# Solve system of equations related to perspective projection

I'm trying to solve system of equations. I have 4 points before perspective projection and after and I whant to obtain matrix of perspective projection $m_{i,j}$ ($x,y,u,v$ are known , $m_{i,j}$ unknown)

$u_{1}= (m_{13} + m_{11}*x_{1} + m_{12}*y_{1})/(1 + m_{31}*x_{1} + m_{32}*y_{1})$ $v_{1}= (m_{23} + m_{21}*x_{1} + m_{22}*y_{1})/(1 + m_{31}*x_{1} + m_{32}*y_{1})$ $u_{2}= (m_{13} + m_{11}*x_{2} + m_{12}*y_{2})/(1 + m_{31}*x_{2} + m_{32}*y_{2})$ $v_{2}= (m_{23} + m_{21}*x_{2} + m_{22}*y_{2})/(1 + m_{31}*x_{2} + m_{32}*y_{2})$ $u_{3}= (m_{13} + m_{11}*x_{3} + m_{12}*y_{3})/(1 + m_{31}*x_{3} + m_{32}*y_{3})$ $v_{3}= (m_{23} + m_{21}*x_{3} + m_{22}*y_{3})/(1 + m_{31}*x_{3} + m_{32}*y_{3})$ $u_{4}= (m_{13} + m_{11}*x_{4} + m_{12}*y_{4})/(1 + m_{31}*x_{4} + m_{32}*y_{4})$ $v_{4}= (m_{23} + m_{21}*x_{4} + m_{22}*y_{4})/(1 + m_{31}*x_{4} + m_{32}*y_{4})$

I tried Solve but it seems to work infinite time.

my input in mathematica

Solve[u1== (m13 + m11*x1 + m12*y1)/(1 + m31*x1 + m32*y1)&&
v1== (m23 + m21*x1 + m22*y1)/(1 + m31*x1 + m32*y1)&&
u2== (m13 + m11*x2 + m12*y2)/(1 + m31*x2 + m32*y2)&&
v2== (m23 + m21*x2 + m22*y2)/(1 + m31*x2 + m32*y2)&&
u3== (m13 + m11*x3 + m12*y3)/(1 + m31*x3 + m32*y3)&&
v3== (m23 + m21*x3 + m22*y3)/(1 + m31*x3 + m32*y3)&&
u4== (m13 + m11*x4 + m12*y4)/(1 + m31*x4 + m32*y4)&&
v4== (m23 + m21*x4 + m22*y4)/(1 + m31*x4 + m32*y4),{m11,m12,m13,m21,m22,m23,m31,m32}]


UPDATE:

another question rises when I try to solve the same equation but when $x,y$ not given, but I have some restrictions that give me more equations. What I trying to do: I have points $u,v$ (distorted points) and I don't have original points $x,y$ but I can foto the same segment in different places and I know that distance $dx=x_{2}-x_{1}$ and $dy=y_{2}-y_{1}$ is always the same,so I want to reconstruct perspective transform that I don't know.

my input in mathematica

Solve[(m31 x1 + m32 y1 + 1) u1 == (m11 x1 + m12 y1 + m13) &&
(m31 x1 + m32 y1 + 1) v1 == (m21 x1 + m22 y1 + m23) &&
(m31 x2 + m32 y2 + 1) u2 == (m11 x2 + m12 y2 + m13) &&
(m31 x2 + m32 y2 + 1) v2 == (m21 x2 + m22 y2 + m23) &&
(m31 x3 + m32 y3 + 1) u3 == (m11 x3 + m12 y3 + m13) &&
(m31 x3 + m32 y3 + 1) v3 == (m21 x3 + m22 y3 + m23) &&
(m31 x4 + m32 y4 + 1) u4 == (m11 x4 + m12 y4 + m13) &&
(m31 x4 + m32 y4 + 1) v4 == (m21 x4 + m22 y4 + m23) &&
(m31 x5 + m32 y5 + 1) u5 == (m11 x5 + m12 y5 + m13) &&
(m31 x5 + m32 y5 + 1) v5 == (m21 x5 + m22 y5 + m23) &&
(m31 x6 + m32 y6 + 1) u6 == (m11 x6 + m12 y6 + m13) &&
(m31 x6 + m32 y6 + 1) v6 == (m21 x6 + m22 y6 + m23) &&
(m31 x7 + m32 y7 + 1) u7 == (m11 x7 + m12 y7 + m13) &&
(m31 x7 + m32 y7 + 1) v7 == (m21 x7 + m22 y7 + m23) &&
(m31 x8 + m32 y8 + 1) u8 == (m11 x8 + m12 y8 + m13) &&
(m31 x8 + m32 y8 + 1) v8 == (m21 x8 + m22 y8 + m23) &&
(m31 x9 + m32 y9 + 1) u9 == (m11 x9 + m12 y9 + m13) &&
(m31 x9 + m32 y9 + 1) v9 == (m21 x9 + m22 y9 + m23) &&
(m31 x10 + m32 y10 + 1) u10 == (m11 x10 + m12 y10 + m13) &&
(m31 x10 + m32 y10 + 1) v10 == (m21 x10 + m22 y10 + m23) &&
dx == x2 - x1 &&
dy == y2 - y1 &&
dx == x4 - x3 &&
dy == y4 - y3 &&
dx == x6 - x5 &&
dy == y6 - y5 &&
dx == x8 - x7 &&
dy == y8 - y7 &&
dx == x10 - x9 &&
dy == y10 - y9
, {m11, m12, m13, m21, m22, m23, m31, m32}]


but Mathematica give me {} just trivial solution.what is wrong?

maybe it's possible to automatically rewrite system of equations as linear system $AX=B$ ?

-
@Ajasja add mathematica code. –  mrgloom Aug 9 '12 at 8:07
What happened to m33 ? - got eliminated from 9th equation? –  Vitaliy Kaurov Aug 9 '12 at 8:09
@Vitaliy Kaurov m33=1 for perspective transform. –  mrgloom Aug 9 '12 at 8:26
regarding your second question: you basically have an overdetermined system, because you have more measuring points than unknowns. You can use PseudoInverse or LeastSquares for that. I'll try to add an example in my answer. –  Thies Heidecke Aug 9 '12 at 12:16
@mrgloom also you should put x1,x2,... in the list of unknowns for Solve then. For measured values which can have measurement errors or inaccuracies it's usually better to use "soft" methods like PseudoInverse or LeastSquares which also work when contradictions arise. –  Thies Heidecke Aug 9 '12 at 12:32

Question 1

stating the problem as a system of linear equations

vars = {m11, m12, m13, m21, m22, m23, m31, m32};
Solve[
u1 (1 + m31*x1 + m32*y1) == (m13 + m11*x1 + m12*y1) &&
v1 (1 + m31*x1 + m32*y1) == (m23 + m21*x1 + m22*y1) &&
u2 (1 + m31*x2 + m32*y2) == (m13 + m11*x2 + m12*y2) &&
v2 (1 + m31*x2 + m32*y2) == (m23 + m21*x2 + m22*y2) &&
u3 (1 + m31*x3 + m32*y3) == (m13 + m11*x3 + m12*y3) &&
v3 (1 + m31*x3 + m32*y3) == (m23 + m21*x3 + m22*y3) &&
u4 (1 + m31*x4 + m32*y4) == (m13 + m11*x4 + m12*y4) &&
v4 (1 + m31*x4 + m32*y4) == (m23 + m21*x4 + m22*y4),
vars
]


seems to work. On a sidenote, it might be interesting to know, that for doing projective transformations Mathematica supplies the LinearFractionalTransform function. So you could state your transformation like this:

generaltrans = LinearFractionalTransform[
{{m11, m12, m13},
{m21, m22, m23},
{m31, m32, 1}}
]


, get the system of equations

eqs = And @@ Flatten @ Table[
, {i, 4}
] /. (a_ == p_/q_ -> a*q == p)


and solve it

Solve[eqs, vars]


.

Example

Regarding your second question i want to show a more practically relevant example that should give you a general idea on how you could approach problems like this. Let's say we have a bunch of points in a regular cartesian grid

pointsxy = Flatten[Table[{i, j}, {j, 5}, {i, 5}], 1];
Graphics@Point[%]


that got transformed by an unknown projective transformation

unknowntransformation = generaltrans /. Thread[vars -> RandomReal[{1, 3}, 8]];


AddNoise = # + RandomVariate[NormalDistribution[0, 0.001], 2] &;


resulting in something like this

pointsuv = Composition[AddNoise, unknowntransformation] /@ pointsxy;
Graphics@Point[%]


and now we ask the question: What does the transformation look like, that did this projection?

For answering this question our plan is the following: First we create the system of linear equations that determine the projection matrix elements that transform our specific points from xy- to uv-space.

MakeProjectionEquations[ptsxy_List, ptsuv_List, m_] := Module[
{t, n},
t = LinearFractionalTransform[m];
n = Min[Length /@ {ptsxy, ptsuv}];
Flatten @ Table[Thread[ptsuv[[i]] == t[ptsxy[[i]]]], {i, n}] /. (a_ == p_/q_ -> a*q == p)
]


Now we want a way to guess the most probable transformation matrix given the points before and after the transformation. Since we might measure more than the 4 points necessary to uniquely determine the transform, we'll do a general least squares fit that determines the matrix elements which give us the best result averaged over all points.

GuessProjection[ptsxy_List, ptsuv_List] := Module[
{m11, m12, m13, m21, m22, m23, m31, m32, eqs, b, A, m},
(* Get our system of linear equations for projection *)
eqs = MakeProjectionEquations[ptsxy, ptsuv, {{m11, m12, m13}, {m21, m22, m23}, {m31, m32, 1}}];
(* extract the matrix and constant vector fulfilling: A.m + b==0, m={m11,m12,...} *)
{b, A} = CoefficientArrays[eqs, {m11, m12, m13, m21, m22, m23, m31, m32}];
(* find a best guess for the values of the projection matrix *)
m = LeastSquares[A, -b];
LinearFractionalTransform[Partition[Append[m, 1], 3]]
]


Now we can guess both the forward and inverse transformation

GuessProjection[pointsxy, pointsuv]
GuessProjection[pointsuv, pointsxy]


and use it to recover the original grid from the measured points (or other points that we might be interested in).

recoveredpointsxy = GuessProjection[pointsuv, pointsxy] /@ pointsuv;
Graphics@Point[%]


Question 2

For your specific problem, where you have the same points in xy-space seen from different perspectives, you can take a similar approach, but in your case the unknowns are not the transformation matrix elements but the coordinates of the points in xy-space, so you would have to adapt the code accordingly (you could take GuessProjection as a guideline). Also instead of giving extra constraints in terms of differences dx, dy, ... you could arbitrarily designate one point as the (0, 0)-point in xy-space and let the LeastSquares-routine figure out the rest.

-
yes it works, but can you explain why mathematica can't get rid of the division by default on it's own? –  mrgloom Aug 9 '12 at 8:34
Because of the divisions, Mathematica treats the problem as general nonlinear, which makes it a lot harder, whereas for systems of linear equations it can use standard linear algebra methods for solving. So whenever you know the problem can be stated as a linear problem, that is usually solved more easily. –  Thies Heidecke Aug 9 '12 at 8:38
@mrgloom Mathematica doesn't get rid of the division because in general it would be plain wrong! a==p/q doesn't imply a*q==p because q could be zero: Resolve@ForAll[{p, q}, Implies[a == p/q, a*q == p]]. Only if we know that q is nonzero that would be alright to do: Resolve@ForAll[{p, q}, q != 0, Implies[a == p/q, a*q == p]] –  Thies Heidecke Aug 9 '12 at 11:56
@mrgloom That's true. Maybe Mathematica might express a solution as a ConditionalExpression then. But it might need a whole lotta time to figure out all the cases and their conditions. That's why it's in general more advisable to put your problem in the easiest form first before letting Mathematica solve it. –  Thies Heidecke Aug 9 '12 at 12:24
@mrgloom Regarding your question if it's possible to automatically rewrite a set of equations into a matrix and constant vector describing a linear system, have a look at the CoefficientArrays function in GuessProjection. It does exactly that tedious rewriting stuff. –  Thies Heidecke Aug 10 '12 at 7:54