Maxime Soler
email: soler.maxime@total.com

S A Total

M Élanie Plainchault
email: melanie.plainchault@total.com

S A Bruno Total

Conche

S A Julien Total

Tierny

Lifted Wasserstein Matcher for Fast and Robust Topology Tracking

Keywords: Topological Data Analysis, Optimal Transport, Feature Tracking

This paper presents a robust and efficient method for tracking topological features in time-varying scalar data. Structures are tracked based on the optimal matching between persistence diagrams with respect to the Wasserstein metric. This fundamentally relies on solving the assignment problem, a special case of optimal transport, for all consecutive timesteps. Our approach relies on two main contributions. First, we revisit the seminal assignment algorithm by Kuhn and Munkres which we specifically adapt to the problem of matching persistence diagrams in an efficient way. Second, we propose an extension of the Wasserstein metric that significantly improves the geometrical stability of the matching of domain-embedded persistence pairs. We show that this geometrical lifting has the additional positive side-effect of improving the assignment matrix sparsity and therefore computing time. The global framework computes persistence diagrams and finds optimal matchings in parallel for every consecutive timestep. Critical trajectories are constructed by associating successively matched persistence pairs over time. Merging and splitting events are detected with a geometrical threshold in a post-processing stage. Extensive experiments on real-life datasets show that our matching approach is up to two orders of magnitude faster than the seminal Munkres algorithm. Moreover, compared to a modern approximation method, our approach provides competitive runtimes while guaranteeing exact results. We demonstrate the utility of our global framework by extracting critical point trajectories from various time-varying datasets and compare it to the existing methods based on associated overlaps of volumes. Robustness to noise and temporal resolution downsampling is empirically demonstrated.

INTRODUCTION

Performing feature extraction and object tracking is an important topic in scientific visualization, for it is key to understanding timevarying data. Specifically, it allows to detect and track the evolution of regions of interest over time, which is central to many scientific domains, such as combustion [START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topology-based data segmentation[END_REF], aerodynamics [START_REF] Hannon | Aircraft wake vortex detection and measurement with pulsed solid-state coherent laser radar[END_REF], oceanography [START_REF] Ringler | A multi-resolution approach to global ocean modeling[END_REF] or meteorology [START_REF] Zhang | Exploring scalar fields using multiple sensor platforms: Tracking level curves[END_REF]. With the increasing power of computational resources and resolution of acquiring devices, efficient methods are needed to enable the analysis of large datasets.

The emergence of new paradigms for scientific simulation, such as in-situ and in-transit [START_REF] Bennett | Combining in-situ and in-transit processing to enable extreme-scale scientific analysis[END_REF][START_REF] Moreland | Examples of in transit visualization[END_REF][START_REF] Rasquin | Electronic poster: Co-visualization of full data and in situ data extracts from unstructured grid cfd at 160k cores[END_REF][START_REF] Rivi | In-situ visualization: State-of-the-art and some use cases[END_REF][START_REF] Yu | In situ visualization for large-scale combustion simulations[END_REF], clearly exhibits the ambition to reach toward exascale computing [START_REF] Committee | Synergistic challenges in data-intensive science and exascale computing[END_REF] in the forthcoming years. In this context, as both spatial and temporal resolutions of acquired or simulated datasets keep on increasing, understanding the evolution of features of interest throughout time proves challenging.

Topological data analysis has been used in the last decades as a robust and reliable setting for hierarchically defining features in scalar data [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF]. In particular, its successful application to time-varying data [START_REF] Bremer | Analyzing and tracking burning structures in lean premixed hydrogen flames[END_REF][START_REF] Sohn | Time-varying contour topology[END_REF] makes it a prime candidate for tracking. Both topological analysis and feature tracking have been applied in-situ [START_REF] Landge | In-situ feature extraction of large scale combustion simulations using segmented merge trees[END_REF][START_REF] Zhang | Insitu feature-based objects tracking for large-scale scientific simulations[END_REF], which demonstrates their interest in the context of large-scale data. Nonetheless, major bottlenecks of state-of-the art topology tracking methods are still the high required computation cost as well as the need for high temporal resolution.

In this paper, we propose a novel feature-tracking framework, which correlates topological features in time-varying data in an efficient and meaningful way. It is the first approach, to the best of our knowledge, combining the setting of topological data analysis with optimal transport for the problem of feature tracking. More precisely, the key idea is to use combinatorial optimization for matching topological structures (namely, persistence diagrams) according to a fine-tuned metric. After exposing our formal setting (Sec. 2), we introduce an extension of the exact assignment algorithm by Kuhn and Munkres [START_REF] Kuhn | The hungarian method for the assignment problem[END_REF][START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF] that we adapt in an efficient way to the case of persistence diagrams (Sec. 3). We highlight the issues raised by the classical Wasserstein metric between diagrams, and propose a robust lifted metric that overcomes these limitations (Sec. 4). We then present the detailed tracking framework (Sec. 5). Extensive experiments demonstrate the utility of our approach (Sec. 6).

Related work

Our framework encompasses the definition, correlation and tracking of topological features in scalar fields. As such, it is related to topological data analysis of scalar fields, tracking techniques, the definition of metrics and combinatorial optimization. Topological analysis techniques [START_REF] Floriani | Morse complexes for shape segmentation and homological analysis: discrete models and algorithms[END_REF][START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF][START_REF] Heine | A survey of topologybased methods in visualization[END_REF][START_REF] Pascucci | Topological Data Analysis and Visualization: Theory, Algorithms and Applications[END_REF] have demonstrated their ability over recent years to capture features of interest in scalar data in a generic, robust [START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF][START_REF] Edelsbrunner | Topological persistence and simplification[END_REF] and efficient manner, for many applications as turbulent combustion [START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topology-based data segmentation[END_REF], computational fluid dynamics [START_REF] Favelier | Visualizing ensembles of viscous fingers[END_REF], material sciences [START_REF] Gyulassy | Interstitial and interlayer ion diffusion geometry extraction in graphitic nanosphere battery materials[END_REF], chemistry [START_REF] Guenther | Characterizing molecular interactions in chemical systems[END_REF], astrophysics [START_REF] Sousbie | The persistent cosmic web and its filamentary structure: Theory and implementations[END_REF], medical imaging [START_REF] Bock | Topoangler: Interactive topology-based extraction of fishes[END_REF][START_REF] Carr | Simplifying flexible isosurfaces using local geometric measures[END_REF]. One reason for their success in applications is the possibility for domain experts to easily translate high level structural notions into topological abstractions, such as contour trees [START_REF] Carr | Computing contour trees in all dimensions[END_REF], Reeb Graphs [START_REF] Biasotti | Reeb graphs for shape analysis and applications[END_REF][START_REF] Pascucci | Robust online computation of Reeb graphs: simplicity and speed[END_REF], Morse-Smale complexes [START_REF] Gyulassy | A practical approach to morse-smale complex computation: Scalability and generality[END_REF]. For instance, in astrophysics the cosmic web can be extracted by querying the most persistent 1-separatrices of the Morse-Smale complex connected to maxima of matter density [START_REF] Sousbie | The persistent cosmic web and its filamentary structure: Theory and implementations[END_REF]. Similar domain-specific notions are translated into topological terms in the above examples. Feature tracking: Topology has been used for feature extraction and tracking in the context of vector fields [START_REF] Post | The state of the art in flow visualisation: Feature extraction and tracking[END_REF][START_REF] Reininghaus | Efficient computation of combinatorial feature flow fields[END_REF][START_REF] Tricoche | Topology tracking for the visualization of time-dependent two-dimensional flows[END_REF], mostly relying on stream lines, path lines [START_REF] Shi | Path Line Attributes -an Information Visualization Approach to Analyzing the Dynamic Behavior of 3D Time-Dependent Flow Fields[END_REF][START_REF] Shi | Path line oriented topology for periodic 2D time-dependent vector fields[END_REF][START_REF] Theisel | Feature flow fields[END_REF][START_REF] Theisel | Stream line and path line oriented topology for 2d time-dependent vector fields[END_REF][START_REF] Theisel | Topological methods for 2d time-dependent vector fields based on stream lines and path lines[END_REF], or tracking punctual singularities [START_REF] Klein | Scale-space tracking of critical points in 3d vector fields[END_REF]. For the latter, a forward streamline integration of critical points is performed in a specific scale space, which adapted for time-varying data would require knowledge about the evolution of the field, and for instance to compute time-derivatives.

For scalar data, features are defined based on attributes that are either geometric (isosurfaces, thresholded regions), or topological (contour trees, Reeb graphs). Similarly, tracking approaches either rely on geometrical (volume overlaps, distance between centers of gravity) or topological extracts (Jacobi set, segment overlap).

Geometrical approaches are based on thresholded connected components [START_REF] Silver | Visualizing evolving scalar phenomena[END_REF], glyphs [START_REF] Reinders | Visualization of timedependent data with feature tracking and event detection[END_REF], cluster tracking [START_REF] Grottel | Visual verification and analysis of cluster detection for molecular dynamics[END_REF], petri nets [START_REF] Ozer | Activity detection in scientific visualization[END_REF], or propose a hierarchical representation [START_REF] Gu | Transgraph: Hierarchical exploration of transition relationships in time-varying volumetric data[END_REF]. Similarly, the core methodology for associating topological features for tracking is often based on overlaps of geometrical domains along with other attributes [START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topology-based data segmentation[END_REF][START_REF] Bremer | Analyzing and tracking burning structures in lean premixed hydrogen flames[END_REF][START_REF] Saikia | Global feature tracking and similarity estimation in time-dependent scalar fields[END_REF][START_REF] Silver | Object-oriented visualization[END_REF][START_REF] Silver | Volume tracking[END_REF][START_REF] Silver | Tracking and visualizing turbulent 3d features[END_REF][START_REF] Silver | Tracking scalar features in unstructured data sets[END_REF][START_REF] Sohn | Time-varying contour topology[END_REF], on tracking Jacobi sets [START_REF] Edelsbrunner | Jacobi sets of multiple morse functions[END_REF][START_REF] Edelsbrunner | Timevarying reeb graphs for continuous space-time data[END_REF], or matching isosurfaces in higher-dimensional spaces [START_REF] Ji | Efficient Isosurface Tracking Using Precomputed Correspondence Table[END_REF][START_REF] Ji | Volume tracking using higher dimensional isosurfacing[END_REF].

Such approaches usually test features in two consecutive timesteps against one another for potential overlaps, then draw the best correspondence between features according to some criterion. Typical criteria include optimizing the overlapping volume, mass, distance between centroids, or a combination of these [START_REF] Reinders | Visualization of timedependent data with feature tracking and event detection[END_REF][START_REF] Samtaney | Visualizing features and tracking their evolution[END_REF]. For this to work, the temporal sampling rate of the underlying data must be such that features in two consecutive timesteps effectively overlap. This first criterion is thus not very robust to temporal downsampling.

Other approaches rely on global optimization [START_REF] Ji | Feature tracking using earth mover's distance and global optimization[END_REF], using the Earth Mover's distance [START_REF] Levina | The earthmover's distance is the mallows distance: some insights from statistics[END_REF] between geometrical features with various attributes such as centroid position, volume and mass. This does not, however, benefit from the natural definition of features offered by topological data analysis, nor from the possibility to simplify features in a hierarchical way. This is a real drawback in the context of noisy data as it implies dealing with large, computation-intensive optimization problems between every pair of timesteps.

Once features have been defined, and a methodology established to associate them in consecutive timesteps, the tracking representation is quite independent of whether geometrical or topological arguments have been used. Most often, graphs are used [START_REF] Laney | Understanding the structure of the turbulent mixing layer in hydrodynamic instabilities[END_REF][START_REF] Robbins | Visualization of splitting and merging processes[END_REF][START_REF] Widanagamaachchi | Tracking features in embedded surfaces: Understanding extinction in turbulent combustion[END_REF], such as Reeb graphs [START_REF] Edelsbrunner | Timevarying reeb graphs for continuous space-time data[END_REF][START_REF] Weber | Feature Tracking Using Reeb Graphs[END_REF] and nested tracking graphs [START_REF] Lukasczyk | Nested tracking graphs[END_REF]. Many popular graph structures are accounted for in [START_REF] Wang | Graphs in scientific visualization: A survey[END_REF]. An inconvenience of extracting rich tracking structures such as these without taking careful attention to potential noise is that it makes the interpretation quite difficult. In [START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topology-based data segmentation[END_REF], the tracking graph is dense and intricate, making exploration impractical. It is therefore mandatory to do filtering and simplification in a post-processing stage, or to cleverly discard noisy events beforehand. Assignment problems: Since we revisit the original algorithm by Kuhn and Munkres, we discuss here some related work in combinatorial optimization. The assignment problem is the discrete optimization problem consisting of finding a perfect matching of optimal cost in a weighted bipartite graph [START_REF] Bertsekas | Network optimization: Continuous and discrete models[END_REF][START_REF] Burkard | Assignment Problems[END_REF][START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF]. In other terms, the problem is to find an optimal one-to-one correspondence between discrete entities (such as singularities in a scalar field at two different timesteps), with a cost associated to each possible correspondence. It can be solved with the seminal Kuhn-Munkres algorithm [START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF]. The auction algorithm [START_REF] Bertsekas | Network optimization: Continuous and discrete models[END_REF][START_REF] Bertsekas | The auction algorithm for transportation problem[END_REF][START_REF] Kerber | Geometry helps to compare persistence diagrams[END_REF] is another popular approach for solving the assignment problem with a user-defined error threshold on the resulting assignment cost. In practice, this threshold is often set to 1% of the scalar range. A more general, continuous formulation of this problem is at the heart of Transportation theory [START_REF] Kantorovich | On the translocation of masses[END_REF][START_REF] Monge | Mémoire sur la théorie des déblais et des remblais[END_REF][START_REF] Villani | Optimal Transport: Old and New[END_REF]. Modern techniques [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] have attracted acute interest for making this theoretical setup central to shape correlation [START_REF] Solomon | Entropic metric alignment for correspondence problems[END_REF] and interpolation [START_REF] Solomon | Convolutional wasserstein distances: efficient optimal transportation on geometric domains[END_REF], which do bear resemblance to feature tracking. Metrics: Since we introduce a new metric for the matching of persistence diagrams, we discuss in the following existing metrics traditionally used in topological data analysis. The Bottleneck and the interleaving distances [START_REF] Chazal | Proximity of Persistence Modules and their Diagrams[END_REF][START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF] have been widely investigated to study the stability of persistence diagrams. These metrics have been notably adapted in the context of kernel methods [START_REF] Carrière | Sliced wasserstein kernel for persistence diagrams[END_REF][START_REF] Reininghaus | A stable multi-scale kernel for topological machine learning[END_REF] in machine learning. In particular, the Bottleneck distance is a special case of the more general Wasserstein metric [START_REF] Kantorovich | On the translocation of masses[END_REF][START_REF] Monge | Mémoire sur la théorie des déblais et des remblais[END_REF] applied to diagram points, also known as the Earth Mover's Distance [START_REF] Levina | The earthmover's distance is the mallows distance: some insights from statistics[END_REF]. The standard approach for computing the discrete Wasserstein metric relies on solving the associated assignment problem, either with an exact Kuhn-Munkres approach [START_REF] Morozov | Dionysus[END_REF][START_REF] Weaver | An implementation of the kuhn-munkres algorithm[END_REF] or with an auction-based approximation [START_REF] Kerber | Geometry helps to compare persistence diagrams[END_REF]. However, as discussed in Sec. 4, when these methods (metric-based [START_REF] Chazal | Proximity of Persistence Modules and their Diagrams[END_REF][START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF] or kernel-based [START_REF] Carrière | Sliced wasserstein kernel for persistence diagrams[END_REF][START_REF] Reininghaus | A stable multi-scale kernel for topological machine learning[END_REF]) are applied as-is for tracking purposes, a high geometrical instability occurs which impairs the tracking robustness, as already observed in the case of vineyards [START_REF] Cohen-Steiner | Vines and vineyards by updating persistence in linear time[END_REF]. Our work (see Sec. 4) addresses this issue.

Contributions

This paper makes the following new contributions:

1. Approach: We present a sound and original framework, which is the first combining topology and transportation for feature tracking, comparing favorably to other state-of-the-art approaches, both in terms of speed and robustness.

2. Metric: We extend traditional topological metrics, for the needs of time-varying feature tracking, notably enhancing geometrical stability and computing time.

Algorithm:

We extend the assignment method by Kuhn and Munkres to solve the problem of persistence matchings in a fast and exact way, taking advantage of our metric.

PRELIMINARIES

This section describes our formal setting and presents an overview of our approach. It contains definitions that we adapted from Tierny et al. [START_REF] Tierny | The Topology ToolKit[END_REF]. An introduction to topology can be found in [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF].

Background

Input data: Without loss of generality, we assume that the input data is a piecewise linear (PL) scalar field f : M → R defined on a PL d-manifold M with d ≤ 3. Values are given at the vertices of M and linearly interpolated on higher dimension simplices.

Critical points: Given an isovalue i ∈ R, the sub-level set of i,

noted f -1 -∞ (i), is the pre-image of the open interval (-∞, i) onto M through f : f -1 -∞ (i) = {p ∈ M | f (p) < i}. The sur-level set is symmetrically given by f -1 +∞ (i) = {p ∈ M | f (p) > i}.
These two objects serve as segmentation tools in many analysis tasks [START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topology-based data segmentation[END_REF].

The points p ∈ M where the topology of f -1 -∞ ( f (p)ε) differs from that of f -1 -∞ ( f (p) + ε) are the critical points of f and their values are called critical values. Critical points can be classified with their index I , which is 0 for minima, 1 for 1-saddles, d -1 for (d -1)-saddles and d for maxima, with d the dimension of M . Persistence diagrams: The distribution of critical points of f can visually be represented by a topological abstraction called the persistence diagram [START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF][START_REF] Edelsbrunner | Topological persistence and simplification[END_REF] (Fig. 1). By applying the Elder Rule [START_REF] Edelsbrunner | Computational Topology: An Introduction[END_REF], critical points can be arranged in a set of pairs, such that each critical point appears in only one pair (c i , c j ) with f (c i ) < f (c j ) and I (c i ) = I (c j ) -1. More precisely, the Elder Rule states that as the value i increases, if two topological features of f -1 -∞ (i), for instance two connected components, meet at a given saddle c j of f , the youngest of the two (the one with the highest minimal value, c i ) dies at the advantage of the oldest (the one with the lowest minimal value). Critical points c i and c j then form a persistence pair.

The persistence diagram D( f ) embeds each pair (c i , c j ) in the plane such that its horizontal coordinate equals f (c i ), and the vertical coordinate of both c i and c j is f (c i ) and f (c j ), corresponding respectively to the birth and death of the pair. The height of the pair P(c i , c j ) = | f (c j )f (c i )| is called the persistence and denotes the life-span of the topological feature created at c i and destroyed at c j . In three dimensions, the persistence of the pairs linking critical points of index (0, 1), (2, 3) and (1, 2) denotes the life-span of connected components, voids and non-collapsible cycles of f -1 -∞ (i). In practice, persistence diagrams serve as an important visual representation of the distribution of critical points in a scalar data-set. Small oscillations in the data due to noise are typically represented by critical point pairs with low persistence, in the vicinity of the diagonal. In contrast, topological features that are the most prominent in the data are associated with large vertical bars (Fig. 1). In many applications, persistence diagrams help users as a visual guide to interactively tune simplification thresholds in multi-scale data segmentation tasks based on the Reeb graph [START_REF] Carr | Simplifying flexible isosurfaces using local geometric measures[END_REF][START_REF] Gueunet | Contour forests: Fast multi-threaded augmented contour trees[END_REF][START_REF] Pascucci | Robust online computation of Reeb graphs: simplicity and speed[END_REF][START_REF] Reeb | Sur les points singuliers d'une forme de Pfaff complètement intégrable ou d'une fonction numérique[END_REF][START_REF] Tierny | Jacobi fiber surfaces for bivariate Reeb space computation[END_REF][START_REF] Tierny | Loop surgery for volumetric meshes: Reeb graphs reduced to contour trees[END_REF] or the Morse-Smale complex [START_REF] Gyulassy | A practical approach to morse-smale complex computation: Scalability and generality[END_REF][START_REF] Gyulassy | Conforming morse-smale complexes[END_REF][START_REF] Robins | Theory and algorithms for constructing discrete morse complexes from grayscale digital images[END_REF]. Wasserstein distance: Metrics have been defined to evaluate the distance between two scalar fields f , g : M → R. The L-norm || f -g|| ν , is a classical example. In the context of topological data analysis, multiple metrics [START_REF] Chazal | Proximity of Persistence Modules and their Diagrams[END_REF][START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF] have been introduced in order to compare persistence diagrams. In our context, such metrics are key to identifying zones in the data which are similar to one another. Persistence diagrams can be associated with a pointwise distance, noted d ν inspired by the L-norm. Given persistence pairs a = (a x , a y ) ∈ D( f ) and b = (b x , b y ) ∈ D(g), d ν is given by Eq. 1:

d ν (a, b) = (|a x -b x | ν + |a y -b y | ν ) 1/ν (1)
The Wasserstein distance [START_REF] Kantorovich | On the translocation of masses[END_REF][START_REF] Monge | Mémoire sur la théorie des déblais et des remblais[END_REF], sometimes called the Earth Mover's Distance [START_REF] Levina | The earthmover's distance is the mallows distance: some insights from statistics[END_REF], noted d W ν , between the persistence diagrams D( f ) and D(g) is defined in Eq. 2:

d W ν D( f ), D(g) = min φ ∈Φ ∑ a∈D( f ) d ν a, φ (a) ν 1/ν ( 2
)
where Φ is the set of all possible bijections φ mapping each critical point a of D( f ) to a critical point b of the same index I in D(g) or to the diagonal, noted diag(a) -which corresponds to the removal of the corresponding feature from the assignment, with a cost d(a, diag(a)).

Assignment problem

The assignment problem is the problem of choosing an optimal assignment of n workers w ∈ W to n jobs j ∈ J, assuming numerical ratings are given for each worker's performance on each job. Given ratings r(w x , j y ) are summed up in a cost matrix (r xy ), finding an optimal assignment means choosing a set of n independent entries of the matrix so that the sum of these elements is optimal. Independent means than no two such elements should belong to the same row or column (i.e. no two workers should be assigned to the same job and no worker should be given more than one job). In other words, one must find a map σ : W → J of workers and jobs for which the sum ∑ x (r(w x , σ (w x ))) is optimal. There are n! possible assignments, of which several may be optimal, so that an exhaustive search is impractical as n gets large.

Similarly, the unbalanced assignment problem is the problem of finding an optimal assignment of n workers to m jobs, where some jobs or workers might be left unassigned. This is the case of assignments between sets of persistence pairs; where costs are defined for leaving specific pairs unassigned.

The Hungarian algorithm [START_REF] Kuhn | The hungarian method for the assignment problem[END_REF][START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF] is the first polynomial algorithm proposed by Kuhn to solve the assignment problem. It is an iterative algorithm based on the following two properties: Theorem 1 If a number is added or subtracted from all the entries of any one row or column of a cost matrix, then an optimal assignment for the resulting cost matrix is also an optimal assignment for the original cost matrix. Subtracting the minimum element from each of the n rows and columns might not be sufficient to make a set of n independent zeros appear. In the above example, initially detected independent zeros are first starred. All columns containing a 0* are then covered (left). An uncovered zero which has a 0* in its row is found and primed; its row is covered and the column of the 0* is uncovered (center). At this point, all zeros are covered by construction. Let ε be the smallest uncovered value. Add ε to every covered row; subtract ε from every uncovered column. This amounts to decreasing uncovered elements by ε and increasing twice-covered elements by ε. The sum of the elements of the matrix has been decreased and a new zero has appeared in an uncovered zone.

Figure 3: Augmenting path phase. After the first non-covered zero (left) is primed and covers updated, there is one non-covered zero Z 1 in the matrix (center), which is then primed. Let Z 2 be the 0* in the column of Z 1 (if any), let Z 3 be the 0' in the row of Z 2 (there is one). Consider the series consisting of 0* (Z 2i ) and 0' (Z 2i+1 ) until it ends at a 0' that has no 0* in its column. Unstar each 0*, star each 0' of the series. The number of starred zeros has increased by one.

This means that the cost matrix (r i j ) can be replaced with (r i j )u iv j where u i (resp. v j ) is an arbitrary number which is fixed for the i th row (resp. the j th column).

Theorem 2 If R is a matrix and m is the maximum number of independent zeros of R (i.e. number of entries valued at 0), then there are m lines (row or columns) which contain all the zeros of R.

This allows to determine whether an optimal assignment has been found and thus constitutes the stop criterion.

The algorithm iteratively performs additions and subtractions on lines and columns of the cost matrix, in a way that globally decreases the matrix cost, until the optimal assignment has been found, that is, until the matrix contains a set of min(m, n) independent zeros.

In the remainder we consider the O(min(m, n) 2 max(m, n)) unbalanced Kuhn-Munkres algorithm [START_REF] Bourgeois | An extension of the munkres algorithm for the assignment problem to rectangular matrices[END_REF][START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF], an improvement over Kuhn's original version which follows the same principles, with an enhanced strategy for finding independent elements. The goal is always to reduce the cost matrix and find a maximal set of independent zeros. These independent zeros are marked with a star: they are candidates for the optimal assignment. Zeros which are candidates for being swapped with a starred zero are marked with a prime. Throughout the algorithm, rows and columns of the matrix are marked as covered to restrict the search for candidate zeros.

The algorithm can be seen as two alternating phases: a matrix reduction phase (Fig. 2) which makes new zeros appear, and an augmenting path phase (Fig. 3) which augments the number of marked (starred) independent zeros.

Persistence assignment problem

The assignment problem for persistence pairs is very similar to the standard unbalanced assignment problem, except additional costs are defined for not assigning elements (i.e. matching persistence pairs with the diagonal). The assignment between diagrams P and Q then involves r i j the numerical rating associated with assigning p i ∈ P with q j ∈ Q, along with r i,-1 (resp. r -1, j ) the numerical rating associated with matching p i (resp. q j ) with the diagonal.

If P and Q are sets of persistence pairs such that card(P) = n and card(Q) = m, then it is possible to solve the persistence assignment problem using the standard Kuhn-Munkres algorithm with the (n + m) × (n + m) cost matrix described by Eq. 3, as proposed in [START_REF] Morozov | Dionysus[END_REF]:

r i j =          d ν (p i , q j ) if 0 < i ≤ n, 0 < j ≤ m d ν (p i , diag(p i )) if n < i ≤ m + n, 0 < j ≤ m d ν (q j , diag(q j )) if 0 < i ≤ n, m < j ≤ m + n 0 if n < i ≤ m + n, m < j ≤ n + m (3)
The first line corresponds to matching pairs from P to pairs from Q; the second one corresponds to the possibility of matching pairs from P to the diagonal; the third one is for matching pairs of Q to the diagonal and the last one completes the cost matrix. The drawback of this approach is that it requires to solve the assignment problem on a (n + m) 2 cost matrix (that potentially contains two non-sparse blocks where persistence elements are located, see Fig. 4), though the number of distinct elements is at most (n + 1) × (m + 1). As seen in Sec. 3, our algorithm addresses this issue.

Overview

This section presents a quick overview of our tracking method, which is illustrated in Fig. 9. The input data is a time-varying PL scalar field f defined on a PL d-manifold M with d ≤ 3.

1. First, we compute the persistence diagram of the scalar field for every available timestep. 2. Next, for each pair of two consecutive timesteps t and t + 1, we consider the two corresponding persistence diagrams D( f t ) and D( f t+1 ). For each couple of persistence pairs (p i , q j ) ∈ D( f t ) × D( f t+1 ), we define a distance metric corresponding to the similarity of these pairs: d ν (p i , q j ) (see Sec. 4). 3. For each pair of consecutive timesteps, we compute a matching function M. Every persistence pair p i of D( f t ) is associated to to M(p i ), which is either a persistence pair q j in D( f t+1 ) or diag(p i ) so as to minimize the total distance ∑ i d(p i , M(p i )).

Finding the optimal M involves solving a variant of the classical Assignment Problem, as presented in Sec. 2.3. Only persistence pairs involving critical points of the same index are taken into account. 4. We compute tracking trajectories starting from the first timestep. If at timestep t the matching associates p i with M t (p i ) = q j , then a segment is traced between p i and q j . If M t (p i ) = diag(p i ), the current trajectory ends. Trajectories are grown following this principle throughout all timesteps. Properties are associated to trajectories (time span, critical index), and to trajectory segments (matching cost, scalar value). 5. Finally, trajectories are post-processed to detect feature merging or splitting events with a user-defined geometric threshold.

OPTIMIZED PERSISTENCE MATCHING

This section presents our novel extension of the Kuhn-Munkres algorithm, which has been specifically designed to address the computation time bottleneck described in Sec. 2.2.

Reduced cost matrix

The classical persistence assignment algorithm based on Kuhn-Munkres considers R, a (n + m) 2 cost matrix. We propose to work instead with R , a reduced (n + 1) × m matrix defined in Eq. 4, where every zero appearing in the last row is considered independent. This amounts to considering that persistence pairs corresponding to rows are not assigned by default. Fig. 4 summarizes the matrices considered by each assignment method.

r i j = d ν (p i , q j ) -diag(q j ) if 0 < i ≤ n, 0 < j ≤ m diag(p i ) if i = n + 1, 0 < j ≤ m (4)
This last row, emulating the diagonal blocks of Fig. 4-b requires a specific handling in the optimization procedure. In particular, it requires the first step of the algorithm to subtract minimum elements from columns (and not rows) so as not to have negative elements in the matrix. As a reminder, the original algorithm proceeds iteratively in two alternating phases: matrix reduction that makes new zeros appear, and augmenting path that finds a maximal set of independent zeros. At the i th iteration, the current maximal set of independent zeros is made of starred zeros. After a matrix reduction, new zeros appeared that can potentially belong to the new maximal set of independent zeros. Such candidates are primed. A single augmenting path (as in Fig. 3) replaces a set of n starred zeros with n + 1 primed zeros, forming a new set of independent zeros with one more element. Rows and columns of the matrix are marked as covered to restrict the search for candidates in the augmenting path phase. Blue blocks of Algorithm 1 indicate our extension of the Kuhn-Munkres algorithm.

In this novel extension, an augmenting path constructed in the corresponding phase can start from a starred zero in the last row (and then potentially find a primed zero in its column), but such a path can never access a starred zero in the last row at another step, for the corresponding column would have been covered prior to this (and thus cannot contain a primed zero, see Algorithm 1). A starred zero in the last row can then never be unstarred.

The Kuhn-Munkres approach has the property to only increase row values (and only decrease column values). When our algorithm working on the reduced matrix R ends, it is therefore not possible that the elements on the top-right corner of the corresponding full matrix R be negative. Furthermore, given Theorem 1, the resulting matrix corresponds to the same assignment problem.

Optimality

Working with the reduced matrix R , however, does not necessarily yield an optimal assignment. When assignments are found in the bottom row, if there has been additions to the matrix rows, then the corresponding R matrix would contain a top-right block that is not zero, and a top-left block that is not zero either. Thus, the stop criterion stated by Theorem 2 may not be respected when k = min(m, n) lines are covered (as the real number of independent zeros in R is m + n). Moreover, in our setup, a starred zero in the last column can never be unstarred; this is allowed in the approach on R, owing to the bottom-right block, initially filled with zeros.

We therefore use the criterion stated in Eq. 5 to ensure that if, at any given iteration of the algorithm, a zero is starred in the last row of column j, the cost of assigning the corresponding persistence pair to any other pair is higher than the cost of leaving both unassigned In our setup, every element in the last row is considered independent, so that it can contain multiple starred zeros (left). This emulates the behavior of the bottom-left matrix block in the classical approach. During an ε-reduction phase (center), we keep track of the (always positive) quantities that were added to matrix rows, hence increasing the top-right block in the classical approach, initially filled with only zeros. If a zero is starred in the last row and j th column, let ρ i be the sum of quantities added to row i throughout the algorithm (right). If for all i, r i j > ρ i , then the persistence pair associated with column i is assigned to the diagonal. If not (which never happened in our experiments), row residuals ρ i and the equivalent residuals for columns ρ j are used to report the partial optimization onto the matrix of the exact classical approach.

(0 for the j th -column pair and the residual value ρ i for i th -row pairs -see Algorithm 1). This specificity is illustrated in Fig. 5.

∀i ∈ 0, n , r i, j > ρ i ⇒ r n+1, j = 0 * (5)

The Eq. 5 criterion is checked whenever a zero appears on the last row after a subtraction is performed on a column by the algorithm. If it is observed, the corresponding column is removed from the problem and the persistence pair is set unassigned.

If the criterion is not respected, we have to report back the reduced problem onto the full matrix (missing banned columns and with reported found residuals ρ). For this, we need to keep track of residuals, that is, values that have been added or subtracted from each row and column throughout the course of the algorithm. Once these residuals have been reported onto the full matrix, there can be no negative element, and all of the optimization work has been reported (so that we do not start all over again from the beginning, but we start from the optimized output of the first phase).

In practice for persistence diagrams, we always observed that the first phase is sufficient to find an optimal assignment. Using this approach prevents from working with two potentially large blocks of persistence elements, typically occurring with the complete matrix for i ∈ n + 1, m + n and j ∈ m + 1, m + n . This property is further motivated by the use of geometrical lifting (Sec. 4). The approach is detailed in Algorithm 1.

Sparse assignment

In practice, it is often observed that some assignments are not possible, and that reordering columns in the associated cost matrix would enable faster lookups and modifications [START_REF] Cui | Solving large-scale assignment problems by kuhn-munkres algorithm[END_REF], using sparse matrices. With persistence diagrams, the following simple criterion (Eq. 6) can be used to discard lookups for potential matchings.

d ν (p, q) > d ν (p, diag(p)) + d ν (q, diag(q)) (6)
Working with our version of the Kuhn-Munkres algorithm then becomes interesting for many assignments verify Eq. 6 (Fig. 6), hence greatly reducing the lookup time for zeros, minimal elements, and the access time for operations performed on rows or columns.

On the contrary, the original full-matrix version of Kuhn-Munkres deals with non-sparse blocks which have to be accessed and modified constantly throughout the course of the algorithm.

LIFTED PERSISTENCE WASSERSTEIN METRIC

This section highlights the limitations of the natural Wasserstein metric applied to time-varying persistence diagrams and presents

Data: R = (r i j ), an (n + 1) × m persistence cost matrix, R the full (n + m) 2 matrix with non-sparse blocks. Result: S a set of starred independent zeros ∀i, ρ i ← 0 // row residuals ∀ j, ρ j ← 0 // column residuals B ← / 0 // banned columns Subtract the persistence element from every row and ρ i Transpose R if n > m and let k = min(m, n) Subtract the min element from every column of R and ρ j Star independent zeros and cover their columns while number of covered columns < k do Find a non-covered zero Z 1 and prime it if Z 1 is in the last row or there is no 0* in its row then Augmenting path phase (Fig. 3) Erase all primes, reset all covers Cover each column of containing a starred zero else Let Z 1 be the 0* in the row of Z 1 Cover this row and uncover the column of Z 1 end if there is no uncovered zero left then Matrix ε-reduction phase (Fig. 2) ρ i ← ρ i + ε for modified rows i ρ j ← ρ jε for modified columns j if ∃ j|r n+1, j = 0 and ∀i ∈ 1, n , r i, j > ρ i then r n+1, j is starred

B ← B ∪ j end end end if ∃ j /
∈ B|r n+1, j = 0 * and ∃i|r i, j < ρ i then Kuhn-Munkres(R i j = R i j + ρ i + ρ j ) with j / ∈ B. end Algorithm 1: Our algorithm for sparse persistence matching. Blue sections allow to emulate the behavior of the three original non-sparse blocks on one single row, while ensuring optimality thanks to the residuals column. Black sections are common with the unbalanced Kunk-Munkres algorithm.

an extension that enhances its geometrical stability. Geometrical considerations are motivated, in terms of accuracy and performance.

Persistence diagrams can be embedded into the geometrical domain (Fig. 1). Doing so, one easily sees how different embeddings can correspond to similar persistence diagrams in the birth-death space. Working in this 2D space does yield irrelevant matchings: as can be seen in Fig. 7, when only the birth-death coordinates of persistence pairs are considered, a matching can be optimal even if it happens between geometrically distant zones. As a consequence, the distance metric between persistence pairs must be augmented with geometrical considerations.

To address this, we propose instead of d ν (Eq. 1) to use the distance defined in Eq. 7:

d li f t,ν (p, q) = (αδ ν birth + β δ ν death + γ 1 δ ν x + γ 2 δ ν y + γ 3 δ ν z ) 1/ν (7)
where δ x , δ y and δ z correspond to geometric distances between the extrema involved in the persistence pairs on a given axis. We process diagonal projections as follows (Eq. 8):

d li f t,ν (p, diag(p)) =(α |p x | ν + β p y ν + γ 1 (δ p x ) ν + γ 2 (δ p y ) ν + γ 3 (δ p z ) ν ) 1/ν (8)
where the terms δ p x , δ p y and δ p z correspond to the geometric distance between the critical points of a given pair p. Intuitively, it accounts for the distance between the critical points to cancel.

A lifted distance is considered by augmenting the geometric distance with coefficients α, β , γ i . This aims at giving more or less showing in color small persistence pairs that will never be assigned in an optimal matching. The light blue pair p ∈ D 2 is such that d(p, diag(p)) + d(q, diag(q)) < d(p, q) for any q ∈ D 1 which is neither light blue nor the first large persistence pair. This results in the cost matrix (right, D 1 pairs are rows and D 2 pairs are columns), where gray elements correspond to pairs (p, q) s.t. d(p, q) > d(p, diag(p)) + d(q, diag(q)). importance to the birth, death or some of the x, y, z coordinates during the matching, depending on applicative contexts. For instance, in practice it is desirable to give less importance to the birth coordinate when dealing with d-(d -1) persistence pairs (in other words, for tracking local maxima, see Fig. 8). For the remainder of the paper and the experiments, we used (α, β , γ i ) = (0.1, 1, 1) for maxima and (α, β , γ i ) = (1, 0.1, 1) for minima, for normalized geometrical extent and scalar values. We observed that using a lifted metric further increases the cost matrix sparsity, resulting in extra speedups.

FEATURE TRACKING

This section describes the four main stages of our tracking framework, relying on the discussed theoretical setup. Without loss of generality, we assume that the input data is a time-varying 2D or 3D scalar field defined on a PL-manifold. Topological features are extracted for all timesteps (Fig. 9, a-b), then matched (Fig. 9,c); trajectories are built from the successive matchings (Fig. 9,d) and post-processed to detect merging and splitting events.

Feature detection

First, we compute persistence diagrams for each timestep. We propose using the algorithm by Gueunet et al. [START_REF] Gueunet | Task-based augmented merge trees with Fibonacci heaps[END_REF], in which only 0- 1 andd-(d -1) persistence pairs are considered.

When the data is noisy, it is possible to discard pairs of low persistence (typically induced by noise) by applying a simple threshold. In practice, this amounts to only considering the most prominent features. Using such a threshold accelerates the matching process, where for approaches based on overlaps, removing topological noise would require a topological simplification of the domain (for example using the approach in [START_REF] Tierny | Generalized topological simplification of scalar fields on surfaces[END_REF]), which is computationally expensive.

Feature matching

If P 1 , P 2 are two sets of persistence pairs taken at timesteps t and t +1, then we use the algorithm described in Sec. 3, with the appropriate distance metric, as discussed in Sec. 4, to associate pairs in P 1 and P 2 . A given pair p 1 ∈ P 1 might be associated to one pair p 2 ∈ P 2 at most, or not associated, and symmetrically.

Trajectory extraction

Trajectories are constructed by simply attaching successively matched segments. For all timesteps t, if the feature matching associates p i with M t (p i ) = q j , then a segment is traced between q j and p i , and is potentially connected back to the previous segment of p i 's trajectory. If M t (p i ) = diag(p i ), the current trajectory ends. Properties are associated to trajectories (time span, critical point index) and to trajectory segments (matching cost, scalar value, persistence value, embedded volume).

Handling merging and splitting events

Given a user-defined geometrical threshold ε, we propose to detect events of merging or splitting along trajectories in the following manner. If T 1 , T 2 : I ⊂ N → R 3 are two trajectories spanning throughout [t i ,t i+n ] and [t j ,t j+m ] respectively, and if for some

k ∈ [i, i + n] ∩ [ j, j + m], d li f t,ν (T 1 (t k ), T 2 (t k )) < ε,
where d li f t,ν is a lifted distance, then an event of merging (or splitting) is detected. We consider that a merging event occurs between T 1 and T 2 at time k, when neither T 1 nor T 2 start at t k . We then consider that the oldest trajectory takes over the youngest. For example, T 1 and T 2 meet (according to the ε criterion) at t k the last timestep of T 2 , and T 2 started before T 1 , then we disconnect the remainder of T 1 from the trajectory before t k and we connect it so as to continue T 2 until T 1 's original end. Similarly, a splitting event occurs between T 1 and T 2 at time k, when neither T 1 nor T 2 end at t k . The process is illustrated in Fig. 10. It is done separately for distinct critical point types: minima, maxima and saddles are not mixed.

RESULTS

This section presents experimental results obtained on a desktop computer with two Xeon CPUs (3.0 GHz, 4 cores each), with 64 GB of RAM. We report experiments on 2D and 3D time-varying datasets, that were either simulated with Gerris [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible euler equations in complex geometries[END_REF] (von Kármán Vortex street, Boussinesq flow, starting vortex), or acquired (Sea surface height, Isabel hurricane). Persistence diagrams are computed with the implementation of [START_REF] Gueunet | Task-based augmented merge trees with Fibonacci heaps[END_REF] available in the Topology ToolKit [START_REF] Tierny | The Topology ToolKit[END_REF]; the tracking is restricted to 0-1 and (d -1)-d pairs. We implemented our matching (Sec. 3) and tracking approaches (Sec. 5) in C++ as a Topology ToolKit module.

Application to simulated and acquired datasets

We applied our tracking framework to both simulated and acquired time-varying datasets to outline specific phenomena.

In Fig. 11, we present the results of the tracking framework applied to an oceanographic dataset. The scalar field (sea surface height) is defined on 365 timesteps on a triangular mesh. We can see interesting trajectories corresponding to well-known oceanic currents. Drifting (a, b) and turbulent behaviors (c) of local extrema are highlighted. In Fig. 12, tracking is performed on the vorticity of highly unstable Boussinesq flow. Thanks to our analysis, trajectories can be filtered according to their temporal lifespan, revealing clearly different trajectory patterns among the turbulent features. This kind of analysis may be easily performed based on other trajectory attributes, depending on applicative contexts. In Fig. 13, we show our approach on a 3D hurricane dataset whose temporal resolution is such that a method based on overlaps of split-tree leaves (see Sec. 6.3) could not extract trajectories. In Fig. 14, our tracking framework correctly follows local extrema of the vorticity field in a simulated vortex street.

Tracking robustness

In the following two sections, we demonstrate the robustness and performance of our tracking framework.

We compare to the greedy approach based on the overlap of volumes [START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topology-based data segmentation[END_REF][START_REF] Bremer | Analyzing and tracking burning structures in lean premixed hydrogen flames[END_REF][START_REF] Saikia | Global feature tracking and similarity estimation in time-dependent scalar fields[END_REF][START_REF] Sohn | Time-varying contour topology[END_REF] of split-tree leaves, which amounts to tracking local maxima. In this approach, for every pair of consecutive timesteps (t,t + 1), split-tree segmentations S t and S t+1 are computed (these are a set of connected regions). Overlap scores are then computed for every pair of regions (s i , s j ) ∈ S t × S t+1 , as the number of common vertices between s i and s j . Scores are sorted and s i is considered matched to the highest positive scoring s j such that s j has not been matched before. Trajectories are extracted by repeating the process for all timesteps.

The robustness of our tracking framework is first assessed on a synthetic dataset consisting of whirling gaussians, on which we applied noise (Fig. 15). Identified trajectories are sensibly the same with a perturbation of 10% of the scalar range. The 75% most important features are still correctly tracked after a 25% random perturbation has been applied to the data. In Fig. 14, our method is compared to the greedy approach, based on overlaps, while decreasing the temporal resolution. The overlap approach yields trajectories corresponding to noise (Fig. 14-e), which can be filtered by applying topological simplification beforehand (this would have a significant computational cost as it requires to modify the original function), or by associating the scalar value of the function to every point in the trajectory and then filtering the trajectory in a post-process step. In our setup, is is much simpler to discard this noise, by using a threshold for discarding small persistence pairs before the matching (implying a faster matching computation). When downsampling the temporal resolution to only 20% of the timesteps, our approach still gives the correct results (Fig. 14,d vs. e). With 15% of the timesteps, our approach (Fig. 14f) still agrees with the tracking on the full-resolution data (Fig. 14-b), until preceding features begin to catch up, resulting in a zig-zag pattern. By comparison, the overlap method fails to correctly track meaningful regions from the beginning of the simulation to its end; it is indeed dependent on the geometry of overlaps, which is unstable. It can be argued that the locality captured by overlaps is emulated in our framework by embedding and lifting the Wasserstein metric, when the overlap method does not take persistence into account when matching regions. Also note that if the saddle component of persistence pairs associated to maxima is ignored (i.e. if α = 0 in Eq. 7 and Eq. 8) during the matching, then the geometrical distance can be insufficient for correctly tracking these persistence pairs (c). Therefore, the problem of matching persistence pairs for tracking topological features cannot be reduced to the (unbalanced) problem of assigning critical points in 4 dimensions (3 for the geometrical extent, one for the scalar value).

Fig. 13 further illustrates the robustness of our approach when downsampling the data temporal resolution. In hurricane datasets, local maxima can be displaced to geometrically distant zones between timesteps if those are taken at multiple-day intervals. This unstable behavior and the absence of obvious overlaps makes it particularly difficult to track extrema; nevertheless, our framework managed to track them at a very low time-resolution.

Tracking performance

We then compare our framework with our implementation of the approach based on overlaps [START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topology-based data segmentation[END_REF] on the ground of performance. Figures are given in Table 1. Note that our approach has the advantage of taking persistence diagrams as inputs, so it can be applied to unstructured or time-varying meshes, for which overlap computations are not trivial. Our approach is also relatively dimension-independent: though in 3D, computing overlaps is very time-consuming (Fig. 1-Isabel), the complexity of the Wasserstein matcher, which only takes With the same temporal resolution, the overlap-based approach (e) does capture small trajectories corresponding to noise, displayed with thinner lines, that have to be filtered for instance using topological simplification [START_REF] Tierny | Generalized topological simplification of scalar fields on surfaces[END_REF]. Considering 1 frame every 7 timesteps (f) still yields correct trajectories up to the point where, every other frame, optimal matchings for the metric are between a feature and the preceding one, due to features traveling fast. The overlap approach (g) is less stable in this case as it extracts erroneous trajectories from the very first stages of the simulation to the end. embedded persistence diagrams as inputs, for a given number of persistence pairs is sensibly equivalent. For both Isabel and Sea surface height datasets, we applied a 4% persistence filtering on input persistence diagrams. As the experiments show, our approach is faster in practice than the overlap method with best-match search.

Matching performance

Next, we compare the performance of the matching method we introduced in Sec. 3 to two other state-of-the art algorithms.

We compare it to the reference approach for the exact assignment problem [START_REF] Weaver | An implementation of the kuhn-munkres algorithm[END_REF] based on the Kuhn-Munkres algorithm, and to our implementation of the approximate approach based on the auction algorithm [START_REF] Bertsekas | The auction algorithm for transportation problem[END_REF][START_REF] Kerber | Geometry helps to compare persistence diagrams[END_REF], on the ground of performance.

Table 2 shows that our new assignment algorithm is up to two orders of magnitude faster than the classical exact approach [START_REF] Weaver | An implementation of the kuhn-munkres algorithm[END_REF]. 16: Running times in seconds of different matching approaches, for decreasing persistence thresholds expressed in percentage of the scalar range. The initial two diagrams containing 14,082 and 14,029 pairs are filtered to remove pairs whose persistence is less than the defined threshold, then a matching is computed with our new method, the reference exact method [START_REF] Weaver | An implementation of the kuhn-munkres algorithm[END_REF], the approximate method [START_REF] Bertsekas | The auction algorithm for transportation problem[END_REF], first with 1% accuracy, then with an accuracy of 10 -4 % of the scalar range.

In particular, the best speedups occur for the larger datasets which indicates that our approach also benefits from an improved scaling.

It is often useful in practice to discard low-persistence pairs prior to any topological data analysis as these correspond to noise. Fig. 16 compares the running times of our approach, [START_REF] Weaver | An implementation of the kuhn-munkres algorithm[END_REF] and [START_REF] Bertsekas | The auction algorithm for transportation problem[END_REF] as more and more low-persistence pairs are taken into account. When removing pairs whose persistence is below 5% of the scalar range, which is commonly accepted as a conservative threshold, approach is faster than all competing alternatives. When considering more low-persistence features, below 4%, our approach is competitive with the approximated auction approach with 1% error. Below 2%, only noise is typically added in the process. The performance of our algorithm becomes comparable to that of the high-precision auction approximation although our approach guarantees exact results.

Limitations

As we described, our framework enables the tracking of 0-1 and d-(d -1) persistence pairs. It would be interesting to extend it to support the tracking of saddle-saddle pairs (in 3D) and see its application to meaningful use cases.

Besides, the lifting coefficients proposed in our metric (Eq. 7) might be seen as supplementary parameters that have to be tuned according to the dataset and applicative domain. Nonetheless, we observed in our experiments that these parameters do not require fine-tuning to produce meaningful tracking trajectories. The extent to which these can be enhanced by fine-tuning is left to future work.

The lifted distance can be generalized to take other parameters, such as the geometrical volume, mass, feature speed, into account, and be fine-tuned to answer the specificity of various scientific domains. Merging and splitting might also be enhanced, or given Table 1: Time performance comparison (CPU time in seconds) between the approach based on overlaps of volumes [START_REF] Bremer | Interactive exploration and analysis of large scale simulations using topology-based data segmentation[END_REF] and our lifted Wasserstein approach. Tracking is performed over 50 timesteps, on structured 2D (Boussinesq, Vortex street), structured 3D (Isabel), and unstructured 2D (Sea surface height) meshes. The preprocessing step (FTM) computes the topology of the dataset. The post-processing step extracts the tracking mesh, computes its attributes, and handles splitting and merging events. We observe a parallel speedup ranging from 4 to 6 for our approach on 8 threads (FTM and matching phases).

Data-set

Pre-proc (s) Matching (s) Post more flexibility, for instance with additional criteria. We also believe that the performance of the post-processing phase can be improved. Additionally, we believe that the approximate auction algorithm can also take the lifted persistence metric into account by performing Wasserstein matchings between persistence pairs in 5 dimensions, and possibly benefit from geometry-based lookup accelerations, as suggested in lower dimension in [START_REF] Kerber | Geometry helps to compare persistence diagrams[END_REF]. It remains to be clarified how the quality of the matchings is affected in practice by using an approximate matching method, and how sparsity can enhance the research phase for the auction algorithm.

We note that the theoretical complexity of our matching method is, as the Munkres method, cubic; however, the two orders of magnitude speedups demonstrated in our experiments allow to study more challenging datasets. For very large case studies, the use of persistence thresholds could prove quite helpful for controlling the computing time of matchings. Among other non-trivial tracking methods, some graph matching methods are based on graph-edit distances [START_REF] Beketayev | Measuring the distance between merge trees[END_REF][START_REF] Gao | A survey of graph edit distance[END_REF]. Their adaptation to the case of persistence diagrams or other topological structures (such as contour trees and Reeb graphs) may enable an additional structural regularization, this ought to be investigated in future work.

CONCLUSION

In this paper, we presented an original framework for tracking topological features in a robust and efficient way. It is the first approach combining topological data analysis and transport for feature tracking. As the kernel of our approach, we proposed a sparsecompliant extension of the seminal assignment algorithm for the exact matching of persistence diagrams, leveraging in practice important speedups. We introduced a new metric for persistence diagrams that enhances geometrical stability and further improves computation time. Overall, in comparison with overlap-based techniques, our approach displays improved performance and robustness to temporal downsampling, as experiments have shown.

We plan to release the implementation of our tracking framework open-source as a part of TTK [START_REF] Tierny | The Topology ToolKit[END_REF] in the near future; we hope that it will be useful to the community with an interest for efficient tracking methods. We look forward to adapting it to tracking phenomena in in-situ contexts, where the large-scale time-varying data is accessed in a streaming fashion. As we are also interested in larger datasets, we are currently carrying out scaling tests on complex physical case studies available at Total S.A., for which one needs specifically adapted rendering techniques [START_REF] Lukasczyk | Nested tracking graphs[END_REF] to apprehend the resulting graphical complexity of the topology evolution.

We also believe that the application potential of our matching framework can be studied for tasks other than time-tracking, for instance, self-pattern matching and symmetry detection [START_REF] Thomas | Multiscale symmetry detection in scalar fields by clustering contours[END_REF], or feature comparisons in ensemble data.

Figure 1 :

 1 Figure 1: Four gaussians defined on a 2D plane (top left) and the associated persistence diagram (right); three gaussians defined on a 3D volume (bottom left) with persistence diagram (right). The most prominent topological features are those which have the longest lifespan, and correspond to the longest vertical bars in persistence diagrams.

Figure 2 :

 2 Figure2: Matrix reduction phase. Subtracting the minimum element from each of the n rows and columns might not be sufficient to make a set of n independent zeros appear. In the above example, initially detected independent zeros are first starred. All columns containing a 0* are then covered (left). An uncovered zero which has a 0* in its row is found and primed; its row is covered and the column of the 0* is uncovered (center). At this point, all zeros are covered by construction. Let ε be the smallest uncovered value. Add ε to every covered row; subtract ε from every uncovered column. This amounts to decreasing uncovered elements by ε and increasing twice-covered elements by ε. The sum of the elements of the matrix has been decreased and a new zero has appeared in an uncovered zone.

Figure 4 :

 4 Figure4: Cost matrices for a balanced assignment problem (left, n × n elements); for a persistence assignment problem with[START_REF] Morozov | Dionysus[END_REF] (center, R with 2n × 2n elements -Eq. 3); and for the same persistence assignment problem with our proposed approach (right, R with (n + 1) × n elements -Eq. 4). Persistence elements in R induce two redundant non-sparse blocks (top-right and bottom-left).

Figure 5 :

 5 Figure5: In our setup, every element in the last row is considered independent, so that it can contain multiple starred zeros (left). This emulates the behavior of the bottom-left matrix block in the classical approach. During an ε-reduction phase (center), we keep track of the (always positive) quantities that were added to matrix rows, hence increasing the top-right block in the classical approach, initially filled with only zeros. If a zero is starred in the last row and j th column, let ρ i be the sum of quantities added to row i throughout the algorithm (right). If for all i, r i j > ρ i , then the persistence pair associated with column i is assigned to the diagonal. If not (which never happened in our experiments), row residuals ρ i and the equivalent residuals for columns ρ j are used to report the partial optimization onto the matrix of the exact classical approach.

Figure 6 :

 6 Figure 6: Persistence diagrams D 1 and D 2showing in color small persistence pairs that will never be assigned in an optimal matching. The light blue pair p ∈ D 2 is such that d(p, diag(p)) + d(q, diag(q)) < d(p, q) for any q ∈ D 1 which is neither light blue nor the first large persistence pair. This results in the cost matrix (right, D 1 pairs are rows and D 2 pairs are columns), where gray elements correspond to pairs (p, q) s.t. d(p, q) > d(p, diag(p)) + d(q, diag(q)).

Figure 7 :

 7 Figure 7: Scalar field f with persistence diagram D( f ) (left), matched with a scalar field g with a similar persistence diagram D(g), but an embedding that swaps the position of the light blue pair with that of the dark blue pair. Matched pairs are displayed with the same color using the non-geometric (center) and geometric Wasserstein metric (right). The latter takes the geometrical embedding into account, preventing similar pairs (regarding persistence) to be assigned if they are geometrically distant.

Figure 8 :

 8 Figure 8: Lifting the birth coordinate. 2D scalar fields with two gaussians (a, b), where the bottom (resp. top) gaussian has the maximum value (a) (resp. b). Using the geometrical metric alone (c) is not sufficient, as the birth coordinate p x misleadingly equalizes the persistence term of pairs of the same color in (a, b): δ a p,yellow = δ b p,yellow , δ a p,blue = δ b p,blue , potentially overcoming the geometrical factor. Lifting the birth coordinate with a small coefficient for associating maxima yields the correct matching (d).

Figure 9 :

 9 Figure 9: Overview of our tracking approach on a dataset consisting of eight whirling gaussians: persistence diagram computations for two consecutive timesteps (a) and (b); matching of persistence pairs of two timesteps (c), propagation of matchings and construction of a trajectory (d).

Figure 10 :

 10 Figure 10: Merging process. Tracking is performed on two gaussians moving from left to right (a). The post-process detects a merging event with a geometrical threshold, then propagates the component identifier of the oldest component (b) and properly reconnects the matching segment. If the oldest component has already the right identifier (c) nothing is done. This process is proposed by analogy with building persistence diagrams.

Figure 11 :

 11 Figure 11: Sea surface height (SSH) captured over 365 days, 1 timestep every day. Local maxima are tracked in the region corresponding to the Agulhas Current, near South Africa (a); it is observed they are slowly drifting towards the west. SSH maxima are also drifting west in the less contrasted zone of the West Australian Current (b). Tracking in the region of the Kuroshio Current, near Japan (c), demonstrates a whirling behavior of local maxima.

Figure 12 :

 12 Figure 12: Boussinesq flow generated by a heated cylinder (a). Feature tracking is performed (b) on the fluid vorticity. Some vortices exist over a long period of time (c), as others vanish more rapidly (d), sometimes akin to noise (e). Feature trajectories can easily be filtered from their lifespan.

Figure 13 :

 13 Figure 13: Tracking performed on the wind velocity on a 3D Isabel hurricane dataset before (left), and after (right) temporal downsampling (1 frame every 5 timesteps). The global maximum is tracked successfully despite the high instability displayed by the scalar field.

Figure 14 :

 14 Figure 14: Simulated von K árm án vortex street (a), on which minima and maxima of the vorticity are tracked with our approach and 1% persistence filtering (b). Only taking the geometry and scalar value into account while doing the matchings (i.e. completely ignoring the birth in the lifted metric), is not sufficient to correctly track features (c). Maxima only are tracked considering 1 frame every 5 timesteps (d).With the same temporal resolution, the overlap-based approach (e) does capture small trajectories corresponding to noise, displayed with thinner lines, that have to be filtered for instance using topological simplification[START_REF] Tierny | Generalized topological simplification of scalar fields on surfaces[END_REF]. Considering 1 frame every 7 timesteps (f) still yields correct trajectories up to the point where, every other frame, optimal matchings for the metric are between a feature and the preceding one, due to features traveling fast. The overlap approach (g) is less stable in this case as it extracts erroneous trajectories from the very first stages of the simulation to the end.

Figure 15 :

 15 Figure 15: Lifted Wasserstein tracking performed on a set of whirling 2D gaussians (a). With noise accounting for 10% of the scalar range (b), feature trajectories are still correctly detected. For 25% noise (c), 75% of the features (namely, the 6 most prominent out of the initial 8) are still correctly tracked despite heavy perturbations.

Figure

  Figure16: Running times in seconds of different matching approaches, for decreasing persistence thresholds expressed in percentage of the scalar range. The initial two diagrams containing 14,082 and 14,029 pairs are filtered to remove pairs whose persistence is less than the defined threshold, then a matching is computed with our new method, the reference exact method[START_REF] Weaver | An implementation of the kuhn-munkres algorithm[END_REF], the approximate method[START_REF] Bertsekas | The auction algorithm for transportation problem[END_REF], first with 1% accuracy, then with an accuracy of 10 -4 % of the scalar range.

Table 2 :

 2 Time performance comparison between the state-of-the-art Munkres-based approach[START_REF] Weaver | An implementation of the kuhn-munkres algorithm[END_REF], and our modified sparse approach.

					-proc
		FTM	[8]	ours	(s)
	Boussinesq	116	75	18	4.7
	Vortex street	45	23	18	2.8
	Isabel (3D)	863	>3k	17	162
	Sea height	568 N.A.	277	113

ACKNOWLEDGMENTS

This work is partially supported by the Bpifrance grant "AVIDO" (Programme d'Investissements d'Avenir, reference P112017-2661376/DOS0021427) and by the French National Association for Research and Technology (ANRT), in the framework of the LIP6 -Total SA CIFRE partnership reference 2016/0010. The authors would like to thank the anonymous reviewers for their thoughtful remarks and suggestions.