5
$\begingroup$

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). 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.

Improve this question
$\endgroup$
5
  • $\begingroup$ I'm curious how you plan to detect if an arbitrary collection of quads or triangles forms a closed surface. $\endgroup$
    – Bill Barth
    Commented Apr 3, 2014 at 21:55
  • $\begingroup$ Unless I'm mistaken, every edge being incident to exactly two faces should be a sufficient condition for closure. Equivalently, I think, you could check that the faces incident to each vertex form a closed fan. $\endgroup$
    – RGMyr
    Commented Apr 3, 2014 at 22:36
  • $\begingroup$ If I do in fact end up using VTK, it looks like this is pretty easy to check with the vtkFeatureEdges class: cmake.org/Wiki/VTK/Examples/Cxx/PolyData/ClosedSurface $\endgroup$
    – RGMyr
    Commented Apr 3, 2014 at 22:42
  • $\begingroup$ You could easily have a self-intersecting topology which satisfies your edge condition but has effectively no inside or outside. $\endgroup$
    – Bill Barth
    Commented Apr 3, 2014 at 23:05
  • $\begingroup$ An additional requirement, that perhaps I didn't spell out, is that the surface may enclose only one local minimum. Please do correct me if I'm wrong, but doesn't a self-intersecting surface imply that some line drawn from the minimum could intersect the surface twice (or more)? For that to be the case, there would need to be another local minimum between the two intersection points, which will be prohibited. The proposed method for finding the surface is actually to increase the contour value until it encloses another minimum, and then step back to the previous contour value. $\endgroup$
    – RGMyr
    Commented Apr 3, 2014 at 23:35

1 Answer 1

Reset to default
3
$\begingroup$

I would still try to see if it is possible to extract the information from (or possibly to extend) VTK. VTK is such a widely used and high quality library that you would probably spend the better part of a year just reimplementing the pieces you need from VTK, and you'd still be only a fraction as fast.

At the very least, you should try to understand the code in VTK that computes isosurfaces.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.