Structure from Motion

Structure from Motion is the problem of estimating, from a set of 2D image correspondences across $n$ images, the 3D positions $\mathbf{X}_j \in \mathbb{R}^3$ of $m$ scene points and the pose $(\mathbf{R}_i, \mathbf{t}_i)$ (and often intrinsics $\mathbf{K}_i$) of each camera $i$, such that projecting each point into each camera that observes it reproduces the measured 2D location as closely as possible. Formally, if $\pi(\cdot)$ is the pinhole projection function, we seek parameters minimizing the reprojection error

where $\mathcal{O}$ is the set of observed (camera, point) pairs and $\mathbf{x}_{ij}$ is the measured pixel location of point $j$ in image $i$. This joint nonlinear least-squares problem, solved over all cameras and points at once, is called bundle adjustment (BA) — named for the "bundle" of light rays from each 3D point through the camera centers that must be adjusted to intersect correctly. Two structural facts define the character of the problem. First, without external metric information (a known baseline, GPS, or a calibration object), the entire reconstruction is recoverable only up to a similarity transform — a global rotation, translation, and uniform scale — because scaling the whole scene and all camera translations together leaves every image projection unchanged. This is the scale ambiguity inherited directly from monocular Epipolar Geometry, where the essential matrix determines relative translation only up to scale. Second, the optimization landscape of $E$ is highly nonconvex in the camera rotations and 3D point positions, so bundle adjustment requires a good initialization; this is why SfM is built as an incremental pipeline that grows a partial, already-good reconstruction one camera at a time rather than trying to solve for everything from a cold start. SfM comes in two major flavors. Incremental SfM (the classical and still most common approach, exemplified by Bundler and COLMAP) starts from a single well-conditioned image pair, triangulates an initial point cloud, then repeatedly registers new cameras via Camera Calibration-style pose estimation (PnP) against the existing 3D points, triangulates new points, and periodically re-runs bundle adjustment. Global SfM instead estimates all relative rotations and translations simultaneously via rotation averaging and translation averaging before a single final bundle adjustment; it is faster and more resistant to drift but more fragile to outlier relative poses. This article focuses on the incremental pipeline, since it is the clearer pedagogical path and the historically dominant method.

The incremental SfM pipeline has five stages, executed in this order for each new image, with global refinement interleaved.

Stage 1 — Correspondence and tracks. Detect and describe features in every image (e.g., SIFT), match pairwise using nearest-neighbor plus ratio test, as covered in Feature Descriptors and Matching, then verify each pairwise match geometrically with RANSAC over the fundamental or essential matrix, discarding image pairs with too few inlier matches. Chains of pairwise matches across three or more images are linked into tracks: a track is a single 3D point's set of observations across all images that see it.

Stage 2 — Initialization. Pick the image pair with the largest number of well-distributed inlier matches and a sufficiently wide baseline (narrow-baseline pairs give poorly conditioned triangulation). Estimate the essential matrix (if intrinsics are known) or fundamental matrix, decompose it into a relative rotation and translation (using Epipolar Geometry), and triangulate the shared tracks to seed the initial 3D point cloud and set the world coordinate frame.

Stage 3 — Incremental registration. For each remaining image, find its 2D-3D correspondences against already-triangulated points (via feature matches to already-registered images), and solve the Perspective-n-Point (PnP) problem inside a RANSAC loop to get that camera's pose robustly. This is the same pose-from-known-points problem used in Camera Calibration, just applied per new image instead of to a calibration target.

Stage 4 — Triangulation and extension. With the new camera's pose fixed, triangulate any track that is now observed by two or more registered cameras but has no 3D position yet (using the multi-view triangulation techniques from Stereo Vision and Depth, extended to $\ge 2$ views via linear least squares or DLT).

Stage 5 — Bundle adjustment. Periodically (e.g., every few registered cameras, and always at the end) run bundle adjustment over all currently registered cameras and triangulated points to jointly minimize total reprojection error, correcting the drift that accumulates from doing steps 3–4 one camera at a time. ``` function IncrementalSfM(images): features[i] = DetectAndDescribe(images[i]) for all i pair_matches = { (i,j): MatchAndGeometricVerify(features[i], features[j]) for all pairs i=2 registered images, not yet in points3D for t in new_tracks: points3D[t] = TriangulateTrack(t, poses, registered) FilterOutlierPoints(points3D, poses) # reprojection error / cheirality checks if len(registered) % K == 0: BundleAdjust(poses[registered], points3D) # local, then eventually global BundleAdjust(poses[registered], points3D) # final global refinement return poses, points3D ``` Bundle adjustment itself is a structured nonlinear least-squares solve. Stack all camera parameters $\boldsymbol\theta_i$ (rotation, translation, possibly intrinsics) and point parameters $\mathbf{X}_j$ into a parameter vector $\mathbf{p}$, and let $\mathbf{r}(\mathbf{p})$ be the vector of all reprojection residuals. Levenberg–Marquardt solves

iteratively, where $\mathbf{J}$ is the Jacobian of residuals with respect to $\mathbf{p}$, $\lambda$ is a damping parameter, and $\mathbf{D}$ is typically $\mathrm{diag}(\mathbf{J}^\top\mathbf{J})$. The crucial structural fact is that any residual involving camera $i$ and point $j$ has nonzero derivative only with respect to camera $i$'s parameters and point $j$'s parameters — no residual couples two different cameras directly, or two different points directly. This gives $\mathbf{J}^\top\mathbf{J}$ (and hence the normal equations) an arrowhead / bordered block-diagonal sparsity pattern, exploited by the Schur complement trick: ``` function BundleAdjust(cameras, points, observations): p = pack(cameras, points) lambda = 1e-3 repeat: r, J = ComputeResidualsAndJacobian(p, observations) # sparse: camera & point blocks only

# Partition normal equations by camera block U and point block V (block-diagonal)

# and camera-point coupling block W:

# [ U W ] [dc] [-b_c]

# [ W^T V ] [dp] = [-b_p] U' = U + lambda diag(U) V' = V + lambda diag(V) # V is block-diagonal: one small block per point S = U' - W @ inv(V') @ W.T # Schur complement, size = #camera-params only dc = solve(S, -b_c + W @ inv(V') @ b_p) dp = inv(V') @ (-b_p - W.T @ dc) # back-substitute, cheap: block-diagonal solve p_new = p + pack(dc, dp) if Cost(p_new) < Cost(p): p = p_new; lambda /= 10 else: lambda *= 10 until converged return unpack(p) ```

Feature matching and geometric verification. With $n$ images and up to $F$ features each, brute-force pairwise matching costs $O(n^2 F^2)$ in the worst case; in practice, approximate nearest-neighbor structures (KD-trees, vocabulary trees) reduce per-pair matching to roughly $O(F \log F)$, and vocabulary-tree image retrieval prunes the $O(n^2)$ pair set down to a much sparser candidate set for large unordered collections. Each geometric verification (RANSAC on the fundamental/essential matrix, 8-point or 5-point minimal solver) costs $O(k \cdot F)$ for $k$ RANSAC iterations, and $k$ itself is chosen adaptively from the estimated inlier ratio $w$ via the standard RANSAC sample bound $k = \log(1-p)/\log(1-w^s)$ for confidence $p$ and sample size $s$ (4 or 5 for these models).

Incremental registration. Registering one new camera by PnP-RANSAC costs $O(k' \cdot |\text{2D-3D matches}|)$ where $k'$ is the PnP RANSAC iteration count; triangulating a new track observed by $v$ views costs $O(v)$ per track via linear least squares (solving a small $2v \times 4$ homogeneous system by SVD). Over the whole pipeline, incremental registration is $O(n)$ camera-addition steps, each cheap in isolation — the real cost is in bundle adjustment.

Bundle adjustment: the crux. A naive dense solve of the normal equations for $n$ cameras (say 6 parameters each) and $m$ points (3 parameters each) requires inverting a $(6n+3m)\times(6n+3m)$ matrix, costing $O((6n+3m)^3)$ — for even a modest reconstruction with $m$ in the hundreds of thousands, this is completely infeasible. The Schur complement trick reduces this dramatically. Because each 3D point is observed by only a small number of cameras, the point-point block $\mathbf{V}$ of the normal equations is exactly block-diagonal (one small $3\times 3$ block per point), so $\mathbf{V}^{-1}$ costs only $O(m)$ (invert $m$ independent $3\times 3$ blocks). Eliminating the point variables leaves the reduced camera system $\mathbf{S} = \mathbf{U} - \mathbf{W}\mathbf{V}^{-1}\mathbf{W}^\top$, a $6n \times 6n$ matrix that is itself sparse (camera $i$ and camera $i'$ are coupled only if they share an observed point) and typically solved with sparse Cholesky or a preconditioned conjugate-gradient solver. This brings one Levenberg–Marquardt iteration down to roughly $O(n^3 + m)$ in the worst case for the sparse camera system (often far better in practice with good sparse solvers and problem structure), instead of $O((n+m)^3)$ — the difference between reconstructing thousands of cameras and reconstructing a handful.

Why bundle adjustment is the right objective (correctness argument). If the true 2D measurement noise is independent, zero-mean, isotropic Gaussian, then the reprojection error $E(\mathbf p)$ is (up to constants) exactly the negative log-likelihood of the observed image points given the camera and structure parameters. Minimizing $E$ is therefore the maximum-likelihood estimator for this noise model, which is the formal justification for treating bundle adjustment — rather than, say, minimizing algebraic epipolar error or 3D point distances — as the gold-standard refinement step: it directly optimizes the quantity that was actually measured (pixel locations) rather than a proxy. Two-view methods (essential matrix estimation, linear triangulation) only approximately optimize this likelihood and are used purely as cheap, well-conditioned initializers for the nonconvex BA landscape; without a good initializer, LM iteration converges to an arbitrary local minimum or diverges.

Honest tradeoffs. LM guarantees only convergence to a local minimum (the reprojection error surface is genuinely nonconvex due to rotation parametrization and the perspective projection nonlinearity), so incremental SfM's dependence on a good starting reconstruction is not a minor implementation detail — it is load-bearing. Global SfM avoids incremental drift accumulation by solving rotation and translation averaging problems that are closer to convex, but is more sensitive to outlier relative poses since there is no incremental error-checking step. Both approaches remain fundamentally limited by the similarity ambiguity: absolute scale must come from an external source — a known object size, stereo baseline, GPS/IMU data, or a calibrated pair of cameras with known baseline as in classic Stereo Vision and Depth rigs.

Consider three images of a small textured box, observing four tracked corner points $A, B, C, D$. Assume intrinsics $\mathbf{K}$ are known from prior Camera Calibration, all cameras have unit focal length in normalized coordinates for simplicity, and matches have already been verified by RANSAC.

Step 1 — Tracks after matching.

Step 2 — Initialize from pair (1,2). This pair has the most shared tracks (A, B, C) with a wide-enough baseline. Estimate the essential matrix $\mathbf{E}_{12}$ via the 5-point algorithm inside RANSAC, decompose into four $(\mathbf{R},\mathbf{t})$ candidates, and pick the one where triangulated points A, B, C all lie in front of both cameras (the cheirality check). Suppose this yields:

(translation direction only — scale is set arbitrarily to $\|\mathbf{t}_2\| \approx 1.005$, fixing the global scale of the whole reconstruction). Triangulating A, B, C from images 1 and 2 by linear least squares (solving the $4\times4$ homogeneous DLT system built from the two projection equations per view) gives initial 3D estimates, e.g. $\mathbf{X}_A \approx (0.6, 0.3, 3.1)$.

Step 3 — Register camera 3 via PnP. Camera 3 shares tracks A, B (with known 3D positions from step 2) plus D (not yet triangulated, so unusable for PnP) — only A and B give 2D–3D correspondences, which is below the minimum of 3 needed for a non-degenerate PnP solve. The pipeline therefore defers camera 3 and instead checks whether any other 2D–3D correspondence exists; finding none sufficient, it falls back to also triangulating D using only the (1,2)/(2,3) overlap once camera 2's pose is fixed, or — more realistically in a bigger dataset — simply waits until enough tracks accumulate. This illustrates a genuine pitfall: sparse early overlap can stall registration, which real pipelines handle by trying several candidate next-images and picking the one with the most 2D–3D correspondences, not a fixed order. Suppose in this toy case a fourth image (not shown) actually supplies a third correspondence, so camera 3 registers successfully via PnP-RANSAC, yielding $\mathbf{R}_3, \mathbf{t}_3$ with a reprojection RMS error of 1.3 pixels (equivalent) on inliers — acceptable, so camera 3 is accepted.

Step 4 — Triangulate D. With camera 2 and camera 3 both posed, triangulate D from its two observations $(0.20,0.29)$ and $(0.13,0.28)$, giving $\mathbf{X}_D \approx (0.9, 1.1, 3.4)$.

Step 5 — Bundle adjustment. Before BA, total reprojection error across all 10 observations (A,B,C,D across their respective images) is, say, $E_0 = 4.7\text{ px}^2$ summed squared error, concentrated mostly in point C (only seen from 2 relatively close viewpoints, poorly conditioned). Running LM for a handful of iterations jointly adjusts $\mathbf{R}_2,\mathbf{t}_2,\mathbf{R}_3,\mathbf{t}_3$ and all four 3D points, redistributing error and typically reducing $E$ to well under 1 px² RMS per observation for a well-conditioned scene — the joint adjustment can improve point C's depth estimate specifically because errors that were previously "hidden" in a single pairwise triangulation are now visible as inconsistencies across the whole observation set. This trace makes the two central lessons concrete: geometry from a single wide-baseline pair seeds an arbitrary but internally consistent scale, and every subsequent camera and point addition should be understood as fitting new evidence into that fixed frame, with periodic bundle adjustment as the only step that reconciles all evidence jointly rather than sequentially.

The mathematical roots of Structure from Motion trace to photogrammetry and projective geometry, but its formulation as a computer vision algorithm — and the recognition that bundle adjustment could be scaled to internet-sized unordered photo collections — is a late-20th/early-21st-century achievement built on decades of multi-view geometry work.

Structure from Motion is the geometric backbone underneath a wide range of production systems. Photogrammetry software (Agisoft Metashape, Pix4D, RealityCapture) runs SfM on drone and handheld imagery to produce surveyed 3D models for construction, agriculture, and archaeology. Google Street View, mapping companies, and autonomous-driving perception stacks use SfM-style bundle adjustment to fuse imagery with GPS/IMU data into consistent maps. Visual SLAM systems used in robotics and augmented reality (ARKit, ARCore, ORB-SLAM) run an online, incremental cousin of this same pipeline, adding local bundle adjustment ("pose-graph" or "windowed" BA) to keep tracking drift-free in real time. Most strikingly, the modern wave of neural 3D reconstruction — NeRF, 3D Gaussian Splatting, and their many derivatives — almost universally takes calibrated camera poses produced by classical SfM (typically COLMAP) as their required input, meaning that this decades-old pipeline remains the unglamorous but indispensable first stage of state-of-the-art neural rendering research.

Ask the tutor: How does the choice of rotation parametrization (quaternion vs. axis-angle vs. rotation matrix) affect the convergence and numerical stability of bundle adjustment?