why perspectivity can be represented by a projective transformation with 6 degrees of freedom?

Question

Original problem comes from R.Hartley & A.Zisserman Multiple View Geometry in Computer Vision at page 34:

So why dof(degree of freedom) of a projective transformation: $$x'=\begin{bmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & h_{33} \end{bmatrix}x=Hx$$ decrease from 8 to 6 in a perspectivity case? How to embody so-called restrictions here?

Answer 1

Summary

Based on strong experimental evidence I would say the statement is incorrect. Essentially what I did was create a bunch of random $3\times3$ matrices, which presumably exercise all the degrees of freedom of a projective transformation, and then try to obtain a 3d embedding for this. This succeeded more often than could plausibly be explained by them accidentially satisfying two more constraints.

Setup

Here is my setup. Without loss of generality, use the coordinate system of the target plane as the frame of reference. So the target point will have coordinates

$$p'=\begin{pmatrix}x'\\y'\\0\end{pmatrix}\;.$$

In relation to this the source plane is rotated and translated. According to the problem statement it might also be scaled, but for every setup with a scaled version there is one with a non-scaled plane at a different position, so only considering rotation and translation is enough. I describe the rotation using a quaternion $q_r+q_ii+q_jj+q_kk$ with $q_r^2+q_i^2+q_j^2+q_k^2=1$ corresponding to a rotation matrix. Using two vectors of that matrix as unit vectors and a point $o$ as the origin of the source plane, we get the source point as

$$p=\begin{pmatrix}o_x\\o_y\\o_z\end{pmatrix}+ x\begin{pmatrix} q_r^2 + q_i^2 - q_j^2 - q_k^2 \\ 2(q_i q_j + q_r q_k) \\ 2(q_i q_k - q_r q_j) \end{pmatrix}+ y\begin{pmatrix} 2(q_i q_j - q_r q_k) \\ q_r^2 - q_i^2 + q_j^2 - q_k^2 \\ 2(q_r q_i + q_j q_k) \end{pmatrix}\;.$$

Also assume the center of the perspectivity to be some point

$$c=\begin{pmatrix}c_x\\c_y\\c_z\end{pmatrix}\;.$$

Then the line $pc$ can be intersected with the plane $z=0$ to obtain $p'$. Formulating that transformation in homogeneous coordinates the transformation matrix describing this operation is

$$B=\scriptstyle\begin{pmatrix} 2(q_iq_k-q_rq_j)c_x - (qr^2+qi^2-qj^2-qk^2)cz & 2(q_rq_i+q_jq_k)c_x + 2(q_rq_k-q_iq_j)c_z & o_zc_x - o_xc_z \\ 2(q_iq_k-q_rq_j)c_y - 2(q_iq_j+q_rq_k)c_z & 2(q_rq_i+q_jq_k)c_y - (q_r^2-q_i^2+q_j^2-q_k^2)c_z & o_zc_y - o_yc_z \\ 2q_iq_k - 2q_rq_j & 2q_jq_k + 2q_rq_i & o_z - c_z \end{pmatrix}$$

unless I made a mistake copying this from the Sage output quoted below. Now pick a matrix $M$ with $3\times3$ random entries. I used integers with a magnitude of at most $1000$, equidistributed. This should arguably exercise all the degrees of freedom a projective transformation matrix has. We have a 3d setup corresponding such a matrix $M$ if the resulting matrices are multiples of one another, i.e. $sM=B$ for some $s\in\mathbb\setminus\{0\}$.

At this point we have $11$ variables: $q_r,q_i,q_j,q_k,o_x,o_y,o_z,c_x,c_y,c_z,s$. We have $10$ non-linear constraints on these: $9$ from the above matrix equation and one from the condition for a unit quaternion. In general this will lead to an infinite number of solutions. So one can pick one additional constraint and then attempt to find a finite number of solutions. Picking $c_z=-1$ I managed to find real solutions for bery roughly half the randomized matrices I tried. For the others some parameters would have been complex-valued. But by choosing any of

$$c_z=-1\quad c_z=-1000\quad c_z=-1/1000\quad o_z=1\quad o_z=1000\quad o_z=1/1000$$

I was able to find a 3d embedding for all $1000$ random matrices that I tried. If there really were only 6 real degrees of freedom, then there would be two more constraints, and me hitting matrices satisfying these constraints a thounsand times in a row is too implausible to believe, I think. Obviously this is not a proof, though.

Code

Here is my Sage code:

PR1.<qr,qi,qj,qk,ox,oy,oz,cx,cy,cz,s> = QQ[]
PR2.<x,y> = PR1[]
qv = vector([qi, qj, qk])
qv.cross_product_matrix()
A = ((qr^2-qv*qv)*identity_matrix(3)
   + 2*qv.column()*qv.row()
   + 2*qr*qv.cross_product_matrix())  # rotation matrix
C = vector([cx,cy,cz])  # center of perspectivity
O = vector([ox,oy,oz])  # origin of source plane
P1 = O + x*A.column(0) + y*A.column(1)
C1 = vector(C.list() + [1])  # homogeneous form
P11 = vector(P1.list() + [1])  # likewise
P2 = P11[2]*C1 - C1[2]*P11  # intersection with z=0
B1 = matrix(PR1, [[
  t.coefficient(x), t.coefficient(y), t.constant_coefficient()] for t in P2])
assert B1 * vector([x,y,1]) == P2  # sanity check
B2 = B1.matrix_from_rows([0,1,3])  # omit zero row for z coordinate
B2.list()
#  [-2*qr*qj*cx + 2*qi*qk*cx - qr^2*cz - qi^2*cz + qj^2*cz + qk^2*cz,
#   2*qr*qi*cx + 2*qj*qk*cx - 2*qi*qj*cz + 2*qr*qk*cz,
#   oz*cx - ox*cz,
#   -2*qr*qj*cy + 2*qi*qk*cy - 2*qi*qj*cz - 2*qr*qk*cz,
#   2*qr*qi*cy + 2*qj*qk*cy - qr^2*cz + qi^2*cz - qj^2*cz + qk^2*cz,
#   oz*cy - oy*cz,
#   -2*qr*qj + 2*qi*qk,
#   2*qr*qi + 2*qj*qk,
#   oz - cz]

def random_experiments(n, max_magnitude=1000):
    solved = []
    not_solved = []
    wrong_dimension = []
    for i in xrange(n):
        M = random_matrix(ZZ, 3, 3, x=-max_magnitude, y=max_magnitude+1)
        dim0 = False
        for eq in [cz+1, 1000*cz+1, cz+1000, oz-1, 1000*oz-1, oz-1000]:
            id = ideal((s*M - B2).list() + [qr^2+qi^2+qj^2+qk^2-1, eq])
            if not id.dimension():
                dim0 = True
                v = id.variety(AA)
                if v:
                    solved.append((M, v))
                    break
        else:
            if dim0:
                not_solved.append(M)
            else:
                print("dimension {}".format(id.dimension()))
                wrong_dimension.append(M)
    res = [solved, not_solved, wrong_dimension]  # Lists of matrices
    res.append(list(map(len, res)))  # result[3] is lengths of lists
    res.append(float(len(solved)/n))  # result[4] is fraction solved
    return tuple(res)

res = random_experiments(1000)
print res1[3:]
# ([1000, 0, 0], 1.0)

Numeric example

To give you an idea of what the randomized examples and their 3d embedding look like, here is one. I've reduced the maximal magnitude of the random matrix elements to $9$ in order to keep expressions a bit simpler. The random matrix for this example is

$$M=\begin{pmatrix} -3 & 7 & 0 \\ 2 & 7 & 8 \\ 9 & 2 & 7 \end{pmatrix}\;.$$

It allowed for $8$ choices of parameters with $c_z=-1$. Since negating the quaternion represents the same rotation, this only represents $4$ actual embeddings in 3d. The first of them was

\begin{align*} -0.29331 \approx q_r &= {\scriptstyle -\frac{1}{2} \, \sqrt{\frac{1179}{700570} \, \sqrt{113} \sqrt{85} - 2 \, \sqrt{\frac{1179}{1401140} \, \sqrt{113} \sqrt{85} + \frac{120447}{1401140}} + 1}} \\ 0.77963 \approx q_i &= {\scriptstyle \frac{1}{2} \, \sqrt{\frac{2007}{700570} \, \sqrt{113} \sqrt{85} + 2 \, \sqrt{\frac{2007}{1401140} \, \sqrt{113} \sqrt{85} + \frac{267003}{1401140}} + 1}} \\ -0.26040 \approx q_j &= {\scriptstyle -\frac{1}{2} \, \sqrt{-\frac{2007}{700570} \, \sqrt{113} \sqrt{85} - 2 \, \sqrt{-\frac{2007}{1401140} \, \sqrt{113} \sqrt{85} + \frac{267003}{1401140}} + 1}} \\ -0.48819 \approx q_k &= {\scriptstyle -\frac{1}{2} \, \sqrt{-\frac{1179}{700570} \, \sqrt{113} \sqrt{85} + 2 \, \sqrt{-\frac{1179}{1401140} \, \sqrt{113} \sqrt{85} + \frac{120447}{1401140}} + 1}} \\ 0.15549 \approx o_x &= {\scriptstyle -\frac{7}{53890} \, \sqrt{26945} \sqrt{26} + \frac{1}{2} \, \sqrt{\frac{49562478}{122699581225} \, \sqrt{26945} \sqrt{26} + \frac{43922976039}{122699581225}} - \frac{13}{85}} \\ -0.65619 \approx o_y &= {\scriptstyle -\frac{207}{700570} \, \sqrt{4121} \sqrt{226} - \frac{1}{2} \, \sqrt{\frac{1321488}{1443524485} \, \sqrt{4121} \sqrt{226} + \frac{390212685}{288704897}} + \frac{32}{85}} \\ -1.7109 \approx o_z &= {\scriptstyle -\frac{7}{8242} \, \sqrt{700570} - 1} \\ 0.090883 \approx c_x &= {\scriptstyle \frac{177}{700570} \, \sqrt{931346} - \frac{13}{85}} \\ 0.091320 \approx c_y &= {\scriptstyle -\frac{207}{700570} \, \sqrt{931346} + \frac{32}{85}} \\ -1.0000 = c_z &= {\scriptstyle -1} \\ -0.10155 \approx s &= {\scriptstyle -\frac{1}{8242} \, \sqrt{8242} \sqrt{85}}\end{align*}

So for example for $x=2, y=5$ you get the preimage

$$p\approx\begin{pmatrix}0.15549\\-0.65619\\-1.7109\end{pmatrix} +2\begin{pmatrix}0.38772\\-0.11964\\-0.91398\end{pmatrix} +5\begin{pmatrix}-0.69241\\-0.69232\\-0.20311\end{pmatrix}\approx \begin{pmatrix}-2.5311\\-4.3571\\-4.5544\end{pmatrix}$$

and its image

\begin{align*} p'&=\left(1+\frac{\sqrt{700570}}{2975}\right)c-\frac{\sqrt{700570}}{2975}p\\& \approx1.2813\begin{pmatrix}0.090883\\0.091320\\-1.0000\end{pmatrix} -0.2813\begin{pmatrix}-2.5311\\-4.3571\\-4.5544\end{pmatrix} \approx\begin{pmatrix}0.82857\\1.3429\\0\end{pmatrix} \end{align*}

which maches the projective transformation

$$\begin{pmatrix} -3 & 7 & 0 \\ 2 & 7 & 8 \\ 9 & 2 & 7 \end{pmatrix}\cdot \begin{pmatrix}2\\5\\1\end{pmatrix} = \begin{pmatrix}29\\47\\35\end{pmatrix} \approx 35\begin{pmatrix}0.82857\\1.3429\\1\end{pmatrix} \;.$$

Answer 2

Assume that the projection plane is fixed wrt the reference frame. Then the projected plane has the six degrees of freedom of any solid (three translational and three rotational), each influencing the projection independently of the others.

Even though the matrix has nine coefficients, they are no independent. Some compatibility relations are required to express a rigid motion rather than an affinity. Should you allow arbitrary coefficients, the projected plane will appear deformed. (For instance you could see a rectangle parallel to the image plane projected like a parallelogram.)