Marching Cubes

How marching cubes works

You have a scalar field: a function f(x, y, z) that returns a number at every point in space. It might be CT-scan density, a 3D noise function for terrain, or a signed distance field where the value is "how far am I from the surface." You want the surface where f = c for some chosen isovalue c — the isosurface — turned into triangles a GPU can draw.

The field is continuous, so there are infinitely many points to check. Marching cubes makes the problem finite by sampling f on a regular 3D grid and then walking — "marching" — through every cube cell formed by 8 neighboring grid samples. For each cube it asks one question per corner: is this corner inside the surface (f < c) or outside (f ≥ c)?

Eight corners, each in one of two states, gives a single 8-bit number from 0 to 255. That number is the cube index. The surface must pass between any two adjacent corners that disagree, so it crosses certain edges of the cube. A precomputed table maps each of the 256 cube indices to exactly which edges are crossed and how to connect the crossing points into triangles:

1. Sample. Read f at the 8 cube corners.
2. Classify. Build the 8-bit cube index from which corners are inside.
3. Look up edges. The edgeTable[cubeIndex] entry tells you which of the 12 edges the surface intersects.
4. Interpolate. On each crossed edge, find where the surface crosses by interpolating between the two corner values.
5. Triangulate. The triTable[cubeIndex] entry lists triples of edges, each triple becoming one triangle.

Two cube indices are trivial: index 0 (all corners outside) and index 255 (all inside) emit no triangles — the surface doesn't pass through. Everything in between emits 1 to 5 triangles. Run this over the whole grid and the triangles from neighboring cubes stitch into a watertight mesh.

Why interpolation matters

The naive version places each surface vertex at the midpoint of a crossed edge. This works, but produces a blocky, staircased mesh — every vertex sits on a grid line, so the surface looks like LEGO. The fix is one line of arithmetic.

On an edge between corner p0 (value v0) and corner p1 (value v1), the surface crosses where the value equals the isolevel c. Assuming the field is linear along the edge:

t = (c - v0) / (v1 - v0)
vertex = p0 + t * (p1 - p0)

If c is exactly halfway between v0 and v1, t = 0.5 and you get the midpoint. But if one corner is barely inside and the other deeply outside, the vertex slides toward the barely-inside corner — and the mesh becomes smooth. This is the single highest-impact line in the whole algorithm. Vertex normals come almost for free: estimate the gradient of f at each corner by central differences, interpolate it the same way as the position, and you have smooth shading without ever computing a triangle face normal.

When to use it — and when not to

• Medical and scientific volumes. CT and MRI scanners output dense 3D density grids; marching cubes is the standard way to extract bone or organ surfaces. This is literally what it was invented for.
• Procedural and destructible terrain. Games like No Man's Sky-style planets and voxel sandboxes store a density field and re-mesh chunks when you dig. Marching cubes gives smooth terrain that a pure cube/Minecraft voxel mesher cannot.
• Metaballs and fluid surfaces. Sum of radial falloff functions, thresholded — marching cubes turns the implicit blob into renderable geometry every frame.
• Reconstructing surfaces from point clouds after fitting a signed distance field (e.g. Poisson reconstruction outputs an indicator field you then march).

Reach for something else when you need sharp creases (use dual contouring or dual marching cubes), when you must guarantee a manifold, hole-free mesh for 3D printing (use marching cubes 33 or a topology-correct variant), or when the field is adaptive / sparse and a uniform grid would waste memory (use an octree-based or out-of-core scheme). If your data is a height field rather than a true volume, you don't need marching cubes at all — a 2D grid triangulation is far cheaper.

Marching cubes vs other isosurface methods

	Marching cubes	Marching tetrahedra	Dual contouring	Surface nets	Marching cubes 33
Vertex placement	On cube edges	On tet edges	One per cell (QEF solve)	One per cell (averaged)	On cube edges
Sharp features	Rounded off	Rounded off	Preserved	Mostly rounded	Rounded off
Topology correct	No (ambiguous faces)	Yes (simplices)	Yes with care	No guarantee	Yes
Cases / table	15 base, 256 entries	3 base, 16 entries	gradient-driven, no table	none	33 base, 256 entries
Triangle count	baseline	~2–3× more	quads, fewer	fewest	baseline
Needs gradients	No (only for normals)	No	Yes (Hermite data)	No	No
Patent history	GE 1985, expired 2005	Patent-free workaround	2002, free	1998, free	free

The headline tradeoff: marching cubes is dead simple and fast (one table lookup per cell) but smears sharp edges and can produce cracks on ambiguous faces. Marching tetrahedra dodges the ambiguity entirely by splitting each cube into 5 or 6 tetrahedra — a tetrahedron has no ambiguous faces — at the cost of 2–3× more triangles. Dual contouring is the modern choice when sharp features matter, but it needs surface normals (Hermite data) and a small least-squares solve per cell.

What the numbers actually say

• Cost is O(N³) to scan, but O(N²) in output. A grid of N³ cells visits every cell once, yet only cells straddling the surface emit triangles, and those form a 2D shell — about N² of them. A 256³ field of a sphere visits 16.7 million cells but emits on the order of a few hundred thousand triangles.
• 1–5 triangles per active cell. The 256-entry triangle table emits at most 5 triangles for any single cube. Most active cells emit 1 or 2.
• The table is 256 × 16 bytes. The triangle table stores up to 15 edge indices plus a terminator per case — about 4 KB total, trivially cache-resident.
• Embarrassingly parallel. Every cube is independent, so the algorithm maps cleanly to the GPU. A compute-shader implementation can mesh a 128³ field in well under a millisecond on modern hardware — fast enough to re-mesh destructible terrain every frame.
• 15 cases → 256 entries. The 2⁸ sign patterns collapse to 15 distinct shapes under rotation, reflection, and complement; the 256-entry table is those 15 expanded back out so no per-cell symmetry math is needed at runtime.

JavaScript implementation

The full tables are long, so this shows the core loop and edge interpolation. The two lookup tables, EDGE_TABLE (256 entries) and TRI_TABLE (256 rows), are the canonical Paul Bourke tables — paste them in from any reference implementation.

// 8 corner offsets, 12 edges (each edge = pair of corner indices)
const CORNERS = [[0,0,0],[1,0,0],[1,1,0],[0,1,0],[0,0,1],[1,0,1],[1,1,1],[0,1,1]];
const EDGES = [[0,1],[1,2],[2,3],[3,0],[4,5],[5,6],[6,7],[7,4],[0,4],[1,5],[2,6],[3,7]];

// Linear interpolation: where does the isolevel cross this edge?
function vertexOnEdge(c, p0, p1, v0, v1) {
  if (Math.abs(c - v0) < 1e-6) return p0;
  if (Math.abs(c - v1) < 1e-6) return p1;
  if (Math.abs(v0 - v1) < 1e-6) return p0;
  const t = (c - v0) / (v1 - v0);
  return [p0[0] + t*(p1[0]-p0[0]), p0[1] + t*(p1[1]-p0[1]), p0[2] + t*(p1[2]-p0[2])];
}

// March one cube whose minimum corner is at (x,y,z), spacing s
function marchCube(field, x, y, z, s, c, out) {
  const val = [], pos = [];
  for (let i = 0; i < 8; i++) {
    const [dx,dy,dz] = CORNERS[i];
    pos[i] = [x + dx*s, y + dy*s, z + dz*s];
    val[i] = field(pos[i][0], pos[i][1], pos[i][2]);
  }
  // Build the 8-bit cube index
  let cubeIndex = 0;
  for (let i = 0; i < 8; i++) if (val[i] < c) cubeIndex |= (1 << i);

  const edges = EDGE_TABLE[cubeIndex];
  if (edges === 0) return;            // fully in or fully out — no surface

  // Interpolate a vertex on every crossed edge
  const vert = [];
  for (let e = 0; e < 12; e++) {
    if (edges & (1 << e)) {
      const [a,b] = EDGES[e];
      vert[e] = vertexOnEdge(c, pos[a], pos[b], val[a], val[b]);
    }
  }
  // Emit triangles per the triangle table (triples of edge indices, -1 terminates)
  const tri = TRI_TABLE[cubeIndex];
  for (let i = 0; tri[i] !== -1; i += 3) {
    out.push(vert[tri[i]], vert[tri[i+1]], vert[tri[i+2]]);
  }
}

// Scan the whole grid
function marchingCubes(field, min, max, s, c) {
  const out = [];
  for (let z = min; z < max; z += s)
    for (let y = min; y < max; y += s)
      for (let x = min; x < max; x += s)
        marchCube(field, x, y, z, s, c, out);
  return out;   // flat array of triangle vertices
}

// Example field: a sphere of radius 1 (isolevel 0 of a signed distance field)
const sphere = (x,y,z) => Math.sqrt(x*x + y*y + z*z) - 1;
const tris = marchingCubes(sphere, -1.5, 1.5, 0.1, 0);

Note that the corner ordering of CORNERS, the bit each corner sets in cubeIndex, and the edge numbering in EDGES must all agree with the table you paste in. Mismatched conventions are the number-one source of garbled output — the algorithm runs cleanly but the triangles connect the wrong points.

Python implementation

In practice nobody hand-codes this in Python for production — you call a library. The reference one-liner is scikit-image's marching_cubes, which implements the topology-correct Lewiner variant:

import numpy as np
from skimage import measure

# Build a 3D scalar field: signed distance to a sphere of radius 30
n = 100
x, y, z = np.mgrid[:n, :n, :n]
center = n // 2
field = np.sqrt((x-center)**2 + (y-center)**2 + (z-center)**2) - 30.0

# Extract the isosurface at f = 0
verts, faces, normals, values = measure.marching_cubes(
    field,
    level=0.0,            # the isovalue c
    spacing=(1.0, 1.0, 1.0),
    gradient_direction='descent',
)

print(verts.shape, faces.shape)   # (~28000, 3) vertices, (~56000, 3) faces

# Save as an STL you can 3D-print or load in Blender
from stl import mesh                      # numpy-stl
m = mesh.Mesh(np.zeros(len(faces), dtype=mesh.Mesh.dtype))
for i, f in enumerate(faces):
    for j in range(3):
        m.vectors[i][j] = verts[f[j]]
m.save('sphere.stl')

And the bare algorithm, mirroring the JavaScript above, if you want to understand every step rather than call a library:

CORNERS = [(0,0,0),(1,0,0),(1,1,0),(0,1,0),(0,0,1),(1,0,1),(1,1,1),(0,1,1)]
EDGES   = [(0,1),(1,2),(2,3),(3,0),(4,5),(5,6),(6,7),(7,4),(0,4),(1,5),(2,6),(3,7)]

def vertex_on_edge(c, p0, p1, v0, v1):
    if abs(v0 - v1) < 1e-9:
        return p0
    t = (c - v0) / (v1 - v0)
    return tuple(p0[k] + t * (p1[k] - p0[k]) for k in range(3))

def march_cube(field, x, y, z, s, c, out):
    pos = [(x + dx*s, y + dy*s, z + dz*s) for dx, dy, dz in CORNERS]
    val = [field(*p) for p in pos]
    cube_index = 0
    for i, v in enumerate(val):
        if v < c:
            cube_index |= (1 << i)
    if EDGE_TABLE[cube_index] == 0:
        return
    vert = [None] * 12
    for e in range(12):
        if EDGE_TABLE[cube_index] & (1 << e):
            a, b = EDGES[e]
            vert[e] = vertex_on_edge(c, pos[a], pos[b], val[a], val[b])
    tri = TRI_TABLE[cube_index]
    for i in range(0, len(tri), 3):
        if tri[i] == -1:
            break
        out.append((vert[tri[i]], vert[tri[i+1]], vert[tri[i+2]]))

Variants worth knowing

Marching tetrahedra. Split each cube into 5 or 6 tetrahedra before marching. A tetrahedron has only 16 cases and no ambiguous faces, so the mesh is always topologically correct and watertight. The cost is 2–3× more triangles. Historically popular because it sidestepped GE's marching cubes patent.

Marching cubes 33 (Chernyaev). Replaces the original 15-case table with 33 cases that correctly resolve every ambiguous face, guaranteeing a manifold, hole-free surface. The variant most people actually want when correctness matters.

Lewiner marching cubes. A 2003 implementation (the one in scikit-image) that combines topological correctness with efficient table generation. The pragmatic default in scientific Python.

Dual contouring. Places one vertex inside each cell and uses the surface normal at each edge crossing (Hermite data) to solve a small quadratic for the vertex position. This lets the vertex land exactly on a sharp corner or crease — the one thing marching cubes structurally cannot do.

Surface nets / naive dual. Like dual contouring but averages edge crossings instead of solving a QEF. Faster and simpler, produces fewer triangles, but rounds sharp features. A great default for smooth organic surfaces.

Adaptive / octree marching cubes. Run on an octree so flat regions use big cells and detailed regions subdivide. Cuts triangle count dramatically but requires careful crack-patching at level-of-detail boundaries where a coarse cell abuts a fine one.

Common bugs and edge cases

• Mismatched corner/edge/table conventions. The most common failure. If your corner numbering or bit assignment disagrees with the table's, you get a mesh that looks like noise. Copy all three — corners, edges, and tables — from one source.
• Cracks on ambiguous faces. The 15-case table resolves diagonally-split faces inconsistently between neighbors, leaving holes. Use marching cubes 33 or a topology-correct variant for watertight output.
• Division by zero in interpolation. When both corner values equal the isolevel, v1 - v0 is zero. Guard it (return either endpoint) before dividing.
• Sampling on the isovalue exactly. If a grid corner's value equals c exactly, the inside/outside test is ambiguous. Use a strict < consistently, and nudge the isovalue by an epsilon if degenerate triangles appear.
• Wrong winding order. Triangles emitted with inconsistent vertex order produce backwards normals and a surface that lights from inside. Reverse the triple if your renderer culls the visible side.
• Forgetting the grid only sees what it samples. Features thinner than the cell spacing vanish entirely — a wall one voxel thin can disappear or fragment. Sample finely enough to resolve your thinnest feature, or the mesh will silently omit it.
• Recomputing shared edge vertices. Adjacent cubes share edges; the naive loop interpolates the same vertex up to four times, bloating the vertex buffer 3–4×. Cache edge vertices in a hash keyed by edge to weld the mesh.