Computing Geomorphologically Salient Networks via Discrete Morse Theory

In this paper we give a mathematically sound formalization of volume persistence and its use for network computation. In Section 2 we provide background on discrete Morse theory. Section 3 describes our method to construct a network based on volume persistence. By design these networks are stable under perturbations of small volume, however, like the channels that stem from striations, they connect one side of the terrain to the other. In Section 4 we develop the additional theory necessary to detect salient fingers. The major challenge here is to characterize the area and corresponding terrain volume that defines a finger. It is clear which area lies between two channels that first bifurcate and then merge again. Correspondingly, it is straightforward to compute the volume of terrain between two channels. But what area determines a finger? When are two fingers similar? Which fingers should be retained in a network given a volume threshold for persistence? To address these challenges, we introduce the notion of “spurs” on the terrain and define salient fingers as those that bound spurs of sufficient volume. Finally, in Section 5 we showcase the results of our method as implemented in our tool TopoTide (https://github.com/tue-alga/topotide).

Traditional Morse theory enables us to study a manifold by defining a Morse function, a differentiable function on the manifold, and studying its features such as critical points and stable and unstable manifolds. As the data describing our terrain is not continuous, but consists of discrete points, we assume the terrain is a two-dimensional cell complex $T$: a surface defined by a set of zero-, one-, and two-cells so that the intersection of any two cells in $T$ is in $T$ itself. We then turn to discrete Morse theory [11, 20].

For cells $\tau$ and $\sigma$, we say that $\tau$ is a face of $\sigma$ (or $\sigma$ is a coface of $\tau$) if $\tau \subseteq \sigma$, and we say that $\tau$ is a facet of $\sigma$ (or $\sigma$ is a cofacet of $\tau$) if $\tau$ is a face of $\sigma$ and $dim(\tau )=dim(\sigma )-1$. A discrete Morse function (DMF) is a map $f:T\to ℝ$ such that every cell has a lower (higher) value than its cofacets (facets) with at most one exception. It may be helpful to visualize $f$ as pairing facets to cofacets where these exceptions occur, in an analog to gradient flow from the smooth setting. These pairings are called gradient pairs and we may illustrate them using arrows pointing from facets to cofacets, as in Figure 4. A cell $\sigma$ is critical if $\sigma$ is not paired with any of its facets or cofacets. A critical zero-cell is called a minimum, a critical one-cell is called a saddle, and a critical two-cell is called a maximum.

Given a DMF $f$, an ascending path is a sequence that alternates between zero- and one-cells, “swimming upstream” along the pairings induced by $f$. Similarly, a descending path is a sequence that alternates between two- and one-cells, following the pairings induced by $f$. The union of all ascending paths forms a spanning forest, and so does the union of all descending paths; these spanning forests are dual, with the exception of saddles. For each facet of a saddle, there exists a unique ascending path ending at that zero-cell and originating at a minimum; we define the ascending Morse–Smale complex (or MS-complex for brevity) as a graph $G=(V,E)$ with nodes $V$ corresponding to the minima of $f$, and arcs $E$ corresponding to the pairs of ascending paths connected via saddles between them. See Figure 4. Similarly, for each cofacet of a saddle, there is a unique descending path ending at that two-cell and originating at a maximum. Note that, while a node $v\in V$ corresponds to a zero-cell, an arc $e\in E$ represents the union of potentially many zero- and one-cells. Furthermore, the highest zero-cell along $e$ is the highest facet of its saddle. Each face of $G$ contains a single maximum. For a subset $E^{\prime }\subseteq E$, we write $G-E^{\prime }$ to denote the arc-induced subgraph on $E\setminus E^{\prime }$, i.e., the graph containing only the arcs $E\setminus E^{\prime }$ and their endpoints.

The input to our methods is a digital elevation model (DEM), a set of points and associated elevations. These points and elevations correspond to the elevations of geographic locations in a real-world or simulated terrain. DEMs come in many varieties, for instance triangular irregular networks (TINs) or other point configurations. In this paper, we suppose our DEM data consists of a regular rectilinear grid of measurement points. However, we emphasize that discrete Morse theory along with many of the methods we describe, generalizes to a wide variety of DEM types and structures.

We assume that our input DEM corresponds to a region homeomorphic to a disk, i.e., there are no arches or more complex structures. If any two zero-cells have identical elevation, we perform a slight perturbation so that general position is satisfied [8]. We add three additional zero-cells to an input DEM. First, we add two zero-cells at elevation $-\mathrm{\infty }$ to represent a global source and sink, and connect to the appropriate regions of the DEM (conceptually, the source and sink along this boundary). We then add a zero-cell at elevation $+\mathrm{\infty }$ to represent the land that is not part of the terrain of interest (e.g., the bank of a river), which we assign elevation $+\mathrm{\infty }$, and connect to the entire (new) boundary, including the minimum points. Thus, $T$ is topologically a sphere (see Figure 5). Note that $T$ contains non-square two-cells adjacent to the infinite points, so $T$ does not have an entirely grid-like structure; however, the region of interest, namely, points that correspond to data measurements, is endowed with a grid structure, and it suffices to conceptualize $T$ as cubulated. We call a DEM adapted in this way a terrain, which is the starting point for our work.

A discrete Morse function (DMF) $f:T\to ℝ$ assigns a real value to every cell of $T$: not only the zero-cells, but also the one- and two-cells. On the other hand, our terrain data consists of elevation values for the zero-cells of $T$, but not the one- and two-cells. Thus, we have to extend these measurements on the zero-cells to every cell of $T$ in such a way that we obtain a discrete Morse function, whose critical points and gradient pairs reflect the input measurements. Several choices have been proposed in the literature, see for example [21, 12, 13, 9, 22, 19, 10]. In this paper we use the method proposed in [22]; we sketch the relevant gradient pairings in the following construction.

A main contribution of [22] shows the gradient pairings described above are valid, and do not contradict given elevation data.

Since a set of valid gradient pairings is equivalent to a DMF inducing those pairings, we generally proceed by discussing gradient pairings rather than an explicit function.

Our filtrations depend on a volume-based parameter. Thus, we first need to define the volume of a region in $T$. Conceptually, we think of the elevation of each zero-cell as volume being divided among its four adjacent two-faces, so that each zero-cell contained in a region contributes zero, one, two, three, or four quarters of its volume (height value) to the total, depending on how many of its adjacent two-faces are in the region.

Note that Definition 3 could easily be adapted for more general settings, e.g., if faces do not necessarily have exactly four nodes, or if the size or shape of the faces is salient. We want to know the volume of regions of the terrain that are above a particular elevation, but as elevation is only defined on zero-cells, we introduce the following concept.

We are generally interested in measuring volumes that correspond to continuous regions of the terrain. This is formalized in the following definition.

Next, we introduce an idea that allows us to consider lexicographically lowest paths.

Roughly speaking, Construction 1 defines ascending paths of an MS-complex by introducing zero-/one-cell pairings that are lexicographically low. This has the following consequence.

Consider a splitting arc $e\in E(G_{{\delta }^{\prime }})$. Then the volume persistence of $e$ is greater than or equal to ${\delta }^{\prime }\ge \delta$, meaning $e\in E(G_{\delta })$. Since our choice of $e$ was arbitrary, every splitting arc of $G_{\delta }^{\prime }$ is a splitting arc of $G_{\delta }$, and $V(G_{{\delta }^{\prime }})\subseteq V(G_{\delta })$. Then the lowest paths between two nodes of $G_{{\delta }^{\prime }}$ are also in $G_{{\delta }^{\prime }}$. Thus, $G_{{\delta }^{\prime }}\subseteq G_{\delta }$. $◀$

We also take the following general position assumptions. Let $f$ be a DMF of a terrain $T$. We assume $f$ is injective on the zero-cells of $T$ with a finite height, i.e., they each have a unique height. Furthermore, given $G$, the MS-complex associated to $T$, each arc of $G$ has a unique volume persistence in the sense of Definition 9.

In Section 3, we associate volumes to splitting arcs of $G$, based on the minimum volume of either of the two hills that the arc separates. Here, we take a finer-grained approach and, instead, associate volume to individual one-cells in a red tree. While a one-cell generally does not separate hills like a splitting arc of $G$ does, we take inspiration from this sort of behavior and instead, consider the volumes of the two subtrees formed by the one-cell’s removal, above a given height. Note that it makes sense to discuss “volume of a subtree” in this setting, as the subtrees are collections of two- and one-cells, and thus, the quarter pillar volume (Definition 3) may be nonzero. More formally, we have the following.

See Figure 8. If the volume of a one-cell at a particular height is low, it does not represent a location in the terrain with significant “separation.” Since we eventually use the structure of one-cells that do have significant separation in order to identify the locations of our finger channels, we filter the red tree in the following way.

Note that we now have both a volume and height parameter. If $\delta$ varies within a range that does not change the volume filtered subnetwork, $G_{\delta }$, then decreasing $\delta$ or $h$ means more one-cells become part of $R_{\delta ,h}$, i.e., we have a bifiltration. However, if $\delta$ varies within a range such that $G_{\delta }$ changes, we are no longer guaranteed this property. In particular, if $\delta$ increases and an arc of $G_{\delta }$ is removed, two previously disconnected red trees merge, potentially causing the involved one-cells to have a jump in their volumes. In other words, one-cells are monotonically included only within intervals of $\delta$ for which $G_{\delta }$ does not change.

First, we define a spur, which is a structure defined for a fixed volume parameter and, conceptually, is determined by noting the order in which various red one-cells have sufficient volume as height varies. It may be helpful to recall that $R$ is a tree made up of one- and two-cells.

See Figure 9. The following lemma gives us insight into the structure of spurs of $R$.

When $i=j$, the claim holds trivially. Then consider some $i\ne j$ and suppose that the path from $b_i$ to $t_i$ intersects the path from $b_j$ to $t_j$, i.e., there is some last two-cell, $p$, contained in both paths. Towards a contradiction, suppose $p$ is not either $t_i$ or $t_j$. Then we can denote $w_i^1$ and $w_i^2$ as the one-cells on the path from $b_i$ to $t_i$ just before and just after $p$ respectively: define $w_j^1$ and $w_j^2$ similarly. Since $p$ is the last two-cell in the intersection, $w_i^2\ne w_j^2$.

Without loss of generality, suppose the significance of $w_i^1$ is no higher than that of $w_j^1$. Note that $w_i^2$ and $w_j^2$ are both adjacent to $w_i^1$, and only $w_i^2$ can have a higher significance than $w_i^2$, otherwise, $w_i^1$ would be adjacent to $t_i$, contradicting $p\ne t_i$. But then both $w_j^2$ and $w_i^2$ would have higher significance than $w_j^2$, contradicting that $w_j^2$ is not adjacent to $t_j$. $◀$

Since a spur is a connected collection of one- and two-cells, we measure the volume of a spur using quarter-pillar volume, just as we did previously for other types of superlevel sets.

A spur with volume persistence of at least $\delta$ corresponds to a region of the terrain with enough volume to separate a channel on either side of it, at least up to its flanking height. Specifically, we place fingers along the boundaries of spurs, defined below.

We choose fingers that lie on spur boundaries; our final lemma helps justify this choice.

Suppose, towards a contradiction, that $w$ is a one-cell in the boundary of $S$, not adjacent to $t$, that does not lie on an ascending path or saddle of the terrain. Recall that every one-cell is either part of an ascending path, is a saddle, or is part of a descending path. I.e., it must be that $w$ is part of a descending path, and thus has a gradient pair to some adjacent two-cell, $q$. But, since $S$ is the entire subtree below its top, we have $q\in S$, which contradicts our assumption that $w$ is on the boundary of $S$. $◀$

We are now ready to compute fingers. We summarize our entire procedure as follows:

1. 1.

   Extend given elevation data to a DMF on all cells (Construction 1).

2. 2.

   Compute $G_{\delta }$, i.e., channels that separate hills of volume no less than $\delta$ (Section 3), for an appropriate choice of $\delta >0$ (in practice, this value is determined by domain experts).

3. 3.

   In each face of $G_{\delta }$, create a red tree (Definition 15).

4. 4.

   For such a red tree, determine its set of spurs (Definition 18).

5. 5.

   Consider the subset, ${𝒮}_{\delta }$, of spurs with volume persistence no less than $\delta$ (Definition 20).

6. 6.

   For each spur $S\in {𝒮}_{\delta }$, consider the two walks (clockwise and counterclockwise) around the boundary of $S$ starting at the top of $S$. In each direction, note the first zero-cell we encounter with height no more than the flanking height of $S$. Such zero-cells define the degree-one end of a finger, and we complete the finger by flowing down via the lowest lexicographical path (along the boundary of $S$, and then possibly through additional ascending paths and saddles) until we reach a minimum of $G_{\delta }$.

7. 7.

   Output $G_{\delta }$ along with all fingers as the resulting network.

We end with a remark that highlights connections between this section and Section 3.

Figure 10 shows the input and output for the Western Scheldt dataset. In the resulting network (Figure 10f) the color of each arc encodes its volume persistence: from dark purple (high volume) to bright pink (low volume). Arcs with volume persistence less than $\delta$ are not shown.

In practice, we found that the method is quite stable under small changes in the DEM. To demonstrate this stability, in Figure 11 we show the networks extracted from historical DEMs from 1999, 2004, 2009, and 2014. Though some areas of the estuary underwent rapid changes during this period, the other parts of the extracted network stay similar.

Figure 12 shows the input and output for the Wadden Sea dataset. This dataset spans the same area as depicted in Figure 1. Figure 12b shows the result of our network extraction method. The method detects one significant bar between two channels, but as expected, does not detect the other features, as they are all fingers. Figure 12d then shows the result of our finger extraction method. The spurs (Figure 12c) represent areas of sufficient volume between fingers, and the extracted fingers capture the tree-like structure of the tidal basin well. In most cases, the flanking height of each spur is a good choice for determining how high the starting point of each finger lies. However, some of the extracted fingers are very short because the flanking height is quite low. This may be alleviated by filtering out fingers in their entirety if their length is below a given threshold.

The networks found by the implementation show that our algorithm is indeed able to identify geomorphologically salient networks, including fingers, in low-relief terrains. Domain experts have confirmed the high quality of the networks and thereby validated our mathematical modeling.