2
$\begingroup$

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?

enter image description here

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)

$\endgroup$
  • $\begingroup$ Even if you fix the "target" rectangle of assigned height and width, the mapping is not unique because of the rotational symmetry of the rectangle, e.g. it could be rotated through an angle $\pi$ without change. $\endgroup$ – hardmath May 20 '18 at 18:29
  • $\begingroup$ Already answered at math.stackexchange.com/q/296794/1257 $\endgroup$ – brainjam May 20 '18 at 18:37
  • $\begingroup$ @brainjam it's related but not exactly the same (see my last remark in question) $\endgroup$ – Basj May 20 '18 at 18:41
  • $\begingroup$ Sorry, I didn't notice that you had 8 source points. Generally 4 is sufficient, along with 4 target points. Since you are posing an overdetermined problem you probably want something like a least squares solution, like pdfs.semanticscholar.org/ea51/… The problem is that you don't know what the target points would be that correspond to any of the source points other than $(x_1,y_1)$ and $(x_8,y_8)$. $\endgroup$ – brainjam May 21 '18 at 15:37
  • $\begingroup$ @brainjam Yes, it's overdetermined and I was also thinking about least squares, but not sure how to do it in this context. See my updated Note 1 in the question for more informations: I do know the target points for all points $A_i$, for $i \in {1, 2, ..., 8}$. $\endgroup$ – Basj May 21 '18 at 15:57
2
$\begingroup$

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 ] enter image description here

enter image description here

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

enter image description here

enter image description here

$\endgroup$
  • $\begingroup$ Thanks! Does this assume it's a parallelogram? In real situation it's not, the opposite sides are not necessary parallel. $\endgroup$ – Basj May 20 '18 at 20:16
  • $\begingroup$ Not at all. The skewed figure can be distorted as a general quadrilateral. $\endgroup$ – Cesareo May 20 '18 at 20:18
  • $\begingroup$ $M$ copes for the scales at any axis and $R$ copes with the rotation and $p_0$ with an eventual origin translation. $\endgroup$ – Cesareo May 20 '18 at 20:19
  • $\begingroup$ Doesn't really address how to do the minimization or limit the space of M $\endgroup$ – mathreadler May 20 '18 at 20:38
  • $\begingroup$ Is a standard minimization procedure that can be handled easily with NMinimize in MATHEMATICA for instance. $\endgroup$ – Cesareo May 20 '18 at 21:51
1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ it's related but not exactly the same (see my last remark in question) $\endgroup$ – Basj May 20 '18 at 18:41
1
$\begingroup$

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:

enter image description here

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()

enter image description here

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.