cgg_vijay

Valerio Pascucci
Lawrence Livermore National Laboratory, California, USA.
pascucci1@llnl.gov

The 3D Morse-Smale complex is a fundamental topological construct that partitions the domain of a real-valued function into regions having uniform gradient flow behavior. In this paper, we consider the construction and selective presentation of cells of the Morse-Smale complex and their use in the analysis and visualization of scientific datasets. We take advantage of the fact that cells of different dimension often characterize different types of features present in the data. For example, critical points pinpoint changes in topology by showing where components of the level sets are created, destroyed or modified in genus. Edges of the Morse-Smale complex extract filament-like features that are not explicitly modeled in the original data. Interactive selection and rendering of portions of the Morse-Smale complex introduces fundamental data management challenges due to the unstructured nature of the complex even for structured inputs. We describe a data structure that stores the Morse-Smale complex and allows efficient selective traversal of regions of interest. Finally, we illustrate the practical use of this approach by applying it to cryo-electron microscopy data of protein molecules.

Keywords: Morse theory, Morse-Smale complexes, 3D scalar field visualization, computational geometry, computational topology, data structures.

Topological methods are being used increasingly for the exploration and analysis of scalar fields. Many of these methods are based on Morse theory, which studies the relationship between the shape of a space and scalar functions defined over it. One such method uses Morse-Smale complexes, which partitions the domain space by clustering regions having uniform gradient flow. Efficient algorithms for computing the Morse-Smale complex in 2D and 3D have been proposed recently [10,12]. These algorithms employ combinatorial methods and are hence stable, in contrast to numerical methods that perform integration to trace the flow lines. Besides providing a visually compact presentation of the scalar function, the Morse-Smale complex also allows us to abstract features in the data. This abstraction has been used to automatically identify important regions, i.e. features, in 2D scalar data and preserve them in a multiresolution representation [4]. Extending this to 3D is made nontrivial by the presence of new types of critical points.

In this paper, we describe different ways to use substructures of the Morse-Smale complex for analyzing 3D scalar fields. These substructures give the skeletal shape of isosurfaces and are therefore useful both for visualization as well as for further analysis of the shape. We illustrate these ideas using cryo-electron microscope images of two biological macromolecules. The 1-dimensional cells of the Morse-Smale complex extract the ring-like shape of the macromolecules. Since storage of the complex is nontrivial, we first describe a data structure that stores the Morse-Smale complex before discussing how individual substructures can be visualized.

The family of isosurfaces of the scalar function sweep out each cell in the Morse-Smale complex starting at a minimum and proceeding toward the maximum while crossing the boundary orthogonally. The Reeb graph is another structure that encodes the evolution of the topology of isosurfaces as they sweep the domain space [20]. Efficient methods have been developed to compute the Reeb graph [6] with special emphasis on the case when it is a tree, namely the contour tree [5,19]. The Reeb graph gives a compressed representation of the isosurface components and can therefore be used for modeling the shape of objects [14,23]. However, it has no geometric information as expressed in the Morse-Smale complex and hence not directly usable for representing geometric features.

Another concept related is the medial axis of a shape in three-dimensional Euclidean space. As introduced by Blum to describe the shape of biological structures, it is the set of centers of spheres that touch the boundary of the shape in at least two points without crossing it [2]. Medial axes are used in a wide variety of applications besides shape representation. The medial axis of isosurfaces does contain geometric information as opposed to Reeb graphs. However, the medial axis is sensitive to small changes in the geometry of the isosurface resulting in possibly big differences between the representation of two topologically equivalent isosurfaces. Working with a medial axis representation can therefore be time consuming. On the other hand, a simple threshold based filtering applied on substructures of the Morse-Smale complex leads to an efficient representation of these features. There is an interesting relationship between medial axis and the Morse-Smale complex, which we do not address further in this paper: if the boundary is an orientable 2-manifold embedded in 3D Euclidean space, we may define the signed distance as a function over the space. The medial axis then consists of 1D and 2D cells in the Morse-Smale complex.

We apply Morse theoretic ideas to study scalar data, usually available as sampled functions. Morse theory was, however, originally developed for smooth functions and therefore, as a first step, we extend it to the piecewise-linear domain. In this section, we give a brief description of this extension while reviewing Morse functions, critical points, lines of steepest ascent and Morse-Smale complexes for three-dimensional manifolds. The paper by Edelsbrunner et al. [12] gives a more complete description of how to transport Morse theoretic ideas from smooth to the piecewise-linear domain in order to define Morse-Smale complexes. We refer the reader to the books by Matsumoto [15] and Milnor [16] for further background on Morse theory.

Smooth manifolds. Let $\mathbb{M}$ be a smooth compact 3-manifold without boundary. A simple example is the unit 3-sphere, which consists of all points at unit distance from the origin in $\mathbb{R}^4$ . A smooth map $f: \mathbb{M}\rightarrow \mathbb{R}$ is a Morse function if none of its critical points are degenerate (i.e. the Hessian matrix at all critical points is non-singular) and no two critical points have the same function value. The Morse Lemma states that a Morse function has quadratic behavior within a local neighborhood of every non-degenerate critical point

. This immediately gives a characterization of the critical points. Figure 1 shows local pictures for the four types of non-degenerate critical points. The criticality of

is characterized by the structure of the oceans, consisting of points

on the sphere around

with

, and continents, consisting of points

on the sphere with

**Fig. 1.** Local pictures, with shaded oceans and white continents, of a regular point and the four types of critical points: minimum, 1-saddle, 2-saddle, and maximum.
$\begin{figure}\vspace*{0.1in} \centering \includegraphics[height=1.2in]{Figs/spheres.eps} \end{figure}$

An integral line of

is a maximal path in $\mathbb{M}$ whose tangent vectors agree with the gradient of

at every point of the path. Each integral line has a natural origin and destination at critical points of

where the gradient becomes zero. The Morse-Smale complex partitions $\mathbb{M}$ into regions by clustering integral lines based on their origin and destination. For example, the 3-dimensional cells of the Morse-Smale complex cluster integral lines that originate at a given minimum and terminate at an associated maximum. The cells of different dimensions are called crystals, quads, arcs, and nodes. Ascending and descending manifolds are obtained as clusters of integral lines having common origin and destination respectively. The Morse-Smale complex can also be viewed as an overlay of the ascending and descending manifolds. Note that the ascending manifolds of

are the descending manifolds of

and vice versa. This duality between the two types of manifolds implies that they have the same structural properties.

Piecewise-linear manifolds. Smooth manifolds are represented by triangulations in the PL setting. Let K be a triangulation of the given 3-manifold $\mathbb{M}$ . The underlying space of K is homeomorphic to $\mathbb{M}$ . K consists of simplices of dimension 0 to 3, which we refer to as vertices, edges, triangles and tetrahedra. The spherical neighborhood of a vertex

is represented by a triangulation of the vertex neighbors of

. Edges and triangles in this triangulation, called the link of

, are exactly the faces of tetrahedra that contain

Let $f: \mathbb{M}\rightarrow \mathbb{R}$ be a continuous function that is linear within each simplex of K. Assuming

is given at the vertices, we can compute it over a simplex by linear interpolation. Oceans consist of simplices in the link where

takes values lower than at

. We count the number of components and tunnels in the ocean, using their reduced Betti numbers, to distinguish regular from critical vertices and to classify the latter. Gradients and hence integral lines are not well defined for the PL domain. However, monotonic curves and surfaces corresponding to the arcs and quads of the Morse-Smale complex can be constructed using the idea of simulated disjointness. When curves and surfaces merge, an infinitesimal separation is maintained between them in order to avoid any illegal crossings. This separation is purely symbolic and no geometric perturbation is actually applied [12]. The function

has its critical points at the nodes of this complex and is monotonic within all the arcs, quadrangles and crystals.

The Morse-Smale complex is constructed during two sweeps of the 3-manifold. The first sweep is in order of decreasing values of the function

and computes the descending manifolds. The ascending manifolds are computed during the second sweep, which is in the order of increasing function value. The structure of the descending manifolds is used while computing the ascending manifolds in order to maintain structural integrity of the Morse-Smale complex. We require special data structures for storing these manifolds and for resolving various kinds of degeneracies that arise during the construction. In an effort to achieve modularity, we store the topology and geometry of the Morse-Smale complex separately. We maintain two data structures, one for the triangulation K of the 3-manifold and the other for the Morse-Smale complex Q of the function

. The two consist of various pieces and are connected to each other as shown in Figure 2.

**Fig. 2.** The data structures used to represent the triangulation K and the Morse-Smale complex Q.
$\begin{figure}\vspace*{0.1in} \centering \includegraphics[height=2.8in]{Figs/datastructure.eps} \end{figure}$

Triangulation. The triangulation of $\mathbb{M}$ is a 3-dimensional simplicial complex consisting of vertices, edges, triangles, and tetrahedra. Miscellaneous information about the simplices is stored in the arrays

, and

. Each simplex is identified by its dimension and its index in the corresponding array. The connectivity between the simplices is represented by another array, K, whose elements are sextets of pointers that connect the triangles in rings about the edges. This data structure is akin to the edge-facet structure introduced in [7]. As illustrated in Figure 3, a sextet is made up of the six ordered versions of a triangle.

**Fig. 3.** Algebra of ordered triangles.
$\begin{figure}\includegraphics[height=1.5in]{Figs/triangles.eps} \end{figure}$

To illustrate the functionality of this data structure, consider the computation of the link of a vertex $p = p_{i}$ . Letting

be one of the triangles that share that vertex, we use depth-first search to traverse all triangles in the star. For each visited triangle

, the edge

belongs to the link of

and so do the triangles that precede and succeed

in the ring around

. Given the initial triangle

, the search takes time proportional to the number of edges in the link. With an additional test of the vertex heights, we can identify the lower link as a subcomplex of the link.

Morse-Smale Complex. The Morse-Smale complex of

consists of simple open cells of dimensions 0, 1, 2, and 3, which we refer to as nodes, arcs, quadrangles, and crystals. Miscellaneous related information is stored in arrays

, and

, as shown in Figure 2. Each cell is identified by its dimension and its index in the corresponding array. The connectivity is again stored in another array, Q, whose elements are octets of pointers, as illustrated in Figure 4.

**Fig. 4.** Algebra of ordered quadrangles.
$\includegraphics[height=1.5in]{Figs/quadrangles.eps}$

A node in the Morse-Smale complex is just one of the vertices of the triangulation, but arcs and quadrangles are more complicated objects that need elaborate representations. Each arc is an open sequence of edges and, as indicated in Figure 2, it is stored as a sequence of oriented edges or pointers into K. Similarly, each quadrangle is an open patch of triangles and, as indicated in Figure 2, it is stored as a 2-manifold of oriented triangles or pointers into

. The connectivity can be recovered from the connectivity information in K, so we can get by with simple and compact representations. It is important that each arc and quadrangle has its own fixed orientation, which is then used in the algebra illustrated in Figure 4.

Normal structures. It is possible that multiple quadrangles merge and pass through the same triangle and multiple arcs pass through a common edge of K. Furthermore, a quadrangle in the Morse-Smale complex could have boundary arcs glued to each other. In order to maintain connectivity information in these cases, we need some additional data structures. For each edge in K, we store a snapshot of the Morse-Smale complex restricted to its neighborhood. This snapshot, called the normal disk of the edge is needed to maintain the order of arcs and quads passing through the edge and the triangles containing it.

The Morse-Smale complex has a rich structure because it encodes topological as well as geometric features of the data. Visualizing it presents the same problems as in the display of the raw data, namely that of visual clutter. However, substructures of the Morse-Smale complex can be displayed individually and these provide specific information about the data. In this section, we describe these substructures and illustrate their use by constructing and visualizing them for a particular scientific dataset. We begin with a description of the data.

Data. Cryo-electron microscopy is a technique used to determine the structure of biological molecules. The biological material is suspended in vitreous ice for imaging. The vitrified samples are transferred to an electron microscope where they are imaged using low electron doses. The cryo-EM datasets are available as density maps over a 3D space from the macromolecular structure database maintained by the European Bioinformatics Institute [17]. We use two datasets: the first one is a DNA helicase called dnaB and the second is that of a complex ¹ between dnaB and another protein called dnaC .

Ascending/Descending arcs. The ascending arcs are piecewise-linear curves between 2-saddles and maxima lying along edges of the input mesh. So, they can be displayed as a collection of edges. Since it is difficult to perceive depth in this display, we render the arcs as thin tubes. In many datasets, there is typically a region of interest, either spatial or in function space. We provide some filters that allow the user to select the arcs to be displayed. For example, the user can clip the arcs to those that originate from saddles with function values in a given range. A threshold arc length can be chosen to clip away short arcs that typically represent less significant topological features. The user can also prune away the arcs that span a short interval in the function space. This is akin to smoothing of the function to remove minor perturbations.

**Fig. 5.** The enzyme DNA helicase breaks hydrogen bonds to unravel the DNA double helix during the replication process. `dnaB` is the major replicative helicase in E-coli and is known to have a ring-shaped hexamer structure. One face of this ring has threefold symmetry and the other has sixfold symmetry. `dnaC` transports `dnaB` to the site of action by binding with it to form a complex. The shape of the `dnaB-dnaC` complex is similar to that of `dnaB` , with a tunnel at the center. Each one of the six subunits of `dnaC` interacts with two subunits of `dnaB` and vice-versa. *Top*: A translucent isosurface extracted from the `dnaB` dataset (left) and from the `dnaB-dnaC` complex dataset (right). Ascending arcs clipped to values above the corresponding isovalue have been overlaid on the isosurface. The ascending arcs trace the ring-like shape of the `dnaB` molecule. The `dnaC` molecule attaches itself to the face of `dnaB` that has sixfold symmetry. The two isosurfaces are different but the tunnel still remains as can be seen in the ascending arcs. *Bottom*: Ascending arcs trace the skeletal shape of isosurfaces for `dnaB` . Note how the loops in the ascending arcs break along with the isosurface.
$\includegraphics[height=2in]{Figs/dnab.ps}$ $\includegraphics[height=2.2in]{Figs/dnabc.ps}$ $\includegraphics[height=2.7in]{Figs/dnabsixholes.ps}$

We use the ascending arcs and isosurfaces to visualize the structure of dnaB and the dnaB-dnaC complex (see Figure 5). The higher density regions in the data correspond to the biological molecule whereas the lower density regions represent the solution. So, larger isovalues are more significant and their corresponding isosurfaces give the 3D structure of the molecule. Each component of the isosurface encloses at least one local maximum. We filter the ascending arcs using the given isovalue as a threshold. These filtered arcs terminate at the enclosed maxima and trace the skeletal shape of the isosurface. The ring like structure is persistent over a sizeable range of density values for both cryo-EM datasets. The corresponding isosurfaces do indeed represent the structure of the molecules. Ascending arcs can be ordered based on the 2-saddle where they originate. Therefore, changing the threshold corresponds to simply choosing the 2-saddles that lie in the chosen range. We found it useful to work with this skeletal shape while exploring the data because the ascending arcs can be computed and rendered much faster than the isosurfaces. The skeleton shape of isosurfaces given by ascending arcs have a natural geometric embedding as opposed to Reeb graphs. In the case of the latter, the nodes of the graph correspond to vertices of the triangulation. However, there is no natural way of embedding the edges of the graph to lie along edges of the triangulation.

Descending arcs behave similarly in other datasets where local minima capture important features. We envision the use of these 1-manifolds for the exploration of volumetric data. This complements existing approaches like the contour spectrum [3]. For example, the skeletal shape can be generated for isovalue ranges determined from the contour spectrum.

Ascending/Descending disks. The ascending and descending disks are piecewise-linear surfaces passing through the triangles of the given mesh. The boundary of each ascending/descending disk consists of a collection of ascending arcs and their corresponding nodes. The ascending disks can be used to annotate individual isosurfaces as well. The intersection of these disks with an isosurface is computed efficiently as an isocontour lying within the disk. Ascending disks can be extracted from the Morse-Smale complex as a collection of quadrangles that contain the source 1-saddle and a maximum as two nodes. Given our data structure that allows efficient traversal from a quadrangle to neighboring quadrangles, this extraction can be performed efficiently.

**Fig. 6.** Visualizations of a function defined as the minimum distance from a set of three points , , and within a unit cube. Assume a minimum at infinity connected to all points on the boundary. The critical points are displayed as small spheres in blue, cyan, green, and red for minima, 1-saddles, 2-saddles, and maxima respectively. The edges in the interior are the descending arcs going between 1-saddles and minima. The ascending arcs are on the boundary connecting 2-saddles to maxima. The figure on the left shows an isosurface with three components about to merge at two 1-saddles. The figure on the right shows an ascending disk originating at one of these 1-saddles.
$\includegraphics[height=2.2in]{Figs/anal_d1_a1_iso.ps}$ $\includegraphics[height=2.2in]{Figs/anal_d1_a1_a2.ps}$

Critical points. The nodes of the Morse-Smale complex represent the critical points of the function containing information about the local features in the data. We display the four types of critical points as color-coded spheres. The number of critical points is large in the presence of multiple features or flat regions. We use symbolic perturbation [9, Section 1.4] to remove the flat regions and obtain a Morse function and filter out the unimportant critical points using a cancellation procedure that first ranks them based on their topological persistence [11]. A geometric realization of this cancellation has been accomplished recently for surfaces [4,10]. Performing a cancellation in 3D is a challenging problem that has been addressed recently [13].

**Fig. 7.** Screenshot of the Morse-Smale complex visualization tool. The panel on the right can be used to set various parameters for displaying arcs, disks, and isosurfaces.
$\includegraphics[height=2.5in]{Figs/screenshot.ps}$

We develop a visualization tool using VTK's class libraries that displays the above-mentioned substructures of the Morse-Smale complex. In addition to various substructures, our visualization tool also allows the user to extract isosurfaces. . Figure 7 shows a screen shot of the software. Figure 6 shows the visualization of a synthetic function using isosurfaces and substructures of the Morse-Smale complex. The function value at each point in a unit cube is computed as the minimum distance from a set of three specified interior points. The isosurfaces in the figure indicate the minima corresponding to these three points. The edges in the interior are descending arcs and the ones on the boundary are ascending arcs. The surface is an ascending disk origination at a 1-saddle lying between two of the minima. The ascending disk can be extracted from the Morse-Smale complex as a collection of quadrangles containing the 1-saddle and some maximum as one of its nodes. Given our data structures, this extraction can be done efficiently.

We illustrate how substructures of the Morse-Smale complex, called ascending arcs, trace the shape of isosurfaces using data from two application domains. These ascending arcs are rendered efficiently and hence used as an exploratory tool to identify interesting isosurfaces. Our results for the cryo-EM data are consistent with the structures of the macromolecules determined by the scientists who performed the electron microscope imaging [21]. Similar structures have also been reported by De-Alarcón et al. who use semi-automatic methods to model the shape of biological macromolecules from low resolution density maps [8]. Descending arcs also play a similar role as the ascending arcs in the visualization and analysis process. 2D substructures, namely ascending and descending disks, can be used to locate atoms in molecules from a given electron density distribution. Such a topological approach is proposed by the Atoms in Molecules (AIM) theory [1]. Morse-Smale complexes, constructed using existing combinatorial algorithms, can be used to determine the topologically defined atoms in a stable manner. We might however require the use of adaptive meshing as proposed by Takahashi et al. in order to automatically remove flat regions from consideration during the computation [22]. The adaptivity also makes the computation efficient because it reduces the number of expensive quantum mechanical calculations involved in the evaluation of the function. A promising approach toward the removal of flat regions is to construct a fair Morse function as proposed by Ni et al [18].

We thank Carmen San Martin and Jose Maria Carazo for their expert comments regarding the structure of the DNA helicase. Initial discussions with Herbert Edelsbrunner and John Harer on constructing Morse-Smale complexes were very useful. Research by the first author was supported by NSF under grant CCR-00-86013. This work was performed under the auspices of the U.S. Department of Energy by University of California Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48.

A Topological Method for Analysis of 3D Scalar Functions

Footnotes