# Application of the projection matrix: Removing perspective distortion

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:

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 insigniﬁcant 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]}]]


1. 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 insigniﬁcant 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!!

• Have you seen ImageTransformation[], ImageForwardTransformation[], and ImagePerspectiveTransformation[]? (Maybe also FindGeometricTransform[] and TransformationFunction[].) Jul 30 '16 at 14:08
• @J.M. Yes I tried it, but I must have done a mistake somewhere, because it is not working... Jul 30 '16 at 14:12
• I am also getting imaginarry matrix coefficients which is a bit strange... Jul 30 '16 at 14:17
• Thank you, but I would like to try to solve it "manually" with those 4 point + assumption Jul 30 '16 at 15:10
• That's OK. Have you seen this? Jul 30 '16 at 15:18

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.

• thank you very much !! This is exactly what I was looking for! ...but how did you get this image ? Did you use my data ? ... because I am receiving a very distorded image. Jul 31 '16 at 6:51
• Also thank you very much for the additional comments !! Yes I have confused homogenous and inhomognous coordinates. I have erased this sentence from my question, but added additional imaformation so that one can much better understand the "physics". Jul 31 '16 at 6:56
• Hm, I've just re-entered the code from your question (starting at "Find two points witth should have the same x-coordinates:") and I get the same result image. Try quitting your kernel and re-evaluate everything. Sometimes old definitions are influencing calculations if symbols aren't reset using Clear Jul 31 '16 at 7:08
• Or maybe NSolve doesn't find the same solution every time? The output I get from NSolve is: {{h01 -> 0.608299 + 0.000282722 I, h02 -> -499.388 - 0.179458 I, h10 -> 0.0931617 - 0.0000789878 I, h11 -> 1.25351 + 0.000203753 I, h12 -> -325.443 - 0.0785295 I, h20 -> 0.000174956 - 4.79868*10^-8 I, h21 -> 0.000757305 + 4.26176*10^-7 I, h22 -> 0.237246 - 0.00024059 I}} Jul 31 '16 at 7:09
• Different NSolve results could be from using different versions. This issue has been raised a few times on MSE. Jul 31 '16 at 15:38

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]}


• Very nice ! Looks great ! Thanks a lot ! +1 Jan 23 '17 at 17:47