Game Dev Mechanics: Marching Cubes — How It Works

What Is Marching Cubes?

Imagine a three-dimensional grid of numbers, a scalar field, where every point in space holds a single floating-point value. You want to draw a surface wherever that value crosses a specific threshold. That surface is called an isosurface, and the Marching Cubes algorithm is how you turn it into triangles. Published in 1987 by William Lorensen and Harvey Cline at General Electric, Marching Cubes remains one of the most important algorithms in real-time graphics. It powers metaballs, fluid surface rendering, procedural cave systems, volumetric destruction, and medical CT scan visualization.

The core idea is simple: divide space into a grid of small cubes, visit each cube one at a time, and decide, based solely on the field values at its eight corners, which triangles to place inside it. Assembled together, those triangles form a continuous, smooth surface. The algorithm's name describes exactly what it does: the surface extraction marches through every cube in the grid.

Scalar Fields and Isosurfaces

A scalar field is a function $f(x, y, z) \to \mathbb{R}$ that assigns a real number to every point in 3D space. The isosurface at threshold $\tau$ is the set of all points where the field equals that value:

$$S_\tau = \{ (x,y,z) \mid f(x,y,z) = \tau \}$$

Think of it as a 3D contour line. Common scalar field functions used in games include:

• Metaball fields: Each ball $i$ at center $\mathbf{c}_i$ with radius $r_i$ contributes $f_i(\mathbf{p}) = r_i^2 / \|\mathbf{p} - \mathbf{c}_i\|^2$. Summing these gives the total field. When two balls approach each other, their combined field exceeds the threshold between them, causing the surface to merge organically.
• Signed Distance Fields: $f$ returns signed distance to the nearest surface, negative inside and positive outside. This enables powerful CSG boolean operations between shapes.
• 3D Noise terrain: A 3D Simplex or Perlin noise function creates caves, overhangs, and natural-looking geology that flat heightmaps cannot represent.

Marching Cubes is agnostic to what the field represents. It only cares whether each corner's value is above or below $\tau$.

The Algorithm Step by Step

Step 1: Sample the Grid

Discretize space into a regular grid with $N^3$ sample points, forming $(N-1)^3$ cubes. Evaluate $f$ at every vertex and cache the results. Each cube shares its eight corner values with neighboring cubes, so each value is computed once and reused by up to eight cubes.

Step 2: Compute the Cube Index

For a given cube, compare each of the eight corner values to the threshold $\tau$. Encode the result as an 8-bit integer:

$$\text{cubeIndex} = \sum_{i=0}^{7} \mathbf{1}[f(v_i) \geq \tau] \cdot 2^i$$

The cube index ranges from 0 (all corners below threshold, no surface intersection) to 255 (all corners above threshold, also no intersection). Only the 254 intermediate cases produce geometry.

Step 3: Look Up the Triangle Configuration

Lorensen and Cline pre-computed the triangulation for all 256 configurations. Two lookup tables encode the result:

• edgeTable[256]: A 12-bit bitmask. Bit $k$ is set if edge $k$ of the cube is intersected by the isosurface. A cube has 12 edges; an edge is active when one endpoint is inside and the other is outside.
• triTable[256][16]: For each cube index, a list of up to five triangles encoded as edge index triples, terminated by -1. Each triangle is three edge indices; the actual vertex positions are computed by interpolating along those edges.

By symmetry arguments, the 256 cases reduce to only 15 topologically distinct configurations. The authors exploited rotational and reflective symmetry to derive the full 256-entry table from those 15 base cases.

Step 4: Interpolate Edge Vertices

For each active edge connecting corners $A$ and $B$ with values $f_A$ and $f_B$, find the exact crossing point using linear interpolation:

$$\mathbf{p} = \mathbf{A} + \frac{\tau - f_A}{f_B - f_A}(\mathbf{B} - \mathbf{A})$$

This interpolation is what gives Marching Cubes its smooth appearance. Without it, if you snapped vertices to grid corners, the result would look like blocky voxels. With interpolation, the surface slides smoothly across the edge as the field changes.

Step 5: Emit Triangles and Compute Normals

Use triTable to assemble interpolated edge vertices into triangles. For normals, you have two options. The simpler approach averages face normals at shared vertices (what Three.js's computeVertexNormals() does). The better approach computes the gradient of the scalar field at each vertex using finite differences:

$$\nabla f \approx \left(\frac{f(x+\epsilon,y,z) - f(x-\epsilon,y,z)}{2\epsilon},\; \frac{f(x,y+\epsilon,z) - f(x,y-\epsilon,z)}{2\epsilon},\; \frac{f(x,y,z+\epsilon) - f(x,y,z-\epsilon)}{2\epsilon}\right)$$

Gradient normals produce smooth shading because they reflect the true shape of the underlying field, regardless of triangulation resolution.

Ambiguous Cases and Fixes

The original algorithm has a subtle flaw: certain configurations are ambiguous. When a cube face has two inside and two outside corners in a diagonal pattern, there are two equally valid ways to triangulate it, and neighboring cubes may choose differently, leaving cracks in the mesh. Chernyaev's 1995 Marching Cubes 33 paper extended the base cases to 33 distinct topological configurations, resolving all ambiguities. For game terrain and metaballs, the original table works fine. For watertight medical meshes, MC33 or Dual Contouring is preferred.

Dual Contouring (Ju et al., 2002) takes a fundamentally different approach: instead of placing vertices on cube edges, it places one vertex per cube at the position that best fits the local surface normals using a Quadric Error Function. This preserves sharp edges and corners that Marching Cubes rounds off, making it popular for architectural and hard-surface procedural generation.

Applications in Games

Marching Cubes appears throughout the industry in both obvious and surprising ways:

• Destructible voxel terrain: Games like 7 Days to Die, Teardown, and No Man's Sky store world data as voxel density grids and use Marching Cubes (or Transvoxel) to skin the surface. Removing voxels updates the density field; the mesh regenerates automatically.
• Fluid surfaces: SPH (Smoothed Particle Hydrodynamics) fluid simulators splat particle density into a grid each frame and extract the water surface with Marching Cubes. Games like Noita use variants of this for liquid simulation.
• Metaball enemies and VFX: Organic, morphing blob enemies and slime effects are trivially implemented as animated metaball fields.
• Procedural planets: No Man's Sky generates planet surfaces by evaluating a noise field in 3D (not 2D), enabling caves, arches, and floating islands impossible with heightmaps.

Performance Considerations

Naive Marching Cubes is $O(N^3)$, visiting every cube regardless of whether it contains surface geometry. For large grids or real-time updates, several optimizations matter:

• Active cell skipping: Skip cubes with all corners inside or all outside. For sparse surfaces, only a small fraction of cubes contribute triangles.
• Shared vertex tables: Use a hash map from edge ID to vertex index to deduplicate vertices shared between adjacent cubes. This halves vertex count and enables proper normal averaging.
• GPU Compute Shaders: Each cube is independent, which makes it well-suited for GPU parallelism. A compute shader pass uses prefix-sum compaction to collect active cells and write triangles into a compact buffer in parallel.
• Transvoxel LOD: Eric Lengyel's Transvoxel algorithm adds transition cells at LOD boundaries so coarse and fine grid regions stitch together without cracks, enabling view-distance-scaled terrain.

Code Walkthrough

// Core march loop — iterates over all (N-1)^3 cubes
function march(field, N, threshold) {
  const verts = [];

  for (let iz = 0; iz < N - 1; iz++)
  for (let iy = 0; iy < N - 1; iy++)
  for (let ix = 0; ix < N - 1; ix++) {

    // Sample 8 corner values
    const v = [
      field[idx(ix,   iy,   iz,   N)],  // v0
      field[idx(ix+1, iy,   iz,   N)],  // v1
      field[idx(ix+1, iy+1, iz,   N)],  // v2
      field[idx(ix,   iy+1, iz,   N)],  // v3
      field[idx(ix,   iy,   iz+1, N)],  // v4
      field[idx(ix+1, iy,   iz+1, N)],  // v5
      field[idx(ix+1, iy+1, iz+1, N)],  // v6
      field[idx(ix,   iy+1, iz+1, N)]   // v7
    ];

    // Build 8-bit cube index
    let ci = 0;
    for (let i = 0; i < 8; i++)
      if (v[i] >= threshold) ci |= (1 << i);

    if (edgeTable[ci] === 0) continue; // no crossing

    // Corner world positions
    const s = 1 / (N - 1); // grid step
    const p = [
      [ix*s, iy*s, iz*s], [(ix+1)*s, iy*s, iz*s],
      [(ix+1)*s, (iy+1)*s, iz*s], [ix*s, (iy+1)*s, iz*s],
      [ix*s, iy*s, (iz+1)*s], [(ix+1)*s, iy*s, (iz+1)*s],
      [(ix+1)*s, (iy+1)*s, (iz+1)*s], [ix*s, (iy+1)*s, (iz+1)*s]
    ];

    // Interpolate vertex on each active edge
    const ev = new Array(12);
    const em = edgeTable[ci];
    const lerp = (a, b) => {
      const t = (threshold - v[a]) / (v[b] - v[a]);
      return p[a].map((val, k) => val + t * (p[b][k] - val));
    };
    if (em & 1)    ev[0]  = lerp(0,1);
    if (em & 2)    ev[1]  = lerp(1,2);
    if (em & 4)    ev[2]  = lerp(2,3);
    if (em & 8)    ev[3]  = lerp(3,0);
    if (em & 16)   ev[4]  = lerp(4,5);
    if (em & 32)   ev[5]  = lerp(5,6);
    if (em & 64)   ev[6]  = lerp(6,7);
    if (em & 128)  ev[7]  = lerp(7,4);
    if (em & 256)  ev[8]  = lerp(0,4);
    if (em & 512)  ev[9]  = lerp(1,5);
    if (em & 1024) ev[10] = lerp(2,6);
    if (em & 2048) ev[11] = lerp(3,7);

    // Emit triangles
    const tri = triTable[ci];
    for (let i = 0; tri[i] !== -1; i += 3) {
      verts.push(...ev[tri[i]], ...ev[tri[i+1]], ...ev[tri[i+2]]);
    }
  }
  return new Float32Array(verts);
}