1. Introduction

The Marching Squares algorithm is a computer graphics algorithm introduced in the 1980s that can be used for contouring. We can use the Marching Squares algorithm to draw the lines of constant:

  • altitude on a topographical map
  • temperature on a temperature heat map
  • pressure on a contour map for a pressure field

In this tutorial, we’ll learn how the Marching Squares algorithm determines the contour from a grid of sample points of the image.

2. Contour Curves and Level Curves

Let’s spend a brief moment understanding the mathematically precise definition of contour lines.

2.1. Mathematical Definition

We have been drawing graphs of functions of one-variable such as f(x)=\sin x or implicit functions like x^2 +y^2 = 1 for so long, that we give the matter little thought. Plotting functions of two variables is a much more difficult task.

The trick to putting together a reasonable graph of a function of two or more variables is to find a way to cut down on the dimensions involved. One way to achieve this is to draw certain special curves that lie on the surface z=f(x,y). These special curves, called contour curves, are obtained by intersecting the surface with horizontal planes z = c for various values of the constant c.

These curves in the xy-plane are called the level curves of the original function f.

The level curve at height c of any function f is the curve in \mathbf{R}^2 defined by the equation f(x,y)=c, where c is a constant. In mathematical notation,

    \[\text{Level curve at height c} = \{(x,y) \in \mathbf{R}^2| f(x,y) = c\}\]

2.2. Examples

Consider the function f:\mathbf{R}^2 \to \mathbf{R} defined as:

    \[f(x,y)=4 - x^2 - y^2\]

By definition, the level curve at height c is:

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

Thus, we see that the level curves for c < 4 are circles centered at the origin of radius \sqrt{4 - c}. They are summarized in the table below:

c Level\ curve\ x^2 + y^2 = 4 - c
-21 x^2 + y^2 = 25
-12 x^2 + y^2 = 16
-5 x^2 + y^2 = 9
0 x^2 + y^2 = 4
3 x^2 + y^2 = 1
4 x^+y^2 = 0 \implies x = 0, y = 0

The family of level curves, the “topographic map” of the surface z = 4- x^2-y^2 is shown in the figure below:

2dfuncplot 02

If we stack all level curves together, we can get a feeling for the complete graph of z=4-x^2 - y^2. It’s a surface that looks like an inverted dish and is called a paraboloid:

2dfuncplot 01

For our discussion, we require the output contour to possess the following properties:

  • It must separate the interior points from the exterior points.
  • Thus, it must be an orientable and closed curve.

3. Essentials of the Algorithm

Consider the simplistic image of a cloud given below. We are interested to trace the contours of the cloud:


For simplicity, let’s assume that the image of a cloud is represented by a square matrix A of 100 \times 100 pixels. Each pixel value a_{ij} represents a single sample of light, carrying brightness information. It ranges between 0 and 255. Black is represented by 0 and white by 255.

3.1. Sampling the Image

We divide the entire image domain D into 4 \text{ rows }\times 4\text{ columns } = 16 squares. Each square has dimensions by 5 \times 5. We sample the image at each of the 5^2 = 25 grid points:

diagram 20220522

Next, we iterate through the grid G to determine the vertex state in comparison to an iso-value v. Each grid vertex is now assigned a binary state (0 or 1). An iso-value \sigma serves as a threshold for determining the vertex state. Suppose we choose \sigma = 200 as the threshold. The vertex state is determined as follows:

  • ON(1)– Pixel value less than the iso-value \sigma
  • OFF(0) – Pixel value greater than the iso-value \sigma

The resulting binary grid F, thus obtained is:

diagram 20220522-1

3.2. Configuration of Each Sub-Grid

Next, a 5 \times 5 subgrid can be associated with a configuration. We iterate through G. At each 5 \times 5 sub-grid, we form a binary code based on the vertex value by iterating clockwise from the top-left corner to the bottom-left corner.

For example, if a given subgrid has four corners:

diagram 20220522-3

its binary code is said to be (0110)_2. The equivalent of the binary code in the decimal number system is called the configuration of the sub-grid. In the above example, the sub-grid has the configuration 6 Thus, we have a configuration associated with each sub-grid as follows:

diagram 20220522-2

3.3. Lookup Table for the Contours

Imagine a single square with the contours of an object cutting through this square. That is, the square does not lie wholly in the interior or exterior of the object. A square has 4 edges and for each edge, we have 2 possible choices – either the contour(boundary) intersects this edge, or it does not. This results in 2^4 = 16 possible configurations.

The configuration of each sub-grid is matched with one entry in the contours lookup table below:

diagram 20220522-4

The contours lookup table is intuitive. As mentioned before, there are 2^4 = 16 possible configurations.

Let’s consider the configuration 7 = (0111)_2:

diagram 20220522-5

We think of the purple region as the interior of an object, and the blue region as the exterior of the object. We denote the ON(1) state by a black vertex and the OFF(0) state by a white vertex. We define a bipolar pair to be two vertices in opposite states.

Intuitively, the contour of an object must intersect the edge of a bipolar pair. Thus, the configuration 7 = (0111)_2 must correspond to an object whose contour intersects the top and left edges of the square. In this fashion, the contours in each sub-grid are approximated, as shown in the figure below:

diagram 20220522-6

4. Algorithm

The Marching Squares algorithm can be divided into two steps:

  1. Sampling the grid points G and converting them to a binary matrix F.
  2. March the squares in F and build a set of edges \gamma.

Let’s see the pseudocode:

Rendered by QuickLaTeX.com

4.1. Sampling the Grid

The input to the Marching Squares algorithm is some sample of the original matrix G. The parameters \text{stepX} and \text{stepY} control the resolution of the sample. For example, if \text{stepX}=10 and \text{stepY}=15, we select the sequence of pixel-values:

    \[\begin{bmatrix} g_{0,0} & g_{0,10} & g_{0,20} & g_{0,30} & \ldots \\ g_{15,0} & g_{15,10} & g_{15,20} & g_{15,30} & \ldots\\ g_{30,0} & g_{30,10} & g_{30,20} & g_{30,30} & \ldots\\ \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix}\]

The result of this operation should yield a matrix of size p=\frac{m}{\text{stepX}} rows and q=\frac{m}{\text{stepY}} columns.

Each selected pixel value is compared with the iso-value \sigma and converted to a bit-value 0 or 1. F is used to hold the bit matrix. The actual (x,y)-coordinates of each pixel in the sample F, concerning the original matrix G is stored as a named tuple in a matrix \text{coord}.

For example, if G[15][20] = 100, \sigma = 200, we would find F[1][2]=1, \text{coord}[1][2]=(x=15,y=20).

Let’s see the pseudocode now:

Rendered by QuickLaTeX.com

4.2. Marching the Squares In \mathbf{F}

Next, we march the squares in the binary grid F. Each square in F, beginning with its top-left corner, has vertices, in clockwise order : F[i,j], F[i,j+1], F[i+1,j+1] and F[i+1,j]. We compute the reference value \kappa of each square. We also compute the mid-points of each square edge: north, east, south, and west. Finally, we determine the contour S by searching the reference value in the contours lookup table:

Rendered by QuickLaTeX.com

Here’s a topographical map of the Mont Blanc region in the Swiss Alps. The Marching Squares algorithm when executed on this image produces nice contours around the mountain peaks:

marching squares

5. Conclusion

In this article, we learned about an interesting Computer Graphics algorithm called the Marching Squares. To summarize, it is based on two core ideas. First, the contours can be constructed piecewise within each grid square, without reference to other squares. Second, each isocontour segment in a grid square can be retrieved from a lookup table.