Marching Cubes Tutorial

This is a short tutorial on the Marching Cubes Algorithm with C++ source code.

If you want a fuller tutorial on this subject please read an article by Paul Bourke.

There are a number of other tutorials on the web:
http://www.cfxweb.net/article.php?sid=199
http://www.cosc.brocku.ca/Project/info/grossi/detmarch.html
http://www.cs.wpi.edu/~matt/courses/cs563/talks/march_cub.html
http://www.siggraph.org/education/materials/HyperVis/vistech/volume/surface4.htm
http://www.whisqu.se/per/docs/math29.htm

Marching Cubes algorithm, developed by Lorensen and Cline in 1987 is used to approximate an isosurface by subdividing a region of space into 3D array of rectangular cells. Each of the eight vertices of a cell is assigned a value, comparing which against the minimum determines if that cell is intersected by the surface. Any vertex with a value less than or equal to the minimum is defined to be inside the surface. For each inside vertex a bit is inserted into a designated index (called cubeIndex in my source code) at a position representing that vertex. The index is used with edgeTable (defined in MCTable.h) to return the edges of the cube that intersect the surface, again as bits of an integer. If none of the edges are intersected the process goes back to the beginning and tests another cell. Otherwise that number is checked bit by bit and a point of intersection is determined (usually by Linear Interpolation, see this page for details) for each intersected edge. That point is saved into an array at the index corresponding to its edge (see Paul Bourke's article for numbering). Finally triangles are created by connecting points on intersected edges of the cell in the order returned by triTable at index formed earlier (cubeIndex). triTable is 256x16 table of integer values, which are used as indices for the array of 12 points of intersection. It defines the right order to connect the intersected edges to form triangles. The process for one cell stops when index of -1 is returned from the table, forming a maximum of 5 triangles. Here are the step in Marching Cubes Algorithm, each one is also shown in the source code (MarchingCubes.cpp):

1. Define space: (minx, maxx, miny, maxy, minz, maxz).
2. Subdivide space: either number of cells or the step size on each axis.
3. Assign a value to each vertex of each box.
   (the above is done outside of the algorithm, the steps below are what the algorithm actually does)
   
4. Compare the value at each vertex of each box to the designated minimum value. This will determine if vertex is inside or outside. For each vertex that is inside a bit is inserted into the index (cubeIndex). Again for numbering of the vertices consult Paul Bourke's article.
5. We index the edgeTable array using cubeIndex from step 4. That in turn returns a number, bits of which represent the 12 edges that intersect the cube (there are 256 possibilities). Zero is returned if all are inside or outside, at which case we don't go any further.
6. If cube is intersected by the surface we assign a point of intersection for edge that is intersected, as determined by the bits in the number returned from edgeTable (using Linear Interpolation). All of the points are put in an array of size 12, each representing its edge.
7. cubeIndex is used again with triTable which returns indexes for the array of points defined in step 6. triTable is a double 256x16 array, representing 256 possibilities of intersection and up to 15 vertices for the triangles formed in each possibility. The points are in order, so for example triTable[3][0] will be the index of the first vertex of a first triangle of the third possibility, and so on. Those indexes are used with array defined in step 6 to group vertices in appropriate triangles.
8. (Optional) Normal vector is defined as the cross product of two vectors for each triangle.
   
   

This is all for the basic Marching Cubes algorithm. Once each step is understood it can be coded in no time. I have provided 2 functions, one that runs Marching Cubes on an array of already computed points and the other which first computes points using some given formula and then executes the first function with that data.

First function (version 1) takes the number of cells to subdivide the space on each axis, the minimum value to be used with during linear interpolation, the array of points (type mp4Vector) with already defined coordinates and values (mp4Vector contains 4 values, see mpVector.h for details). The size of the array is (ncellsX+1)(ncellsY+1)(ncellsZ+1) where ncells is the number of subdivisions on each axis. Final argument is integer in which number of triangles formed is saved.

The second function does more work for your. It starts with step one and goes down to step 8, while the first function needs steps 1 to 3 done before hand. It accepts maximum and minimum values and a number of subdivisions for each axis, a minimum value used during linear interpolation, a function which computes value based on the point passed to its only argument (type float(mpVector)), and the integer in which the number of triangles is saved.

Download here (534 kb): MCExample.zip

Source Code

Please email me with suggestions and/or bugs at:
myp@andrew.cmu.edu or mikepolyakov@hotmail.com