Marching cubes is an algorithm to create a triangle mesh from an isosurface. It is used to render terrain or fluids in 3D for example.

Let's have an example in 2D. We start from a 2D grid and each point has a value. The goal is to separate points that are below and above a certain value, called isolevel. In the example, the isolevel is 0.

It's the same idea in 3D. To simplify the process, we divide the grid in cubes and create the isosurface in each of them. There are 2⁸ = 256 possibilities, but it goes down to only 15 if we allow rotations, plane symmetries and flipping the values (for example, a grid full of negative numbers is the same as a grid full of positive numbers).

Actually, if look closely at those cases from the original publication, case 14 and 11 are the same, if we flip along the X plane. So in fact there are 14 cases.

So, the goal is to identify in which case we are for each for all of the cubes and add the triangle to the mesh with the right transformation. It is guaranteed that adjacent cubes will have a surface that connects.

But it takes too much time to find the best case for each cube, so instead we can precompute all the 256 cases with the correct transformation, so we just have to look into a table to get correct surface.

What people do generally is to find the table on the internet and paste it in their code. But I had spare time so I decided to create mine.

I first associated a number to each corner of the cube (in black), and a number for each of the segments (in green). Our table will associate a configuration to a list of triangles. Each triangle is composed of points at the center of a segment. For example, in case 1, we have one triangle with points 2, 3 and 11.

We can represent the corners by binary numbers (above or below the isolevel), so each configuration can be represented by an 8-bits number, between 0 and 255. For example, in case 1, the corner 3 is activated, so the binary number is 00010000, which is the number 16.

Then I created all the base cases by hand and create a base table. The next step is to generate all the remaining cases. For this, I created a Python script that will brute force every transformation for all the configurations to find what base case it is associated to. A transformation is just a permutation of the corners. Here is a rotation around the Y axis:

The transformations we need to consider are the rotations around each 3 axis and one symmetry on the YZ axis. We can go from any configuration to a base case by combining those transformations: up to 3 rotations and up to one symmetry. Also, when there are more than four corners that are activated, we can invert them, because we get the same surface. We just have to flip the normal in that case. We need to flip it too when there is a symmetry.

Once we have our transformation, we can apply it in reverse to the segments and get the surface. After a lot of debugging, here is what I obtained.

The source code can be found here.

Now, let's see the results. I can fill a 3D grid with values from a Perlin Noise, and I obtain this.

But the result is quite blocky. The first thing we can do is to interpolate the point on the segment: we should put the point closer to a corner if it is closer to the isosurface than the opposite. This create a smoother shape:

But the normals are not good. We should compute a normal at each point instead of doing it for the face. The problem is that my mesh is not optimized for the moment. I just add the triangles for each cubes, without using the ones that are already placed. There are 3 points for each face, although some could be used for several faces.

I fixed this by using a hash map that associates an edge to an index. In my example, I went from 195 000 vertices to 33 000.

Also, when interpolating, some points from adjacent edges are too close, and this creates a micro face, which is bad when we compute the normal. So I also merge those points. To do this, I used a spatial hashing grid to accelerate the search. I went from 33 000 vertices to 30 000.

Now my mesh is correctly built, I have no duplicate points, so I know which faces are attached to which point. Then I can average the normals of the adjacent faces to get one normal per vertex.

Finally, I used this algorithm to render fluids. I had already implemented a 3D fluid simulation before, and I can build an isosurface from the fluid density.

I had to write a simple shader to add a skybox and some reflexions, but here is the final result: