Multi-View Stereo

Multi-view stereo (MVS) computes dense 3D geometry of a scene from many photographs whose cameras are already calibrated and posed. It is the step that follows structure from motion: SfM recovers the cameras and a sparse cloud of feature points, and MVS uses those cameras to estimate a surface point for nearly every pixel. Two-view stereo matching is its special case; with more views, each surface point is seen from several baselines, and a point hidden in one image can still be matched in others.

Problem Definition

Seitz et al. (2006) state the goal as reconstructing a complete 3D model from images taken from known viewpoints. The input is a set of images $I_0, \dots, I_N$ with camera matrices $P_i = K_i[R_i \mid \mathbf{t}_i]$, usually from SfM and bundle adjustment. The output is a surface, or samples of it, that is photo-consistent: each surface point should look the same in every image in which it is visible, up to noise and lighting.

Because the cameras are known, each pixel’s correspondence search runs along its viewing ray, parameterized by depth; with many cameras in arbitrary positions, images are warped onto hypothesized surfaces instead of being rectified. The result inherits the cameras’ scale, which images alone fix only up to a global factor (multi-view geometry). Some prior on the scene’s extent is also needed: Seitz et al. note that otherwise a different flat picture placed in front of each camera would explain the images. Algorithms take a bounding volume, silhouettes, a depth range per view, or matched feature points.

Input and Output

• Input: images with intrinsics and poses, often with the sparse SfM point cloud, which supplies depth ranges and shows which images overlap.
• Output: a depth map per view, a dense point cloud or set of oriented patches, a triangle mesh, or a volume (occupancy or signed distance on a 3D grid). Pipelines often pass through several, for example depth maps fused into points and then meshed.

Variants

• Objects versus scenes. Single objects seen from all around, scenes that extend beyond the images and contain clutter, and crowded scenes with moving obstacles (Furukawa and Ponce, 2010).
• Controlled versus unstructured capture. Internet photo collections, the target of Schönberger et al. (2016), vary in resolution, illumination, and viewpoint, and some views are misregistered.

Challenges

• Occlusion. A point is visible in only some views, and which ones depends on the unknown geometry. Seitz et al. (2006) group visibility models into geometric reasoning from the current shape, heuristics such as using only nearby cameras, and treating occluded views as outliers.
• Weak texture. Aanæs et al. (2016) found lack of texture, with meshing, the most pressing issue on DTU; specularity and lighting changes caused fewer problems than they expected.
• Non-Lambertian surfaces look different from each viewpoint; the Middlebury evaluation used nearly Lambertian objects because most algorithms assume them.
• Scale. Volumes and cost volumes grow with the cube of their resolution.

Approaches

Photo-consistency

Methods score a hypothesized surface point by projecting it, or a small patch around it, into the images. Measures include the variance of the projected colors and window-based sum of squared differences or normalized cross-correlation (NCC) between image pairs (Seitz et al., 2006), defined under cost volume; NCC ignores gain and offset changes in brightness. Per-view costs such as $1 - \mathrm{NCC}$ are combined over the views in which the point is believed visible. A plain average lets one occluded view spoil a correct hypothesis, so robust variants average only the best views, drop views that score poorly, or select views per pixel.

Taxonomy

Seitz et al. (2006) classified MVS algorithms by six properties, among them the reconstruction algorithm, of which they described four classes:

1. Volumetric cost, then surface extraction, as in voxel coloring, which sweeps once through a grid keeping low-cost voxels, or graph cuts on a volumetric random field.
2. Surface evolution by removing voxels, level sets, or mesh deformation. Space carving (Kutulakos and Seitz, 2000) belongs here: it removes voxels that are not photo-consistent from an enclosing volume. Kutulakos and Seitz proved that, within a given initial volume, the union of all photo-consistent shapes, the photo hull, contains every other one, and that carving away only inconsistent voxels converges to it. Because carving keeps what it cannot reject, the result bulges outward where texture is weak (Seitz et al., 2006).
3. Depth-map methods, which compute a depth map per view and enforce consistency between them or merge them.
4. Feature-based methods, which match and triangulate features and fit a surface to them.

Plane-sweep stereo

Collins (1996) swept a plane through space and counted, in each cell of the plane, the back-projected feature rays from all images passing through it; cells with enough votes became 3D features. Each image maps to the plane by a homography, and the cost grows linearly with the number of images, one of Collins’s requirements for true multi-image matching.

Dense plane sweeping fixes a reference camera $K[I \mid \mathbf{0}]$ and sweeps fronto-parallel planes $\mathbf{n}^\top\mathbf{X} = z$, $\mathbf{n} = (0, 0, 1)^\top$. A source camera $K_k[R_k \mid \mathbf{t}_k]$ sees the plane at depth $z$ through

$$H_k(z) = K_k\left(R_k + \frac{\mathbf{t}_k\,\mathbf{n}^\top}{z}\right)K^{-1},$$

which maps reference pixels to source pixels (the homography article’s plane $\mathbf{n}^\top\mathbf{X} + d = 0$ with $d = -z$, which turns its minus sign into a plus). Warping each source image by $H_k(z)$ and comparing it with the reference gives one depth slice of the cost volume. Image shifts are proportional to inverse depth, as with disparity, so planes are spaced uniformly in $1/z$. The volume is read out by winner-take-all or regularized, for example with semi-global matching. Fronto-parallel planes favor fronto-parallel surfaces; several plane orientations or per-pixel normals reduce the bias.

Patch-based MVS

PMVS (Furukawa and Ponce, 2010) represents the surface by small oriented rectangles, each with a center, a normal, and a reference image. A patch’s discrepancy is the average of one minus NCC between its reference image and the other images that see it, sampled on a 5 × 5 or 7 × 7 grid on the patch; images that disagree strongly with the reference are dropped, which handles specular highlights and obstacles. The algorithm is match, expand, filter: matched Harris and difference-of-Gaussian features seed patches, patches spread into neighboring image cells, and filtering removes patches that violate visibility constraints. It needs no bounding box or depth range, and its oriented points can be meshed by Poisson surface reconstruction (point cloud).

PatchMatch MVS

PatchMatch Stereo (Bleyer, Rhemann and Rother, 2011) propagates random per-pixel slanted planes to neighbors (stereo matching). In MVS, each reference pixel holds a depth and a normal whose plane induces a homography into every source image. Schönberger et al. (2016), building on a probabilistic framework by Zheng et al. and released in COLMAP, estimate per pixel the depth, the normal, and which source images to trust, guided by priors on triangulation angle, incidence angle, and resolution. A geometric consistency term checks each depth against those estimated in the source views, and fusion merges pixels whose depths, normals, and reprojection errors agree into an oriented point cloud.

Depth-map fusion and meshing

Curless and Levoy (1996) fused range images volumetrically. Each depth map becomes a signed distance to the surface along the lines of sight, truncated to a band around the surface and added into a cumulative weighted signed distance function on a voxel grid; the surface is its zero-level isosurface, which they argued is least-squares optimal under their assumptions. This truncated signed distance function (TSDF) is widely used to fuse depth maps from stereo or depth cameras.

Learned MVS

MVSNet (Yao et al., 2018) makes plane sweeping learnable: it warps learned features by differentiable homographies, combines any number of views by their variance, and regularizes the volume with 3D convolutions; depth is the expectation over hypotheses (soft argmin), refined with the reference image. Because the 3D volume dominates memory, R-MVSNet (Yao et al., 2019) regularizes 2D cost maps sequentially along depth with a gated recurrent unit, and the cascade cost volume of Gu et al. (2020) narrows the depth range at each finer level of a feature pyramid. These networks still need known cameras, and their per-view depth maps are still filtered and fused.

Relation to neural rendering

Neural radiance fields and 3D Gaussian splatting also fit a scene to calibrated images, usually with SfM cameras, but aim at novel-view synthesis rather than an explicit surface.

Example

The script ray-casts a reference view and four source views, 0.4 m to its left, right, top, and bottom, of a textured slanted plane ($Z = 4 + 0.3X$ m) behind a textured box. It sweeps 160 fronto-parallel planes from 2 m to 6 m, uniform in inverse depth, warps each source image with $H_k(z)$, computes a $7 \times 7$ NCC cost, aggregates it over views, and picks depth by winner-take-all. Errors are measured away from the image border, and separately on pixels whose surface point the box hides in at least one source view.

import warnings
import numpy as np
from scipy import ndimage

rng = np.random.default_rng(0)
H, W, f = 120, 160, 160.0
K = np.array([[f, 0, W / 2], [0, f, H / 2], [0, 0, 1]])
Kinv = np.linalg.inv(K)

# Scene in reference-camera coordinates: a slanted plane Z = 4 + 0.3 X behind
# a box spanning X in [-0.5, 0.3], Y in [-0.35, 0.35], Z in [2.5, 3.0].
box_lo, box_hi = np.array([-0.5, -0.35, 2.5]), np.array([0.3, 0.35, 3.0])
freq = rng.normal(0, 12, (40, 3))           # random 3D sinusoids = texture
phase = rng.uniform(0, 2 * np.pi, 40)

def texture(P):
    return np.sin(P @ freq.T + phase).mean(axis=-1)

def look_at(C, target=np.array([0.0, 0.0, 3.5])):
    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.stack([x, np.cross(z, x), z])     # rows: camera axes in world frame
    return R, -R @ C                         # X_cam = R X + t

def render(R, t):
    """Ray-cast the scene; return the image and the depth along the optical axis."""
    u, v = np.meshgrid(np.arange(W), np.arange(H))
    rays = np.stack([u, v, np.ones_like(u)], -1).reshape(-1, 3) @ Kinv.T
    d = rays @ R                              # ray directions in the world (R^T r)
    C = -R.T @ t
    s_plane = (4 + 0.3 * C[0] - C[2]) / (d[:, 2] - 0.3 * d[:, 0])
    with np.errstate(divide="ignore"):        # slab test for the box
        t1, t2 = (box_lo - C) / d, (box_hi - C) / d
    near = np.minimum(t1, t2).max(1); far = np.maximum(t1, t2).min(1)
    s_box = np.where((near < far) & (near > 0), near, np.inf)
    s = np.minimum(s_plane, s_box)
    P = C + s[:, None] * d
    return texture(P).reshape(H, W), (s * rays[:, 2]).reshape(H, W)

ref, gt = render(np.eye(3), np.zeros(3))
sources = [look_at(np.array(c)) for c in
           [(0.4, 0, 0), (-0.4, 0, 0), (0, 0.4, 0), (0, -0.4, 0)]]
src_imgs = [render(R, t)[0] for R, t in sources]

depths = 1 / np.linspace(1 / 2.0, 1 / 6.0, 160)   # uniform in inverse depth
u, v = np.meshgrid(np.arange(W), np.arange(H))
pix = np.stack([u.ravel(), v.ravel(), np.ones(H * W)])
n = np.array([0.0, 0.0, 1.0])                     # fronto-parallel planes Z = z
win = lambda a: ndimage.uniform_filter(a, 7)       # 7 x 7 window

def ncc_cost(a, b):
    ma, mb = win(a), win(b)
    cov = win(a * b) - ma * mb
    va, vb = win(a * a) - ma ** 2, win(b * b) - mb ** 2
    return 1 - cov / np.sqrt(va * vb + 1e-6)       # 0 = perfect match

# cost[k, j] = matching cost of source view k at depth plane j
cost = np.full((len(sources), len(depths), H, W), np.nan)
for k, ((R, t), img) in enumerate(zip(sources, src_imgs)):
    for j, z in enumerate(depths):
        Hz = K @ (R + np.outer(t, n) / z) @ Kinv   # reference pixel -> source pixel
        q = Hz @ pix
        x, y = q[0] / q[2], q[1] / q[2]
        warped = ndimage.map_coordinates(img, [y, x], order=1, cval=np.nan).reshape(H, W)
        valid = win(np.isfinite(warped).astype(float)) > 0.999
        c = ncc_cost(ref, np.nan_to_num(warped))
        cost[k, j] = np.where(valid, c, np.nan)

# Ground-truth visibility: is the reference pixel's surface point occluded in source k?
X = (pix * gt.ravel()).T @ Kinv.T                   # reference pixels back-projected
occluded = np.zeros((len(sources), H, W), bool)
for k, (R, t) in enumerate(sources):
    Xc = X @ R.T + t
    q = Xc @ K.T
    zs = ndimage.map_coordinates(render(R, t)[1], [q[:, 1] / q[:, 2], q[:, 0] / q[:, 2]],
                                 order=0, cval=np.inf)
    occluded[k] = (Xc[:, 2] > zs + 0.05).reshape(H, W)
inner = np.zeros((H, W), bool); inner[12:-12, 12:-12] = True
occ_any = occluded.any(0) & inner
print(f"evaluated pixels: {inner.sum()}, occluded in at least one view: {occ_any.sum()}")

def errors(agg):
    est = depths[np.argmin(np.where(np.isnan(agg), np.inf, agg), axis=0)]
    return np.abs(est - gt)

def report(name, agg):
    e = errors(agg)
    print(f"{name:<22s} median {np.median(e[inner]) * 100:4.1f} cm | "
          f"within 5 cm: all {np.mean(e[inner] < 0.05) * 100:5.1f}%, "
          f"occluded-in-some-view {np.mean(e[occ_any] < 0.05) * 100:5.1f}%")

with warnings.catch_warnings():   # all-NaN slices: no valid view
    warnings.simplefilter("ignore", RuntimeWarning)
    for m in range(1, len(sources) + 1):
        report(f"{m} view(s), mean", np.nanmean(cost[:m], axis=0))
    # Robust aggregation: drop the worst view at each hypothesis, so a view in
    # which the pixel is occluded does not spoil the correct depth.
    report("4 views, best 3 of 4", np.nanmean(np.sort(cost, axis=0)[:3], axis=0))
    report("4 views, best 2 of 4", np.nanmean(np.sort(cost, axis=0)[:2], axis=0))

Output:

evaluated pixels: 13056, occluded in at least one view: 1722
1 view(s), mean        median  1.7 cm | within 5 cm: all  75.9%, occluded-in-some-view  52.7%
2 view(s), mean        median  1.0 cm | within 5 cm: all  87.1%, occluded-in-some-view  47.6%
3 view(s), mean        median  1.2 cm | within 5 cm: all  84.4%, occluded-in-some-view  39.3%
4 view(s), mean        median  1.0 cm | within 5 cm: all  86.1%, occluded-in-some-view  36.2%
4 views, best 3 of 4   median  1.3 cm | within 5 cm: all  88.1%, occluded-in-some-view  74.4%
4 views, best 2 of 4   median  1.7 cm | within 5 cm: all  80.4%, occluded-in-some-view  69.7%

A second view raises the share of pixels within 5 cm from 76% to 87%: a wrong match that fits one image pair rarely fits another. The median error is set mainly by the plane spacing, about 1.3 cm at the box and 2.6–4.4 cm on the evaluated background. Further views barely help a plain mean, and they hurt the 13% of pixels that some view cannot see, from 53% to 36% correct, because each added view is another chance to average an occluded cost into the true depth. Dropping the worst view at each hypothesis lifts these pixels to 74% and gives the best overall result; keeping only two of four discards too much evidence.

Datasets and Benchmarks

• Middlebury multi-view stereo (Seitz et al., 2006): two objects, temple and dino, imaged at 640 × 480 from a spherical gantry in hemisphere, ring, and sparse-ring sets, with laser-scanned ground truth. It is separate from the two-view Middlebury benchmarks.
• DTU (Aanæs et al., 2016): 80 scenes seen from 49 or 64 robot-positioned cameras under seven lighting conditions, with structured-light reference scans. MVSNet was trained and evaluated on it.
• Tanks and Temples (Knapitsch et al., 2017): large indoor and outdoor scenes captured as high-resolution video, with industrial laser-scanner ground truth, in an intermediate and an advanced group; it evaluates complete SfM and MVS pipelines.
• ETH3D (Schöps et al., 2017): high-resolution DSLR images and multi-camera videos of indoor and outdoor scenes, registered to laser scans.

Evaluation Metrics

Reconstructions $\mathcal{R}$ are compared with ground truth $\mathcal{G}$ by nearest-neighbor distances in both directions.

• Accuracy uses the distances $e_{r \to \mathcal{G}} = \min_{g \in \mathcal{G}} \lVert r - g \rVert$ from reconstructed points. Middlebury reports the distance within which 90% of them lie; DTU reports the mean and median after discarding distances above 20 mm.
• Completeness uses the reverse distances $e_{g \to \mathcal{R}}$. Middlebury reports the percentage of ground-truth points within 1.25 mm of the reconstruction; DTU again reports mean and median.
• Precision, recall, and F-score. Tanks and Temples defines precision $P(d)$ as the percentage of reconstructed points within $d$ of the ground truth, recall $R(d)$ as the percentage of ground-truth points within $d$ of the reconstruction, and their harmonic mean

$$F(d) = \frac{2\,P(d)\,R(d)}{P(d) + R(d)},$$

with a per-scene threshold $\tau$ as the main measure. A few precise points give low recall and filling space gives low precision, so neither raises $F$. ETH3D calls the same quantities accuracy and completeness, reports them at thresholds from 1 cm to 50 cm, ranks by F1, averages over small voxels so that point density cannot game the scores, and ignores points in space its scanner never observed.

Applications

• 3D models from photographs for film, games, and metrology (Furukawa and Ponce, 2010), and aerial mapping, the setting of Collins’s (1996) example.
• Landmarks and cities from internet photos, and dense models for image-based rendering, localization, and classification (Schönberger et al., 2016).

Open Problems

• Textureless and non-Lambertian surfaces, such as white walls, glass, and mirrors, where photo-consistency is uninformative or wrong.
• Meshing noisy, incomplete point sets into surfaces without smoothing away detail.
• Generalization of learned methods to scenes and cameras unlike their training data.