A Computer Algorithm for Reconstructing a Scene from Two Projections

“A Computer Algorithm for Reconstructing a Scene from Two Projections” is a three-page letter by H. Christopher Longuet-Higgins of the Laboratory of Experimental Psychology at the University of Sussex, published in Nature on 10 September 1981. It showed that if eight points can be matched between two perspective views, the relative orientation of the two viewpoints and the three-dimensional positions of the points can be computed by solving linear equations. The $3 \times 3$ matrix at the center of the method, which the paper calls $Q$, is the matrix now called the essential matrix, and the method is the origin of the eight-point algorithm. The letter was reprinted in Fischler and Firschein’s 1987 collection Readings in Computer Vision.

Longuet-Higgins was a theoretical chemist by training: he held the chair of theoretical chemistry at Cambridge from 1954 and left it in 1967 for Edinburgh to work on artificial intelligence, moving to Sussex in 1974 (University of Edinburgh notice). The IEEE PAMI Technical Committee’s Longuet-Higgins Prize, presented each year at CVPR to CVPR papers from ten years earlier that have had significant impact, is named after him.

Problem

The paper presents the problem as relevant to photographic surveying, to binocular vision, where the observer’s non-visual information about each eye’s orientation and focal length is far less accurate than the retinal images themselves, and to perceiving motion with one eye. Following Marr and Poggio, it separates two processes, finding corresponding points and using them to recover depth, and assumes the first has been solved. What remains is to find the relative orientation of the two viewpoints.

Photogrammetrists already knew that five matched points determine the relative orientation in general. Computing it, however, meant solving five simultaneous third-order equations iteratively; for this the paper cites E. H. Thompson’s 1959 formulation of relative orientation. Longuet-Higgins’s point was that with eight points the problem needs no iteration.

Contribution

• A bilinear constraint. Each match gives one equation that is linear in the nine entries of a single matrix built from the rotation and translation, so eight matches determine the matrix up to scale.
• Closed-form recovery. The translation direction, the rotation, and every point’s depth follow from the matrix by elementary algebra, with the sign ambiguity resolved by requiring points to lie in front of both viewpoints.
• Degenerate configurations, listed without proof, in which the eight equations are not independent.

Method

The paper assumes calibrated views: each point has image coordinates $x_\mu = X_\mu / X_3$, where $(X_1, X_2, X_3)$ are its coordinates relative to the first viewpoint and $X_3$ is the forward coordinate, and a dummy $x_3 = 1$ is appended. The second view is defined likewise, with primed coordinates related by

$$X'_\mu = R_{\mu\nu}(X_\nu - T_\nu),$$

with summation over repeated indices, $R$ a rotation, and the length of $T$ taken as the unit of distance. The paper defines

$$Q = RS, \qquad S = \begin{bmatrix} 0 & T_3 & -T_2 \\ -T_3 & 0 & T_1 \\ T_2 & -T_1 & 0 \end{bmatrix},$$

shows that $X'_\mu Q_{\mu\nu} X_\nu = 0$ because $S$ is antisymmetric, and divides by $X_3 X'_3$ to obtain the constraint between image coordinates, $x'_\mu Q_{\mu\nu} x_\nu = 0$.

In the site’s notation. The paper’s $x$ and $x'$ are the normalized coordinates $\hat{\mathbf{x}} = K^{-1}\mathbf{x}$ and $\hat{\mathbf{x}}' = K'^{-1}\mathbf{x}'$ of the epipolar geometry articles; the cameras are $P = K[\,I \mid \mathbf{0}\,]$ and $P' = K'[\,R \mid \mathbf{t}\,]$ with $\mathbf{t} = -RT$, so $T$ is the second camera’s center in the first camera’s frame. Since $S = -[T]_\times$ and $[RT]_\times = R[T]_\times R^\top$,

$$Q = -R[T]_\times = [\mathbf{t}]_\times R = E,$$

and the paper’s constraint is $\hat{\mathbf{x}}'^\top E\,\hat{\mathbf{x}} = 0$. With pixel coordinates, the same matrix appears as $F = K'^{-\top} E K^{-1}$, the fundamental matrix.

The algorithm, as the paper lists it, has six steps:

1. Write the equation $x'_\mu Q_{\mu\nu} x_\nu = 0$ for each of eight points and solve for the ratios of the $Q_{\mu\nu}$.
2. Scale $Q$ so that $\operatorname{tr}(Q^\top Q) = 2$, which holds for $|T| = 1$.
3. Read $T$ off $Q^\top Q = I - TT^\top$: the diagonal gives the magnitudes of its components, the off-diagonal entries their relative signs.
4. Recover the rotation row by row. Each row of $Q$ is $Q_\alpha = T \times R_\alpha$; with $W_\alpha = Q_\alpha \times T$, the rows are $R_\alpha = W_\alpha + W_\beta \times W_\gamma$ for $(\alpha, \beta, \gamma)$ a cyclic permutation of $(1, 2, 3)$.
5. Compute each point’s depth $X_3$ from its first image coordinate in the second view, $x'_1$, then $X_1 = x_1 X_3$, $X_2 = x_2 X_3$, and the primed coordinates from $X' = R(X - T)$.
6. Check that $X_3$ and $X'_3$ are positive. If both are negative, reverse $T$ and return to step 5; if they have opposite signs, reverse $Q$ and return to step 4.

The paper notes that the signs of $T$ and $Q$ give four distinct solutions, and that the positive-depth condition holds only when both signs are chosen correctly. This is the four-fold ambiguity, resolved by cheirality, of the essential matrix article.

The code below follows the six steps in the paper’s notation, using NumPy, on eight noise-free matches. It checks the recovered motion and the identity $Q = [\mathbf{t}]_\times R$, then computes the rank of the $8 \times 9$ equation matrix for two configurations the paper calls degenerate.

import numpy as np
from scipy.spatial.transform import Rotation

np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(0)

def skew(v):  # skew(v) @ u == np.cross(v, u)
    return np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])

# Longuet-Higgins's setup: X' = R (X - T), with |T| = 1 (eqs. 4-6)
R_true = Rotation.from_rotvec([0.05, -0.30, 0.10]).as_matrix()
T_true = np.array([1.0, 0.1, 0.2]) / np.linalg.norm([1.0, 0.1, 0.2])
X = rng.uniform([-2, -2, 4], [2, 2, 8], size=(8, 3))       # eight scene points
Xp = (X - T_true) @ R_true.T
x, xp = X / X[:, 2:], Xp / Xp[:, 2:]                      # image coordinates, x3 = 1

def eight_point(x, xp):
    """The 1981 procedure, steps 1-6, in the paper's notation."""
    A = np.einsum("ni,nj->nij", xp, x).reshape(-1, 9)     # x'_mu x_nu Q_mu,nu = 0 (eq. 13)
    Q = np.linalg.svd(A)[2][-1].reshape(3, 3)
    Q /= np.sqrt(0.5 * np.trace(Q.T @ Q))                 # step 2
    TT = np.eye(3) - Q.T @ Q                              # step 3: Q~Q = I - T T^T (eq. 17)
    k = np.argmax(np.diag(TT))
    T = TT[:, k] / np.sqrt(TT[k, k])
    def depths(R, T):                                     # step 5 (eqs. 31, 32, 4)
        r = R[0] - xp[:, :1] * R[2]
        X3 = (r @ T) / np.sum(r * x, axis=1)
        return X3, ((x * X3[:, None] - T) @ R.T)[:, 2]
    for _ in range(2):
        W = np.cross(Q, T)                                # step 4: W_a = Q_a x T (eq. 20)
        R = W + np.cross(W[[1, 2, 0]], W[[2, 0, 1]])      # R_a = W_a + W_b x W_g (eq. 27)
        X3, X3p = depths(R, T)
        if np.all(X3 < 0) and np.all(X3p < 0):            # step 6: both negative,
            T, Q = -T, -Q                                 # reverse T (Q = RS keeps its form)
            X3, X3p = depths(R, T)                        # and redo step 5
        if np.all(X3 > 0) and np.all(X3p > 0):
            return Q, R, T
        Q = -Q                                            # opposite signs: reverse Q, step 4

Q, R, T = eight_point(x, xp)
ang = np.degrees(np.linalg.norm(Rotation.from_matrix(R @ R_true.T).as_rotvec()))
print("rotation error (deg):", f"{ang:.1e}", "  |T - T_true|:", f"{np.linalg.norm(T - T_true):.1e}")

# The paper's Q is the essential matrix E = [t]x R of P' = K'[R | t], with t = -R T
t = -R_true @ T_true
print("max |Q - [t]x R|:", f"{np.abs(Q - skew(t) @ R_true).max():.1e}")

# Degenerate configurations named in the paper: rank of the 8 x 9 system
def rank(Xs):
    Xps = (Xs - T_true) @ R_true.T
    a, b = Xs / Xs[:, 2:], Xps / Xps[:, 2:]
    return np.linalg.matrix_rank(np.einsum("ni,nj->nij", b, a).reshape(-1, 9), tol=1e-9)

plane = np.c_[rng.uniform(-2, 2, (7, 2)), np.full(7, 6.0)]           # seven in a plane
cube = np.array([[i, j, k] for i in (-1, 1) for j in (-1, 1) for k in (-1, 1)]) + [0, 0, 6]
print("rank, general position:", rank(X))
print("rank, seven coplanar:  ", rank(np.vstack([plane, [0.5, -0.3, 7.5]])))
print("rank, cube vertices:   ", rank(cube))

Output:

rotation error (deg): 2.4e-13   |T - T_true|: 1.6e-15
max |Q - [t]x R|: 3.8e-15
rank, general position: 8
rank, seven coplanar:   7
rank, cube vertices:    7

With exact data the published steps recover the motion to rounding error, and the paper’s $Q$ coincides with $[\mathbf{t}]_\times R$. When seven of the points are coplanar, or the eight lie at the corners of a cube, the eight equations have rank seven and no longer fix $Q$. Step 1 here takes the null vector from a singular value decomposition, which with exactly eight points and exact data gives the same answer as solving the eight equations directly.

Results

The letter contains no experiments on images. Its results are the derivation and the algorithm, together with an accuracy estimate: if image directions are accurate to a few seconds of arc, the paper states, forward coordinates can be estimated with great accuracy out to about ten times the distance between the two centers of projection, and to about a hundred times if the points are well spread in depth.

On degeneracy, the paper states that the method fails when the eight equations are not independent, which happens if four of the points lie on a line or seven lie in a plane, and, unexpectedly, if the configuration includes six points at the vertices of a regular hexagon or consists of the eight vertices of a cube. It says the explanation is too long to give, that these cases are unrelated to any genuine ambiguity in interpreting the two projections, and that moving one offending point slightly removes the degeneracy.

Impact

The paper calls its matrix $Q$ and gives it no name. Huang and Faugeras (1989) wrote the matrix as $E$, described it as central to the eight-point linear algorithm, and proved that a $3 \times 3$ matrix factors as a skew-symmetric matrix times a rotation if and only if one singular value is zero and the other two are equal. By 1990 Horn referred to the so-called essential matrix, and Hartley (1997) credits Tsai and Huang with the step, now standard, of enforcing rank 2 by zeroing the smallest singular value.

The same linear method applies unchanged to pixel coordinates from uncalibrated cameras. Hartley (1997) notes that the matrix computed this way came to be called the fundamental matrix, and that it determines the reconstruction only up to a projective transformation.

The algorithm acquired a reputation for sensitivity to noise; the prevailing view, Hartley (1997) wrote, was that it was of little practical use. Horn (1990) had pointed out that a linear estimate does not satisfy the constraints of an essential matrix unless the data are exact, so it cannot be decomposed consistently. Hartley (1997) traced much of the poor performance to conditioning rather than to the method: translating each image’s points to their centroid and scaling them to an average distance of $\sqrt{2}$ from the origin before solving made the normalized eight-point algorithm, in his experiments on five image pairs, nearly as accurate as the best iterative methods while running about twenty times faster.

Limitations

The paper treats exact data. It solves eight equations exactly, gives no least-squares formulation for more points, and does not address noise in the matches or outliers. Its rotation formula assumes that $Q$ has exactly the structure of $RS$; with noisy data the estimated $Q$ does not, which is the inconsistency Horn (1990) discusses. It needs calibrated views, eight matches where five suffice in principle, and a baseline that is not small compared with the scene depth.

What Came After

The calibrated case moved toward minimal solvers. Nistér (2004) reduced the five-point problem to finding the roots of a tenth-degree polynomial and used the solver inside a hypothesize-and-test loop, the RANSAC setting in which relative pose is now usually estimated. The decomposition of $E$ into rotation and translation is now usually done with a singular value decomposition, and the choice among the four candidates by triangulating points and checking their depths, which is the paper’s step 6 in modern form.