Marching Cubes: A High Resolution 3D Surface Construction Algorithm

“Marching Cubes: A High Resolution 3D Surface Construction Algorithm” is a seven-page paper by William E. Lorensen and Harvey E. Cline of General Electric Corporate Research and Development, presented at SIGGRAPH ’87 in Anaheim and published in Computer Graphics, volume 21, number 4, in July 1987. It turned a stack of medical slices into a triangle model of one constant-density surface, one cube of eight neighbouring samples at a time. The method as used today, with its ambiguities and their repairs, is covered in marching cubes; this article is about the 1987 paper.

Problem

The motivation was clinical. Computed tomography (CT), magnetic resonance (MR) and single-photon emission computed tomography (SPECT) produce stacks of 2D slices that radiologists are trained to read but referring physicians struggle to picture in 3D. The authors split the pipeline into acquisition, image processing, surface construction (the paper’s subject, with the surface chosen by a density value) and display (Figure 1).

The related work names four families of methods and what each throws away. Connecting contours on consecutive slices, as in Keppel’s triangulation, is ambiguous when a slice holds several contours, and user interaction, which partly helps, should be minimal in a clinic; it also discards the connectivity between slices. Herman and colleagues’ cuberille approach thresholds the volume into blocks and tries to recover shading from them. Ray casting by Farrell and by Höhne shades by depth alone or with an unnormalized gradient. The Mayo Clinic’s volume display shows every value and relies on motion for the 3D impression. The aim was to take connectivity, surface position and gradient from the original data and output triangles for conventional renderers.

Contribution

• A case table built by divide and conquer. The surface is located in one logical cube at a time, formed by four pixels on each of two adjacent slices (Figure 2). A table indexed by the 8-bit corner pattern (Figure 4), $2^8 = 256$ cases, lists the edges the surface crosses and their triangles.
• Two symmetries for building the table. Swapping inside and outside leaves the topology of the triangulated surface unchanged, and rotations of the cube map cases onto each other. The authors triangulated a few basic patterns by hand, each non-empty one with one to four triangles (Figure 3), and generated the 256 entries from them.
• Linear interpolation along edges. Higher-degree interpolation, also tried, helped little with at most four triangles per cube.
• Normals from the data gradient. A constant-density surface has no gradient component along its tangent, so the normalized gradient is the normal. It is estimated at cube corners by central differences (Eqs. 2–4) and interpolated to vertices for Gouraud shading.
• Efficiency and solid modelling. Edge reuse between neighbouring cubes, and Boolean operations built on the cube index.

Method

The seven-step summary in §4 keeps four slices in memory: two for the cubes and two more for the central differences.

Coherence. For a cube not on the first slice, line or pixel, only three of its twelve edges are new (edges 6, 7 and 12 in Figure 5). Saving only the previous pixel’s and line’s intersections, since a slice’s needed too much memory, gave a threefold speed-up. Averaging four pixels into one cut the triangle count for high-resolution slices, at some loss of detail.

Solid modelling. The index classifies a cube as outside (0), inside (255) or on the surface, and any surface with an inside–outside function, such as a plane $ax + by + cz - d$ or a connectivity mask, classifies cubes the same way. A truth table over the two classifications gives Boolean intersection (Figures 6 and 7): a cube on one surface and inside the other keeps that surface’s triangles, a cube outside either gives none, and a cube on both has its triangles clipped against each other with the Sutherland–Hodgman algorithm. An “or” table extracts several surfaces. This yields cutting, likened to surgery, and capping; texture mapping a cut plane, which clinicians call reformatting, attenuates each normal’s length by the slice data.

Inconsistencies in the printed text

Three passages disagree with each other or with the figures.

• 14 or 15 patterns. The text says the symmetries reduce 256 cases “to 14 patterns”, but Figure 3 draws 15 cubes, numbered 0 to 14, of which pattern 0 is empty. The companion patent describes 14 cases needing manual analysis plus one trivial case (Cline & Lorensen, 1987). Later papers count 15 (Nielson & Hamann, 1991; Chernyaev, 1995; Lewiner et al., 2003).
• Labels. On p. 164 a corner at or above the surface value is assigned one and called inside (or on), a corner below it zero and outside; the next sentence reverses the labels, calling the one outside and the zero inside. An edge is intersected when its labels differ either way, so the case table is unaffected (our reasoning).
• Central differences. Eqs. 2–4 divide $D(i+1, j, k) - D(i-1, j, k)$ by the edge length $\Delta x$, not by $2\Delta x$, and likewise for $y$ and $z$. The missing factor $\tfrac12$ is the same on all three axes, so the estimate is twice the gradient, and normalizing removes it; what matters for the direction is dividing by the separate $\Delta x$, $\Delta y$, $\Delta z$, since the CT slices were 1.5 mm apart with 0.8 mm pixels (our reasoning).

A fourth passage is loose but not wrong: complement symmetry, the paper says, leaves cases with zero to four vertices above the value, 128 of them. There are 163 such cases, since a four-vertex case and its complement both have four; 128 is the number of complement pairs (code below).

Checking the symmetry count

The rotational reduction was done “by inspection”. The code applies the 24 rotations of the cube to all 256 cases, with and without complement, and repeats the count with all 48 symmetries, reflections included.

import itertools
import numpy as np

# Cube corners as (x, y, z) in {0, 1}^3; corner c has bit c in the 8-bit case index.
corners = np.array([[c & 1, c >> 1 & 1, c >> 2 & 1] for c in range(8)])
centred = 2 * corners - 1                                  # coordinates in {-1, 1}^3

# The 48 symmetries of the cube are the signed permutation matrices;
# the 24 with determinant +1 are rotations.
mats = []
for perm in itertools.permutations(range(3)):
    for signs in itertools.product((1, -1), repeat=3):
        M = np.zeros((3, 3), int)
        M[range(3), perm] = signs
        mats.append(M)
rotations = [M for M in mats if round(np.linalg.det(M)) == 1]

def corner_map(M):
    """Permutation of the 8 corners induced by the matrix M."""
    image = centred @ M.T
    return [int(np.flatnonzero((centred == p).all(1))[0]) for p in image]

def apply(perm, case):
    return sum(1 << perm[c] for c in range(8) if case >> c & 1)

def classes(group, complement):
    perms = [corner_map(M) for M in group]
    seen, reps = set(), []
    for case in range(256):
        if case in seen:
            continue
        orbit = {apply(p, case) for p in perms}
        if complement:
            orbit |= {255 - o for o in orbit}
        seen |= orbit
        reps.append(case)
    return reps

n_in = lambda case: bin(case).count("1")
print("rotations:", len(rotations), " all symmetries:", len(mats))
print("cases with 0-4 vertices at or above the value:", sum(n_in(c) <= 4 for c in range(256)))
print("complement pairs:", len({min(c, 255 - c) for c in range(256)}))
for name, group in (("rotations", rotations), ("rotations + reflections", mats)):
    for comp in (False, True):
        reps = classes(group, comp)
        print(f"{name:24s} complement={comp!s:5s} classes={len(reps):2d}  "
              f"by vertex count {np.bincount([min(n_in(r), 8 - n_in(r)) if comp else n_in(r) for r in reps]).tolist()}")

Output (NumPy 2.3.5):

rotations: 24  all symmetries: 48
cases with 0-4 vertices at or above the value: 163
complement pairs: 128
rotations                complement=False classes=23  by vertex count [1, 1, 3, 3, 7, 3, 3, 1, 1]
rotations                complement=True  classes=15  by vertex count [1, 1, 3, 3, 7]
rotations + reflections  complement=False classes=22  by vertex count [1, 1, 3, 3, 6, 3, 3, 1, 1]
rotations + reflections  complement=True  classes=14  by vertex count [1, 1, 3, 3, 6]

Rotations and complement give 15 classes, the empty one included, with one, one, three, three and seven classes having zero to four marked corners. In our count from the page image, Figure 3 has the same distribution (pattern 0 empty, 1 with one marked corner, 2–4 with two, 5–7 with three, 8–14 with four). Reflections would merge two mirror-image four-corner classes and leave 14, so the figure’s 15 cubes are the classes under rotations plus complement, without reflections.

Results

The implementation was written in C and ran on Sun workstations under Unix, VAX computers under VMS and an IBM 3081 under IX/370. Models were displayed with an in-house z-buffer program or a General Electric Graphicon 700, at 10,000 triangles per second. On a VAX 11/780, model creation took from 100 seconds for 64 × 64 × 48 SPECT data to 30 minutes for 260 × 260 × 93 CT studies; the IBM 3081 was twelve times faster. Triangle counts are proportional to surface area.

Three case studies were rendered at 512 × 512 without antialiasing:

• CT. A head study of a twelve-year-old boy with a hole in the skull near the nose, 93 axial slices 1.5 mm thick with 0.8 mm pixels, provided by D. C. Hemmy of the Medical College of Wisconsin. The bone surface has 550,000 triangles and the soft-tissue surface 375,000 (Figures 8–11). A pierced ear and the impression of adhesive tape show the resolution; Figure 11 adds a texture-mapped sagittal cut.
• MR. An adult male volunteer, 128 coronal slices 1.9 mm thick from a flow-compensated 3D Fourier-transform sequence acquired in 9 minutes. Figure 12 shows sagittal cut planes, texture-mapped from the coronal data and intersected with the skin surface; the largest model has 330,000 triangles.
• SPECT. 29 coronal slices of 64 × 64 pixels of the heart, a low-resolution test; the surface of the blood pool in the diastolic heart has 5,000 triangles (Figure 13). The §6 timing quotes 64 × 64 × 48 SPECT data, so it may refer to a different study.

Data from analytic functions are mentioned but not shown. The evaluation is visual: no geometric accuracy is measured, no earlier method is compared, and timings are tied to data sizes rather than to the figures.

Impact

The cube-by-cube table became a common way to mesh a sampled scalar field, first in visualization and later in 3D vision:

• 1987. General Electric held U.S. Patent 4,710,876 on the method, filed by Cline and Lorensen on June 5, 1985 and granted on December 1, 1987. Patents filed before June 8, 1995 run for the greater of twenty years from filing and seventeen from grant (USPTO, MPEP §2701), so its term ended by June 2005 at the latest (our computation).
• 1988. Levoy’s direct volume rendering paper (technical-report version) cites marching cubes as a local surface detector applied at every sample, yielding many voxel-sized polygons, and argues that such binary surface decisions create spurious surfaces or holes where features are small or poorly defined.
• 1995. Chernyaev describes marching cubes as well known and usually considered the basic method for surface rendering in medical applications.
• 1996. Curless and Levoy meshed their fused range data with marching cubes and a disambiguating table, thanking Lorensen for his tables (§5.2.3, acknowledgements); TSDF fusion still ends this way.
• 1998. ACM reprinted the paper in Seminal Graphics: Pioneering Efforts That Shaped the Field.
• 2011. KinectFusion named marching-cubes-type extraction as one of two ways to view a signed distance surface, choosing raycasting for its live loop (Newcombe et al., 2011).
• 2019. Occupancy Networks meshed a learned occupancy function with marching cubes after hierarchical refinement (Mescheder et al., 2019), as is common for neural implicit representations.
• 2020. Lorensen published a history of the algorithm in IEEE Computer Graphics and Applications.

Limitations

The authors named the output size, over 500,000 triangles, reduced with cut planes, connectivity and filtering at some loss of resolution; that MR data need image processing to select a structure; and that the triangle count approaches the number of display pixels as resolution grows.

The best-known flaw was found later: on a face whose inside corners lie on one diagonal, the table built by complement symmetry can join the crossings differently in the two cubes sharing the face, leaving a hole. Nielson and Hamann (1991) showed such a hole, credited Dürst with pointing out the problem, proposed the asymptotic decider, and found ambiguous configurations in a CT volume, though rarely (their Table 1). Chernyaev (1995) needed 33 configurations instead of 15 to match the trilinear interpolant, and Lewiner et al. (2003) completed his table with a manifold guarantee. The marching cubes article explains these fixes and the method’s other weaknesses.

What Came After

The conclusion already pointed to a successor by the same authors, dividing cubes, which outputs points instead of triangles at a density chosen to cover the raster display; they argued that both give the same image quality because shading governs perceived quality. Direct volume rendering dropped the intermediate geometry altogether (Levoy, 1988). Surface extraction became the last stage of reconstruction, meshing fused signed distance volumes, Poisson indicator functions and learned fields into polygon meshes.