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

% A description assuming 4 correspondences
% 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 ; 
%      0  0  0 x1 y1 1 ;
%      x2 y2 1 0  0  0 ;
%      0  0  0 x2 y2 1 ;
%      x3 y3 1 0  0  0 ;
%      0  0  0 x3 y3 1 ;
%      x4 y4 1 0  0  0 ;
%      0  0  0 x4 y4 1 ;
%      ];
%  
% B = [x1p y1p x2p y2p x3p y3p]';
%  

%populates A matrix and b vector
%A matrix and b vector are calculated by writing out the matrix
%multiplication of the original linear algebra problem, placing all
%unknowns in the x vector, and populating the corresponding A matrix and b
%vector which should be in terms of the known points.  
A = [];
b = [];
for index = 1:size(point,1)
    xCurrent = point(index,1);
    xpCurrent = pointp(index,1);
    yCurrent = point(index,2);
    ypCurrent = pointp(index,2);
    
    A = [A; 
         xCurrent yCurrent 1 0 0 0;
         0 0 0 xCurrent yCurrent 1];
     b = [b; xpCurrent; ypCurrent];
end

% least squares solution
 x = (A'*A)^-1 * A'*b;
%  x = A\B;   This is an alternative way to get the result.

 a = x(1);
 b = x(2);
 c = x(3);
 d = x(4);
 e = x(5);
 f = x(6);
 g = 0;
 h = 0;

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