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

%% 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

iniMM=[];
fgMap = zeros(bR, bC);

%% get input noisy motion field
[npx, npy] = readMVs(mv_path,nFrm, bC, bR, HALFPIX);

%% GME with MRF classification out of the GD loop
mrfInGme(npx, npy, coorX, coorY, bR, bC, blkSiz, iniMM, fgMap, REJRATE);
    

function mrfInGme(npx, npy, coorX, coorY, bR, bC, blkSiz, iniMM, fgMap, REJRATE)

GM_PERS = 4;  % perspective model
MRF_ENA = 1;

% GME based on all inliers using Newton-Raphson method  
iMap = ones(bR, bC);
mmPers = mvGME_NR(GM_PERS, npx(:), npy(:), iMap(:), ...
    coorX(:), coorY(:), 6, REJRATE, iniMM);

% get motion parameters of perspective motion model
iniMM = mmPers(6+1,:);

% object segmentation
fgMap = resSeg(npx, npy, iniMM, coorX, coorY, blkSiz, 0, MRF_ENA, 2);
        

% 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
mvQuant(resMVFx, resMVFy, thX, thY);

        