A Volumetric Method for Building Complex Models from Range Images

“A Volumetric Method for Building Complex Models from Range Images” by Brian Curless and Marc Levoy of Stanford University appeared at SIGGRAPH ‘96 in August 1996. It merged many aligned range scans into a single surface by averaging signed distances along the sensor’s lines of sight in a voxel grid, extracting the zero level set, and plugging holes with surfaces placed where space observed to be empty meets space never observed. The authors reported models of up to 2.6 million triangles from as many as 70 scans, including the dragon and the Happy Buddha later released in the Stanford 3D Scanning Repository. The averaging scheme, now called TSDF fusion, became the mapping step of KinectFusion fifteen years later. Section, equation and figure numbers below are from the authors’ PDF on the Stanford paper page, paginated 1 to 10.

Problem

By the mid-1990s, fast range scanners produced range images, depth values on a regular lattice, for medicine, reverse engineering and film-making (§1). The paper’s scanner was a laser-stripe optical triangulation system: a laser sheet sweeps the object and an off-axis CCD camera turns the reflected stripe into depth profiles (Fig. 1b). The error of such a measurement is ellipsoidal about each range point, with its long axis along a line of sight (§3, citing earlier analyses): the laser’s, or the camera’s with the authors’ spacetime analysis (§5.1). One range image never covers a whole object, so many must be merged.

The paper’s review (§2) sorts earlier methods as follows.

• Unorganized points. Delaunay triangulation, alpha shapes and Hoppe et al.‘s signed distance, in the authors’ assessment, discard normals and reliability estimates and are fragile at high curvature and with outliers or systematic distortions.
• Polygonal merging. Turk and Levoy’s zippering (1994) eroded redundant triangles, stitched the remaining boundaries and moved vertices to a consensus of the original data; Soucy and Laurendeau and Rutishauser et al. merged meshes otherwise. All could still fail catastrophically at high curvature.
• Volumes. Octree and binary-voxel methods had not produced surfaces and cannot represent uncertainty or average noise; stereo-based averaging gave sparse, low-resolution results. Hilton et al. (ECCV, April 1996) also merged weighted signed distance functions but, in the authors’ view, did not address directional uncertainty, incremental updating, space efficiency or hole filling.
• Occupancy grids. Elfes and Matthies’s probabilistic grids put the surface on a ridge of the probability field, and ridge finding lacked robust solutions, a main reason the authors extract an isosurface instead.

Contribution

The introduction lists seven properties a merging algorithm should have, and the paper claims all of them, where earlier algorithms had subsets (§1, abstract):

1. Representation of range uncertainty, which is directional for range scanners.
2. Use of all range data, so that redundant observations average out noise.
3. Incremental and order-independent updating: a model after any scan, no bias from earlier scans, easy parallelization.
4. Time and space efficiency for the many scans a complex object needs.
5. Robustness to outliers and systematic distortions, without holes or self-intersecting surfaces.
6. No restriction on topological type.
7. Ability to fill holes with plausible surfaces, giving a watertight model.

Method

Cumulative signed distance. Each range image is triangulated without bridging depth discontinuities, its vertices weighted as Soucy and Turk did, by the dot product of normal and viewing direction, less near mesh boundaries. A ray from the sensor through each voxel near the surface meets this mesh, giving a signed distance and an interpolated weight (Fig. 5) for a running weighted average (Eqs. 1–4). Values are stored in fixed point between $D_{\min} < 0$ and $D_{\max} > 0$ (§3). The TSDF fusion article covers the update, the 1D weighted mean (Eq. 5), the falloff behind the surface and Appendix A’s least-squares sketch (Eq. 9), whose full proof the paper defers to Curless’s 1996 PhD thesis.

Sign convention. The paper’s ramps increase away from the sensor (Fig. 3): empty space holds $D_{\min}$ and unseen space $D_{\max}$ (§4, Fig. 6b). This site and KinectFusion put positive values in front of the surface, and KinectFusion (§2.4) and Nießner et al. (§3) describe the Curless and Levoy representation with that sign, the reverse of the paper’s figures. The zero set is unaffected.

Hole filling. Voxels are unseen, empty (carved along lines of sight in front of a range surface) or near the surface (§4, Fig. 6); the TSDF fusion article explains how one extraction across zero-weight voxels then yields hole fillers on the empty–unseen frontier, joined to the observed surface. The binary transition aliases, so only the hole-fill triangles are smoothed, by weighted averaging of neighbouring vertices. Scans of a backdrop placed outside the grid carve more space without adding geometry.

Implementation. The volume is run-length encoded with empty, unseen and varying runs, saving 10:1 to 20:1, usually less memory than the final mesh (§5.2.1). Likening fusion to inverse volume rendering, the authors resample each range image so that its scanlines align with the voxel grid and stream through both into a double-buffered volume (§5.2.2, Fig. 7); a binary tree of depth bounds skips voxels far from the surface, a typical speed-up of 15 without carving and 5 with it. The surface is extracted with marching cubes (Lorensen and Cline, 1987), using Montani et al.’s lookup table that resolves the ambiguous cases, only in varying voxels and at empty–unseen boundaries (§5.2.3).

Example

The code checks the paper’s 1D argument (Fig. 3, Eq. 5) along one line of sight: 12 scans of a surface at 100 mm, noisier at grazing incidence and weighted by the cosine, fused as truncated ramps on 0.35 mm voxels (the dragon’s voxel size). The ramps use this site’s sign.

import numpy as np

# One line of sight through a voxel column: the surface is at range 100 mm.
# Each of n scans measures it with noise sigma0 / cos(theta) (grazing views are noisier)
# and gets the paper's weight, the cosine of the angle between normal and view (Sec. 3).
# Ramps use this site's sign (positive in front of the surface); the paper's is opposite.
TRUE, H, SIGMA0, N = 100.0, 0.35, 0.10, 12          # mm; 0.35 mm voxels as for the dragon
x = np.arange(98.0, 102.0, H)                        # voxel centres along the ray


def fuse(r, w, delta, order=None):
    """Eqs. 3-4: running weighted average of ramps d_i = r_i - x, kept where |d_i| <= delta."""
    D, W = np.zeros_like(x), np.zeros_like(x)
    for i in (range(len(r)) if order is None else order):
        d = r[i] - x
        wi = np.where(np.abs(d) <= delta, w[i], 0.0)   # ramp and weight cut off at +-delta
        D = np.where(wi > 0, (W * D + wi * d) / np.maximum(W + wi, 1e-30), D)
        W = W + wi
    return D, W


def zero_crossing(D, W):
    """First + to - sign change between observed voxels, linearly interpolated."""
    for j in range(len(x) - 1):
        if W[j] > 0 and W[j + 1] > 0 and D[j] > 0 >= D[j + 1]:
            return x[j] + H * D[j] / (D[j] - D[j + 1])
    return np.nan


def scans(seed):
    rng = np.random.default_rng(seed)
    cos = rng.uniform(0.2, 1.0, N)
    return TRUE + rng.normal(0, SIGMA0 / cos), cos


print("ramp half-width  max |R - weighted mean| (mm), seeds 0-99  failures")
for delta in (0.1, 0.2, 0.5, 1.0, 2.0):
    dev, fail = [], 0
    for s in range(100):
        r, w = scans(s)
        R = zero_crossing(*fuse(r, w, delta))
        fail += np.isnan(R)
        if not np.isnan(R):
            dev.append(abs(R - (w @ r) / w.sum()))      # Eq. 5
    print(f"{delta:10.1f} mm  {max(dev):30.4f}  {fail:8d}")

r, w = scans(0)
D1, _ = fuse(r, w, 2.0)
D2, _ = fuse(r, w, 2.0, order=np.random.default_rng(1).permutation(N))
print(f"\nshuffled-order check, seed 0: max |D - D_shuffled| = {np.abs(D1 - D2).max():.1e} mm")

err = np.array([[r[0] - TRUE, r.mean() - TRUE, zero_crossing(*fuse(r, w, 2.0)) - TRUE]
                for r, w in map(scans, range(1000))])
print("\nRMS error over seeds 0-999 (mm): first scan alone, unweighted mean, fused zero crossing")
print("  ".join(f"{e:.4f}" for e in np.sqrt((err ** 2).mean(0))))

Output:

ramp half-width  max |R - weighted mean| (mm), seeds 0-99  failures
       0.1 mm                          0.5042        38
       0.2 mm                          0.4246         3
       0.5 mm                          0.0507         0
       1.0 mm                          0.0459         0
       2.0 mm                          0.0000         0

shuffled-order check, seed 0: max |D - D_shuffled| = 4.4e-16 mm

RMS error over seeds 0-999 (mm): first scan alone, unweighted mean, fused zero crossing
0.2312  0.0650  0.0489

With ramps 2 mm long on each side, every scan contributes at the crossing, the fused field is linear there, and the zero crossing equals the weighted mean of Eq. 5 for all 100 seeds. With 0.5 or 1.0 mm, a few grazing scans miss the crossing and shift it by up to 0.05 mm; with 0.1 mm, 38 of 100 seeds give no crossing at all: the trade-off of §3. Shuffling the scans changes only rounding. Over 1000 seeds the fused surface has an RMS error of 0.049 mm, against 0.231 mm for one scan and 0.065 mm for an unweighted mean; the weights help because this simulation gives the noisiest views the lowest weights.

Results

All scans came from a Cyberware 3030 MS scanner, using either its hardware triangulation or the spacetime analysis, with distances laid along the matching lines of sight (§5.1). Inputs reached 70 scans, 12 million vertices and grids of up to 160 million voxels (§5.2). Figure 8 gives the statistics:

Model	Scans	Input triangles	Voxel size (mm)	Volume dimensions	Time (min)	Output triangles	Holes
Dragon	61	15 M	0.35	712 × 501 × 322	56	1.7 M	324
Dragon + fill	71	24 M	0.35	712 × 501 × 322	257	1.8 M	0
Buddha	48	5 M	0.25	407 × 957 × 407	47	2.4 M	670
Buddha + fill	58	9 M	0.25	407 × 957 × 407	197	2.6 M	0

• Timing and accuracy. Merging the observed surfaces took under an hour on a 250 MHz MIPS R4400 processor; carving and hole filling, not fully optimized, took 3–5 hours, less than scanning and alignment. The RMS distance between range points and the reconstruction was about 0.1 mm for both models, roughly the scanner’s accuracy (§6).
• Dragon (Fig. 11). Hole filling on the 61 scans closed some holes plausibly but left uncarved regions that produced extraneous surfaces. Ten backdrop scans carved that space, and after smoothing the fillers the 1.8 million polygon model was watertight.
• Happy Buddha (Fig. 12). The 20 cm plastic and rosewood statuette, painted matte grey, merged from 48 scans into a mesh less noisy than a single scan, with a hole in its base that ten backdrop scans filled. A stereolithography hardcopy was built in 500 layers 150 µm apart.
• Drill bit (Fig. 9). A 1.6 mm drill bit scanned from 12 orientations with conventional triangulation, whose false edge extensions confuse merging, gave a watertight volumetric model where zippering failed catastrophically.

Released models. As of October 2026, the Stanford 3D Scanning Repository (page last modified 28 May 2026) lists the Happy Buddha, the dragon and the drill bit as first appearing in this paper, with raw range data (backdrop scans included for the Buddha and dragon) and reconstructions made with the method, whose implementation Stanford distributes as VripPack. Its figures differ from the paper’s: about 60 Buddha and about 70 dragon scans, and conservatively decimated reconstructions of about 1.09 and 1.13 million triangles; it describes the released dragon as containing numerous small holes. The page warns that the merged models hide outliers, noise and misalignment, making them unrealistic input for evaluating surface reconstruction.

Impact

By the repository’s account, the Digital Michelangelo Project’s Scanalyze software invoked this volumetric merging, which also reconstructed the Lucy statue at 0.5 mm with hole filling. KinectFusion (Newcombe et al., 2011) fused every frame of a handheld Kinect on a GPU with the paper’s weighted running average (its Eqs. 11–12), citing it for the representation and the truncation (§2.4, §3.3); it also records that Rusinkiewicz et al.’s 2002 real-time scanning system optimized its final models offline with this method (§2.5). Nießner et al. (2013) present voxel hashing explicitly as an extension of the Curless and Levoy method, keeping its cumulative average and adding sparse storage so that it scales (§3). DeepSDF (Park et al., 2019) credits the paper, with KinectFusion, with pioneering TSDF fusion of noisy depth maps (§2).

Limitations

The paper’s discussion (§7) names algorithmic limitations. Sharp corners are hard to bridge when no scan sees both faces. Behind each surface the ramp interferes with scans of the opposite side, so thin surfaces thicken and sharp corners round off. Other limits come from optical triangulation: only external surfaces are seen, complex objects need very many scans, laser and camera must both see each point, and shiny, dark or bright surfaces cause trouble, so objects were often painted flat grey.

Hole fillers are plausible bounds, not measurements, and need enough carving to avoid false surfaces. The repository adds that the carved Buddha contains small bridges, so its genus is larger than it looks. Kazhdan et al. (2006, §5.2) found “lipping” at sharp creases in this method’s output, which they called VRIP, on the Buddha and a fragment of the Forma Urbis Romae, attributing it to the distance function being grown, as they put it, perpendicular to the view direction rather than the surface normal. The method also assumes aligned scans, treating the remaining error as the scanner’s (§3).

What Came After

The averaging carried over; the sensor, speed and storage changed. KinectFusion projected voxels into each depth map instead of casting rays into a triangulated range image, and tracked the camera with iterative closest point against the fused model. Voxel hashing removed the dense grid’s memory cost, the role run-length encoding played in 1996 (our comparison). Poisson surface reconstruction offered an alternative for oriented points without lines of sight, and DeepSDF moved signed distance functions from grids into networks (neural implicit representations).