
function iter_seg_LS(px, py, para, MAX_ITER, nH, nW, coorX, coorY, blkSiz)

% 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


fgMV = -5; % MVx, MVy --> 3; emulating a moving object (3, 3)

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

stdDev = para(1);   % the standard deviation of Gaussian noise
fgSiz = para(2);  % the foreground size

fgN=round((fgSiz*nH*nW).^0.5);
strPixH=0;  % starts from upper-left corner
strPixW=0;
% initialize foreground map
fgMap = zeros(size(px));

ITER_RANSAC = 500;


% add moving object on the top of MV field by GM model.
strBlkH = max(round(strPixH/blkSiz), 1);
strBlkH = min(strBlkH, nH-fgN+1);
strBlkW = max(round(strPixW/blkSiz), 1);
strBlkW = min(strBlkW, nW-fgN+1);

endBlkH = strBlkH+fgN-1;
endBlkW = strBlkW+fgN-1;

fgMVF = zeros(size(px));
fgMVF(strBlkH:endBlkH,strBlkW:endBlkW) = fgMV;
opx = px+fgMVF;
opy = py+fgMVF;

% object moving from upper left to bottom right
strPixH = strPixH - fgMV;  
strPixW = strPixW - fgMV;

% 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

MAXITER = MAX_ITER(1);
[iterMM] = SEGGME(npx, npy, coorX, coorY, nH, nW, blkSiz, MAXITER);


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

MAXITER = 1;
REM_MODE = 3;
MM=[];
PREV_MM = [];
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);
dM_th = [1e-5, 1e-5, 1e-3, 1e-5, 1e-5, 1e-3, 1e-5, 1e-5];

TEST_SEG = 1:20;

ENA_SEG_EVA = 1;

for ii=1:nITER+1
    
    if ii == 1
        % initializing the parameter
        MM = zeros(1,8);
        MM(1) = 1;
        MM(5) = 1;
        MM(3) = sum(npx(:))/(bR*bC);
        MM(6) = sum(npy(:))/(bR*bC);
        iterMM(ii, :)= MM;
    elseif ii == 2
        MM = [];
        % the MV outliers removal cascade
        iMap = type1Map;
        
        MM = gmeLsme(coorX, coorY, npx, npy, iMap, MM);
        
        % PSNR for perspective motion model
        iterMM(ii, :)= MM;
        
        % object segmentation
        fgMap = resSeg(npx, npy, MM, coorX, coorY, blkSiz, PO2, MRF_ENA);
    elseif ii == nITER+1
%         iMap = ones(bR, bC);
        iMap = 1-fgMap;
        MM = gmeLsme(coorX, coorY, npx, npy, iMap, MM);
        iterMM(ii, :) = MM;
    else
        % the MV outliers removal cascade
        if ENA_SEG_EVA
            if find(TEST_SEG==ii)
                colorSegMap(type1Map, fgMap, 16);
            end
        end

        iMap = 1-fgMap;
        
        MM = gmeLsme(coorX, coorY, npx, npy, iMap, MM);
        
        % PSNR for perspective motion model
        iterMM(ii, :) = MM;
                
        % object segmentation
        fgMap = resSeg(npx, npy, MM, 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

function colorSegMap(t1Map, t2Map, ext)
[m, n] = size(t1Map);
d1 = ones(m, n).*255; % three color R, G, B,
d2 = ones(m, n).*255; % three color R, G, B,
d3 = ones(m, n).*255; % three color R, G, B,

id = find(t1Map==0);
d1(id) = 0;
d2(id) = 0;
id = find(t2Map==1);
d1(id) = 0;
d3(id) = 0;

RGB(:,:,1) = kron(d1, ones(ext));
RGB(:,:,2) = kron(d2, ones(ext));
RGB(:,:,3) = kron(d3, ones(ext));

figure, imshow(uint8(RGB));
