function [psnr] = gme_cas_psnr_ls(seq_in, mv_path, seg_path, pC, pR, nFrm, REJRATE, MAXITER)

%% input file
[fid_seq_in message]= fopen(seq_in,'rb');

%% parameters of a synthetic motion field
%
blkSiz = 8; % block size
bC=pC/blkSiz; % number of column in blocks
bR=pR/blkSiz; % number of row in blocks

% get the coordinates
[coorBlkY,coorBlkX]=ndgrid(1:bR,1:bC);
coorX = coorBlkX.*blkSiz-blkSiz/2;
coorY = coorBlkY.*blkSiz-blkSiz/2;

HALFPIX=2; % Half pixel

for ii = 1:nFrm
    % get one frame from YUV sequence (YUV: 411 format) 
    [yCurr, uCurr, vCurr] = readOneFrame(fid_seq_in, pR, pC);
    if ii == 1
        iniH=[];
        fgMap = zeros(bR, bC);
    else
        %% get input noisy motion field
        [npx, npy] = readMVs(mv_path,ii, bC, bR, HALFPIX);

        %% GME with MRF classification out of the GD loop
%         [fgMap, iniMM] = mrfOutGme(npx, npy, coorX, coorY, ...
%             bR, bC, blkSiz, iniMM, fgMap, REJRATE, MAXITER, ii);
        [fgMap, MM, iniH] = mrfInGme(npx, npy, coorX, coorY, ...
            blkSiz, iniH, fgMap, REJRATE, MAXITER, ii, seg_path, yCurr, uCurr, vCurr);
    
        psnr(ii-1)=psnrGMC(yPrev, yCurr, MM); 
%         psnr(ii-1) = 0;
    end
    yPrev = yCurr;  % copy it as a reference frame
end

fclose(fid_seq_in);

% the PSNR after applying global motion compensation.
function psnr = psnrGMC(yPrev, yCurr, M)
yWarp=gmeTF(yCurr, M, 1);
D=yWarp-yPrev;
[pR, pC] = size(D);
mse=sum(D(:).*D(:))/(pC*pR);
psnr=10*log10(255^2/mse); % PSNR for affine motion model


function [fgMap, iniMM] = mrfOutGme(npx, npy, coorX, coorY, ...
    bR, bC, blkSiz, iniMM, fgMap, REJRATE, MAXITER, ii)

REM_MODE = 3;

GM_TRAN = 1;  % translational model
GM_ISOT = 2;  % isotripic model
GM_AFFI = 3;  % affine model
GM_PERS = 4;  % perspective model
REM_MODE = 3;

PO1 = 0.125;      % type-1 outlier --> 10% of MVF
PO2 = 0.125;      % type-2 outleir --> 20% of MVF
PO = PO1 + PO2;  % the overall outleirs --> 30% of MVF

PO_NR = 0.1;      % Newton-Raphson --> 10% of MVF


% GME based on all inliers using Newton-Raphson method  
if ii == 2  % for the first frame, we use more iterations
    iMap = ones(bR, bC);
    mmPers = mvGME_NR_test(GM_PERS, npx(:), npy(:), iMap(:), ...
        coorX(:), coorY(:), 6, REJRATE, iniMM);

    % PSNR for perspective motion model
    iniMM = mmPers(6+1,:);
else
    if REM_MODE == 4
        % the MV outliers removal cascade
        iMap = OutlierRem(npx, npy, PO1, PO2, fgMap, REM_MODE);

        lenPO12 = length(find(iMap(:)==0));
        len = bC*bR*(1-REJRATE)- lenPO12;
        lenNRPO = bC*bR*PO_NR;

        if len > lenNRPO
            nrRejectRate = 1-PO_NR-lenPO12/(bC*bR);
        elseif len < 0
            nrRejectRate = REJRATE;
        else
            nrRejectRate = REJRATE;
        end

        mmPers = mvGME_NR_test(GM_PERS, npx(:), npy(:), iMap(:), ...
            coorX(:), coorY(:), MAXITER, nrRejectRate, iniMM);
    else
        % the MV outliers removal cascade
        iMap = OutlierRem(npx, npy, PO1, PO2, fgMap, REM_MODE);
        mmPers = mvGME_NR_test(GM_PERS, npx(:), npy(:), iMap(:), ...
            coorX(:), coorY(:), MAXITER, REJRATE, iniMM);
    end

    % PSNR for perspective motion model
    iniMM = mmPers(MAXITER+1,:);
end

if REM_MODE == 3 || REM_MODE == 2 || REM_MODE == 4
    [mmpx, mmpy] = mvGen_f(iniMM, coorX, coorY);% SNR check-up
    resMVFx = npx - mmpx;
    resMVFy = npy - mmpy;

    % MV quantization
    [mvSidX, mvSidY, segMap] = initCluster(resMVFx, resMVFy);

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

    % 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);
    end
end


function [fgMap, MM, iniH] = mrfInGme(npx, npy, coorX, coorY, ...
    blkSiz, iniH, fgMap, REJRATE, nITER, IDFRM, SEG_PATH, img_y, img_u, img_v)

REM_MODE = 3;

GM_TRAN = 1;  % translational model
GM_ISOT = 2;  % isotripic model
GM_AFFI = 3;  % affine model
GM_PERS = 4;  % perspective model
REM_MODE = 3;
MRF_ENA = 1;

PO1 = 0.05;      % type-1 outlier --> 10% of MVF
PO2 = 0.40;      % type-2 outleir --> 20% of MVF
PO = PO1 + PO2;  % the overall outleirs --> 30% of MVF

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

if IDFRM == 2
    nITER = 4;  % more iterations for the first frame
end

for ii=1:nITER
   
    % iteratively refine GME
    if ii == 1
        % GME based on all inliers using Newton-Raphson method  
        if IDFRM == 2  % for the first frame, we use more iterations
            % the MV outliers removal cascade
            iMap = OutlierRem(npx, npy, PO1, PO2, fgMap, blkSiz, 3);
            [MM, iniH] = gmeLs(coorX, coorY, npx, npy, iniH, iMap, REJRATE, 0.001);
        else
            % the MV outliers removal cascade
            iMap = 1-fgMap;
            [MM, iniH] = gmeLs(coorX, coorY, npx, npy, iniH, iMap, REJRATE, 0.001);
        end
        
        % object segmentation
        fgMap = resSeg(npx, npy, MM, coorX, coorY, blkSiz, PO2, MRF_ENA, 2);
        
    elseif ii == nITER
        % the MV outliers removal cascade
        iMap = 1-fgMap;
        [MM, iniH] = gmeLs(coorX, coorY, npx, npy, iniH, iMap, REJRATE, 0.001);
        
        % object segmentation
        fgMap = resSeg(npx, npy, MM, coorX, coorY, blkSiz, PO2, MRF_ENA, 1);
        
        
    else
        % the MV outliers removal cascade
        iMap = 1-fgMap;
        [MM, iniH] = gmeLs(coorX, coorY, npx, npy, iniH, iMap, REJRATE, 0.001);
        
        % object segmentation
        fgMap = resSeg(npx, npy, MM, coorX, coorY, blkSiz, PO2, MRF_ENA, 2);
        
    end
end

function [rcoorX, rcoorY] = coorEst(m, coorX, coorY)       
factor = m(7).*coorX+m(8).*coorY+1;
rcoorX = (m(1).*coorX+m(2).*coorY+m(3))./factor;
rcoorY = (m(4).*coorX+m(5).*coorY+m(6))./factor;


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, mode)

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

% moving region threshold
t = 1.2;
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

% field
maxReg = max(segMap(:));
if maxReg > 1
    % if so, then, predict the foreground map
    fgMap = predFgMap(npx, npy, segMap, blkSiz, PO2, mode);
else
    fgMap = zeros(size(segMap));
end
