
function [iMap] = segMRF(MVFx, MVFy, iMap)
[m,n]=size(iMap);
maxReg=max(iMap(:));
beta = 3.5;
for iter = 1:5  % 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



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

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;
