# Extract Surface (Marching Cubes)

The Extract Surface algorithm described in this section is based on the Marching Cubes algorithm, which was designed by William E. Lorensen and Harvey E. Cline to extract surface information from a 3D array of values. The applications of this algorithm are mainly concerned with medical visualizations such as CT and MRI scan data images.

### Background

The basic notion of the algorithm is that we can define a voxel by the pixel values at the eight corners of the cube. We also can define an isosurface which intersects the cube, refer to Figures 1 and 2. If one or more pixels of a voxel have values less than the user-specified isovalue, and one or more have values greater than this value, we know that the voxel must contribute some component of the isosurface.

By determining which edges of the voxel are intersected by the isosurface, we can create triangular patches that divide the voxel between regions within the isosurface and regions outside. By connecting the patches from all voxel on the isosurface boundary, we get a surface representation.

The indexing convention for vertices and edges used in the algorithm is shown in Figure 1 below. Example: If for example, the value at vertex 3 is below the isosurface value and all the values of all the other vertices are above the isosurface value, then we would create a triangular facet which cuts through edges 2,3, and 11 as shown in Figure 2. The exact position of the vertices of the triangular facet depends on the relationship of the isosurface value to the values at the vertices 3-2, 3-0, 3-7 respectively. Therefore, we must specify a threshold value, first. For this value, some voxels will be entirely inside or outside the corresponding isosurface and some voxels will be intersected by the isosurface. In the first pass of the algorithm, the voxels that are intersected by the threshold value are identified. In the second pass, these voxels are examined and a set of one or more polygons is produced. These polygons are then output for rendering.

Each of the voxel's 8 vertices can be either inside or outside of the isosurface value, thus, there are 2'8' = 256 possible ways how an isosurface can intersect the voxel. To simplify the algorithm, we can reduce this complexity by assuming that cell combinations are duplicated under the following conditions:

• Rotating by any degree over any of the 3 primary axis (X,Y, or Z);
• Mirroring the isosurface's shape across any of the 3 primary axis;
• Inverting all corners and flipping the normals of the relating polygons.

By taking this symmetry into account, we can reduce the original 256 combinations to a total of 15 combinations, refer to Figure 3. Now we can use a table (called edgeTable) which maps the vertices under the isosurface to the intersecting edges. An 8 bit index can be formed from the table where each bit corresponds to a vertex, refer to [MarchingCubes.html#wp1780799 TableÂ 1].

Table 1. edgeTable
cubeindex = 0;
if (grid.val < isolevel) cubeindex |= 1;
if (grid.val < isolevel) cubeindex |= 2;
if (grid.val < isolevel) cubeindex |= 4;
if (grid.val < isolevel) cubeindex |= 8;
if (grid.val < isolevel) cubeindex |= 16;
if (grid.val < isolevel) cubeindex |= 32;
if (grid.val < isolevel) cubeindex |= 64;
if (grid.val < isolevel) cubeindex |= 128;

Looking up the edgeTable returns a 12 bit number, each bit corresponding to an edge. 0 if the edge isn't cut by the isosurface and 1 if the edge is cut by the isosurface. If none of the edges are cut the table returns zero. This occurs when the cubeindex is 0, which means that all vertices are below the isosurface, or 0xff (all vertices above the isosurface).

Using the example shown in Figure 1, where only vertex 3 was below the isosurface, cubeindex would equal 0000 1000 or 8. The edgeTable{8} equals 1000 0000 1100. This means that edge 2,3, and 11 are intersected by the isosurface, refer to Figure 2.

Now the intersection points can be calculated by linear interpolation. If P1 and P2 are the vertices of a cut edge and V1 and V2 are the scalar values of each vertex, the intersection point P is given by the following equation:

The last part of the algorithm involves forming the correct triangular facets from the positions where the isosurface intersects the edges of the grid cell. In order to do that, another table (also called triTable) is used. This table uses the same cubeindex from edgeTable, but allows the vertex sequence to be searched for as many triangular facets as necessary to represent the isosurface within the grid cell. Each of the non-trivial configurations results in between 1 and 4 triangles is added to the isosurface. The actual vertices themselves can be computed either by interpolating along the edges, or, by default their location to the middle of the edge.

Since we can create surface patches for a single voxel, we can apply this algorithm to the entire volume.We can either treat each cube independently, or we can propagate edge intersections between the cubes which share the edges. The sharing of edge/vertex information also results in a more compact model.