Deskew and rotate a photographed rectangular image (aka "perspective correction")

Question

I photographed a rectangular shape, but:

• the camera angle is not perfect
• the image is possibly rotated.

Given the coordinates of 8 points $ A_1 (x_1, y_1), ..., A_8 (x_8, y_8)$, how to transform the input coordinates into perfectly rectangular output coordinates?

Note:

1) Additional information: the deskewed rectangle has a width $W = 21\ cm$ and a height $H = 29.7\ cm$. Also $A_1 A_2 A_4 A_3$ is a square of 1 cm x 1 cm, and $A_5 A_6 A_8 A_7$ too.

2) The problem will probably have multiple solutions, up to rotations of 90 or 180 degrees

3) It's probably possible with only 6 points (or even less?) but I'd like to make use of all the 8 points to improve quality (it's an overdetermined problem, so using least-squares could probably help, but I don't know how to do it in this context)

Answer 1

1) Only four points (like four corners of a sheet of paper) are enough to do the deskewing, the transform is known as a "homography" (as stated in another answer). See also this answer about these transformations. Here is how to do it on an image with only a few lines of Python + OpenCV:

2) When we want to use 8 points to do it, the problem is overdetermined, and we want to find the optimal solution by minimising an error (least-squares style). It's exactly what findHomography from the library OpenCV does (see remarks there about minimization of error). Here is how to do it with Python + OpenCV:

import cv2 
import numpy as np
import matplotlib.pyplot as plt

img = cv2.imread('IMG_20180522_211242_2.jpg')
pts1 = np.float32([[58,642],[147,627],[83,733],[168,716],[2320,2654],[2291,2566],[2238,2675],[2211,2588]])
pts2 = np.float32([[50,50],[150,50],[50,150],[150,150],[2100-50, 2970-50],[2100-50,2970-150],[2100-150,2970-50],[2100-150,2970-150]])
size = (2100,2970)

M, mask = cv2.findHomography(pts1,pts2)
dst = cv2.warpPerspective(img,M,size)

plt.subplot(121),plt.imshow(img),plt.title('Input')
plt.subplot(122),plt.imshow(dst),plt.title('Output')
plt.show()

Note: I now realize this is a programming+math question/answer (if so, feel free to migrate it to StackOverflow), or maybe half math/half programming, and then another more math-oriented-answer is welcome.

Answer 2

You need to compute a projective transformation (aka homography, collineation). Details at Finding the Transform matrix from 4 projected points (with Javascript)

If you want to pursue least squares methods you might find this Mathematica answer helpful: https://mathematica.stackexchange.com/questions/9244/solve-system-of-equations-related-to-perspective-projection

Answer 3

QuadrilateralGrid[q_List, nx_, ny_] := Module[{u, i, grf, graphics = {}}, For[i = 1, i <= nx, i++, u = i/nx; AppendTo[graphics, ParametricPlot[(q[[1]] u + q[[2]] (1 - u)) v + (q[[4]] u + q[[3]] (1 - u)) (1 - v), {v, 0, 1}]]]; For[i = 1, i <= ny, i++, u = i/ny; AppendTo[graphics, ParametricPlot[(q[[2]] u + q[[3]] (1 - u)) v + (q[[1]] u + q[[4]] (1 - u)) (1 - v), {v, 0, 1}]]]; AppendTo[graphics, ListLinePlot[{q[[1]], q[[2]], q[[3]], q[[4]], q[[1]]}, PlotStyle -> Red]]; Return[graphics] ]

QuadrilateralCoords[q_List, p_List] := Module[{coords = {}, i, x, y, sol, c0 = q[[3]], c1 = q[[4]] - q[[3]], c2 = q[[2]] - q[[3]], c3 = q[[1]] - q[[4]] - q[[2]] + q[[3]], s1, s2, s0}, For[i = 1, i <= Length[p], i++, sol = Quiet[ Solve[Thread[c0 + x c1 + y c2 + x y c3 == p[[i]]], {x, y}]]; If[Length[sol] < 2, s0 = ({x, y} /. sol)[[1]], s1 = {x, y} /. sol[[1]]; s2 = {x, y} /. sol[[2]]; s0 = s2; If[s1[[1]] > 0 && s1[[1]] < 1 && s1[[2]] > 0 && s1[[2]] < 1, s0 = s1]]; AppendTo[coords, s0] ]; Return[coords] ] QuadrilateralPlotPoints[q_List, coords_] := Module[{graphics = {}, nc = Length[coords], i, u, v, p1, p2, p0, , diag1 = Norm[q[[1]] - q[[3]]], diag2 = Norm[q[[2]] - q[[4]]], diag, rad}, diag = Max[diag1, diag2]; rad = 0.008*diag; For[i = 1, i <= nc, i++, u = coords[[i, 1]]; v = coords[[i, 2]]; p1 = q[[1]] u + q[[2]] (1 - u); p2 = q[[4]] u + q[[3]] (1 - u); p0 = p1 v + p2 (1 - v); AppendTo[graphics, Graphics[{Red, Disk[p0, rad]}]]]; AppendTo[graphics, ListLinePlot[{q[[1]], q[[2]], q[[3]], q[[4]], q[[1]]}, PlotStyle -> Red]]; Return[graphics] ] quadrilateral = {{0, 0}, {2, 2}, {1, 5}, {-3, 2}}; base = {{0, 0}, {2, 0}, {2, 3}, {0, 3}}; Show[QuadrilateralGrid[quadrilateral, 10, 10], PlotRange -> All ] Show[QuadrilateralGrid[base, 10, 10], PlotRange -> All ]

basepoints = Table[{0.5 Cos[t] + 1, 0.5 Sin[t] + 1}, {t, 0, 2 Pi, 0.3}]; coordsbase = QuadrilateralCoords[base, basepoints]; Show[QuadrilateralPlotPoints[quadrilateral, coordsbase]] Show[QuadrilateralPlotPoints[base, coordsbase]]