%returns the H matrix given point and pointp (transformed coordinates)
function out = findHMatrix (point, pointp)

% x1 = point1(1);
% x2 = point2(1);
% x3 = point3(1);
% x4 = point4(1);
% y1 = point1(2);
% y2 = point2(2);
% y3 = point3(2);
% y4 = point4(2);
% 
% x1p = point1p(1);
% x2p = point2p(1);
% x3p = point3p(1);
% x4p = point4p(1);
% y1p = point1p(2);
% y2p = point2p(2);
% y3p = point3p(2);
% y4p = point4p(2);
% 
% A = [x1 y1 1 0  0  0  -x1*x1p -y1*x1p; 
%      0  0  0 x1 y1 1  -x1*y1p -y1*y1p;
%      x2 y2 1 0  0  0  -x2*x2p -y2*x2p;
%      0  0  0 x2 y2 1  -x2*y2p -y2*y2p;
%      x3 y3 1 0  0  0  -x3*x3p -y3*x3p;
%      0  0  0 x3 y3 1  -x3*y3p -y3*y3p;
%      x4 y4 1 0  0  0  -x4*x4p -y4*x4p;
%      0  0  0 x4 y4 1  -x4*y4p -y4*y4p;
%      ];
%  
% B = [x1p y1p x2p y2p x3p y3p x4p y4p]';
%  
A = [];
b = [];
for index = 1:size(point,1)
    xCurrent = point(index,1);
    xpCurrent = pointp(index,1);
    yCurrent = point(index,2);
    ypCurrent = pointp(index,2);
    
    %calculates A matrix and B vector
    A = [A; 
         xCurrent yCurrent 1 0 0 0 -xCurrent * xpCurrent -yCurrent * xpCurrent;
         0 0 0 xCurrent yCurrent 1  -xCurrent * ypCurrent -yCurrent * ypCurrent;];
     b = [b; xpCurrent; ypCurrent];
end

%least squares solution
 x = (A'*A)^-1 * A'*b;
%  x = A\B;  This is an alternative way to obtain the solution.
 
 a = x(1);
 b = x(2);
 c = x(3);
 d = x(4);
 e = x(5);
 f = x(6);
 g = x(7);
 h = x(8);

 H = [a b c;
      d e f;
      g h 1];
 out = H; 

end