Structure from Motion

Structure from motion (SfM) recovers, from a set of photographs alone, the pose of each camera, usually its focal length and distortion, and a sparse cloud of 3D points on the scene. Each stage has its own article: feature detection and matching, two-view geometry with the essential or fundamental matrix, camera pose by perspective-n-point, triangulation, and bundle adjustment. This article shows how they combine and compares the two main strategies, incremental and global.

Problem Definition

Given $m$ images of a static scene, find cameras $P_i = K_i[\,R_i \mid \mathbf{t}_i\,]$ and points $\mathbf{X}_j$ that minimize the reprojection error over all observations $\mathbf{x}_{ij}$ of point $j$ in image $i$:

$$\min_{\{K_i, R_i, \mathbf{t}_i\},\, \{\mathbf{X}_j\}} \; \sum_{(i,j)} \rho\Big(\big\lVert \mathbf{x}_{ij} - \pi(K_i, R_i, \mathbf{t}_i, \mathbf{X}_j) \big\rVert^2\Big),$$

where $\pi$ projects a point into an image and $\rho$ is a robust loss. This is the bundle adjustment cost, but SfM must also find what bundle adjustment assumes: which observations belong together, and a starting point close enough to the minimum for a local method to reach it.

The solution is never unique. With known or self-calibrated intrinsics, cameras and points are determined only up to a 3D similarity transformation (7 degrees of freedom); with unconstrained intrinsics, up to a projective one (see multi-view geometry). Nothing in the images fixes the scale, so metric scale and a world frame need outside information, such as a known distance, GPS, or ground control points; Photo Tourism aligned its models to maps by an interactively chosen similarity transformation.

Input and Output

Input: an unordered image collection, such as internet photos of a landmark, or video frames, optionally with priors: EXIF focal lengths, calibrated intrinsics, GPS positions.

Output: the pose and intrinsics of every registered image, and for every reconstructed point its 3D position and track, the images and keypoints that observe it. Not every image is registered: on an internet collection of Rome, Schönberger and Frahm’s (2016) system registered about 21,000 of 74,000 images in its largest model. The points form a sparse point cloud; dense surfaces are computed afterwards from the recovered cameras by multi-view stereo.

Variants

• Unordered collections versus sequences. Video is ordered and moves little between frames, which eases matching, but short baselines mean tracks must span many frames to triangulate well (Hartley and Zisserman, 2004, §18.6). Real-time sequential variants are visual odometry and SLAM.
• Calibrated, uncalibrated, or self-calibrated cameras, single cameras or rigs, rigid or non-rigid scenes, and incremental, global, or hierarchical strategies.

Challenges

• Wrong matches that pass verification. Repetitive facades and symmetric buildings yield matches consistent with a wrong two-view geometry; Pan et al. (2024) report rotation averaging collapsing on a rotationally symmetric ETH3D scene. Watermarks, timestamps, and frames can link internet photos of different landmarks (Schönberger and Frahm, 2016), and one mismatch merges the tracks of two points.
• Drift. Each new camera is estimated from points that earlier cameras estimated, so errors accumulate; without repeated refinement, incremental SfM soon drifts into an unrecoverable state (Schönberger and Frahm, 2016).
• Degenerate geometry. Matches of a planar scene, or between images related by a pure rotation, fit a single homography and leave the epipolar geometry ill-determined, so the initial pair is chosen so that no homography explains its matches (Snavely et al., 2006). Nearly collinear motion, common in sequential data, leaves camera positions poorly determined by translation directions (Pan et al., 2024), and self-calibration leaves an ambiguity for critical motion sequences, such as pure translation or rotation about a single axis (Hartley and Zisserman, 2004, §19.11).
• Size. Exhaustive matching is quadratic in the number of images, and bundle adjustment dominates run time.

Approaches

Correspondence search

Features such as SIFT are detected and described in every image and matched between candidate pairs; for large collections, image retrieval chooses the candidates, as when Schönberger and Frahm matched each image against its 100 nearest neighbors from a vocabulary tree. Each pair is geometrically verified by fitting an essential matrix, fundamental matrix, or homography with RANSAC. Verified pairs are the edges of the view graph (or scene graph), and chaining matches through it produces the tracks.

Incremental SfM

Incremental SfM grows one reconstruction from a seed:

1. Initialization. Choose a pair with many matches and a wide baseline, estimate its relative pose from the essential matrix, for example with the five-point algorithm inside RANSAC, and triangulate the inliers. A bad seed may never be recovered from (Schönberger and Frahm, 2016).
2. Registration. Choose the next image, traditionally the one seeing the most triangulated points (COLMAP also rewards their even spread), and find its pose by PnP inside RANSAC. For uncalibrated images, Photo Tourism estimated the full camera matrix by the direct linear transform, which also gives an initial focal length.
3. Triangulation. Tracks now seen by two registered images become points, if their rays meet at a large enough angle and the points lie in front of the cameras.
4. Bundle adjustment and filtering. Refine, remove observations with large reprojection error, and repeat. COLMAP runs local bundle adjustment after each registration and global adjustment whenever the model has grown by a set percentage.

Because each new measurement is checked against an already reliable model, incremental SfM rejects outliers well; it pays with repeated bundle adjustments, dependence on the seed and registration order, and drift.

Global SfM

Global SfM estimates all cameras at once from the view graph’s two-view geometries, usually in two steps (Moulon et al., 2013). Rotation averaging finds absolute rotations $R_i$ consistent with the relative rotations, $R_{ij} \approx R_j R_i^\top$, by minimizing a robust sum of rotation distances over the edges. Translation averaging then finds camera centers $\mathbf{c}_i$ whose directions $(\mathbf{c}_j - \mathbf{c}_i)/\lVert \mathbf{c}_j - \mathbf{c}_i \rVert$ agree with the measured ones. This step is harder: relative translations are known only in direction, decomposing a two-view geometry needs accurate intrinsics, and collinear motion is degenerate (Pan et al., 2024). Points are then triangulated, and one global bundle adjustment refines everything.

Using the whole graph spreads error over loops instead of accumulating it along a chain, and one bundle adjustment replaces many. The weakness is outliers: no reliable partial model exists to check a wrong edge against. Wilson and Snavely (2014) therefore project the translation problem onto many 1D directions, where it becomes an ordering problem, minimum feedback arc set; edges that contradict the best orderings are removed as outliers (1DSfM). GLOMAP (Pan et al., 2024) replaces translation averaging with a joint global positioning of cameras and points, and its authors report accuracy on par with or better than COLMAP’s at a much lower run time.

Intrinsics, self-calibration, and priors

In principle, rigid motion alone determines constant intrinsics and yields a metric reconstruction, the auto-calibration of camera calibration (Hartley and Zisserman, 2004, ch. 19), but systems rely on simple camera models and priors. Photo Tourism modeled each camera by rotation, center, and one focal length, with the principal point at the image center, and initialized the focal length from EXIF tags when that value lay between 0.7 and 1.4 times the DLT estimate: EXIF values are usually good but occasionally off by more than a factor of two. COLMAP adds one radial distortion coefficient for internet photos, fixes the principal point because estimating it is ill-posed, and discards cameras whose field of view or distortion become implausible.

History

Ullman (1979) analyzed the problem as one of visual perception, built on what his abstract calls the structure-from-motion theorem: three orthographic views of four non-coplanar points determine their structure. Longuet-Higgins (1981) gave a linear algorithm for two calibrated views (his paper), and Tomasi and Kanade (1992) recovered all orthographic cameras and points at once by factorization (see multi-view geometry).

Snavely, Seitz and Szeliski’s Photo Tourism (2006) applied incremental SfM to photos from Flickr and Google, which its authors described as the first successful demonstration on such collections; their 2,635-photo Notre Dame set took about two weeks, mostly bundle adjustment. Building Rome in a Day (Agarwal et al., 2009) distributed matching over a cluster and reconstructed a reduced skeletal set of images before adding the rest; its 150,000 Rome images took 21 hours on 496 cores. COLMAP (Schönberger and Frahm, 2016) revised each incremental stage and was released as open source; Pan et al. (2024) call it the most widely used incremental system.

Learned components

Learned features and matchers can replace SIFT in the correspondence stage, leaving the geometry unchanged (see feature matching). VGGSfM (Wang et al., 2024) instead makes the whole pipeline differentiable: a deep point tracker produces tracks, a network predicts all cameras at once, and a differentiable bundle adjustment layer refines them. DUSt3R (Wang et al., 2024) regresses, for a pair of uncalibrated images, a 3D point for every pixel in a common frame and globally aligns many pairs; focal lengths and poses are read off the aligned points. Learning also depends on SfM: VGGSfM’s landmark model was trained on MegaDepth, labeled by SfM tools.

Example

The script reconstructs six synthetic cameras moving sideways past 400 points, each seen with 0.5-pixel noise by the cameras whose image it falls in. It initializes from cameras 0 and 1 with a linear eight-point essential matrix and the cheirality test, registers each further camera by linear PnP (the DLT), triangulates newly seen tracks, and runs a full bundle adjustment, with the Jacobian’s sparsity pattern, after each step. The data has no outliers, so RANSAC is omitted. Finally it aligns the result to the ground truth by a least-squares similarity (see iterative closest point).

import numpy as np
from scipy.optimize import least_squares
from scipy.sparse import lil_matrix
from scipy.spatial.transform import Rotation as Rot

rng = np.random.default_rng(3)
f, w, h, sigma = 800.0, 640, 480, 0.5       # focal length, image size, pixel noise

# Ground truth: a strip of points, and 6 cameras moving sideways past it while turning slightly.
X_true = rng.uniform([-1.5, -1.0, 4.0], [4.0, 1.0, 7.0], (400, 3))
R_true = [Rot.from_euler("yxz", [-0.05 * i, 0.02 * (-1) ** i, 0.0]).as_matrix() for i in range(6)]
C_true = [np.array([0.5 * i, 0.05 * np.sin(i), 0.1 * i]) for i in range(6)]
t_true = [-R @ C for R, C in zip(R_true, C_true)]

def project(R, t, X):                        # normalized image coordinates (K removed)
    Xc = X @ R.T + t
    return Xc[:, :2] / Xc[:, 2:]

# Observations: a point is seen where it projects inside the image, with noise.
obs = {}                                     # (camera, point) -> normalized observation
for i in range(6):
    xn = project(R_true[i], t_true[i], X_true)
    px = f * xn + [w / 2, h / 2]
    for j in np.flatnonzero((px[:, 0] > 0) & (px[:, 0] < w) & (px[:, 1] > 0) & (px[:, 1] < h)):
        obs[i, j] = xn[j] + sigma / f * rng.standard_normal(2)
seen = lambda i: {j for (c, j) in obs if c == i}

def triangulate(cams, j):                    # linear (DLT) triangulation of point j
    A = []
    for i in cams:
        R, t = pose[i]
        P = np.c_[R, t]
        x, y = obs[i, j]
        A += [x * P[2] - P[0], y * P[2] - P[1]]
    X = np.linalg.svd(np.array(A))[2][-1]
    return X[:3] / X[3]

# 1. Two-view initialization: essential matrix from cameras 0 and 1 (8-point, no outliers here).
common = sorted(seen(0) & seen(1))
a = np.array([np.r_[obs[0, j], 1] for j in common])
b = np.array([np.r_[obs[1, j], 1] for j in common])
E = np.linalg.svd(np.einsum("ni,nj->nij", b, a).reshape(-1, 9))[2][-1].reshape(3, 3)
U, _, Vt = np.linalg.svd(E)
U, Vt = U * np.linalg.det(U), Vt * np.linalg.det(Vt)
W = np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]])
best = -1
for R1 in (U @ W @ Vt, U @ W.T @ Vt):        # four (R, t) candidates: keep the one with
    for t1 in (U[:, 2], -U[:, 2]):           # most points in front of both cameras
        pose = {0: (np.eye(3), np.zeros(3)), 1: (R1, t1)}
        Xs = np.array([triangulate([0, 1], j) for j in common])
        n_front = np.sum((Xs[:, 2] > 0) & ((Xs @ R1.T + t1)[:, 2] > 0))
        if n_front > best:
            best, init = n_front, (R1, t1)
pose = {0: (np.eye(3), np.zeros(3)), 1: init}   # |t| = 1: the first baseline is the unit
points = {j: triangulate([0, 1], j) for j in common}

def bundle_adjust():
    cams, pts = sorted(pose), sorted(points)
    ci = {c: k for k, c in enumerate(cams)}
    pj = {p: k for k, p in enumerate(pts)}
    ob = [(ci[c], pj[p], x) for (c, p), x in obs.items() if c in ci and p in pj]
    ic, ip, xo = (np.array(v) for v in zip(*ob))
    def unpack(v):                           # camera 0 fixed, x of camera 1's t fixed (gauge)
        c = np.r_[np.zeros(6), v[:3], pose[1][1][0], v[3:6 * len(cams) - 7]].reshape(-1, 6)
        return c, v[6 * len(cams) - 7:].reshape(-1, 3)
    def res(v):
        c, X = unpack(v)
        Xc = Rot.from_rotvec(c[ic, :3]).apply(X[ip]) + c[ic, 3:]
        return (f * (Xc[:, :2] / Xc[:, 2:] - xo)).ravel()   # residuals in pixels
    c0 = np.concatenate([np.r_[Rot.from_matrix(pose[c][0]).as_rotvec(), pose[c][1]] for c in cams])
    v0 = np.r_[c0[6:9], c0[10:], np.array([points[p] for p in pts]).ravel()]
    rms0 = np.sqrt(np.mean(res(v0) ** 2) * 2)
    n = len(cams)                            # Jacobian sparsity: a residual depends only on
    S = lil_matrix((2 * len(ob), len(v0)), dtype=int)   # its camera's 6 and its point's 3 values
    for o in range(len(ob)):
        k = ic[o]
        cols = [] if k == 0 else [0, 1, 2, 3, 4] if k == 1 else range(6 * k - 7, 6 * k - 1)
        cols = [*cols, *range(6 * n - 7 + 3 * ip[o], 6 * n - 4 + 3 * ip[o])]
        S[2 * o: 2 * o + 2, cols] = 1
    sol = least_squares(res, v0, x_scale="jac", jac_sparsity=S)
    c, X = unpack(sol.x)
    for k, cam in enumerate(cams):
        pose[cam] = (Rot.from_rotvec(c[k, :3]).as_matrix(), c[k, 3:])
    for k, p in enumerate(pts):
        points[p] = X[k]
    return rms0, np.sqrt(np.mean(sol.fun ** 2) * 2)   # RMS per point (2D distance)

print(f"init  cameras 0,1: {len(points)} points, RMS before/after BA "
      "{:.3f} / {:.3f} px".format(*bundle_adjust()))

# 2. Incremental registration: PnP (linear DLT on normalized points), triangulation, BA.
for i in range(2, 6):
    js = sorted(seen(i) & set(points))
    A = []
    for j in js:
        Xh, (x, y) = np.r_[points[j], 1], obs[i, j]
        A += [np.r_[Xh, 0 * Xh, -x * Xh], np.r_[0 * Xh, Xh, -y * Xh]]
    P = np.linalg.svd(np.array(A))[2][-1].reshape(3, 4)
    U, s, Vt = np.linalg.svd(P[:, :3])
    P = P * np.sign(np.linalg.det(U @ Vt)) / s.mean()   # scale so that P = [R | t]
    pose[i] = (np.sign(np.linalg.det(U @ Vt)) * U @ Vt, P[:, 3])
    new = [j for j in seen(i) - set(points) if sum((c, j) in obs for c in pose) >= 2]
    for j in new:
        points[j] = triangulate([c for c in pose if (c, j) in obs], j)
    rms = bundle_adjust()
    print(f"add camera {i}: PnP from {len(js)} points, {len(new)} new points, "
          f"RMS before/after BA {rms[0]:.3f} / {rms[1]:.3f} px")

# 3. Compare with ground truth after a least-squares similarity alignment of the points.
pts = sorted(points)
A, B = np.array([points[j] for j in pts]), X_true[pts]
mu_a, mu_b = A.mean(0), B.mean(0)
U, s, Vt = np.linalg.svd((B - mu_b).T @ (A - mu_a))
D = np.diag([1, 1, np.sign(np.linalg.det(U @ Vt))])
Ra, scale = U @ D @ Vt, np.trace(np.diag(s) @ D) / np.sum((A - mu_a) ** 2)
align = lambda P: scale * (P - mu_a) @ Ra.T + mu_b          # reconstruction -> world
C_est = align(np.array([-R.T @ t for R, t in (pose[i] for i in range(6))]))
pos_err = np.linalg.norm(C_est - np.array(C_true), axis=1)
rot_err = [np.degrees(Rot.from_matrix(R_true[i] @ Ra @ pose[i][0].T).magnitude()) for i in range(6)]
n_obs = sum((c, j) in obs for c in range(6) for j in pts)
d = 6 * 6 - 7 + 3 * len(pts)                                  # free parameters
print(f"{len(pts)} of {len(X_true)} points, {n_obs} observations; final RMS expected "
      f"{sigma * np.sqrt(2 * (1 - d / (2 * n_obs))):.3f} px")
print(f"scale: 1 unit = {scale:.4f} m (true distance between cameras 0 and 1: "
      f"{np.linalg.norm(C_true[1] - C_true[0]):.4f} m)")
print("camera position error (mm):", np.array2string(1000 * pos_err, precision=2))
print("camera rotation error (deg):", np.array2string(np.array(rot_err), precision=3))
print(f"point error: median {1000 * np.median(np.linalg.norm(align(A) - B, axis=1)):.1f} mm")

Output:

init  cameras 0,1: 263 points, RMS before/after BA 1.406 / 0.351 px
add camera 2: PnP from 219 points, 45 new points, RMS before/after BA 1.239 / 0.464 px
add camera 3: PnP from 208 points, 57 new points, RMS before/after BA 1.990 / 0.511 px
add camera 4: PnP from 190 points, 23 new points, RMS before/after BA 0.714 / 0.538 px
add camera 5: PnP from 163 points, 2 new points, RMS before/after BA 1.208 / 0.553 px
390 of 400 points, 1560 observations; final RMS expected 0.555 px
scale: 1 unit = 0.5098 m (true distance between cameras 0 and 1: 0.5116 m)
camera position error (mm): [13.57 10.6  10.61 10.94 13.17 14.62]
camera rotation error (deg): [0.069 0.013 0.025 0.046 0.087 0.131]
point error: median 11.8 mm

Each linear step leaves up to two pixels of error, which bundle adjustment reduces to the noise level: the final 0.553 pixels matches the 0.555 expected at the optimum for this many observations and parameters (see reprojection error). The reconstruction’s unit is the first baseline, which only the alignment with ground truth reveals to be 0.51 m. Passing the sparsity pattern cuts the run time here from about 14 seconds to 1.

Datasets and Benchmarks

Ground truth for photo collections is scarce. The 1DSfM collections (Wilson and Snavely, 2014) are Flickr photos of landmarks such as Alamo, Notre Dame, and Piccadilly; their reference is a sequential reconstruction based on Bundler, the Photo Tourism system, placed roughly in meters with geotags, so they measure agreement, not true error. Schönberger and Frahm used them alongside Quad, which has ground-truth camera locations. ETH3D (Schöps et al., 2017) registers high-resolution DSLR images and multi-camera videos to laser scans of indoor and outdoor scenes. The PhotoTourism benchmark behind the Image Matching Challenge (Jin et al., 2021) evaluates the camera poses recovered from landmark photos against pseudo ground truth built by COLMAP.

Evaluation Metrics

• Mean reprojection error after the final bundle adjustment: necessary but not sufficient, since a wrong, self-consistent model can fit its own inliers well.
• Completeness: registered images, reconstructed points, and mean track length.
• Camera position error against ground truth after a similarity alignment, for example the median distance.
• Relative pose AUC: for every image pair, the larger of the rotation and translation-direction angular errors, summarized by the area under its cumulative curve up to thresholds such as 3°, 5°, and 10°. Pan et al. (2024) use it for unordered images; it needs no alignment.
• Run time, which differs by orders of magnitude between systems.

Applications

• Photo browsing: Photo Tourism used the recovered cameras to move between photos of a landmark in 3D.
• Cameras for dense reconstruction: multi-view stereo and novel-view synthesis take their cameras from SfM.
• Mapping and visual localization: new images are localized against an SfM model by PnP.
• Photogrammetry of buildings, terrain, and heritage sites, and training data for learned 3D models.

Open Problems

• Symmetric and repetitive scenes, where verified matches can be consistently wrong.
• Robust global positioning: translation averaging is still the weak step of global SfM.
• Hard imagery: low overlap, extreme viewpoint and illumination change, and forward motion.
• Learned reconstruction at scale: the VGGSfM authors note that their method cannot yet process thousands of images, as classical pipelines do.
• Dynamic scenes, where the static-scene assumption fails.