The algorithm was developed by William E. Lorensen and Harvey E. Cline as a result of their research for General Electric. At General Electric they worked on a way to efficiently visualize data from CT and MRI devices. The premise of the algorithm is to divide the input volume into a discrete set of cubes. By assuming linear reconstruction filtering, each cube, which contains a piece of a given isosurface, can easily be identified because the sample values at the cube vertices must span the target isosurface value. For each cube containing a section of the isosurface, a triangular mesh that approximates the behavior of the trilinear interpolant in the interior cube is generated. The first published version of the algorithm exploited rotational and reflective symmetry and also sign changes to build the table with 15 unique cases. However, due to the existence of ambiguities in the trilinear interpolant behavior in the cube faces and interior, the meshes extracted by the Marching Cubes presented discontinuities and topological issues. Given a cube of the grid, a face ambiguity occurs when its face vertices have alternating signs. That is, the vertices of one diagonal on this face are positive and the vertices on the other are negative. Observe that in this case, the signs of the face vertices are insufficient to determine the correct way to triangulate the isosurface. Similarly, an interior ambiguity occurs when the signs of the cube vertices are insufficient to determine the correct surface triangulation, i.e., when multiple triangulations are possible for the same cube configuration. The popularity of the Marching Cubes and its widespread adoption resulted in several improvements in the algorithm to deal with the ambiguities and to correctly track the behavior of the interpolant. Durst in 1988 was the first to note that the triangulation table proposed by Lorensen and Cline was incomplete, and that certain Marching Cubes cases allow multiple triangulations. Durst's 'additional reference' was to an earlier, more efficient isosurface polygonization algorithm by Wyvill, Wyvill and McPheeters. Later, Nielson and Hamann in 1991 observed the existence of ambiguities in the interpolant behavior on the face of the cube. They proposed a test called Asymptotic Decider to correctly track the interpolant on the faces of the cube. In fact, as observed by Natarajan in 1994, this ambiguity problem also occurs inside the cube. In his work, the author proposed a disambiguation test based on the interpolant critical points, and added four new cases to the Marching Cubes triangulation table. At this point, even with all the improvements proposed to the algorithm and its triangulation table, the meshes generated by the Marching Cubes still had topological incoherencies. The Marching Cubes 33, proposed by Chernyaev in 1995, is one of the first isosurface extraction algorithms intended to preserve the topology of the trilinear interpolant. In his work, Chernyaev extends to 33 the number of cases in the triangulation lookup table. He then proposes a different approach to solve the interior ambiguities, which is based on the Asymptotic Decider. Later, in 2003, Nielson proved that Chernyaev's lookup table is complete and can represent all the possible behaviors of the trilinear interpolant, and Lewiner et al. proposed an implementation to the algorithm. Also in 2003 Lopes and Brodlie extended the tests proposed by Natarajan. In 2013, Custodio et al. noted and corrected algorithmic inaccuracies that compromised the topological correctness of the mesh generated by the Marching Cubes 33 algorithm proposed by Chernyaev.
Algorithm
The algorithm proceeds through the scalar field, taking eight neighbor locations at a time, then determining the polygon needed to represent the part of the isosurface that passes through this cube. The individual polygons are then fused into the desired surface. This is done by creating an index to a precalculated array of 256 possible polygon configurations within the cube, by treating each of the 8 scalar values as a bit in an 8-bit integer. If the scalar's value is higher than the iso-value then the appropriate bit is set to one, while if it is lower, it is set to zero. The final value, after all eight scalars are checked, is the actual index to the polygon indices array. Finally each vertex of the generated polygons is placed on the appropriate position along the cube's edge by linearly interpolating the two scalar values that are connected by that edge. The gradient of the scalar field at each grid point is also the normal vector of a hypothetical isosurface passing from that point. Therefore, these normals may be interpolated along the edges of each cube to find the normals of the generated vertices which are essential for shading the resulting mesh with some illumination model.
