Solve system of equations related to perspective projection

Question

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$ ?

Answer 1

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[
        Thread[{u[i], v[i]} == generaltrans[{x[i], y[i]}]]
        , {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]];

with added noise

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.