Detecting Critical Regions in Multidimensional Data Sets

Revised version of the Research Report No. 2007-52

Madjid Allili^{†}, David Corriveau, Sara Derivière^{∗}, Tomasz Kaczynski^{†}

January 13, 2009

Abstract

We propose a new approach, based on the Conley index theory, for the detection and classification of critical regions in multidimensional data sets. The use of homology groups makes this method consistent and successful in all dimensions and allows to generalize visual classification techniques based solely on the notion of connectedness which may fail in higher dimensions.

1 Introduction

In the geometric study of scalar functions f : X → ℝ on a multidimensional data set X, the goal is to extract features that enhance our understanding of the measured data and allow us to represent the data in a structure that is easy to manipulate for subsequent processing. In the case of a twodimensional data set X, f can be geometrically interpreted as a height field. The features of interest are critical points of f, that is, the points where the gradient ∇f is nul. They are interpreted as peaks, pits, or mountain passes, also called saddles. Once the critical points are extracted, different techniques are used to analyze relationships between them and to trace structures such as ridge and ravine lines, and isolines. In the case of data of higher dimensions, the geometric interpretation of critical points is more complex but those points play equally important role in further investigation, such as the construction of the level sets given by f = c. These level sets are interpreted as isolines in dimension 2, isosurfaces in dimension 3 and so on.

∗supported by the Tomlinson Visiting Scholarship of Bishop’s University

†supported by the NSERC of Canada discovery grant

The mathematical tool commonly used to detect and classify critical points is the Morse theory [14], [6, Chapter 7]. This theory relies on very strong smoothness and non-degeneracy assumptions, described in Section 2, which appear not realistic in discrete models. Among the non-degeneracy assumptions is one, that critical points are isolated and have distinct critical values. This assumption fails when considering, for example, a height field of a terrain with water surfaces as minimum plateaus, sandbars near a seashore may become saddle isolines at a low tide, and ridges of volcano craters as maximum closed isolines.

In [2], we have introduced the notion of critical regions in the case of a height field for a 2-dimensional square lattice. In that case, the detection and classification of a critical region C can be done by studying the number of connected components of its upper and lower wraps, \overline{wrap} (C) and \underline{wrap} (C). These are, roughly speaking, the subsets of their immediate neighborhoods where the values of the function are respectively greater and smaller than f(C). Our approach is to some extent influenced by a classification of critical regions for continuous but not necessarily smooth functions in ℝ^{3} given by Weber et al. [20].

In many applications of the classical Morse theory to imaging science, one spends a lot of effort on forcing data so to obey the assumptions of the theory. This practice is validated by the argument, that a small perturbation of a given function can bring it to the generic case where all critical points are non-degenerated, isolated, and correspond to distinct critical values. Among most systematic studies of that kind one could point the work of Edelsbrunner et al. [9, 10]. This approach may be natural in some problems but there are many other situations where it is not applicable. One, addressed by Arnold [4, 1984], arises in the study of parameterized family of equilibrium states. A small perturbation can remove a degeneracy at one parameter value but it would create a new one at another value.

Another arises in the singularity theory [15], [4, 1993]. Critical regions of polynomial functions of several variables are extensively studied there, under the name of singular curves and surfaces. The Whitney-Cayley umbrella [4, 1993] given by

x^{2} − y^{2}z = 0

is a simple classical example of an algebraic level surface in R^{3} self-intersecting at a singular curve (the z-axis). Any level set g(x) = c, where x ∈ ℝ^{d}, can be turned into a critical region of

f(x) = (g(x) − c)^{2}h(x).

Thus one can easily construct critical regions with a non-trivial topology such as a circle or a wedge of several circles. Such sets arise for example in the study of Milnor fibers [11] of complex functions.

In computational geometry, there are well-developed methods for construction of non-singular level surfaces f = c. However, the mesh generation algorithms break near the singular zone ∇f = 0. A particular attention is required for critical regions of saddle-type, because these are the curves where the level surface self-intersects. In singular zones, one has to investigate the local topology nearly on a pixel level, in order to determine how to locally complete the construction of the studied surface so to preserve its global properties.

Since we now work on multidimensional data, we have to be cautious about carrying over intuitions from planar geometry to higher dimensions: There is a tendency in geometric visualisation research to overlook, for example, the fact that the number of connected components may no longer be sufficient to describe changes in topology of studied sets when parameters change. One has to take into account such topological invariants as homotopy type, homotopy groups, or homology groups. We focus on homology descriptors, because there is a vast library of convenient programs computing homology such as [7], while the homotopy type is not constructive, and the computability of homotopy groups is problematic. In Section 2, we recall the basic terminology concerning the Morse and Conley indices and we give examples motivating the use of homology descriptors.

In Section 3, we propose a method for detecting and classifying critical regions of discrete functions inspired by the Conley index theory. This theory is a generalization of the Morse theory in the following sense: The Morse index discussed in Section 2, classifies the geometric nature of isolated critical points of a smooth non degenerated function f. The gradient field ∇f of that function generates a flow whose equilibrium points are the critical points of f. The Conley index applies not only to gradient flows but to arbitrary continuous flows as well as to discrete dynamical systems, and it classifies not only equilibrium points but also isolated invariant sets. In particular, a critical region can be presented as an isolated invariant set and its Conley index is a natural substitute of the Morse index. We give algorithmic definitions of regular and critical regions inspired by the Conley index theory for functions defined on a unit-size cubical grid in ℝ^{d}. This choice is well justified in practice, because many applications in 3D imaging science produce regularly gridded data by sampling scalar fields at uniform intervals of time and space. This data is given in the form of scalar functions defined on the vertices of a hypercubic decomposition of space. Regularly gridded data does not require an explicit storage of cell adjacency information, which would be necessary, if we wanted to subdivide each cube into tetrahedral grids. Thus it results in lower complexity and lower storage requirements. Note that irregularly gridded data is typically organized into tetrahedral cells using techniques such as the ones based on Voronoi diagrams and Delaunay triangulations [6]. Our method is valid in all dimensions. It requires verifying whether or not relative homology groups of certain sets are trivial. In our case, the verification is implemented at a low cost, since it reduces to so-called elementary collapses, see [12, Sec. 2.4].

In Section 4, we derive conclusions related to the concept of the Morse connection graph (MCG) introduced in [3]. This is a graph whose nodes are critical components and edges display the existence of trajectories connecting them. An new algorithm building MCG for planar images and its implementation are presented in [2]. The method is based on the construction of two multivalued maps F_{+} and F_{−}, defined on a set of pixels whose values are subsets of neighboring pixels. The iterations of F_{+} and F_{−} imitate instantaneous movements of a point along a trajectory of a flow. There are some works somewhat similar in spirit to our approach e.g. Edelsbrunner [8] or Boczko at al. [5] on polygonal flows. The second work is motivated by applications to dynamical systems rather than imaging. Both are based on building adapted triangulations, while we are using a previously mentioned fixed cubical grid. In this paper, we reinforce the MCG Algorithm of [2], by presenting its improved and simplified formulation, extended to arbitrary dimensions.

In Section 5, we discuss some open questions and prospects for future applications to construction of isosurfaces.

2 Motivating examples from the Morse and Conley theories

Let M be a smooth manifold and f : M → ℝ a function of class C^{2}. At this point we do not assume that M is compact, so we may consider M = ℝ^{d}. A point p ∈ M is critical if the gradient ∇f vanishes at p and it is called regular otherwise. The function f is called a Morse function if all of its critical points p are non degenerate. This condition implies in particular that p is an isolated critical point.

Given a Morse function f, the index of any critical point p, denoted by λ(p), is the number of negative eigenvalues of the Hessian H_{f}(p). If λ(p) = 0, p is a local minimum and if λ(p) = d, it is a local maximum. The values 0 < λ(p) < d, characterize saddles which may topologically differ from each other when d > 3.

Consider a function whose critical points are still isolated but possibly degenerate. Then there may occur a critical point p of “inflection type”. which is inessential in the sense that a perturbation of f in a small isolating neighborhood N of p may cause the removal of singularity from N. One of possible topological definitions may be formulated in terms of the sublevel and superlevel sets

N_{n} = {x ∈ N | f(x) < f(p)},
N_{p} = {x ∈ N | f(x) > f(p)},

The intersection of their closures is the level set

N_{z} = {x ∈ N \ {p} | f(x) = f(p)}.

Whether p is a regular point or an inessential critical point, if N is a disc, the sets N_{n} and N_{p} should be its topological half-discs. At a maximum or minimum point, one set is N \ {p} and another one is empty. At a critical saddle such as p = 0 of the function f(x, y) = x^{2}−y^{2}, each of the sets N_{n} and N_{p} consist of two disjoint wedges limited by the crossing lines x^{2} = y^{2}. We may be tempted to use the numbers n_{n}, n_{p} and n_{z} of connected components of N_{n}, N_{p}, and N_{0} respectively, to distinguish various types of topological critical points. This approach seems to work well in dimensions 2 and 3 but it fails to characterize critical points in dimensions higher than 3 even in the classical Morse setting, as the following example illustrates.

Example 2.1 Let f : ℝ^{4} → ℝ be given by

f(x, y, z, t) = x^{2} + y^{2} − z^{2} − t^{2}.  (1)

Then p = 0 is an isolated non-degenerate critical point with the value 0. It is a saddle point and its Morse index is λ(p) = 2. Consider its neighborhood N = \overline{B}(0, r) the closed ball of radius r centered at 0. Since f is radially homogeneous, the choice of r is not important so we may assume r = 1, For the same reason, the radial projection of N \ {0} onto the unit sphere S^{3} given by

x^{2} + y^{2} + z^{2} + t^{2} = 1

is a homotopy deformation. Consider the discussed sets N_{n}, N_{p} and N_{z} deformed to S^{3}:

L_{n} = N_{n} ∩ S^{3} = {(x, y, z, t) ∈ S^{3} | x^{2} + y^{2} < z^{2} + t^{2}},
L_{p} = N_{p} ∩ S^{3} = {(x, y, z, t) ∈ S^{3} | x^{2} + y^{2} > z^{2} + t^{2}},

and

L_{z} = N_{z} ∩ S^{3} = {(x, y, z, t) ∈ S^{3} | x^{2} + y^{2} = z^{2} + t^{2}}.

Since L_{z} can be expressed by x^{2}+y^{2} = 1/2 = z^{2}+t^{2}, it is a torus (product of two circles). Its complements L_{n} and L_{p} in S^{3} are “donuts” that is, products of a circle by a 2-disc. All three sets are connected. The same is true for a regular point, hence the number of connected components does not permit to distinguish a saddle from a regular point in ℝ^{4}.

In the study of critical regions or, more explicitly, connected regions consisting of non-isolated critical points, the three-dimensional space is sufficient to construct examples showing that the connected components of sublevel and superlevel sets do not provide sufficient information.

Example 2.2 Let M be a “donut-shape” solid in ℝ^{3} obtained by revolution of the disc D_{r_{1},r_{2}} given by

(x − r_{2})^{2} + z^{2} ≤ r_{1}^{2}, y = 0

about the z-axis, where 0 < r_{1} < r_{2}. Its boundary ∂M is a torus. We use toric coordinates (r, θ, ϕ) in M, that is, (r, θ) are polar coordinates for D_{r_{1},r_{2}} and ϕ is the angle of revolution. First, we want to define f(r, θ, ϕ) at ϕ = 0 by

f(r, θ, 0) = (x − r_{2})^{2} − z^{2} = r^{2}(1 − 2 sin^{2} θ).

The center of D_{r_{1},r_{2}} is a simple saddle of this function on D. We extend this definition to M by the formula

f(r, θ, ϕ) = r^{2} (1 − 2 sin^{2}(θ − ϕ/2)) .  (2)

The revolution of the center of D produces a circle S of radius r_{2} given by r = 0. Both f and ∇f vanish on S, so this circle is a singular curve of f forming the self-intersection of the isosurface f = 0. Moreover, all points on it are simple saddles of the function f restricted to the plane perpendicular to the circle.

Let now N be a neighborhood of S in M defined in toric coordinates by r ≤ δ, for any chosen 0 < δ < r_{1}. Consider the torus T := ∂N. At ϕ = 0, the boundary of the disc D_{δ,r_{2}} is composed of four arcs limited by isolines at θ = ±π/4, ±3π/4, with f > 0 on the left and right arcs A_{L}, A_{R}, and f < 0, on the upper and lower arcs A_{+}, A_{−}. The arc A_{L} is revolved to A_{R} forming a band on T and, similarly, A_{+} is revolved to A_{−}. Thus each of the sets L_{p} := N_{p} ∩ T and L_{n} := N_{n} ∩ T is a band wrapped twice around the torus T. The set L_{z} := N_{z} ∩ T consists of two circles separating the band L_{p} from L_{n}. By the same argument on radial projection as in the previous example, each of the sets L_{p}, L_{n} and L_{z} is a deformation retract of N_{p}, N_{n}, and N_{z} respectively, so it has the same homotopy type. Thus N_{p} and N_{n} are connected. If we wanted to use this information only, S would be classified as a regular region. In the next section, we discretise this example and show that our homology criterion classifies it as a saddle-type region.

Below we give a brief overview of the Conley index [17, 16] in the context of a flow generated by the gradient field of a C^{2} function f.

A subset S ⊂ X is called an invariant set of φ, if φ(t, x) ∈ S for all x ∈ S and all t ∈ ℝ. It is called isolated if there exists a neighborhood N such that S is the maximal invariant set in N. The main purpose of the Conley index theory of flows is to describe isolated invariant sets.

Let N be the closure of a bounded region in X and ∂N its boundary. A set L ∈ N is called the exit set of N if given x ∈ N and t_{0} > 0 such that φ(t_{0}, x) ∉ N, then there exists 0 ≤ t_{1} ≤ t_{0} such that φ(t_{1}, x) ∈ L. Obviously, we must have L ⊂ ∂N. Here is a homological version of the classical Ważewski principle [19] which inspired the definition of the Conley index.

Theorem 2.3 [12, Prop. 10.40] Suppose that L is a closed subset of ∂N. If H_{∗}(N, L) is non trivial then N contains a non-empty invariant set in its interior.

The original Ważewski principle formulated in terms of deformation retraction gave rise to the Conley theory. Given a maximal invariant set S isolated by N, the homotopical Conley index C(S) is the relative homotopy type [N, L] of the pair (N, L). Since the homotopy type is not constructive, in modern application-oriented formulations of the Conley index theory [16], C(S) tends to be replaced by computable descriptors such as the homological Conley index given by CH_{∗}(S) = H_{∗}(N, L). Additional assumptions on N and L are imposed in order to make the index definition dependent only on S, and not on the choice of N, and in order to achieve desired topological properties, such as stability of the index with respect to small perturbations of the flow. We avoid introducing those assumptions here because, in our work, we will not use the Conley index literarily but only as an inspiration for the digital setting. The following well known example shows how the homological Conley index generalizes the Morse index.

Example 2.4 Let p be a non-degenerate critical point of the Morse index λ(p) = m of a function f defined on a d–dimensional manifold. Consider the flow of the differential equation ẋ = −∇f(x). Let N be an admissible isolating neighborhood for the Conley index of {p} and L its exit set. Then

CH_{k}(\{p\}) = H_{k}(N, L) ≅ H_{k}(S^{m}, \{s_{0}\}) ≅ \{ \begin{matrix}ℤ & if k = m, \\ 0 & otherwise, \end{matrix}

where S^{m} is the m-dimensional sphere and s_{0} is a chosen base point in S^{m}.

3 Critical regions of discrete functions

Our goal in this section is to introduce a discrete analogy of critical points and critical regions of functions f : ℝ^{d} → ℝ.

The geometric structure we are going to work with is the cubical grid K defined in [12, Sec. 2.1]. Recall that K is the collection of elementary cubes of the form

Q = I_{1} × I_{2} × ⋯ I_{d}

where I_{j} = [k, k + 1] or I_{j} = {k}, k ∈ ℤ (the set of integers), including the empty set. We denote by K^{d} the subset of K consisting of d–dimensional cubes, also called full cubes. These are the cubes which have no degenerate intervals I_{j} = {k} in their expression, equivalently, which are not proper faces of other cubes in K. Note that K is a particular case of a regular cellular complex. A set X ⊂ ℝ^{d} is called a cubical set, if it is a finite union of elements of K, and it is called full if it is a finite union of elements of K^{d}. We should distinguish between a finite set X ⊂ K and the support or polytope of X which is the cubical set X ⊂ ℝ^{d} given by

X = |X | = ⋃ X .

Given a cubical set X, K(X) is the restriction of the grid K to the cubes contained in X.

Would chosen units appear too coarse, one may apply the rescaling isomorphism Λ^{k} : ℝ^{d} → ℝ^{d},

Λ^{k}(x_{1}, x_{2}, …, x_{d}) = (kx_{1}, kx_{2}, …, kx_{d}),

where k ∈ ℤ is called a scaling factor. The corresponding refined grid is the image of K under the inverse rescaling Λ^{1/k} = (Λ^{k})^{−1}:

Λ^{1/k}K = {Λ^{1/k}Q | Q ∈ K}.

Here is the natural extension of the concept of neighborhood to the combinatorial setting.

Definition 3.1 Let A be a bounded set in ℝ^{d}. The combinatorial unit–scale wrap of A is a subset of K defined by

wrap(A) = {P ∈ K^{d} | P ∩ A ≠ ∅}.

It is a finite set, so its support denoted by wrap(A) = |wrap(A)| is a cubical set. The factor k scaled wrap of A is a subset of Λ^{1/k}K defined by

wrap^{1/k}(A) = {P ∈ Λ^{1/k}K^{d} | P ∩ A ≠ ∅}.

Its support with respect to the refined grid is denoted by wrap^{1/k}(A).

The study of scaled wraps is motivated by the following theorem.

Theorem 3.2 (Allili [1]) Given a cubical set A and a scaling factor k ≥ 3, the inclusion A ↪ wrap^{1/k}(A) induces an isomorphism in homology.

Definition 3.3 The factor k scaled exterior boundary of a cubical set A is the topological boundary ∂ wrap^{1/k}(A) of A.

In practice, the values of functions f : ℝ^{d} → ℝ are only known at finitely many grid points. Our approach, inspired by the combinatorial multivalued maps used in [13, 12, 18] is to define combinatorial maps on full elementary cubes rather that on their vertices. Thus we consider discrete functions^{1}

f : X^{d} → ℝ, where X^{d} ⊂ K^{d}.

This approach is very natural if we view the smallest objects in an image, n-pixels, as d-cells of a cubical grid K. The passage from our definition to a vertex definition is by simple translation of coordinates x_{i} ↦ x_{i} + 1/2, i = 1, 2, … d. Conversely, if a function f is given on a cubical set X = |X^{d}| ⊂ ℝ^{d}, we may define the discretization of f on X^{d} by taking its values at the center of each elementary full cube.

1As a matter of fact writing f : X^{d} → ℝ is a slight abuse of notation. In fact, we consider f : X^{d} → 10^{−k}ℤ for a fixed k ∈ ℤ, ℤ integers.

Definition 3.4 Let X^{d} ⊂ K^{d}. Consider a function f : X^{d} → ℝ and an elementary cube Q ∈ X^{d} such that wrap(Q) ⊂ X^{d}. The upper wrap of Q is defined by

\overline{wrap} (Q) := \{P ∈ wrap(Q) | f(P) > f(Q)\} .

Analogously, the lower wrap of Q is given by

\underline{wrap} (Q) := \{P ∈ wrap(Q) | f(P) < f(Q)\} ,

and the level wrap of Q is given by

wrap_{z}(Q) := {P ∈ wrap(Q) | f(P) = f(Q)} .

The extension of the upper and lower wraps to cubical sets A = |A| is more delicate. We assume that wrap(A) ⊂ X. Since the values of f may vary, we want the inequalities in the earlier formulas to be satisfied locally:

\overline{wrap} (A) := \{P ∈ wrap^{∗}(A) | f(P) > f(Q) for all Q ∈ A ∩ wrap^{∗}(P)\} ,

\underline{wrap} (A) := \{P ∈ wrap^{∗}(A) | f(P) < f(Q) for all Q ∈ A ∩ wrap^{∗}(P)\} .

wrap_{z}(A) := wrap^{∗}(A) \ (\underline{wrap} (A) ∪ \overline{wrap} (A))

where wrap^{∗}(A) = wrap(A) \ A. As for wraps, the notation \overline{wrap} and \underline{wrap} is used for the supports of upper and lower wraps \overline{wrap} and \underline{wrap} .

Analogous terminology and notation is carried over to the scaled wraps introduced in definition 3.1.

The upper and lower wraps are analogies of the exit set and entrance set from the Conley index theory. In [2], the algorithms detecting, first, critical pixels in ℝ^{2} and, secondly, their critical components C ⊂ K^{2}, where based on the number of edge-connected components of combinatorial upper and lower wraps \overline{wrap} (C) and \underline{wrap} (C). We did not need to consider scaled wraps introduced above, so we need to explain why do we use them now. As we showed in the previous section, the connectivity concept is not sufficient in dimensions higher than 2. In order to compute homology groups, we need to pass to geometric carriers. Therefore, following the approach chosen in [2], the set

N := |wrap(C)| = \overline{wrap} (C)

would be a natural candidate for the isolating neighborhood of C ⊂ K^{2}, and \overline{wrap} (C) and \underline{wrap} (C) respectively, for its entrance and exit sets. However, there is a problem related to the fact that while the combinatorial sets Figure 1: The region C of pixels with the value 0 is classified accordingly to the algorithm in [2] as regular, because both \overline{wrap} (C) and \underline{wrap} (C) are non-empty and edge-connected. However \overline{wrap} (C) contains ∂wrap(C), as in the case of a minimum.

\overline{wrap} (C) and \underline{wrap} (C) are disjoint, their carriers may intersect, resulting in misleading topological information. This is illustrated in Figure 3. One possible remedy to the described problem is to consider open wraps instead of the closed ones. However, computing homology of open polyhedra is a substantially more complex task than that of compact polyhedra. We shall achieve the same effect, while staying in the class of compact cubical sets, by considering a scaled wrap bounded by the contour around C marked in Figure 3.

Definition 3.5 Let X^{d} ⊂ K^{d} and k ≥ 3. Consider a function f : X → ℝ and an elementary cube Q ∈ X^{d} such that wrap(Q) ⊂ X^{d}. Define

N := wrap^{1/k}(Q), L_{n} := \underline{wrap} (Q) ∩ ∂N L_{p} := \overline{wrap} (Q) ∩ ∂N,

and

L_{z} := wrap_{z}(Q) ∩ ∂N.

The elementary cube Q is called ordinary if H_{∗}(N, L_{p}) = 0 and H_{∗}(N, L_{n}) = 0. Otherwise, it is called singular.

Here is the related algorithm.

Algorithm 3.6 Detecting singular cubes For each elementary full cube Q build N, ∂N, L_{p}, L_{n}, L_{0} H := H_{∗}(N, L_{p}) = 0 and H_{∗}(N, L_{n}) = 0 if H = TRUE then Q is ordinary

else Q is singular endif

As it was pointed out in [2] in the planar case, the adjacent singular cells may, in a sense, cancel each other. In other words, the type of criticality of a singular cube cannot be decided without looking at neighboring cubes. This is taken care of by the following definitions.

Definition 3.7 Let X^{d} ⊂ K^{d}, k ≥ 3 and consider a function f : X^{d} → ℝ. A set C of singular elementary cubes in K^{d} is called an isolated singular component if

(a) wrap(C) ⊂ X;

(b) C = |C| is connected;

(c) C is isolated in the sense that any P ∈ wrap^{∗}(C) is ordinary.

The isolating neighborhood of C = |C| is defined by

N = wrap^{1/k}(C).

The lower, upper, and level sets of C are respectively defined by

L_{n} = \underline{wrap} (C) ∩ ∂N, Ł_{p} = \overline{wrap} (C) ∩ ∂N, L_{z} = wrap_{z}(C) ∩ ∂N

Definition 3.8 An isolated singular component C is called regular if H_{∗}(N, L_{p}) = 0 and H_{∗}(N, L_{n}) = 0. Otherwise it is called a critical component. A singular cube Q is called critical, if it belongs to a critical component, otherwise, it is called regular. An ordinary cube Q is, by definition, regular.

In a sense, singular-regular cubes are an analogy of an inessential or removable singularities in mathematical analysis. Given any singular cell Q ⊂ K^{d}, the component C(Q) = C is constructed by the following algorithm:

Algorithm 3.9 Sorting components in ℝ^{d} For each singular cube Q, C := {Q} while P ∈ wrap^{∗}(C) ∩ X is singular

C := C ∪ {P }

build wrap^{∗}(C) endwhile build N, ∂N, L_{p}, L_{n} do

H := H_{∗}(N, L_{p}) = 0 and H_{∗}(N, L_{n}) = 0 if H = TRUE then C is a regular component else if L_{n} = ∂N then C is a maximum component else if L_{p} = ∂N then C is a minimum component else C is a saddle component

endif break

A practical implementation of the Algorithms 3.9 and 3.6 is a work in progress. In order to illustrate the ideas behind these algorithms, we shall construct discrete analogies of the Examples 2.1 and 2.2, in the setting of cubical grids.

Example 3.10 We come back to Example 2.1 with the unit euclidean sphere in ℝ^{4} replaced by the unit sphere with respect to the maximum norm ||x|| = max_{i=1...4} |x_{i}|. For the convenience of presentation, we consider a cubical grid shifted so that the origin of coordinates is the center of a full elementary cube σ, and rescaled so that N = wrap(σ) = [−1, 1]^{4}. By discretizing f so that the values of f in the centers of the 3–dimensional faces in ∂N determineed L_{p}, L_{n} and L_{z}, we obtained:

L_{p} = ({−1, 1} × [−1, 1] × [−1, 1]^{2}) ∪ ([−1, 1] × {−1, 1} × [−1, 1]^{2})
= ∂[−1, 1]^{2} × [−1, 1]^{2},
L_{n} = ([−1, 1]^{2} × {−1, 1} × [−1, 1]) ∪ ([−1, 1]^{2} × [−1, 1] × {−1, 1})
= [−1, 1]^{2} × ∂[−1, 1]^{2},
L_{z} = ∅.

We used the software CHomP [7] and obtained the result:

H_{∗}(N, L_{p}) ≃ H_{∗}(N, L_{n}) ≃ \{ \begin{matrix}ℤ & if k = 2 \\ 0 & otherwise \end{matrix}

the Algorithms 3.6 and 3.9, the cell σ is classified as a saddle.

Example 3.11 We present a cubical analogy of Example 2.2 and use the homology program of [7] to classify the criticality of the studied region. We choose the radius r_{2} = 2. The circle S discussed in Example 2.2 is inscribed in the cubical donut-shape region

C = ([−3, 3]^{2} \ [−1, 1]^{2}) × [−1, 1].

The values of the function f given by (2) at the centers of the cubes of C = K^{3}(C) are either equal or close to 0. By means of a suitable tresholding, we obtained its discretization f̄ with the value 0 on all cubes of C. Next, the space has been rescaled by the factor of 10 in order to exclude the z-axis and to obtain a discrete approximation f̄ that reflects the behavior of f in the neighborhood N = wrap^{1/10}(C) of C. To avoid working with fractional grids, we redefined C as

C := ([−30, 30]^{2} \ [−10, 10]^{2}) × [−10, 10]

and N as N := wrap(C). The direct implementation of the Algorithm 3.9 has confirmed that C is an isolated singular component and it generated the sets L_{p} and L_{n}. Using the software CHomP [7], we obtained the result:

H_{∗}(N, L_{p}) ≃ H_{∗}(N, L_{n}) ≃ \{ \begin{matrix}ℤ & if k = 2 \\ ℤ & if k = 1 \\ 0 & if k = 0 \end{matrix}

which shows that C is a saddle component. The computed homology is, in fact, the homology of the wedge of the circle and the 2D sphere relatively the wedge point. The snapshots of C and of the solid-cube sets \overline{wrap} (C) and \underline{wrap} (C) giving rise to L_{p} and L_{n} are shown in Figure 3.11.

Figure 2: The sets C, \overline{wrap} (C) and \underline{wrap} (C) in Example 3.11.

We finish this section with several remarks.

(1) In the definition of the homological Conley index for flows, the condition H_{∗}(N, L_{p}) ≠ 0 would be sufficient for criticality, one would not need to take N_{n} into account. This is due to the fact, that the flow of −∇f has the same trajectories as that of ∇f oriented in the opposite direction. In the digital case we study here, there is no such symmetry, so we need to take both L_{p} and L_{n} into consideration. Such differences between the smooth and digital case are discussed in grater detail in [2].

(2) It is possible to modify the definition of the sets L_{p} and L_{n} so to build them of full cubes rather than (d − 1)–dimensional faces in ∂N. Namely, one can define N as the two-layer wrap of C:

N := wrap^{1/k}(wrap^{1/k}(C)),

while L_{p} and L_{n} can be defined as the intersections of \overline{wrap} (C) and, respectively, \underline{wrap} (C) with the second layer of this scaled wrap. If C consists of a single cube, the scaling by k = 3 is sufficient but in general, one should use scaling by k = 5 for topology preservation. This approach may be convenient for the uniformity of data structure since it does not require generating faces of pixels within the program that builds the pairs (N, Lp) and (N, L_{n}). Note that generating faces takes place implicitly in the homology computation.

(3) The condition H^{∗}(N, L_{p}) = 0 (and the same for L_{n}) in our definitions and algorithms can be replaced by a stronger condition that the homotopy type [N, L_{p}] of the pair of spaces (N, L_{p}) is non-trivial. However, as already commented in the introduction, such a condition could not be verified in practice. A yet stronger but combinatorial condition would be “N collapses to L_{p}”, which means that N can be deformed to N_{p} by a chain of elementary collapses described in [12, Sec. 2.4]. Elementary collapses preserve not only homology but also homotopy type and the algorithm Collapse actually is the fastest method of verifying that H_{∗}(X, L_{p}) = 0.

(4) Note that there are many topologically distinct types of saddle components. In the case of a single elementary cube, one can characterize it as a minimum, maximum, or a k–fold saddle, k ∈ ℕ, by means of H_{∗}(N, L_{p}) and H_{∗}(N, L_{p}). For larger components C, the topology of C may be nontrivial and one has to take into account H ∗(C). An attempt of classifying all cases that may occur would be equivalent to solving most of the challenges of the singularity theory.

4 Dynamics of a multidimensional scalar field and associated Morse Connections Graph

In [2], we gave a detailed description of the dynamics of a given height field and we associated a structure called the Morse connections Graph whose nodes are critical components of the height field and the edges show different connecting trajectories between the nodes. This graph encodes a summary of topological features of the data and makes it easy to manipulate in many imaging applications such as shape representation and retrieval. We present here a very simple reformulation of these concepts in the context of a multidimensional scalar field.

Let now K be the cubical grid in ℝ^{d} and X^{d} ⊂ K^{d}. Consider a function f : X^{d} → ℝ. We want to introduce a multivalued map on X^{d} such that the dynamics of its iterates in the sense of [13] would exhibit the dynamical properties analogous to those of the gradient flow of a Morse function of f. First, we define the distance dist (Q, P) between two elementary cubes as the euclidean distance between their centers. Note that, if Q, P ∈ X^{d} intersect at a common face of dimension k, then

dist (Q, P) = \sqrt{d − k}.

The directional derivative of f at Q in the direction of an adjacent cell P is defined by

∂_{P} f(Q) = \frac{f(P) − f(Q)}{dist (Q, P)}.

By convention, ∂_{Q}f(Q) = 0. Given a critical cell Q, we put

m(Q) = min_{P ∈\underline{wrap} (Q)} ∂_{P} f(Q) and M(Q) = max_{P ∈\overline{wrap} (Q)} ∂_{P} f(Q).

Definition 4.1 Given f : X^{d} → ℝ, the ascending system is the semidynamical system of iterates of the map F : X^{d} → Pℝ, which is defined on cells Q ∈ X^{d} as follows. If Q is critical, we put F_{+}(Q) = Q. If Q is regular, we put

F_{+}(Q) = \{P ∈ \overline{wrap} (Q) | ∂_{P} f(Q) = M(Q)\} .

The descending system is generated by the map F_{−} : X^{d} → Pℝ defined as before with \overline{wrap} (Q) replaced by \underline{wrap} (Q) and M(Q) replaced by m(Q).

Note that the inverses of F_{+} and F_{−} may have empty values (and these two maps are not mutual inverses), hence the prefix ”semi” in the name for the system.

Definition 4.2 Let F : X → PX be anyone of the two maps introduced in Definition 4.1. Consider a critical component C ⊂ X and let C = |C|. The stable and unstable manifolds of C relatively to F are defined by

\begin{aligned}W^{u}(C, F) & = & F(wrap(C), ℤ_{+}) & = & ⋃_{Q∈wrap(C),n≥1} F^{n}(Q); \\ W^{s}(C, F) & = & F(wrap(C), ℤ_{−}) & = & ⋃_{Q∈wrap(C),n≥1} F^{−n}(Q). \end{aligned}

If F is clear from the context, we will just write W^{s}(C), W^{u}(C).

The following theorems proved in [2] remain valid in arbitrary dimension:

Theorem 4.3 Let X^{d} be finite and let F be anyone of the two maps introduced in Definition 4.1. If C and D are critical components such that W^{u}(C) ∩ W^{s}(D) ≠ ∅, then there exists a trajectory connecting C to D in the sense that it is issued at a cell in wrap(C) and ends at a cell in wrap(D).

Theorem 4.4 Let X^{d} be finite and let F be anyone of the two maps introduced in Definition 4.1. Let {C_{i}}_{i=1,2,...k} be the set of all critical components in X^{d}. Then

⋃_{i=1}^{k} W^{s}(C_{i}) = X^{d} = ⋃_{i=1}^{k} W^{u}(C_{i}).

Let C and D be two critical components of f : X → ℝ. We say there is an upward connection from C to D, denoted C ↗ D, if

W^{u} (C, F_{+}) ∩ W^{s} (D, F_{+}) ≠ ∅.

There is a downward connection from C to D, denoted C ↘ D, if

W^{u} (C, F_{−}) ∩ W^{s} (D, F_{−}) ≠ ∅.

We can now define the Morse Connection Graph as follows:

Definition 4.5 The Morse Connections Graph MCG_{f} = (V_{f}, E_{f}) is a graph whose vertices (or nodes) V_{f} and edges E_{f} are defined as follows:

V_{f} = {critical components of f};
E_{f} = {(C_{i}, C_{j}) ∈ V_{f} × V_{f} | C_{i} ↗ C_{j} or C_{i} ↘ C_{j}}

Equivalently, (C_{i}, C_{j}) is an edge of the graph if

W^{u}(C_{i}, F_{+}) ∩ W^{s}(C_{j}, F_{+}) ≠ ∅ or W^{u}(C_{i}, F_{−}) ∩ W^{s}(C_{j}, F_{−}) ≠ ∅.

The algorithm constructing Morse Connections Graph based on the described multivalued maps, and experimentation on two-dimensional image data are presented in [2]. The Morse Connections Graph algorithm itself is dimension independent but an experimentation in higher dimensions has not been done yet at the time of writing this paper.

5 Directions of the future work

A natural question to rise is, weather or not our construction is stable with respect to perturbation of f or with respect to rescaling (that is, with respect to a change of resolution). Since our model assumes a discrete function as an input, and stability concerns continuous data, these questions cannot be answered within the framework of this paper, because that would involve different methods of discretizing, that is interpolating a continuous fonction. Thus interpolation methods should be integrated into the study.

Another interesting direction is to study extensions of the Euler formula in the context of our work. This formula is used (sometimes enforced) as a criterion of correctness of critical points extracted from discrete data. Our two-dimensional experiments test positively when the isolation and nondegeneracy assumptions hold and when there are no critical pixels on the boundary of the region. But the Euler formula cannot hold in presence of critical components. In fact, a critical component itself may have an Euler characteristics different than that of a point. Thus the goal is to find a suitable generalization of the Euler formula that would hold for the degenerate case.

Yet another direction is an application of our results to the problem of isosurfaces which has been the main motivation for our work. Isosurface extraction deals with the problem of generating a surface defined by the preimage of a scalar value c, called an isovalue, under a scalar function f of several variables. It is a powerful tool for visualizing and investigating volumetric scalar fields. Indeed, by varying the isovalue c and studying each corresponding isosurface, one can completely explore a given volume data, at any appropriate level of detail. Moreover, in domains such as medical image processing, many imaging modalities give images (such as CT and MRI scan data images) that are piecewise constant, because the intensity is related to tissue type. It follows that isosurfaces correspond approximately to tissue regions.

Studying the representation of all possible isosurfaces is not realistic in large scale data sets. Therefore, one needs to limit the study to isovalues where “interesting” behavior occurs. However, it is very difficult to determine which isovalues are of interest, since failing to consider certain isovalues may result in missing important features of the volume data. The analysis of the topological properties of a given discrete scalar field with appropriate methods such as the ones inspired by Morse and Conley index theories as is carried out in this work makes it possible to determine critical regions that correspond to the topological changes for which relevant isosurface behavior occurs. By determining the locations of critical components, one can track all the fundamental changes in the scalar data and reconstruct an appropriate isosurface for any given isovalue using meshing techniques and marching cubes algorithms (see [20] and the references therein).

Although a panoply of methods for reconstruction of non-singular isosurfaces are available, the techniques utilized break near the singular zones where one has to investigate the local topology nearly on a pixel level, in order to determine how to locally complete the construction of the studied surface so to preserve its global properties.

References

[1] M. Allili, Une approche algorithmique pour le calcul de l’homologie de fonctions continues, Ph. D. Thesis, Dep. de math. et d’info., Univ. de Sherbrooke (1998).

[2] M. Allili, D. Corriveau, S. Derivière, T. Kaczynski, and A. Trahan, Discrete dynamical system framework for construction of connections between critical regions in lattice height data, Journal of Math. Imaging and Vision 28 2 (2007), 99–111.

[3] D. Corriveau, M. Allili, and D. Ziou, Morse connections graph for shape representation, In Lecture Notes in Computer Science 3708, Springer-Verlag 2005, 219–226.

[4] V. I. Arnold, Catastrophe Theory, Springer-Verlag 1984; —, The Theory of Singularities and Its Applications, Acad. Naz. dei Lincei, Pisa 1993.

[5] E. Boczko, W.D. Kalies, and K. Mischaikow, Polygonal approximation of flows, Topology and Its Appl. 154 (2007), 2501–2520.

[6] J. D. Boissonnat and M. Teillaud (editors), Effective Computational Geometry for Curves and Surfaces, Math. and Visualization Series, Springer-Verlag 2007.

[7] CHomP, Computational Homology Project, http://chomp.rutgers.edu.

[8] H. Edelsbrunner, Surface reconstruction by wrapping finite point sets in space, Discr. and Comp. Geom., in Ricky Pollack and Eli Goodman Festschrift (eds. B. Aronov, S. Basu, J. Pach, M. Sharir), Springer-Verlag 2003, 379–404.

[9] H. Edelsbrunner, J. Harer, and A.J. Zomorodian, Hierarchical Morse-Smale complexes for piecewise linear 2-manifolds, In AMC Symposium on Computational Geometry, Discrete & Compu. Geom. 30 1 (2003), 87–107.

[10] H. Edelsbrunner, J. Harer, V. Natarajan, and V. Pascucci, Morse– Smale complexes for piecewise linear 3-manifolds, In Proc. IEEE Conf. Visualization (2004), 275–280.

[11] K. Houston, Topology of differentiable mappings, in Handbook of Global Analysis, 1213, Elsevier Amsterdam 2008, 493–532.

[12] T. Kaczynski, K. Mischaikow, and M. Mrozek, Computational Homology, Appl. Math. Sci. Series 157, Springer-Verlag 2004.

[13] T. Kaczynski and M. Mrozek, Conley index for discrete multivalued dynamical systems, Topology and Its Appl. 65 (1995), 83–96.

[14] Y. Matsumoto, An Introduction to Morse Theory, Mathematical Monographs 208, Amer. Math. Soc. 2002.

[15] J. Milnor, Singular Points of Complex Hypersurfaces, Annals of Mathematics Studies 61, Princeton Univ. Press 1968.

[16] K. Mischaikow and M. Mrozek Conley Index Theory, Handbook of Dynamical Systems II: Towards Applications, (ed. B. Fiedler), North-Holland 2002.

[17] J. Robbin and D. Salamon, Dynamical systems, shape theory, and the Conley index, Ergodic Theory and Dynamical Systems, 8* (1988), 375–393.

[18] A. Szymczak, Index pairs: From dynamics to Combinatorics and Back, Ph.D. thesis, Georgia Inst. Tech., Atlanta 1999.

[19] T. Ważewski, Une méthode topologique de l’examen du phénomène asymptotique relativement aux équations différentielles ordinaires, Rend. Accad. Nazionale dei Lincei, Cl. Sci. fisiche, mat. e naturali, Ser. VIII, III (1947), 210–215.

[20] G. H. Weber, G. Scheuermann, and B. Hamann, Detecting critical regions in scalar fields, In IEEE TCVG Symposium on Visualization (2003), 1–11.

Department of Computer Science Bishop’s University Sherbrooke (Québec), Canada J1M 1Z7 mallili@ubishops.ca

Département de mathématiques Université de Sherbrooke, Sherbrooke (Québec), Canada J1K 2R1 {d.corriveau, s.deriviere, t.kaczynski}@usherbrooke.ca