function [wSNR, woSNR, fltSNR, ranSNR] = ...
    testNoiObjGME_GD_multiple(opx, opy, stdDev, m_in, MAX_ITER, nH, nW, coorX, coorY, blkSiz)

ENA_PROPOSED = 1;
ENA_PURE_GD = 1;
ENA_FLT_GD = 1;
ENA_RANSAC = 1;

[px, py] = mvGen_f(m_in, coorX, coorY);
figure, quiver(px, py)
figure, quiver(opx, opy)
% The outlier removal mode:
%   1: removing type-1 outlier
%   2: removing type-2 outlier
%   3: removing both type-1 and type-2 outlier
%
REM_MODE = 3;

PO1 = 0.05;      % type-1 outlier --> 10% of MVF
PO2 = 0.20;      % type-2 outleir --> 20% of MVF
PO = PO1 + PO2;  % the overall outleirs --> 30% of MVF
rOUTLIERS = 0.65; % estimate of outliers in a motion field.
REJRATE = 0.65; % the ratio of inliers to overall MVs

% the global motion model
GMMODE = 4; % 1: Translational, 2: Isotropic, 3: Affine, 4: Perspective

% initialize foreground map
fgMap = zeros(size(opx));

ITER_RANSAC = 500;

for runs = 1:1:50

    % Gaussian noise on X/Y channel, and rounding it to integer
    % pixel
    nX = stdDev.*randn(nH, nW);
    npx = round(opx+nX);
    nY = stdDev.*randn(nH, nW);
    npy = round(opy+nY);

    %% GME with the proposed outlier removal
    % inlier map initialization
    if ENA_PROPOSED == 1
        MAXITER = MAX_ITER(1);
        [iterMM] = SEGGME(npx, npy, coorX, coorY, nH, nW, blkSiz, MAXITER, REJRATE);
        for ii = 1:MAXITER+1
            [mmpx, mmpy] = mvGen_f(iterMM(ii,:), coorX, coorY);% SNR check-up
            wSNR(runs, ii) = mvSNR(px, py, mmpx, mmpy, ones(nH, nW));
        end
    else
        wSNR(runs, 1:MAXITER+1) = 0;
    end


    %% pure GD GME 
    if ENA_PURE_GD == 1
        MAXITER = MAX_ITER(2);
        iMap = ones(size(px));
        mm1 = mvGME_NR(GMMODE, npx(:), npy(:), iMap(:), ...
            coorX(:), coorY(:), MAXITER, REJRATE, []);
        for ii = 1:MAXITER+1
            [mmpx, mmpy] = mvGen_f(mm1(ii,:), coorX, coorY);% SNR check-up
            woSNR(runs, ii) = mvSNR(px, py, mmpx, mmpy, ones(nH, nW));
        end
    else
        woSNR(runs, 1:MAXITER+1) = 0;
    end


    %% GME with filter from Ref. [6]
    % inlier map initialization
    if ENA_FLT_GD == 1
        MAXITER = MAX_ITER(3);
        iMap = MVRemFlt(npx, npy);
        % check if at least 4 MVs obtained after filtering
        if sum(iMap(:)) > 4

            % GME based on all inliers using Newton-Raphson method
            mm0 = mvGME_NR(GMMODE, npx(:), npy(:), iMap(:), ...
                coorX(:), coorY(:), MAXITER, REJRATE, []);
            for ii = 1:MAXITER+1
                [mmpx, mmpy] = mvGen_f(mm0(ii,:), coorX, coorY);% SNR check-up
                fltSNR(runs, ii) = mvSNR(px, py, mmpx, mmpy, ones(nH, nW));
            end
        else
            if runs == 1
                fltSNR(runs, 1:MAXITER+1) = 0;
            else
                fltSNR(runs, :) = fltSNR(runs-1, :);
            end
        end
    else
        fltSNR(runs, 1:MAXITER+1) = 0;
    end

    %% GME with RANSAC
    % get the reference coordinates
    if ENA_RANSAC == 1
        MAXITER = MAX_ITER(4);
        
        rcoorX = coorX + npx;
        rcoorY = coorY + npy;
        
        if MAXITER > 100
                [ranM, nIter, inliers] = ranGME(GMMODE, rcoorX(:), ...
                    rcoorY(:), coorX(:), coorY(:), ii);
                iMap = zeros(size(npx));
                iMap(inliers) = 1;

                % check if at least 4 MVs obtained after filtering
                if sum(iMap(:)) > 4

                    % GME based on all inliers using Newton-Raphson method
                    % using 3 iterations
                    ITER_GD = 3;
                    mm0 = mvGME_NR(GMMODE, npx(:), npy(:), iMap(:), ...
                        coorX(:), coorY(:), ITER_GD, REJRATE, []);

                    [mmpx, mmpy] = mvGen_f(mm0(ITER_GD+1 ,:), coorX, coorY);% SNR check-up
                    rsnr = mvSNR(px, py, mmpx, mmpy, ones(nH, nW));
                    ranSNR(runs, 1) = max(0, rsnr);
                else
                    if runs == 1
                        ranSNR(runs, 1) = 0;
                    else
                        ranSNR(runs, 1) = ranSNR(runs-1, 1);
                    end
                end
        else
            for ii = 1:MAXITER+1
                [ranM, nIter, inliers] = ranGME(GMMODE, rcoorX(:), ...
                    rcoorY(:), coorX(:), coorY(:), ii);
                iMap = zeros(size(npx));
                iMap(inliers) = 1;

                % check if at least 4 MVs obtained after filtering
                if sum(iMap(:)) > 4

                    % GME based on all inliers using Newton-Raphson method
                    % using 3 iterations
                    ITER_GD = 3;
                    mm0 = mvGME_NR(GMMODE, npx(:), npy(:), iMap(:), ...
                        coorX(:), coorY(:), ITER_GD, REJRATE, []);

                    [mmpx, mmpy] = mvGen_f(mm0(ITER_GD+1 ,:), coorX, coorY);% SNR check-up
                    rsnr = mvSNR(px, py, mmpx, mmpy, ones(nH, nW));
                    if ii == 1
                        ranSNR(runs, ii) = max(0, rsnr);
                    else
                        ranSNR(runs, ii) = max(ranSNR(runs,ii-1), rsnr);
                    end
                else
                    if runs == 1
                        ranSNR(runs, 1:MAXITER+1) = 0;
                    else
                        ranSNR(runs, :) = ranSNR(runs-1, :);
                    end
                end
            end
        end
    else
        ranSNR(runs, 1:MAXITER+1) = 0;
    end
end

mPropSnr = mean(wSNR);
mPureGDSnr = mean(woSNR);
mRANSnr = mean(ranSNR);
mFltGDSnr = mean(fltSNR);

function [iterMM] = SEGGME(npx, npy, coorX, coorY, bR, bC, blkSiz, nITER, REJRATE)

MAXITER = 3;
REM_MODE = 3;
iniMM=[];
fgMap = zeros(bR, bC);
PO1 = 0.05;      % type-1 outlier --> 10% of MVF
PO2 = 0.20;      % type-2 outleir --> 20% of MVF
PO = PO1 + PO2;  % the overall outleirs --> 30% of MVF
REJRATE = 0.65; % the ratio of inliers to overall MVs
MRF_ENA = 1;
% the global motion model
GMMODE = 4; % 1: Translational, 2: Isotropic, 3: Affine, 4: Perspective

TYPE1PO = 0.15;
type1Map = OutlierRem(npx, npy, TYPE1PO, PO2, fgMap, blkSiz, 1);

for ii=1:nITER
   
    % iteratively refine GME
    if ii == 1
        % the MV outliers removal cascade
%         iMap = OutlierRem(npx, npy, PO1, PO2, fgMap, blkSiz, REM_MODE);
        iMap = type1Map;
        mmPers = mvGME_NR(GMMODE, npx(:), npy(:), iMap(:), ...
            coorX(:), coorY(:), MAXITER, REJRATE, iniMM);
        % PSNR for perspective motion model
        iniMM = mmPers(MAXITER+1,:);
        iterMM = mmPers;
                
        % object segmentation
        fgMap = resSeg(npx, npy, iniMM, coorX, coorY, blkSiz, PO2, MRF_ENA);
        
    elseif ii == nITER
%         iMap = ones(bR, bC);
        iMap = 1-fgMap;
        mmPers = mvGME_NR(GMMODE, npx(:), npy(:), iMap(:), ...
            coorX(:), coorY(:), MAXITER, REJRATE, iniMM);
        % PSNR for perspective motion model
        iniMM = mmPers(MAXITER+1,:);
        iterMM(ii+1, :) = iniMM;
    else
        % the MV outliers removal cascade
%         iMap = OutlierRem(npx, npy, PO1, PO2, fgMap, REM_MODE);
        iMap = 1-fgMap;
        iMap = type1Map.*iMap;
        mmPers = mvGME_NR(GMMODE, npx(:), npy(:), iMap(:), ...
            coorX(:), coorY(:), MAXITER, REJRATE, iniMM);
        % PSNR for perspective motion model
        iniMM = mmPers(MAXITER+1,:);
        iterMM(ii+1, :) = iniMM;
                
        % object segmentation
        fgMap = resSeg(npx, npy, iniMM, coorX, coorY, blkSiz, PO2, MRF_ENA);
    end
end

function segMap = preproMap(npx, npy, segMap, thX, thY)
maxId = max(segMap(:));

if maxId == 1
    return;
end

for jj = 1:maxId
    [idr, idc]=find(segMap(:) == jj);
    regInfo(jj,1) = length(idr);
    regInfo(jj,2) = mean(npx(idr));
    regInfo(jj,3) = mean(npy(idr));
end

% for test
regID = 2;
mapBak = segMap;

magSid = regInfo(:,2).^2+regInfo(:,3).^2;
absSiz = abs(regInfo(:,2))+abs(regInfo(:,3));
for jj = 1:maxId
    if regInfo(:,1) == 0
        continue;
    end
%         if magSid(jj) < 8
    if abs(regInfo(jj,2)) < thX
%             if absSiz(jj) < 4
        if abs(regInfo(jj,3)) < thY
            segMap(mapBak==jj) = 1;
        else
            segMap(mapBak==jj) = regID;
            regID = regID + 1;
        end
    else
        segMap(mapBak==jj) = regID;
        regID = regID + 1;
    end
end

% segmentation performed on residual MV field
function fgMap = resSeg(npx, npy, MM, coorX, coorY, blkSiz, PO2, MRF_ENA)

% object segmentation
[mmpx, mmpy] = mvGen_f(MM, coorX, coorY);% SNR check-up
resMVFx = npx - mmpx;
resMVFy = npy - mmpy;

% moving region threshold
t = 1.5;
thX = t*(mean(resMVFx(:).^2)^0.5);
thY = t*(mean(resMVFy(:).^2)^0.5);

% MV quantization
[segMap] = initCluster(resMVFx, resMVFy, thX, thY);

% segmentation using Markov Random Field classificaiton
if MRF_ENA
    [segMap] = segMRF(resMVFx, resMVFy, segMap);
end

% test map
segMap = preproMap(resMVFx, resMVFy, segMap, thX, thY);

% check whether there is a foreground present in MV
% field
maxReg = max(segMap(:));
if maxReg > 1
    % if so, then, predict the foreground map
    fgMap = predFgMap(npx, npy, segMap, blkSiz, PO2, 2);
else
    fgMap = zeros(size(segMap));
end