function gme_synth_mvf_ls();

clear all
close all
warning off all

fprintf(1,'****************************************************\n');
fprintf(1,'*      GME performance analysis on synthetic MVF   *\n');
fprintf(1,'****************************************************\n');
fprintf(1,'       Designed by Yue-Meng Chen, Jul,  2009      \n\n');

%% parameters of a synthetic motion field
%
blkSiz = 8; % block size
MVRES = 1;  % integer-pixel: 1, half-pixel: 2
height = 288; nH = height/blkSiz;
width =352; nW = width/blkSiz;

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

% four global motion models to be tested!
m(1, :) =   [0.95     0        10.4238    0      0.95      5.7927    0           0];
m(2, :) =   [0.9964   -0.0249  1.0981     0.0856 0.9457    -7.2      0           0];
m(3, :) =   [0.9964   -0.0249  6.0981     0.0249 0.9964    2.5109    -2.7e-5     1.9e-5];
m(4, :) =   [1        0        4.4154     0      1         0         -1.1263e-4  0];


% constructing a object map
% three objects with sizes: 10%, 20%, and 5%
objF = ones(nH, nW);
fgSiz = [0.10, 0.20, 0.05];  % object with size up to 25% of overall MV field
fgN=round((fgSiz(1)*nH*nW).^0.5);

% first object -- square region
strBlkH = 6;
strBlkW = 12; %(nW-fgN*2)/2;
endBlkH = strBlkH+fgN-1;
endBlkW = strBlkW+fgN-1;
id = 2;
objF(strBlkH:endBlkH,strBlkW:endBlkW) = id;

% second object -- Circular region
r = round((fgSiz(2)*nH*nW/pi).^0.5);
x = endBlkW+r;
y = floor(nH/2);
id = 3;
for jj = 1:nH
    for ii = 1:nW
        dis = round(((jj-y)^2+(ii-x)^2)^0.5);
        if dis < r
            objF(jj, ii) = id;
        end
    end
end

% the third object -- triangular region
id = 4;
tri_strBlkH = endBlkH+1;
tri_strBlkW = strBlkW;
tri_fgSiz = round((fgSiz(3)*nH*nW).^0.5);
for ii = 1:tri_fgSiz
    objF(tri_strBlkH+ii-1,tri_strBlkW-(ii-1):tri_strBlkW+(ii-1)) = id;
end
testSegMap(objF, 8, 1);

% foreground obj motion
mvObj = [3, -7; -6, 4; 5, 5];
GMID = 1;
% generate motion field in floating point
[MVFx, MVFy] = mvGen_f(m(GMID, :), coorX, coorY);

for REGID = 1:3
    
    id = find(objF==REGID+1);
    
    MVFx(id) = MVFx(id) + mvObj(REGID, 1);
    MVFy(id) = MVFy(id) + mvObj(REGID, 2);
end

stdDev = 2.2;
nX = stdDev.*randn(nH, nW);
npx = round(MVFx+nX);
nY = stdDev.*randn(nH, nW);
npy = round(MVFy+nY);

figure, quiver(1:nW, nH:-1:1, npx, npy);
xlim([0 45]), ylim([0 37]);

% verify the size of foreground/background objects
for ii = 1:4
    id = find(objF==ii);
    perc(ii) = round((length(id)/(nH*nW))*100);
end

sd = [0.7, 1.5, 2.0, 3.0];  % stdDev: 2.2, fgSiz: 0.10
         
% the parameters for the cascade
MAXITER = [20 20 20 500];
TEST_SCENARIO = [2];  % change to different scenarios
ENA_FIGURE = 0;

for ii = 1:4 % Four types of outlier sizes

    if find(TEST_SCENARIO==ii)

        disp(sprintf('scenario: %d,', ii));

        stdDev = sd(ii);
        nX = stdDev.*randn(nH, nW);
        npx = round(MVFx+nX);
        nY = stdDev.*randn(nH, nW);
        npy = round(MVFy+nY);
        
        % test 4 algorithms with different noise level and object size
        iter_seg_LS_mvf(npx, npy, MAXITER, nH, nW, coorX, coorY, blkSiz);
    end
end

disp('done!')

function plotGmeLs(propSNR, gdSNR, fltSNR, ranSNR, GMID, para)

figure, % Sigma: 0.7
plot(1:50, propSNR(:,end), 'r-', 'LineWidth',2)
hold on
plot(1:50, gdSNR(:,end), 'b:', 'LineWidth',2)
hold on
plot(1:50, fltSNR(:,end), 'k-.', 'LineWidth',2)
hold on  
plot(1:50, ranSNR(:,end), 'g--', 'LineWidth',2)
hold off
legend('SEG\_LS','LSS-MS', 'FLT\_LS', 'RAN\_LS')
box on
grid off
xlabel('Frame index','FontSize',18)
ylabel('Average SNR (dB)','FontSize',18)
title(sprintf('GM %d, sd=%1.1f, outliers: %d%%', ...
    GMID,para(1), para(2)*100) ,'FontSize',18)
set(gca,'FontSize',18)
