Local bilinear computation of Jacobi sets

We propose a novel method for the computation of Jacobi sets in 2D domains. The Jacobi set is a topological descriptor based on Morse theory that captures gradient alignments among multiple scalar fields, which is useful for multi-field visualization. Previous Jacobi set computations use piecewise linear approximations on triangulations that result in discretization artifacts like zig-zag patterns. In this paper, we utilize a local bilinear method to obtain a more precise approximation of Jacobi sets by preserving the topology and improving the geometry. Consequently, zig-zag patterns on edges are avoided, resulting in a smoother Jacobi set representation. Our experiments show a better convergence with increasing resolution compared to the piecewise linear method. We utilize this advantage with an efficient local subdivision scheme. Finally, our approach is evaluated qualitatively and quantitatively in comparison with previous methods for different mesh resolutions and across a number of synthetic and real-world examples.


Introduction
There exist numerous methods for the analysis and visualization of multi-fields. One particular approach is the comparison of scalar fields via topological descriptors. The understanding, interpretation, and computation of topological descriptors remain a challenging task [30]. One of these tools is the Jacobi set. The Jacobi set is based on Morse theory and captures gradient alignments among multiple scalar fields. Currently, Jacobi set computation is based on piece- University of Utah (SCI), Salt Lake City, United States 3 Linköping University (ITN), Linköping, Sweden wise linear approximations on triangulations [8]. Piecewise linear interpolation typically produces non-smooth and inaccurate results. In the case of edge-based Jacobi set extraction, it leads to zig-zag patterns and discretization artifacts as illustrated in the top row of Fig. 1. These zig-zag patterns are even preserved for higher resolutions since the traditional method only extracts the edges of the underlying triangulation and does not utilize the cell interior. To address these problems, previous works explore the simplifications of Jacobi sets by either simplifying the scalar functions (indirect methods) or the extracted Jacobi sets (direct methods). While these approaches may improve the Jacobi set representation, they all rely on extracting edges that form the Jacobi set using a piecewise linear approach.
In this paper, our main goal is to extend the piecewise linear computation to a bilinear one to arrive at a more accurate geometry of the extracted Jacobi set. To achieve this, we introduce an equivalent formulation of the piecewise linear method and then extend it to our higher-order approach. Our contributions can be summarized as follows: • We provide an equivalent formulation for the computation of Jacobi sets and embed the piecewise linear method into this framework. • Due to the reformulation, we introduce a novel Jacobi set extraction that depends on bilinear interpolation, Fig. 1 Comparison of the traditional piecewise linear method (top) and our local bilinear method (bottom) for a fixed resolution and different zoom levels on a rectangular grid. The Jacobi set is marked in black and the color-coding corresponds to the ground-truth gradient alignment field (of two analytic functions). The color bar is ordered from positive (top) to negative (bottom), where white indicates the zero. It can be observed that the piecewise linear method results in zig-zag patterns and discretization artifacts, whereas our method leads to a smoother representation of the Jacobi set leading to a smoother, non-edge-based Jacobi set representation. • We present an efficient local subdivision method that further improves the extracted Jacobi set representation.
We compare our method to the piecewise linear method for three different datasets: an analytic dataset, a dataset with a von Kármán vortex street, and the Hurricane Isabel dataset. The paper is structured as follows: After the introduction, we review related work in Sect. 2. To lay the foundations, Sect. 3 provides the theoretical background of Jacobi sets as well as the piecewise linear method to calculate them. In Sect. 4, we introduce the relevant concepts for our local bilinear approach. Subsequently, we derive our approach in three steps and present an efficient local subdivision scheme. Next, Sect. 5 provides a thorough evaluation of our method. Finally, the conclusion with future research directions is presented in Sect. 6.

Related work
The context of this paper is the analysis and visualization of multi-fields via the comparison of scalar fields. In the following, we will summarize previous works that address this topic by extracting joint topological structures [17,30], some of which are closely related to the notion of Jacobi sets Simultaneous analysis of isosurfaces of several scalar fields using a variation density function has been proposed by Nagaraj et al. [23]. Tools studying multivariate data topologically go one step further. Such methods are, in particular, the piecewise linear computation of Jacobi set [8], the Reeb space [11] and its variants [6], Pareto sets [19], and multivariate persistent homology [5]. The Reeb space compresses the components of the level sets of a multivariate mapping and obtains a summary representation of their relationships. Instead of exploring Jacobi sets in the domain, Chattopadhyay et al. [7] studied multi-field topology by considering projections of the Jacobi set into the Reeb space. Thus, these two concepts are shown to be related as the image of the Jacobi sets under the mapping corresponds to certain singularities in the Reeb space [6].
Another related concept is that of Pareto sets [19]. Here, a set of scalar fields is analyzed with respect to consensus or disagreement among ascending and descending manifolds, critical points, and their connectivity. A follow-up paper [18] investigates the relation to Jacobi sets and derives subset and equivalence relations between Jacobi sets and Pareto sets. In addition, the mapper construction [28] and the joint contour nets [6] can be considered as discrete variants of the Reeb space.
All these previous works either introduce a different topological descriptor or establish connections to the concept of Jacobi sets. In contrast, our approach aims to improve the geometrical representation of Jacobi sets with a higher-order interpolation scheme.
As part of topological data analysis tools, the Jacobi set can be used to compare multiple scalar fields [10], track critical points [4,8], define silhouettes of objects [12], and extract ridges in image analysis [26]. It has also been employed in various scientific applications [1,2,8,20,22]. However, in practice, the Jacobi set contains degeneracies and discretization artifacts that can be reduced by simplification of the Jacobi set. Bremer et al. [4] perform persistent simplification of the scalar field to remove small loops of the Jacobi set that lie entirely within successive time steps. Luo et al. [21] compute the Jacobi set by approximating the gradients of a pair of scalar fields. Nagaraj and Natarajan [24] simplify the Jacobi set by explicitly reducing the number of components. Bhatia et al. [3] identify and remove unimportant portions of the Jacobi set by determining how the pair of scalar fields must be modified locally. Our proposed work utilizes bilinear interpolation to obtain a smoother and more accurate geometric representation of the Jacobi set (compared to the piecewise linear approach), while preserving its topology. Therefore, our work can be seen as a completely different approach. To the best of our knowledge, the computation of Jacobi sets under higher-order interpolation schemes has not been explored and understood, and our investigation using bilinear interpolation is the first step in this direction.
Another related concept is the general search for parallel vectors in two fields, e.g., for the computation of ridgelines or the detection of vortex cores in the flow. An introduction to the parallel vectors operator has been given by Peikert and Roth [27], formulated as a search for roots of the crossproduct of the two vector fields. As a computational solution, they suggest an extraction via zero-level isolines for twodimensional vector fields, e.g., by using marching squares. Further improvements for higher-dimensional fields were presented by Theisel et al. [29] and van Gelder et al. [14], who express the problem as a streamline integration problem and make use of tracing algorithms.
While our approach also works with cross-products and an extraction via zero-level isolines, we focus on an edge-based identification of gradient alignments. To be more precise, our Jacobi set extraction starts with scalar fields, computes edge-based gradients in a local bilinear fashion, and identifies the gradient alignments. This means that we do not extract a global gradient field, which is necessary for the parallel vectors operator. Furthermore, due to the preservation of the topology of the piecewise linear method [8], our approach guarantees several topological properties; in particular, it satisfies the Even Degree Lemma [8], which is at the core of our connection method.

Background
In this section, we review the underlying theory of Jacobi sets and outline the piecewise linear computation of Jacobi sets by Edelsbrunner and Harer [8].

Jacobi sets and Morse theory
The Jacobi set is defined for a pair of Morse functions on a smooth d-manifold M. For the purpose of this paper, we restrict our attention to a subset M ⊂ R 2 of the Euclidean space R 2 . Given a smooth function f : M → R, the gradient ∇ f (x) at the point x ∈ M is well-defined by its partial derivatives. The point x is a critical point if the gradient vanishes; otherwise it is a regular point. A critical point x is non-degenerate if the Hessian (i.e., the matrix of secondorder derivatives) is invertible. The function f is a Morse function if every critical point is non-degenerate and the critical points have distinct function values [9, page 128].
Given two Morse functions f , g : M → R, the Jacobi set J( f , g) is defined as the set of points where the gradients of both functions are linearly dependent [8], that is In other words, points in the Jacobi set are referred to as the simultaneous critical points of f and g [8]. For any value t ∈ R, we consider the restriction of f to the level sets can be alternatively formulated as the closure of the set of critical points of f t , that is The Jacobi set is symmetric, i.e., J( f , g) = J(g, f ). Furthermore, it is a smoothly embedded 1-manifold in M [8].

Piecewise linear computation of Jacobi sets
We now describe the traditional piecewise linear computation of Jacobi sets following the work of Edelsbrunner and Harer [8]. Let K be a triangulation of M and let φ, ψ : Vert(K ) → R be functions defined on the vertices of K , Vert(K ). They consider the functions f and g to be the piecewise linear extensions of φ and ψ, respectively, f , g : |K | → R. Then, they model both of these functions as the limits of a series of smooth functions, lim n→∞ f n = f and lim n→∞ g n = g. The Jacobi set is constructed as the limit, J( f , g) := lim n→∞ J( f n , g n ).
The computation of J( f , g) is performed by tracing the critical points of the 1-parameter family of functions h λ = f + λg along the edges of K . To be more precise, the algorithm determines for each edge ab the value λ, such that h λ (a) = h λ (b). This leads to λ = λ ab = ( f b − f a )/(g a −g b ). To determine whether such an edge ab is critical (i.e., belonging to the Jacobi set), the algorithm examines the link of ab. Recall that the star of a simplex τ ∈ K is the collection of simplices that contains τ ; and the link includes simplices in the closure of the star, i.e., St(τ ) = {σ ∈ K | τ ≤ σ } and Lk(τ ) = {ν ∈ Cl(St(τ )) | ν ∩ τ = ∅}. The lower star of a vertex a consists of all simplices in the star that have a as their highest vertex, while the lower link is the portion of the link that bounds the lower star. The lower link of ab is the subcomplex of vertices v i whose h λ ab function values are smaller than the ones from a and b.
Let l be the restriction of h λ ab to the link of ab. The computation relies on the following lemma.
Critical Edge Lemma An edge ab belongs to J if and only if vertex a is a critical point of l. Moreover, the multiplicity of ab in J is the multiplicity of a as critical point of l.
An edge ab with the vertices v 1 and v 2 in its link Lk(ab) (see Fig. 2) is a Jacobi edge or critical if and only if [25] We differentiate between the criticality of the edges: if Eq. 1 holds, the edge is maximal; if Eq. 2 holds, the edge is minimal. These inequalities may lead to degenerate cases, which are resolved by using the Simulation of Simplicity [13]. We refer to Natarajan [25] for details.
Recall that the degree of a vertex is the number of edges that contain the vertex. After computing a set of critical edges, the next step of the algorithm is to construct the 1-manifold J. Based on the Even Degree Lemma by Edelsbrunner and Harer [8], stated below, the algorithm concatenates the edges in pairs, thus forming the 1-manifold J. Our proposed approach builds upon the above edge-based Jacobi set extraction under bilinear interpolation.
Even Degree Lemma The degree of a vertex in J is even.

Our approach
In this section, we first outline the motivation for our approach by establishing a new perspective on the piecewise linear method. This lays the foundation for our proposed local bilinear method that preserves the topology of Jacobi set from [8] and extracts the Jacobi set locally with a bilinear interpolation scheme. Thus, the restriction to an edge-based Jacobi set representation is avoided and the geometry is enhanced. Finally, an appropriate local subdivision method for an efficient Jacobi set extraction is presented.

Motivation and derivation
The piecewise linear computation of Jacobi sets by Edelsbrunner and Harer [8] has several limitations. One of these is that only the edges of the underlying triangulation can be identified as components of the Jacobi set. This leads to zig-zag patterns, poor convergence for increasing mesh resolution, and an imprecise Jacobi set representation. To avoid these issues, the natural way is to use a higher-order interpolation scheme, which is the core of our approach.
To apply this, we propose to use an equivalent formulation of Jacobi sets, that is In our setting, where M is a smooth 2-manifold embedded in R 3 , the cross-product of ∇ f (x) and ∇g(x) can be expressed with the help of the unit vectors e i and the determinant Assuming that M lies in the x y-plane, the cross-product is We define κ x , called gradient alignment value, as The gradient alignment value κ x characterizes the linear independence of the gradients at the point x ∈ M. This leads to the scalar field κ x : M → R as a function of x, which we call the gradient alignment field. According to the Jacobi set definition, the point x belongs to the Jacobi set J( f , g) if and only if κ x ( f , g) = 0. Our next step is to examine how the (smooth) gradient alignment field κ x ( f , g) behaves for both linearly and bilinearly interpolated functions f and g, the corresponding gradient alignment fields are denoted κ li x and κ bi x , respectively. We assume f to be linear in each of the triangles {abv 1 } and {av 2 b} (for the notation, see Fig. 3). For the triangle {abv 1 }, we use the formula f abv 1 (x, y) := a f 00 +a f 10 x +a f 01 y, which leads to the following system of linear equations: where The analytical solution leads to the following gradient of f abv 1 : describes the area of the parallelogram spanned by the edges v 1 a and v 1 b. An analogous calculation can be done for f av 2 b , g abv 1 , and g av 2 b . According to Equation 4, the gradient alignment value at the vertex v 1 is given by Analogously, for the triangle {av 2 b}, we compute the gradient alignment value at the vertex v 2 via These calculations show that the resulting values κ li v 1 ( f , g) and κ li v 2 ( f , g) stem from the gradient alignment field κ li x , that is constant on each triangle (because the gradients ∇ f and ∇g are constant on each triangle). With these observations and the fact that the Jacobi set is equal to the zero-level set of κ x , it is possible to identify critical edges via a change of sign with respect to κ li v 1 ( f , g) and κ li v 2 ( f , g). Comparing these considerations with the approach of Edelsbrunner and Harer [8], we will notice that the two approaches are equivalent. To show this, we consider an edge ab with Lk(ab) = {v 1 , v 2 }. For the edge ab to be critical, either Eq. 1 or Eq. 2 needs to be true. Let us assume that Eq. 1 holds. We can reformulate it in terms of κ li v 1 and κ li v 2 : An analogous calculation can be done for Eq. 2 with swapped signs for both conditions. As a result, the following relations, which explicitly show the connection between the approach of Edelsbrunner and Harer and the gradient alignment field, can be established: where the term g a − g b determines whether an edge ab is maximal or minimal. At this point, the consequences are twofold. Firstly, the gradient alignment values κ li v 1 ( f , g) and κ li v 2 ( f , g) characterize a new Jacobi set condition for the edge ab via Secondly, we observe that the piecewise linear approach by Edelsbrunner and Harer is also equivalent to the condition in Eq. 8. In fact, the conditions only differ in a scaling factor: whereas the gradient alignment field for linearly interpolated functions f and g (on each triangle) encodes both the topology and geometry in the condition (due to the scaling A abv 1 and A av 2 b ), the approach by Edelsbrunner and Harer solely relies on topological relations. Now, the next step is to investigate how the gradient alignment field κ x ( f , g) behaves for bilinearly interpolated functions f and g. To this end, let f and g be given on the The coefficients can be analogously computed to Eq. 5 via a system of linear equations. However, in this case the gradients are not constant and are given by Hence, the gradient alignment field κ x is linear as the following general calculation shows: This aspect enables the computation of the zero-level set in the quadrilateral, instead of identifying critical edges of the underlying triangulation. In fact, the method by Edelsbrunner and Harer can be extended using a bilinear interpolation scheme.
To be more precise, the gradient alignment values κ li v 1 and κ li v 2 (see Eq. 6 and Eq. 7) are used for the topological identification of critical edges, whereas for the geometrical identification a bilinear interpolation scheme is used.
To this end, the gradient alignment values κ bi v 1 and κ bi v 2 in the case of bilinearly interpolated f and g are given by the following formulas (we skip the technical details): The values κ li a and κ li b are given by with the area for counter-clockwise sorted p 1 ,p 2 , and p 3 : The denominator D is computed via

Algorithm
An overview of our local bilinear method is given in Fig. 4, and the algorithm can be divided into three steps:

Computation of linear gradient alignment values κ li
x with the identification of critical edges. 2. Computation of Jacobi set points via the bilinear gradient alignment values κ bi x . 3. Generation of the Jacobi set representation J( f , g).

Identification of critical edges
In the previous subsection, we formulated a method (compare Eq. 8), that is topologically equivalent to the one of Edelsbrunner and Harer [8] in terms of the (linear) gradient alignment values κ li v 1 and κ li v 2 . We use these values to identify critical Jacobi set edges. This is the first part of our method and also the initial step in Algorithm 1. To be more precise, given edges E = (e i ), we iterate over each edge e i (line 2), identify v 1 and v 2 as the lower link of e i (line 3), and compute the gradient alignment values κ li v 1 and κ li v 2 via Eq. 6 and Eq. 7 (line 4). Afterward, the critical edge test via the change of sign method is performed (line 5).
While the piecewise linear approach stops here (identifying critical edges of the underlying triangulation), the next step of our method is to determine a geometrically more precise Jacobi set representation via Jacobi set points.

Computation of Jacobi set points
In contrast to the piecewise linear computation by Edelsbrunner and Harer, we avoid marking critical edges of the underlying triangulation by introducing Jacobi set points. As illustrated in Fig. 4a, a Jacobi set point is conceptually assigned to a critical edge ab and lies on the line segment that connects the vertices v 1 and v 2 via m (compare Fig. 2). Find v 1 , v 2 ∈ Lk(e i ) 4 Compute κ li v 1 , κ li v 2 // Eq. 6 and Eq. 7 Compute κ bi v 1 , κ bi v 2 // Eq. 10 and Eq. 11 add (e i , p i ) to J E Given the fact that bilinearly interpolated functions f and g result in a linear gradient alignment field on the quadrilateral {av 1 bv 2 }-cell, there is only one zero on this line segment if the signs of the values κ v 1 and κ v 2 differ. In this context, the zero of the line segment is given by where λ is the weighting coefficient that characterizes the location of the zero via Thus, the main idea is that the collection of all Jacobi set points approximates the Jacobi set (see Fig. 4b and c). The computation of Jacobi set points is the second part of Algorithm 1 (starting with a critical edge e i in line 6) and will be explained in the following in more detail.
To compute a Jacobi set point, we compute the bilinear gradient alignment values κ bi v 1 and κ bi v 2 (line 6). Then, we perform the change of sign method for κ bi v 1 and κ bi v 2 (line 7). In the case of different signs, the weighting coefficient λ characterizing the Jacobi set point is computed via bilinear gradient alignment values κ bi v 1 and κ bi v 2 (line 8). Otherwise, linear gradient alignment values κ li v 1 and κ li v 2 are used (line 10). Finally the exact location of the Jacobi set point is computed via Eq. 12 (line [11][12][13][14][15] and added to the list of pairs (e i , p i ) of Jacobi set points p i and critical edges e i (line 16).

Connection of Jacobi set points
Given a collection of Jacobi set points computed by Algorithm 1, our goal is to connect them in a way such that they geometrically approximate the Jacobi set J( f , g). To achieve this, we utilize the preserved topological properties of the (linear) gradient alignment values (from Edelsbrunner and Harer, see Eq. 8) to apply the Even Degree Lemma. According to the lemma, the vertices a and b of a Jacobi edge ab have an even degree. It implies that a Jacobi edge ab is connected to an even number of other Jacobi edges. This is due to the fact that the degrees of a and b are even and since both degrees include the edge ab as well, the number of Jacobi edges connected to ab equals (degr(a) − 1) + (degr(b) − 1). The important consequence of this is: Each Jacobi set point is connected to an even number (greater than zero) of Jacobi set points and, in particular, there are no isolated Jacobi set points.
Based on this result, we can describe a method to connect Jacobi set points. An illustration of this method is presented in Fig. 5, and a detailed description is given by Algorithm 2, which will be explained in the following. Given the list of pairs (e i , p i ) of Jacobi set points p i that are related to critical edges e i , we iterate over each critical edge e i (line 2) to find all critical edges e i j that share a vertex with e i (line 3). These are used to construct line segments l = (p i , p j ) with the respective Jacobi set points p i and p j (line [4][5]. This method produces a collection of line segments that geometrically approximates the Jacobi set while preserving the topology. Algorithm 2 could be considered as a "dualization" process where each critical edge is replaced by its corresponding Jacobi set point (as its "dual"); in fact, two Jacobi set points are connected by a line segment if their corresponding critical edges share the same vertex. This aspect is visualized in Fig. 5, where we use a color-coding and numbering to highlight the connectivity. In the top, the pairs of critical edges and Jacobi set points are represented (the dotted lines emphasize the pairs). The bottom image shows the Jacobi set representation with our connection method. The effect of our connection method is that the adjacency of critical edges (e.g., edge 1 is connected to the edges 0, 2, 3 and 5) is reflected by the connectivity of Jacobi set points (see the respective Jacobi sets points). In this way, the topology of the piecewise linear method is preserved. We provide a formal proof of the homotopy equivalence in Appendix.

Local subdivision method
One of the advantages of our local bilinear approach is the usage of a geometrically more precise extraction of the Jacobi set. This typically leads to better convergence properties compared to piecewise linear schemes. Therefore, we propose a local subdivision scheme that makes use of the already identified critical edges by only doing computations in these locations. Given a critical edge e i = ab and the related Jacobi set point p i within a {av 1 bv 2 }-cell (see Fig. 6), we now describe a single subdivision step for the edge ab:  This procedure is done for each critical edge, which means a local efficient subdivision scheme is used. An illustration is presented in Fig. 4d. It can be observed that T-nodes occur, which violate the rules of a triangulation. However, for our purposes, we do not need a consistent global triangulation for further computations. The new critical edges and Jacobi set points will only be computed in the local neighborhood, and therefore, the T-nodes do not affect the repeated application of the local subdivision method.

Results
To evaluate our local bilinear method qualitatively and quantitatively, we apply it to three different datasets: an analytical example, a simulated von Kármán vortex street, and the multi-field Hurricane Isabel dataset. Additionally, it is compared to the piecewise linear method by Edelsbrunner and Harer [8]. The methods are performed using MATLAB (R2022a) on a MacBook Pro with an Intel Dual-Core i5 CPU @ 3.1 GHz and 8 GB of RAM. The computation times for each of the datasets are given in Table 1.
For the quantitative analysis, we introduce the gradient alignment measure μ that measures the accuracy of the In short, the measure computes the deviation of the Jacobi set with regard to the globally computed gradient alignment field κ x ( f , g). The background of this measure is that the gradient alignment value is zero for x ∈ J( f , g).

Analytic dataset
The first dataset consists of two artificially generated scalar fields. Both are constructed by the sum of three bivariate normal distributions on the unit square. To get an impression of the two functions, a visualization of the contour lines is presented in Fig. 7. The figure also shows the Jacobi set of these two functions (black line) given by an analytical formula. Therefore, we can compare the Jacobi set extracted by our local bilinear approach (or by the piecewise linear approach) with the true Jacobi set.
In Fig. 1, the Jacobi sets are computed by the piecewise linear approach (top row) and our local bilinear method (bottom row) for a fixed resolution and different zoom levels on a (triangulated) rectangular grid. First of all, both methods identify the fundamental topological structure correctly. For each zoom level, the piecewise linear approach produces a solution that highly suffers from the edge-based geometric representation of the Jacobi set (which stems from the underlying triangulation), leading to zig-zag patterns. Furthermore, the zoomed-in image (top right) reveals that the Jacobi set representation consists of several ("glued") triangles of the underlying triangulation, resulting in misleading topological patterns. While our local bilinear method obtains the same topology, these triangles are merged geometrically, leading to a clearer geometric representation. In sum, our local bilinear method results in a smoother and more precise Jacobi set representation (e.g., see third column).
In the following, we also evaluate the proposed local subdivision method. It can be applied to both our local bilinear method and the piecewise linear method. The qualitative results are presented in Fig. 8, and the corresponding gradient alignment measure is depicted in Fig. 9.
In general, both Jacobi set representations produce incorrectly detected regions and artifacts for the initial configuration, especially at the bottom right. This is due to the coarse irregular mesh and resulting degenerate areas. With the application of subdivision, both methods produce better results; however, our local bilinear method results in much smoother representations. This observation can also be noticed in the error plot in Fig. 9. The gradient alignment measure shows that our approach results in less error and has a faster convergence than the piecewise linear method.

Von Kármán vortex street dataset
While the analytic dataset provides first insights into our approach, we now apply our method to a simulated von Kármán vortex street [15]. The dataset consists of 1501 time steps describing the velocity and has a resolution of 640 × 80. We use the magnitudes of the flow velocity for two consecutive time steps as the two scalar fields for the Jacobi set computation. For our examples, we start at time step 1500, where the vortex street is fully emerged.
A visualization of the von Kármán vortex street for this time step is depicted in Fig. 10 top. The figure also shows the Jacobi set computed by the piecewise linear approach (middle) and our approach (bottom). In general, we observe that the Jacobi set captures the topological structure of each Kármán vortex. However, similar to previous observations, the zig-zag patterns of the piecewise linear approach make it hard to distinguish between different topological structures.
For a zoomed-in view, the top rows of Fig. 11 show the framed area in Fig 10 in more detail. We observe that both  The gradient alignment measure for different subdivision steps (analogously to Fig. 8) for the piecewise linear (red) and our local bilinear (blue) method. The vertical dashed lines correspond to the snapshots in Fig. 8 methods fail to produce the correct topological structure at the top left (compare the gradient alignment field in the background). However, due to our connection method, our Jacobi set representation suggests that this region is problematic and may contain wrongly detected areas. In the piecewise linear case, the regular identification of Jacobi set patterns via the underlying triangulation makes it hard to distinguish between correctly and incorrectly detected areas. Besides of that, at the middle and top right areas, our local bilinear method captures the geometry better due to the smooth Jacobi set representation.
Another comparison of these two methods is done by using the local subdivision method (rows 2-4 in Fig. 11). We observe for both methods that the incorrect topological structures in the top left disappear after one subdivision step. Another important aspect is shown in the third row (with two subdivision steps). Topologically identified patterns, such as the isolated triangle at the bottom left, can disappear visually for our bilinear method due to geometrical contraction. This can never happen for the piecewise linear approach. Applying subdivision four times, we observe that the representation of topological structures is significantly improved. While the piecewise linear approach identifies the circular shape rather poorly (due to the underlying grid), the geometry of our method captures the circularity more clearly (in the limit, the shape is nearly an ellipse).

Hurricane Isabel dataset
Jacobi sets are motivated as a topological descriptor of multiple scalar fields with different quantities. Here, we perform  Fig. 10 (see the white box). As scalar functions, flow magnitudes of two con-secutive time steps are used. Furthermore, the subdivision method is applied in the zoomed-in areas for both methods Fig. 12 Comparison of the piecewise linear (a) and local bilinear (b) computation of Jacobi sets for the Hurricane Isabel dataset. For the Jacobi set, a pressure and temperature scalar field is used. Furthermore, our local bilinear method with two subdivision steps is illustrated (c) the Jacobi set extraction on the multi-field real-world dataset Hurricane Isabel 1 . The dataset consists of multiple 3D scalar fields and a 3D velocity vector field that have a resolution of 500 × 500 × 100 and 48 time steps.
We use the pressure and temperature scalar field at time step 30 and height 50 to compare the Jacobi set computation methods in Fig. 12. At this height, the hurricane has a large spatial expansion and reveals many characteristic structures.
Comparing the piecewise linear method (Fig. 12a) with our method (Fig.12b), we observe that many areas are marked 1 Hurricane Isabel data produced by the Weather Research and Forecast (WRF) model, courtesy of NCAR and the U.S. National Science Foundation (NSF) in the top left region. Whereas in the piecewise linear case it is nearly impossible to track topological structures, our method (without subdivision) already has a clearer representation due to geometrical enhancements. However, it can be observed that the preservation of topology leads to a folding of line segments, which can clutter areas in the representation (see exemplary Fig. 5). Our subdivision method (Fig 12c) solves this issue by thinning out those areas. This leads to a much clearer representation.
To further demonstrate that our proposed local subdivision method improves the results, the other zoomed-in area at the center shows that the identified topological structures match more clearly the zero-level set (white areas) of the gradient alignment field in the background. In fact, several connected topological patterns are separated correctly.

Conclusion
In this paper, we recapped the piecewise linear computation of Jacobi sets and embedded its theory into our new formulation. This view enables a new formulation that extends the piecewise linear computation to be based on a bilinear interpolation scheme. While preserving the topological structures, our bilinear interpolation results in a smoother and geometrically more accurate Jacobi set representation. In addition, our efficient local subdivision method tailored to our approach further improves the Jacobi set representation.
In the future, we could extend our local bilinear method to higher dimensions as well as more than two scalar functions. Additionally, other interpolation methods could be used that are based on finite elements or higher interpolation schemes such as biquadratic or bicubic interpolation.
Another idea is to evaluate our local bilinear method in terms of detecting new features. This not only includes the raw results but also the representation via connected line segments. In this regard, the appearance and the spatial spread of the connected line segments (representing a band-like structure) could contain more information about the underlying data, e.g., uncertainty or other features. In addition, since our method significantly facilitates the investigation of Jacobi sets, another research direction is to find new interpretations and applications of Jacobi sets across diverse research areas, such as image processing or computer vision.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Fig. 13
Applying Nerve Lemma to a piece of the Jacobi set

Appendix: Proof of homotopy equivalence
In Algorithm 2, each critical edge is replaced by its corresponding Jacobi set point, where two Jacobi set points are connected by a line segment if their corresponding critical edges share the same vertex. We claim that such a mapping preserves the topology due to the Nerve Lemma (see [16, Corollary 4G.3, page 459] and [9, Nerve Theorem, page 71]). In our context, as shown in Fig. 13 left, the input to Algorithm 2 is a set of critical edges e i , which form a simplicial complex K together with their vertices. Recall the closure of a simplex contains the set of its faces. We form a cover U of K by treating the closure of each critical edge e i as a cover element U i := Cl(e i ), where K ⊆ U i ∈U U i (see Fig. 13 middle). Each cover element U i is a convex and closed set, and the intersection of U i and U j is contractible. Therefore, the Nerve Lemma applies and the nerve of U, which contains a 3-simplex (a yellow tetrahedron) as shown in Fig. 13 right, is homotopy equivalent to K . The output of Algorithm 2 is by construction precisely the 1-skeleton of the nerve of U.
Ingrid Hotz Ingrid Hotz is currently a Professor in Scientific Visualization at the Linköping University in Sweden. Her research interests lie in data analysis and scientific visualization, ranging from basic research questions to effective solutions to visualization problems in applications.
Bei Wang Bei Wang is an Assistant Professor at the School of Computing and a faculty member at the SCI Institute, University of Utah. She received her PhD in Computer Science from Duke University. Her research interests include topological data analysis, data visualization, computational topology, computational geometry, machine learning, and data mining.