%
% t = .001;  % Distance threshold for deciding outliers
%
function [mm, iniH] = gmeLs(coorX, coorY, npx, npy, iniH, iMap, STEP, t)

% one step to further remove outliers according to last least square
% estimate

rcoorX = coorX + npx;
rcoorY = coorY + npy;

if ~isempty(iniH)

    

    % get the inliers from segmentation
    oriIn = find(iMap==1);
    N = length(iMap(:));
    
    if length(oriIn) > N*STEP
        x1 = [coorX(oriIn)'; coorY(oriIn)'; ones(1,length(oriIn))];
        x2 = [rcoorX(oriIn)'; rcoorY(oriIn)'; ones(1,length(oriIn))];  
    %     x1 = [coorX(:)'; coorY(:)'; ones(1,N)];
    %     x2 = [rcoorX(:)'; rcoorY(:)'; ones(1,N)];  
        [x1, T1] = normalise2dpts(x1);
        [x2, T2] = normalise2dpts(x2);

        % Calculate, in both directions, the transfered points    
        Hx1    = iniH*x1;
        invHx2 = pinv(iniH)*x2;

        % Normalise so that the homogeneous scale parameter for all coordinates
        % is 1.

        x1     = hnormalise(x1);
        x2     = hnormalise(x2);     
        Hx1    = hnormalise(Hx1);
        invHx2 = hnormalise(invHx2); 

        d2 = sum((x1-invHx2).^2)  + sum((x2-Hx1).^2);
        IXIN = find(abs(d2) < t);

        if length(IXIN) < N*STEP
            [B,IX] = sort(abs(d2), 'ascend');
            IXIN = IX(1:round(N*STEP));
        end

        inliers = oriIn(IXIN); %#ok<FNDSB> 
    else
        inliers = oriIn;
    end

%     [mcoorX, mcoorY] = coorEst(iniMM, coorX, coorY);  % estimate reference coordinate
%     emvx = npx-mcoorX+coorX;
%     emvy = npy-mcoorY+coorY;
%     emv = emvx.^2+emvy.^2;
% 
%     oriIn = find(iMap==1);
% 
%     N = length(emv);
% 
%     % for the first iteration, we only use pre-filtered MVs
%     in_emv = emv(oriIn);
%     [B,IX] = sort(in_emv);
% 
%     if length(oriIn) > STEP*N
%         IXIN = IX(1:round(N*STEP));
%         inliers = oriIn(IXIN);
%     else
%         inliers = oriIn;
%     end
else
    inliers = find(iMap(:)==1);
end
% 
% rcoorX = coorX + npx;
% rcoorY = coorY + npy;

[H, iniH] = homography_ls(rcoorX(:), rcoorY(:), coorX(:), coorY(:), inliers);

H=H';
ranM=H(:);
ranM = ranM./H(3,3);

mm = ranM(1:end-1)';
