In Defense of the Eight-Point Algorithm

“In Defense of the Eight-Point Algorithm” is a paper by Richard I. Hartley of GE Corporate Research and Development in Schenectady, New York, published in IEEE Transactions on Pattern Analysis and Machine Intelligence in June 1997. It grew out of “In Defence of the 8-Point Algorithm”, presented at the 1995 International Conference on Computer Vision (ICCV) in Cambridge, Massachusetts, which the journal paper cites for part of its comparison; this article treats the two as one publication. The paper argued that the linear eight-point algorithm for the fundamental matrix, then widely regarded as too sensitive to noise to be useful, failed because of poor numerical conditioning rather than any flaw in the method, and that translating and scaling the image coordinates before solving fixes it. The fix, often called Hartley normalization, became standard practice in linear estimation for multiple-view geometry, and the 1995 conference paper received the Helmholtz Prize in 2013. The fundamental matrix article covers the method as used today; this article covers the paper.

Problem

The eight-point algorithm goes back to Longuet-Higgins (1981), for the essential matrix of calibrated cameras; applied to pixel coordinates it computes the fundamental matrix. Each match gives one equation $\mathbf{x}'^\top F \mathbf{x} = 0$, linear in the entries of $F$; stacking them gives $A\mathbf{f} = \mathbf{0}$, and with noise $\mathbf{f}$ is the unit vector minimizing $\|A\mathbf{f}\|$, the least eigenvector of $A^\top A$.

The algorithm had a poor reputation. The abstract describes the prevailing view that it was extremely susceptible to noise and so “virtually useless for most purposes”. Many alternatives had been proposed, among them iterative methods minimizing distances in the image; Hartley cites INRIA comparisons by Luong, Deriche, Faugeras, Papadopoulo, and Zhang, and notes that all were considerably more complicated.

Contribution

• A diagnosis: the poor results came from implementations that ignored the conditioning of the linear system.
• A fix: in each image, translate the points so their centroid is at the origin and scale them so their average distance from it is $\sqrt{2}$, solve for $F$, and transform the result back.
• A justification from the condition number, covering both the linear solution and the rank-2 correction.
• Experiments on five real image pairs against the unnormalized algorithm and several iterative methods, mostly using implementations supplied by other researchers.

Method

The algorithm is not invariant to coordinates. If points are transformed as $\tilde{\mathbf{x}} = T\mathbf{x}$ and $\tilde{\mathbf{x}}' = T'\mathbf{x}'$, the matrix $T'^{-\top} F T^{-1}$ fits the transformed points with exactly the same algebraic residuals. But the constraint $\|\mathbf{f}\| = 1$ is not preserved by the change of variables, so the least-squares solution in the new coordinates, transformed back, is a different $F$. The result depends on the coordinate frame, even on whether the origin is at the corner or the center of the image, a problem Hartley compares to one Bookstein observed in conic fitting; the points should therefore be expressed in a canonical frame. (The calibrated coordinates $K^{-1}\mathbf{x}$ of the essential matrix are one such transformation; Hartley’s $T$ is computed from the points instead.)

Conditioning of the linear system. The relevant quantity is $\kappa = d_1/d_8$, the ratio of the largest to the second-smallest eigenvalue of $A^\top A$; the smallest, $d_9$, would be zero for exact data. Hartley shows that a perturbation of the size of $d_8$ can mix the last two eigenvectors and so change the solution drastically, while being a tiny relative change to $A^\top A$ when $\kappa$ is large. For a $200 \times 200$ image with points around $(100, 100, 1)^\top$, each row of $A$ has entries from about $10^4$ down to 1, so the diagonal of $A^\top A$ ranges from about $10^8$ down to 1. Using the interlacing property of eigenvalues of symmetric matrices, he bounds $\kappa \geq 10^8 / (10^4 + 1)$ and notes that in practice it is far larger. Scaling so that the homogeneous coordinates are on average about one reduces this imbalance. Translation helps because, with a large common offset, the informative part of each coordinate sits in its third or fourth significant figure.

Two scalings. Isotropic scaling applies one scale factor to both axes after centering, chosen so that the average distance from the origin is $\sqrt{2}$, which makes the average point $(1, 1, 1)^\top$. Non-isotropic scaling instead maps the centered points to a cloud whose two principal moments both equal one, using a Cholesky factorization of $\sum_i \mathbf{x}_i \mathbf{x}_i^\top$.

The rank-2 correction. The estimated $F$ is replaced by the closest singular matrix in Frobenius norm, by zeroing its smallest singular value, a step Hartley credits to Tsai and Huang (1984). Normalization matters here too. In pixel coordinates the entries of $F$ differ by orders of magnitude, about $10^{-4}$ in the top-left $2 \times 2$ block for points of size $(100, 100, 1)^\top$, and the Frobenius correction perturbs all entries by similar absolute amounts. The smallest entries, which change the most relative to their size, are the ones multiplied by the largest coordinates when computing epipolar lines. After normalization all entries have similar magnitude and are treated alike.

Results

The experiments used five image pairs, of houses, a statue, a museum in Grenoble, a calibration jig, and a corridor, with epipoles ranging from far from the image center to inside the image. Other researchers supplied the matches, found automatically and cleaned of outliers. For each number of points $N$, from eight up to three-quarters of the matches, each algorithm was run 100 times on random subsets. The error was the average distance, in both images, from each matched point to the epipolar line of its partner, over all matches rather than only those used in the fit. The reported results:

• Conditioning. On the house images, normalization improved the condition number by about $10^8$.
• Stages. Normalization helped the linear solution most, but helped the rank-2 correction significantly as well; the best results normalized both. For $N = 8$ it has no effect on the linear solution, since eight equations are solved exactly.
• Normalized against unnormalized. Better on all five pairs: slightly on the statue and corridor, where all methods performed similarly, and dramatically where the matches were less accurate.
• Against the optimal algorithm. Almost as good on all pairs as Hartley’s own “optimal” method, which minimizes the displacement of the matched points needed to satisfy a rank-2 $F$ exactly, the maximum-likelihood estimate under Gaussian noise.
• Against iterative methods. Zhengyou Zhang’s implementations of an iterative epipolar-distance method and a gradient-based method performed almost identically. Against the gradient-based method, initialized with the unnormalized algorithm, the normalized algorithm was better on the houses, worse on the museum, slightly worse on the calibration jig, and indistinguishable on the other two pairs. Initialized with the normalized algorithm, the iterative methods were not significantly different from the optimal one, which Hartley read as evidence that initialization matters more than the choice of cost function.
• Isotropic against non-isotropic. Almost indistinguishable, so the paper recommends the simpler isotropic scaling.
• Reconstruction. On the calibration jig, projective reconstructions from the normalized algorithm were almost as accurate as the optimal algorithm’s and far better than the unnormalized one’s.

The paper concludes that the normalized algorithm runs about 20 times faster than the best iterative methods and is far easier to code, while without normalization errors can reach 10 pixels. It names other linear problems where the same normalization applies directly: computing projective transformations between point sets, which includes the homography, estimating the trifocal tensor, and computing a projective camera matrix by the direct linear transform.

The NumPy code below repeats the comparison on synthetic data: $512 \times 512$ images with the origin at the top-left corner, 200 points, and 1 pixel of Gaussian noise. For $N$ = 8, 15, and 50 it fits $F$ 200 times with each version, all with the rank-2 correction, and reports the medians of $\log_{10}\kappa$ and of the paper’s error measure over all 200 matches.

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

rng = np.random.default_rng(0)

# Two cameras, 512 x 512 images, pixel origin at the top-left corner
K = np.array([[700.0, 0, 256], [0, 700.0, 256], [0, 0, 1]])
R = Rotation.from_euler("xyz", [2, -10, 3], degrees=True).as_matrix()
t = np.array([1.0, 0.2, 0.1])
X = rng.uniform([-2, -2, 5], [2, 2, 10], size=(200, 3))

def project(R, t):
    x = (X @ R.T + t) @ K.T
    return x / x[:, 2:]

x1_true, x2_true = project(np.eye(3), np.zeros(3)), project(R, t)

def isotropic(x):
    """Hartley (1997), Sec. 5.1: centroid to origin, mean distance sqrt(2)."""
    c = x[:, :2].mean(0)
    s = np.sqrt(2) / np.linalg.norm(x[:, :2] - c, axis=1).mean()
    return np.array([[s, 0, -s * c[0]], [0, s, -s * c[1]], [0, 0, 1]])

def non_isotropic(x):
    """Sec. 5.2: centroid to origin, both principal moments equal to one."""
    P = np.eye(3)[::-1]                           # sum u u^T = N U U^T with U upper
    U = P @ np.linalg.cholesky(P @ (x.T @ x / len(x)) @ P) @ P    # triangular
    return np.linalg.inv(U)                       # affine, since U[2] = (0, 0, 1)

def eight_point(x1, x2, normalize=None):
    T1 = normalize(x1) if normalize else np.eye(3)
    T2 = normalize(x2) if normalize else np.eye(3)
    y1, y2 = x1 @ T1.T, x2 @ T2.T
    A = np.einsum("ni,nj->nij", y2, y1).reshape(-1, 9)   # rows from x'^T F x = 0
    w = np.linalg.eigvalsh(A.T @ A)[::-1]                # eigenvalues, decreasing
    F = np.linalg.svd(A)[2][-1].reshape(3, 3)            # least eigenvector of A^T A
    U, S, Vt = np.linalg.svd(F)
    F = U @ np.diag([S[0], S[1], 0]) @ Vt                # closest rank-2 matrix
    return T2.T @ F @ T1, w[0] / w[7]                    # denormalize; d1/d8

def epipolar_error(F, x1, x2):
    """Mean point-to-epipolar-line distance in pixels, both images."""
    l2, l1 = x1 @ F.T, x2 @ F
    r = np.abs(np.sum(x2 * l2, 1))
    return np.mean(np.r_[r / np.hypot(l2[:, 0], l2[:, 1]), r / np.hypot(l1[:, 0], l1[:, 1])])

methods = {"unnormalized": None, "isotropic": isotropic, "non-isotropic": non_isotropic}
sigma = 1.0                                              # pixel noise
noise = lambda: np.c_[rng.normal(0, sigma, (len(X), 2)), np.zeros(len(X))]
print("   N  method          median log10(d1/d8)  median error (px)")
for n in (8, 15, 50):
    stats = {m: ([], []) for m in methods}
    for trial in range(200):
        x1, x2 = x1_true + noise(), x2_true + noise()
        idx = rng.choice(len(X), n, replace=False)       # N matches used to fit F
        for m, f in methods.items():
            F, k = eight_point(x1[idx], x2[idx], f)
            stats[m][0].append(np.log10(k))
            stats[m][1].append(epipolar_error(F, x1, x2))  # error over all 200 matches
    for m, (k, e) in stats.items():
        print(f"  {n:2d}  {m:<14}  {np.median(k):19.1f}  {np.median(e):17.2f}")

skew = lambda v: np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])
F_true = np.linalg.inv(K).T @ skew(t) @ R @ np.linalg.inv(K)     # F = K'^-T [t]x R K^-1
print(f"  true F error on noisy matches: {epipolar_error(F_true, x1_true + noise(), x2_true + noise()):.2f} px")

Output:

   N  method          median log10(d1/d8)  median error (px)
   8  unnormalized                   11.5              13.93
   8  isotropic                       4.6               6.53
   8  non-isotropic                   4.4               6.49
  15  unnormalized                   10.1               4.61
  15  isotropic                       3.1               1.66
  15  non-isotropic                   3.0               1.64
  50  unnormalized                    9.6               2.26
  50  isotropic                       2.7               1.19
  50  non-isotropic                   2.6               1.19
  true F error on noisy matches: 1.13 px

Normalization lowers $\kappa$ by about seven orders of magnitude. With 50 matches the normalized estimates come within 0.06 pixels of the true $F$’s error and the unnormalized one has nearly twice their error; the two scalings agree to a few hundredths of a pixel, as in the paper. At $N = 8$ the linear solutions with and without normalization coincide to rounding error (checked separately), so the halved error there comes entirely from the rank-2 correction.

Impact

Hartley and Zisserman’s Multiple View Geometry in Computer Vision (2004) credits Hartley (1997) with making normalization public in the computer vision literature (§4.9.1), calls normalization an essential step of the direct linear transform that must not be considered optional (§4.4.4), and applies it to homographies (§4.4) and camera resection (§7.1). Its normalized 8-point algorithm (§11.2) is the paper’s method, its comparison of algorithms (§11.5) uses the paper’s five image pairs, and its first recommendation (§11.5.1) is not to use the unnormalized version.

The Helmholtz Prize, awarded by the IEEE PAMI Technical Committee at each ICCV to ICCV papers from about ten years earlier that have had a significant impact on computer vision, went in 2013 to the 1995 conference paper “In Defence of the 8-Point Algorithm”, one of eleven papers so honored that year.

Limitations

• An algebraic cost. The algorithm still minimizes an algebraic residual, not a geometric error, and the paper does not claim it matches the best iterative methods.
• The rank-2 correction. Hartley acknowledges that minimizing the Frobenius distance to a singular matrix has little theoretical justification, citing his own work on minimizing algebraic distance as a better-founded alternative.
• An informal argument. The case for scaling and translation is given only an intuitive justification; the paper says a complete analysis of the condition number is too complex to undertake.
• Clean data. The matches were already cleaned of outliers; robust estimation, as in a RANSAC loop, is not addressed.
• Uncalibrated only. Evaluation in the calibrated case, for the essential matrix, is left as possible future work.

What Came After

Chojnacki, Brooks, van den Hengel, and Gawley (2003) gave a different explanation of why normalization helps. They showed that the normalized algorithm’s estimate, before the rank-2 correction, is the minimizer of a particular cost function, computable from a generalized eigenvalue problem, and argued that under a statistical model of the data this cost is better founded than the unnormalized algorithm’s, which augments Hartley’s conditioning argument. They place this in a line of statistical analyses by Torr, by Mühlich and Mester, and by Torr and Fitzgibbon.

In practice the normalized algorithm took the role the paper argued for: in the usual fundamental matrix pipeline it refits RANSAC inliers and initializes nonlinear refinement. For calibrated cameras, the minimal five-point algorithm became the usual solver for relative pose.