Marching Cubes and the Ambiguity at a Cell Face

A scalar field sampled on a grid is easy to store and easy to evaluate pointwise. A surface extraction algorithm converts that grid into triangles approximating the level set at a chosen value. Marching cubes does this by processing each cell of the grid independently, looking at the sign of the field at its eight corners, and emitting triangles from a precomputed table. The design is local, which is what makes it fast and parallel, and locality is also the origin of its well-known topological failure.

The case table and why it is smaller than it looks

Each of the eight corners of a cell is either inside or outside the level set. That gives $2^8 = 256$ configurations. Enumerate them literally and the table is 256 entries, each listing up to five triangles, with each triangle’s vertices on specific cell edges. The original 1987 paper presents it that way.

Most of the 256 are redundant. Inverting all eight bits — swapping inside and outside — produces the same surface with reversed orientation, so 256 collapses to 128. Rotating a cell by the 24 symmetries of the cube groups the remaining cases into 15 topologically distinct configurations, of which one is the empty case with no triangles. So a hand-written implementation can hold 15 base cases plus a table of geometric transforms, and the Lookup tables in most libraries are 256 entries only because the transform machinery costs more instructions than the table costs bytes.

The geometry inside each case comes from linear interpolation along the cell’s twelve edges. A corner at position $p_a$ with value $f_a$ and a neighbour at $p_b$ with value $f_b$ straddling the level $L$ has a crossing point at $t = (L - f_a)/(f_b - f_a)$ along the edge. That linear interpolation is the only approximation in the algorithm: the surface is piecewise linear within each cell, so a curved isosurface becomes a facet where the true surface would be curved, and the error is bounded by the cell size times the field’s second derivative on that cell.

# canonical corner numbering: corner i has bits (x, y, z) = (i & 1, i >> 1 & 1, i >> 2 & 1)
EDGES = [(0, 1), (2, 3), (4, 5), (6, 7), (0, 2), (1, 3),
         (4, 6), (5, 7), (0, 4), (1, 5), (2, 6), (3, 7)]

def cell_triangles(vals, corners, level):
    """Triangles for one cell. TABLE maps the 8-bit case to edge-id triples."""
    case = sum(1 << i for i, v in enumerate(vals) if v < level)
    pts = {}
    for eid, (a, b) in enumerate(EDGES):
        if (vals[a] < level) != (vals[b] < level):        # edge straddles the level
            t = (level - vals[a]) / (vals[b] - vals[a])
            pts[eid] = tuple(corners[a][k] + t * (corners[b][k] - corners[a][k])
                             for k in range(3))
    return [[pts[e] for e in tri] for tri in TABLE[case]]

The table is the substance; the loop is a few lines. Note that the table must be orientation-consistent as well as positionally correct, because a triangle wound the wrong way produces an inward-facing normal and a hole in any downstream signed volume computation.

Vertices are shared, faces are not

A naive implementation produces each cell’s triangles independently and appends them to a global list. Adjacent cells compute the same crossing point along their shared face independently, from their own corner values, and the two computations agree exactly because they use the same two endpoint values and the same interpolation formula. In exact arithmetic the duplicate vertices coincide; in floating point they are bitwise identical only if the arithmetic order is identical. Since it is, the vertices coincide, and welding them by exact position is safe.

That is an unusual property. Most mesh welding needs a tolerance; here the interpolated coordinate along a shared edge is produced by the same expression from the same inputs, so the position hash matches exactly. Emitting an index per (edge, cell) pair rather than per position triples the vertex count for no benefit, and the standard fix is to key the vertex cache on the grid edge identity — the pair of lattice indices at its endpoints — rather than on geometry. A grid of $N^3$ cells has roughly $3N^3$ interior edges, and every vertex produced lies on one of them.

The triangle count that a surface extraction produces is not controlled by the input. A sphere at radius $r$ sampled on a grid of spacing $h$ has an area $4\pi r^2$, and marching cubes covers it with triangles of area roughly $h^2/2$, so the count is around $8\pi r^2 / h^2$. Halving the grid spacing quadruples the triangle count, and the extraction becomes a memory bottleneck long before the sampling does.

The ambiguous face and the two surfaces it allows

The independence of cells is exactly where the algorithm breaks. Consider one face of a cell with four corners, two inside diagonally opposite each other and two outside. On that face, the contour separates the inside corners. There are two ways to connect them, and both are consistent with the corner signs. The surface passing through the cell can either connect the two inside corners through a saddle, or keep them separated.

Marching cubes makes this decision implicitly, once per case in the table, and the decision is made independently for each cell. Two adjacent cells sharing that face may choose differently. Each cell’s surface is individually valid, and they meet along the shared face’s four edges. But the connection pattern inside one cell can contradict the pattern inside the other, and the resulting mesh has a topological defect: a hole, or a pair of intersecting sheets, depending on the case. The surface can even change genus in a way the true isosurface does not.

The ambiguity is not a rare corner case. It appears at every place where the field has a saddle point inside a cell face, which is exactly the neighbourhood of the medial axis between two close surfaces. A thin wall separating two blobs, a nearly-touching pair of spheres, the gap between fingers: all produce ambiguous configurations.

Resolving it consistently

Two classes of fix exist. The first keeps the cell-local structure and makes the decision consistently. The asymptotic decider computes, from the four corner values, where the bilinear interpolant on the face has its saddle and whether that value is inside or outside the level. That test decides which pairing of inside corners the true bilinear surface has, and both cells sharing the face make the same decision because they have the same four corner values. The decision separates the cases where a tunnel passes through the face from those where a wall does.

The second approach abandons the cell table entirely. Marching tetrahedra splits each cube into six tetrahedra and applies a trivial three-case rule inside each. A tetrahedron has four corners and therefore sixteen sign patterns, two of which place all four corners on one side and produce nothing. No remaining configuration is ambiguous, because a tetrahedron has no face with diagonally opposite same-sign corners, so the three topological cases cover it completely. The price is more triangles: a cube split into six tetrahedra needs several triangles to cover a surface that one cube’s table might cover with two, and the added diagonals produce slivers along the split.

Dual contouring takes a different route. It places one vertex per cell that straddles the surface, positioned by minimizing a quadratic error against the cell’s edge crossings, and connects those vertices around each grid edge. The topology of the output is then determined by the grid edges rather than by the cell faces, so the ambiguity has no place to arise: each grid edge is either crossed or not. The output is a quad mesh with better-shaped elements than marching cubes’ triangles, and it reproduces sharp features more accurately because the vertex position is not constrained to lie on an edge. The cost is a solve per cell and a mesh whose topology must be checked for manifoldness separately, since a cell with two disjoint surface components produces one vertex that both components share.

Sign conventions and the practical edge cases

The extracted mesh is only as good as the sign convention of the input. If the field is a sampled signed distance grid stored with an implicit truncation band, cells outside the band have no sign information, and extraction at a level within the band is fine while extraction at a level outside it invents surfaces from garbage. Truncated fields store values only within a narrow band around the zero level, and any isovalue outside that band requires resampling.

Samples exactly equal to the level value are another edge. The comparison f < L versus f <= L changes which corners are inside, and a grid plane coincident with a flat surface can produce zero-area triangles or no triangles at all depending on the choice. The usual convention is to perturb the level by a tiny amount, which is honest but leaves the choice visible in the output at exactly the flat regions where it matters least.

Finally, the topology of the result is not guaranteed to be that of the true isosurface, even with a consistent ambiguity rule. The resolution of the grid sets a scale below which features are invisible, and two surfaces closer than one cell diagonal merge into one. That is a sampling limitation rather than an algorithmic one, and it has no fix other than refining the grid.