function [segMap] = initCluster(resMVFx, resMVFy, thX, thY)

% for test
[m, n] = size(resMVFx);
cX = resMVFx(:);
cY = resMVFy(:);
W = ones(size(cX));

% set the centroid of first cluster as (0, 0)
mvSidX(1) = 0;
mvSidY(1) = 0;
D_ALL=W.*(cX.^2+cY.^2);

% get the overall distortion of first cluster, which includes all MVs
ii = 1;
E(ii)=sum(D_ALL(:));

% initialize the segmetnation map
segMap = ones(size(W));
prevMap = segMap;

% set target region to split
regSplit = 1;

while (1)
    
    % find out the regions of intereste
    rM = find(segMap(:) == regSplit);
    
    % sort it out in descent order
    [B, IX] = sort(D_ALL(rM), 'descend');
    
    % split target cluster into two cluster, and initialize their centroid
    
    sizNewCluster = round(length(rM)/2);
    
    % first half
    inId = IX(1:sizNewCluster);
    rM1 = rM(inId);
    newSidX=sum(W(rM1).*cX(rM1))/sum(W(rM1));
    newSidY=sum(W(rM1).*cY(rM1))/sum(W(rM1));
    
    % second half
    if regSplit == 1
        % we fix the centroid of the first cluster as (0, 0) throughout
        % the MV quantization.
    else
        inId = IX(sizNewCluster:end);
        rM2 = rM(inId);
        mvSidX(regSplit)=sum(W(rM2).*cX(rM2))/sum(W(rM2));
        mvSidY(regSplit)=sum(W(rM2).*cY(rM2))/sum(W(rM2));
    end
        
    
    % obtain new set of centroid for weighted vector quantization
    sidX=[mvSidX, newSidX];
    sidY=[mvSidY, newSidY];
    
    ii=ii+1;
    [D_REG, D_ALL, mvSidX, mvSidY, segMap]  = vqMVF(cX, cY, W, sidX, sidY, 5);
    E(ii)=sum(D_ALL(:));
    
    % check whether the Vector Quantization should be ended
    stop = 0;

    % condition - 1,
    % the region centroid is too close to go further on the split
    for jj = 1:length(mvSidY)
        [idr, idc]=find(segMap(:) == jj);
        regInfo(jj,1) = length(idr);
        if jj == 1
            regInfo(jj,2) = 0;
            regInfo(jj,3) = 0;
        else
            regInfo(jj,2) = mean(cX(idr));
            regInfo(jj,3) = mean(cY(idr));
        end
    end

    % condition - 2, 
    % the small distoriton --> the E_n < 0.05 * E_1
    if ii > 2
        if abs(E(ii)-E(ii-1)) < 0.05*abs(E(2)-E(1))
            break;
        end
    end
    
    % condition - 3,
    % the to-be-split region is smaller than 10% of overall region, 
    % or any other region is smaller than 5% of overall MVF size
    
	[C, I] = max(D_REG); % find the index of the cluster with max distortion.
    
    stop = 0;
    for jj = 1:length(mvSidY)
        regSiz=regInfo(jj,1);
        if jj == I
%             TH = 0.1.*m.*n;
            TH = 0.05.*m.*n;
        else
%             TH = 0.05.*m.*n;
            TH = 0.025.*m.*n;
        end
        
        if regSiz<TH
            stop = 1;
            break;
        end
    end
    
    if stop == 1
        break;
    else
        regSplit = I;
    end
    
    
    if ii == 6
        break;
    end
end

% for test
regID = 2;
mapBak = segMap ;
maxRegId = max(segMap(:));
if maxRegId > 1
    magSid = regInfo(:,2).^2+regInfo(:,3).^2;
    absSiz = abs(regInfo(:,2))+abs(regInfo(:,3));
    for jj = 2:length(mvSidY)
        if regInfo(jj,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
end

segMap = reshape(segMap, m, n);

% plot energy curve interms of iterations
% plotEng(E);

function plotEng(E)

figure,plot(E,'b*-', 'linewidth', 2)
box on
grid off
xlabel('the number of motion clusters','FontSize',16)
ylabel('Quantization error','FontSize',16)
title('step-wise MV quantization: Garden','FontSize',16)
set(gca,'FontSize',16)

% the returning map has the index same as the input MV centroid.
function [D_REG, D_ALL, sidX, sidY, segMap]  = vqMVF(resMVFx, resMVFy, W, sidX, sidY, numIter)

% the quantization noise
for k=1:length(sidX)
    wdist(k, :) = W(:).*((resMVFx(:)-sidX(k)).^2+(resMVFy(:)-sidY(k)).^2);
end
[D_ALL, segMap] = min(wdist);

% Iterations
for ii=1:numIter
    % Update MVs
    for k=1:length(sidX)
        [idr, idc]=find(segMap == k);
        % we start k equal to 2, since the centroid of the first cluster 
        % will be fixed to (0, 0).
        if k > 1
            sidX(k) = sum(resMVFx(idc).*W(idc))/sum(W(idc));
            sidY(k) = sum(resMVFy(idc).*W(idc))/sum(W(idc));
        end
        D_REG(k)=sum(D_ALL(idc));
    end
    
    if ii == numIter
        break;
    end
    
    % update the quantization noise       
    for k=1:length(sidX)
        wdist(k, :) = W(:).*((resMVFx(:)-sidX(k)).^2+(resMVFy(:)-sidY(k)).^2);
    end
    [D_ALL, segMap] = min(wdist);
    
end
% segMap = b;

% plot(1:numIter+1, snr),
% hold on

