Find perspective projection matrix under constraints

Question

I have taken a photo (with a rectilinear lens) of a house's outside wall that sits on top a roof, from which I wish to take measurements from. However, to do so I would need to correct the perspective so the wall appears parallel to the image

To assist with this, I have some known properties of the wall, which I have labelled above:

• ∠172 = 40° (the roof is 40° from the vertical)
• ∠756 = 90°, ∠734 = 90° (the white "slats" are horizontal)
• Points 1,3,5,7 are colinear, and divide the segment into thirds (the slats have identical heights)
• Point 6 is the midpoint of points 4 and 7 (similar triangles)
• Points 2,4,6,7 are colinear

---

First of all, lets make some definitions:

$$ \begin{align} n &= \textrm{Point number as indicated on the diagram} \\ \textrm{image}[n] &= \mathbb{R}^2 \textrm{ coordinates of points in the original (left) image} \\ \textrm{error}[n] &= \mathbb{R}^2 \textrm{ measurement error of the above coordinates} \\ \textrm{base}[n] &= \mathbb{R}^2 \textrm{ coordinates of points in the projected (right) image} \\ \textrm{projection} &= \begin{pmatrix} p[1,1] & p[1,2] & p[1,3] \\ p[2,1] & p[2,2] & p[2,3] \\ p[3,1] & p[3,2] & 1 \end{pmatrix} \textrm{ arbitrary projection using }\mathbb{R}^3\textrm{ homogeneous coordinates} \\ \textrm{cosAngle}(\textrm{a, b, c}) &= \frac{(\textrm{a} - \textrm{b})\cdot(\textrm{c} - \textrm{b})}{|\textrm{a} - \textrm{b}| |\textrm{c} - \textrm{b}|} \textrm{ cosine of angle ∠abc} \end{align} $$

So this problem simply becomes an optimisation problem where we want to minimise the $\textrm{error}[n]$ terms, and then extract the $\textrm{projection}$ or $\textrm{base}$ terms, under the following constraints:

$$ \begin{align} & \textrm{(base is a perspective transform of image)} \\ \textrm{base}[n] &= \frac{\textrm{projection} \begin{pmatrix}\textrm{image}[n] + \textrm{error}[n] \\ 1\end{pmatrix}}{\begin{pmatrix}p[3,1] \\ p[3,2] \\ 1\end{pmatrix} \cdot \begin{pmatrix}\textrm{image}[n] + \textrm{error}[n] \\ 1\end{pmatrix}} \\ & \textrm{(arbitrary reference coordinates to absorb rotation, scale, and translation degrees of freedom)} \\ \textrm{base}[7] &= \begin{pmatrix}0 \\ 0\end{pmatrix} \\ \textrm{base}[1] &= \begin{pmatrix}0 \\ -1\end{pmatrix} \\ & \textrm{(colinearity and proportionality)} \\ \textrm{base}[5] &= \frac{1}{3} \textrm{base}[1] \\ \textrm{base}[3] &= \frac{2}{3} \textrm{base}[1] \\ \textrm{base}[6] &= \frac{1}{2} \textrm{base}[4] \\ \cos 0 &= \textrm{cosAngle}(\textrm{base}[4], \textrm{base}[7], \textrm{base}[2]) \\ & \textrm{(angles)} \\ \cos 40° &= \textrm{cosAngle}(\textrm{base}[1], \textrm{base}[7], \textrm{base}[2]) \\ \cos 90° &= \textrm{cosAngle}(\textrm{base}[7], \textrm{base}[5], \textrm{base}[6]) \\ \cos 90° &= \textrm{cosAngle}(\textrm{base}[7], \textrm{base}[3], \textrm{base}[4]) \end{align} $$

---

Here is my attempt to perform least squares minimisation with Mathematica 12.0

image = {
    {2871.87,1626.04},
    {3240.15,1722.65},
    {2863.93,1794.74},
    {3035.29,1953.99},
    {2855.83,1966.57},
    {2951.74,2048.33},
    {2846.37,2167.32}
};
error = Table[{ex[n], ey[n]}, {n, 1, 7}];
base = Table[{bx[n], by[n]}, {n, 1, 7}];
projection = Table[If[x == 3 && y == 3, 1, p[x,y]], {x, 1, 3}, {y, 1, 3}];
cosAngle[a_, b_, c_] := (a-b).(c-b) / Norm[a-b] / Norm[c-b];

constraints = {
    ((projection.Join[#, {1}])/(projection[[3]].Join[#, {1}]) & /@ (image + error))[[All, 1;;2]] == base,
    (*arbitrary reference coordinates*)
    base[[7]] == {0, 0},
    base[[1]] == {0, -1},
    (*colinearity and proportionality*)
    base[[5]] == base[[1]] / 3,
    base[[3]] == 2 base[[1]] / 3,
    base[[6]] == base[[4]] / 2,
    cosAngle[base[[4]],base[[7]],base[[2]]] == Cos[0],
    (*angles*)
    cosAngle[base[[1]],base[[7]],base[[2]]] == Cos[40 Degree],
    cosAngle[base[[7]],base[[5]],base[[6]]] == Cos[90 Degree],
    cosAngle[base[[7]],base[[3]],base[[4]]] == Cos[90 Degree]
};

NMinimize[{
    Total[#.# & /@ error],
    And@@constraints
}, Select[Flatten@Join[projection, error, base], !SameQ[#, 1] &]]

But Mathematica is unable to find a solution, erroring with NMinimize::nosat: Obtained solution does not satisfy the following constraints... with a list of 13 constraints afterwards

Am I approaching this problem correctly? Have I made some incorrect assumptions?
Why is Mathematica unable to find a solution?