PROPOSAL FOR RESEARCH TO BE CONDUCTED FOR
Naval Research Laboratory
In Response to BAA96-015

Advanced Information Technology Branch
800 North Quincy Street, Arlington, VA, 22217-5660
c/o Dr. Larry Rosenblum, ONR 311 rosenblum@ait.nrl.navy.mil
tel: (703) 696-0990

Novel Morse Theoretic Tools for
Virtual Reality Exploration of Volumetric Datasets:
Graphing the Visible Human Project Dataset

Project Period: 2 + 3 years, from 1/97 to 12/99 , 1/00 to 12/02.
Amount Requested: $621,168.53.
($ 121,228.80) ($ 124,865.66) ($ 128,611.63) ($ 132,469.98) ($ 136,444.08)

Principal Investigator:
Daniel B. Karron, Ph.D. e-mail: karron@casi.net
Computer Aided Surgery, Inc.

James Cox, Ph.D. e-mail: karron@casi.net
City University of New York Graduate Center and Brooklyn College

Officer To Whom Award Documents Should Be Mailed:
Daniel B. Karron, Ph.D. e-mail: karron@casi.net
President and Chief Technical Officer, Computer Aided Surgery, Inc.
300 East 33rd Street, Suite 4N, New York, NY 10016,
tel: (212) 86-8748 fax: (212) 685 - 0736

Abstract

We propose to create the technology to generate an ``Object Spectrum Graph'' (OSG) (a generalized Reeb Graph, applied to the hypersurface defined by the entire scalar function) of all of the available separable objects in a sample data set. We will use as our sample dataset the Visible Human Project data. The OSG will be a valuable digraph index of each and every separable object, their parantage and children. The OSG will enable us to construct a visual graph representation of all possible topologically distinct objects over all threshold isovalues. For each family of isosurfaces (for a given isovalue threshold range) each bounded object includes every detail, but can be drawn as a single shell with as much or as little Level of Detail (LOD) tileing as required from any particular Point of View (POV).

Should we find a particular detail worth closing in on from our POV, we can recurse into only that part of the object that is visible in our viewing frustrum or cube, and only as deeply as our pixilated screen space warrents. We can vary the isovalue and display the changes to the associated isosurfaces, thus allowing for both spatial and function value navigation of the dataset. We have written a 2D version of the OSG code and desire funding in order to implement this theory as a new collection of algorithms for higher dimensions.

We will break new ground in many fundamental geometric and topological theoretic problems as we proceed with the implementation of our prototype system. We will be using the theory of (Marsdon) Morse Functions recast in a digital context (Digital Morse Theory) as a foundation for a new LOD management technology.

In our Phase 1 effort we will write indexing and viewing software in 3 dimensions available in 2 years that will run via VRML from a SGI supercomputer host on the internet.

We wish to mount a follow up phase 2 effort wherein we request an additional 3 years funding so that we may implement a viewer capable of managing up to 6 dimensions of information (using dimensions of time, or RGB Color) that will have the capability to track data objects changing over time (such as the beating heart), or enable any linear (or other) combination of multi-spectral information. We also will develop application of the OSG for registering as well as understanding multispectral data as separate or combined dimensions of information.

Introduction

This project will investigate and apply an extension of Morse theory to the analysis of volumetric data: data that represents the value of some real-valued or real vector-valued function sampled at discrete points in space. We call this extension Digital Morse Theory, and in the case of a real valued function it allows one to identify and render (using the SpiderWeb Algorithm) each topologically distinct isosurface bounded object over the entire range of threshold values. That is, one can completely explore, apprehend and navigate the volume data, at any appropriate level of detail. We will apply the theoretical and practical insights gained by the investigators to extend the results to vector valued functions. In particular we will apply our results to the visual human data set, which is a vector-valued (RGB) volume dataset. We will produce an object spectrum graph (which generalizes the notion of Reeb graph) with associated topological information, representative surfaces for each topologically distinct object, and the means for navigating the dataset in a variety of ways. For a specific linear combination of RGB values and isovalue we can visually navigate the volume spatially and at changing resolution scales. We can vary the isovalue and see the resulting changes to the isosurfaces (at critical points). Finally we propose to allow navigation in RGB space, that is, varying the coefficients of the convex combination of RGB values and displaying the topological changes to the isosurface bounded objects.

There have been various similar approaches in the published literature. The key to our approach is that we can correctly identify the critical isovalues at which the topology changes in the corresponding isosurfaces in discrete data. This applies the notion of Morse theory to criticalities as the isovalue ranges, as opposed to applications of Morse theory on individual isosurfaces, as has been done in the literature to date[Kergosien92][Itoh96][Shinagawa92]. This virtual reality navigation of volume data might seem an ambitious visualization project but it is precisely because the Digital Morse Theoretic approach is so powerful that it enables such general apprehension of volume and visual data. It should be noted that the insights of Morse theory have been applied by researchers, but always to the topologically correct rendering of a single object or several objects for fixed isovalues or for the static segmentation of a picture into contiguous similiarly (greyscale) valued regions of pixels ( [Gauch92] has a good review).

In contrast, over the last few years we have been applying Morse theoretic ideas in a higher dimension: Regard the data as defining a single hypersurface; varying the isovalue results in specific slices through this surface represented by the isosurfaces for that value. By identifying criticalities and associated changes to the number and topology of objects as the isovalue varies, we can characterize all possible isosurfaces and hence all possible segmentations of the scene into separable topologically distinct objects. We can still view the data at multiple resolutions, without building models at multiple resolutions, but rather by filling in the details at a given resolution, as needed. The overall structure of all possible isosurfaces is known and represented in the OSG, and the sample isosurfaces for each critical value bounded interval of function values.

This work represents a fundamental contribution to the understanding and analysis of data representing function values, sampled at discrete points in space (volumetric data). Such data is ubiquitous, occurring in many applications areas. Many of the current techniques multiresolution are tailored for the particular case of segmenting visual data into objects defined by similar pixel intensity and employ the usual gradient edge detection, Gaussian blurring, and other techniques drawn from computer vision and signal processing. Our techniques extend to other situations where segmentation based on visual cues may not be appropriate. For some examples, one may be looking for isothermal patterns and structures embedded in atmospheric data, or looking to identify regions of similar velocity in data representing a fluid flow.

A system for the critical point analysis of two dimensional data has been developed and the theory and algorithms have been extended (in a simple and natural way) to three dimensional volumetric data. The analysis seems to extend quite naturally to higher dimensions, and investigator Karron is employing these ideas in the four dimensional analysis of heart data taken over time. The investigation will implement the software tools to solve a variety of problems: segmentation, visualization, navigation at varying scales of the data space, and a priori (prior to growing any surface) simplification of the data. This is novel; while the computer science literature is brimming with data compression and simplification schemes, this is an attempt to simplify spatial functional data in a manner that maintains the essential topology of the function's graph. This is to contrasted with maintaining the topology of one or several surfaces for a fixed function isovalue, which may be thought of as a 3D slice through the function's graph. By characterizing all possible topologically distinct and consistent objects over the entire range of isovalues, it can assist in the aforementioned and more standard image processing techniques of segmentation and multiresolution shape description. More precisely the analysis we propose identifies all topologically distinct objects, at all possible thresholds, and thus subsumes all topologically consistent segmentations of the data.

Some segmentations (for example by gradient edge detection) do not produce objects bounded by isovalued surfaces. By mapping the basins of attraction of DMT criticalities (as described in the section on applications of the theory) we can assist in this type of spatial division of the volumetric data (for volume rendering).

The analysis that leads to these new and powerful results represents a novel synthesis of techniques and theory drawn from discrete and combinatorial mathematics (especially computational geometry), topology (including digital topology) and the more classical analysis of continuum mathematics, such as the calculus of multivariate functions and the calculus of variations. Potential applications are in areas as diverse as scientific imaging and visualization, solid and geometric modeling, computational fluid mechanics, geophysical modeling (including climatology), geometric and topological analysis of morphogenesis, as well as the imaging problems of object registration, object matching, segmentation (what is the correct isovalue), computing important topological invariants and variances, surface simplification (reducing the number of polygons in a tiling), and data simplification (removing lattice points without changing the topology of any isosurface for any isovalue). The fundamental analysis and categorization of all possible level set topologies can aid in the search for geometric patterns in the data. It also provides a good starting point for model based and statistical reasoning about discrete data.

Theoretical Background

There is much common ore to be mined in the disciplines of solid and geometric modeling for CAD/CAM, scientific visualization, and robotics. The recent work of investigator Cox has been in the areas of sensory based robot navigation or motion planning, [Cox89], [CoxYap91], [CoxWang92], complexity theoretic analysis [ChenCoxMishra91], [CoxMishraEricson96] and constraint programming [CoxMcAloonTretkoff92][CoxMcAloonTretkoff90][CoxMcAloon93]. Investigator Cox, in collaboration with Dan Karron and Bud Mishra of NYU, has worked in the related areas of scientific imaging and solid modeling [KarronCox95][KarronCoxMishra94][CoxKarron93][Karron92ParisPoster].

Let us consider first the basic problem of a real-valued function of three spatial variables. Stated briefly one is given a set of real numbers at discrete points in space, where the numbers represent a real-valued function, such as the density of some material. We assume that the numbers are given at the points of the integer lattice, but this is for simplicity of presentation only. The data can in fact be anisotropically sampled. The basic problem is to segment the points (or pixels) into contiguous objects of similar density, bounded by an interface surface of some specific threshold value. One approach segments the data and then attempts to recover a consistent bounding surface for each object. Alternatively, one can choose a specific function value and plot the level sets (isosurfaces) for this isovalue, with an object defined as the topological interior of this surface. The problem then is to interpolate a surface representing the points of some particular function value, e.g. the density of the skin surface of a scalp or the surface of a brain, or the points in a fluid shear surface. This is known as the isosurface construction problem.

Many algorithms have been proposed for this problem, and most desire desire to construct an isosurface that is an oriented manifold, see for example [Lorensen87], [Cline87], [Cline88], [Kalvin91]. The problems encountered by most isosurface construction algorithms occur because of what researchers interpret as noisy data. This seeming noise is manifested in anomalous situations where portions of the data will not yield a surface with well defined interior and exterior regions. Most surface construction algorithms attempt to grow oriented surfaces by using complicated case table driven techniques, that are often ad hoc. Even the best and most widely used of these algorithms [Lorensen87], the so-called Marching Cubes algorithm, can break down in the presence of noisy or anomalous data.

By contrast the recently discovered SpiderWeb algorithm [Karron92ParisPoster][CoxKarron93] will, under the assumption of a bounded region of nonzero density, always produce a collection of consistently oriented manifold surfaces, irrespective of noise in the data. The algorithm has been applied to medical imaging from CT and MRI at NYU Medical Center by the investigator and colleagues Dan Karron and Bud Mishra, and they have applied for a patent. Since the surfaces are given in triangulated form, we can compute the topological genus of each constructed manifold surface, with little additional overhead, using the Euler-Poincare statistic.

The alternate approaches to surface construction (mentioned above) include volume visualization and digital contouring. The dissertation [Kalvin91] contains a good survey of the area. The limitations of the pure voxel renderings approach have motivated the digital surface approach. In this approach, each lattice point or vertex is considered to be an atomic entity and is regarded as inside or outside some object. Adjacency relationships are defined between the entities, and a surface is assumed to exist ``between" any pair of adjacent vertices for which one is in, and the other is out. There is extensive work on developing a theory of digital topology using the natural generalization of this discrete rubric [Herman93][Udupa82][Udupa91][Udupa92].

The SpiderWeb produces an interpolated surface, consisting of a triangle tiling (simplicial complex). Note that the well-known Marching Cubes [Lorensen87][Cline87][Cline87a][Cline88] also produces a triangle tiling. Many of these algorithms work by means of a unspecified rules encoded into a look-up table [Lorensen87][Cline87][Cline88][Baker88][Artzy81][Bloomenthal88][Wallin91][Nagae93]. There has yet to emerge a consensus as to how to define topological correctness [Ning93][VanGelder90][Wilhelms90], and various approaches to solve the remaining problems all seem to work.

Various schemes to disambiguate [Kalvin91][Kalvin91a] or re-tile with triangles and/or tetrahedra [Frey91] apparently ``rare''[Nielson91] ambiguous circumstances have been implemented. Arbitrary but consistent decision algorithms nevertheless generate orientable manifold surfaces. Investigation of these issues has led to the broader understanding embodied in the emerging Digital Morse Theory.

NOTE: We wish to emphasize that the correct knitting of a single manifold isosurface has been acomplished in a variety of ways and is now considered a minor problem. The important point is that the examination of these ambiguous situations and the SpiderWeb algorithm provided the starting point of a much deeper analysis based on Digital Morse Theory. It is the understanding of the role of the underlying Morse Singularities that cause tileing difficulties that opens up new approaches to old problems.

The Basic SpiderWeb Algorithm

The data is assumed given at lattice points in space which we organize into adjacent cubels (the cubic volume of the lattice), see figure 2. Isovalue hits are interpolated on the cubel edges (figure 1), wherever the endpoint edges straddle the isovalue. A simply proved property is that each face can have only an even number of hits. We connect these hits to form a triangulated surface.

Here is a brief summary of the algorithm. Two hits in a cubel are considered adjacent if they are on the same 2-hit face, or they are adjacent on a 4-hit face. The hits on a 4-hit face are divided into 2 adjacent pairs according to a disambiguation rule given in the next section, see figure 4 . A connected system of hits within a cubel (usually all the hits if there are no 4-hit faces) is the transitive closure of the adjacency relation within that cubel. For each connected system in a cubel, the algorithm selects (interpolates) an interior articulating point (AP) (center of mass of the hits) and makes it the common vertex of all triangles grown among the connected system of hits, each triangle connecting the AP with a pair of adjacent hits in a face, as shown in figure 3 . Note that in some cases a single cubel may have more than one connected system of hits and thus more than one articulating point. The triangles constructed in a cubel are automatically connected to the triangles grown in adjacent cubels as each pair of adjacent hits is shared by two adjacent cubels. One can grow all surface portions independently in each cubel and track individual surfaces by following the triangle adjacencies.

Remarkably, this simple rule produces a collection of oriented manifolds, without use of complex case tables[CoxKarron93]. A given surface can be tracked in depth-first or breadth-first fashion, in time proportional to number of cubels that the surface intersects, a sparse set in the data lattice. Further, since the representation of each surface is given as a simplicial complex, typical geometric and solid modeling operators [Hoffmann89], [Mishra91] are easily accomplished.

The analysis of ambiguities in the data (what we call 4-hit faces) led to the recognition the relation of these ambiguities to topology-changing saddle critical points in the data and the development of DMT. DMT gives insight into the nature of the criticalities present. This work uses Morse theory in a different setting then Shinawaga et. al., [Shinagawa91a][Shinagawa91b], in the sense that our analysis of the density function is independent of the choice of isovalue; the height function is obtained by varying the isovalue. In the following section a brief introduction to Digital Morse Theory is given.

Figure 1: Illustration of an threshold crossing, or ``hit'' along an edge.

Figure 2: Cubel Nomenclature.

Figure 3: Illustration of a simple triangular webbing inside a cubel.

Figure 4: Four hit ambiguity nomenclature. Illustration of ``knit In'' and ``knit Out''
rules, where In is the higher value, ``knit high'', Out is the lower value, ``knit low''.

Figure 5: Triangle hnitting ambiguity creates a potential new point for resolution
of the ambiguity. We will show later that this is evidence for the presence of a
singularity inside the cubel face/body.

Digital Morse Theory

Recall a function f defined on a compact smooth manifold is Morse if it has a finite discrete set of critical points (points where the differential vanishes) and that at each of these critical points the Hessian is non-singular. The critical values are the values of the function at critical points. In the case considered in this proposal f is assumed to be Morse on a finite open region of 3-dimensional Euclidean space (containing the data lattice), and is identically 0 in small neighborhoods of the boundary of that region.

We are guaranteed to have no topological changes in the level sets { p | f(p) = tau } or { p | f(p) <= tau } for tau varied in the open interval between any pair of critical values. Further, the topological changes at a critical point can be completely characterized by the number of negative eigenvalues of the Hessian at the point.

These two facts constitute the fundamental result of Morse theory [Goresky88][Milnor63][Morse34]. We apply these results to the volumetric data.

Digital Morse Criticalities

The essential properties of any volume data set can be described by a hierarchy of critical points. Lattice points that have respectively greater and lesser sampled values than their lattice neighbors (defined under certain adjacency conditions given in the axioms) are respectively maxima and minima. This is because, for example, a maximal lattice vertex will spawn the creation of a new object when the isovalue is dropped below its value. This is true not only for SpiderWeb, but for all the isosurface algorithms cited above. Lattice vertices that are maximal with respect to one direction with respect to the two adjacent lattice points in one direction but minimal in another direction form saddles. Further, as we show below, some saddle points are implied by the selection of a disambiguation rule for 4 hit faces, and are interior cubel points.

NOTE: This holds regardless of the disambiguation rule chosen. As the isovalue is decreased through a range containing the penultimate high and the penultimate low on the face, a local topological change occurs. Between these two values the face will have 4 hits. From examining the pattern of readings on the face's vertices, one can determine if 4 hits can occur and the range of isovalues that engendered the hits.

In lattice data, the maxima or minima do not have to be a points. We allow for connected sets of lattice points forming plateaus, flat basins, flat straight ridges, and circular ridges (cirques). These type of digital (and hence non-Morse) criticalities are generated by binary scenes or connected sets of identically valued lattice points. Nevertheless, in passing through the data we can identify all critical isovalues and the above critical point structures.

One can characterize the importance of a critical point by the size of the region that it owns or inherits (its basin of attraction) and thereby eliminate critical points that represent small innesential or textural features.

The key point is that if one identifies all critical points and associated critical isovalues, one can then divide the range of function values into critical value-free intervals. The set of isosurfaces for any value in such an interval will be topologically equivalent to the set of surfaces for any other value in the same interval. Size and geometry will continuosly vary as the isovalue is varied, but the number and topological genus of the surfaces is invarient within the interval.

Further, the index of a critical point (the so-called Morse data) determines what type of change occurs as the isovalue is varied through the corresponding critical value. For volume data this is quite easy to characterize. For example as the isovalue decreases through a local maximum value, a new object is spawned surrounding this maximal critical point. As the value decreases through a saddle value, either distinct objects merge or portions of a single object merge, producing a change in topology. A vertex saddle is a point that is locally maximal (resp. minimal) with respect to two adjacent lattice vertices in a single direction and locally minimal (resp. maximal) to the four other adjacent vertices. A saddle point can also be embedded in a cubel face, in the case of the pattern producing 4 hit faces, and we will discuss this criticality below. The importance of these critical points can be demonstrated in a simple thought experiment that helps to make intuitively obvious the insight gained from Morse theory and its application to volume data.

Visualizing Morse Theory: The Fish Tank Story

A good intuative illustration of the application of (Stratified) Morse Theory to do a little thought experiment.

Let us consider a fish tank completely filled with water, and also containing two completely submerged smooth sand hills, one higher than the other. Connecting these two hills is submerged pass in between (two peaks and one pass)

Consider what happens as the water level is slowly lowered. As the water level is above the highest peak there are zero island objects in the scene. At the precise water height where the water level equals the highest peak, an interesting phenomena occurs. At this water level, we have a degenerate contour centered at the highest peak, but this island contour has zero area, zero perimeter, infinite curvature, and no gradient. Infinitesimally lowering the water causes this singularity to develop into a non singular contour ring island. There is now one island object in our scene. The peak is a local maxima, and from the peak, traveling in all directions a sufficiently short distance, all heights are less then or equal to the starting point. Should we be on a plateau, we must continue to travel until we see a change in height.

Figure 6: A completely submerged topology, where every point is below
the water level, or threshold.

Figure 7: An emerging topology, where the first object in the scene emerges.

Figure 8: A second emerging object. Between this value and the previous
threshold value, there is no topological changes, only geometric shape changes.

Figure 9: The two objects merge into one at a saddle point.

As the water level continues to drop, it will just touch or kiss the penultimate peak. Again it will produce a degenerate, singular contour before the point expands into a real island. There is now two island objects in the scene, as we can see in figure 8.

Continuing to lower the water level brings a pair of peninsulas closer and closer to each other till they just touch. This is shown in figure 9. The point where they touch is a point of non manifold contact where two contour rings touch. Most contouring algorithms never permit this osculating to occur. At this point, the gradient vanishes, and one point is shared by two sets of contours. Infinitesimally lower water levels cause the osculating point to expand into a path of finite width (an isthmus), and our two island objects merge into one object.

The final type of critical point can be represented by a lagoon, where decreasing water levels produce a shrinking lake. Again, we see as the water level just equals the nadir of the lagoon, we obtain a degenerate contour lake of zero perimeter, zero arena, no gradient, and infinite curvature. Similar to the definition of a peak, a local minima is a point at which traveling in all directions a sufficiently short distance will give a height that is greater to or equal to the starting point.

The basic insight from this thought experiment is that for any spectral sweep of tau varying from infinity to 0, each positive maxima (peak) adds one object to the scene, each minima (pit) removes an object from the scene, and each saddle point (pass) (locally) joins objects from the scene. The topological changes are completely described by the index of the criticality; the number of negative eigenvalues of the Hessian at the critical value. We are also guaranteed to have no topological changes for tau varied between any pair of critical values. There are some additional subtleties to this experiment, but the essential properties of any volume data set can be described by a hierarchy of these critical points.

These two facts are the fundamental result of Morse theory [Goresky88][Milnor63][Morse34]. The three dimensional density function defines a three dimensional manifold embedded in four dimensional space (the fourth axis is tau), and we are looking at the criticalities of the height function, i.e., projecting onto the tau axis. A projection onto a specific value of tau yields the isosurfaces for that particular threshold.

Toward Defining Critical Points: Axioms for Correct Isosurfaces

In order to justify our application and extension of Morse theoretic ideas to the digital domain, it is important to show that our assumptions are in fact shared (implicitely) by other practioners in the field, independent of the segmentation, object detection, or surface construction technique that they may employ. However they construct or classify their objects, these objects will be bounded by component isosurfaces that satisfy some basic properties.

First let us give the definitions of connected sets of lattice points or cubel vertices. These definitions are reminiscent of both graph connectivity and digital topology [Herman93].

Two vertices are edge adjacent if they share a cubel edge, and thus differ in one lattice coordinate. These are just the standard 4 and 6 adjacencies in two and three dimensional binary images or digital spaces. Path connectivity is the transitive closure of the edge adjacency relation. A set of cubel vertices is connected if for each pair of vertices in the set there is an (edge) connected path of vertices from the set, between the two vertices. We call this omega-connected since it is reminiscent of the similar notion of digital topology. If two vertices differ by at most one in each coordinate direction (but are not necessarily edge-adjacent) and the total difference is at most 2, we will say that they are alpha adjacent and define alpha-connectivity accordingly.

Recall that the set of lattice points with non-zero data value is bounded. Toward developing a digital Morse theory, let us make the following assumptions (or axioms) on a correct set of isosurfaces for a given volumetric data set and fixed isovalue tau.

Firstly we assume that the set of proposed isosurfaces for k-dimensional lattice data, at a non-critical isovalue, (at critical values, the isosurfaces are allowed to touch) are a collection of disjoint compact oriented k-1-manifolds.

Secondly we assume that there is one and only one surface intersection on a cubel edge (interior) if and only if the incident lattice vertex pair is labeled In and Out (respectively).

Thirdly we assume that the isosurfaces intersections with each k-dimensional cubel is a k-1 dimensional set that divides each k-dimensional cubel into a finite set of simply connected k-dimensional regions. More formally, the intersection of the complement of all the isosurfaces with any k dimensional cubel V is a finite collection of simply connected k-dimensional sets.

Fourthly we give the disambiguation rule for 4-hit faces. If a cubel has no interpolated saddle critical point, each such connected region of a cubel V contains a maximal omega-connected (within V) set of In vertices or Out vertices. If the cubel has an interpolated saddle critical point of taucritical and tau > taucritical then each region contains a maximal omega-connected set of In vertices or a maximal alpha-connected set of Out's (knit high, or knit In). If tau > taucritical then each region contains a maximal omega-connected set of Out vertices or a maximal alpha-connected set of In (knit low, or knit Out).

The justification for the first assumption comes from the basic goals of isosurface construction or digital imaging. The second assumption is simply an Occam's razor assumption that says that there should be an interface surface between any two adjacent data values that cross the isovalue threshold, and that without further information (more refined physical sampling) we assume a single crossing. The third assumption says that a correct set of surfaces should intersect the cubels in a Morse-like way and that we don't allow extraneous handles in the interior of cubels.

The fourth assumption embodies what is called the knit high (in) above the critical value and knit low below the critical value rule. This serves to disambiguate the 4-hit faces. Now this fourth axiom may seem kind of technical, but what is says it quite simple. If a cubel face has diagonally opposite pairs of In and Out readings then we disambiguate them, as the threshold value varies, in the following manner.

The critical value is the value, subject to the form of hit interpolation (linear in this case), below which the pairs of hits that share the low valued cubel vertices become closer than the pairs sharing high valued vertices. Above the critical value the diagonally opposed In vertices are part of locally (relative to the cubel) separate objects. At the critical value the two local objects merge. Below the critical value the In vertices are part of the same local object interior. Now these four assumptions serve to uniquely (up to homeomorphism) define a set of isosurfaces.

The fact is that one can take the fourth axiom in several different ways. In the next subsection and [KarronCoxMishra94] evidence is given to suggest why the above choice of axiom 4 is the ``right'' choice. That is, assuming that there is a single internal saddle critical point in this ambiguous cubel face, axiom 4 leads to the simplest and most reasonable set of local topological transitions as the threshold is varied.

Examples of Knitting Rules for Axiom 4

Here we compare knitting rules. The first, shown in figure 10, is ``Knit High above taucrit, Knit Low below'', and the second, shown in figure 11, is ``Knit Low above taucrit, Knit High below''. We can see that the first rule is appropriate, since it corresponds to the local topological changes that one would expect if there were only one saddle critical point. By contrast, we get a number of extra local transitions in figure 11 . In particular, in 4 and 8, we get splitting and merging. There are similar problems with other rule choices.

Figure 10: ``knit high above knit low below'' rule. Note that there is only one
discontinuity, at the critical point, at 6, which is as expected by Digital Morse Theory.

Figure 11: ``knit low above knit high below'' rule. Note that there are three
discontinuous transitions, from 3 to 4, at 6, and from 8 to 9. Using the wrong
knitting rule identified incorrect criticalities.

Identifying Criticalities Required for Building Graph

Regardless of the choice of knitting rule a consistent segmentation will imply a critical point in this pattern of readings. That is as we vary the isovalue through the critical value almost all practitioners regardless of preferred technique, will produce images whose topology changes (for example two objects merge) as the isovalue falls below the critical value. However, the choice of a knitting rule is critical for correct identification of criticalties for creating the Morse decomposition to give the Reeb Graph and Object Spectra Graph.

Applications of the Theory to (Simple) Volume Data

Identifying Critical Points and Building the Object Spectra

One first identifies all locally maximal lattice points or connected sets of locally maximal critical points and records the critical values. Then one identifies the minimal critical points. Finally one identifies the saddle points, including those implied by the 4 hit faces. One divides the range of function values into intervals bounded by these critical values. One also stores information on the associated critical points at each interval boundary (including the Morse data).

Once one has completely characterized all of the possible Digital Morse criticalities, we can now use the natural global relationship between criticalities to build the object spectrum (with the associated Morse data). One traverses the range of values, by intervals, starting with the largest value and records the number of objects as each is spawned and the topological changes at each saddle point extremal point.

An object spectrum, such as indicated in figure 12, can provide a guide to what and where objects are at any threshold. Additionally one plots a sample isosurface for a sample isovalue in each critical value-free range. By plotting at each saddle and minimum one can respectively determine if two objects merge (resp. split) or a single object increases (resp. decreases) genus in the former case and whether an object vanishes (resp. appears) or a hole (or handle) vanishes (or appears) in the latter case. Each criticality can trace a lineage back to a parent criticality, or groups of criticalities. This provides a powerful way to organize features in a particular geography.

We trace this ancestory of an objects associated with a criticality, since new objects can result from the splitting of an object at a saddle point, and at a threshold above the largest data value we have no objects (all lattice points are exterior). Further we can run simultaneous best-first scanning (grassfire) algorithms from the criticalities to label each region of influence. Each peak or pit can be considered to ``own'' the region surrounding it, bounded by the largest connected set of lattice points that contain no saddle point, and are not owned by another criticality. This is well known in the theory of dynamical systems, where the owned region is known as the critical point's ``basin of attraction". The grassfire algorithm labels each point according to the critical point that owns it. For example we begin a wavefront propogation from each local maximum, and label in stages all lattice points of lesser value than the maximum, the waves propogating according to decreasing isovalue. In particular we place neighbors on a priority queue (max heap) for continued evaluation. If there is an interval bounded by critical values 60 and 50, with corresponding critical points p and q, the labeling would propogate

from p and label a contiguous set of lattice vertices of value greater than 50, in stages. When the value of 50 is reached a wavefront would propogate from q in stages, again by decreasing values until the wavefronts meet at a saddle point. This type of algorithm has been used in discrete potential function calculation in robotics, for example by Dr. Latombe[Latombe91]. This results in a nice segmentation of the data by critical points.

After we partition the data space by this criticality parentage, the relative importance of a criticality, or group of criticalities can be measured by the area and perimeter of the regions it controls (the size of the basin of attraction). We can eliminate or ignore small criticalities, and choose to build hierarchically more detailed models in this theory. Small noisy features, only a few boxels in extant, can be eliminated, and the regions condensed into broad, monotonic rectlinear regions, called super-cubels. A similar super cubel concept was developed by Schroeder[Schroeder92] and Kalvin[Kalvin93], but was applied to a surface, and not applied to a volume. The SpiderWeb can knit triangles spanning many cubels on such a region. These triangles will preserve the topology of the region, but will only approximate the curvature in this region.

Figure 12: For every valid isovalue, tau, we can look up what objects exist
there, and how many distinct objects are present. This can be an important
guide to apprehending a volume. This is a kind of Reeb Graph.

Figure 13: Monotonic ``Super Boxel'' hits behave the same way as regular, simple
boxels, with only 2 hits possible along entire perimeter, the dashed line.

Implications for Data Space Navigation

Using the above ideas one can approximately render a region (using super cubels), and then render (locally) increasing detail as needed. Thus the interior of the super cubel is retiled only as needed. Again, it should be emphasized that super cubels are constructed independently of the particular isovalue.

One can also navigate in threshold space along the tau-axis. For example as the isovalue tau is varied in a critical point-free interval the triangulated SpiderWeb surfaces continually changes as the hit points defining the triangles slide continuously along the lattice (cubel) edges. Since there is no essential topological change, we will develop efficient ways to grow and change our representation as the isovalue is varied. The principal investigator has already implemented this ability for two dimensional datasets.

This presents exciting possibilities for the ``virtual reality" visual navigation of volumetric data. In addition, by completely characterizing all possible topologies that exist in the dataset (over all threshold values) one obtains a powerful tool for visual pattern recognition, e.g., in atmospheric or geophysical data. This will improve greatly on ``hunt and peck" techniques applied to segmentation (selection of the appropriate isovalue range to reveal the structure desired) and aid in application of the more analytic techniques, as well as the various object recognition techniques applied from computer vision.

In some cases it should be apparent that one is not looking for objects as such, but isovalued structures not necessarily bounding regions of smoothly varying greyscale intensity. In these cases visual object recognition techniques (gradient edge detection, for example) may not really apply. DMT data space analysis and navigation would really show its unique utility in these cases. The investigators believe that there may be applications for these tools in many areas.

Visual Navigation of Vector
Valued Volume Data

It now remains to extend the ideas to allow the data space navigation of higher dimensional data, using the visual human dataset as our for example. For fixed linear combination of the RGB value at each point we can identify the criticalities and navigate the dataset as above. The criticalities vary as we change the coefficients of the combination, that is, move within RGB space. That is we can define a scalar function for each convex combination of coefficients. We can regard ourselves as navigating within this convex set as moving within RGB space. But again since each point in RGB space represents a 4D slice through a five dimensional manifold in six dimensional space (that is yields a single volume data set as defined above) the insights of Morse theory still apply. In particular there will be regions (cells) of RGB space where the set of Morse criticalities in the projected volume data set remain essentially the same, that is, yield the same or topological equivalent object spectra. Now the goal is to graph this higher dimensional object spectra, as a family of object spectra, with associated RGB regions. This set of object spectra can be represented by a graph, with properties similar to the object spectral graph (which is itself a generalization of a Reeb graph) for a simple volume data set. Dividing RGB spaces into these cells will have to be accomplished by extending the techniques used in cylyndrical algebraic decomposition to the digital (volumetric data) domian. We are currently investigating this idea and others for dealing with data representing vector fields.

We contend that this higher dimensional extension of Digital Morse theory represents a ``right'' way to apprehend complex multidimesional volume data, and the true way to apply the insights of Morse Theory and variational calculus to the data analysis and visualization problem.

Research Objectives

Digital Morse Theory (DMT) will provide the ability to completely apprehend volume data lattices independent of any particular isovalue. In particular, this gives us the ability to identify, up to topological equivalence, all possible isosurfaces that can be grown for any threshold value.

Develop Digital Morse Theory

The investigation will continue to develop the theoretical foundations of Digital Morse Theory, contributing to the fundamental information theoretic understanding of volumetric data.

Develop New Algorithms

The investigation will apply this theory in the continued implementation of software systems, which will be made available to the researcher community via Internet. The investigation proposes to develop:

Develop New Software

Develop Applications in Other Visualization
Applications Other Than Medical

The investigation will seek to apply these techniques in other applications areas, emphasizing scientific visualization and visual navigation of data.

Statement of Work

Two Phases

Our work will be divided into two phases. It is expected that we will have available to this project an SGI ONYX multiprocessing supercomputer funded by our other to be funded (by Rick Satava) DARPA virtual reality project.

Combined Pure Theory and Practical Software Effort

This effort will consist of a combined computer science theory effort and a software effort, with no resources allocated to hardware. We are asking for long term funding (5 years aggregate), contingent on promising results and of the Phase 1 deliverables, of the two investigators to do the programming work.

Qualifications of Investigators

Both Karron and Cox are world class computer programmers. Karron has years of experience in computer graphics, and Cox has years of experience in designing and programming basic computer algorithms for academic as well as commercial applications.

Phase 1 Description Summary

Phase 1 lasts for two years and will consist of writing C++, SGI Performer (or COSMO), Tk/Tcl code to compile the OSG, interactively and selectively display portions of the OSG, and a feedback POV/LOD limiviewer to display salient features in the viewing frustrum, and cull features outside our POV and to recursively draw those features with only as much detail as the screen pixilation warrents. Our sample dataset will consist of the visible human dataset, where the RGB pixel data will be projected into a grey scale scalar black and white image.

Phase 2 Description Summary

Phase 2 will build on our experience using DMT OSG in three dimensions and we will extend our OSG compiler, display, and POV/LOD to multidimensional / multispecral data sets. We will explore the full diminsionality of the visible human dataset in RGB spac.

This is expected to take 3 years and depends on sucessful implementation of the core code classes in three dimensions in the phase 1. Specific work items are detailed below.

Phase 1 Statement of Work

Individual aspects of the construction and refinement of the Phase 1 OSG Compiler are as follows:

Specific work tasks in the Phase 1 viewer are development of a feedback viewer that can navigate in

Phase 2 Statement of Work

Phase 2 is an optional follow on project that depends on our demonstration of mastery of the Phase 1 deliverables. We will continue to refine our software techniques with new methods and classes to extend our visualization to surfaces and criticalites in 6 dimensions. We will use the additional dimensions of RGB space.

We will develop each of the above tasks in the Phase 2 to include additional dimensions. Our overall goal is to develop code to visualize 6 dimensional criticalties as projections of RGB space in the visible human dataset.

Deliverables

Development and progress on each individual deliverable will be simultaneous and incrimental for each bulleted item.

Common Phase 1 and 2 Deliverables

In general, Phase 1 and Phase 2 deliverables will be:

Phase 1 Deliverables

Phase 2 Deliverable

Our phase 2 deliverables will be essentially the same the the phase 1 deliverables, only extended into three additional dimensions of RGB space.

Time of Delivery of Deliverables

As this is a research effort, precise achievement of each deliverable will be incrimental, over the life of the project. We will incrimentally improve each component deliverable from the first buggy prototype to published public methods. We expect a significant period of public release of the code, and use by various other researchers, and feedback as to bugs, desired new features, and deletion of unnecessary features. Each deliverable item will evolve and will improve over the entire life of the project.

References

[Artzy81] E. Artzy, G. Frieder, and G. T. Herman. The theory, design, implementation and evaluation of a three-dimensional surface detection algorithm. Computer graphics and Image Processing, 15:1-24, January 1981.

[Baker88] H. Harlyn Baker. Building, Visualizing, and Computing on Surfaces of Evolution. IEEE Computer Graphics and Applications, 8(4):31-41, July 1988.

[Bloomenthal88] Jules Bloomenthal. Polygonization of Implicit Surfaces. Computer Aided Geometric Design, 5:341-355, 1988.

[ChenCoxMishra91] Jianer Chen, James Cox, and Bud Mishra. NL Heirarchy. Information Processing Letters, 39, 1991.

[Cline87] Harvey E. Cline, C. L. Dumoulin, H. R. Hart Jr., and William E. Lorensen. 3D Reconstruction of the Brain from Magnetic Resonance Images Using a Connectivity Algorithm. Magnetic Resonance Imaging, 5:345-352, July 1987.

[Cline87a] Harvey E. Cline and William E. Lorensen. System and method for the display of surface structures contained within the interior region of a solid body. United States Patent Number 4,710,876, 1987.

[Cline88] Harvey. E. Cline, William. E. Lorensen, S. Ludke, and Bruce. C. Teeter C. R. Crawford. Two Algorithms for the Reconstruction of Surfaces from Tomograms. Medical Physics, 15(3):320-327, May 1988.

[Cox89] James Cox. On-Line Motion Planning. Ph.D. dissertation, New York University, Courant Institute, Department of Computer Science, February 1989.

[CoxKarron93] James L. Cox, Daniel B. Karron, and B. Mishra. The SpiderWeb Algorithm for Surface Construction from Medical Volume Data: Geometric Properties of its Surface. Innovations Et Technologie en Biologie et Medecine, 14(6):634-656, November 1993.

[CoxMcAloon93] James Cox and Ken McAloon. Decision Procedures for Constraint-based Extensions of Datalog. Constraint Logic Programming: Selected Research. MIT Press, Cambridge, Massachusetts, 1993.

[CoxMcAloonTretkoff90] James Cox, Ken McAloon, and Carol Tretkoff. Computational Complexity and Constraint Logic Programming Languages : Extended Abstract, chapter 2. Logic Programming: Proceedings of the 1990 North Americal Conference. MIT Press, Cambridge, Massachusetts, 1990.

[CoxMcAloonTretkoff92] James Cox, Ken McAloon, and Carol Tretkoff. Computational complexity and constraint logic programming languages. Annals of Mathematics and Artificial Intelligence, 5, 1992.

[CoxMishraEricson96] James Cox and Lars Ericson and Bud Mishra. The Average Case Complexity of Multilevel Syllogistic. Communications in Pure and Applied Mathematics, 1996. Accepted for publication subject to revisions.

[CoxWang92] James Cox and Xiaozhaung Wang. Adaptive potential function calculation for sensory based robot navigation. Research report, Brooklyn College of City University of New York, Computer Science Department, October 1992.

[CoxYap91] James Cox and Chee Yap. On-Line Motion Planning: Case of a Planar Rod. Annals of Mathematics and Artifical Intelligence, 1(3), 1991.

[Frey91] Pascal Frey and Michel Gautherie. Generation Automatique D'un Maillage 3D Dans Un Ensemble De Voxels: Application à la modèlisation thermique numèrique 3D en thermothèrapie ultrasonore. Innovation Et Technologie En Biologie Et Medecine, 12(4):428-442, April 1991.

[Gauch92] John M. Gauch. Multiresolution Image Shape Description. Springer Series In Perception Engineering. Springer, New York, 1992.

[Goresky88] Mark Goresky and Robert MacPherson. Stratified Morse Theory, volume 14 of Series 3, A Series of Modern Surveys In Mathematics. Springer-Verlag, Berlin, 1988.

[Herman93] Gabor T. Herman. Oriented Surfaces in Digital Spaces. CVGIP: Graphical Models and Image Processing, 55(5):381-396, September 1993.

[Hoffmann89] Christoph Martin Hoffmann. Geometric and Solid Modeling: An Introduction. Morgan Kaufmann Publishers, Inc., San Mateo, California, 94403, 1989.

[Itoh96] Takayuki Itoh and Koji Koyamada. Automatic Isosurface Propagation Using an Extrema Graph and Sorted Boundary Cell Lists. IEEE Transactions on Visualization and Computer Graphics, 1(4):319 - 327, December 1995.

[Kalvin91] Alan D. Kalvin. Segmentation and Surface-Based Modeling of Objects in Three-Dimensional Biomedical Images. Ph.D. dissertation, New York University, Computer Science Department, 1991.

[Kalvin91a] Alan D. Kalvin, Court B. Cutting, Betsy Haddad, and Marlyn E. Noz. Constructing Topologically Connected Surfaces for the Comprehensive Analysis of 3D Medical Structures. In SPIE Medical Imaging V: Image Processing, volume 1445, pages 247-258. Society for Photo-Optical Instrumentation Engineers, 1991.

[Kalvin93] Alan D. Kalvin and Russell H. Taylor. SuperFaces: Polyhedral Approximation with Bounded Error. Research Report RC19135 (82286), IBM Research Division, IBM Research Division, T. J. Watson Research Center, Yorktown Heights, New York 10598, April 1993.

[Karron92ParisPoster] Daniel Karron and James Cox. Mathematical Analysis of ``SpiderWeb'' Surface Construction Algorithm]. In J. P. Morucci, R. Plonsey, J. L. Coatrieux, and S. Laxminarayan, editors, Proceedings of the 14th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, volume 14, pages 1250-1252, Piscataway, NJ, 1992. IEEE Engineering in Medicine and Biology Society.

[KarronCox95] Daniel Karron and James Cox. The spiderweb algorithm for extracting 3D objects from volume data. In Nicholas Ayache, editor, Computer Vision, Virtual Reality and Robotics in Medicine, number 905 in Lecture Notes in Computer Science: Proceedings of the First International Conference, CVRMed'95, pages 481-486, Berlin, 1995. INRIA-Sophia Antipolis, Springer-Verlag.

[KarronCoxMishra94] Daniel B. Karron, James Cox, and Bhubaneswar Mishra. New findings from the spiderweb algorithm: Toward a digital morse theory.In Richard A. Robb, editor, Visualization in Biomedical Computing - '94, number 2359 in SPIE Proceedings, pages 643-657. SPIE, SPIE, October 1994.

[Kergosien92] Yannick L. Kergosien. Topology and Visualization: From Generic Singularities to Combinatorial Shape Modelling. In Tosiyasu L. Kunii and Yoshihisa Shinagawa, editor, Modern Geometric Computing for Visualization, pages 31-54. Springer, 1992.

[Latombe91] J. C. Latombe. Robot Motion Planning. Kluwer Publishers, 1991.

[Lorensen87] William E. Lorensen and Harvey E. Cline. Marching Cubes: A High Resolution 3D Surface Construction Algorithm. ACM Computer Graphics, 21(4):163--169, July 1987.

[Milnor63] J. Milnor, M. Spivak, and R. Wells. Morse Theory, volume 51 of Annals of Mathematics Studies. Princeton University Press, Princeton, New Jersey, 5th edition, 1963.

[Mishra91] Bud Mishra. Notes on Solid Modeling. Lecture notes, New York University, Computer Science Department, 1991.

[Morse34] Marston. Morse. The calculus of variations in the large, volume 18 of Colloquium Publications Series. American Mathematical Association, Providence, R.I., 8th printing edition, 1934.

[Nagae93] Takanori Nagae, Takeshi Agui, and Hiroshi Nagahashi. Orientable Closed Surface Construction from Volume Data. IEICE Transaction Inf. and Systems, E76D(2):269-273, February 1993.

[Nielson91] Gregory M. Nielson and Bernd Hamann. The Asymptotic Decider: Resolving the Ambiguity in Marching Cubes. In Gregory M. Nielson and Larry Rosenblum, editors, Proceedings of Visualization 91, volume 2, pages 83-91, 1991.

[Ning93] Paul Ning and Jules Bloomenthal. An Evaluation of Implicit Surface Tilers. IEEE Computer Graphics and Applications, 13(6):33-41, November 1993.

[Schroeder92] William J. Schroeder, Jonathan A. Zarge, and William E. Lorensen. Decimation of triangle meshes. ACM Computer Graphics, 26(2):65-70, July 1992.

[Shinagawa91a] Yoshihisa Shinagawa, Tosiyasu L. Kunii, and Yannick L. Kergosien. Surface Coding Based on Morse Theory. IEEE Computer Graphics and Applications, 11(5):66-78, September 1991.

[Shinagawa91b] Yoshihisa Shinagawa and Tosiyasu L. Kunii. Constructing a Reeb Graph Automatically from Cross Sections. IEEE Computer Graphics and Applications, 11(6):44-51, November 1991.

[Shinagawa92] Yoshihisa Shinagawa and Tosiyasu L. Kunii. Using Surface Coding to Detect Errors in Surface Reconstruction. In Tosiyasu L. Kunii and Yoshihisa Shinagawa], editor, Modern Geometric Computing for Visualization, pages 227 - 240. Springer, 1992.

[Udupa82] Jayaram K. Udupa, Sargur N. Srihari, and Gabor T. Herman. Boundary Detection in Multidimensions. IEEE Transactions on Pattern Analysis and Machine Intelligence, 4(1):41-50, January 1988.

[Udupa91] Jayaram K. Udupa, Hsiu-Mei Hung, and Keh-Shih Chuang. Surface and Volume Rendering in Three-Dimensional Imaging: A Comparison. Journal of Digital Imaging, 4(3):159-168, August 1991.

[Udupa92] Jayaram K. Udupa and Gabor T. Herman. Boundaries in Multidimensional Digital Scenes: Theory and Algorithms. Technical Report MIPG182, Medical Image Processing Group, University of Pennsylvania, Medical Image Processing Group, Department of Radiology, University Of Pennsylvania, Brockley Hall, Fourth Floor, 418 Service Drive, Philadelphia, Pennsylvania, 19104-6021, January 1992.

[VanGelder90] Allan Van Gelder and Jane Wilhelms. Topological Ambiguities in Isosurface Generation. Technical report, University of California at Santa Cruz, 1990.

[Wallin91] Ake Wallin. Constructing Isosurfaces from CT Data. IEEE Computer Graphics and Applications, 11(6):28-33, November 1991.

[Wilhelms90] Jane Wilhelms and Allan Van Gelder. Topological considerations in isosurface generation: Extended abstract. ACM Computer Graphics, 24(5):79-86, November 1990.