Marching Cubes — Extracting Isosurfaces from Volumetric Data

MRI scans, fluid level-set simulations, metaballs, and procedural terrain all share a common representation: a three-dimensional scalar field where the interesting surface is the set of points where the field equals a threshold value (an isosurface). Marching Cubes, invented by Lorensen and Cline at GE in 1987, converts any such implicit surface into a triangular mesh suitable for rasterized rendering — and despite being nearly 40 years old, it remains the algorithm of choice in medical imaging software, 3D game engines, and real-time fluid renderers worldwide.

1. Isosurfaces and Implicit Representations

A scalar field F: ℝ³ → ℝ assigns a real value to every point in 3D space. An isosurface at isovalue c is the set S = {(x, y, z) : F(x,y,z) = c}. Examples:

Metaball: F(x,y,z) = Σ rᵢ² / (|p − pᵢ|²) with isovalue c = 1; the surface is where contributions sum to 1.
Signed distance field (SDF): F gives the signed distance to the nearest surface; isosurface at c = 0 is the surface itself.
MRI tissue density: F(x,y,z) = Hounsfield unit of voxel; isosurface at c ≈ 400 HU extracts bone.
Fluid level set: ϕ = 0 is the air-water interface; negative ϕ is liquid interior.

All of these produce isosurfaces that cannot be directly rasterized — they must first be converted to an explicit triangle mesh.

2. The Voxel Grid

Marching Cubes samples the scalar field on a regular 3D grid of voxels (volumetric pixels). Each voxel is a cube with 8 corner samples. The algorithm processes one cube at a time — "marching" through each cube in the grid:

For a cube with corners c_0 \dots c_7 and scalar values v_0 \dots v_7: 1. Classify each corner as INSIDE (v_i < isovalue) or OUTSIDE (v_i > isovalue) 2. Form an 8-bit index: cubeIndex = Σ (v_i < isovalue ? 1 : 0) << i 3. Look up edgeTable[cubeIndex]: which of the 12 edges are intersected 4. For each intersected edge, interpolate the exact crossing position 5. Look up triTable[cubeIndex]: how to connect edge midpoints into triangles Result: 0 to 5 triangles per cube, seamlessly joined across cube boundaries.

3. The 256-case Lookup Table

With 8 corners each classifiable as inside or outside, there are 2⁸ = 256 unique corner configurations. However, complementarity and rotational symmetry reduce these to only 15 fundamentally distinct cases:

All outside: no triangle.
One corner inside: 1 triangle.
One edge inside: 1 triangle (different orientation).
Two adjacent corners: 2 triangles.
Two opposite corners on same face: 2 triangles (S-shape).
Two diagonally opposite corners: 2 triangles (crossing).
Three corners on same face: 2 triangles + 1.
Three corners L-shape: 3 triangles.
Complete face inside: 2 triangles forming a quad.
Face + adjacent corner: 3 triangles.
Two faces: 4 triangles.
S-shape diagonal: 4 triangles.
Three faces: 4 triangles.
Almost full (7 corners): 1 triangle.
All inside: no triangle.

All 256 cases are generated by applying the 15 base cases under rotations and reflections. The resulting triTable[256][16] encodes which edges form each triangle (edge indices 0–11; −1 terminates the list).

4. Trilinear Interpolation

Marching Cubes places triangle vertices not at grid corners but at interpolated crossing points on each edge. For an edge between corners i and j with scalar values vᵢ and vⱼ:

t = (isovalue - v_i) / (v_j - v_i) // linear parameter t \in [0, 1] p = p_i + t * (p_j - p_i) // interpolated 3D position No interpolation (t = 0.5) produces "blocky" Minecraft-like surfaces. Trilinear interpolation produces smooth curved surfaces matching the field.

function interpolateEdge(p1, p2, v1, v2, iso) {
  if (Math.abs(v1 - iso) < 1e-5) return p1;
  if (Math.abs(v2 - iso) < 1e-5) return p2;
  const t = (iso - v1) / (v2 - v1);
  return {
    x: p1.x + t * (p2.x - p1.x),
    y: p1.y + t * (p2.y - p1.y),
    z: p1.z + t * (p2.z - p1.z)
  };
}

5. Ambiguous Cases and MC33

The original 1987 marching cubes has ambiguous face configurations: when opposite corners of the same face are both inside (or both outside), there are two topologically distinct ways to triangulate the cube, and the naïve lookup table may choose inconsistently — causing holes in the final mesh.

6 ambiguous face cases out of 15 base configurations. If adjacent cubes resolve the same ambiguous face differently, the resulting mesh has a seam — a topological hole that breaks watertightness and causes rendering artefacts.

MC33 (Chernyaev, 1995) extends the lookup table to 33 distinct cases by specifying a consistent face-disambiguation strategy. It guarantees a topologically correct, hole-free mesh. Most production implementations (OpenVDB, VTK, Three.js MarchingCubes) use MC33 tables.

6. Computing Surface Normals

Smooth normals are critical for shading. Two approaches:

Gradient-based normals

The normal to an isosurface at point p is proportional to the gradient of the scalar field: n̂ = ∇F / |∇F|. Approximate the gradient by central differences on the voxel grid:

\partialF/\partialx \approx (F(x+h, y, z) - F(x-h, y, z)) / (2h) Per-corner gradients are then interpolated to triangle vertices using the same t parameter as position interpolation, producing smooth normals independent of triangle density.

Flat facet normals

Take the cross product of two triangle edges: n = (v₁ − v₀) × (v₂ − v₀). Fast but produces faceted appearance on low-resolution grids.

7. Applications

Medical imaging: CT/MRI volumetric datasets reconstructed to 3D meshes for surgical planning and visualization (ITK-SNAP, 3D Slicer, VTK).
Fluid simulation: Level-set fluid simulators (OpenVDB, Houdini FLIP) run marching cubes each frame to extract the liquid surface mesh for rendering.
Terrain generation: Voxel-based games (Minecraft-like engines) use modified marching cubes for smooth terrain with overhangs and caves impossible in heightfield terrain.
Metaballs: Organic, blobby shapes for character animation rigs and special effects — the classic demo-scene effect of merging spheres.
Neural SDFs: NeRF and neural implicit surface networks (NeuS, VolSDF) output signed distance fields that marching cubes converts to meshes for mesh-based rendering pipelines.

Performance: A 256³ voxel grid contains 16 million cubes. A simple CPU Marching Cubes implementation processes roughly 50–100 million cubes per second on a modern core. GPU implementations (compute shaders with workgroup prefix sums for compacting output triangles) achieve 10–100× speedup, enabling real-time marching cubes in game engines.

💧 Try the Fluid Simulation →