Triangulation

Triangulation recovers the position of a point in space from its images in two or more views whose cameras are known. Each image point defines a ray from its camera’s center through the pixel, and the 3D point lies where the rays meet. It is the step that turns matched pixels into geometry.

With perfect measurements the problem is trivial. With real measurements the rays almost never meet, so triangulation is an estimation problem: which point best explains the noisy image positions? Hartley and Sturm (1997) showed that the answer depends on what “best” means, and that the most common methods do not minimize the error that actually occurs, which is in the image.

Problem Definition

A point $\mathbf{X}$ is observed in $n \ge 2$ images. The camera matrices $P_i = K_i R_i\,[\,I \mid -\mathbf{C}_i\,]$ are known (see camera parameters), and the image measurements $\mathbf{x}_i$ are the true projections plus noise. The task is to estimate $\mathbf{X}$ .

The back-projected ray of measurement $\mathbf{x}_i$ , written in homogeneous coordinates as $\tilde{\mathbf{x}}_i$ , is

\mathbf{X}(\lambda) = \mathbf{C}_i + \lambda\, \mathbf{d}_i, \qquad \mathbf{d}_i = (K_i R_i)^{-1}\, \tilde{\mathbf{x}}_i, \quad \lambda > 0.

Two rays in space generically do not intersect. Equivalently, in two views the noisy pair does not satisfy the epipolar constraint $\mathbf{x}'^\top F\, \mathbf{x} = 0$ of the fundamental matrix, where $\mathbf{x} = \tilde{\mathbf{x}}_1$ and $\mathbf{x}' = \tilde{\mathbf{x}}_2$ , and the corresponding rays are skew.

Input and Output

Input: for each view, a camera matrix $P_i$ and an image point $\mathbf{x}_i$ , already corrected for lens distortion; typically one correspondence track per point, produced by feature matching or tracking.
Output: a 3D point, either Euclidean $(X, Y, Z)$ or homogeneous $(X_1, X_2, X_3, X_4)$ , optionally with a covariance and the residual reprojection errors.

Variants

Two-view and multi-view. Stereo rigs give two observations per point; a structure-from-motion track may have dozens.
Calibrated and uncalibrated. With known intrinsics, a reconstruction is Euclidean up to a similarity (an unknown scale if poses came from images alone). With only a fundamental matrix, the cameras are known only up to a projective transformation, and the triangulation method should give the same answer in every projective frame.

Challenges

Noise and parallax. Depth is determined by the angle between the rays, the parallax angle. When it is small, a tiny angular error moves the intersection a long way along the viewing direction.
Degenerate geometry. Points on the baseline, the line joining two camera centers, project to the epipoles and their depth cannot be determined from two views. Points near it, as in forward camera motion, are poorly constrained.
Outliers. A wrong correspondence gives a confident but meaningless point.
Cheirality. A valid point must lie in front of every camera that sees it, which algebraic methods do not enforce.

Approaches

Midpoint method

The classical geometric answer is to find the shortest segment between the two rays, by solving a small least-squares problem for $\lambda_1$ and $\lambda_2$ in $\mathbf{C}_1 + \lambda_1 \mathbf{d}_1 \approx \mathbf{C}_2 + \lambda_2 \mathbf{d}_2$ , and return its midpoint. It is cheap and intuitive, but it measures error in 3D space rather than in the images, so a camera close to the point and a camera far from it are treated as if their errors were the same size. Hartley and Sturm also noted that it is neither affine nor projectively invariant, since the common perpendicular is not an affine concept.

Linear triangulation (DLT)

Writing $\lambda \tilde{\mathbf{x}}_i = P_i \mathbf{X}$ and eliminating the scale with the third row gives two linear equations per view, with $\mathbf{p}_i^{k\top}$ the $k$ -th row of $P_i$ :

x_i\, \mathbf{p}_i^{3\top} \mathbf{X} - \mathbf{p}_i^{1\top} \mathbf{X} = 0, \qquad y_i\, \mathbf{p}_i^{3\top} \mathbf{X} - \mathbf{p}_i^{2\top} \mathbf{X} = 0.

Stacking them gives $A\mathbf{X} = \mathbf{0}$ with $A$ of size $2n \times 4$ , the direct linear transform used throughout multiple-view geometry (Hartley and Zisserman, 2004, §12.2). Hartley and Sturm name two ways to solve it in the least-squares sense:

Homogeneous (their “Linear-Eigen”): minimize $\lVert A\mathbf{X} \rVert$ subject to $\lVert \mathbf{X} \rVert = 1$ ; the solution is the right singular vector of $A$ with the smallest singular value. It can represent points at or near infinity, where $X_4 \approx 0$ .
Inhomogeneous (their “Linear-LS”): set $X_4 = 1$ and solve the resulting $2n \times 3$ system by ordinary least squares. It cannot represent points at infinity.

Both minimize an algebraic residual with no direct geometric meaning: scaling a row of $A$ changes the answer. Hartley and Sturm showed that the inhomogeneous form is affine invariant, since an affine transformation preserves $X_4 = 1$ , while the homogeneous form is neither affine nor projectively invariant. Their iterative variant reweights each view’s equations by $1/(\mathbf{p}_i^{3\top}\mathbf{X})$ from the previous estimate, so that the residuals approach image distances; it improved on the linear methods but occasionally failed to converge near the epipoles.

Minimizing reprojection error

If image noise is independent, isotropic, and Gaussian, the maximum-likelihood estimate minimizes the reprojection error:

\hat{\mathbf{X}} = \arg\min_{\mathbf{X}} \sum_{i=1}^{n} \big\lVert \mathbf{x}_i - \pi(P_i \mathbf{X}) \big\rVert^2,

where $\pi$ divides by the third coordinate. This cost is measured in the images, so it is invariant to any projective transformation applied to the cameras and the point together.

For two views, Hartley and Sturm (1997) found the global minimum in closed form. They recast the problem as finding the pair of corrected points $\hat{\mathbf{x}}_1, \hat{\mathbf{x}}_2$ closest to the measurements that satisfy the epipolar constraint exactly. Parametrizing the pencil of epipolar lines by $t$ , they write the summed squared distance from each measured point to its corresponding epipolar line as a rational function of $t$ . Setting its derivative to zero gives a polynomial of degree 6, whose real roots, together with $t = \infty$ , are the candidates. Once the points are corrected, their rays meet exactly and any method gives the point. Their paper also shows that the cost can have up to three local minima, one of their examples having a spurious minimum even for a perfect match, so a local iterative search from a poor start can stop in the wrong one. Hartley and Zisserman (2004, §12.4–12.5) present the same method alongside the cheaper Sampson approximation, a first-order correction of the points toward the epipolar constraint.

With more than two views there is no comparable closed form, and in practice the linear estimate initializes a Gauss–Newton or Levenberg–Marquardt minimization of the reprojection error. With three unknowns per point, each iteration solves a $3 \times 3$ system, and a few iterations suffice from a good start. In structure from motion the same cost is minimized jointly over all points and cameras by bundle adjustment, for which per-point triangulation provides the starting points.

L∞ triangulation

The sum of squared reprojection errors can have several local minima in multiple views. Hartley and Schaffalitzky (2004) proposed minimizing instead the maximum reprojection error over the views, the $L_\infty$ norm, and showed that for multiview triangulation this cost has a single local minimum on a convex domain, so its global optimum can be found reliably with standard constrained optimization. Kahl (2005) cast triangulation, camera resectioning, and homography estimation in this form as quasi-convex problems and solved them by bisection with second-order cone programming. The price is sensitivity to outliers: the maximum error is set by the single worst measurement, so $L_\infty$ methods need clean data or an explicit outlier-removal step.

Accuracy and Uncertainty

For a rectified stereo pair with focal length $f$ in pixels and baseline $B$ , a point at depth $Z$ has disparity $d = fB/Z$ (see stereo vision). Differentiating gives the first-order depth error

\delta Z \approx \frac{Z^2}{f B}\, \delta d.

Depth uncertainty grows with the square of the distance. A one-pixel disparity error that is negligible at 2 m can be meters at 50 m. More generally, the 3D uncertainty is an ellipsoid elongated along the viewing direction, and its elongation is governed by the parallax angle between the rays; Hartley and Zisserman (2004, §12.6) discuss the probability distribution of the triangulated point. A first-order covariance follows from the Jacobian $J$ of the reprojection residuals at the solution, $\Sigma_{\mathbf{X}} \approx \sigma^2 (J^\top J)^{-1}$ . When the parallax is very small, depth error is no longer approximately linear in the image noise, and the first-order prediction underestimates it.

Cheirality and Robustness

The homogeneous DLT solution is defined only up to sign, and noise can place a point behind a camera. Checking that the depth $\mathbf{p}_i^{3\top}\mathbf{X}$ , with consistently signed $\mathbf{X}$ and camera matrices, is positive in every view is the cheirality test. It rejects bad points, and it is also how relative pose is disambiguated: decomposing an essential matrix yields four candidate camera pairs, and only one places triangulated points in front of both cameras.

Outlier correspondences are normally removed beforehand, during geometric verification with RANSAC. Afterwards, points are commonly filtered by maximum reprojection error, minimum parallax angle, and cheirality; for long tracks, triangulating from subsets of views and keeping the views consistent with the best estimate removes residual mismatches.

Example

The NumPy code below triangulates 500 noisy points with the midpoint method, both DLT variants, and five Gauss–Newton iterations on the reprojection error, started from the homogeneous DLT. The first uses a converging pair with the second camera closer to the scene; the second uses a parallel pair and shrinks the baseline, comparing the root-mean-square depth error with the first-order prediction $Z^2 \sqrt{2}\,\sigma / (fB)$ , since the disparity combines the independent noise of two image coordinates.

import numpy as np

rng = np.random.default_rng(0)
f, sigma = 800.0, 1.0                       # focal length (px), pixel noise (px)
K = np.array([[f, 0, 320], [0, f, 240], [0, 0, 1.0]])

def camera(C, target=None):
    """P = K R [I | -C]; R = I, or R turns the optical axis toward target."""
    R = np.eye(3)
    if target is not None:
        z = (target - C) / np.linalg.norm(target - C)
        x = np.cross([0.0, 1.0, 0.0], z); x /= np.linalg.norm(x)
        R = np.vstack([x, np.cross(z, x), z])
    return K @ R @ np.hstack([np.eye(3), -C[:, None]])

def project(P, X):
    x = np.hstack([X, np.ones((len(X), 1))]) @ P.T
    return x[:, :2] / x[:, 2:]

def midpoint(xs, Ps):
    """Midpoint of the shortest segment between the two back-projected rays."""
    C = [-np.linalg.solve(P[:, :3], P[:, 3]) for P in Ps]          # centres
    d = [np.linalg.solve(P[:, :3], np.append(x, 1.0)) for x, P in zip(xs, Ps)]
    a = np.linalg.lstsq(np.column_stack([d[0], -d[1]]), C[1] - C[0], rcond=None)[0]
    return 0.5 * (C[0] + a[0] * d[0] + C[1] + a[1] * d[1])

def dlt(xs, Ps, homogeneous=True):
    """Linear triangulation: two equations per view, A X = 0."""
    A = np.vstack([[x[0] * P[2] - P[0], x[1] * P[2] - P[1]] for x, P in zip(xs, Ps)])
    if homogeneous:                          # min ||A X|| subject to ||X|| = 1
        X = np.linalg.svd(A)[2][-1]
        return X[:3] / X[3]
    return np.linalg.lstsq(A[:, :3], -A[:, 3], rcond=None)[0]   # fix X4 = 1

def gauss_newton(X, xs, Ps, iters=5):
    """Minimize the sum of squared reprojection errors over the 3D point."""
    for _ in range(iters):
        r, J = [], []
        for x, P in zip(xs, Ps):
            h = P @ np.append(X, 1.0)
            r.append(h[:2] / h[2] - x)
            J.append((P[:2, :3] * h[2] - np.outer(h[:2], P[2, :3])) / h[2] ** 2)
        X = X - np.linalg.lstsq(np.vstack(J), np.concatenate(r), rcond=None)[0]
    return X

def experiment(Ps, n=500):
    X = np.column_stack([rng.uniform(-1, 1, n), rng.uniform(-1, 1, n),
                         rng.uniform(4, 6, n)])          # points 4-6 m away
    obs = [project(P, X) + rng.normal(0, sigma, (n, 2)) for P in Ps]
    est = {"midpoint": [], "DLT (SVD)": [], "DLT (X4 = 1)": [], "Gauss-Newton": []}
    for xs in zip(*obs):
        est["midpoint"].append(midpoint(xs, Ps))
        est["DLT (SVD)"].append(dlt(xs, Ps))
        est["DLT (X4 = 1)"].append(dlt(xs, Ps, homogeneous=False))
        est["Gauss-Newton"].append(gauss_newton(est["DLT (SVD)"][-1], xs, Ps))
    return X, obs, {k: np.array(v) for k, v in est.items()}

# 1. Converging cameras at different distances from the scene
Ps = [camera(np.zeros(3)), camera(np.array([1.5, 0.0, 3.0]), np.array([0.0, 0.0, 5.0]))]
X, obs, est = experiment(Ps)
print("converging pair, noise 1 px, 500 points")
print(f"{'method':14s} {'mean 3D error (mm)':>19s} {'RMS reprojection (px)':>22s}")
for k, E in est.items():
    rep = np.concatenate([project(P, E) - x for P, x in zip(Ps, obs)])
    print(f"{k:14s} {1000 * np.linalg.norm(E - X, axis=1).mean():19.1f}"
          f" {np.sqrt(2 * np.mean(rep ** 2)):22.3f}")

# 2. Parallel (rectified) pair: depth error versus baseline B
print("\nparallel pair, Gauss-Newton estimate")
print(f"{'baseline (m)':>12s} {'RMS depth error (mm)':>21s} {'predicted (mm)':>15s}")
for B in [1.0, 0.5, 0.2, 0.1, 0.05]:
    Ps = [camera(np.zeros(3)), camera(np.array([B, 0.0, 0.0]))]
    X, obs, est = experiment(Ps)
    dZ = est["Gauss-Newton"][:, 2] - X[:, 2]
    pred = np.sqrt(np.mean((X[:, 2] ** 2 * np.sqrt(2) * sigma / (f * B)) ** 2))
    print(f"{B:12.2f} {1000 * np.sqrt(np.mean(dZ ** 2)):21.1f} {1000 * pred:15.1f}")

Output:

converging pair, noise 1 px, 500 points
method          mean 3D error (mm)  RMS reprojection (px)
midpoint                      10.8                  0.979
DLT (SVD)                     10.8                  0.960
DLT (X4 = 1)                  10.8                  0.960
Gauss-Newton                  10.5                  0.731

parallel pair, Gauss-Newton estimate
baseline (m)  RMS depth error (mm)  predicted (mm)
        1.00                  45.1            46.1
        0.50                  88.4            91.5
        0.20                 237.6           229.1
        0.10                 458.6           457.8
        0.05                1083.2           933.0

Gauss–Newton cuts the reprojection error substantially, to near the value expected at the optimum: four image coordinates constrain three unknowns, leaving one degree of freedom of residual per point, so the root-mean-square image residual is about $\sigma/\sqrt{2} \approx 0.71$ pixels. The 3D error improves only slightly: minimizing image error is not the same as minimizing 3D error, and in a well-conditioned configuration the simple methods are already close. The second table shows the baseline effect: halving $B$ roughly doubles the depth error, as $Z^2/(fB)$ predicts, until at 5 cm the disparity of about 8 pixels is so small that the error becomes visibly nonlinear and exceeds the first-order prediction.

Evaluation Metrics

Synthetic experiments measure 3D error against ground truth; real data, which rarely has it, is judged by reprojection error in pixels, often with the parallax angle and number of views per point. The two can disagree: in Hartley and Sturm’s synthetic Euclidean experiments, the midpoint and linear methods had much larger image residuals than the polynomial method yet comparable or slightly smaller 3D errors; after a projective transformation of the cameras, the same methods were far worse in both measures.

Applications

Stereo depth. Every disparity map from stereo vision is a dense triangulation in the rectified geometry.
Structure from motion and multi-view reconstruction. New points are triangulated as tracks grow, then refined with the cameras by bundle adjustment.
Visual SLAM and odometry. Map points are triangulated from keyframes, with parallax and cheirality checks deciding when a point is reliable enough to add.
Motion capture and photogrammetry. Calibrated multi-camera systems triangulate markers or targets in every frame.

Open Problems

Globally optimal $L_2$ triangulation from many views is far more expensive than the two-view polynomial method, so practical systems rely on a linear start and local refinement, which assumes the start lies in the right basin. Other concerns are reliable uncertainty at low parallax, robust triangulation of long tracks with some mismatches, and the fact that the cameras are themselves estimates.