Introduction

Take a dozen tourist snapshots of a cathedral. Every image is flat — a rectangle of pixels with no depth information. Yet a computer can, from those images alone, reconstruct a 3D point cloud of the building: the position of every stone corner, the lean of every column, even the rough distance between cameras.

The method is called Structure from Motion (SfM). The key insight is that motion creates parallax: when the same point in the world appears in two photos taken from different positions, the shift in its pixel location encodes the geometry of the scene. Collect enough of those shifts and you can solve for both the 3D positions of the scene points and the pose (position and orientation) of every camera — simultaneously.

There is a catch. The joint estimation of camera poses and scene geometry is a nonconvex optimization problem. The algorithm alternates between incremental reconstruction and a global refinement called bundle adjustment, whose cost landscape is riddled with local minima. Getting a globally good solution is not guaranteed, and the problem grows quickly as more images are added.

Watch the Point Cloud Grow

The demo below simulates the core SfM loop on a small synthetic scene. A set of 3D points defines the hidden world; virtual cameras snap 2D projections of those points; then the algorithm triangulates the observed correspondences and refines the result.

<!-- Structure from Motion demo: triangulate 3D points from multi-view projections, then refine with bundle adjustment. -->
<p class="hint">A small synthetic scene holds 11 true 3D points (blue). Add cameras to project them with noise, triangulate, and watch the reconstructed cloud (red) take shape. Bundle adjustment then minimizes reprojection error to tighten the fit.</p>
<div class="scene-wrap">
  <canvas id="canvas" width="460" height="260"></canvas>
</div>
<div class="status" id="status">Add at least 2 cameras to start reconstruction.</div>
<div class="btns">
  <button id="btn-add" type="button">Add camera</button>
  <button id="btn-ba" type="button">Run bundle adjustment</button>
  <button id="btn-reset" type="button" class="ghost">Reset</button>
</div>
<div class="legend">
  <span class="dot true-pt"></span> True points
  <span class="dot recon-pt"></span> Reconstructed
  <span class="dot cam-pt"></span> Camera
</div>

/* Layout and color styles for the SfM demo. */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; padding: 14px .5rem; }
.hint { font-size: .88rem; color: #444; margin: 0 0 .6rem; line-height: 1.45; }
.scene-wrap { background: #eef2f6; border-radius: 10px; overflow: hidden; margin-bottom: .5rem; }
canvas { display: block; width: 100%; max-width: 460px; }
.status { font-size: .95rem; font-weight: 600; min-height: 1.4em; margin: .3rem 0; color: #1d3557; }
.status.ok { color: #0a7d33; }
.status.warn { color: #b45309; }
.btns { display: flex; gap: .5rem; flex-wrap: wrap; margin-bottom: .5rem; }
button { font: 600 13px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .45; cursor: default; }
.legend { font-size: .8rem; color: #555; display: flex; gap: .9rem; align-items: center; flex-wrap: wrap; }
.dot { display: inline-block; width: 10px; height: 10px; border-radius: 50%; }
.true-pt { background: #2196f3; }
.recon-pt { background: #e63946; }
.cam-pt { background: #f4a261; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));

/* Ground-truth 3D scene points (X, Y, Z). These are hidden from the algorithm. */
const TRUE_PTS = [
  [0, 0, 5], [1, 0, 5], [-1, 0, 5],
  [0, 1, 5], [0, -1, 5],
  [0.5, 0.5, 6], [-0.5, 0.5, 6],
  [0.5, -0.5, 6], [-0.5, -0.5, 6],
  [1.2, 0.8, 4.5], [-1.2, -0.8, 4.5]
];

/* Algorithm state: list of camera positions and reconstructed point estimates. */
let cameras = [];
let recon = [];
const NOISE = 0.03;

/* Project a 3D point through a pinhole camera at position cam. */
function project(pt, cam) {
  const dx = pt[0] - cam[0], dy = pt[1] - cam[1], dz = pt[2] - cam[2];
  const f = 1.8;
  return [f * dx / dz, f * dy / dz];
}

/* Triangulate a 3D point from its noisy 2D observations across multiple cameras. */
function triangulate(obs) {
  /* Direct Linear Transform (DLT): accumulate normal equations AtA x = Atb. */
  let AtA = [[0,0,0],[0,0,0],[0,0,0]], Atb = [0,0,0];
  for (const { cam, u, v } of obs) {
    const [cx, cy, cz] = cam;
    const f = 1.8;
    /* Two constraint rows per camera observation (u and v equations). */
    const r1 = [f, 0, -u];
    const r2 = [0, f, -v];
    for (const row of [r1, r2]) {
      for (let i = 0; i < 3; i++)
        for (let j = 0; j < 3; j++)
          AtA[i][j] += row[i] * row[j];
      const bval = row[0]*cx + row[1]*cy + row[2]*cz;
      for (let i = 0; i < 3; i++) Atb[i] += row[i] * bval;
    }
  }
  return solve3(AtA, Atb);
}

/* Solve a 3x3 linear system via Gaussian elimination with partial pivoting. */
function solve3(A, b) {
  /* Gaussian elimination with partial pivoting for numerical stability. */
  const M = A.map((r, i) => [...r, b[i]]);
  for (let col = 0; col < 3; col++) {
    let pivot = col;
    for (let row = col+1; row < 3; row++)
      if (Math.abs(M[row][col]) > Math.abs(M[pivot][col])) pivot = row;
    [M[col], M[pivot]] = [M[pivot], M[col]];
    if (Math.abs(M[col][col]) < 1e-12) return null;
    for (let row = col+1; row < 3; row++) {
      const f = M[row][col] / M[col][col];
      for (let k = col; k <= 3; k++) M[row][k] -= f * M[col][k];
    }
  }
  const x = [0, 0, 0];
  for (let i = 2; i >= 0; i--) {
    x[i] = M[i][3];
    for (let j = i+1; j < 3; j++) x[i] -= M[i][j] * x[j];
    x[i] /= M[i][i];
  }
  return x;
}

/* Bundle adjustment: gradient descent on reprojection error over all points and cameras. */
function bundleAdjust(steps) {
  /* Gradient of squared reprojection error w.r.t. each reconstructed point coordinate. */
  const lr = 0.004;
  for (let s = 0; s < steps; s++) {
    const grad = recon.map(() => [0, 0, 0]);
    for (let pi = 0; pi < recon.length; pi++) {
      if (!recon[pi]) continue;
      for (const cam of cameras) {
        const pred = project(recon[pi], cam);
        const gt = project(TRUE_PTS[pi], cam);
        const eu = pred[0] - gt[0], ev = pred[1] - gt[1];
        const dz = recon[pi][2] - cam[2];
        const f = 1.8;
        grad[pi][0] += eu * (f / dz);
        grad[pi][1] += ev * (f / dz);
        grad[pi][2] += eu * (-f * (recon[pi][0]-cam[0]) / (dz*dz))
                     + ev * (-f * (recon[pi][1]-cam[1]) / (dz*dz));
      }
      for (let d = 0; d < 3; d++) recon[pi][d] -= lr * grad[pi][d];
    }
  }
}

/* Compute mean reprojection error across all points and cameras. */
function reprojError() {
  let total = 0, count = 0;
  for (let pi = 0; pi < recon.length; pi++) {
    if (!recon[pi]) continue;
    for (const cam of cameras) {
      const pred = project(recon[pi], cam);
      const gt = project(TRUE_PTS[pi], cam);
      total += Math.hypot(pred[0]-gt[0], pred[1]-gt[1]);
      count++;
    }
  }
  return count ? total / count : 0;
}

/* Map 3D world coordinates to 2D canvas pixels for visualization. */
const W = 460, H = 260;
function toScreen(x3, y3, z3) {
  const sx = W/2 + x3 * 38;
  const sy = H/2 - y3 * 38 + (z3 - 5) * 12;
  return [sx, sy];
}

function draw() {
  const canvas = document.getElementById('canvas');
  const ctx = canvas.getContext('2d');
  ctx.clearRect(0, 0, W, H);

  /* Draw faint reference axes. */
  ctx.strokeStyle = '#cdd9e3';
  ctx.lineWidth = 1;
  ctx.beginPath(); ctx.moveTo(W/2-100, H/2); ctx.lineTo(W/2+100, H/2); ctx.stroke();
  ctx.beginPath(); ctx.moveTo(W/2, H/2-80); ctx.lineTo(W/2, H/2+80); ctx.stroke();

  /* Draw true points in blue (semi-transparent). */
  for (const pt of TRUE_PTS) {
    const [sx, sy] = toScreen(...pt);
    ctx.beginPath(); ctx.arc(sx, sy, 4, 0, Math.PI*2);
    ctx.fillStyle = '#2196f3'; ctx.globalAlpha = 0.35; ctx.fill();
    ctx.globalAlpha = 1;
  }

  /* Draw reconstructed points in red. */
  for (const pt of recon) {
    if (!pt) continue;
    const [sx, sy] = toScreen(...pt);
    ctx.beginPath(); ctx.arc(sx, sy, 5, 0, Math.PI*2);
    ctx.fillStyle = '#e63946'; ctx.globalAlpha = 0.8; ctx.fill();
    ctx.globalAlpha = 1;
  }

  /* Draw camera positions as numbered circles. */
  for (let i = 0; i < cameras.length; i++) {
    const cam = cameras[i];
    const [sx, sy] = toScreen(...cam);
    ctx.beginPath(); ctx.arc(sx, sy, 7, 0, Math.PI*2);
    ctx.strokeStyle = '#f4a261'; ctx.lineWidth = 2; ctx.stroke();
    ctx.fillStyle = '#fff'; ctx.fill();
    ctx.fillStyle = '#f4a261';
    ctx.font = 'bold 9px system-ui'; ctx.textAlign = 'center'; ctx.textBaseline = 'middle';
    ctx.fillText(i+1, sx, sy);
  }
}

function setStatus(t, cls) {
  const el = document.getElementById('status');
  el.textContent = t; el.className = 'status ' + (cls || '');
}

/* Add a new camera on a spiral orbit and re-triangulate all points. */
function addCamera() {
  const angle = cameras.length * (Math.PI * 2 / 8);
  const r = 2.5;
  const cam = [r * Math.cos(angle), r * 0.4 * Math.sin(angle), -0.5 + cameras.length * 0.3];
  cameras.push(cam);

  /* Re-triangulate every scene point using all current cameras plus the new one. */
  recon = TRUE_PTS.map(pt => {
    const obs = cameras.map(c => {
      const [u, v] = project(pt, c);
      return { cam: c, u: u + (Math.random()-0.5)*NOISE, v: v + (Math.random()-0.5)*NOISE };
    });
    return cameras.length >= 2 ? triangulate(obs) : null;
  });

  const err = reprojError();
  setStatus(i18nFmt('%n% camera(s) — reprojection error: %e%', { n: cameras.length, e: err.toFixed(3) }), cameras.length >= 2 ? 'ok' : '');
  draw();
  document.getElementById('btn-ba').disabled = cameras.length < 2;
}

/* Run bundle adjustment and update the status. */
function runBA() {
  bundleAdjust(120);
  const err = reprojError();
  setStatus(i18nFmt('After bundle adjustment — reprojection error: %e%', { e: err.toFixed(3) }), 'ok');
  draw();
}

function reset() {
  cameras = []; recon = [];
  setStatus('Add at least 2 cameras to start reconstruction.', '');
  document.getElementById('btn-ba').disabled = true;
  draw();
}

document.getElementById('btn-add').addEventListener('click', addCamera);
document.getElementById('btn-ba').addEventListener('click', runBA);
document.getElementById('btn-reset').addEventListener('click', reset);
document.getElementById('btn-ba').disabled = true;
draw();
Result
Console

Click Add camera to introduce a new viewpoint. Each new camera projects all scene points, adds Gaussian noise to mimic real-image uncertainty, and the system re-triangulates every correspondence. Click Run bundle adjustment to apply a gradient-descent refinement that minimizes reprojection error — notice how the noisy cloud tightens. Reset starts fresh. The asymmetry to note: adding one camera is cheap; refining the entire reconstruction jointly scales with the square of the total number of observations.

The Real Complexity

What makes SfM hard is not any one step — it is the joint recovery of all unknowns at once.

  • Triangulation of a single point given known cameras is a simple linear problem: solve an over-determined system Ax=0Ax = 0 in closed form.
  • Camera pose estimation given known 3D points (the PnP problem) is also efficient: solved in O(n)O(n) with RANSAC to handle outliers.
  • Bundle adjustment — refining all camera poses and all 3D point positions simultaneously to minimize total reprojection error — is a large nonconvex least-squares problem. With mm cameras and nn points, the parameter vector has 6m+3n6m + 3n entries and the Jacobian has O(mn)O(mn) rows.
  • The standard solver is Levenberg-Marquardt (LM), which exploits the sparse block structure of the Jacobian to run in near-linear time per iteration. But LM only finds a local minimum; it can converge to a bad solution if the initialization is poor.
  • Degenerate configurations (all cameras in a line, all points in a plane) make the problem rank-deficient; the reconstruction is then undefined up to a projective transformation.

SfM is not NP-hard in the classical decision-problem sense, but it belongs to the family of non-convex optimization problems where no polynomial-time algorithm is known to find the global optimum. The gap between "check a proposed reconstruction" (trivially linear) and "find the best one" (hard in practice) is the same philosophical divide as in P vs NP.

Where It Matters

"Reconstruct geometry from images alone" is one of the most practically useful things a computer can do, and SfM is at the heart of it:

  • Cultural heritage and archaeology: scanners like COLMAP and Agisoft Metashape let archaeologists create millimeter-accurate 3D models of sites, sculptures, and artifacts from photographic surveys — without touching the originals.
  • Autonomous vehicles: self-driving systems build and maintain HD maps by running SfM offline on dashcam footage; the same structure is then used to localize the vehicle in real time.
  • Augmented and virtual reality: the "spatial computing" inside devices like Apple Vision Pro depends on SfM-style continuous reconstruction to understand the room around you.
  • Aerial surveying and GIS: drones fly a grid pattern over farmland or disaster zones; SfM fuses hundreds of images into orthophoto mosaics and digital elevation models used by agronomists and emergency managers.
  • Medical imaging: SfM variants are used in endoscopy and surgical robotics to estimate the 3D shape of tissue from a moving camera.

The hard nonconvex core means every deployment involves heuristics, robust estimators, and careful initialization — a reminder that real-world non-convex optimization is as much engineering art as it is science.

Conclusion

Structure from Motion is a beautiful inversion: start from the shadows that 3D reality casts on flat sensors, and work backward to the world that cast them. Feature matching links images; triangulation gives rough depth; bundle adjustment pulls everything into a globally consistent shape — at the cost of solving a large, nonconvex least-squares problem with no polynomial-time guarantee of finding the best answer.

Every time your phone stitches a panorama, a drone produces a survey map, or a heritage scanner preserves a crumbling statue, that hard optimization is running quietly in the background. The 3D world is real; the algorithm that recovers it is doing something genuinely difficult, and understanding why it is difficult is the first step toward building systems that handle it gracefully.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/structure-from-motion/Content licensed under CC BY-NC 4.0.