Isosurface extraction has become a ubiquitous problem in the scientific visualization literature. We are all familiar with Lorensen and Cline's Marching and Dividing Cubes algorithms [6, 1], Wilhelm and Van Gelder's hierarchical algorithms [15], and some of the more recent improvements based on locating the critical points of the field or reparameterizing the domain based on value with active lists [4], span filters [3], both active lists and span filters [13], or span-space [12]. While all of these algorithms fall under the broad umbrella of ``isosurfacing'' methods, they clearly differ in their domains of applicability. In this paper we will focus on those methods that generate a connected, topological surface while extracting. Such a surface is often stored as a list of vertex positions, and a list of triangle triplets (the three vertex indices which compose each triangle). This representation is distinct from the ``geometric primitives'' format, where each triangle contains its three vertex locations, rather than indices into a vertex list. The geometric representation is useful when the surface's only purpose is to be rendered. In this case, the triangle list can be sent directly to a graphics pipeline for rendering. In contrast, the topological surface is more useful when the surface will be operated on for other purposes such as model editing, boundary condition application, and decimation/refinement. In all of these cases, it is not just important that the surface ``look right'' but that it be topologically connected.
Isosurface extraction is not the only method used for constructing topological surfaces - indeed, the vision literature is rich with methods for reconstructing models from range data [2, 5, 14]. For the purpose of this paper, we will restrict our interest to topological surface generation from scalar field data. Additionally, while there are many algorithms for extracting isosurfaces from data mapped to unstructured meshes, we will not be discussing those here. For the moment, we restrict our domain of interest to topological surface extraction from scalar field data on regular grids. Such data is prevalent in medical applications, as data from CT and MR scanners is most often comprised of parallel slices (which can then be stacked up to build a regular lattice) or a volume of regular hexahedral elements (voxels).
Within the domain of surface extraction from regularly gridded scalar data, several algorithms have been published. The most intuitive of these, an extension of Marching Cubes, is the seed-growth algorithm [11]. This algorithm was not motivated by the need for topological surfaces (that result was merely a by-product); rather, it was introduced as a way of extracting a single surface connected to a starting node in order to restrict the surface domain. This algorithm's strength also turns out to be a hindrance, though, as the user often wants more than a single connected component as output. Furthermore, this algorithm turns out to be exceedingly difficult to parallelize because the relevant data can not be easily decomposed into discrete regions.
Another topological, regular grid extraction method is the ``Bin and Coalesce'' method - a generalization of Rossignac and Borrel's algorithm for surface decimation [10]. This algorithm takes the geometric data discussed above and, for each triangle, bins its vertices by location. When every triangle has been processed, vertices located in the same bin are coalesced into a single vertex, and a topologically connected surfaces is output. The drawbacks to this method are that it requires a fair amount of memory overhead, does not necessarily preserve the topology of the initial surface, and must be applied as a post-processing step to an already-extracted surface.
Montani et al. implemented a discrete surface construction method that bisects, rather than linearly interpolates, the intersected edges [7]. The advantage to their algorithm is that the output surface can be easily decimated, as the resultant vertices lie within a finite number of indices. However, as with the original Marching Cubes algorithm, they must postprocess their data in order to obtain a connected, topological surface.
Finally, there is the edge-hashing method, as described at the end of Wilhelm and Van Gelder's octree isosurfacing paper [15]. This method uses hashing to create connected surfaces during extraction, rather than after. It works from the observation that if linear interpolation is used to locate surface-lattice intersection, then each edge in the lattice contains at most one vertex of the isosurface. Additionally, every isosurface vertex located on an interior (non-boundary) edge of the volume, will be found by exactly three other voxels - namely those voxels which share that edge. Working from these observations, Wilhelm and Van Gelder add each vertex location to their list the first time they encounter it, and then store the new vertex index in the hash table based on its edge. When this index is needed by the adjoining cells in the future, it is retrieved from the hash table. Finally, in order to conserve memory, they remove each interior vertex from the hash table after the fourth (and, by definition, final) time it is accessed. This method is very efficient, and succeeds in building a topologically correct surface as output. Furthermore, it can be parallelized in a straightforward manner - the only complication coming from sharing information about vertices on the boundaries between processors.
While these algorithms capitalize on much of the continuity of the data, what they fail to exploit is the predictability with which the data is accessed. In short, none of the above algorithms has very good cache performance, and hashing, in particular, has notoriously bad cache performance. As memory performance becomes an increasingly important factor in overall algorithm performance, it becomes necessary to give high priority to memory performance when designing algorithms. This is precisely what we have done in implementing our new topological isosurface extraction algorithm: recognizing the possible structure of memory accesses for topological isosurface construction, we have built a data structure and algorithm to match these access patterns.