Jacobi Set Simplification for Tracking Topological Features in Time-Varying Scalar Fields

Abstract

The Jacobi set of a bivariate scalar field is the set of points where the gradients of the two constituent scalar fields align with each other. It captures the regions of topological changes in the bivariate field. The Jacobi set is a bivariate analog of critical points, and may correspond to features of interest. In the specific case of time-varying fields and when one of the scalar fields is time, the Jacobi set corresponds to temporal tracks of critical points, and serves as a feature-tracking graph. The Jacobi set of a bivariate field or a time-varying scalar field is complex, resulting in cluttered visualizations that are difficult to analyze. This paper addresses the problem of Jacobi set simplification. Specifically, we use the time-varying scalar field scenario to introduce a method that computes a reduced Jacobi set. The method is based on a stability measure called robustness that was originally developed for vector fields and helps capture the structural stability of critical points. We also present a mathematical analysis for the method, and describe an implementation for 2D time-varying scalar fields. Applications to both synthetic and real-world datasets demonstrate the effectiveness of the method for tracking features.

1 Introduction

The relationship between scalar fields defined over a spatial domain of interest corresponds to interesting phenomena or captures application-specific features. For example, a water mass in oceanography is defined by the temperature-salinity (T-S) curve helland1916nogenTSDiag , fields such as pressure and wind speed are used to identify structural characteristics of a hurricane doleisch2004IsabelScivis . This motivates the study of techniques for visualizing multifield or multivariate data hansen2014scientific . Various structures have been introduced to analyze the relationship between multiple fields and specifically for bivariate fields, such as continuous scatterplot Bachthaler2008CSP , fiber surface carr2015fiber , Reeb space Edelsbrunner2008ReebSpace , Pareto set huettenberger2014decomposition , and Jacobi set edelsbrunner2004jacobi . We restrict our attention to the Jacobi set, a topological descriptor that is a generalization of the notion of critical points to the multivariate setting. For a bivariate field, the Jacobi set is defined as the set of points in the domain where the gradients of the two fields align with each other.

This paper focuses on time-varying scalar fields that may be represented as a bivariate field over space-time, where the time function is one of the two scalar fields. The Jacobi set of a time-varying scalar field represents the tracks of critical points of the scalar field over time. Critical points may correspond to interesting features, such as vortices in fluid flows or eddies in the ocean, and hence their tracks are of interest. The Jacobi set is, in this case, a feature-tracking graph. The Jacobi set computed from piecewise linear (PL) scalar functions is sensitive to noise, dense, and cluttered, making it difficult to analyze its structure. We propose a method that computes a simplified Jacobi set, resulting in a smaller number of clutter-free tracks while retaining important features that are represented by the original Jacobi set. The intuition from time-varying fields is central to the development of our simplification method and we consider this as a first step towards the development of methods that are applicable to generic bivariate fields.

We present a Jacobi set simplification method that considers both temporal and spatial coherence of features. It formulates critical points of a scalar field as zeros of the corresponding gradient field thereby enabling a bottom-up simplification directed by a stability measure called robustness that was introduced for vector fields Primoz2015RobustnessSimplification4VF . Clusters of critical points that are determined by the robustness parameter are tracked over time. The track of a cluster of critical points replaces the tracks of all critical points within the cluster thereby reducing the number of tracks while retaining their spatial and temporal structure. The resulting collection of tracks constitutes the simplified Jacobi set. Key contributions of this paper include

1. A robustness-based Jacobi set simplification algorithm for time-varying scalar fields,

2. A mathematical analysis supporting the algorithm that guarantees the existence of a corresponding simplified vector field, and

3. An implementation of the method for 2D time-varying scalar fields.

The effectiveness of the method is demonstrated via experiments on a synthetic dataset and multiple fluid flow datasets.

2 Related work

This paper focuses on simplifying the Jacobi set for time-varying scalar fields, resulting in reduced clutter in the visualizations of feature tracks. This section summarizes previous work on structures, including the Jacobi set, that have been employed for bivariate field visualization, different applications of the Jacobi set, methods for simplification, and feature tracking.

Multivariate field visualization. Different concepts and structures have been introduced to study the relationship between individual fields in a multivariate dataset. Continuous scatterplot (CSP) Bachthaler2008CSP , a generalization of scatterplot to continuous fields, provides a dense visualization in range space. This visual representation helps identify regions of interest that may in turn be mapped to the spatial domain using fiber surfaces carr2015fiber . CSPs and fiber surfaces have been used for studying bivariate and multivariate fields from multiple application domains Blecha2019Nuclear ; Raith2019Tensor ; sharma2021segmentation ; Sharma2023CSPOperators . Nagaraj and Natarajan Nagaraj2011isoExtractMultifield introduce a variation density function that captures the relationship between multiple scalar fields over isosurfaces of a given scalar field. This profile helps identify interesting isovalues of scalar fields in multivariate scenarios. Chattopadhyay et al. Chattopadhyay2014Jacobist introduce the Jacobi structures as the critical features in the Reeb space Edelsbrunner2008ReebSpace , establishing a relationship between the Jacobi set and singular fibers. The Reeb space is the bivariate analog of the Reeb graph. A fiber component, the bivariate counterpart of an isocontour component, maps to a point in the Reeb space. The joint contour net Carr2014JCN is a discrete representation of the Reeb space. The representation and computation of Reeb space is challenging, and hence it is often not computed explicitly in practice. The Jacobi set is a subset, specifically the image of singular fibers, of the Reeb space and is amenable to fast computation. The notion of Pareto set huettenberger2014decomposition helps identify consensus regions using Pareto optimality and dominance in multifields.

Jacobi set applications. The Jacobi set has been applied for many tasks in image analysis and visualization. Edelsbrunner et al. Edelsbrunner2004LnGcompJS introduce local and global comparison measures based on the Jacobi set to support comparative visualization. In image processing, Norgard and Bremer Norgard2013ridgeValley utilize the Jacobi set to extract image ridges, proposing a combinatorial algorithm to create a ridge-valley graph that captures information about all ridges in a given image. The Jacobi set has been used in geophysics for estimating relationships between geopotential height and total ozone column Artamonova2017physicsJS and in forestry for automatic tree ring detection Makela2020AutomaticTR . Tierny and Carr use the Jacobi fiber surfaces, the preimage of a Jacobi edge, to compute the Reeb space of a bivariate field tierny2016jacobi . Sharma and Natarajan sharma2022FF propose an output-sensitive algorithm for fiber surface computation that further facilitates the computation of individual fiber surface components in the vicinity of a Jacobi edge selected by the user. The performance of the algorithm depends on the size of the Jacobi set, and it will benefit from a simplified Jacobi set.

Jacobi set simplification. In practical scenarios, the Jacobi set consists of a complex network of edges, including noisy artifacts that result from the PL representation. Analyzing it can be challenging and will benefit from a simplification method. Simplifying the bivariate field will, indirectly, result in a simpler Jacobi set. The bivariate field may be simplified by employing topological simplification directed by persistence on the individual scalar fields ELZ02 ; bremer2007topological ; lukasczyk2020localized . However, in a time-varying scenario, such an approach may introduce temporal inconsistencies. Post simplification, a critical point may no longer have a correspondence in adjacent time steps, breaking the temporal continuity. The Jacobi set is not a stable structure. Small changes in the underlying bivariate field can produce significant changes in the Jacobi set. This instability poses challenges for such indirect simplification, especially in time-varying fields.

Direct simplification methods aim to directly remove noisy edges or sections of Jacobi curves that are deemed to be less significant. These approaches are directed by a metric that measures the significance of edges of the Jacobi set. Suthambhara and Natarajan Suthambhara2011jacobisetsimplification propose a method to reduce the number of connected components in the Jacobi set by posing the simplification problem as the extraction of level sets and offset contours. Bhatia et al. bhatia2015localjacobisetsimplification generalize the concept of critical point cancellation to Jacobi sets and simplify ’Jacobi regions’. The method identifies less significant portions of the Jacobi set and removes them by modifying the individual scalar fields locally. Recently developed methods klotzl2022local ; klotzl2022reduced compute a clutter-free visual representation of the Jacobi set by considering a bilinear interpolant. The method removes zig-zag artifacts while maintaining the topological structure.

Our proposed method may be classified under the direct simplification category. It reduces the number of Jacobi edges by clustering edges that potentially represent a common feature for a given approximation threshold. It directly computes a simplified version of the Jacobi set as opposed to simplifying the noisy Jacobi set.

Feature tracking. In a time-varying scalar field, the Jacobi edges represent tracks of critical points. Since critical points are good representatives of topological features, critical point tracking is a preferred approach for topological feature tracking. Feature tracking has been studied widely, but it has not been formulated using the Jacobi set except for a few instances edelsbrunner2008time . Reininghaus et al. reininghaus2011efficient use ideas from discrete Morse theory and combinatorial feature flow fields to track critical points. They introduce the notion of integrated persistence, combining spatial persistence of critical points with their temporal evolution to determine the significance of a particular track as opposed to direct elimination of short-lived features. Other works lukasczyk2017nested ; saikia2017global ; Lukasczyk2019DNTG ; Yan2021Survey focus on tracking regions of interest, such as predefined connected components instead of critical points. Feature tracking has also been studied for vector fields TRICOCHE2002topologytracking ; Post2003SOFlowVis ; theisel2003feature ; Tino2011StableFeatureFF ; Bujack2020TimeFlowSTA .

Our method is inspired by the notion of tracking graphs that represent the correspondences and tracks of individual connected components of spatial features. We formulate the critical point tracking problem as tracking zeros of a time-varying vector field. To achieve this, we use intuition, definitions, and analyses from previous work on tracking critical points in vector fields using robustness chazal2011computing ; Primoz2015RobustnessSimplification4VF .

Figure 1: Jacobi set of a 2D synthetic bivariate field. (a) Scalar field $f(x,y)=x^{2}+y^{2}$ (b) Scalar field $g(x,y)=x(y-8)$ (c) Level sets of $f$ (green) and $g$ (red). The Jacobi set $\mathbf{J}(f,g)$, shown in black, passes through the points where gradients of $f$ and $g$ align.

3 Background

This section presents the necessary definitions, notations, and mathematical preliminaries for describing Jacobi set simplification. It introduces critical points of a scalar field and the notion of Jacobi set for a bivariate and a time-varying scalar field, both in the smooth and PL setting edelsbrunner2022computational .

Scalar field. A scalar field $f:\mathbb{M}\rightarrow\mathbb{R}$ is a real-valued function defined on an $n$-dimensional domain $\mathbb{M}$ and maps a point on the domain to a single scalar value. The scalar values usually represent physical quantities such as pressure, temperature, speed, or height, although they may also be synthetically generated. Figure 1(a) shows a 2D scalar field, $f(x,y)=x^{2}+y^{2}$. A point $p\in\mathbb{M}$ is a critical point of $f$ if and only if the gradient of $f$ at $p$ vanishes. Else, $p$ is a regular point. The scalar value $f(p)$ is called a critical or regular value, respectively. A critical point $p$ is said to be non-degenerate if the Hessian at $p$ is non-singular. If all the critical points of $f$ are non-degenerate and have distinct values, then $f$ is called a Morse function. A key result from Morse theory states that all smooth functions are either Morse or can be perturbed using an infinitesimally small quantity into a Morse function. All scalar functions considered in this paper are assumed to be Morse functions.

Bivariate field. A multivariate function or a multifield is a collection of scalar functions / scalar fields defined over a common spatial domain. A bivariate field $\{f,g\}:\mathbb{M}\rightarrow\mathbb{R}^{2}$ is a specific case where two scalar fields are defined over a domain $\mathbb{M}$. Bivariate field analysis becomes beneficial when important data features depend on the interaction between two individual fields as opposed to being characterized independently by the two fields. Figure 1 shows two scalar fields that can either be studied independently or as a bivariate field $\{f,g\}$.

Jacobi set. Given a bivariate field $\{f,g\}:\mathbb{M}\rightarrow\mathbb{R}^{2}$ such that the intersection of the sets of the critical points of $f$ and $g$ is a null set, the Jacobi set $\mathbb{J}=\mathbb{J}(f,g)=\mathbb{J}(g,f)$ is the set of points where gradients of $f$ and $g$ are linearly dependent.

$$\begin{split}\mathbb{J}(f,g)=\{x\in\mathbb{M}\ |\ \nabla f+\lambda\nabla g=0\text{ or}\\ \lambda\nabla f+\nabla g=0\}\end{split}$$ (1)

for some $\lambda\in\mathbb{R}$. In other words,

$$\begin{split}\mathbb{J}(f,g)=\{x\in\mathbb{M}\ |\ \text{$x$ is a critical point of }f+\lambda g\text{ or}\\ \lambda f+g\}\end{split}$$ (2)

Equation 2 essentially characterizes the Jacobi set as the set of critical points of a collection of scalar fields parameterized by $\lambda\in\mathbb{R}$. Let $k\in\mathbb{R}$ be a regular value of $g$. The scalar field $f_{k}:g^{-1}(k)\rightarrow\mathbb{R}$ is called the restriction of the scalar field $f$ to the level set $g^{-1}(k)$. The Jacobi set is equivalently defined as

$$\mathbb{J}(f,g)=cl(\{p\in\mathbb{M}\ |\ p\mbox{ is a critical point of }f_{k},k\in\mathbb{R}\})$$ (3)

The Smooth Embedding Theorem edelsbrunner2004jacobi states that, under generic conditions, the Jacobi set of two Morse functions is a smoothly embedded $1$-manifold in $\mathbb{M}$. Figure 1(c) shows the level sets of two scalar fields $f$ and $g$ using isocontours. The Jacobi set (black) of the bivariate field $\{f,g\}$ is a curve on the plane. The individual scalar fields are shown in Figure 1(a) and Figure 1(b).

Figure 2: Identifying a Jacobi edge. (Left) $h_{\lambda_{e}}(v_{1})<h_{\lambda_{e}}(a)$ and $h_{\lambda_{e}}(v_{2})>h_{\lambda_{e}}(a)$, which implies that $ab$ is a regular edge. (Middle) Both $h_{\lambda_{e}}(v_{1})$ and $h_{\lambda_{e}}(v_{2})$ are smaller than $h_{\lambda_{e}}(a)$, which implies that $ab$ is a maximum. (Right) Edge $ab$ is a minimum Jacobi edge.

Piecewise linear function. In practice, the scalar fields are often available as a sample over a spatial domain. The spatial domain is represented as a mesh, either a structured grid or a triangulation. The samples of the scalar field are available at the vertices of the mesh. A piecewise linear (PL) approximation of the scalar field is obtained via linear interpolation in the interior of cells of all dimensions.

Let $f:|K|\rightarrow\mathbb{R}$ be a PL function defined on a triangulation $K$ of a manifold $\mathbb{M}$. Critical vertices of $f$ are characterized by the Betti numbers of its lower link delfinado1993incremental ; edelsbrunner2022computational . Non-Morse functions can be transformed into Morse functions by addressing degeneracies through the simulation of simplicity edelsbrunner1990simulation . In the following, we always assume that the scalar field is a Morse function.

The Jacobi set of two PL functions may be computed as the set of critical points of the family of functions $h_{\lambda}=f+\lambda g$, following Equation 2. According to the Critical Edge Lemma edelsbrunner2004jacobi , the Jacobi set lies on the edges of the triangulation $K$. If an edge $e=(a,b)\in K$ belongs to the Jacobi set then the function $h_{\lambda}$ is necessarily flat along the edge, i.e., $h_{\lambda_{e}}(a)=h_{\lambda_{e}}(b)$. This value $\lambda_{e}$ for which $h_{\lambda}$ is flat can be computed as $\lambda_{e}=(f(b)-f(a))/(g(a)-g(b))$. Further, the edge $e$ is classified as critical based on a characterization of the lower link of $e$ with respect to $h_{\lambda_{e}}$. Figure 2 shows the categorization of an edge $ab$ as regular, maximum, or minimum. If $h_{\lambda_{e}}(v_{1})<h_{\lambda_{e}}(a)$ and $h_{\lambda_{e}}(v_{2})>h_{\lambda_{e}}(a)$ then the edge $ab$ is a regular edge. If both $h_{\lambda_{e}}(v_{1})$ and $h_{\lambda_{e}}(v_{2})$ are lower or higher than $h_{\lambda_{e}}(a)$, then $ab$ is an extremum edge. The Jacobi set is obtained by stitching the critical edges at their endpoints. The resulting set is a 1-manifold due to the Even Degree Lemma, which states that the degree of each vertex in the Jacobi set is even.

Figure 3: Jacobi set of a time-varying scalar field. (a) A synthetic 2D time-varying scalar field obtained by applying incremental rotations over time. Select time steps within the range [1..299]. (b) All time steps of the scalar field. stacked together to show the spiral path followed by the two pairs of maxima. (c) The Jacobi set of the time-varying scalar field is noisy. It is difficult to identify the two primary data features due to clutter. (d) Edges of the Jacobi set filtered based on the scalar value. Edges with endpoints having scalar value less than $1$ are discarded, leaving only a few Jacobi edges that are in the spatial proximity of important features.

Time-varying scalar field. Let $f:\mathbb{M}\times\mathbb{R}\rightarrow\mathbb{R}$ be a time-varying scalar function. The second parameter of $f$ is considered as the time axis. Consider a time function $g:\mathbb{M}\times\mathbb{R}\rightarrow\mathbb{R}$ defined as $g(p,t)=t$. Following the formulation in Equation 3, the Jacobi set $\mathbb{J}(f,g)$ is the collection of tracks of critical points of $f$ at every time step. This is equal to the set of critical points of $f$ at every time step and the temporal edges that represent the correspondences between critical points lying on two consecutive time steps.

In this paper, we consider 2D PL time-varying scalar fields $\{f_{t}\}$, where $f_{t}$ is a 2D PL scalar field defined in time step $t$. The domain $\mathbb{M}\times\mathbb{R}$ is triangulated by extending a triangulation of the spatial domain. We include temporal edges between corresponding vertices in two consecutive time steps and insert diagonals to decompose the resulting prism into three tetrahedra. The Jacobi set $\mathbb{J}(\{f_{t}\})$ of the time-varying scalar field is computed over this triangulated domain using the standard PL algorithm edelsbrunner2004jacobi .

Figure 3(a) shows a subset of $300$ time steps of a synthetic time-varying scalar field. The scalar field contains two pairs of Gaussians undergoing rotation over time. In Figure 3(b), all time steps are stacked together, highlighting the rotational pattern formed by these pairs. Figure 3(c) shows the Jacobi set, considering time as the second field. However, the Jacobi set is cluttered, and the noisy Jacobi edges occlude the prominent rotational pattern. Figure 3(d) shows Jacobi edges close to prominent features by applying a manual threshold. Edges with endpoints having scalar field value less than $1$ are discarded. In this paper, our goal is to simplify the Jacobi set such that the simplified Jacobi set captures the important data patterns.

4 Robustness-based simplification of the Jacobi set

In this section, we introduce the notion of a simplified Jacobi set for time-varying scalar fields based on the robustness of critical points of a vector field. Critical points of a scalar field occur where the corresponding gradient field vanishes. Extremum critical points of the scalar field correspond to sources and sinks in the gradient field. A saddle corresponds to a saddle in the associated gradient field. Hence, the Jacobi set, consisting of the tracks of critical points of a time-varying scalar field, is the same as the track of zeros of the associated gradient field. We first briefly introduce the notion of robustness of critical points in a vector field and then describe our adaptation to gradient fields followed by a method for Jacobi set simplification.

Robustness. In a 2D vector field, the critical point is a point where the vector vanishes and is characterized by its Poincaré index. Sources and sinks have an index of +1, and a saddle has an index of -1. Figure 4(a) shows a 2D synthetic vector field. A subset of its critical points, source $M$ (red), sinks $m_{i}$ (blue), and saddles $s_{i}$ (green) are highlighted. Robustness is a measure that quantifies the stability of a critical point in a vector field. Intuitively, it represents the difficulty or effort required for perturbing away the critical point.

The formal definition of robustness depends on a topological structure called the merge tree. The merge tree of a scalar field $f$ is an abstract representation of the evolution of sublevel set components of $f$. Leaves of the merge tree correspond to the minima of the underlying scalar field and represent the birth of a component. The internal nodes correspond to saddle critical points and represent the merging of sublevel set components.

Consider a 2D vector field $V:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$. The robustness of a critical point of $V$ is determined by the critical points that lie within a specific spatial neighborhood of the critical point. This spatial neighborhood is in turn determined by the merge tree of the vector magnitude function $\|V\|$ defined on the plane. A connected component of the sublevel set of $\|V\|$ contains one or more nodes of the merge tree $\mathcal{MT}(\|V\|)$, see Figure 4. The collection of nodes contained within $C$ together with their incident arcs constitutes a subtree of $\mathcal{MT}(\|V\|)$. Given such a path connected component $C$, the degree $deg(C)$ is defined as the sum of indices of critical points of the vector field $V$ that are enclosed within the component. The degree $deg(C)$ is assigned to the root of this subtree. Note that critical points of the vector field appear as leaf nodes (minima) of $\mathcal{MT}(\|V\|)$. The degree of these leaf nodes is equal to their Poincaré index. The degree of an internal node of $\mathcal{MT}(\|V\|)$ is equal to the sum of degrees of all leaves of the subtree rooted at the node.

The static robustness $\delta$ of a critical point is defined as the vector magnitude of the lowest ancestor node in the merge tree with degree 0 chazal2011computing . This paper utilizes only the notion of static robustness and we often refer to it as robustness, without the qualifier. In the vector field shown in Figure 4, the source $m_{4}$ has a lower robustness value as compared to the other sources. Note that static robustness of zeros in vector fields is theoretically different from the concept of persistence of critical points in scalar fields chazal2011computing ; edelsbrunner2011quantifying ; Primoz2015RobustnessSimplification4VF . We discuss the difference between the two approaches towards Jacobi set simplification in section 6.

Figure 4: (a) A synthetic 2D vector field visualized using a LIC image and annotated with a subset of its critical points: sources (blue), sink (red), and saddles (green). Isocontours of the vector magnitude field are shown in different colors. The corresponding sublevel sets contain subsets of critical points. All contours corresponding to the sublevel set of a critical point may not be visible. Contours at smaller isovalues may lie close to the critical point and hence are occluded by the glyph used to show the critical point. (b) The merge tree captures the evolution of the connected components of the sublevel set components. Nodes are annotated with the degree. The degree of the leaf nodes is equal to their Poincaré index, +1 for sources/sinks and -1 for saddles. The degree of an internal node is equal to the sum of degrees of all leaf nodes within the subtree rooted at the node.

Jacobi set and critical point tracks. Given a 2D PL time-varying scalar field $\{f_{t}\}$, where $f_{t}:\mathbb{R}^{2}\rightarrow\mathbb{R}$ is the scalar field in the $t^{th}$ time step, we derive the Jacobi set $\mathbb{J}(\{f_{t}\})$ by introducing the corresponding time-varying gradient field $\{\nabla f_{t}\}$, where $\nabla f_{t}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$. The critical points of $f_{t}$ are equal to the points where $\nabla f_{t}$ vanishes. Minima, maxima, and saddles of the scalar field behave as sources, sinks, and saddles of the corresponding gradient field, respectively. These critical points are the 0-value minima of the gradient magnitude field $G_{t}=\|\nabla f_{t}\|$. Hence, the set of tracks of critical points of $\{f_{t}\}$ is the same as the set of tracks of zeros of $\{G_{t}\}$. We denote the original Jacobi set of the time-varying scalar field either as $\mathbb{J}(\{f_{t}\})$ or $\mathbb{J}(\{G_{t}\})$.

Simplifying the Jacobi set. A $\delta$-sublevel set of the gradient magnitude field $G_{t}$ is the preimage $G_{t}^{-1}(-\infty,\delta]$. Let $\mathcal{C}(G_{t},\delta)$ denote the set of connected components of the $\delta$-sublevel set. Each connected component in $\mathcal{C}(G_{t},\delta)$ contains a collection of critical points of the corresponding gradient vector field $\nabla f_{t}$, similar to the illustration in Figure 4 for a generic vector field. These critical points correspond to critical points of the scalar field $f_{t}$ for all $\delta\geq 0$ and for all time $t$. Consequently, all edges of the Jacobi set $\mathbb{J}(\{f_{t}\})$ pass through the collection $\mathcal{C}(G_{t},\delta)$ for all $t$. The collection of edges of $\mathbb{J}(\{f_{t}\})$ passing through a pair of spatially overlapping components from adjacent time steps may be represented by a single Jacobi edge, resulting in a reduction in the size of the Jacobi set. Figure 1 in the supplementary material illustrates the simplification. Increasing $\delta$ results in larger components of $\mathcal{C}(G_{t},\delta)$, which in turn implies a larger number of critical points clustered together into a single component and represented by a single Jacobi edge. This implies that $\delta$ is a natural parameter for driving Jacobi set simplification.

Given a vector-valued function $V:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$, the function $\widetilde{V}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ is called a $\delta$-perturbation of $V$ if

$$\sup_{p\in\mathbb{R}^{2}}\lVert V(p)-\widetilde{V}(p)\rVert\leq\delta.$$

The notion of $\delta$-perturbation aids in the numerical realization of the process that removes critical points with low robustness. Specifically, it quantifies the effort required to remove the critical points. We transport the notion of robustness to time-varying scalar functions and to their Jacobi set by studying and characterizing $\delta$-robust critical points of the associated gradient vector field.

A valid simplification $\mathbb{J}^{*}$ of $\mathbb{J}(\{G_{t}\})$ is a Jacobi set that is derived from $\mathbb{J}(\{G_{t}\})$ via a $k\delta$-perturbation of $\nabla f_{t}$ and satisfies the following conditions:

1. $\mathbb{J}^{*}=\mathbb{J}(\{\widetilde{G}_{t}\})$, where $\widetilde{G}_{t}$ is the gradient magnitude field of a $k\delta$-perturbation ($k$ is a small positive integer) of $\nabla f_{t}$ that consists of only $\delta$-robust critical points.

2. $\mathcal{C}(G_{t},\delta)=\mathcal{C}(\widetilde{G}_{t},\delta)$ for all $t$.

3. $\phi_{t}(\widetilde{G}_{t})(x)\leq\phi_{t}(G_{t})(x)$ for all $t$, where the function $\phi(G_{t})$ maps a pair of $\delta$-sublevel set components to the number of Jacobi edges between them, $\phi_{t}(G_{t}):\mathcal{C}(G_{t},\delta)\times\mathcal{C}(G_{t+1},\delta)\rightarrow\mathbb{Z}^{+}\cup\{0\}$.

The time-varying gradient magnitude field $G_{t}$ and its perturbation $\{\widetilde{G}_{t}\}$ are close to each other, resulting in a Jacobi set $\mathbb{J}^{*}$ that is close to $\mathbb{J}(\{G_{t}\})$. The first condition also ensures that noisy critical points, that are not $\delta$-robust, are perturbed away and the stable critical points form the simplified Jacobi set. The second condition ensures that $G_{t}$ is not altered significantly and the geometry of $\mathbb{J}(\{G_{t}\})$ is closely approximated. Finally, the third condition ensures a reduction in the number of Jacobi set edges by considering each pair of $\delta$-sublevel set components from adjacent time steps that contributes an edge to the Jacobi set.

5 Algorithm

We now present a method to compute a valid simplification of a given Jacobi set $\mathbb{J}(\{G_{t}\})$. The input scalar field $\{f_{t}\}$ at every time step is defined on the vertices of a grid and is interpolated within the interior of cells. The algorithm consists of the following three steps:

Gradient computation. A gradient field is computed for each PL function $f_{t}$ via central differences. At a grid vertex, the partial derivative along the X-axis is computed as the difference between the scalar values of its left and right neighbors. Similarly, neighbors along the Y-axis are used to compute the partial derivative along the Y-axis. The two partials constitute the gradient vector. The sequence of gradient fields $\{\nabla f_{t}\}$ and the corresponding gradient magnitude fields $\{G_{t}\}$ are computed in the first step. For smooth functions, the critical points of a scalar field $f_{t}$ correspond to a zero of the gradient vector field $\nabla f_{t}$. However, such a clear correspondence may not exist for a PL scalar field. So, we assign a zero value for $\nabla f_{t}$ and $G_{t}$ at vertices that are identified as critical points of $f_{t}$. $G_{t}$ is evaluated at all grid vertices and linearly extended in the interior of the cells.

$\delta$-sublevel set computation. In the second step, we compute the connected components of the $\delta$-sublevel sets of $\{G_{t}\}$. section 7 provides the implementation details.

Simplified Jacobi set computation. The simplified Jacobi set is computed as the tracking graph of the $\delta$-sublevel set components of $\{G_{t}\}$. Correspondence between two $\delta$-sublevel set components from adjacent time steps is determined via spatial overlap. If $C_{1}$ and $C_{2}$ are two components that have a non-empty spatial overlap and further the degrees $deg(C_{1})\not=0$ and $deg(C_{2})\not=0$, then we insert an edge of the Jacobi set between $C_{1}$ and $C_{2}$. A geometric embedding of this edge is computed by using the centroid of the component as a representative and inserting the edge connecting the two representatives. If the component $C_{1}$ from the earlier time step has a non-empty overlap with multiple components from the next time step, then we choose the component with the largest overlap in terms of area.

The Jacobi set is hence computed as a collection of edges between representatives of $\delta$-sublevel set components. The resulting set may contain short or broken tracks, due to short-lived critical points caused by noise. A post-processing step addresses this issue as discussed in section 7.

6 Analysis

In this section, we present a theoretical analysis of the algorithm based on a study of the stability of the zeros in the gradient field. We begin by stating two results regarding Riemann maps conway2012functions ; conway2012functions2 .

Theorem 6.1 (Riemann Mapping Theorem conway2012functions )

Let $A\subset\mathbb{C}$ be a simply connected region which is not the whole plane and let $a\in A$. There exists a unique analytic function $R:A\rightarrow\mathbb{C}$ satisfying the following properties:

• $R(a)=0$ and the derivative $R^{\prime}(a)>0$

• $R$ is injective

• $R(A)=\mathbb{D}$, the unit open disk in the plane

Theorem 6.2 ( conway2012functions2 )

If $A$ is a bounded simply connected region and $R^{-1}:\mathbb{D}\rightarrow A$ is the Riemann map with $R^{-1}(0)=0$ and $(R^{-1})^{\prime}(0)>0$, then $R^{-1}$ extends to a homeomorphism of $\overline{\mathbb{D}}$ onto $\overline{A}$ iff the boundary of $A$, $Bd(A)$, is a Jordan curve. Further, $Bd(A)$ is homeomorphic to $Bd(\mathbb{D})$.

Given a connected component of a $\delta$-sublevel set of the gradient magnitude field $G=\|\nabla f\|$ whose degree is non-zero, we will use the Riemann mapping theorem to show the existence of a simplified vector field. This simplified vector field will contain a single zero within each component of the $\delta$-sublevel set $\mathcal{C}(G,\delta)$. We  prove the following result:

Theorem 6.3

Let $X\in\mathcal{C}(G,\delta)$ be a simply connected component of the $\delta$-sublevel set of the gradient magnitude field G, with non-zero degree and its boundary is a Jordan curve. There exists a continuous vector field $\widetilde{V}:\overline{X}\rightarrow\mathbb{R}^{2}$ that satisfies the following properties:

• $\widetilde{V}$ is a $2\delta$-perturbation of $\nabla f_{|\overline{X}}$,

• $\widetilde{V}_{|Bd(X)}=\nabla f_{|Bd(X)}$, and

• $\widetilde{V}$ has exactly one singularity.

Proof

We first reduce the problem onto a simpler domain, namely the disk, rather than dealing with a generic simply connected Jordan region $X$. The proof consists of two steps.

Step 1: Reduce the problem to a disk. Theorem 1 and 2 guarantee the existence of a homeomorphism $R:\overline{X}\rightarrow\overline{\mathbb{D}}$, the extended Riemann map. Now, define a vector field on the closed unit disk $U:\overline{\mathbb{D}}\rightarrow\mathbb{R}^{2}$ as a composition $U=\nabla f_{|\overline{X}}\circ R^{-1}$, which allows us to transfer a simplified vector field designed within $\overline{\mathbb{D}}$ onto $\overline{X}$ as shown in Figure 5.

Step 2: Construct a simplified vector field on the unit disk. We construct a continuous vector field $\widetilde{V}_{D}:\overline{\mathbb{D}}\rightarrow\mathbb{R}^{2}$ on the closed disk as $\widetilde{V}_{D}(r,\theta)=rU(1,\theta)$. This function is composed with $R$ to give $\widetilde{V}=\widetilde{V}_{D}\circ R:\overline{X}\rightarrow\mathbb{R}^{2}$, see Figure 5

Figure 5: A simplified gradient vector field is constructed within the unit disk $\overline{\mathbb{D}}$ and then transferred to the region $\overline{X}$.

The function $\widetilde{V}$ satisfies the three properties from Theorem 3. Since $U$ has constant magnitude $\delta>0$ on $Bd(X)$ and hence does not have any zeros, it follows that $\widetilde{V}_{D}$ and $\widetilde{V}$ have exactly one singularity. Further, it is evident from the construction that $\widetilde{V}_{D|Bd(\mathbb{D})}=U_{|Bd(\mathbb{D})}$ and hence $\widetilde{V}$ and $\nabla f_{|\overline{X}}$ are equal on $Bd(X)$. Finally, since $X$ is a $\delta$-sublevel set, the following is true for all values of $0\leq r\leq 1$ and $0\leq\theta\leq 2\pi$

$$\|U(1,\theta)\|=\delta\implies\|rU(1,\theta)\|\leq\delta\implies\|\widetilde{V}_{D}\|\leq\delta$$

Applying the triangle inequality, we have

$$\|\widetilde{V}_{D}(r,\theta)-U(r,\theta)\| \leq 2\delta,\mbox{and}$$

$$\|\widetilde{V}(x)-\nabla f(x)\| \leq 2\delta$$

Hence, the simplified vector field $\widetilde{V}$ within $X$ is a $2\delta$-perturbation of the gradient field, and its magnitude is smaller than $\delta$.

Theorem 3 guarantees the existence of a simplified continuous vector field that differs from the gradient of the input scalar field only within the $\delta$-sublevel set components. Further, each zero of the simplified vector field within a $\delta$-sublevel set component is at least $\delta$-robust.

Next, we consider the case where the connected component of the $\delta$-sublevel set has degree-0. We use an adaptation of a theorem developed in the study of robustness for vector fields chazal2011computing to show that the gradient field can be simplified within such components to remove all critical points.

Theorem 6.4 ( chazal2011computing )

Let $X\in\mathcal{C}(G,\delta)$ and let $\nabla f$ have degree-0 in $X$. There exists a continuous vector field $\widetilde{V}:\overline{X}\rightarrow\mathbb{R}^{2}$ that satisfies the following properties:

• $\widetilde{V}$ is a $\delta$-perturbation of $\nabla f_{|\overline{X}}$,

• $\widetilde{V}_{|Bd(X)}=\nabla f_{|Bd(X)}$, and

• $\widetilde{V}$ has no singularity.

Theorem 4 implies that perturbing $\nabla f$ within all such regions $X$ results in a simplified continuous vector field that preserves the $\delta$-sublevel sets. The interpretation of static robustness of a critical point as its stability or the perturbation required to remove it also follows from Theorem 4.

The above results guarantee that the zeros of $\nabla f$ can be perturbed away to retain either a single stable zero (Theorem 3) or remove all zeros (Theorem 4). The perturbation results imply that the zeros (critical points of $f$) within the $\delta$-sublevel set are less significant, need not be tracked individually, and can be replaced by a single stable zero. This corollary drives our simplification algorithm. The algorithm tracks connected components of $\mathcal{C}(G_{t},\delta)$ across time $t$ if their degree is non-zero. A single edge replaces the collection of Jacobi edges that are incident on the critical points lying within these $\delta$-sublevel set components. The simplified Jacobi set computed by our algorithm adheres to the valid simplification constraints outlined in section 4 and represents the tracks of zeros of the perturbed vector fields described in Theorems 3 and 4. However, the perturbed fields are not computed explicitly.

The algorithm reduces the number of edges and hence computes a simplified Jacobi set. However, it captures the macroscopic trends of movement of the critical points. Birth and death of less significant critical points do not affect the result and do not appear in the simplified Jacobi set because they occur in the interior of $\delta$-sublevel sets. The $\delta$-sublevel sets provide a clustering of critical points that are in geometric proximity to each other and evolve in a coherent manner over time. Thus, the algorithm ensures temporal and spatial consistency of the simplified Jacobi set while reducing the number of tracks. In contrast, tracking methods that rely on standard persistence-based filtering do not require persistence pairs to be consistently paired over time, nor do the pairs have to be in geometric proximity to each other. Also, the persistence values of critical points may vary significantly over time.

7 Implementation

This section describes details of an implementation of the Jacobi set simplification algorithm, optimizations that improve the runtime performance and the visual clarity of the simplified Jacobi set. All visualizations in this paper are generated using ParaviewAyachit2015paraview and TTK Tierny2018ttk . In the following discussion, $T$ denotes the total number of time steps and $n$ is the total number of vertices in the input 2D grid.

Persistence simplification. Simulation of simplicity, a symbolic perturbation of the scalar field, is used to handle degeneracies when classifying critical points. However, this results in the creation of spurious critical points within regions where the scalar function is flat. These critical points of $f_{t}$ are removed within each time step using a persistence-driven topological simplification ELZ02 . A small persistence threshold is sufficient to remove these critical points. This helps reduce unnecessary processing of these spurious critical points in subsequent steps.

Critical point computation. The critical points and their types (maxima, minima, or saddle) are computed and stored. This requires $O(nT)$ time. Storing the type of critical point helps optionally filter the Jacobi set to show edges corresponding to a specific type.

Gradient field computation. Computing the gradient magnitude field $G_{t}$ at each time step $t$ and explicitly setting the gradient magnitude of critical points to zero requires a runtime of $O(nT)$.

Figure 6: $\delta$-sublevel set computation. The $\delta$-isocontour (red) of the gradient magnitude field has two components and the $\delta$-sublevel set contains three critical points $a,b$, and $c$. The algorithm returns the collection of blue cells, a conservative approximation of the $\delta$-sublevel set.

$\delta$-sublevel set computation. We compute the $\delta$-sublevel sets at each time step $t$ in a single pass over the list of all vertices. If $G_{t}(p)\leq\delta$ for a vertex $p$ then all cells containing $p$ are inserted into the reported $\delta$-sublevel set. This reported $\delta$-sublevel set is a superset of the true $\delta$-sublevel set. The additional cells ensure that a connected component consisting of a single point is also represented by a spatial region and may be tracked via spatial overlap. Figure 6 illustrates the computation with a small example consisting of two components and three critical points. The running time for this step is $O(nT)$.

Connected component computation. Connected components of the $\delta$-sublevel sets are computed using the connectivity filter of TTK. Further, all degree-0 components are removed. The filter uses a flood fill algorithm and has a runtime of $O(nT)$. Removing all components with degree 0 also takes $O(nT)$ time.

Simplified Jacobi set. The spatial overlap based tracking filter in TTK, called “Tracking from Overlap” lukasczyk2017nested , is used to compute a tracking graph. Edges of the simplified Jacobi set are essentially the edges of this tracking graph. This computation has a worst case runtime of $O(nT)$. Further details are available in Section 1 of the supplementary material.

Post-processing. Noise in the data could result in short-lived or broken tracks. We identify and resolve such issues in a post-processing step to produce the simplified Jacobi set. This step analyzes the proximity of the start and end points of different tracks and determines if they may be merged. Specifically, it checks two criteria as illustrated in Figure 7. First, for a track $\tau_{1}$ that begins at time $t$ and at spatial location $(x,y)$, a search within the time interval $[t-\epsilon_{t},t]$ determines if there exists another track whose end point lies within an $\epsilon_{s}$-neighborhood of $(x,y)$ i.e., within a ball of radius $\epsilon_{s}$ centered at $(x,y)$. If such a track $\tau_{2}$ is found, its endpoint is connected to the start of $\tau_{1}$, as shown in Figure 7(a). If multiple such tracks are found, the one with the closest end time is declared as $\tau_{2}$. Second, for a track $\tau_{1}$ that ends at time $t$ and at spatial location $(x,y)$, a search is performed to determine if there exists another track that ends later than time $t$ and whose start time is within the interval $[t-\epsilon_{t},t]$ and start point lies within an $\epsilon_{s}$-neighborhood of $(x,y)$. If such a track $\tau_{2}$ is found, the end point at time step $t$ in $\tau_{1}$ is connected to the start point at time step $t+1$ in $\tau_{2}$ and the segment of $\tau_{2}$ up to time $t+1$ is discarded. Figure 7 (b) illustrates this case. Significant critical points are expected to correspond to long tracks. Based on this intuition, a length threshold $l$ is applied to discard all short tracks.

In the above implementation, a total of 5 parameters may be tuned to achieve an appropriately simplified Jacobi set. A persistence threshold $\epsilon_{p}$, the parameter $\delta$, the maximum length $\epsilon_{t}$ of recoverable temporal breaks in tracks, the size $\epsilon_{s}$ of the spatial neighborhood used to merge tracks, and a track length threshold $\epsilon_{l}$. Table 1 in the supplementary material lists the parameter values used in our experiments. Combining the running time of the individual steps, the total runtime of the algorithm is $O(nT)$.

Figure 7: Connecting broken tracks based on spatial and temporal proximity. (a) Tracks $\tau_{2}$ and $\tau_{3}$ end within a time interval $\epsilon_{t}$ of the start of track $\tau_{1}$ and their end points lie within an $\epsilon_{s}$-ball (orange) centered at the start point of $\tau_{1}$. The temporally closer track $\tau_{2}$ is merged with $\tau_{1}$ (green edge). (b) End time of track $\tau_{1}$ start time of $\tau_{2}$ are closer than the spatial ($\epsilon_{s}$) and temporal ($\epsilon_{t}$) thresholds. They are merged (green edge) and the initial segment of $\tau_{2}$ (dashed) is discarded.

8 Experimental results

In this section, we present results on four datasets – a synthetic rotating Gaussians dataset, 2D von Kármán vortex street, Boussinesq flow, and sea surface height. We provide qualitative insights into the reduction of the number of Jacobi edges and showcase the benefit of identifying prominent patterns and tracking features corresponding to critical regions. Table 1 in the supplementary material lists the statistics on the size of the simplified Jacobi set, demonstrating the effectiveness and utility of the method in practice. It also lists the parameter values used in the experiments for completeness.

The simplified Jacobi sets reported here and in subsequent results include tracks of all critical point types – maxima, minima, and saddles. The implementation may be extended to filter and selectively display tracks that correspond to a specific critical type. All tracks are rendered using cubic spline interpolation. The focus of the prototype implementation in Python is functionality and not efficiency. Code optimization using efficient libraries is likely to result in significant runtime improvement.

Rotating Gaussians. A synthetic dataset generated as a time-evolving sum of Gaussians helps analyze the correctness of the algorithm and its implementation. Two pairs of maxima, each consisting of two Gaussian centers (maxima) rotate about the origin over time as shown in Figure 8. The Jacobi set is noisy and cluttered as discussed earlier. We observe that upon filtering the Jacobi set (orange) to retain edges with scalar value greater than $1$, we obtain the spiral pattern followed by the rotating Gaussians. Tracks in the simplified Jacobi set shown in Figure 8(c) effectively capture this pattern, indicating that the simplified Jacobi set is a geometrically accurate representation of the feature (Gaussian center) evolution. Furthermore, a single track for each critical region shows the successful clustering of critical points. The straight track in the center corresponds to the global minimum. No post-processing was required to produce the results.

Figure 8: (a) A time-varying scalar field represented as a stack of color-mapped images. The vertical direction corresponds to time. (b) Edges of the original Jacobi set (orange) after filtering out edges where the scalar value is smaller than 1. (c) The simplified Jacobi set consists of three components. Two tracks capture the spiral pattern of rotating Gaussians. The straight track corresponds to the global minimum.

Figure 9: (a) Three time steps from the 2D von Kármán vortex street dataset color-mapped with $z$-component of vorticity. (b) Jacobi set is large and noisy. (c) Simplified Jacobi set is clutter-free. The tracks capture the flow pattern of vortices from left to right.. (d) A zoomed-in view shows that each track represents an individual critical region around a maximum or minimum, namely a vortex. The vertical track represents a critical region near the domain boundary and corresponds to a saddle.

2D von Kármán vortex street. This dataset of the flow of a fluid around a cylindrical obstacle is generated using Gerris flow solver gerrisflowsolver ; Guenther17 . The vorticity field is computed using the 2D velocity field available at each grid vertex over 1501 time steps. Since the velocity is a 2D field, only the $z$-component of vorticity is non-zero and is used for the analysis. The dataset is available as a single 3D grid where the $z$-direction denotes time.

Figure 9(a) shows three selected time steps of the vorticity field. A persistence-driven simplification with a small threshold of 0.12% successfully removes spurious critical points from each time step, but the Jacobi set remains cluttered as shown in Figure 9(b). A persistence-directed simplification of individual scalar fields may not necessarily produce good results. The simplified Jacobi set Figure 9(c) clearly captures the movement of vortices from left to right. Similar patterns are reported in earlier work, in particular using the lifted Wasserstein matcher (Soler2018LiftedWasserstein, , Figure 14). Figure 9(d) provides a zoomed-in view showing tracks passing through critical regions of maxima and minima, implying that the geometry of the key tracks is well preserved. The simplified Jacobi set is a clean visual representation with reduced zigzag and minimal clutter, thanks to a substantial reduction (400$\times$) in size. The vertical track passing through the center of the domain represents a cluster of saddles that constitute a large critical region near the domain boundary.

Figure 10: (a) Select time steps of the Boussinesq flow dataset showing the $z$-component of vorticity. (b) Tracks of all critical regions. (c-e) Tracks filtered based on their length ($>750,50-750,<50$). (f) A close-up view of the long tracks shows that several critical regions are tracked despite the turbulent behavior.

Boussinesq flow. This dataset representing 2D flow due to a heated cylinder is generated by a simulation using the Gerris flow solver gerrisflowsolver ; Guenther17 . The vorticity field is computed from the velocity field, as described above, over 2001 time steps and is represented as a scalar field on a 3D grid. Figure 10(a) shows the field in select time steps.

The simplified Jacobi set represents tracks of critical points of all types. The tracks are classified based on their length and shown in Figure 10(b-e). A visual inspection indicates a significantly larger number of long tracks identified by our algorithm compared to the ones reported by the lifted Wasserstein matcher (Soler2018LiftedWasserstein, , Figure 12), a method that relies on persistence pair matching for tracking. Despite the innately unsteady nature of this dataset, our algorithm successfully identifies 38 long tracks that last longer than 750 time steps. Figure 10(f) presents a zoomed-in view of the long tracks passing through important extrema regions, emphasizing the validity of the computed tracks. The number of medium length tracks shown in Figure 10(d), indicates that many critical regions have a life span between 50 and 750 time steps. Experimental results on this turbulent dataset demonstrate the utility of the post-processing heuristic. Tracks in proximity are effectively merged, resulting in a 15$\times$ reduction in the total number of tracks, from 4597 to 312. The number of short tracks is small.

Figure 11: Long tracks ($>$300 days) in the sea surface height dataset. Extrema tracks exhibit a westward drift due to wind. The red track depicts unstable saddles near the coastline.

Sea surface height. This dataset from Marine Copernicus copernicus measures the sea surface height from Aug 2021 to Aug 2022. The analysis aims to track significant extrema of the sea surface height over one year to identify features that persist for a long period of time. A length threshold $\epsilon_{l}=300$ days is chosen to identify such long tracks. The tracks of extrema shown in Figure 11 capture the westward drift. This drift is possibly due to the direction of the wind. The simplified Jacobi set is able to capture additional tracks when compared to those extracted using LWM, a direct tracking of critical points (Soler2018LiftedWasserstein, , Figure 11). We use data from a different year due to unavailability of data for the year presented in LWM. Nevertheless, it is evident via a visual comparison that LWM misses multiple prominent features. The noisy red track corresponds to the inconsistent and unstable saddles near the boundary or the coastline. Our algorithm successfully clusters and represents all such saddles with a single track.

9 Conclusions

We described a novel Jacobi set simplification method in the context of time-varying scalar fields, where the Jacobi set represents tracks of critical points. The simplification is based on an analysis of the associated gradient field and an adaptation of the notion of robustness for vector fields. Experiments show that the method produces clutter-free tracks and can potentially identify a larger number of feature tracks when compared with a method for tracking critical points. We also presented a theoretical guarantee of the existence of a simplified vector field to support the algorithm. A generalization of the theoretical analysis to 3D and to regions that are not Jordan is non-trivial, because there is no direct extension of the Riemann mapping theorem.

See pages - of supp.pdf