Rendering volumes using Marching Cubes

 

Marching Squares (2D case):

Given a function f(x,y), where (x,y) are pixels in an image, marching squares is a way to approximate the curve along f(x,y) = 0.

For example, consider the function below (which you can edit), evaluated over the unit square:

To the right you can see a very low resolution (10×10) rendering of this function. Suppose we want to know the shape of the curve where this function has its roots (that is, where f(x,y) = 0).

Ideally we'd like to know this without having to evaluate the function at more samples.


Marching squares provides a way to get a sense of what a level-set curve of a function looks like, without taking more samples.

The key insight is that the curve can be approximated just by looking at those pixels bounded by corner points (i,j),(i+1,j),(i+1,j+1),(i,j+1) for which the signs of f at the four corners are not all the same. If the signs of f are different at two adjoining corner points of a pixel's square, that means the curve will cut the edge which connects those two corners.

One thing we need to figure out is where this transition happens along each such edge.

Given a value of A at corner a, and a value of B at adjoining corner b, we can compute the proportional distance t of the transition point along the edge [a,b] by observing, by similar triangles:

     t/A = (1-t)/-B

     -Bt = (1-t)A

     -Bt = A - tA

     (A-B)t = A

     t = A / (A-B)
Each corner can have two states: f < 0 or f ≥ 0, so in general, there are sixteen cases, as shown in the diagram to the right. Consider the second case along the top row of the diagram, where f at the top left corner (i,j) of a pixel is positive, but is negative at the other three corners of the pixel.

In this case, there is a transition point p along the top edge -- between (i,j) and (i+1,j), and another transition point q along the left edge -- between (i,j) and (i,j+1). Within this pixel, we can approximate the f(x,y)==0 curve by the line segment [p,q].

So that for any pixel we need to do three things:

  1. Figure out which edges, if any, of the pixel contain transition points
  2. Compute the locations of these points;
  3. Draw line segments between transition points, to approximate pieces of the curve.

Marching Cubes (3D case):

Marching cubes is the 3D equivalent of marching squares. Rather than approximate a closed curve where f(x,y)=0 via small straight edges inside square pixels, as in marching squares, the marching cubes algorithm approximates a closed surface where f(x,y,z)=0 via small triangles inside cubic voxels. The technical paper describing this algorithm, published by Lorensen and Cline in 1987, has been cited more often than any other paper in the field of computer graphics.

Each voxel cube has eight corners, which can be numbered as follows:

0 x=0 y=0 z=0
1 x=1 y=0 z=0
2 x=0 y=1 z=0
3 x=1 y=1 z=0
4 x=0 y=0 z=1
5 x=1 y=0 z=1
6 x=0 y=1 z=1
7 x=1 y=1 z=1

Because the value of f(x,y,z) at each of these eight corners can be either positive or negative, there are 28 or 256 cases to consider. These are shown in the figure to the right.

The key is to put all of this information into this table. The table has 256 entries, one for each of the 256 cases. Each entry contains between 0 and 4 triangles, which is the number of triangles that will be produced by the marching cube algorithm for a voxel of that type.

Each triangle is described by the three edges of the cube that contain its respective vertices, and each vertex is described by identifying one cube corner as well as the orientation of the cube edge that contains that vertex.

For example, a particular vertex of a triangle in the table may be described by the number sequence 0,1, indicating that this vertex lies on edge [0,1] of the cube. This is the edge that connects the x=0 y=0 z=0 corner of the cube and the x=1 y=0 z=0 corner of the cube.

Marching Tetrahedra (simpler to implement, somewhat less efficient):

To avoid the big table look-up of Marching Cubes, the technique I generally use is to split up each voxel into six tetrahedra. Given the same corner numbering we used for Marching Cubes, we can partition the voxel cube by "turning on" the binary bits of the numbered corners in different orders, giving the six tetrahedra:

[0,1,2,7] , [0,1,5,7] , [0,2,3,7] , [0,2,6,7] , [0,4,5,7] , [0,4,6,7]
Since a tetrahedron has only four edges, there are only two non-trivial boundary cases: (1) the boundary is a single triangle, or (2) the boundary is a four sided shape, which can be split into two triangles.

This algorithm is a little less efficient than Marching Cubes, because it generally produces more triangles for each boundary cube. However it requires much less code, and therefore is easier to program and to debug.