Five-Point Algorithm

The five-point algorithm computes the relative pose of two calibrated cameras from the smallest possible number of point matches. Given five correspondences in normalized image coordinates, it finds every essential matrix $E$ with $\hat{\mathbf{x}}'^\top E\,\hat{\mathbf{x}} = 0$ for all five: up to ten real solutions, each giving a rotation and a translation direction. Nistér (2003, 2004) gave an efficient method that reduces the problem directly to the roots of a single polynomial of degree ten, matching the number of solutions. Run inside RANSAC, it is the standard way to estimate the motion between two images with known intrinsics, in visual odometry, structure from motion, and SLAM initialization. The essential matrix article covers $E = [\mathbf{t}]_\times R$ and pose recovery; this article covers the solver.

Problem

The relative pose of two calibrated views, a rotation and a translation whose length cannot be recovered, has five degrees of freedom. Each match gives one equation, $\hat{\mathbf{x}}'^\top E\,\hat{\mathbf{x}} = 0$, so five matches are the minimum that leaves finitely many solutions, and the five-point problem asks for all of them. In photogrammetry, where it is called relative orientation, iterative adjustment needs an approximate starting configuration; a hypothesize-and-test loop needs a direct solver that returns every solution (Stewénius et al., 2006).

Inputs and Outputs

• Input: five correspondences $\hat{\mathbf{x}}_i \leftrightarrow \hat{\mathbf{x}}'_i$ in normalized coordinates, $\hat{\mathbf{x}} = K^{-1}\mathbf{x}$, with lens distortion removed.
• Output: up to ten essential matrices, each up to scale, and for each the one $(R, \mathbf{t})$ candidate, if any, that puts the points in front of both cameras.

Intuition

The five epipolar equations are linear in the nine entries of $E$, so they leave a four-dimensional subspace of $3 \times 3$ matrices. Every matrix in it satisfies the five matches, but almost none is an essential matrix; the essential matrices form a curved set defined by cubic equations, and the solutions are where the two meet. Fixing the scale leaves three unknowns and ten cubic equations. The rest is elimination, until a single polynomial in one unknown remains.

Algorithm

Following Nistér (2004):

1. Compute a basis $X, Y, Z, W$ of the null space of the $5 \times 9$ linear system.
2. Write $E = xX + yY + zZ + W$, fixing the scale.
3. Substitute $E$ into the ten cubic constraints on essential matrices.
4. Reduce the resulting $10 \times 20$ coefficient matrix by Gauss–Jordan elimination.
5. Form a $3 \times 3$ matrix $B(z)$ of polynomials in $z$ from its rows; $\det B(z)$ has degree ten.
6. For each real root $z$, recover $x$ and $y$ from $B(z)$ and form $E$.
7. Decompose each $E$ into four pose candidates and keep the one that passes the cheirality test, triangulating the points and checking that they lie in front of both cameras.

Mathematical Formulation

The linear step. Each match gives $\tilde{\mathbf{q}}^\top \tilde{\mathbf{e}} = 0$, where $\tilde{\mathbf{e}}$ holds the entries of $E$ row by row and $\tilde{\mathbf{q}} = \hat{\mathbf{x}}' \otimes \hat{\mathbf{x}}$. For points in general position the five rows have rank five, so the null space is four-dimensional and

$$E = xX + yY + zZ + W .$$

The cubic constraints. A real nonzero $3 \times 3$ matrix is an essential matrix if and only if

$$2\,EE^\top E - \operatorname{tr}\!\left(EE^\top\right) E = 0,$$

which is the form used by Stewénius et al. (2006) and Li and Hartley (2006); Nistér writes the same constraint with a factor $\tfrac{1}{2}$ on the trace term. Together with the rank condition

$$\det E = 0,$$

these are ten cubic equations in $(x, y, z)$. For real matrices the trace constraint alone already characterizes essential matrices, rank included; the determinant adds a tenth equation that the elimination below uses.

Elimination. Each cubic equation combines the 20 monomials of degree at most three. Nistér orders the columns so that ten monomials, $x^3, y^3, x^2y, xy^2, x^2z, x^2, y^2z, y^2, xyz, xy$, come first, and the rest are grouped as $x$, $y$, and $1$ times polynomials in $z$. Gauss–Jordan elimination makes the first ten columns the identity. Six of the reduced rows, labeled $\langle e \rangle$ to $\langle j \rangle$, have leading monomials $x^2z$, $x^2$, $y^2z$, $y^2$, $xyz$, and $xy$, so the differences

$$\langle k \rangle = \langle e \rangle - z\langle f \rangle, \quad \langle l \rangle = \langle g \rangle - z\langle h \rangle, \quad \langle m \rangle = \langle i \rangle - z\langle j \rangle$$

cancel the leading terms and contain only $x$, $y$, and $1$ multiplied by polynomials in $z$ of degree three, three, and four. Written as $B(z)\,(x, y, 1)^\top = \mathbf{0}$, the $3 \times 3$ matrix $B(z)$ must be singular at a solution, so

$$\langle n \rangle = \det B(z) = 0$$

is a polynomial of degree ten in $z$. For each real root, $(x, y, 1)$ is proportional to the cross product of two rows of $B(z)$, which gives $x$ and $y$ in closed form. The conference version (Nistér, 2003) used the nine trace equations only and reached degree ten by combining two $4 \times 4$ determinants of degree eleven; the journal version adds $\det E = 0$ and needs only the $3 \times 3$ matrix.

How many solutions. Nistér (2004) summarizes the history: Kruppa showed in 1913 that there are at most eleven solutions, and later work, which he cites to Demazure (1988), Faugeras and Maybank (1990), and others, showed that there are at most ten, and ten in general counting complex ones. Kruppa’s procedure required intersecting two sextic curves, and Nistér describes his method as a refinement of Philip’s (1996) derivation, which led to a polynomial of degree thirteen.

Complexity

The solver has no tuning parameters, and its work per sample is fixed: a small null-space computation, elimination on a $10 \times 20$ matrix, and the roots of a tenth-degree polynomial, which Nistér brackets with Sturm sequences and a companion-matrix eigenvalue solve finds more simply. He reported about 110–140 μs per call for two views on a 550 MHz machine, including overhead.

The payoff is in RANSAC. With confidence $p$ and inlier fraction $w$, the number of samples is $N = \log(1-p) / \log(1-w^s)$. At $p = 0.99$ and $w = 0.3$, this gives 1,893 samples for $s = 5$, against 21,055 for the seven-point algorithm and 70,188 for eight points (see the table in the RANSAC article). Each sample yields several hypotheses to score: in Nistér’s random trials the tenth-degree polynomial had 4.55 real roots on average, and 2.74 hypotheses remained after the cheirality test.

Implementation

The NumPy code below follows steps 1–7, with polynomial products as 3D convolutions of coefficient arrays and roots from np.roots, which uses the companion matrix. It solves from five noise-free matches, keeps the candidates that put those points in front of both cameras, and scores each on 25 further matches, for a general and a planar scene.

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

rng = np.random.default_rng(1)

def pmul(a, b):
    """Product of two polynomials in x, y, z stored as 4x4x4 coefficient arrays."""
    return convolve(a, b)[:4, :4, :4]

# Columns of Nister's 10 x 20 matrix: ten leading monomials, then x[2], y[2], [3] in z
MONOS = ["xxx", "yyy", "xxy", "xyy", "xxz", "xx", "yyz", "yy", "xyz", "xy",
         "xzz", "xz", "x", "yzz", "yz", "y", "zzz", "zz", "z", ""]
IDX = [tuple(m.count(c) for c in "xyz") for m in MONOS]

def five_point(x1, x2):
    """Essential matrices with x2^T E x1 = 0 for five normalized matches (Nister 2004)."""
    Q = np.einsum("ni,nj->nij", x2, x1).reshape(5, 9)         # 5 x 9 linear system
    X, Y, Z, W = np.linalg.svd(Q)[2][5:].reshape(4, 3, 3)       # basis of its null space
    E = np.zeros((3, 3, 4, 4, 4))                               # E = xX + yY + zZ + W
    E[..., 1, 0, 0], E[..., 0, 1, 0], E[..., 0, 0, 1], E[..., 0, 0, 0] = X, Y, Z, W
    det = sum(pmul(pmul(E[0, (i + 1) % 3], E[1, (i + 2) % 3])
                   - pmul(E[0, (i + 2) % 3], E[1, (i + 1) % 3]), E[2, i]) for i in range(3))
    EEt = np.array([[sum(pmul(E[i, k], E[j, k]) for k in range(3)) for j in range(3)]
                    for i in range(3)])
    L = 2 * EEt - np.trace(EEt) * np.eye(3)[..., None, None, None]  # 2EE^T - tr(EE^T) I
    cubics = [sum(pmul(L[i, k], E[k, j]) for k in range(3)) for i in range(3) for j in range(3)]
    A = np.array([[c[m] for m in IDX] for c in [det] + cubics])  # 10 equations, 20 monomials
    A = np.linalg.solve(A[:, :10], A)                             # Gauss-Jordan elimination
    # Rows e..j lead with x^2 z, x^2, y^2 z, y^2, xyz, xy: k = e - z f, l = g - z h, m = i - z j
    def coeffs(r):  # coefficients of x, y and 1 in row r, as polynomials in z
        return [A[r, 10:13], A[r, 13:16], A[r, 16:20]]
    B = [[np.polysub(p, np.polymul([1, 0], q)) for p, q in zip(coeffs(r), coeffs(r + 1))]
         for r in (4, 6, 8)]                                      # 3 x 3 matrix B(z)
    p1 = np.polysub(np.polymul(B[0][1], B[1][2]), np.polymul(B[0][2], B[1][1]))
    p2 = np.polysub(np.polymul(B[0][2], B[1][0]), np.polymul(B[0][0], B[1][2]))
    p3 = np.polysub(np.polymul(B[0][0], B[1][1]), np.polymul(B[0][1], B[1][0]))
    n = sum(np.polymul(p, b) for p, b in zip((p1, p2, p3), B[2]))  # det B(z), degree 10
    roots = np.roots(n)
    Es = []
    for z in roots[np.abs(roots.imag) < 1e-9].real:
        x, y = np.polyval(p1, z) / np.polyval(p3, z), np.polyval(p2, z) / np.polyval(p3, z)
        e = x * X + y * Y + z * Z + W
        Es.append(e / np.linalg.norm(e))
    return len(np.trim_zeros(n, "f")) - 1, Es

def decompose(E):
    """Four (R, t) candidates from the SVD of E."""
    U, _, Vt = np.linalg.svd(E)
    U, Vt = U * np.sign(np.linalg.det(U)), Vt * np.sign(np.linalg.det(Vt))
    D = np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]])
    return [(R, s * U[:, 2]) for R in (U @ D @ Vt, U @ D.T @ Vt) for s in (1, -1)]

def depths(R, t, a, b):
    """Depths of the triangulated point in both cameras (linear triangulation)."""
    P1, P2 = np.eye(3, 4), np.c_[R, t]
    A = np.array([a[0] * P1[2] - P1[0], a[1] * P1[2] - P1[1],
                  b[0] * P2[2] - P2[0], b[1] * P2[2] - P2[1]])
    X = np.linalg.svd(A)[2][-1]
    X = X[:3] / X[3]
    return X[2], (R @ X + t)[2]

def poses(Es, x1, x2):
    """Keep, for each E, the candidate that puts all five points in front of both cameras."""
    return [(R, t) for E in Es for R, t in decompose(E)
            if all(min(depths(R, t, a, b)) > 0 for a, b in zip(x1, x2))]

def sampson(E, x1, x2):
    """Squared Sampson distance of each match to the epipolar constraint."""
    Ex1, Etx2 = x1 @ E.T, x2 @ E
    return np.sum(x2 * Ex1, 1) ** 2 / (Ex1[:, 0]**2 + Ex1[:, 1]**2 + Etx2[:, 0]**2 + Etx2[:, 1]**2)

# Ground truth: X' = R X + t, normalized coordinates x_hat = K^-1 x
R_true = Rotation.from_rotvec(np.radians(10) * np.array([0.3, 0.9, 0.3])).as_matrix()
t_true = np.array([1.0, 0.1, 0.2])

def errors(R, t):
    rot = np.degrees(np.linalg.norm(Rotation.from_matrix(R.T @ R_true).as_rotvec()))
    cos = t @ t_true / np.linalg.norm(t_true)
    return rot, np.degrees(np.arccos(np.clip(cos, -1, 1)))

for name in ("general scene", "planar scene"):
    XY = rng.uniform(-2, 2, (30, 2))
    Zw = 6 - 0.3 * XY[:, 0] if name == "planar scene" else rng.uniform(4, 8, 30)
    Xw = np.c_[XY, Zw]                          # 30 points in front of camera 1
    x1 = Xw / Xw[:, 2:]
    Xc2 = Xw @ R_true.T + t_true
    x2 = Xc2 / Xc2[:, 2:]
    degree, Es = five_point(x1[:5], x2[:5])
    cands = poses(Es, x1[:5], x2[:5])
    print(f"{name}: degree {degree}, {len(Es)} real roots, {len(cands)} pass cheirality")
    for R, t in cands:
        E = np.cross(t, R.T).T                  # [t]x R
        d = np.sqrt(np.median(sampson(E, x1[5:], x2[5:])))   # on the 25 unused points
        front = sum(min(depths(R, t, a, b)) > 0 for a, b in zip(x1, x2))
        print("  Sampson %.0e, in front %2d/30, R error %5.2f deg, t error %6.2f deg"
              % (d, front, *errors(R, t)))

Output:

general scene: degree 10, 2 real roots, 1 pass cheirality
  Sampson 4e-13, in front 30/30, R error  0.00 deg, t error   0.00 deg
planar scene: degree 10, 6 real roots, 5 pass cheirality
  Sampson 3e-03, in front 30/30, R error 10.02 deg, t error 100.28 deg
  Sampson 2e-14, in front 30/30, R error  0.00 deg, t error   0.00 deg
  Sampson 8e-16, in front 30/30, R error  7.77 deg, t error  49.48 deg
  Sampson 3e-03, in front 30/30, R error 11.68 deg, t error  84.53 deg
  Sampson 2e-03, in front 30/30, R error  7.48 deg, t error  42.54 deg

For the general scene, two of the ten roots are real, one survives the cheirality test, and it recovers the pose to rounding error. For the planar scene, five candidates pass the test on the sample. The other 25 matches reject three, with Sampson distances near $3 \times 10^{-3}$ in normalized units (about 2 pixels at a focal length of 800 pixels), but two fit every point exactly: the true pose and a second one, 7.8° off in rotation and 49° in translation direction, with all 30 points in front of both cameras. This is the planar two-fold ambiguity; here every point is closer to the first camera than to the second.

Properties and Behavior

• Up to ten physical solutions. Nistér found configurations with ten distinct real solutions, all with the points in front of both cameras, about twice in $10^7$ random trials. Extra matches, a third view, or the RANSAC score choose among the solutions.
• Planar scenes. Coplanar points are degenerate for the fundamental matrix, leaving a family of solutions. For the five-point method they are not: Nistér (2004) states that two calibrated views of a plane give two essential matrices, or one when the baseline is perpendicular to the plane. Cheirality removes the false one unless all visible points are closer to one camera than the other, and a third view generally resolves it. Stewénius et al. (2006) note that the method likewise keeps working for other critical surfaces: the true solution is among those returned.
• Accuracy under noise. In Nistér’s synthetic experiments, the five-point method gave much better translation directions than the six-, seven-, and eight-point methods for sideways motion with a short baseline and a shallow scene. For forward motion the others improved greatly, the eight-point method most. The five-point method’s rotation errors were larger than the others’, though only hundredths of a degree.
• Numerical accuracy. With exact data, Nistér reported median errors in $E$ near $10^{-10}$ with Sturm bracketing and $10^{-14}$ with the companion matrix, plus a small fraction of much larger errors.

Limitations

• Calibration is required. Nistér’s experiments suggest the focal length should be within about 10% for any benefit over uncalibrated methods.
• Pure rotation. If $\mathbf{t} = \mathbf{0}$, then $\hat{\mathbf{x}}' \propto R\hat{\mathbf{x}}$, and every matrix $[\mathbf{t}]_\times R$, for any $\mathbf{t}$, satisfies all the epipolar constraints. The solutions form a continuum and the translation direction is undefined; very short baselines approach this case. A homography or a rotation-only model fits such motion instead.
• Minimal samples amplify noise. A solution from five noisy matches can be far from the best fit to all inliers, so it is refined on the inlier set, typically by minimizing a geometric error and later by bundle adjustment.

Variants

• More than five points. Nistér proposed using the four right singular vectors with the smallest singular values as $X, Y, Z, W$, a least-squares variant that still enforces the constraints exactly.
• Gröbner-basis solver. Stewénius, Engels, and Nistér (2006) showed that Gauss–Jordan elimination of the same $10 \times 20$ matrix, with the monomials in graded lexicographic order, yields a Gröbner basis. From it they build a $10 \times 10$ action matrix for multiplication by $x$, whose eigenvectors give all solutions at once. They reported that this removed the tail of large numerical errors of the original.
• Hidden-variable resultant. Li and Hartley (2006) treat $z$ as a parameter, write the ten equations as $C(z)\,\mathbf{m}(x, y) = \mathbf{0}$ over the ten monomials in $x$ and $y$ up to degree three, and solve $\det C(z) = 0$, which they found to be exactly of degree ten. The null vector of $C(z)$ then gives $x$ and $y$.
• Three views. Nistér (2004) solves on the first and last of three views, triangulates the five points, computes the middle camera with a three-point perspective-n-point solver, and keeps the hypothesis that best reprojects the five points, a unique solution per sample even for a planar scene.

Minimal solvers in general. The five-point algorithm is one instance of a broader approach: solve each geometric problem from the smallest data set that gives finitely many solutions, so that RANSAC samples as few points as possible. Other examples are the seven-point fundamental matrix and P3P, and relative pose with an unknown common focal length from six points, which has 15 solutions in general and for which Stewénius, Nistér, Kahl, and Schaffalitzky gave an action-matrix solver in 2005. Kukelova, Bujnak, and Pajdla (2008) automated the construction, with a generator that produces Gröbner-basis solvers from the polynomial equations.