1.4. Isosurface Generation

One technique (among many others) of visualizing a 3D volume is to extract isosurfaces and render those as a regular polygonal surface. It can be used to extract the surfaces of bones or organs in medical scans , such as MRIs and CTs.

Assume for the moment that we are given a scalar field f : ³ ’ . Then the task of finding an isosurface would just be to find all solutions (i.e., all roots) of the equation .

Since we live in a discrete world (at least in computer graphics), the scalar field is usually given in the form of a curvilinear grid : the vertices of the cells are called nodes , and we have one scalar and a 3D point stored at each node (see Figure 1.8). Such a curvilinear grid is usually stored as a 3D array, which can be conceived as a regular 3D grid (here, the cells are often called voxels ).

Figure 1.8: A scalar field is often given in the form of a curvilinear grid. By doing all calculations in computational space, we can usually save a lot of computational effort.

The task of finding an isosurface for a given value t in a curvilinear grid amounts to finding all cells of which at least one node (i.e., corner) has a value less than t and one node has a value greater than t . Such cells are then triangulated according to a look-up table (see Figure 1.9). So, a simple algorithm works as follows [Lorensen and Cline 87]: compute the sign for all nodes , and then considering each cell in turn , use the eight signs as an index into the look-up table, and triangulate it (if at all).

Figure 1.9: Cells straddling the isosurface are triangulated according to a look-up table. In some cases, several triangulations are possible, which must be resolved by heuristics.

Notice that in this algorithm, we have used only the 3D array; we have not used the information about exactly where in space the nodes are (except when actually producing the triangles ). We have, in fact, made a transition from computational space (i.e., the curvilinear grid) to computational space (i.e., the 3D array). So in the following, we can, without loss of generality, restrict ourselves to consider only regular grids, that is, 3D arrays.

The question is: how can we improve the exhaustive algorithm? One problem is that we must not miss any little part of the isosurface. So, we need a data structure that allows us to discard large parts of the volume where the isosurface is guaranteed to not be. This calls for octrees.

The idea is to construct a complete octree over the cells of the grid [Wilhelms and Gelder 90] (for the sake of simplicity, we will assume that the grid s size is a power of two). The leaves point to the lower-left node of their associated cell (see Figure 1.10). Each leaf ½ stores the minimum ½ _min and the maximum ½ _max of the eight nodes of the cell. Similarly, each inner node of the octree stores the minimum/maximum of its eight children.

Figure 1.10: Octrees offer a simple way to compute isosurfaces efficiently ..

Observe that an isosurface intersects the volume associated with a node ½ (inner or leaf node) if and only if ½ _min ‰ t ‰ ½ _max . This already suggests how the algorithm works: start with the root and visit recursively all the children where the condition holds. At the leaves, construct the triangles as usual.

This can be accelerated further by noticing that if the isosurface crosses an edge of a cell, that edge will be visited exactly four times during the complete procedure. Therefore, when we visit an edge for the first time, we compute the vertex of the isosurface on that edge, and store the edge together with the vertex in a hash table. So, whenever we need a vertex on an edge, we first try to look up that edge in the hash table. Our observation also allows us to keep the size of the hash table fairly low. When an edge has been visited for the fourth time, we know that it cannot be visited anymore; therefore, we remove it from the hash table.

Figure 1.11: Volume data layout should match the order of traversal of the octree.