Marching Cubes and the Case Table That Cannot Be Interpolated

An isosurface is the set of points where a scalar field equals a chosen level. For a field given only at grid samples, that set is defined implicitly through the interpolant, and extracting it means finding the polygonal approximation of the zero set of $f - c$. Marching cubes does this one cell at a time, exploiting the fact that inside a cell the interpolant is a trilinear polynomial whose behaviour is completely determined by the signs of the eight corner values relative to the level.

That reduction is the whole design: the geometry inside a cell depends only on which corners lie above the level, never on the magnitudes of the values. There are $2^8 = 256$ sign patterns, so a table indexed by the pattern can hold the complete answer, and extraction becomes a table lookup per cell plus a linear interpolation per crossed edge. The cost of the reduction is that the table necessarily discards magnitude information, and one consequence of that discard is an ambiguity that no table indexed only by sign can resolve.

From eight signs to a topology

Enumerate the corners of a cube in a fixed order, for example with corner 0 at the origin and the numbering advancing along $x$, then $y$, then $z$. Build an index byte by testing each corner against the level: bit $i$ set if the value at corner $i$ exceeds $c$. A cube is crossed if its index is neither 0 nor 255; a cube entirely on one side contains no surface and is skipped.

An edge of the cube connects two corners that differ in exactly one coordinate. An edge is crossed when its two endpoints have opposite sides, which is an XOR of their bits. There are twelve edges and at most twelve crossed, so the index byte’s bits determine the crossed-edge set. The table entry for each of the 256 indices is a list of triangles, each given as three edge identifiers in the range 0..11.

The number of distinct configurations is far smaller than 256. Under the symmetries of the cube (24 rotations, reflection, and complementation of the index) the configurations collapse to 15 classes. This is why the table can be small and why it can be checked by hand: fourteen base patterns and their complements, with rotations and reflections mapping every index onto one of them. Up to five triangles are needed for a single cell.

Within a cell, the vertices lie on crossed edges, and their positions come from a one-dimensional root find of the trilinear interpolant along the edge. Since an edge has only two samples and the trilinear interpolant restricted to an edge is the linear interpolant of those two samples, the root is exact:

$$t = \frac{c - f_a}{f_b - f_a}, \qquad \mathbf{p} = \mathbf{p}_a + t,(\mathbf{p}_b - \mathbf{p}_a).$$

When $f_b - f_a$ is small the division is ill-conditioned and $t$ jumps wildly for tiny changes in the level, which produces the characteristic needle artefacts on near-degenerate edges. The practical guards are to reject edges whose gradient magnitude is below a tolerance, or to clamp $t$ into $[\epsilon, 1-\epsilon]$.

CORNERS = [(x, y, z) for z in (0, 1) for y in (0, 1) for x in (0, 1)]
EDGES = [(a, b) for a in range(8) for b in range(a + 1, 8)
         if bin(a ^ b).count("1") == 1]          # the twelve cube edges

def cell_index(vals, level):
    idx = 0
    for i, v in enumerate(vals):
        if v > level:
            idx |= 1 << i
    return idx

def edge_point(vals, coords, edge, level):
    a, b = EDGES[edge]
    fa, fb = vals[a], vals[b]
    denom = fb - fa
    t = 0.5 if abs(denom) < 1e-12 else (level - fa) / denom
    t = min(max(t, 0.0), 1.0)
    pa, pb = coords[a], coords[b]
    return tuple(pa[k] + t * (pb[k] - pa[k]) for k in range(3))

The EDGES construction by popcount of the XOR is a way of writing the twelve edges without transcribing a list, and it makes the correspondence between corner numbering and edge numbering explicit: edge $e$ is the pair whose index bytes differ in one bit.

What the sign-only entry cannot determine

Consider a face of the cube with two diagonally opposite corners above the level and two below. The surface must cross that face, and the trace of the surface on the face is a curve between two crossed edges, but there are two topologically distinct ways to connect them: the curve can cut off one corner or the other. Both are consistent with the interpolant restricted to the face, because the signs alone do not determine which.

The face-restricted trilinear interpolant is a bilinear function, and its zero set is a pair of hyperbolic arcs. In the alternating configuration those arcs can separate the two corners above the isovalue from each other, or the two below, and both are consistent with the signs. The deciding quantity is the value of the bilinear function at its interior critical point, the saddle. Write the four face corners as $f_{00}, f_{10}, f_{01}, f_{11}$ in local face coordinates and let $g = f - c$ be the offset from the isovalue. The bilinear form is $g(u,v) = a + bu + cv + duv$ with $a = g_{00}$, $b = g_{10} - g_{00}$, $c = g_{01} - g_{00}$, $d = g_{00} - g_{10} - g_{01} + g_{11}$. Its critical point is at $u = -c/d$, $v = -b/d$, and the saddle value there is

$$S = a - \frac{bc}{d} = \frac{g_{00},g_{11} - g_{10},g_{01}}{g_{00} + g_{11} - g_{10} - g_{01}}.$$

The sign of $S$ decides the connectivity: if the saddle lies above the isovalue, the region above it is joined through the middle of the face and the two corners below the isovalue are cut off separately. Note that the numerator alone is not the test, because the denominator can be negative; the ratio is what carries the sign. Nielson and Hamann’s asymptotic decider is this comparison, and it needs the magnitudes, not only the signs.

A table indexed only by the corner signs cannot make this decision at all. The classic symptom in a mesh built without the test is a hole: neighbouring cells disagree about how the surface passes through their shared face, and the triangles do not meet. Visibility of the defect depends on the geometry. In a rendered image a hole may read as a background color showing through from behind the surface, which at the far side of an object looks like a dark spot rather than an obvious gap.

The complementary symptom is a false connection, where the mesh joins two lobes of the surface that the interpolant keeps separate. Both symptoms change the genus of the extracted surface, so any subsequent computation that depends on topology, such as computing a volume, a connectivity count, or an Euler characteristic, inherits the error.

The two standard repairs are to add the asymptotic decider to the case table, effectively making the entry a function of one extra bit per ambiguous face, and to switch to a decomposition that has no ambiguity by construction. Marching tetrahedra splits each cube into six tetrahedra and contours each one; a tetrahedron has four corners and the interpolant restricted to it is linear, so the zero set is a plane and the case is unique. The price is more triangles for the same surface, roughly a factor of two to three in triangle count, and a surface whose facets follow the tetrahedral decomposition rather than the cube’s, which changes the appearance of the extracted mesh. A cube-based table with all 256 entries verified against the decider avoids both the extra triangles and the holes, at the cost of a larger table.

Welding, normals, and the cost of the mesh

The naive output is a triangle soup: each cell emits its own vertices and neighbouring cells emit duplicates on the shared edges. A cube edge is shared by four cells, so a vertex on a regular interior edge appears up to four times. Deduplicating with a hash keyed by the edge identity, a pair of global corner indices, is the usual approach; the alternative is to reuse the vertex index cached from a previous cell in scan order, which works for a structured grid and not for an unstructured one.

Vertex normals for shading come from the gradient of the scalar field, evaluated at the vertex position by central differences on the grid. The gradient of the trilinear interpolant has a different value on each side of a cell face, so normals computed from it show the cell structure. Smoothing the field, increasing the sampling density, or averaging the gradients of the incident cells all reduce the crease, and all change the surface that is being described. If the field is analytic, evaluating the analytic gradient at the interpolated position is more accurate than differencing the samples.

The output size grows as the surface area divided by the square of the cell size, so refining a grid by a factor of two along each axis roughly quadruples the triangle count for a smooth surface. Memory for the mesh, at 12 bytes per vertex position plus 12 for a normal plus three 4-byte indices per triangle, is often comparable to or larger than the volume it came from. Streaming extraction, where cells are visited in slabs and only the completed triangles are kept, is the common way to keep both the volume and mesh in memory simultaneously; the cost is that a cell’s triangles depend on the sign of its neighbours, so the sweep must retain one layer of corner values beyond the slab being emitted.

Values, not just signs

The case-table reduction is exact only for the topology question it answers. Everything else about the surface, position, orientation, and curvature, depends on the values, and the interpolation step is where all of that information re-enters. A consequence is that two datasets with identical sign patterns produce identical meshes up to the vertex positions, so an extraction from a low-resolution field and from a high-resolution one can have the same connectivity but very different geometry.

The dependence also explains why the level set method and marching cubes are different tools for related problems. Marching cubes extracts one level from a fixed grid. Level set methods advect the implicit function and re-extract it at each step, so the surface can move between cells without a new mesh being generated from scratch. The extraction step inside such a solver has the same table and the same ambiguity.

The colors applied to the extracted mesh, whether by the field’s own values or by a second quantity sampled onto the surface, carry their own conventions; see Perceptually Uniform Colormaps and the Failure of Rainbow Ramps for the constraints on that mapping. Where the field is not a grid at all but a set of scattered points, the surface question changes completely, as covered in Scattered Data Interpolation and Radial Basis Functions.