function gmc_seg(seq_in, mv_path, ENA_REGION_MERGE, 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) 
    [img_y, img_u, img_v] = readOneFrame(fid_seq_in, pR, pC);
    if ii == 1
        iniMM=[];
        fgMap = zeros(bR, bC);
    else
        % get input noisy motion field
        [npx, npy] = readMVs(mv_path,ii, bC, bR, HALFPIX);

        % global motion compensation
        [resMVFx, resMVFy, segMap, iniMM, fgMap] = ...
            gmc(npx, npy, coorX, coorY, REJRATE, MAXITER, iniMM, fgMap, ii, ENA_REGION_MERGE);

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

[color_seg_map] = colorSegMap(img_y, segMap, 8);
figure, imshow(color_seg_map);

fclose(fid_seq_in);

function [iMap] = segMRF(MVFx, MVFy, iMap)
[m,n]=size(iMap);
maxReg=max(iMap(:));
beta = 3.5;
for iter = 1:6  % feels like 6 iteration is enough
    % initialize the segmentation map.
    [stat] = initSeg(MVFx(:), MVFy(:), iMap(:), maxReg);

    % initialization of the segmentation map
    for ii = 1:maxReg
        eng(ii,:,:) = log10(2*pi*((stat(ii,2)*stat(ii,4)).^0.5))+...
            0.5.*((MVFx-stat(ii,1)).^2)./stat(ii, 2)+...
            0.5.*((MVFy-stat(ii,3)).^2)./stat(ii, 4);
    end
    
    if iter==1
        init_eng=initEng(MVFx,MVFy,iMap,stat,beta,maxReg);
    end
    
    padMap = wextend(2,'sym',iMap,1);
    for jj=1:m
        for ii=1:n
            for rr = 1:maxReg
                % upper
                if rr==padMap(jj,ii+1);
                    douEng = -beta;
                else
                    douEng = beta;
                end
                % bottom
                if rr==padMap(jj+2,ii+1);
                    douEng = douEng-beta;
                else
                    douEng = douEng+beta;
                end
                % left
                if rr==padMap(jj+1,ii);
                    douEng = douEng-beta;
                else
                    douEng = douEng+beta;
                end
                % right
                if rr==padMap(jj+1,ii+2);
                    douEng = douEng-beta;
                else
                    douEng = douEng+beta;
                end

                localEng(rr) = douEng + eng(rr,jj,ii);
            end
            [C,I]=min(localEng);
            padMap(jj+1,ii+1)=I;
            iterEng(jj,ii)=C;
        end
    end
    iMap=padMap(2:end-1,2:end-1);
    EE(iter) = sum(iterEng(:));
end

% disable for Zeng's testing, and enable otherwise
iMap = noiseSup(iMap, 2);

% EE=[sum(init_eng(:)) EE];
% figure,plot(0:length(EE)-1,EE,'b*-', 'linewidth', 2)
% box on
% grid off
% xlabel('Iterations','FontSize',18)
% ylabel('Posterori energy','FontSize',18)
% title('ICM Optimization: Garden','FontSize',18)
% set(gca,'FontSize',18)


function eng=initEng(MVFx,MVFy,iMap,stat,beta,maxReg)

padMap = wextend(2,'sym',iMap,1);
% upper
diff = iMap-padMap(1:end-2,2:end-1);
doubleTon = ones(size(iMap)).*beta;
doubleTon(diff==0)=-beta;
douEng = doubleTon;
% bottom
diff = iMap-padMap(3:end,2:end-1);
doubleTon = ones(size(iMap))*beta;
doubleTon(diff==0)=-beta;
douEng = douEng+doubleTon;
% left
diff = iMap-padMap(2:end-1,1:end-2);
doubleTon = ones(size(iMap))*beta;
doubleTon(diff==0)=-beta;
douEng = douEng+doubleTon;
% right
diff = iMap-padMap(2:end-1,3:end);
doubleTon = ones(size(iMap))*beta;
doubleTon(diff==0)=-beta;
douEng = douEng+doubleTon;

uMVx=zeros(size(iMap));
uMVy=zeros(size(iMap));
sigmaMVx=zeros(size(iMap));
sigmaMVy=zeros(size(iMap));
for ii = 1:maxReg
    uMVx(iMap==ii)=stat(ii,1);
    uMVy(iMap==ii)=stat(ii,3);
    sigmaMVx(iMap==ii)=stat(ii,2);
    sigmaMVy(iMap==ii)=stat(ii,4);
end

sinEng = log10(2*pi*((sigmaMVx.*sigmaMVy).^0.5))+...
    0.5.*((MVFx-uMVx).^2)./sigmaMVx+0.5.*((MVFy-uMVy).^2)./sigmaMVy;

eng = douEng+sinEng;

function [stat] = initSeg(inMVx, inMVy, iMap, maxReg)
for ii = 1:maxReg
    [R, C] = find(iMap(:)==ii); % get the region index
    stat(ii, 1) = mean(inMVx(R)); % mean of MVx;
    err = (inMVx(R)-stat(ii,1)).^2;
    stat(ii, 2) = mean(err); % variance of MVx;
    stat(ii, 3) = mean(inMVy(R)); % mean of MVy;
    err = (inMVy(R)-stat(ii,3)).^2;
    stat(ii, 4) = mean(err); % variance of MVy;
end

% the PSNR after applying global motion compensation.
function [resMVFx, resMVFy, segMap, iniMM, fgMap] = ...
            gmc(npx, npy, coorX, coorY, REJRATE, MAXITER, iniMM, fgMap, ii, ENA_REGION_MERGE)

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.10;      % type-1 outlier --> 10% of MVF
PO2 = 0.15;      % type-2 outleir --> 20% of MVF
PO = PO1 + PO2;  % the overall outleirs --> 30% of MVF

blkSiz = 8;

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

    % PSNR for perspective motion model
    iniMM = mmPers(6+1,:);
else
    % the MV outliers removal cascade
    iMap = OutlierRem(npx, npy, PO1, PO2, fgMap, blkSiz, REM_MODE);
    mmPers = mvGME_NR_test(GM_PERS, npx(:), npy(:), iMap(:), ...
        coorX(:), coorY(:), MAXITER, REJRATE, iniMM);
    % PSNR for perspective motion model
    iniMM = mmPers(MAXITER+1,:);
end

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

% moving region threshold
t = 1.2;
if ENA_REGION_MERGE == 1
    thX = t*(mean(resMVFx(:).^2)^0.5);
    thY = t*(mean(resMVFy(:).^2)^0.5);
else
    thX = 0;
    thY = 0;
end

    % MV quantization
%     [mvSidX, mvSidY, segMap] = initCluster(resMVFx, resMVFy);
    [segMap] = initCluster(resMVFx, resMVFy, 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, 1);
    end
end

function [segMap] = regMerge(MVFx, MVFy, segMap)
[m, n] = size(segMap);
% labeling regions
[C, R] = find(segMap(:)==0);

pNum = m*n-length(C);
regID = 0;
mapID = 1;
outMap = zeros(m, n);
while pNum > 0
    z = zeros(size(segMap));
    z(segMap(:)==mapID) = 1;
    
    mapID = mapID+1;
    
    regZ = z.*regID;
    L = bwlabel(logical(z), 4);
    maxReg = max(L(:));
    L = L+regZ;
    regID = regID + maxReg;
    pNum = pNum-sum(z(:));
    outMap = outMap + L;
end

for ii = 1:regID
    [R, C] = find(outMap(:)==ii);
    if length(C) < 0.01*m*n
        outMap(R) = 0;
    end
end

% re-ordering all regions
[segMap] = regOrdering(outMap, MVFx, MVFy);
segMap=segMap+1;
numTH = [0.3, 0.4, 0.6, 0.8];
segMap = regGrowing(MVFx, MVFy, segMap, numTH);
segMap=segMap-1;

[C, R] = find(segMap(:)==0);
if ~isempty(C)
    v = length(C);
end

function [iMap] = regOrdering(iMap, inMVx, inMVy)

% region re-ordering accorindg to its motion size
maxReg = max(iMap(:));
tstMap = iMap;
for ii = 1:maxReg
    [R, C] = find(iMap(:)==ii); % get the region index
    if ~isempty(R)
        mMVx = mean(inMVx(R));
        mMVy = mean(inMVy(R));
        mMVSiz(ii) = mMVx.^2+mMVy.^2;
    else
        mMVSiz(ii) = 0;
    end
end
[b, ix] = sort(mMVSiz,'descend');

for ii = 1:maxReg
    newID = find(ix==ii);% get the region index
    iMap(tstMap==ii) = newID;
end

function [iMap] = regGrowing(MVFx, MVFy, iMap, numTH)

[m, n] = size(iMap);
for nLoop = 1:20
    [cR, cC] = find(iMap==1);
    if isempty(cR)
        break;
    elseif nLoop > length(numTH)
        gTH = 1;
    else
        gTH = numTH(nLoop);
    end
        
    dReg = []; neiID = [];
    
    % calculate the threshold
%     x = norminv([conf (1-conf)],0,1);
%     gTH = x(2);
    
    % get 4-adjacency block address
    [cID, lID, rID, uID, bID] = neiBlkAdd(cR, cC, m);

    % get the statistical information
    stat = regStat(MVFx(:), MVFy(:), iMap(:));
    padMap = wextend(2,'sym',iMap,1);
    mvx = MVFx(cID);
    mvy = MVFy(cID);

    % the distance between input block and the neighboring regions
    [dReg(:,1), neiID(:,1)] = regDist(mvx, mvy, stat, padMap(lID));
    [dReg(:,2), neiID(:,2)] = regDist(mvx, mvy, stat, padMap(rID));
    [dReg(:,3), neiID(:,3)] = regDist(mvx, mvy, stat, padMap(uID));
    [dReg(:,4), neiID(:,4)] = regDist(mvx, mvy, stat, padMap(bID));

    % get the region with minimum distances
    [C, I] = min(dReg, [], 2);
    
    if gTH == 1,
        id = 1:length(cR);
        iMap(cID(id)) = neiID(id'+(I-1).*length(cR));
    else
        [f,x] = ecdf(C(:));
        z = find(f>gTH);
        id = find(C<x(z(1)));
        iMap(cID(id)) = neiID(id+length(cR)*(I(id)-1));
    end
    
    % for testing
%     testSegMap(iMap);
end

function [dist, id] = regDist(mvx, mvy, stat, id)
dist = ((mvx-stat(id,2)).^2+(mvy-stat(id,3)).^2)./(stat(id,4).^2);

function [cID, lID, rID, uID, bID] = neiBlkAdd(cR, cC, m)
cID = cR+(cC-1)*m;
lR = cR+1; lC = cC; lID = lR+(lC-1)*(m+2); % the left block
rR = cR+1; rC = cC+2; rID = rR+(rC-1)*(m+2); % the right block
uR = cR; uC = cC+1; uID = uR+(uC-1)*(m+2); % the upper block
bR = cR+2; bC = cC+1; bID = bR+(bC-1)*(m+2); % the bottom block

function [stat] = regStat(inMVx, inMVy, iMap)
maxReg = max(iMap);
% initialize the null region, set the ungrouped region an extremely large
% MV.
stat(1, 1) = 1;
stat(1, 2) = 2555; % mean MV in horizontal
stat(1, 3) = 2555; % mean MV in vertical
stat(1, 4) = 0;   % the variance
for ii = 2:maxReg
    [R, C] = find(iMap==ii); % get the region index
    stat(ii, 1) = length(R);
    stat(ii, 2) = mean(inMVx(R));
    stat(ii, 3) = mean(inMVy(R));
    err = (inMVx(R)-stat(ii,2)).^2+(inMVy(R)-stat(ii,3)).^2;
    mErr = mean(err);
    % 0.1 is the factor to prevent mErr goes to 0
    stat(ii, 4) = 0.1+mErr.^0.5; 
end
