Computing a matrix which transforms a quadrangle to another quadrangle in 2D

Question

In the figure below the goal is to compute the homography matrix H which transforms the points a1 a2 a3 a4 to their counterparts b1 b2 b3 b4. That is:

[b1 b2 b3 b4] = H * [a1 a2 a3 a4]

What way would you suggest to be the best way to calculate H(3x3). a1...b4 are points in 2D which are represented in homogeneous coordinate systems ( that is [a1_x a1_y 1]', ...). EDIT: For these types of problems we use SVD, So i would like to see how this can be simply done in Matlab.

EDIT:

Here is how I initially tried to solve it using svd (H=Q/P) in Maltlab. Cosider the following code for the given example

px=[0 1 1 0];  % a square
py=[1 1 0 0];

qx=[18 18 80 80];    % a random quadrangle
qy=[-20 20 60 -60];
if (DEBUG)
  fill(px,py,'r');
  fill(qx,qy,'r');
end

Q=[qx;qy;ones(size(qx))];
P=[px;py;ones(size(px))];
H=Q/P;
H*P-Q
answer:
   -0.0000         0         0         0         0
  -20.0000   20.0000  -20.0000   20.0000    0.0000
   -0.0000         0         0         0   -0.0000

I am expecting the answer to be a null matrix but it is not!... and that's why I asked this question in StackOverflow. Now, we all know it is a projective transformation not obviously Euclidean. However, it is good to know if in general care calculating such matrix using only 4 points is possible.

Answer 1

You can try the cp2tform function, which infers spatial transformation from control point pairs. Since parallel is not preserved in your case, you should set the transformtype to be 'projective'. More info is here