Eight-point algorithm

With

${\displaystyle \mathbf {y} ={\begin{pmatrix}y_{1}\\y_{2}\\1\end{pmatrix}}}$   and   ${\displaystyle \mathbf {y} '={\begin{pmatrix}y'_{1}\\y'_{2}\\1\end{pmatrix}}}$   and   ${\displaystyle \mathbf {E} ={\begin{pmatrix}e_{11}&e_{12}&e_{13}\\e_{21}&e_{22}&e_{23}\\e_{31}&e_{32}&e_{33}\end{pmatrix}}}$

the constraint can also be rewritten as

${\displaystyle y'_{1}y_{1}e_{11}+y'_{1}y_{2}e_{12}+y'_{1}e_{13}+y'_{2}y_{1}e_{21}+y'_{2}y_{2}e_{22}+y'_{2}e_{23}+y_{1}e_{31}+y_{2}e_{32}+e_{33}=0\,}$

or

${\displaystyle \mathbf {e} \cdot {\tilde {\mathbf {y} }}=0}$

where

${\displaystyle {\tilde {\mathbf {y} }}={\begin{pmatrix}y'_{1}y_{1}\\y'_{1}y_{2}\\y'_{1}\\y'_{2}y_{1}\\y'_{2}y_{2}\\y'_{2}\\y_{1}\\y_{2}\\1\end{pmatrix}}}$   and   ${\displaystyle \mathbf {e} ={\begin{pmatrix}e_{11}\\e_{12}\\e_{13}\\e_{21}\\e_{22}\\e_{23}\\e_{31}\\e_{32}\\e_{33}\end{pmatrix}}}$

that is, ${\displaystyle \mathbf {e} }$ represents the essential matrix in the form of a 9-dimensional vector and this vector must be orthogonal to the vector ${\displaystyle {\tilde {\mathbf {y} }}}$ which can be seen as a vector representation of the ${\displaystyle 3\times 3}$ matrix ${\displaystyle \mathbf {y} '\,\mathbf {y} ^{T}}$.

Each pair of corresponding image points produces a vector ${\displaystyle {\tilde {\mathbf {y} }}}$. Given a set of 3D points ${\displaystyle \mathbf {P} _{k}}$ this corresponds to a set of vectors ${\displaystyle {\tilde {\mathbf {y} }}_{k}}$ and all of them must satisfy

${\displaystyle \mathbf {e} \cdot {\tilde {\mathbf {y} }}_{k}=0}$

for the vector ${\displaystyle \mathbf {e} }$. Given sufficiently many (at least eight) linearly independent vectors ${\displaystyle {\tilde {\mathbf {y} }}_{k}}$ it is possible to determine ${\displaystyle \mathbf {e} }$ in a straightforward way. Collect all vectors ${\displaystyle {\tilde {\mathbf {y} }}_{k}}$ as the columns of a matrix ${\displaystyle \mathbf {Y} }$ and it must then be the case that

${\displaystyle \mathbf {e} ^{T}\,\mathbf {Y} =\mathbf {0} }$

This means that ${\displaystyle \mathbf {e} }$ is the solution to a homogeneous linear equation.

A standard approach to solving this equation implies that ${\displaystyle \mathbf {e} }$ is a right singular vector of ${\displaystyle \mathbf {Y} }$ corresponding to a singular value that equals zero. Provided that at least eight linearly independent vectors ${\displaystyle {\tilde {\mathbf {y} }}_{k}}$ are used to construct ${\displaystyle \mathbf {Y} }$ it follows that this singular vector is unique (disregarding scalar multiplication) and, consequently, ${\displaystyle \mathbf {e} }$ and then ${\displaystyle \mathbf {E} }$ can be determined.

In the case that more than eight corresponding points are used to construct ${\displaystyle \mathbf {Y} }$ it is possible that it does not have any singular value equal to zero. This case occurs in practice when the image coordinates are affected by various types of noise. A common approach to deal with this situation is to describe it as a total least squares problem; find ${\displaystyle \mathbf {e} }$ which minimizes

${\displaystyle \|\mathbf {e} ^{T}\,\mathbf {Y} \|}$

when ${\displaystyle \|\mathbf {e} \|=1}$. The solution is to choose ${\displaystyle \mathbf {e} }$ as the left singular vector corresponding to the smallest singular value of ${\displaystyle \mathbf {Y} }$. A reordering of this ${\displaystyle \mathbf {e} }$ back into a ${\displaystyle 3\times 3}$ matrix gives the result of this step, here referred to as ${\displaystyle \mathbf {E} _{\rm {est}}}$.

Another consequence of dealing with noisy image coordinates is that the resulting matrix may not satisfy the internal constraint of the essential matrix, that is, two of its singular values are equal and nonzero and the other is zero. Depending on the application, smaller or larger deviations from the internal constraint may or may not be a problem. If it is critical that the estimated matrix satisfies the internal constraints, this can be accomplished by finding the matrix ${\displaystyle \mathbf {E} '}$ of rank 2 which minimizes

${\displaystyle \|\mathbf {E} '-\mathbf {E} _{\rm {est}}\|}$

where ${\displaystyle \mathbf {E} _{\rm {est}}}$ is the resulting matrix from Step 2 and the Frobenius matrix norm is used. The solution to the problem is given by first computing a singular value decomposition of ${\displaystyle \mathbf {E} _{\rm {est}}}$:

${\displaystyle \mathbf {E} _{\rm {est}}=\mathbf {U} \,\mathbf {S} \,\mathbf {V} ^{T}}$

where ${\displaystyle \mathbf {U} ,\mathbf {V} }$ are orthogonal matrices and ${\displaystyle \mathbf {S} }$ is a diagonal matrix which contains the singular values of ${\displaystyle \mathbf {E} _{\rm {est}}}$. In the ideal case, one of the diagonal elements of ${\displaystyle \mathbf {S} }$ should be zero, or at least small compared to the other two which should be equal. In any case, set

${\displaystyle \mathbf {S} '={\begin{pmatrix}s_{1}&0&0\\0&s_{2}&0\\0&0&0\end{pmatrix}},}$

where ${\displaystyle s_{1},s_{2}}$ are the largest and second largest singular values in ${\displaystyle \mathbf {S} }$ respectively. Finally, ${\displaystyle \mathbf {E} '}$ is given by

${\displaystyle \mathbf {E} '=\mathbf {U} \,\mathbf {S} '\,\mathbf {V} ^{T}}$

The matrix ${\displaystyle \mathbf {E} '}$ is the resulting estimate of the essential matrix provided by the algorithm.

The problem is that the resulting ${\displaystyle \mathbf {Y} }$ often is ill-conditioned. In theory, ${\displaystyle \mathbf {Y} }$ should have one singular value equal to zero and the rest are non-zero. In practice, however, some of the non-zero singular values can become small relative to the larger ones. If more than eight corresponding points are used to construct ${\displaystyle \mathbf {Y} }$, where the coordinates are only approximately correct, there may not be a well-defined singular value which can be identified as approximately zero. Consequently, the solution of the homogeneous linear system of equations may not be sufficiently accurate to be useful.

Hartley addressed this estimation problem in his 1997 article. His analysis of the problem shows that the problem is caused by the poor distribution of the homogeneous image coordinates in their space, ${\displaystyle \mathbb {R} ^{3}}$. A typical homogeneous representation of the 2D image coordinate ${\displaystyle (y_{1},y_{2})\,}$ is

${\displaystyle \mathbf {y} ={\begin{pmatrix}y_{1}\\y_{2}\\1\end{pmatrix}}}$

where both ${\displaystyle y_{1},y_{2}\,}$ lie in the range 0 to 1000–2000 for a modern digital camera. This means that the first two coordinates in ${\displaystyle \mathbf {y} }$ vary over a much larger range than the third coordinate. Furthermore, if the image points which are used to construct ${\displaystyle \mathbf {Y} }$ lie in a relatively small region of the image, for example at ${\displaystyle (700,700)\pm (100,100)\,}$, again the vector ${\displaystyle \mathbf {y} }$ points in more or less the same direction for all points. As a consequence, ${\displaystyle \mathbf {Y} }$ will have one large singular value and the remaining are small.

As a solution to this problem, Hartley proposed that the coordinate system of each of the two images should be transformed, independently, into a new coordinate system according to the following principle.

• The origin of the new coordinate system should be centered (have its origin) at the centroid (center of gravity) of the image points. This is accomplished by a translation of the original origin to the new one.
• After the translation the coordinates are uniformly scaled so that the mean of distances from the origin to the points equals ${\displaystyle {\sqrt {2}}}$.

This principle results, normally, in a distinct coordinate transformation for each of the two images. As a result, new homogeneous image coordinates ${\displaystyle \mathbf {\bar {y}} ,\mathbf {\bar {y}} '}$ are given by

${\displaystyle \mathbf {\bar {y}} =\mathbf {T} \,\mathbf {y} }$

${\displaystyle \mathbf {\bar {y}} '=\mathbf {T} '\,\mathbf {y} '}$

where ${\displaystyle \mathbf {T} ,\mathbf {T} '}$ are the transformations (translation and scaling) from the old to the new normalized image coordinates. This normalization is only dependent on the image points which are used in a single image and is, in general, distinct from normalized image coordinates produced by a normalized camera.

The epipolar constraint based on the fundamental matrix can now be rewritten as

${\displaystyle 0=(\mathbf {\bar {y}} ')^{T}\,((\mathbf {T} ')^{T})^{-1}\,\mathbf {F} \,\mathbf {T} ^{-1}\,\mathbf {\bar {y}} =(\mathbf {\bar {y}} ')^{T}\,\mathbf {\bar {F}} \,\mathbf {\bar {y}} }$

where ${\displaystyle \mathbf {\bar {F}} =((\mathbf {T} ')^{T})^{-1}\,\mathbf {F} \,\mathbf {T} ^{-1}}$. This means that it is possible to use the normalized homogeneous image coordinates ${\displaystyle \mathbf {\bar {y}} ,\mathbf {\bar {y}} '}$ to estimate the transformed fundamental matrix ${\displaystyle \mathbf {\bar {F}} }$ using the basic eight-point algorithm described above.

The purpose of the normalization transformations is that the matrix ${\displaystyle \mathbf {\bar {Y}} }$, constructed from the normalized image coordinates, in general, has a better condition number than ${\displaystyle \mathbf {Y} }$ has. This means that the solution ${\displaystyle \mathbf {\bar {f}} }$ is more well-defined as a solution of the homogeneous equation ${\displaystyle \mathbf {\bar {Y}} \,\mathbf {\bar {f}} }$ than ${\displaystyle \mathbf {f} }$ is relative to ${\displaystyle \mathbf {Y} }$. Once ${\displaystyle \mathbf {\bar {f}} }$ has been determined and reshaped into ${\displaystyle \mathbf {\bar {F}} }$ the latter can be de-normalized to give ${\displaystyle \mathbf {F} }$ according to

${\displaystyle \mathbf {F} =(\mathbf {T} ')^{T}\,\mathbf {\bar {F}} \,\mathbf {T} }$

In general, this estimate of the fundamental matrix is a better one than would have been obtained by estimating from the un-normalized coordinates.