Application of the projection matrix: Removing perspective distortion

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.

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