I am trying to implement a perspective projection. For example, I have a square with a given size, and known corners' coordinates (A1 .. A4) . I have an image where this square is displayed at some specific position. I know its position and the coordinates of the corners of the displayed square B1 .. B4. I want to find a transformation matrix M that will
translate B -> A.
I followed method Mathematical Illustrations, chapter 10, and I managed to transform 3 corners well. But I cannot find the matrix to transform correctly 4 corners... I find the matrix as on page 7, but it does not translate all 4 points well.
~
As your reference indicates, the standard way to represent a general projective transformation is to augment your (x,y) 2-D coordinates with a w coordinate (to get homogeneous coordinates [x,y,w]), and use a 3x3 matrix to transform them.
It is not clear from your question exactly how you are currently trying to use your points, but I suspect that there is some confusion with how to use the extra coordinate. Mathematically, the thing to remember is that you can multiply the entire homogeneous vector by an arbitrary nonzero scaling factor (because in the end dividing by the third coordinate will cancel out the scaling factor). However, it is sometimes hard to figure out what this means in a practical sense...
For this problem, set up the system as follows (I am treating homogeneous coordinates as row vectors here, to match your reference):
Find a 3x3 matrix T, such that:
it maps each point b to corresponding point a: a = b T
specifically:
w' * [a_x a_y 1] = w * [b_x b_y 1] * [ T_xx T_xy T_xw ]
[ T_yx T_yy T_yw ]
[ T_wx T_wy T_ww ]
The tricky bit: you can't just set w and w' to 1. This is because a projective transformation T doesn't necessarily leave the third coordinate unchanged! (If it did, you could skip the homogeneous bit altogether, and get an affine transform using only three point pairs...).
One way to express this (in a way that can be solved readily) is to add a parameter, K = w' / w. Then you can express the above as a linear system, as follows:
for each point b and corresponding point a,
[a_x a_y 1] * K = [b_x b_y 1] * [ T_xx T_xy T_xw ]
[ T_yx T_yy T_yw ]
[ T_wx T_wy T_ww ]
Remember, all the a and b values are fixed constants; you're trying to solve for the nine elements of your T matrix. Each point pair adds three constraints (from the vector equality above) and one additional parameter (the K for that equation), for a total of 4*3=12 constraints, and 9+4=13 parameters. So, we're missing one constraint here...
An additional constraint is needed because the matrix T effectively has an arbitrary scaling factor: it is mapping between homogeneous coordinates (which all divide out an arbitrary scaling factor anyway), so there is no inherent way to decide what this scaling factor is. One way to fix this factor is to arbitrarily set T_ww = 1:
for each point pair (a, b):
[a_x a_y 1] * K = [b_x b_y 1] * [ T_xx T_xy T_xw ]
[ T_yx T_yy T_yw ]
[ T_wx T_wy 1 ]
This gives you a fully-determined linear system:
0 = b1_x*T_xx + b1_y*T_yx + 1*T_wx - a1_x*K1
0 = b1_x*T_xy + b1_y*T_yy + 1*T_wy - a1_y*K1
-1*1 = b1_x*T_xw + b1_y*T_yw - 1*K1
0 = b2_x*T_xx + b2_y*T_yx + 1*T_wx - a2_x*K2
0 = b2_x*T_xy + b2_y*T_yy + 1*T_wy - a2_y*K2
-1*1 = b2_x*T_xw + b2_y*T_yw - 1*K2
0 = b3_x*T_xx + b3_y*T_yx + 1*T_wx - a3_x*K3
0 = b3_x*T_xy + b3_y*T_yy + 1*T_wy - a3_y*K3
-1*1 = b3_x*T_xw + b3_y*T_yw - 1*K3
0 = b4_x*T_xx + b4_y*T_yx + 1*T_wx - a4_x*K4
0 = b4_x*T_xy + b4_y*T_yy + 1*T_wy - a4_y*K4
-1*1 = b4_x*T_xw + b4_y*T_yw - 1*K4
You now have 12 linear equations in 12 variables (matrix elements T_** and additional scaling parameters K*). Use a linear algebra library to solve the system and get your results (assuming you don't feel like writing your own solver...).
A better answer that solves an 8x8 system is given in the math stack exchange. https://mathematica.stackexchange.com/questions/9244/solve-system-of-equations-related-to-perspective-projection
