Finding closed equipotential surfaces on a 3D grid
In short, I'm looking for either: (1) Publications or other sources dealing with contour/isosurface finding algorithms, so that I can write my own implementation (and parallelize as best I can), or (2) Suggestions about high-level libraries in Fortran, C/C++, or other languages, that may be useful for finding such surfaces from a scalar 3D grid.
I am currently analyzing data in checkpoint files from a very large hydrodynamics simulation (AMR w/ Lagrangian sink particles -- data in HDF5 format), and would like to find an algorithm for efficiently computing equipotential surfaces (a.k.a. countours, or isosurfaces) from the gravitational scalar potential grid. In particular, I need to be able to find the largest closed equipotential surface around a chosen local minimum (to within some acceptable step error), that may not enclose any other minima. This surface will be used as a selection criteria for further analysis of the volume formed by the enclosed cells. This link (particularly section 3) provides a description of the procedure I'm roughly attempting to replicate, though the IDL code provided takes 2D FITS maps of column density as input, and their data is much smaller: http://www.astro.umd.edu/~hgong/GRID_core.htm
This seems to be very similar to the problem of computing and rendering isosurfaces in graphics (i.e. some variant of the marching cubes algorithm seems relevant), but since I'm interested in using the surface as a boundary for selecting all interior grid cells, rather than visualizing it, most of the libraries I'm familiar with seem ill suited (e.g. VTK). 180
I would happily be proven wrong, however, if someone knows how I might use a graphics algorithm for this purpose -- maybe there's a way to find the isosurface with a graphics library routine and then simply extract information about the surface location or boundary cells?
Any suggestions about how to approach this problem, or ideas about potentially useful libraries are much appreciated. If I write the algorithm more or less from scratch, I will likely use Fortran (easier for me to parallelize), but I would be happy to use C/C++ or other languages if they support libraries that are better suited to the task.
Thanks in advance!
Update: After a bit more research, I'm optimistic that a contour tree will be the perfect data structure for this problem. Here is a rather informal description. It seems to provide a wonderfully abstract graph representation of how the topology of the set of isosurfaces varies with isovalue.
Local maxima and minima are represented as leaf nodes, and are said to create or terminate a "component" (graph edge). Interior vertices/nodes represent the joining or splitting of components, or the changing of genus, at critical points -- topological events, in other words. It is also possible to sweep through the space between the isosurfaces associated with two events, and assign all data points in this region to the component connecting the two events. This more thorough publication calls this an augmented contour tree.
With this in mind, my current strategy is pretty straightforward: once an augmented contour tree has been produced, then the cells I want to associate with a local minimum are simply those associated with the component between the local minimum leaf node and the closest topological event vertex. 180