Application of the projection matrix: Removing perspective distortion

Question

I know that there are many posts about image treatment, but I still have trouble with this little challegene. Baisically I want to remove the perspective distorsion form this image: (source: http://www.canvas-of-light.com/2015/03/how-to-correct-perspective-distortion-in-lightroom-5/)

How I started: An image can be transformed in 2D (planar) using various transformation rules:

Source: https://www.microsoft.com/en-us/research/wp-content/uploads/2004/10/tr-2004-92.pdf

Those transformation can all be described by transformation matrices which map the untransformed image to the transformed one. The most general case (for planar) is the projective transformation.

The projective matrix:

ProjectMat2D[h00_, h01_, h02_, h10_, h11_, h12_, h20_, h21_, 
   h22_] := {{h00, h01, h02}, {h10, h11, h12}, {h20, h21, h22}};
ProjectMat2D[h00, h01, h02, h10, h11, h12, h20, h21, 
  h22] // MatrixForm

The projective matrix has 3x3 components, since it works with homogeneous coordinates. Homogenous coordinates are basically the normale coordinates with an added row of 1. {{x},{y},{1}} ... for more information: https://en.wikipedia.org/wiki/Homogeneous_coordinates

Finding the matrix coeffiecnts:

Now the task is to find the coefficient of the transformation matrix.

Let's first calculate the transformed xy-Vector in homognenous coordinates.

xyTransHomoge = 
 ProjectMat2D[h00, h01, h02, h10, h11, h12, h20, h21, 
   h22].{{x}, {y}, {1}}; xyTransHomoge // MatrixForm

To retreive the homogenous coordinates as output, the last row needs to be 1. Therefore we devide the vector by the value of the last entry (to make it 1)

Last[xyTransHomoge]

{h22 + h20 x + h21 y}

1/(h22 + h20 x + h21 y)*xyTransHomoge // MatrixForm

We drop the {1} entry and receive the projection formulas which map the transformed to the untransformed coordinates:

$x'=\frac{h00*x+h01*y+h02}{h20*x+h21*y+h22}$

$y'=\frac{h10*x+h11*y+h12}{h20*x+h21*y+h22}$

I find four points in the transformed image and the corresponding points in the untransfromed image (by assumption). Basically I say that the four points on the transformed image should be a rectangle on the untransformed image. -> same x Position for two points on a vertical line

We get set of 8 equations, which we can solve to find the matrix coeffiecnts.

We set h00=1, since the matrix needs only to be determined up to a insignificant multiplicative factor.

Okay, enough theory... let's do it:

1. Find two points witth should have the same x-coordinates:

   x11 = 634.0260655160307;
   y11 = 672.1241170908934;
   x12 = 608.4425576092435;
   y12 = 523.0723753730896;
   Show[notreDame, 
    Graphics[{Point[{x11, y11}, VertexColors -> Green], 
      Point[{x12, y12}, VertexColors -> Green]}]]
   

2. Find two other points witth should have the same x-coordinates:

   x21 = 926.9209508015505;
   y21 = 611.3228302929783;
   x22 = 938.2421227197373;
   y22 = 516.5080154781632;
   Show[notreDame, 
    Graphics[{Point[{x21, y21}, VertexColors -> Green], 
      Point[{x22, y22}, VertexColors -> Green]}]]
   

To find the corresponding homogenious coordinates (in undeformed image) let's choose the common x-Value in both sets of points. The inhomognious y-values are considered to be the same as the homogenious y-values (apporx.).

x11h = x11;
y11h = y11;
x12h = x11;
y12h = y12;
x21h = x21;
y21h = y21;
x22h = x21;
y22h = y22;

Now we get a system of equation that we can solve for the matrix coefficient. We set h00=1, since the matrix needs only to be determined up to a insignificant multiplicative factor.

h00 = 1;
coef = NSolve[{
   x11 == (h00*x11h + h01*y11h + h02)/(h20*x11h + h21*y11h + h22), 
   y11 == (h10*x11h + h11*y11h + h12)/(h20*x11h + h21*y11h + h22),
   x12 == (h00*x12h + h01*y12h + h02)/(h20*x12h + h21*y12h + h22), 
   y12 == (h10*x12h + h11*y12h + h12)/(h20*x12h + h21*y12h + h22),
   x21 == (h00*x21h + h01*y21h + h02)/(h20*x21h + h21*y21h + h22),
   y21 == (h10*x21h + h11*y21h + h12)/(h20*x21h + h21*y21h + h22),
   x22 == (h00*x22h + h01*y22h + h02)/(h20*x22h + h21*y22h + h22), 
   y22 == (h10*x22h + h11*y22h + h12)/(h20*x22h + h21*y22h + h22)
   },
  {h01, h02, h10, h11, h12, h20, h21, h22}
  ]

projeMat = 
  Flatten[{{h00, h01, h02}, {h10, h11, h12}, {h20, h21, h22}} /. coef,
    1];

The output has imaginary values ??

To remove the perspective distortion, we need to inverse the transformation matrix and apply it to the image:

invProjeMat = Inverse[projeMat];

ImageTransformation[notreDame, TransformationFunction[invProjeMat]]

...but this is not working?

Thanks for any help!!

Answer 1

You're very close: First, ImageTransformation by default assumes that

the range of the coordinate system for the input image is [...] {{0,1},{0,a}}, where a is the aspect ratio

If you want to work with pixel coordinates, you have to add PlotRange->Full.

Second, the transformation passed to ImageTransformation should transform coordinates from the transformed image to the source image. The small imaginary values you see are inaccuracies due to NSolve - just use the real part, and you get the image you'd expect:

ImageTransformation[notreDame, TransformationFunction[Re@projeMat],
    PlotRange -> Full]

Now, a few further tips:

• You can save yourself a lot of typing if you use arrays instead of separate variables. So e.g. in instead of variables x11, y11, x12, y12,..., you'd use one variable sourceImagePoints = {{624,672},{608,623},...} you could use Dot, Map, Outer and so on to get much simpler code.

• Your system of equations is linear, you don't need NSolve. Solve does just fine, but you have to help it a little, by multiplying the denominators to the left side. You can do this by adding /. (x_ == a_/b_) :> (x*b == a) after your equations:

  Solve[{...)} /. (x_ == a_/b_) :> (x*b == a), {...}]

• I'm pretty sure you're using the terms "homogeneous" and "inhomogeneous" wrong. They don't mean "coordinates in the transformed/untransformed images", respectively. The word homogeneous means (more or less) that if you multiply a vector with a scalar (except 0), that vector refers to the same point in space. It's very elegant mathematical trick that lets you treat translations and projective transformations as simple matrix multiplications. You usually use homogeneous coordinates for both source image coordinates and destination image coordinates, and only convert to euclidean coordinates at the end, when you call some graphics routine.

Answer 2

The workflow can be made easier using the FindGeometricTransform function.

First we find 4 points on the original image that defines a rectangle in the undistorted image. This can be done easily using the "Get Coordinates" in the tool.

pts1 = {{616.597, 569.06}, {634.702, 685.615}, {916.721, 681.76}, {930.127, 566.84}};

We want to map this rectangle in the distorted image to a true rectangle in the undistorted image. To define a rectangle in the undistorted image, we can just use the rectangle defined by the two corners for convenience

pts2 = {{pts1[[1, 1]], pts1[[1, 2]]}, {pts1[[1, 1]], 
   pts1[[3, 2]]}, {pts1[[3, 1]], pts1[[3, 2]]}, {pts1[[3, 1]], 
   pts1[[1, 2]]}}

This is what the two rectangles look like

Graphics[{PointSize[Large], Red, Point[pts1], Green, Point[pts2], 
  Transparent, EdgeForm[Red], Polygon[pts1], EdgeForm[Green], 
  Polygon[pts2]}]

Our task is to transform the red rectangle to the green rectangle. The same transformation will take us from the distorted image the undistorted one.

The transformation can be found by

t = FindGeometricTransform[pts2, pts1][[2]]

We see that it indeed does the job

Graphics[{PointSize[Large], Red, Point[t@pts1], Green, Point[pts2], 
  Transparent, EdgeForm[Red], Polygon[t@pts1], EdgeForm[Green], 
  Polygon[pts2]}]

Now we can undistort the image using the same transformation

Row@{Show[
   norteDame,
   Graphics[{PointSize[Large], Red, Point[pts1], Transparent, 
     EdgeForm[Red], Polygon[pts1]}], ImageSize -> Large],
  Show[
   ImagePerspectiveTransformation[norteDame, t, PlotRange -> Full],
   Graphics[{PointSize[Large], Green, Point[pts2], Transparent, 
     EdgeForm[Green], Polygon[pts2]}], ImageSize -> Large]}