Get Transformation Matrix from Points I have built a little C# application that allows visualization of perpective transformations with a matrix, in 2D XYW space. Now I would like to be able to calculate the matrix from the four corners of the transformed square. Here is a screenshot to help you understand what I am working with:

The idea is to allow the user to move the corners of the square and update the matrix. The corners are currently located at (1,1), (1,-1), (-1,-1), (-1,1). Is there an algorithm that will calculate a 3x3 matrix given four 2D points?
If I understand this correctly, every matrix corresponds to a set of four points, and every set of four points corresponds to one or more equivalent matrices ('equivalent' meaning 'producing identical transformations').
I searched for an algorithm to do this, but didn't have much luck.
I figured out that I could do it by creating eight equations, one for each variable in the four points, and then setting one of the matrix values to one, and solving for the other eight with algebra. However, the equations grow much too complicated to do this all successfully on pencil and paper.

This is the process I used to try to get it working.
So this is the basic matrix transformation formula.
$\begin{pmatrix}a & b & c\\
d & e & f\\
g & h & i
\end{pmatrix}\begin{pmatrix}x\\
y\\
z
\end{pmatrix}=\begin{pmatrix}ax+by+cz\\
dx+ey+fz\\
gx+hy+iz
\end{pmatrix}$
The resulting point is then converted from homogeneous to Euclidean coordinates.
$\begin{pmatrix}x\\
y\\
z
\end{pmatrix}$=$\begin{pmatrix}\frac{x}{z}\\
\frac{y}{z}
\end{pmatrix}$
So given a collection of points, we transform them like this.
$\begin{pmatrix}a & b & c\\
d & e & f\\
g & h & i
\end{pmatrix}\begin{pmatrix}x_{n}\\
y_{n}\\
z_{n}
\end{pmatrix}=\begin{pmatrix}x'_{n}\\
y'_{n}
\end{pmatrix}$
These are the formulas used for the transformation.
$x'_{n}=$$\frac{ax_{n}+by_{n}+cz_{n}}{gx_{n}+hy_{n}+iz_{n}}$
$y'_{n}=$$\frac{dx_{n}+ey_{n}+fz_{n}}{gx_{n}+hy_{n}+iz_{n}}$
We then define four base points, which are the corners of our square.
$x'_{0}=1$
$y'_{0}=1$
$z'_{0}=1$
$x'_{1}=1$
$y'_{1}=-1$
$z'_{1}=1$
$x'_{2}=-1$
$y'_{2}=-1$
$z'_{2}=1$
$x'_{3}=-1$
$y'_{3}=1$
$z'_{3}=1$
This gives us the following system of equations, in the form that lets us determine the transformed points from the matrix.
$x'_{0}=$$\frac{a+b+c}{g+h+i}$
$y'_{0}=$$\frac{d+e+f}{g+h+i}$
$x'_{1}=$$\frac{a-b+c}{g-h+i}$
$y'_{1}=$$\frac{d-e+f}{g-h+i}$
$x'_{2}=$$\frac{-a-b+c}{-g-h+i}$
$y'_{2}=$$\frac{-d-e+f}{-g-h+i}$
$x'_{3}=$$\frac{-a+b+c}{-g+h+i}$
$y'_{3}=$$\frac{-d+e+f}{-g+h+i}$
Now we want to reverse the transformation, and find the matrix that produces the above points. Since we have 9 unknowns and 8 equations, we need to add another equation.
$i=1$
Now all that is left is to solve the system equations to find the formulas for the matrix values. I'm not patient enough nor good enough with algebra to do this myself, so I used an online calculator to solve the system of equations. The formulas it gave almost worked, but had some glitches with y-coordinates.
I think this can be narrowed down to 2 questions:


*

*Are the above calculations wrong, or does the online calculator have a bug?

*Is there an easier algorithm for this?
 A: Let T be the unknown transformation matrix (3X3).
Let A be the matrix of original co-ordinates (4X3).
Let B be the matrix of new co-ordinates (4X3).
Now A*T = B. Solve for 9 unknowns using 12 equations to get the transformation. As you can see you need only 3 points, not 4 to fully define the transformation.
A: I'll use $f(x, y)$ to represent the 2D point $(x, y)$ multiplied by an assumed 3x3 transformation matrix:
$\begin{pmatrix}
x'\\
y'\\
z'
\end{pmatrix}
:=
\begin{pmatrix}
a & b & c\\
d & e & f\\
g & h & i
\end{pmatrix}
\begin{pmatrix}
x\\
y\\
1
\end{pmatrix}$
$f(x, y) :=
\begin{pmatrix}
\frac{x'}{z'}\\
\frac{y'}{z'}
\end{pmatrix}
=
\begin{pmatrix}
\frac{ax + by + c}{gx + hy + i}\\
\frac{dx + ey + f}{gx + hy + i}
\end{pmatrix}$
Given the transformation matrix, the corners of a square centered on $(0, 0)$ can be transformed:
$f(1, 1) = (x_1, y_1)\\
f(1, -1) = (x_2, y_2)\\
f(-1, -1) = (x_3, y_3)\\
f(-1, 1) = (x_4, y_4)$
Forgive me for using slightly different naming and notation than the original question. It was a long time ago.
The problem is to find a transformation matrix (values for $a$ through $i$) given these 8 output variables. I used the following equations to compute a solution.
First, some temporary variables to simplify the process:
$j = x_1 - x_2 - x_3 + x_4\\
k = -x_1 - x_2 + x_3 + x_4\\
l = -x_1 + x_2 - x_3 + x_4\\
m = y_1 - y_2 - y_3 + y_4\\
n = -y_1 - y_2 + y_3 + y_4\\
o = -y_1 + y_2 - y_3 + y_4$
Then values for the matrix elements, calculated in reverse order:
$i = 1\\
h = \frac{(jo - ml)i}{mk - jn}\\
g = \frac{kh + li}{j}\\
f = \frac{y_1(g + h + i) + y_3(-g - h + i)}{2}\\
e = \frac{y_1(g + h + i) - y_2(g - h + i)}{2}\\
d = y_1(g + h + i) - f - e\\
c = \frac{x_1(g + h + i) + x_3(-g - h + i)}{2}\\
b = \frac{x_1(g + h + i) - x_2(g - h + i)}{2}\\
a = x_1(g + h + i) - c - b$
I have a proof of concept implementation: https://jsfiddle.net/kendfrey/oxm5L6q0/5/
For the unhealthily curious, my working out was done in Lean: https://gist.github.com/kendfrey/28c49ebc1c28e543ca6322e167fb8fd8
Notes:

*

*Since there are 9 unknown variables and 8 equations, I'm free to choose a value for $i$.

*This works in the general case but has some degenerate cases where division by zero occurs.

A: First lets the the naming straight. 2D means two-dimensional so this space only needs two axis - X and Y. For transformations in this space you only need a two dimensional matrix, lets call it T. To find T you only need two points (in fact you find out one matrix for two points and then you can check if it also works for the other points you have).
Now lets say you have the two points a and b and their transformed images a' and b'.

From this by matrix multiplication you get four equations:

These are made by simple multiplying T x a and T x b. From these it is fairly easy to derive the value of each number in the T matrix. I started by extracting t11 and t21 from the first two equations and moved on to the next. Here is the result:

This is almost enough for your application to work. The last thing you need is to handle a special case when a1 is zero. In which case some of these formulas do not work. So you have to start from the equations and derive to this:

All of this is enough for simple transformations. You will need to extend this with one more dimension.
A: When the image of a square in space can be an arbitrary quadrilateral (a non-parallelogram) then the map in question cannot be linear and is not obtained by applying a matrix to the coordinates of the original points in space. Instead the mapping functions are of the form
$$\xi={a_1 x+ b_1 y+ c_1\over a_2 x+ b_2 y + c_2}\ ,\qquad \eta={a_3 x+ b_3 y+ c_3\over a_4 x+ b_4 y + c_4}\ .$$
Here $(\xi,\eta)$ are coordinates on the image screen, and $(x,y)$ are coordinates in the plane in which the original square is situated. The coefficients $a_1$, $\ldots$, $c_4$ appearing in these formulae depend on the various choices made: Location of the origins in the original and the image plane, adopted coordinate systems in these planes, position of these planes vs. the center $O$ of the used perspective.
Maybe the following source is of some help: http://www.geometrictools.com/Documentation/PerspectiveMappings.pdf
