Previous : Visualising Scalar Fields : Introduction

The Marching Squares algorithm is derived from the Marching Cubes algorithm, which was used to visualise data from medical scans.

I’m going to use Marching Squares to create 2D isosurfaces from the data in a scalar field, starting with a simple demo that creates a contour map, then a couple of other well known use cases : Metaballs and Tile Mapping.

*I’m generating isosurfaces, not isolines, because it’s easier in Unity. If you’re interested in drawing isolines, I’ve added in some notes at the relevant stages to show you how to do this.*

#### Making a contour map

*You’ll find this example in the ContourMap scene in the demo project (found here).*

To generate the contour maps, we’re going to follow these steps:

- Generate scalar data and build the grid
**Iterate through the grid and draw the isosurfaces****Find the right configuration for each square****Smooth the edges****Draw the contour meshes**

The bold parts (steps 2-4) encompass the marching squares algorithm.

##### Generating Scalar data and Building the grid

To march some squares, we first need the squares. Typically, the grid would usually be determined by the scalar field, with the sample points sitting in the corners of the squares, as shown below.

But we don’t have a scalar field yet…We could arbitrarily make one up and manually input a load of values, but it’s probably easier to switch things around a little and create the grid first, then dynamically determine the values for the scalar field based on the size of the grid we choose. If you already have a scalar field, skip the part where I make up the values for the grid and just use the ones you have.

So, first up we’ll create the grid. We need to decide on the resolution (number of squares) of our grid, which will determine how accurately the isosurfaces/isolines are drawn. I’ve started with 50, as it’s a nice round number, but when we get to the drawing part I’ll show some results with different resolutions so you can see the impact. I’m using a 16:9 screen in Unity, so will get a grid size of around 100 x 57 squares.

Next up, we want to work out our scalar field. Let’s imagine that we have a cone sitting in the middle of the screen. Its base diameter is the height of the screen, and it stretches out towards us with a height equal to half the height of the screen. We’re going to make a height map for this cone, where each point in our grid is the height to the surface at each point (x,y).

The distance of the sample point to the centre of the cone can be determined using pythagoras’ theorem:

`(point.x - coneCentre.x)<sup style="top: -0.5em;">2</sup> + (point.y - coneCentre.y)<sup style="top: -0.5em;">2</sup> = distanceFromCentre<sup style="top: -0.5em;">2</sup>`

Now we work out how far from the peak of the cone the point is.

`distanceFromPeak = distanceFromCentre / angleOfSlope`

We know that the angle of the slope of the curve is 45° because the base radius is equal to the height, so

`distanceFromPeak = distanceFromCentre / tan(45°)`

`distanceFromPeak = distanceFromCentre / 1`

`distanceFromPeak = distanceFromCentre`

Now let’s use this to get the height:

`height = heightOfCone - distanceFromPeak`

#### Finding the right configuration for each square

Now we have our grid, we can create the contour isosurfaces. We are going to create four surfaces, with threshold values of 0, 1, 2, 3 and 4 respectively.

To do this, we look at each square and inspect the corner values. There are two different states, inside our threshold value and outside our threshold value, and four different corners to consider, so we have 16 (2^{4}) possible cases:

A black dot in the corner indicates that the value in that corner is higher than the threshold value and should be inside our line. A white dot indicates that the value is lower than the threshold and outside the line.

We’re going to create a 4-bit case id (0000), with one bit assigned to each corner of the square (going from right to left of the id, we assign the bottom left corner, then bottom right, top right, then top left). We use bitwise OR to set the bit associated with a corner when the corner’s value is over the threshold value.

Values used to set the bits :

0001 : Bottom left corner bit (binary to decimal = 1)

0010 : Bottom right corner bit (binary to decimal = 2)

0100 : Top right corner bit (binary to decimal = 4)

1000 : Top left corner bit (binary to decimal = 8)

Let’s say both the top left and bottom right corners are over the threshold value. Using bitwise OR ( | ) to set the bits:

We start our id equal to 0, with no corners evaluated:

`id = 0000`

We set the bottom right corner bit by ORing it with 0010

`id |= 0010 `

`And we get an `

id equal to 0010 (or 2 in decimal)

Then set the top left corner bit by ORing it with 1000 (8 in decimal).

`id |= 1000`

And we get an id equal to 1010 (or 10 in decimal)

In C# :

```
int GetConfigurationForSquare (GridSquare square)
{
<span style="visibility: hidden;">++++</span>int caseId = 0;
<span style="visibility: hidden;">++++</span>if (square.bottomLeftValue >= threshold)
<span style="visibility: hidden;">++++</span>{
<span style="visibility: hidden;">++++++++</span>caseId |= 1;
<span style="visibility: hidden;">++++</span>}
<span style="visibility: hidden;">++++</span>if (square.bottomRightValue >= threshold)
<span style="visibility: hidden;">++++</span>{
<span style="visibility: hidden;">++++++++</span>caseId |= 2;
<span style="visibility: hidden;">++++</span>}
<span style="visibility: hidden;">++++</span>if (square.topRightValue >= threshold)
<span style="visibility: hidden;">++++</span>{
<span style="visibility: hidden;">++++++++</span>caseId |= 4;
<span style="visibility: hidden;">++++</span>}
<span style="visibility: hidden;">++++</span>if (square.topLeftValue >= threshold)
<span style="visibility: hidden;">++++</span>{
<span style="visibility: hidden;">++++++++</span>caseId |= 8;
<span style="visibility: hidden;">++++</span>}
<span style="visibility: hidden;">++++</span>return caseId;
}
```

#### Ambiguous Cases

There are two cases (Case 5 and Case 10) in the case diagram above that can be drawn in two different ways: two triangles in the corners (black corners and red lines) or a saddle running through the centre (connecting the two black corners with black lines).

The images below show the different ways to draw case 10 (blue is over or equal to the threshold):

Both are valid, so how do we know which to use?

One way to decide is to get a extra value for the centre of the square that takes the average of the square’s corners : if that point is above the threshold value, then you use the saddle configuration (righthand image above), otherwise you use the configuration that keeps the two surfaces separate (lefthand image).

The easy way (and the method I’ve used here) is just to pick one and be done with it – I’ve always defaulted to the configuration that separates the isosurfaces.

We won’t see this issue in this current tutorial as it only presents itself when there is more than one peak in the scalar field (and we’re only showing one cone). The Metaballs post (after this one), does have potential to exhibit this behaviour and I’ll show the results of using both configurations there so you can see the difference it makes.

#### Drawing the cases

Once you have a case that matches the corner configuration of the square, you know what to draw. From the case diagram, if you’re drawing isolines, you draw the lines shown in that case. If you’re drawing an isosurface, you draw the part of the case encompassed by the line and the black corners.

In the image below, you can see the different cases drawn onto a mesh (cyan lines are debug lines so you can see the squares properly), from case 0 in the top left square to case 15 at the bottom, and the triangles needed to draw each case.

There are eight possible vertices on each square that are used for drawing these triangles : 0,0 is bottom left, and 2,2 is top right and midpoints are at 1 (so the mid point on the left edge is 0,1).

You can see case 7 below is made up of three triangles:

To draw this in Unity, we create a mesh, then pass that mesh a list of the vertices (points) and triangles that create our shape.

First, our vertices :

```
Vector2[] vertices = new Vector2[5];
vertices[0] = new Vector2(0, 0);
vertices[1] = new Vector2(0, 1);
vertices[2] = new Vector2(1, 2);
vertices[3] = new Vector2(2, 2);
vertices[4] = new Vector2(2, 0);
```

And now our triangles. Each triangle is a list of three entries – the indices of the points in the triangle – as entered in the vertices array above.

```
int[] triangles = new int[]
{
<span style="visibility: hidden;">++++</span>0, 1, 2, // points (0,0),(0,1),(1,2)
<span style="visibility: hidden;">++++</span>2, 3, 0, // points (1,2),(2,2),(0,0)
<span style="visibility: hidden;">++++</span>0, 3, 4 // points (0,0),(2,2),(2,0)
}
```

You need to make sure that when you define your triangles, you list the vertices in a clockwise order, so the shape is facing the camera. Defining them in an anticlockwise order will mean the shape will face backwards (you won’t be able to see the mesh from the front).

You can see specific vertices and triangles in the MarchingSquaresModel class in the demo project.

#### Results

Now we’ve gone through every grid square, drawn the correct configuration and updated the mesh, we’ve done enough to display our scalar field. For ease, I’ve created four different meshes, one for each threshold value.

So what should we expect to see? If you think back to the cone again; imagine we’ve cut a slices out of the cone where the height equals the threshold values (0, 1, 2, 3 and 4 – the cone’s peak is at a height of 5), and then slapped those slices on top of one another. It should look like five circles, with the radius getting smaller with each circle.

With a resolution of 50 squares, this is what it looks like:

You can definitely make out the different isosurfaces, which is good, but they don’t follow the expected circle shape very well. A higher grid resolution should help, and cranking it up to 100 gives us this:

The isosurfaces are more defined here, but the edges are still pretty lumpy. We need to smooth them out!

#### Smoothing the edges

For our isosurface to be accurate (and therefore smooth), all points on its edge need to be equal to the threshold value. This is not the case with our visualisation so far.

Let’s take another look at the case diagram:

For each case, you can see that when a line passes between two corners (in all but cases 0 and 15) it always passes through the middle of that edge, exactly between the two corners.

The value at that point will be whatever the average value of the two corners is, not necessarily our threshold value. This is why we currently have a bumpy circle and not a smooth one: we’re only approximating the line.

You can see this a bit better on the diagram below:

Diagram A shows our current setup. We have a square that has a case 2 configuration and so we’ve drawn a line connecting edge AB to edge BD, from point S to point T. We’re looking to draw an isoline with a threshold value of 1 and have approximated this by making point S and T sit in the middle of their edges, as outlined in the case diagram. However, as you can see from the values at each corner, the point on edge AB where the value is 1 is not in the centre – the centre has a value of 2.5 ((B – A) / 2) – so the correct point must be closer to corner A. Similarly, the point on the edge BD where the value is 1 is not in the centre of BD (where T currently sits), but closer to point D.

To move S and T to the point on their edge where the value is 1, and create a more accurate isoline, we linearly interpolate between the corners.

There are four scenarios when we need to do this : when y is at an edge and x is a midpoint (i.e. point S in the image above), and when x is at an edge and y is a midpoint (i.e. point T in the image above).

Given a square with corner values A, B, C and D (see image above) and the threshold value, we can find the point between corners where the value is 1:

NB : this is for a scale of y that runs from bottom to top.

```
if (point.<i>x</i> == EDGE && point.<i>y</i> == MIDPOINT)
{
<span style="visibility: hidden;">++++</span>point.<i>y</i> = ((<i>Threshold</i> - square.A) / (square.<i>C</i> - square.<i>A</i>));
}
else if (point.<i>x</i> == EDGE && point.<i>y</i> == MIDPOINT)
{
<span style="visibility: hidden;">++++</span>point.<i>y</i> = ((<i>Threshold</i> - square.<i>B</i>) / (square.<i>D</i> - square.<i>B</i>));
}
else if (point.<i>x</i> == MIDPOINT && point.<i>y</i> == EDGE)
{
<span style="visibility: hidden;">++++</span>point.<i>x</i> = ((<i>Threshold</i> - square.<i>A</i>) / (square.<i>B</i> - square.<i>A</i>));
}
else if (point.<i>x</i> == MIDPOINT && point.<i>y</i> == EDGE)
{
<span style="visibility: hidden;">++++</span>point.<i>x</i> = ((<i>Threshold</i> - square.<i>C</i>) / (square.<i>D</i> - square.<i>C</i>));
}
```

This gives us a new line, shown in diagram B, where points S and T have a value of 1. Notice S has moved closer to corner A (A’s value is closer to our threshold of 1 than B’s), and T is now a lot closer to D.

Implementing interpolation, we get this result (with a resolution of 50) – nice and smooth!

Cranking up to 80 smooths the inside circle a little more:

And that’s all there is to it!

If you’re interested in playing around with Marching Squares some more, I’m going to follow up with a couple more demos showing some well known uses cases – Metaballs and Map Tiles.

Previous : Visualising Scalar Fields : Introduction

Next : Visualising Scalar Fields : Marching Squares Part 2 : Metaballs