%x1=[x1,y1] vector of the points in the first image
%x1=[x1,y1] vector of the corresponding points in the second image
A = zeros(numpoints*2, 8);
b = zeros(numpoints*2,1);
for i=1:numpoints
    A(i*2-1, :) = [x1(i), y1(i), 1, 0, 0, 0, -x1(i)*x2(i), -y1(i)*x2(i)];
    A(i*2, :)   = [0, 0, 0, x1(i), y1(i), 1, -x1(i)*y2(i), -y1(i)*y2(i)];
    b(i*2-1) = x2(i);
    b(i*2)   = y2(i);
end

h = inv(A'*A)*A'*b;
H = [h(1), h(2), h(3);
     h(4), h(5), h(6);
     h(7), h(8), 1 ];