
% REJRATE = 0.7; % the ratio of inliers to overall MVs for GD based approach

function iMap = MVRemFlt(npx, npy)

% the parameters for the cascade
TH1 = 0.08;          % the threshold on magnitude -- 40 percent
TH2 = 0.03;          % the threshold on magnitude -- 40 percent

% symmetric extension
expx = wextend('2D','sym',npx,1);
expy = wextend('2D','sym',npy,1);

[m, n] = size(npx);

% magnitude-based filtering
% get the threshold
THSQ = TH1^2;
magRes1 = magFlt1(expx, expy, THSQ, m, n);
magRes1(magRes1<3)=0;
THSQ = TH2^2;
magRes2 = magFlt2(expx, expy, THSQ, m, n);
magRes2(magRes2<2)=0;

iMap = magRes1+magRes2;
iMap(iMap>0)=1;

function numM = magFlt1(px, py, THM2, m, n)

% calculate the magnitude of input MVs
oriM = (px.^2+py.^2)*THM2;
oriM = oriM(2:end-1, 2:end-1); % removing all image boundary pixels
numM = zeros(m, n);

% horizontal difference
ehpx = px(:, 1:end-1)-px(:, 2:end);
ehpy = py(:, 1:end-1)-py(:, 2:end);
ehM = (ehpx.^2+ehpy.^2);
numM = inliers(oriM, ehM(2:end-1, 2:end), numM);
numM = inliers(oriM, ehM(2:end-1, 1:end-1), numM);

% vertical difference
evpx = px(1:end-1, :)-px(2:end, :);
evpy = py(1:end-1, :)-py(2:end, :);
evM = (evpx.^2+evpy.^2);
numM = inliers(oriM, evM(2:end, 2:end-1), numM);
numM = inliers(oriM, evM(1:end-1, 2:end-1), numM);

% diagonal difference (45C)
ed1px = px(1:end-1, 1:end-1)-px(2:end, 2:end);
ed1py = py(1:end-1, 1:end-1)-py(2:end, 2:end);
ed1M = (ed1px.^2+ed1py.^2);
numM = inliers(oriM, ed1M(2:end, 2:end), numM);
numM = inliers(oriM, ed1M(1:end-1, 1:end-1), numM);

% diagonal difference (135C)
ed2px = px(1:end-1, 2:end)-px(2:end, 1:end-1);
ed2py = py(1:end-1, 2:end)-py(2:end, 1:end-1);
ed2M = (ed2px.^2+ed2py.^2);
numM = inliers(oriM, ed2M(2:end, 1:end-1), numM);
numM = inliers(oriM, ed2M(1:end-1, 2:end), numM);

function numM = magFlt2(px, py, THM2, m, n)

% calculate the magnitude of input MVs
oriM = (px.^2+py.^2)*THM2;
oriM = oriM(2:end-1, 2:end-1); % removing all image boundary pixels
numM = zeros(m, n);

% horizontal difference
mpx = (px(:, 1:end-2)+px(:, 3:end))/2;
mpy = (py(:, 1:end-2)+py(:, 3:end))/2;
ehpx = mpx-px(:, 2:end-1);
ehpy = mpy-py(:, 2:end-1);
ehM = (ehpx.^2+ehpy.^2);
numM = inliers(oriM, ehM(2:end-1, :), numM);

% vertical difference
mpx = (px(1:end-2, :)+px(3:end, :))/2;
mpy = (py(1:end-2, :)+py(3:end, :))/2;
evpx = mpx-px(2:end-1, :);
evpy = mpy-py(2:end-1, :);
evM = (evpx.^2+evpy.^2);
numM = inliers(oriM, evM(:, 2:end-1), numM);

% diagonal difference (45C)
mpx = (px(1:end-2, 1:end-2)+px(3:end, 3:end))/2;
mpy = (py(1:end-2, 1:end-2)+py(3:end, 3:end))/2;
ed1px = mpx-px(2:end-1, 2:end-1);
ed1py = mpy-py(2:end-1, 2:end-1);
ed1M = (ed1px.^2+ed1py.^2);
numM = inliers(oriM, ed1M, numM);

% diagonal difference (135C)
mpx = (px(1:end-2, 3:end)+px(3:end, 1:end-2))/2;
mpy = (py(1:end-2, 3:end)+py(3:end, 1:end-2))/2;
ed2px = mpx-px(2:end-1, 2:end-1);
ed2py = mpy-py(2:end-1, 2:end-1);
ed2M = (ed2px.^2+ed2py.^2);
numM = inliers(oriM, ed2M, numM);

function numM = inliers(oriM, eM, numM)
ee = oriM - eM;
numTmp = zeros(size(numM));
numTmp(ee>0) = 1; % the MVs with a less magnitude change
numM = numM + numTmp;
