NOTE: After I wrote this article I realized that the Marching Cubes is not as simple as most people think. I’m going to write another part about problems I’ve found.

I always wanted to implement my own Marching Cubes algorithm. Last days I had some time so I decided to give it a go…

I started from checking the internet resources, so I’m not reinventing the wheel, and I supposed that there’s many different implementations, in many languages and so on… What I’ve found is that there’re just few good articles about that topic, and almost all of them are based on some ‘magical’ table, that contains all the data you need to render every cube case.

I could copy that *magic* table and use it in my code but it’s obviously not the option, but even ‘creating’ those data on my own doesn’t seem to be fun. So I decided to go another way and to create all the data algorithmically. I tried not to use any static tables, however there’re still some places that I’d like to improve and simplify.

I’m not going to describe here the Marching Cubes algorithm itself, there are many articles about it on the internet, I’ll just write up my approach.

Algorithmic marching cubes

The plan is really simple – I don’t want to create a full-blown Marching Cubes library, at this stage I don’t care about performance too much, I just want to make some simple to use code, that can render isosurface. It should be also easy to understand and use in few lines. I don’t want it to be complex to configure and build and I don’t want to use any dependencies but just old, good OpenGL.

Cube and vertex ordering

We’ll compute a lot of information about our cubes geometry, so it’s really important to arrange cube vertices in a way that will simplify further computations. As we know cube has 8 corners, so we can store corner number using 3 bits. Our cube is also a 3D figure with 3 axes, so let’s assign bits (0,1,2) to axes (x,y,z) – bit value 0 means 0.0f along particular axis, bit 1 means 1.0f.

For example vertex of index (000b) is placed in (0,0,0) while vertex (110b) is (0,1,1) in 3D space.

You can see axes and vertex ordering on the image below.

In most of the code I’ve ever seen, that used/created any cubes, there were typically two loops of vertices – lower and upper, but my code uses two Zorro signs :). Thanks to this corner ordering you can very easily move between corner index, it’s coordinates, it’s neightbours, etc… You’ll see the effects in the code.

At the beginning we want to compute our corners positions so we iterate over all of them and translate their index to 3D vector

for( int v = 0; v < 8; v++ ) {
    for( int compound = 0; compound < 3; compound++ )
        vertexOffset[v].f[compound] = (v&(1<<compound)) ? 1.0f : 0.0f;
}

Edges

The next concept is an edge, that connects two corners. Thanks to our corner ordering we know, that two corners are connected by an edge, if their indexes differ by one and only one bit.

Here I iterate over all corner-corner combinations and check if they form an edge.

int edgesNum = 0;
for( int v1 = 0; v1 < 7; v1++ ){
    for( int v2 = v1+1; v2 < 8; v2++ ){
        if( _oneBitDiff(v1,v2) ) {  //one bit difference
            edgeToVertex[edgesNum][0] = v1;
            edgeToVertex[edgesNum][1] = v2;
            edgesNum++;
        }
    }
}
bool _oneBitDiff( int v1, int v2 )
{
    int diff = v1 ^ v2;    //  XOR
    if( diff == 1 || diff == 2 || diff == 4 )
        return true;
    return false;
}

These were our basic data, but we can also add few helper methods, like finding an edge by two corner indexes:

int MarchingCubes::_findEdge( int v1, int v2 )
    for( int edge = 0; edge < 12; edge++ ) {
        if( v1 == edgeToVertex[edge][0] &&
            v2 == edgeToVertex[edge][1] )
                return edge;
        if( v2 == edgeToVertex[edge][0] &&
            v1 == edgeToVertex[edge][1] )
                return edge;
    }
    return -1;
}

Conversion of 8-bit code to int[8] sign table:

void _codeToSignTable( int code, int* tab )
{
    for( int counter = 0; counter < 8; counter++ ) {
        int i = 1 << counter;
        if( code & i )
            tab[counter] = 1;
        else
            tab[counter] = -1;
    }
}

When I polish more of this code I’ll describe more utility functions.

Let’s make some geometry…

The cube has 8 corners, and each of them can be inside or outside of our surface so it gives us 256 combinations. For (almost)each of those combinations we should render some triangles. But how to know what triangles do you need?

I start by iterating over all of those 256 combinations, so each time I have just one case to analyze. Cases ‘0’ and ‘255’ means cubes that are entirelly inside/outside of the surface so they don’t need any triangles, all other cases will result in at least one triangle to render. Given the one corner combination we can check it over few basic situations. The simplest of them is a “single vertex”.

Single vertex

This is a case, when we have one of the corners on one side of the surface while corner neightbours are on the other side. This case needs just one triangle and you can see it on the image below. Corners marked by black or white dots show what value do we expect on each corner. If the corner isn’t marked by any dot then we don’t care about its value.

We iterate over each corner.
The current corner has some sign (+/-)
Get current corner neightbours
If all neightbours have different sign than current corner – we have a “single vertex” here.
If that’s the case then create one triangle between 3 edges starting from current corner.

Single edge

Slightly more complicated case, where instead of one corner we have two corners that shares an edge. In this case we need two triangles to represent it.

Half Split

This is the case, where the surface goes along one of the axes and half of our corners are inside of surface and others are outside. The code is different, because we need to iterate over axes and not the corners. For each of 3 axes we check if it splits our corners correctly and if we find such case we create 2 triangles along this axis.

At this point I could see that something is actually rendered but most of the corner combinations were still missing. So now is the time for more interesting cases. Let’s call them the “advanced ones”.

Now I could also focus a bit on my rendering code and some scalar field implementation. I created simple animation, tweaked few params, etc. so I could continue my further work.

Advanced cases

I created a bit of code that prints triangle data for all of 256 cases and I looked for the empty ones. At the beginning you have a feeling that there’s so many missing cases that you’ll never finish this work 🙂 but most of those cases are just reflection or rotation of another ones so it’s really not so hard. For each new case I’ve drawn a cube and corresponding corner signs and tried to find some triangle configuration that’s best for particular case. Then I implemented a method that tests for this particular case.

Here’s what I’ve found.

Triple vertex

This is the case, where you have three corners of one cube face on one side of the surface and their neightbours on the other side.

Four vertex

In this case the cube is split in half, but not along any particular axis.

Snake

That was the last case I’ve found and I really couldn’t find better name for it :). I like to imagine it as starting from some corner, then moving along one of the axes, then along another one and the last one. It uses every axis only once and all axes will be used. This forms a line made of 4 corners.

When I started to work on this implementation I assumed that I will probably not find an algorithm for every combination and I’ll eventually need to add some precomputed data, but after I’ve found the last case described here it seems that my cube case table is full. Hovewer there are still some places that I’m not happy of, especially in code that creates triangles. I hope I’ll improve those parts some day.

Despite the fact that I have generated geometry for every cube case, it seems that there’re still some errors in my surface representation. I’ll try to track them down when I have some spare time.

In the next version of the project I plan to add some geometry buffers, smooth normals and some simple caching, so we’ll make the results a bit nicer :).

Here’s the project source: MarchingCubes_1

Algorithmic Marching Cubes