Marching Squares

Level Sets

A level set of a 3D surface is given by the set of points

\[\{(x,y)\quad |\quad f(x,y) = c\}.\]

Informally, this corresponds to slicing the surface with a horizontal plane at height \(c\) and graphing the curve that the intersection traces out onto the 2D plane. They are essentially a 2D projection of a 3D surface where the height information is encoded in the curve itself (i.e. distinct curves have distinct heights).

Figure 1 shows the level sets of the function \(f(x,y) = \sin(\pi x)\sin(\pi y)\).

Figure 1: Level sets of f(x,y) = sin(πx)sin(πy) for several values of c, colored by level value. Each curve is the set of points where the function equals the labeled constant.

This is what is generated by matplotlib’s plt.contour function. But you might think that plotting the level set requires solving non-linear equations like

\[\sin(\pi x)\sin(\pi y)=1/2.\]

It turns out there is a simple and elegant algorithm for drawing level sets without explicitly solving these equations called marching squares.

Paraboloid Level Sets

We work with \(f(x,y) = x^2/4 + y^2\) for the rest of this post because its level sets have a closed form we can compare against. The level sets take the form

\[\{(x,y) | x^2/4 + y^2=c\}\]

which are horizontal ellipses with the major axis twice as long as the minor axis as shown in Figure 2.

Figure 2: Level sets of f(x,y) = x²/4 + y² for several values of c.

We can see that the level sets are ellipses by dividing through by \(c\) to get

\[\left(\frac{x}{2\sqrt{c}}\right)^2+\left(\frac{y}{\sqrt{c}}\right)^2=1.\]

This is an ellipse with semi-major axis length \(2\sqrt{c}\) and semi-minor axis length \(\sqrt{c}\).¹

For a given level set \(c\), the marching squares algorithm begins by evaluating \(f(x,y)\) on a grid and determining whether each point is above or below \(c\).² Figure 3 shows a 4 by 4 grid for the level set \(x^2/4+y^2=1\).

Figure 3: A 4×4 grid covering the domain. Blue corners lie inside the c = 1 ellipse (f(x,y) < 1); grey corners lie outside. The dashed curve is the exact level set.

The Python code below shows how to implement this subroutine.

import numpy as np

def f(x: float, y: float) -> float:
    return x**2 / 4 + y**2

def label_corners(xs: np.ndarray, ys: np.ndarray, c: float) -> np.ndarray:
    """Return a boolean grid where True means the corner is inside (f < c)."""
    return np.array([[f(x, y) < c for x in xs] for y in ys])

xs = np.linspace(-3, 3, 5)   # 5 points → 4×4 grid of cells
ys = np.linspace(-2, 2, 5)

inside = label_corners(xs, ys, c=1)

The next part of the algorithm is to examine each cell and determine if there are any edges where the level set crosses. This occurs precisely when \(f(x,y)\) goes from being above \(c\) to below \(c\) along an edge of the cell (i.e. when an edge connects a blue and gray circle). To determine where the level set crosses, we can use a simple linear interpolation. Supposing \(f(x_a, y_a) > c\) and \(f(x_b, y_b) \le c\),

\[\begin{align*} x_{cross} &= \frac{c-f_b}{f_a - f_b} x_a + \left(1-\frac{c - f_b}{f_a - f_b}\right) x_b\\ y_{cross} &= \frac{c-f_b}{f_a - f_b} y_a + \left(1-\frac{c - f_b}{f_a - f_b}\right) y_b \end{align*}\]

This is a weighted average of the \(x\) and \(y\)-coordinates, weighted by how close \(c\) is to the corresponding values of \(f\) at each corner. We can easily enumerate the crossings for a given cell with the following Python code

from typing import TypeAlias

Point: TypeAlias = tuple[float, float]

def find_crossings(corners: list[Point], fvals: list[float], c: float) -> list[Point]:
    """Return the edge crossing points for one cell of the marching squares grid."""
    crossings = []
    for k in range(4):
        nk = (k + 1) % 4
        fa, fb = fvals[k], fvals[nk]
        if (fa > c) != (fb > c):            # edge straddles the level set
            xa, ya = corners[k]
            xb, yb = corners[nk]
            w = (c - fb) / (fa - fb)        # weight for corner a; (1-w) for corner b
            crossings.append((w * xa + (1 - w) * xb, w * ya + (1 - w) * yb))
    return crossings

If the cell has 2 crossings, then we connect the two points with a line and this becomes our approximation of the level set’s trajectory through the current cell. Figure 4 illustrates approximating the ellipse level sets using a coarse 4 by 4 grid of cells.

4×4 grid of cells — click a cell

Cell detail

Figure 4: Click any cell in the left panel to inspect it. Blue dots are inside the ellipse (f < 1), grey dots are outside. Diamonds mark where the contour crosses each edge; the w value gives the weight assigned to that corner in the linear interpolation.

Notice how each cell has either zero or one line segment. But what about more complicated level sets? What other configurations could occur?

The Saddle Ambiguity

Since each corner can have a value either above or below the level set value (i.e. gray or blue), and there are 4 corners, there are \(2^4=16\) different configurations. However, because of symmetry, there are only 5 distinct cases that need to be enumerated which are shown in Figure 5.

Figure 5: The five crossing cases for a marching squares cell, with one representative shown for each. The 2⁴ = 16 total corner configurations group into these 5 cases: all corners the same (0 crossings), one inside corner (2 crossings), two adjacent inside corners (2 crossings), three inside corners (2 crossings), and two diagonally opposite inside corners — the saddle case (4 crossings).

The first four of these cases are unambiguous because they either require no line segment or there is exactly one line segment that can be drawn based on the two crossing points. However, the fifth case is genuinely ambiguous. There are four crossings and two different ways to pair the points to create line segments.

We can resolve this ambiguity by computing \(f(x,y)\) at the center of the cell. Figure 6 shows how each value determines the correct pairing.

Figure 6: Resolving the saddle ambiguity by sampling f at the cell center. Left: with two diagonally opposite inside corners there are four crossings and two geometrically valid pairings. Middle: if f at the cell center is below c, the two blue corners form one connected interior region and the level set arcs around each grey corner. Right: if f at the cell center is above c, each blue corner is an isolated island and the level set arcs around each separately.

Figure 7 shows an ellipse rotated by 45 degrees whose \(c = 1\) level set requires resolving the saddle ambiguity on a 4×4 grid.

4×4 grid (tilted 45°) — click a cell

Cell detail

Figure 7: The grid is for the tilted ellipse f(x,y) = 5x²/8 − 3xy/4 + 5y²/8 = 1 (the c=1 level set of a 45°-rotated ellipse). The center cell is a saddle case: all four edges are crossed. The star marks where f is sampled to resolve the ambiguity; because f(0,0) < 0 the center is inside, so the two inside corners are connected and the level set arcs around each outside corner separately.

The Python code below implements the full algorithm

Python Code

Segment: TypeAlias = tuple[Point, Point]

def resolve_saddle(
    corners: list[Point], fvals: list[float], crossings: list[Point], c: float
) -> list[Segment]:
    """Resolve the 4-crossing saddle ambiguity by sampling the cell center."""
    cx = sum(p[0] for p in corners) / 4
    cy = sum(p[1] for p in corners) / 4
    center_inside = f(cx, cy) < c
    if (fvals[0] < c) == center_inside:
        # BL and center share status: arc around the two corners of opposite status
        return [(crossings[0], crossings[1]), (crossings[2], crossings[3])]
    else:
        return [(crossings[0], crossings[3]), (crossings[1], crossings[2])]

def marching_squares(xs: np.ndarray, ys: np.ndarray, c: float) -> list[Segment]:
    """Return all contour segments approximating the level set f(x, y) = c."""
    segments = []
    for j in range(len(ys) - 1):
        for i in range(len(xs) - 1):
            corners = [
                (xs[i],     ys[j]    ),  # BL
                (xs[i + 1], ys[j]    ),  # BR
                (xs[i + 1], ys[j + 1]),  # TR
                (xs[i],     ys[j + 1]),  # TL
            ]
            fvals = [f(x, y) for x, y in corners]
            crossings = find_crossings(corners, fvals, c)
            if len(crossings) == 2:
                segments.append((crossings[0], crossings[1]))
            elif len(crossings) == 4:
                segments.extend(resolve_saddle(corners, fvals, crossings, c))
    return segments

xs = np.linspace(-3, 3, 50)
ys = np.linspace(-2, 2, 50)
segments = marching_squares(xs, ys, c=1)

Conclusion

Marching squares answers the question “how do you draw the curve defined by \(f(x, y) = c\) without solving that equation?” Surprisingly, the answer is to never solve it at all. Marching squares proceeds by

• classifying each grid corner as inside or outside the level set
• locating where the boundary crosses each edge by linear interpolation
• connecting the crossings.

Fourteen of the sixteen possible corner configurations are unambiguous. The one exception is the saddle case, where two diagonally opposite corners are inside and four edge crossings admit two geometrically valid pairings. A single extra evaluation of \(f\) at the cell center resolves the ambiguity: whichever pairing connects corners that share interior status with the center is correct.

The curve defined by \(f(x, y) = c\) can be drawn to arbitrary precision using nothing more sophisticated than linear interpolation!

1. The \(c = 1\) ellipse can be parameterized as \(x(\theta) = 2\cos(\theta)\), \(y(\theta) = \sin(\theta)\) for \(\theta \in [0, 2\pi)\); more generally \(x(\theta) = 2\sqrt{c}\cos(\theta)\), \(y(\theta) = \sqrt{c}\sin(\theta)\). Most level sets do not admit such a closed form, which is why marching squares works from function evaluations alone. ↩

2. Checking whether \(f(x,y) > c\) or \(f(x,y) < c\) is equivalent to checking the sign of \(g(x,y) = f(x,y) - c\), so without loss of generality the algorithm tracks the zero level set of \(g\). ↩