Marching Squares

A contour-tracing algorithm that extracts iso-lines from a 2D scalar field by classifying each grid cell into one of 16 cases and looking up the corresponding edge segments.

The 16 cases

A scalar field assigns a numeric value to every grid vertex. Choose a threshold and each vertex is either inside (above threshold — bright) or outside (below — dark). A cell has four corners so there are 2⁴ = 16 possible configurations, indexed by a 4-bit number: TL = bit 3 (8), TR = bit 2 (4), BR = bit 1 (2), BL = bit 0 (1).

Each case maps to zero, one, or two line segments. Cases 5 and 10 are saddle ambiguities — the correct topology is genuinely ambiguous and either of the two valid interpretations can be chosen.

Classify a cell

Click any corner to toggle it inside or outside the threshold. The case index, bit pattern, and contour segments update instantly.

Click a corner to toggle it.

Tracing a scalar field

A 26 × 14 grid of values computed from a sum of Gaussian blobs. Each vertex is classified by the threshold and marching squares runs on every cell. Drag the threshold to watch the contour move through the field.

Threshold

38%

Inside vertex

Outside vertex

Contour (interpolated)

Each contour vertex is placed by linear interpolation along the edge — not simply at the midpoint. The crossing position is t = (threshold − v₀) / (v₁ − v₀), producing smooth contours even on coarse grids.

Midpoint vs linear interpolation

The simplest implementation places each contour vertex at the exact midpoint of the crossing edge, ignoring the actual field values. Linear interpolation uses t = (threshold − v₀) / (v₁ − v₀) to find exactly where the field crosses. On a fine grid the difference is small; on a coarse grid the midpoint contour looks like a staircase of 45° steps.

Threshold

38%

The code

The complete algorithm in three pieces. The surprising thing is how little code it takes — the domain knowledge lives entirely in the lookup table.

// Bit encoding: TL=8  TR=4  BR=2  BL=1
// Edge indices:  0=left  1=top  2=right  3=bottom
const CELL_EDGES = [
  [],             // 0  — all outside, no crossing
  [[0,3]],        // 1  — BL:       left → bottom
  [[3,2]],        // 2  — BR:       bottom → right
  [[0,2]],        // 3  — BL+BR:   left → right
  [[1,2]],        // 4  — TR:       top → right
  [[1,2],[0,3]], // 5  — TR+BL:   top→right  +  left→bottom  ← saddle
  [[1,3]],        // 6  — TR+BR:   top → bottom
  [[0,1]],        // 7  — TR+BR+BL: left → top
  [[0,1]],        // 8  — TL:       left → top
  [[1,3]],        // 9  — TL+BL:   top → bottom
  [[0,1],[3,2]], // 10 — TL+BR:   left→top  +  bottom→right  ← saddle
  [[1,2]],        // 11 — TL+BR+BL: top → right
  [[0,2]],        // 12 — TL+TR:   left → right
  [[3,2]],        // 13 — TL+TR+BL: bottom → right
  [[0,3]],        // 14 — TL+TR+BR: left → bottom
  [],             // 15 — all inside, no crossing
];

Cases 5 and 10 produce two segments each — the only two configurations where the contour is genuinely ambiguous. In case 5 (TR + BL active), two diagonally opposite corners are inside. The contour could either connect top→right and left→bottom (treating the blobs as separate) or top→left and bottom→right (treating them as merged). Both are valid. This implementation picks the first interpretation consistently.

Cases 0 and 15 (all outside / all inside) produce no segments. They are their own trivial complements.

Complement symmetry: case N and case (15 − N) always produce the same set of segments — swapping inside and outside doesn't change where the boundary runs, only which side is which.

// Classify each cell's four corners, then look up segments
for (let cy = 0; cy < ROWS; cy++)
  for (let cx = 0; cx < COLS; cx++) {
    const tl = field(cx,   cy  ) > threshold ? 8 : 0;
    const tr = field(cx+1, cy  ) > threshold ? 4 : 0;
    const br = field(cx+1, cy+1) > threshold ? 2 : 0;
    const bl = field(cx,   cy+1) > threshold ? 1 : 0;

    const caseIndex = tl | tr | br | bl;  // bitwise OR builds the index

    CELL_EDGES[caseIndex].forEach(([a, b]) => drawSegment(a, b, cx, cy));
  }

Bitwise OR as table index. Assigning each corner a distinct power of 2 (8, 4, 2, 1) means the bitwise OR of the four classification results is exactly the case index — no branching, no conditionals. Four comparisons produce the table row directly. This is the key insight that keeps the inner loop tiny and branch-free.

Each cell is fully independent of its neighbours. The algorithm visits every cell exactly once — O(rows × cols) — and each cell only reads four vertex values. This makes it trivially parallelisable; a GPU compute shader assigns one thread per cell and runs the entire grid in a single pass.

// Find where the threshold crosses edge e of cell (cx, cy)
// using linear interpolation between the two endpoint values
function edgeCrossing(e, cx, cy, threshold) {
  // Map edge index to the two grid vertices it connects
  const [gx0,gy0, gx1,gy1] = [
    [cx,   cy,   cx,   cy+1],  // 0: left  — TL→BL
    [cx,   cy,   cx+1, cy  ],  // 1: top   — TL→TR
    [cx+1, cy,   cx+1, cy+1],  // 2: right — TR→BR
    [cx,   cy+1, cx+1, cy+1],  // 3: bottom — BL→BR
  ][e];

  const v0 = field(gx0, gy0);
  const v1 = field(gx1, gy1);
  const t  = (threshold - v0) / (v1 - v0);  // fraction along edge [0, 1]

  return [gx0 + t * (gx1 - gx0),
          gy0 + t * (gy1 - gy0)];
}

Why interpolate? Using only the midpoint places every contour vertex at the exact centre of each edge, regardless of how close the field value is to the threshold. Linear interpolation — t = (threshold − v₀) / (v₁ − v₀) — positions the vertex precisely where the field crosses the threshold, assuming the field varies linearly between vertices. On a fine grid the difference is small; on a coarse grid (like the demo above) it prevents the contour from looking like a staircase of 45° segments.

Edge case: if v₀ === v₁ the denominator is zero (the edge runs exactly along the iso-surface). The result is clamped to t = 0.5.