**The marching squares algorithm**

I recently had to code the marching squares algorithm for an assignment involving implicit surfaces.

The idea behind this algorithm is simple. You have a grid, with a value assigned to each intersection. In the case of implicit surfaces, those values represent the result of evaluating the implicit function at that point in space.

For example, the implicit function for a full unit circle is

So if we were to evaluate it at each point in a grid, we’d get something like this.

Forget about the circle right now. We’re only interested in the values themselves. The point at the center is (0,0), so the value there is 1.00 because

The marching squares algorithm needs that we classify each square in the grid as belonging to one of sixteen classes, depending on the function values in their corners. These values, as far as the algorithm is concerned, are binary. In this case, we will consider negative values to be 0 and positive ones to be 1. The class of your square is determined by the binary sequence that describes its vertices, read top-left top-right, bottom-right, bottom left. The square below, for instance, can be represented by the sequence 0011 (negative, negative, positive, positive).

This sequence means that the line filling this square must start at the edge on the left and finish at the edge on the right. This is because the ends of the line we draw must start and end on edges whose vertices are of a different sign.

So, for every square we have 4 values that correspond to the 4 vertices, and depending on whether they’re negative or positive we must draw a different thing. The naive approach is as simple as it is horrid. It would be close to this.

if(v0 <; 0 && v1 <; 0 && v2 <; 0 && v3 <; 0){ drawNothing(); } else if(v0 <; 0 && v1 <; 0 && v2 <; 0 && v3 > 0){ drawLeftToBottom(); } else if(v0 <; 0 && v1 <; 0 && v2 > 0 && v3 <; 0){ drawRightToBottom(); } //...

And so on for 16 painful condition checks. This is terrible, of course, and it can be avoided.

**The function pointer approach**

First of all, each sequence can be understood as a binary number, and we can calculate its decimal value very easily.

double decimalValue = v0*8 + v1*4 + v2*2 + v3;

Like that, we could replace the flood of if statements with a simple switch block. That, however, is not the approach we’re going for.

Taking advantage of the distinct numerical representation of each sequence, we can put together and array of pointers to functions.

For doing this, C++ requires of us to define the type contained in our array. We want to pass to each function the points of the current square and the implicit function, so that it can calculate the exact locations where the line will begin and end. As you can see, I’m defining this in a class called Grid.

typedef void(Grid::*FunctionPointer)(Point2, Point2, Point2, Point2, Implicit &);

Afterwards, we can declare the array itself.

FunctionPointer drawMSCase[16];

And initialize it, mapping each draw function to a position in the array (the decimal value of its sequence). Notice how the case we saw earlier, 0011, which is 3 in binary, is mapped to drawEdgeH, that is, a function that draws a horizontal line. That is the method that draws a line between the two vertical edges.

drawMSCase[0] = &Grid::drawNothing; drawMSCase[1] = &Grid::drawEdgeLT; drawMSCase[2] = &Grid::drawEdgeTR; drawMSCase[3] = &Grid::drawEdgeH; drawMSCase[4] = &Grid::drawEdgeRB; drawMSCase[5] = &Grid::drawEdgeBRLT; drawMSCase[6] = &Grid::drawEdgeV; drawMSCase[7] = &Grid::drawEdgeBL; drawMSCase[8] = &Grid::drawEdgeBL; drawMSCase[9] = &Grid::drawEdgeV; drawMSCase[10] = &Grid::drawEdgeBLTR; drawMSCase[11] = &Grid::drawEdgeRB; drawMSCase[12] = &Grid::drawEdgeH; drawMSCase[13] = &Grid::drawEdgeTR; drawMSCase[14] = &Grid::drawEdgeLT; drawMSCase[15] = &Grid::drawNothing;

After doing this, our code for identifying the class of each square and drawing the corresponding option will look like this.

int a0 = v0 < 0 ? 0 : 1; int a1 = v1 < 0 ? 0 : 1; int a2 = v2 < 0 ? 0 : 1; int a3 = v3 < 0 ? 0 : 1; int index = a0*8 + a1*4 + a2*2 + a3; (this->*drawMSCase[index])(p0, p1, p2, p3, f) ;

In fact, we should write a function that gives us the decimal value. The code would then end up looking like this

int index = getDecimalValue(a0, a1, a2, a3); (this->*drawMSCase[index])(p0, p1, p2, p3, f) ;

Much cleaner than the 16 if statements we started with, indeed.

In this case, this approach does not provide us with much more than cleaner code, but for other situations, where our data structures do not have a fixed size, mapping functions to lists can be of huge advantage.

**Switching to lambdas**

Function pointers, however, are a little old fashioned. A more modern approach would be to use a map of lambda expressions.

First, we create a structure that maps integers to functions.

std::map<int, std::function<void(Point2, Point2, Point2, Point2, Implicit &)> > drawMSCase;

We then insert each one of the functions.

drawMSCase[0] = [] (Point2, Point2, Point2, Point2, Implicit &) { //Code that does the drawing for case 0 here! }; //.. drawMSCase[15] = [] (Point2, Point2, Point2, Point2, Implicit &) { //Code that does the drawing for case 15 here! };

And now we can call the function doing this.

drawMSCase[index](p0, p1, p2, p3, f);

The code for the function pointer version (a Qt Creator project) can be found here.