FTK: A Simplicial Spacetime Meshing Framework for Robust and Scalable Feature Tracking
Hanqi Guo, David Lenz, Jiayi Xu, Xin Liang, Wenbin He, Iulian R. Grindeanu, Han-Wei Shen, Tom Peterka, Todd Munson, Ian Foster
aa r X i v : . [ c s . G R ] N ov IEEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 1
FTK: A High-Dimensional Simplicial MeshingFramework for Robust and ScalableFeature Tracking
Hanqi Guo,
Member, IEEE,
David Lenz, Jiayi Xu, Xin Liang, Wenbin He, Iulian R. Grindeanu,Han-Wei Shen,
Member, IEEE,
Tom Peterka,
Member, IEEE,
Todd Munson, and Ian Foster
Fellow, IEEE
Abstract —We present the Feature Tracking Kit (FTK), a framework that simplifies, scales, and delivers various feature-trackingalgorithms for scientific data. The key of FTK is our high-dimensional simplicial meshing scheme that generalizes both regular andunstructured spatial meshes to spacetime while tessellating spacetime mesh elements into simplices. The benefits of using simplicialspacetime meshes include (1) reducing ambiguity cases for feature extraction and tracking, (2) simplifying the handling of degeneraciesusing symbolic perturbations, and (3) enabling scalable and parallel processing. The use of simplicial spacetime meshing simplifiesand improves the implementation of several feature-tracking algorithms for critical points, quantum vortices, and isosurfaces. As asoftware framework, FTK provides end users with VTK/ParaView filters, Python bindings, a command line interface, and programminginterfaces for feature-tracking applications. We demonstrate use cases as well as scalability studies through both synthetic data and scientific applications including Tokamak, fluid dynamics, and superconductivity simulations. We also conduct end-to-end performance studies on the Summit supercomputer. FTK is open-sourced under the MIT license: https://github.com/hguo/ftk.
Index Terms —Feature tracking, high-dimensional meshing, distributed and parallel processing, critical points, isosurfaces, vortices. ✦ NTRODUCTION
Feature tracking is a core topic in scientific visualization forunderstanding dynamic behaviors in time-varying simula-tion and experimental data. By tracking features such asextrema, vortex cores, and boundary surfaces, one can high-light key regions in visualization, reduce data to store, and • Hanqi Guo is with the Mathematics and Computer Science Division,Argonne National Laboratory, Lemont, IL 60439, USA.E-mail: [email protected] • David Lenz is with the Mathematics and Computer Science Division,Argonne National Laboratory, Lemont, IL 60439, USA.E-mail: [email protected] • Jiayi Xu is with the Department of Computer Science and Engineering, the Ohio State University, Columbus, OH 43210, USA.
E-mail: [email protected] • Xin Liang is with Oak Ridge National Laboratory, Oak Ridge, TN 37830,USA.E-mail: [email protected] • Iulian R. Grindeanu is with the Mathematics and Computer ScienceDivision, Argonne National Laboratory, Lemont, IL 60439, USA.E-mail: [email protected] • Wenbin He is with Bosch Research North America, Sunnyvale, CA 94085,USA.E-mail: [email protected] • Han-Wei Shen is with the Department of Computer Science and Engi-neering, the Ohio State University, Columbus, OH 43210, USA.E-mail: [email protected] • Tom Peterka is with the Mathematics and Computer Science Division,Argonne National Laboratory, Lemont, IL 60439, USA.E-mail: [email protected] • Todd Munson is with the Mathematics and Computer Science Division,Argonne National Laboratory, Lemont, IL 60439, USA.E-mail: [email protected] • Ian Foster is with the Data Science and Learning Division, ArgonneNational Laboratory, Lemont, IL 60439, USA.E-mail: [email protected] enable further analysis based on the dynamics of features inscientific data.This paper introduces a general framework that deliversa collection of feature-tracking tools to end users, scalesfeature-tracking algorithms in distributed and parallel envi-ronments, and simplifies the development of new feature-tracking algorithms. The motivations for developing thisframework are threefold. First, although feature-trackingcapabilities appear sporadically in today’s data analysis andvisualization tools, a general-purpose toolset is lacking thatwould enable users to track and analyze features in scientificworkflows. In community tools such as VTK [1], VTK-m [2],ParaView [3], VisIt [4], and TTK [5], most algorithms focuson single-timestep data, and only a few filters are providedfor tracking features over time. Object tracking for videos isavailable in computer vision libraries such as OpenCV [6],but scientific data differ significantly from natural videosin their feature definitions and data representation. Second,few existing feature-tracking algorithms are designed forscalability and parallel processing. The advent of exascalecomputing means that data produced by supercomputersneed to be efficiently handled by the same scale of comput-ing resources. In both in situ and post hoc scenarios, thesheer data size and high complexity of tracking algorithmsnecessitate distributing data to many computing nodes andusing GPU accelerators when available. Third, no devel-oper framework exists for eliminating redundant efforts toimplement application-specific feature-tracking algorithms.Implementing feature tracking algorithms from scratch canbe daunting; the management of time-varying data, thehandling of degenerate cases, and the parallelization oftracking algorithms are needed in many applications; suchfeatures do not exist in publicly available software libraries.
EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 2
To these ends, we identify the common ground—high-dimensional meshing—among many tracking algorithmsfor isosurfaces, critical points, and vortex cores. By extrud-ing the spatial mesh into the time dimension, a spacetimemesh connects the cells in the spatial mesh over adjacenttimesteps. For example, in 3D isosurface tracking, marchingcubes [7] are generalized to higher dimensions [8, 9] by iter-ating and classifying 4D spacetime cells with lookup tables.In critical point tracking [10, 11], the movement of criticalpoints can be captured by identifying the spatiotemporalcells that contain critical points. Likewise, in tracking quan-tum vortices in complex-valued scalar fields [12, 13], themoving trajectories of vortex core lines can be reconstructedin spacetime meshes.We present the Feature Tracking Kit (FTK), which in-troduces high-dimensional simplicial meshing for robust andscalable feature tracking. Compared with spacetime meshesused previously, the key difference of our method is thatall mesh elements are simplices . Previous feature trackingmethods extruded 2D triangles into 3D prisms [10, 13] and3D cubes into 4D cubes [8, 12]; but neither a prism nor acube in the extruded mesh is simplicial.Simplicial meshes offer three benefits for feature track-ing: specificity, stability, and scalability. First, simplicialmeshes eliminate ambiguities in feature tracking, similar tohow marching tetrahedra [14] eliminates isosurface ambi-guity. In nonsimplicial cells such as cubes, multiple featuresintersect the same cell, causing ambiguities that require at-tention. We show that with the spacetime piecewise linearity(PL) assumption, no disambiguation is needed for trackingcritical points and isosurfaces in simplicial meshes. Sec-ond, simplicial meshes ease the handling of degeneracies,enabling robust feature tracking. Degeneracies, such as acritical point on an edge or isosurface intersecting a vertex,may lead to loss or duplication of the detection resultsdue to numerical instabilities [15]. Simplicial meshes enablethe use of Simulation of Simplicity (SoS) [16]—a matureprogramming technique to simplify the handling of de-generacies in computational geometry—to generate robust,combinatorial, and consistent tracking results regardlessof numerical instabilities. Third, simplicial meshes makeit straightforward to accelerate feature-tracking algorithmswith both GPU parallelism and distributed parallelism. Incases when the feature detection is independent in each cell,we can distribute the tasks to different computing resourcesfor concurrent and scalable processing.In this study, we design and implement the simplicialsubdivision of two types of high-dimensional meshes— ( n +1) -D simplicial prism and ( n +1) -D regular meshes—toenable robust and scalable feature tracking in both unstruc-tured and regular meshes in R n , in order to support thetracking of critical points (0D features in 2D/3D), quantumvortices (1D features in 3D), and isosurfaces (2D featuresin 3D) for a wide range of applications. The n -dimensionalspace consists of both 2D/3D space and time, and allmesh elements are simplices. Enabled by high-dimensionalsimplicial meshing, each individual tracking algorithm hasnovelties in disambiguation, degeneracy handling, and scal-ability. We also enable efficient mesh element traversal overtime. Considering that time-varying data are large andstreamed from simulations in situ, one can iterate spacetime mesh elements within a sliding window of a few timestepsfor out-of-core and streaming data access.As a software framework, FTK provides ParaView plu-gins, Python bindings, command line interfaces, and pro-gramming interfaces for end users to track a variety offeatures both in situ and post hoc. We demonstrate the useof FTK for fluid dynamics, fusion, and superconductivitysimulations as well as high-speed imaging experiments. Insummary, the novelty of this paper is in its combination ofseveral individual technical contributions: • A high-dimensional simplicial meshing scheme thatgeneralizes and tessellates both regular and unstruc-tured spatial meshes to spacetime simplices (Sec-tion 3) • A robust and scalable critical point tracking algo-rithm that handles degeneracies in a consistent man-ner with no ambiguities (Section 4) • A scalable implementation of quantum vortex track-ing with distributed parallelism (Section 5) • A robust and scalable isosurface tracking algorithmthat avoids ambiguities and handles degeneracies ina consistent manner (Section 6) • A software framework for users to track featureswith distributed and parallel environments both insitu and post hoc (Sections 7 and 8) • Comprehensive performance studies of FTK algo-rithms on both supercomputers and commodityhardware (Section 9)
ELATED W ORK
This section reviews related work in the extraction andtracking of critical points, quantum vortices, and isosur-faces. We also briefly review high-dimensional meshing.In the following, we refer to the process of independentlydetecting features in individual timesteps as feature extrac-tion, as opposed to feature tracking, which is the process ofreconstructing trajectories of features through multiple con-secutive timesteps. For a comprehensive review of featureextraction and tracking, see Post et al. [17]; Heine et al. [18]review topology-based methods for visualization.
In general, a critical point is defined as the location where avector field vanishes. Critical points are also defined in thegradient fields of scalar functions. Below, we review criticalpoint extraction and tracking algorithms.
Numerical methods have been used to locate critical pointswhere all vector components are zero simultaneously, as-suming the vector field can be interpolated based on discreterepresentations. While finding such zero-crossings has beenstudied for bilinear and trilinear schemes [19], the piecewiselinear (PL) interpolation of vector fields is more widely usedin various applications because of its simplicity [20–22]. Tric-oche et al. characterized higher-order critical points in 2DPL vector fields by partitioning neighboring regions basedon different flow behavior [21]. That approach was furthergeneralized to 3D vector fields [23]. Besides PL, extraction
EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 3 of critical points in piecewise-constant vector fields can beachieved by discrete Hodge decomposition [24].A major issue with numerical methods is their sensitivityto numerical instabilities. As illustrated in Figure 1(a) and(b), a critical point may be identified multiple times if thecritical point resides on the boundary of cells. To this end,Bhatia et al. [15] introduced the use of symbolic perturba-tion [16] in testing if a simplicial cell contains critical points,leading to a robust critical point detection in a combina-torial manner, illustrated in Figure 1(c). Our study furthergeneralizes the use of symbolic perturbation to ensure thatthe tracking of critical points is robust and combinatorial, asdemonstrated in following sections. (a) (b) (c)
Fig. 1. Nonrobust (a, b) versus robust (c) critical point extraction. Withnumerical methods, when a critical point resides on an edge, the criticalpoint may or may not be detected by all triangles that share the sameedge; in this case, the number of detected critical points could rangefrom zero to six because of numerical instabilities. With robust criticalpoint extraction and tracking, the single critical point will be detectedand associated with one of the triangles in a combinatorial manner.
Topology methods include the use of the Poincar´e indextheorem, (discrete) Morse theories, and contour trees/Reebgraphs. For example, Poincar´e’s index theorem can be usedto test if critical points exist in 3D regions [25, 26] or PLsurfaces [27]. With Morse decomposition, Chen et al. [28]proposed a vector field topology representation of 2D PLvector fields with graphs, such that critical points can beidentified as part of the vector field topology. For scalarfields, critical points are the key constituents of the scalarfield topology, including Reeb graphs [29] and contourtrees [30] extracted with well-established algorithms.
Spacetime meshing methods perform interpolation overtime to track critical points. Tricoche et al. [10] extruded 2Dtriangular cells into 3D spacetime prisms, detected entriesand exits of singularities on prism faces, and then identifiedpaths of critical points. Garth et al. [11] generalized this ap-proach to 3D by extruding from tetrahedra to 4D tetrahedralprisms. For both methods, the vector field is assumed linearover both space (barycentric interpolation) and time (linearinterpolation). Note that similarly to bilinear and trilinearinterpolations, such spacetime interpolation in prisms is notPL in spacetime, as explained in the next section.
Feature flow field (FFF) methods [31] use a derived FFFvector field to characterize feature movements, such thatfeature trajectories can be computed as FFF streamlines. Forcritical point tracking, one needs to find an appropriateset of critical points in spacetime as the seeds, computestreamlines from the seeds by numerical integration, andthen slice streamlines back into individual timesteps. To address instabilities in the numerical integration, Weinkaufet al. [32] proposed a method to improve convergence. Kleinand Ertl generalized FFF to track critical points in scalespace [33]. Reininghaus et al. [34] proposed a combinatorialversion of FFF based on the discrete Morse theory to trackcritical points in 2D time-varying scalar fields.
Nearest-neighbor and region-overlapping approaches are heuristics to track critical points. For example, Wanget al. [35] reconstructed critical point trajectories by join-ing proximal and same-type critical points in adjacenttimesteps. Skraba and Wang [36] used the closeness ofrobustness, the minimum amount of perturbation neededto cancel features, to link corresponding critical points inadjacent timesteps.
We use quantum vortices as an example of tracking 1Dfeatures. Quantum vortices, or simply vortices, are topo-logical defects in superconductivity [37], superfluidity [38],and Bose-Einstein condensates. Simulations produce 3Dcomplex-valued fields that combine both amplitudes andphase angles. Singularities in phase fields are closed 1Dcurves embedded in 3D Euclidean spaces. By definition, avortex is the locus of points such that − I C ∇ θ ( x ) · d l = 2 kπ, k = 0 , (1)where θ ( x ) is the phase field; C is an infinitesimal contourthat encircles the vortex curve; d l is the infinitesimal arc on C ; and k is a nonzero integer usually equal to ± . Vortex extraction.
Based on the definition, a straight-forward approach to extract vortices in 3D meshes is tofirst check whether the contour integral is nonzero for eachface boundary and then to trace singularity curves alongfaces [12, 37]. Guo et al. [13] proved that a triangulatedmesh cell intersects up to one singularity line, and thusthat simplicial mesh subdivision leads to combinatorial andconsistent extraction results.
Vortex tracking.
A spacetime meshing approachwas proposed to correspond vortex curves in adjacenttimesteps [12, 39] by examining faces in both space and time.As a result, mesh faces testing positive for a singularity formgraphs that characterize the movement of singularities assurfaces. Guo et al. [13] used triangular/tetrahedral prismsas the spacetime cells to extract and track singularities.However, ambiguities still exist because spacetime prismsare nonsimplicial. In Section 5, we demonstrate that a high-dimensional simplicial mesh eliminates ambiguities in aconsistent manner and allows parallel vortex curve trackingin distributed environments.Quantum vortices are fundamentally different from vor-tices in fluid flows [40], which are swirling centers of flowsand have been defined by level sets or extremum lines of λ eigenvalues [41] and vorticity magnitude [42]. Depend-ing on definitions, tracking of fluid flow vortices may beachieved by connected component labeling [43] or FFF [31]. Isosurface extraction —the task of reconstructing polygonsurfaces with a given isovalue in a 3D scalar field—is fun-damental to scientific visualization. The marching cubes [7]
EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 4 algorithm extracts isosurfaces in regular grid data based onlookup tables; disambiguating how surfaces are connectedinside a cube was the key research problem for a decade.A cubic cell has = 256 possible ways to intersect anisosurface, which boil down to 15 unique configurations.Ambiguities exist when vertex values have alternating signson any faces. Marching tetrahedra [14] is a promisingmethod to eliminate ambiguities by tessellating inputs intosimplicial cells; there are only two unambiguous cases ofintersections for each tetrahedron. Isosurface tracking.
Two distinct approaches exist fortracking isosurfaces: volume tracking [43, 44] and higher-dimensional isosurfacing [8]. The former extracts regions ofinterest independently in each timestep and then associatesregions across adjacent timesteps based on overlaps. The lat-ter generalizes marching cubes to 4D spacetime; the outputsare 3D objects embedded in 4D and can be sliced back into2D surfaces in individual timesteps for visualization. Similarto marching cubes in 3D, ambiguities exist in 4D spacetime,and researchers have shown that ione can disambiguate 4Dcases by triangulation [9]. We demonstrate in Section 6, asimplified implementation of 4D isosurface tracking basedon our high-dimensional simplicial meshes.
Work in high-dimensional meshing has focused mostly on4D data because of the dimensionality of the physics, whilemany algorithms such as Delaunay triangulation work inarbitrary dimensions.
High-dimensional meshing for computational sci-ences.
Recently, scientists have started to explore the use of4D meshes [45, 46] to numerically solve time-dependent par-tial differential equations (PDEs) in spacetime as opposedto traditional timestepping approaches. We believe that ourmethod could be directly applied to spacetime mesh data;but because of challenges of increased complexity, memoryfootprint, and cost to converge, the majority of scientific datatoday are still stored and represented in discrete timesteps.
High-dimensional meshing for scientific visualization.
Spacetime meshing approaches, which are limited mostlyto prisms to date, have been successfully used to tracksingularities in vector fields [10, 11] and phase fields [12, 13,37, 39]. Prisms are a straightforward choice but challengesexist in handling ambiguities, as discussed in later sections.
IGH - DIMENSIONAL SIMPLICIAL MESH
We design and implement the simplicial subdivision of ( n + 1) -D simplicial prism and ( n + 1) -D regular meshes,respectively, in order to enable robust and scalable featuretracking in unstructured and regular grid meshes, wherethe dimensionality n is 2 or 3 for the spatial domain. Theadditional dimension is time in this study, and we assumethat the spatial mesh does not change over time. The sim-plicial subdivision of ( n + 1) -D regular meshes, which is aspecial case of the subdivision of ( n + 1) -D prism meshes, isimplemented separately for the efficient handling of images,volumes, and curvilinear grids.For example, in the case of n = 2 , the input is a (time-invariant) triangular mesh (illustrated in Figure 2(a)). One can extrude the mesh into 3D by replicating and elevat-ing vertices in the new dimension, forming 3D triangularprisms (Figure 2(b)). The output 3D mesh is a subdivisionof triangular prisms, and each mesh element in the outputmesh is simplicial (Figure 2(c)). We also categorize andindex simplices in all dimensions for efficient traversal (Fig-ure 2(d)). In the rest of this section, we formalize definitions(Section 3.1), describe the subdivision of ( n +1) -D simplicialprism meshes (Section 3.2), and introduce the subdivisionof ( n + 1) -D regular meshes as a special case of subdividingprism meshes (Section 3.3). Formally, an n -simplex is the convex hull of n + 1 points a , a , . . . , a n in R n . An n -simplicial complex is the set of k -simplices ( k = 0 , , . . . , n ); the intersection of any two k -simplices ( k > ) is either a common ( k − -simplicial face(such as triangles, edges, and vertices) or the empty set. Forthis study, all simplices are nondegenerate; that is, each k -simplex bounds a k -dimensional volume.We define the ( n + 1) -D simplicial prism as the ex-trusion of an n -simplex a a . . . a n to one dimensionhigher. Denoting the R n +1 coordinates of each point a i as x i = ( x i, , x i, , . . . , x i,n ) ⊺ with the identical last compo-nent x i,n for all i , the extruded prism includes anothersimplex b b . . . b n ; the coordinates of each point b i are ( x i, , x i, , . . . , x ′ i,n ) ⊺ , x i,n < x ′ i,n . Note that a i and b i sharethe same first n coordinates and that the last coordinateis different. In addition to the two simplices, the prismincludes n edges a b , a b , . . . , a n b n . Note that an ( n +1) -Dsimplicial prism is nonsimplicial.We further define the ( n + 1) -D simplicial prism mesh as a collection of ( n + 1) -D simplicial prisms obtained byextruding a conforming simplicial mesh into one dimensionhigher. Our goal is to tessellate the ( n + 1) -D simplicialprism mesh into a simplicial complex without adding newvertices. The tessellation must be conformal , such that anytwo k -simplices either share a common face or have nothingin common. We shown examples of conformal and non-conformal subdivisions in Figure 3. n +1)-D simplicial prismmeshes We first review the concept of staircase triangulation [47] andthen generalize the staircase triangulation to the conformalsubdivision of prism meshes. Without loss of generality,we describe the case of n = 3 , the extrusion from anunstructured 2D triangular mesh to a 3D triangular prismmesh, followed by its subdivision into a 3D tetrahedralmesh. Assuming the input is given by a list of triangles(2-simplices), our algorithm extrudes each triangle into aprism in a new dimension; each triangular prism is furthersubdivided into three tetrahedra in a conformal manner. Staircase triangulation of a 3D triangular prism.
Asillustrated in Figure 2(d), a triangular prism consisting oftwo sets of three vertices, each lying in a plane, may besubdivided into three tetrahedra. For example, consider onetrio of points lying in the plane z = 0 and the other lyingin the plane z = 1 . Denoting the “lower” vertices as a a a and the “upper” vertices as b b b , we may subdivide this EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 5 a0 a4a5 a6a7 a8a9a1 a2 a3 (d) (e)(a) (b) (c) t y p e - I t e t ( b o tt o m ) t y p e - II t e t ( c e n t e r ) type-III tet(top)type-I face(prism base)type-II face(prism lower)type-III face(prism lower)type-V face(edge upper)type-IV face(edge lower) t y p e - I e d g e t y p e - III e d g e t y p e - II e d g e a0 a4a5 a6 a7 a8 a9a1 a2 a3b0 b4b5 b6 b7 b8 b9b1 b2 b3 Fig. 2. Extrusion of simplicial mesh: (a) the input 2D simplicial mesh, (b) the extruded 3D simplicial prism mesh, (c) the output 3D simplicial meshas the conformal subdivision of the 3D simplicial prism mesh, (d and e) subdivision of a 3D triangular prism. Unique types of edges, faces, andtetrahedra in the extruded mesh are illustrated in d and e. (a) (b) (c) (d)
Fig. 3. Conformal (a and c) and nonconformal (b and d) subdivisions ofneighboring triangular prisms and 3-cubes, respectively. a a a b b b a a a b b b a a a b b b Fig. 4. All monotone paths from a to b in the triangular prism a a a − b b b ; each path corresponds to a staircase triangulation of the prism. triangular prism into three tetrahedra: a a a b , a a b b ,and a b b b . As documented by DeLoera et al. [47], thevertex list of each tetrahedra corresponds to a monotonestaircase beginning with a and ending with b in Figure 4;each vertex is immediately above or to the right of theprevious vertex in the grid. Staircase triangulation of an ( n + 1) D simplicial prism.
The staircase triangulation can be generalized to higherdimensions, and an ( n + 1) -D simplicial prism may besubdivided into n + 1 ( n + 1) -simplices without the in-troduction of new vertices. First, we impose an orderingon the n + 1 vertices in the lower and upper hyperplanes(and use the same ordering in both hyperplanes). Second,we identify the n + 1) points of the prism with the grid { , , , . . . , n } × { , } . Third, we compute all monotonepaths on the grid. The staircase triangulation of 2-, 3-, and4-simplicial prisms are listed in Table 1. Staircase triangulation of ( n +1)-D simplicial prismmesh. The staircase subdivision method produces confor-mal subdivisions along prism boundaries, given a global
1. Monotone here means that both alphabets and subscripts areascending. For example, an edge like b a or a a cannot appear. TABLE 1Monotone staircases and triangulation of 1D, 2D, 3D (triangular), 4D(tetrahedral) and 5D simplicial prisms.
1D 2D 3D 4D 5D a b a a b a a a b a b b b b a b b b b b a b b a a b b a a b b b a a b b b b a b b b a a a b b a a a b b b a a a a b a a a a b b a a a a a b a b a b a b a a a a b b b b a a b a b b (4-prism figureomitted) vertex ordering on a simplicial prism mesh. In a 3D casein Figure 2(a-c), we assign a global order to each vertexand then subdivide each prism with staircase triangulation.In practice, this approach is equivalent to subdividing eachquadrilateral into 3D triangular prism meshes. For example,the quadrilateral a a b b is subdivided into two triangles a a b and a b b along the monotonous edge a b . Mesh element indexing.
Each k -simplex in the sub-divided ( n + 1) -D prism mesh can be one to one mappedto a tuple of integer ID, type, and timestep for traversingand compact storage. Considering the extrusion along thenew dimension for multiple layers of vertices (e.g., multipletimesteps), we use the same triangulation scheme for eachlayer and design an efficient indexing of simplices in alldimensions in the new mesh. For k = 3 , there are threetypes of 3-simplices: bottom, center, and upper tetrahedra(or type-I, II, and III tetrahedra), such that one can indexeach tetrahedron with a tuple of original triangle ID, type,and timestep. The original triangle ID is the integer indexof the triangle in the original mesh. For k = 2 , to uniquelyindex 2-simplices, we identify five unique types of faces:prism base, prism lower, prism higher, edge lower, andedge upper. For example, the “top” triangle of a prism canbe indexed by the “bottom” triangle of the same prism inthe next timestep; triangles on quadrilaterals can also beindexed by neighboring prisms in the same layer. As such,each 2-simplex can be uniquely indexed by the tuple oforiginal triangle/edge ID, type, and timestep. Likewise, for k = 1 , we identify three unique types of edges; each can beindexed by the original vertex/edge ID, type, and timestep. Mesh element queries.
We provide functions suchas vertices() , sides() , and side_of() that featuretracking algorithms can use to query a mesh element in the EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 6 extruded mesh. The vertices() function returns the list ofvertices of the given mesh element ID. The sides() func-tion provides a list of ( k − -simplicial sides of the given k -simplex, such as the triangular faces of a tetrahedron. Incontrast, the side_of() function gives a list of ( k + 1) -simplices that contains the given k -simplex. For example, atriangular face in 3D simplicial meshes is usually containedby two tetrahedra unless the face is on the boundary of thedomain; likewise, a tetrahedron in a 4D simplicial mesh isusually shared by two pentachora (4-simplices). Ordinal and interval mesh elements.
For ease of fea-ture tracking, we further categorize k -simplices into ordinal and interval types. A simplex is an ordinal type if each of itsvertices resides in the same timestep in the extruded mesh;otherwise it is an interval type. For example, the lower edgetriangle (type-V face in Figure 2(d)) has two vertices in thelower layer and one vertex in the upper layer and is thusan interval type. The type-I edge is ordinal because bothvertices are in the same layer.There are two reasons to distinguish ordinal and intervaltypes. First, this distinction allows feature-tracking algo-rithms to consume data in a streaming and out-of-core man-ner for both in situ and post hoc processing, as discussed inthe next paragraph. Second, the distinction allows efficientslicing of output trajectories. If one needs only individualtimesteps, rather than intervals, it is straightforward toselect features identified in ordinal mesh elements.In a streaming pipeline, assuming data of each timestep , , . . . , n t − are available in the ascending order, n t beingthe number of timestep, we show that one can traverse all k -simplices while keeping a sliding window of two timestepsof data. For each i th timestep, we first traverse ordinal typesand then traverse interval types if i < n t − . Because eachinterval type consists of vertices from adjacent timesteps,both the i th and ( i + 1) th timesteps must be available inmemory. As a result, it streams timesteps and traverses all k -simplices without having spacetime data of all timestepsin memory simultaneously. Complexity.
The space complexity of maintaining aconformal subdivision of an ( n + 1) -simplicial prism meshis the order of the number of simplices in all dimensions inthe original n -simplicial mesh. For example, in the case of n = 2 , we need to maintain lists of all triangles, edges, andvertices of the original triangular mesh. We also maintainlists of sides and parents for all simplices in the originalmesh, in order to accelerate the query of sides and parentsin the extruded mesh. The time complexity of querying asimplex and getting the vertex list of the simplex is constant. n +1)-D regular mesh Subdividing a regular mesh is a special case of that in thepreceding subsection but does not require maintaining amesh data structure (e.g., lists of vertices and triangles). Wedefine the ( n + 1) -D regular simplicial mesh as a simplicialsubdivision of the ( n + 1) -D regular mesh without introduc-ing additional vertices. Recursive subdivision of ( n + 1) -D regular mesh. Conceptually, one can recursively subdivide an ( n + 1) -Dregular mesh based on the simplicial extrusion of an n -D regular simplicial mesh. For n = 0 , the extruded mesh (a) (b)(c) (d) Fig. 5. Regular mesh (a and c) and regular simplicial mesh (b and d),each in 3D and 4D, respectively. In the 4D cases (c and d), we showthe projection of a single 4-cube. Numbers encode both indices andcoordinates of vertices. (1D regular grid) is already simplicial. For n ≥ , one canextrude cells in a n -D regular mesh into prisms and followTable 1 to triangulate the prisms. As a result, each n -D cubeis subdivided in the same way into n ! congruent and disjoint n -simplices, and the subdivision is conformal. Precomputation of the subdivision for n -D unit cube. In practice, we precompute the subdivision of the n -D unitcube for any n , which enables direct access to an n -Dsimplicial mesh without recursive subdivision. For example,the unit 2-cube can be subdivided into two 2-simplices: , , , , , (2)where each 2-digit is the coordinate and ID of the vertex,and each line is a 2-simplex. By extruding the simplices, theconformal simplicial subdivision of the unit 3-cube containssix tetrahedra: , , , , , , , , , , , , , , , , , , , (3)as is also illustrated in Figure 5(a). Likewise, the unit 4-cubecan be subdivided into twenty-four (
4! = 24 ) pentachora, asillustrated in Figure 5(b). , , , , , , , , , , , , . . . , , , , . (4)After the precomputation, each cube in the n -D regularmesh is subdivided in the same way. EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 7 (0,0) (1,0) (2,0)(0,1) (1,1) (2,1)(0,2) (1,2) (2,2) (a) (c)(b) , { } sides )( = side_of )( , { } = Fig. 6. Indexing and querying mesh elements in regular simplicial mesh:(a) a 2D regular simplicial mesh; (b) indexing simplicial elements with k : ( i i ) / type , where k is the dimensionality of the simplex, ( i i ) is the corner coordinates of the cube that contains the simplex, and type is the unique simplex type ID; (c) sides() and side_of() of asimplicial element. Mesh element indexing.
We use the tuple of simplicialdimension k , corner coordinates, and the unique type ID toindex a k -simplex in the n -D simplicial mesh. The cornercoordinates encode the location of the n -cube that containsthe simplex. The unique type ID is designed to encode asimplex within the cube; one cannot use other cubes to indexthe same simplex. In a 2D case in Figure 6(a), although thereare five edges (1-simplices) in each 2-cube, we have onlythree unique types of edges, because horizontal and verticaledges are always shared between neighboring cubes. Forexample, to index the top edge of cube (1, 1), we can use thecube (1, 2) to find the same edge. Figure 6(b) enumeratesall unique simplex types in 2-cubes. Table 2 enumeratesthe number of unique types up to five dimensions. If thedomain is bounded, mesh element IDs can be one to onemapped to integers for iterating and traversing. TABLE 2Number of unique k -simplextypes in an n -D regularsimplicial mesh, for different n . n k Mesh element queries.
We also provide the vertices() , sides() ,and side_of() functionsdefined in the previoussubsection for n -D regularsimplicial meshes. Theresults of each function areprecomputed for each uniquetype. As illustrated in the2D mesh in Figure 6(c), thesides of Type-II 1-simplices(diagonal edges) include two vertices; two triangular cellscontain the same Type-I 1-simplices (horizontal edges). Complexity.
The space complexity of maintaining an n -D regular simplicial mesh is O ( n !) , but n does not exceed4 for tracking features for 3D data. Note that the spacecomplexity does not grow with the size of the regular grid,because precomputed unit-cube subdivisions are stored in-stead of an explicit list of mesh elements. Such implicit meshdata structure allows queries of a simplex in constant time. RACKING FEATURES : CRITICAL POINTS
We describe here the use of high-dimensional simplicialmeshes to track critical points in 2D and 3D vector fields.
We assume that the input n -dimensional time-varying vec-tor field v : R n +1 → R n is piecewise linear (PL). The PL f (a) (b) Fig. 7. Two-pass critical point trajectory reconstruction for 2D and 3Dvector fields with 3D (a) and 4D (b) spacetime simplicial meshes,respectively. The first pass tests whether triangular/tetrahedral sidesintersects the trajectory, and then the second pass associates every pairof intersected sides if they share the same tetrahedron/pentachoron. assumption implies that the vector field is defined on a sim-plicial spacetime mesh; each cell is an ( n + 1) -simplex. Forexample, each cell in a 2D or 3D time-varying vector fieldis a tetrahedron (3-simplex) or pentachoron (4-simplex),respectively. The ( n + 1) -simplicial spacetime mesh can beconstructed based on an existing n -dimensional mesh, asdetailed in Section 3. Thus, v is C continuous along anycombination of spatial and temporal directions; that is, thereexist A ∈ R n × ( n +1) and b ∈ R n for each ( n + 1) -simplex S such that v = Ax + b , x ∈ S .We further assume that the time-varying vector fieldis generic . This means vector values on each vertex i arenonzero ( v i = ), and vectors at verticies of any k -simplex( k = 0 , , . . . , n + 1 ) are affinely independent. Thus, criticalpoints in generic vector fields may be found in the theinterior of n -simplices instead of on cell boundaries. Becausevector fields in scientific applications may not be generic,in the end of this subsection we discuss the relaxation ofthe generic assumption by using the simulation of simplic-ity [16] technique.A (spacetime) critical point x c ∈ R n +1 is the locationwhere the vector value v ( x c ) is zero. We focus on criticalpoints that are nondegenerate; meaning that the (spatial) Ja-cobian J v at x c is nondegenerate. Based on the eigensystemof J v , the critical point x c can be further categorized intovarious types such as sources, sinks, and saddles. In thecase that v is the gradient field of a scalar field, the criticalpoint types are maxima, minima, and saddles.A critical point trajectory (or simply trajectory) is a locusof critical points in space and time, which are 1D curves em-bedded in R n +1 . Because v is PL, critical point trajectoriesare PL parametric curves, and the intersection with each cellis a line segment. Critical point trajectories can be sliced intoa set of critical points at an arbitrary time t by intersectingthe hyperplane t = t . In the following sections we discussmethods for reconstructing critical point trajectories. As illustrated in Figure 7, we use a two-pass algorithm to di-rectly reconstruct critical point trajectories from v . Withoutloss of generality, we describe this algorithm with 2D time-varying vector fields ( v : R → R ). In the first pass, weiterate each triangular face (2-simplex) to determine whethera trajectory intersects the face by solving the inverse linearinterpolation problem: u u u v v v µ µ µ = , (5) EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 8
Algorithm 1
Two-pass algorithm of tracking 0-, 1-, and 2- features—critical points, quantum vortices, and isosurfaces,respectively—with high-dimensional simplicial mesh. S is the set of simplices that test positive, UF being union-find. S ← ∅ , UF ← ∅ for each tri ∈ simplices(2) doif test(tri) then ⊲ Eq.5 S ← S ∪ tri end ifend forfor each tet ∈ simplices(3) do T ← S ∩ tet.sides(2)UF.unite( T ) end for S ← ∅ , UF ← ∅ for each tri ∈ simplices(2) doif test(tri) then ⊲ Eq.1 S ← S ∪ tri end ifend forfor each penta ∈ simplices(5) do T ← S ∩ penta.sides(2)UF.unite( T ) end for S ← ∅ , UF ← ∅ for each edge ∈ simplices(1) doif test(edge) then ⊲ Eq.11 S ← S ∪ edge end ifend forfor each penta ∈ simplices(5) do T ← S ∩ penta.sides(1)UF.unite( T ) end for m a x m a x s a d d l e m a x m a x s a d d l e m a x m a x s a d d l e m a x (a) (b) (c) s i n k s i n k s i n k saddle t i m e Fig. 8. Example of (a) filtering, (b) smoothing, and (c) simplification ofcritical point trajectories. Colors indicate different critical point types. where ( µ , µ , µ ) ⊺ are the normalized barycentric coor-dinates of the trajectory intersection with the face, and ( u , v ) ⊺ , ( u , v ) ⊺ , and ( u , v ) ⊺ are the vector values onthe three vertices of the 2-simplex. If µ j ∈ [0 , for all j ∈ { , , } , then the triangular face is punctured by atrajectory, and the spacetime coordinates of the critical point x c can be interpolated as x c = x c y c t c = x x x y y y t t t µ µ µ , (6)where ( x c , y c , t c ) ⊺ are the spacetime coordinates of thecritical point and ( x j , y j , t j ) ⊺ are the spacetime coordinatesof each vertex ( j ∈ { , , } ). In the second pass, we iterateover each tetrahedral cell (3-simplex) to associate its sidesthat are punctured by trajectories, because one can provethat each tetrahedron intersects up to one trajectory. Com-plete trajectories can be constructed by pairing every twopunctured triangular faces of the same tetrahedron.In general, for arbitrary dimensionality n , the two-pass algorithm to reconstruct critical point trajectories in v : R n +1 → R n is following. The first pass iterates overeach n -simplex to determine whether the simplex intersectsa trajectory based on Equation 5. The second pass iterateseach ( n + 1) -simplex and pairs its sides ( n -simplices) thatintersect a trajectory. The two-pass algorithm can be easilyparallelized with both distributed and GPU parallelism, asdiscussed in the following sections.The output reconstructed critical point trajectories areclosed curves in spacetime; they either end on domainboundaries or are loops. Within each curve, the critical pointtype may alternate, and there may be multiple monotoneintervals with respect to time. For example, as illustrated in Figure 8(a), each loop characterizes a maximum-saddle pairin the gradient field; the maximum-saddle pair establishesand annihilates simultaneously. In Figure 8(c), we see thebirth of a saddle-sink pair; the saddle further merges withanother sink soon after the birth. One can further simplifyand filter trajectories based on their attributes, as discussedin the following sections. The above two-pass algorithm assumes that PL vector fieldsare generic, an assumption that often does not hold for real-world data. For example, gradients of an integer-valuedimage may be exactly zero at vertices; gradients based oncentral-differences, which are rational numbers, may beaffinely dependent, causing nongeneric situations. In fluidflows, nonslip conditions lead to zero velocities on bound-aries. With nongeneric cases, however, there is no guaranteethat each ( n + 1) -simplex has up to one pair of intersectedsides. A critical point may reside on the boundary of the cell,causing numerical instabilities during the test of µ k ∈ [0 , in Equation 5. As a result, the critical point may or maynot be detected in the current n -simplex; chances are thatthe same critical point be detected by neighboring cells,causing nonrobust and noncombinatorial tracking results,as illustrated in Figure 1.We use the concept of simulation of simplicity [16] (SoS)to compute critical point trajectories robustly and combi-natorially in nongeneric vector fields. As proved by Bhatiaet al. [15] with the Brouwer degree theory, a critical pointexists in the n -simplex { x , x , . . . , x n } if and only if liesinterior of the convex hull of { v , v , . . . , v n } , where v j ( j = 0 , , . . . , n ) is the vector value on each vertex x i . Asa result, the critical point test is reduced to the point-in-simplex predicate, which can be determined by the sign ofdeterminant of the matrix ( v , v , . . . , v n , ) . In nongenericcases, such as a critical point on the boundary of the simplex,the SoS prevents the determinant from becoming zero byadding a symbolic perturbation such that the critical point isassociated with one of the simplices in a consistent manner. We provide three postprocessing approaches to help usersfilter, smooth, and simplify trajectories that result fromtracking critical points.
Filtering.
One can filter results based on theattributes—time duration, topology, persistence, scalar
EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 9 (a) (b)x yz x yz
Fig. 9. Nonrobust (a) versus robust (b) tracking of critical points in 20instances of the 3D moving minimum synthetic data, each with the sameorigin but random rational directions. Color encodes time in the figure. value, if applicable—of trajectories. Figure 8(a) is an exam-ple of filtering loops in the gradient of a scalar field. Typ-ically, a loop exists when a small transient bump appears,introducing a saddle-extremum pair. Such a loop may befiltered out based on the time duration or persistence of theloop. Likewise, trajectories can be filtered based on otherattributes, as demonstrated in the following sections.
Simplification.
One can also simplify trajectories thatchange directions frequently in time based on a threshold ofpersistence in time. As illustrated in Figure 8(c), for example,a saddle-sink pair is born right before the saddle mergeswith another sink. Because the saddle may be caused bynoise, we provide the simplification function to eliminatethe saddle and merge the trajectory as a consistent sink type.
Smoothing.
Although our trajectory reconstruction al-gorithm is robust, the evaluation of Jacobians J v may sub-ject to numerical instabilities, causing inconsistent criticalpoint types along trajectories. One can smooth the typesbased on a window-shifting approach. We iterate each pointin the trajectory and check whether the current critical pointtype is consistent with both the precedent and antecedent. Ifan inconsistency is identified, we mark the inconsistency lo-cation and modify the type after the iterations, as illustratedin Figure 8(b). The window size depends on the applicationand analysis needs. For example, we set the half-windowsize to be two consecutive timesteps in our experiments. We validate the effectiveness and evaluate the robustness ofour critical point tracking approach with synthetic data.
3D moving minimum.
We use the following scalarfunction to synthesize 3D time-varying scalar field data withthe known position of the single minimum to evaluate thenumeric robustnesss of our method: f ( x , t ) = k x − ( x + d · t ) k , (7)where ( x , t ) ∈ R n +1 are the spatiotemporal coordinates,and x and d are arbitrary vectors in R n . As a result, thetrajectory of the single minimum x c ( t ) in the data is x c ( t ) = x + d · t. (8)In Figure 9, we synthesized 20 different instances with thesame x = (10 , , ⊺ but with different d . The scalarfunction f is discretized into a × × grid, which (a) (b) σ=0, ground truth σ=0.08, unfilteredσ=0.02, unfiltered σ=0.08, filtered Fig. 11. Critical point tracking and filtering in the 2D synthetic spiralwoven data. (a) left: no noise injection, right: σ = 0.02 without filtering,(b) left: σ = 0.08 without filtering, right: σ = 0.08 with filtering. represents a (0 , , × (20 , , domain. Each componentof the moving direction d is a random rational numbersuch that the trajectory must intersect at least one gridpoint including x , causing degenerate cases similar to thatof Figure 1. In Figure 9(a), because the grid point of x is shared by multiple pentachora in the spacetime mesh,the number of tetrahedra that numerically test positive forcontaining a critical point ranges from zero to tens. Becausethe degenerate cases cause ambiguity in tracing, trajectoriesin the figure appear dashed and isolated. In Figure 9(b),our robust detection approach guarantees that each criticalpoint is exclusively associated with a tetrahedron such thattrajectories can be robustly tracked without any ambiguity. t i m e Fig. 10. Criticalpoint tracking in 2Dunstructured meshdouble gyre flow.
Color encodes ID oftrajectories.
2D double gyre flow.
Wedemonstrate critical point trackingin a 2D unstructured mesh withthe double gyre function, which iswidely used to study Lagrangiancoherent structures [48]. The dou-ble gyre function is defined in thedomain of [0 , × [2 , , and wegenerate a triangular mesh with2,098 triangles and 1,100 vertices todemonstrate the results. Figure 10visualizes critical trajectories in t ∈ [0 , ; the timespan ∆ t betweenadjacent timesteps is 0.1. The colorof each trajectory is categorical andencodes the unique ID of the tra-jectory. Slices at t = 0 and t = 2 . visualize flow directions with lineintegral convolution (LIC), and theslice at t = 1 . visualizes the 2D mesh. As shown in thefigure, two critical points move back and forth along the x -axis periodically, and the saddle points repeatedly enterand exit the 2D domain. Because the double gyre functionis analytical, we verified that the vector field is exactly zeroat all points in the trajectories, and all critical points areidentified by our method.
2D spiral woven with perturbations.
We design thespiral woven function f to evaluate our critical point track-ing method with the presence of noise: f ( x, y ) = cos( x cos t − y sin t ) sin( x sin t + y cos t ) . (9)This function has a finite number of critical points includingminima, maxima, and saddles for a bounded 2D domain.The critical points rotate around the origin at a fixed angular EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 10 speed over time. We discretized the function into a gridand injected Gaussian noise with two different standarddeviations to the data, as shown in Figure 11. Note thatthe range of the data is [ − , , and thus perturbation of σ = 0 . and σ = 0 . introduces up to 3% and 12% relativeerror, respectively, to the data in the 3-sigma limits. Becausethe noise injection produces many artificial bumps in thedata, the output critical point trajectories contain artifacts.Artifacts such as small loops can be removed through tra-jectory filtering, as discussed in the previous subsection. We demonstrate use cases of our critical point trackingmethods in science applications.
Fluid dynamics.
We track vector field critical points tovisualize the vortex streets in a 2D flow-past-cylinder sim-ulation (Figure 12(a)). The data are available on a uniformgrid with the resolution of × for 1,001 timesteps. Weuse the average velocity as the frame of reference and trackall critical over time. Based on the trajectories, we can seethat critical points are created in pairs past the cylinder andmove evenly toward the other side of the domain. Tokamak simulations.
We use trajectories of scalarfield critical points to characterize the dynamics of blobsin Tokamak fusion plasma simulations (Figure 12(b)). Un-derstanding blobs—regions of high turbulence that coulddamage the Tokamak—is critical for reactor design andfuture power generation. Assuming each blob correspondsto a maximum/minimum in the preconditioned and nor-malized electron density field, we reconstruct trajectoriesof critical points in a single poloidal plane consisting of56,980 vertices and 112,655 triangles in an XGC particle-in-cell simulation [50]. These trajectories can enable furtheranalyses of the physical properties of blobs.
Turbulent vortices.
We track and visualize scalar fieldcritical points in the publicly available turbulent vortexdataset [51], which contains 100 timesteps of scalarvalue data characterizing vorticity magnitudes. We visualizeresults for the first 30 timesteps (Figure 12(c)) to avoid over-plotting. There are 549 maximal trajectories out of 85,911trajectories in all types (minima, saddles, and maxima). Inthe visualization, one can see how high-turbulent locationsmove over space and time. Additional applications.
Our method has been ap-plied to other visualization techniques and science appli-cations. For visualization techniques, our method enablesfeature-preserving vector field compression by adaptingsufficient local error bounds for lossy compressors to pre-serve critical points in decompressed data without any falsepositives, false negatives, and false types [22]. For scienceapplications, we also demonstrated the use of our toolsto track microparticle clouds in laboratory plasma experi-ments [52].
First, in the case of critical point tracking in scalar fields, onehas to use the derived gradient fields as the input. However,the spacetime PL gradient field implies C continuity of theoriginal scalar field. The outputs may be distorted because of the smoothness of differentiation kernels used for gradi-ent derivation. That being said, the output trajectories are asaccurate as the quality of the gradient field and how closethe scalar field is to C continuity. In addition, the fidelityof critical point trajectories is related to the spacetime res-olution of the inputs. Downsampling the data may lead toan oversimplified topology of the trajectories. We leave thestudy of topology change due to downsampling to futurework. t i m e m a x m a x s a dd l e Fig. 13. Concurrentmerges and splits ofthree critical points.Color encodes tra-jectory ID.
Second, the determination of criticalpoint types is subject to numerical in-stabilities. Because Jacobians are usuallyestimated numerically, in extreme casesthe signs of eigenvalues may changewhen perturbation is introduced. Wetherefore introduced a smooth filter toheuristically correct classification errors,but one can also remediate the problembased on domain-specific knowledge.Third, the PL assumption preventsthe native identification of higher-ordercritical points [21, 23]. An example ofhigher-order cases is illustrated in Fig-ure 13(b), where two local maxima peri-odically merge and split. However, theresult includes a persisting trajectorythat characterizes one of the maximaand a number of maxima-saddle loopsperiodically. At the time of merging,ideally three critical points (two maximaand one saddle) should merge and thensplit into another three critical points.Because only up to one trajectory is allowed to intersect eachcell, our method cannot identify the “3-in-3-out” event overtime. We leave the study of higher-order dynamics of criticalpoints to future work as well.
Our approach offers three improvements over existing ap-proaches that use spacetime meshing [10, 11]. First, ourmethod uses simplicial meshes instead of prism meshes in4D spacetime. Because each of our simplicial cells intersectsup to one singularity, our method consistently avoids ambi-guity when multiple trajectories intersect a prism. Second,the PL assumption makes it easier to localize zero-crossingsin simplicial cells than in prism cells, a task that involvesiterative numerical algorithms such as Newton’s method.Third, our simplicial mesh enables the robust tracking ofcritical points based on SoS, producing combinatorial andconsistent results when nongeneric cases occur.Compared with approaches based on feature flow fields(FFFs) [31–34], our method is numerically robust and com-putationally scalable. FFFs are vector fields that characterizethe movement of critical points such that the trajectoriescan be computed as the streamlines in FFFs. First, in com-parison with FFF, which requires numerical approximationsof the gradients of vector fields, our method does notrequire such transformation. Numerical integration in FFF,such as Runge–Kutta methods, further introduces error instreamline tracing. Second, our two-pass algorithm can be
EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 11 T i m e (a) (b) (c) Fig. 12. Case studies of critical point tracking. (a) 2D flow-past-cylinder vector field (data courtesy of Weinkauf and Theisel [49]); color encodes theID of each trajectory. The two 2D slices at t = 60 and t = 200 each visualize the velocity direction and magnitude with line integral convolution imageand pseudo colors, respectively. (b) XGC fusion simulation (data courtesy of Dominski et al.); color of trajectories encodes whether a critical pointis a maximum (red) or minimum (blue). (c) Turbulent vortices (data courtesy of Ma); color encodes time (0 ≤ t ≤
30) of all maximal trajectories inthe large image; color encodes critical point types (blue for minima, white for saddles, and red for maxima) in the small image. (b)(a)
3D Case 4D Case I 4D Case II 4D Case III2D Case (c) θ0 mod(θ1-θ0, 2π)+mod(θ2-θ1, 2π)+mod(θ0-θ2, 2π)=2π θ1θ2 Fig. 14. Case(s) of a quantum vortex intersecting a 2-, 3-, and 4-simplices, respectively. embarrassingly parallelized, while computing streamlinesin distributed and parallel environments in general can bedifficult to scale [53, 54].
RACKING FEATURES : QUANTUM VORTICES
We demonstrate the use of high-dimensional simplicialmeshing to track 1D topological defects (quantum vortices)in 3D complex-valued scalar fields produced by supercon-ductivity, superfluidity, and Bose–Einstein condensate sim-ulation data. Unlike critical points in vector fields, vorticesare 1D curves in individual timesteps, and the trajectories ofvortices are 2D surfaces embedded in 4D spacetime.
Vortex tracking can be achieved by a two-pass reconstruc-tion algorithm, as outlined in Algorithm 1 (middle). Wereview here how vortices intersect 2- and 3-simplices andthen introduce the two-pass vortex surface reconstruction.
Contour integral test on triangular faces.
We testwhether each spacetime 2-simplex intersects a vortex basedon the definition in Equation 1. As shown in Figure 14(a), atriangular face intersects a vortex at most once. The contourintegral can be calculated by accumulating phase shiftsalong each edge: I C ∇ θ ( x ) · d l = X i =0 ∆ θ i,j , i = 0 , , , j = i + 1 mod , (10) where the phase shift ∆ θ i,j is the modulo of θ j − θ i and π ; θ i and θ j , respectively, are the phase angle at the i thand j th vertex of the triangular face. The test is positive ifthe contour integral equals ± π . Because the modulo is lessthan π , we assume that the data resolution in both spatialand temporal dimensions is sufficiently fine that the phasedifference between each pair of adjacent vertices is lessthan π . Scientists would need to refine the spatiotemporalresolution of simulations if the discretization is too coarse. Vortex curve reconstruction.
One can reconstruct the1D topology of vortices in the R subspace of individualtimesteps. Because vortices are closed curves In R , eachclosed volume has an even number of intersections onboundaries; each tetrahedral cell has up to one pair ofintersected faces [13], as illustrated in Figure 14(b). As aresult, vortex curves can be reconstructed by associatingintersected faces that share the same tetrahedral cells. Theprocess involves two passes: first iterating all triangularfaces and then scanning tetrahedral cells that have inter-sected faces. Vortex surface reconstruction.
The two-pass curve re-construction can be generalized to 4D spacetime to charac-terize the trajectory of vortex curves as surfaces. In 4D, eachpentachoron has five tetrahedral faces; each tetrahedronmay intersect a vortex. Considering that tetrahedron sharestriangular sides in the pentachoron, the number of triangu-lar sides that test positive may be 3, 4, or 5 (correspondingto Case I, II, or III in Figure 14). Each intersection sits on thesame 2-manifold of the vortex surface. In the reconstructionof the vortex surface, the first pass tests all triangular faces in4D spacetime, and then the second pass joins all triangularfaces that test positive and share the same pentachoron.
Benefits.
The use of spacetime simplicial meshing sim-plifies and improves on previous implementations basedon cubic meshes [12] and triangular prism meshes [13].In the former study, ambiguity exists when two vorticespenetrate the same cube. In the latter study, to eliminateambiguities with cubic cells, each spatial cube is tessellatedinto six tetrahedra such that each tetrahedron intersects upto one vortex line. Although mesh elements in individual
EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 12 timesteps are simplices, simplicial prisms are used in space-time, causing two problems: (1) two different functions areneeded to test triangular and quadrilateral faces of a prism,and (2) ambiguity is still possible with prism cells. Withsimplicial spacetime meshes, we need only one function totest triangular faces, and there is no ambiguity.
Figure 15 demonstrates vortex tracking results of a time-dependent Ginzburg-Landau superconductivity simulation,produced by the finite-difference partial differential equa-tion solver GLGPU [55]. In this case, the resolution of themesh is × × , and the number of the timestepsis 200. As a result, 6.5M out of 251M spacetime triangularfaces tested positive for encircling a vortex, forming a set ofdisjoint surfaces.From the scientific perspective, the dynamics of vorticesdetermine all electromagnetic properties of the supercon-ducting materials, and recombinations of vortices are di-rectly related to energy dissipation [56]. The surface-basedvisualization produced by our tools enables scientists to in-vestigate the time-varying features with a single image. Thereconstructed surface also makes it possible to derive themoving speed of each vortex when there is no topologicalchanges; the moving speed has positive correlation with thevoltage drop, which is critical to the material design. RACKING FEATURES : ISOSURFACES
Isosurface tracking in FTK is also based on high-dimensional simplicial meshes. As demonstrated by Bhani-ramka et al. [8] and Ji et al. [9], isosurfaces can be trackedin 3D scalar fields by extracting and slicing levelsets in R regular meshes. While the major complexities of previousefforts—disambiguation of isosurfaces intersecting the samehypercube—may be resolved by triangulation, our imple-mentation with high-dimensional simplicial meshing intrin-sically avoids ambiguities in a consistent manner and scalesto larger computing resources with the FTK framework. We formalize the isosurface tracking problem as the re-construction of isovolumes in the time-varying scalar field f : R n +1 → R , isovolumes being the solution of f = c ,where c is the isovalue. In general, the isovolume is an ( n − -dimensional object embedded in R n +1 . With n = 3 ,the object can be further sliced into 2D surfaces with fixedtime values. We assume that f is generic; that is, the scalarvalue f i = c for each vertex i , and scalars at vertices ofany k -simplex ( k = 1 , , . . . , n + 1 ) are affinely independent.Thus, each edge (2-simplex) in the spacetime mesh intersectsthe level set at most once; the edge neither resides on thelevel set nor intersects the it at vertices. Similar to the robustcritical point tracking method in Section 4, we use symbolicperturbations [16] to relax the generic assumption and toreconstruct 4D level sets in a robust and consistent mannerfor real application data.As shown in Algorithm 1, the reconstruction consists oftwo passes for isovolume reconstruction in 4D: the edge pass and pentachoron pass. In the first pass, we checkwhether every edge (2-simplex) intersects a level set bysolving the inverse interpolation: µf + (1 − µ ) f = c, (11)where f and f are the scalar values on each vertex ofthe simplex, ( µ, − µ ) ⊺ being the barycentric coordinatesof the intersection point. The edge intersects a level set ifand only if µ ∈ (0 , . In the second pass, we iterate everypentachoron; edges that intersect the level set are associatedand labeled with the same identifier.In the case of a 3D time-varying scalar field, the outputisovolumes are 3D objects and can be represented as atetrahedral grid, which can be further sliced into isosurfacesin individual timesteps for visualization and analysis. Asillustrated in Figure 16(a), the intersection between a pen-tachoron and an isovolume has only two distinct cases.We use the positive and negative sign, respectively, torepresent whether the scalar value on a vertex is greateror less than the threshold. In the 4D case I ( + − − − − or other permutations), one of the vertices has a differentsign than other vertices do; in this case, the isovolumeinside the pentachoron is the single tetrahedron consistingof the four intersections. In the 4D case II ( + + − − − orother permutations), the pentachoron has six intersections,which form a polytope and can be triangulated into threetetrahedra, as explained below.Without loss of generality, we show that the second case + + − − − leads to an 3-polytope that can be tessellated intothree tetrahedra. In the five tetrahedral sides of the penta-choron, we find three tetrahedra with + + −− and the othertwo tetrahedra with + − −− . As illustrated in Figure 16(b),in the 3D case I ( + −−− ), the intersection is a triangle; in the3D case II ( + + −− ), the intersection of the tetrahedron andthe isovolume leads to a coplanar quadrilateral. As a result,we have two triangles and three quadrilaterals, forming aprism-like polytope that resides in the same 3D subspace.With the same staircase triangulation method described inSection 3.2, we can decompose the polytope into three non-overlapping tetrahedra in a combinatorial manner. The robustness of this method is guaranteed by the assump-tion that the function is generic; that is, the scalar valueof each vertex is either greater or less than the isovalue,and Equation 11 has one solution. For real problems thatusually do not observe this assumption, degenerate casesmay appear. For example, should the scalar value of avertex exactly equal the isovalue, the intersection may ormay not be identified by other edges that share the samevertex, leading to unstable results. Should an edge resideon an isosurface, every point on the edge is a solution ofEquation 11, leading to numerical instabilities.In nongeneric cases, we ensure the robustness with thesame SoS programming technique [16] used for robust crit-ical point tracking (described in Section 4.3), and we regardthe test in Equation 11 as a special case of critical pointdetection in a 1D vector field. Thus, if an intersection is onany vertex of the edge, the SoS prevents the divisor frombecoming zero by adding a symbolic perturbation, such that
EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 13 loop annihilation in early stagesrecombination t=132 t=141 0 200time
Fig. 15. Trajectory surfaces of quantum vortices in a 3D superconductivity simulation data with 200 timesteps; color encodes the timestep. Vorticesat timestep 132 are visualized with tubes. The two zoomed subfigures visualize how a pair of vortices exchanges parts before and after therecombination event. - -+ - - ++- + +- + + +- -+- (a)
4D Case I 4D Case II 3D Case I 3D Case II(b)
Fig. 16. Cases of an isovolume intersecting a pentachoron (a) and atetrahedron (b). Signs indicate whether the scalar value on a vertex isgreater or less than the given isovalue. (a) (b) t=1295t=1296t=1297t=1298t=10t=0
Fig. 17. Tracking isosurfaces: (a) projected isovolume (top) and isosur- faces (bottom) of synthetic data and (b) isosurfaces of supernova data.
Colors in (a) and (b) encode time and surface ID, respectively. the intersection is exclusively associated with one of theedges that share the same vertex.
We demonstrate our isovolume reconstruction with bothsynthetic and simulation data in Figure 17 and below.
Synthetic data.
We show in Figure 17(a) the reconstruc-tion of isovolume of f = 0 for the function f ( x, y, z, t ) = x − αt on a × × grid and 12 timesteps, with α = 0 . .Since the closed form of the surface trajectory ( x = αt ) isavailable, we can verify that the reconstructed isovolumeand the isosurfaces sliced from the isovolume are correct. Supernova simulation data.
We reconstruct the iso-volume of a supernova simulation dataset [57] with theisovalue of 0.8 and visualize the sliced isosurfaces of four timesteps in Figure 17(b). The resolution of the data is ,which in turn produces ≈ G spacetime edges to test inter-sections for the 4-timestep data. As a result, the isovolumeconsists of 25M intersection points and 141M tetrahedra,which can be sliced back into individual timesteps forvisualization and analysis.
Compared with existing methods of isosurfacing in higherdimensions [8, 9], the benefits of our method include(1) straightforward implementation, (2) no disambiguationcases, (3) guaranteed robustness, and (4) straightforwardparallelization with GPUs and distributed environments.First, the two-pass algorithm can be written in a few lines ofcode based on FTK’s meshing APIs; there is no need to gen-erate large lookup tables for higher-dimensional marchingcubes. Second, compared with marching cubes [7], there areno ambiguity cases with simplicial meshes. Our method canbe viewed as a generalization of marching tetrahedra [14]in higher dimensions, which automatically eliminates anyambiguities. Third, the use of the simplicial mesh makes itpossible to ensure robustness with symbolic perturbations,leading to consistent tracking results. Fourth, our two-passalgorithm can be easily distributed and computed withparallel computing resources.
LIBRARY DESIGN
We describe FTK’s meshing, numeric, union-find, and I/Ocomponents, the design of which takes into account paral-lelization for both distributed- and shared-memory environ-ments, and the needs of both in situ and post hoc analyses.We also introduce developer-level APIs provided in FTK forimplementing feature-tracking algorithms.
FTK implements the following C++11 APIs to support thedevelopment of feature-tracking algorithms that traversemesh elements along different dimensions of ( n + 1) -Dspacetime simplical meshes. • element_for(int k,function
Union-find is a key algorithm for partitioning a set ofmesh elements into disjoint subsets for feature tracking. Weimplement an asynchronous distributed union-find methodto enable distributed feature tracking. Because the operationof merging sets is order-independent, we use an asyn-chronous message-passing approach to exchange interpro-cess requests; the final results remain unchanged even ifmessages are received in random order. By eliminating fre-quent and expensive global synchronizations, our methodoutperforms existing implementations based on the bulk-synchronous parallel programming model. Details of ourdistributed union-find algorithm can be found in [61]. from netCDF4 import Datasetfrom scipy.ndimage import gaussian_filterimport pyftk as ftkdataset = Dataset('Ni_pol.nc', 'r')ni = dataset.variables['NI_TRZ'] (a)(b)
Fig. 18. Example uses of FTK’s (a) ParaView plugins and (b) Pythonbindings.
Inputs.
FTK libraries can be built with various I/O li-braries, including VTK [1], NetCDF [62], HDF5 [63], andADIOS2 [64], to load data in a variety of formats. Datastreamed in situ can be loaded into memory with ADIOS2.In addition,
NumPy data are supported through PyFTK,described in Section 8.3.
Outputs.
Output formats include VTK, text, andPython objects. With VTK, trajectories and surfaces aretransformed into vtkPolyData ; isovolumes are writtenin vtkUnstructuredGrid formats. FTK also supportshuman-readable text formats. Python objects can be re-trieved and serialized into either JSON or pickle formats.
N SITU AND POST HOC ANALYSIS WITH
FTK
FTK provides four different utilities for end users:VTK/ParaView filters, a command line interface, Pythonbindings, and a C++ programming interface. TheVTK/ParaView filters are designed for interactive visual-ization, with the possibility to couple simulations throughParaView Catalyst and other in situ frameworks. The com-mand line interface can be used for loosely coupled in situprocessing and post hoc analyses. The Python bindingsare designed for post hoc analyses and integration withdata science libraries. The C++ programming interface isdesigned mainly for tightly coupled in situ analyses.
We developed VTK wrappers to use FTK interactively withParaView. The filters are designed for interactive visualiza-
EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 15 tion, but they also can be used with ParaView Catalyst forin situ processing.
Synthetic data sources.
We implemented syntheticdata generators for users to learn FTK filters. Synthetic datainclude double gyre, spiral woven, and merger functions(Section 4.5). Users can tune parameters to control spa-tiotemporal resolutions of the data.
Filters.
We wrapped FTK’s feature-tracking algorithmsinto VTK filters, such as vtkCriticalPointTracker and vtkLevelsetTracker . Currently, the inputs need to beimage data types, and the output vtkPolyData can bedirectly rendered and further processed with other filterswith ParaView. Figure 18(a) demonstrates a critical pointtracking case: a scalar field is generated from a syntheticdata source, and then critical point trajectories are extractedand transformed into tubes for visualization.
We provide executables for tracking critical points, quantumvortices, and isosurfaces in data obtained from files and insitu data streams. For file inputs, FTK supports multiple fileformats, including NetCDF, HDF5, VTK, ADIOS2, and rawbinaries. Users need to explicitly specify variable names forself-described formats. For in situ data streams, users mustspecify ADIOS2 stream sources, variable names, mesh infor-mation, and other needed parameters. The command lineinterface automatically loads and handles data in parallel ifexecuted with MPI. Users are also provided with options touse multithreading and GPU accelerators.
PyFTK, which is a Python binding library for FTK, en-ables post hoc feature tracking and analysis in Python.PyFTK functions take
NumPy arrays as inputs, allowingusers to load data with Python’s netCDF4 , h5py , andother I/O libraries. PyFTK also enables easy integrationwith other libraries such as SciPy and
Matplotlib fordata science workflows. For example, in Figure 18(b), weload BOUT++ [65] fusion plasma simulation data with theNetCDF4 Python module, apply Gaussian smoothing inspacetime with
SciPy , and track critical points with PyFTK.
The FTK’s C++ APIs enable users to couple simulation codeswith FTK and/or customize feature-tracking algorithms.The former can be achieved with user-level APIs; the latterrequires the use of developer-level APIs.
User-level APIs.
APIs are available to users to feedtime-varying data into FTK algorithms. Users can call push_field_data() functions in various tracker classes.Each timestep of the input field data must be convertedinto ftk::ndarray , a dynamic multidimensional arraythat eases the internal data transformation and access. Util-ity functions are provided to convert C-style arrays and vtkDataArray s into ftk::ndarray . Developer-level APIs.
Users can customize feature-tracking algorithms by using element_for() , sides() ,and side_of() functions for access to high-dimensionalmesh elements. One can also use linear algebra functions tohelp localize features on both CPUs and GPUs. T i m e ( s e c o n d s ) Num ber of CPU cores / GPUs I d e a l s c a l i n g s l o p e Critical point, 20 timesteps(256 , double prec.)Isosurface, 5 timesteps(432 , double prec.)Quantum vortex, 10 timesteps(128x512x64, single prec.) (CPU)(GPU) Fig. 19. Timings, in seconds, of critical point tracking, isovolume, andquantum vortices computations for different numbers of CPU cores andGPUs on the Summit supercomputer.
ERFORMANCE E VALUATION
We benchmarked the scalability of FTK on Summit, a 200-petaflop supercomputer in Oak Ridge National Labora-tory. The system consists of 4,608 IBM AC922 computingnodes, each of which has two 22-core Power9 CPUs andsix NVIDIA Tesla V100 GPUs. The clock frequency of eachCPU is 3.07 GHz, and each node is equipped with 1,600 GBof high-bandwidth memory. The interconnection betweennodes is 100Gbps EDR InfiniBand. As shown in Figure 19,we characterize the strong scalability by measuring theexecution time of solving three feature-tracking problemswith different numbers of processes, each of which useseither one GPU or one CPU core.
GPU acceleration.
Compared with the CPU perfor-mance, the GPU acceleration typically ranges from O (100) to O (1000) . The acceleration is expected because there is nosynchronization or communication between GPU threads;each thread independently tests a mesh element. The mag-nitude of GPU acceleration varies, possibly because of dif-ferent complexities of handling mesh elements, precisionof numeric algorithms (double vs. single), and data accesscost (e.g., accessing scalars, vectors, and Jacobians in criticalpoint tracking vs. accessing only scalars in isosurface track-ing). We will study performance variability in future work. Scalability.
CPUs appear to scale better than GPUs.As we continue increase the number of processes by eachfactor of 2, the acceleration of using GPUs saturates fasterthan when using CPUs. Because the iteration over sim-plices are accelerated by GPUs, the non-GPU cost (e.g., datamovement and connected component) dominates when theworkload per process decreases. In our applications such asfusion and superconductivity simulations, which typicallyproduce 3 to 4 timesteps per second [13, 66], our algorithmsare able to keep up with the data-producing rate in situ.We will further investigate the performance and resourceconsumption during in situ processing in the future.
10 C
ONCLUSIONS AND F UTURE W ORK
We have demonstrated the use of spacetime simplicialmeshes for simplifying, scaling, and delivering feature-
EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 16 tracking algorithms. Our contributions include high-dimensional simplicial meshing and feature-tracking algo-rithms over these meshes. We have presented methods fortessellating ( n + 1) -D simplicial prisms and regular meshesin arbitrary dimensions. We track three types of features:critical points (0D features embedded in 2D/3D), quantumvortex curves (1D features embedded in 3D), and isosur-faces (2D features embedded in 3D). We have implementedthese contributions in a suite of feature-tracking tools calledFTK. FTK provides VTK/ParaView plugins, a commandline interface, Python bindings, and C++ programming in-terfaces for both in situ and post hoc analysis and visualiza-tion. Performance studies show that FTK algorithms can beaccelerated by both GPUs and distributed parallelism, withhigh scalability.In the future we plan to investigate four aspects: sup-porting more types of features, adapting to time-varyingmeshes, incorporating scale spaces, and porting to newprocessor and memory technologies. First, we will expandFTK’s support of features including parallel vectors (1Dfeatures in 3D), ridge/valley surfaces (2D features in 3D),and interval volumes (3D features in 3D). Second, we willdevelop spacetime mesh generalizations for time-varyingmeshes such as finite volume and adaptive mesh refinementmeshes for broader scientific applications. Third, we willincorporate scales as the fifth dimension for scale-spacetracking with our high-dimensional meshes. Fourth, we willport FTK to emerging hardware including GPUs, FPGAs,and high-bandwidth memory in order to enable acceleratedcomputing and data handling. A CKNOWLEDGMENTS
This research is supported by the Exascale ComputingProject (ECP), project number 17-SC-20-SC, a collaborativeeffort of Department of Energy Office of Science and theNational Nuclear Security Administration, as part of theCo-design center for Online Data Analysis and Reduction(CODAR) [66]. It is also supported by the U.S. Departmentof Energy, Office of Advanced Scientific Computing Re-search, Scientific Discovery through Advanced Computing(SciDAC) program, and by Laboratory Directed Researchand Development (LDRD) funding from Argonne NationalLaboratory, provided by the Director, Office of Science, ofthe U.S. Department of Energy under Contract No. DE-AC02-06CH11357. R EFERENCES [1] W. Schroeder, K. Martin, and B. Lorensen,
The Visualization Toolkit(4th ed.) . Kitware, 2006.[2] K. Moreland, C. M. Sewell, W. Usher, L. Lo, J. S. Meredith,D. Pugmire, J. Kress, H. A. Schroots, K. Ma, H. Childs, M. Larsen,C. Chen, R. Maynard, and B. Geveci, “VTK-m: Accelerating thevisualization toolkit for massively threaded architectures,”
IEEEComputer Graphics and Applications , vol. 36, no. 3, pp. 48–58, 2016. [3] U. Ayachit,
The ParaView Guide: A Parallel Visualization Application .Kitware, 2015.[4] H. Childs, E. Brugger, B. Whitlock, J. Meredith, S. Ahern, D. Pug-mire, K. Biagas, M. Miller, C. Harrison, G. H. Weber, H. Krishnan,T. Fogal, A. Sanderson, C. Garth, E. W. Bethel, D. Camp, O. R ¨ubel,M. Durant, J. M. Favre, and P. Navr´atil, “VisIt: An end-user tool for visualizing and analyzing very large data,” in
High PerformanceVisualization—Enabling Extreme-Scale Scientific Insight , E. W. Bethel,H. Childs, and C. Hansen, Eds., 2012, pp. 357–372. [5] J. Tierny, G. Favelier, J. A. Levine, C. Gueunet, and M. Michaux,“The Topology ToolKit,”
IEEE Transactions on Visualization andComputer Graphics , vol. 24, no. 1, pp. 832–842, 2018.[6] G. Bradski, “The OpenCV Library,”
Dr. Dobb’s Journal of SoftwareTools , 2000.[7] W. E. Lorensen and H. E. Cline, “Marching cubes: A high resolu-tion 3D surface construction algorithm,” in
Proc. SIGGRAPH ’87 ,1987, pp. 163–169.[8] P. Bhaniramka, R. Wenger, and R. Crawfis, “Isosurfacing in higherdimensions,” in
Proc. IEEE Visualization ’00 , 2000, pp. 267–273.[9] G. Ji, H. Shen, and R. Wenger, “Volume tracking using higherdimensional isocontouring,” in
Proc. IEEE Visualization ’03 , 2003,pp. 209–216.[10] X. Tricoche, T. Wischgoll, G. Scheuermann, and H. Hagen,“Topology tracking for the visualization of time-dependent two-dimensional flows,”
Computers & Graphics , vol. 26, no. 2, pp. 249–257, 2002.[11] C. Garth, X. Tricoche, and G. Scheuermann, “Tracking of vectorfield singularities in unstructured 3D time-dependent datasets,”in
Proc. IEEE Visualization ’04 , 2004, pp. 329–336.[12] H. Guo, C. L. Phillps, T. Peterka, D. Karpeyev, and A. Glatz,“Extracting, tracking and visualizing magnetic flux vortices in 3Dcomplex-valued superconductor simulation data,”
IEEE Transac-tions on Visualization and Computer Graphics , vol. 22, no. 1, pp. 827–836, 2016.[13] H. Guo, T. Peterka, and A. Glatz, “In situ magnetic flux vortex vi-sualization in time-dependent Ginzburg-Landau superconductorsimulations,” in
Proc. IEEE Pacific Visualization Symposium , 2017,pp. 71–80.[14] A. Doi and A. Koide, “An efficient method of triangulating equi-valued surfaces by using tetrahedral cells,”
IEICE Transactions onInformation and Systems , vol. E74-D, no. 1, pp. 214–224, 1991.[15] H. Bhatia, A. Gyulassy, H. Wang, P. Bremer, and V. Pascucci,“Robust detection of singularities in vector fields,” in
TopologicalMethods in Data Analysis and Visualization III, Theory, Algorithms,and Applications , P. Bremer, I. Hotz, V. Pascucci, and R. Peikert,Eds. Springer, 2014, pp. 3–18.[16] H. Edelsbrunner and E. P. M ¨ucke, “Simulation of simplicity: atechnique to cope with degenerate cases in geometric algorithms,”
ACM Transactions on Graphics , vol. 9, no. 1, pp. 66–104, 1990.[17] F. H. Post, B. Vrolijk, H. Hauser, R. S. Laramee, and H. Doleisch,“The state of the art in flow visualisation: Feature extraction andtracking,”
Computer Graphics Forum , vol. 22, no. 4, pp. 775–792,2003.[18] C. Heine, H. Leitte, M. Hlawitschka, F. Iuricich, L. D. Floriani,G. Scheuermann, H. Hagen, and C. Garth, “A survey of topology-based methods in visualization,”
Computer Graphics Forum , vol. 35,no. 3, pp. 643–667, 2016.[19] A. L. Haynes and C. E. Parnell, “A trilinear method for findingnull points in a three-dimensional vector space,”
AIP Physics ofPlasmas , vol. 14, no. 082107, pp. 1–8, 2007.[20] J. Helman and L. Hesselink, “Representation and display of vectorfield topology in fluid flow data sets,”
IEEE Computer , vol. 22,no. 8, pp. 27–36, 1989.[21] X. Tricoche, G. Scheuermann, and H. Hagen, “A topology simpli-fication method for 2D vector fields,” in
Proc. IEEE Visualization2000 , 2000, pp. 359–366.[22] X. Liang, H. Guo, S. Di, F. Cappello, M. Raj, C. Liu, K. Ono,Z. Chen, and T. Peterka, “Toward feature-preserving 2D and3D vector field compression,” in
Proc. IEEE Pacific VisualizationSymposium 2020 , 2020, pp. 81–90.[23] T. Weinkauf, H. Theisel, K. Shi, H. Hege, and H. Seidel, “Extractinghigher order critical points and topological simplification of 3Dvector fields,” in
Proc. IEEE Visualization ’05 , 2005, pp. 559–566.[24] K. Polthier and E. Preuß, “Identifying vector field singularitiesusing a discrete Hodge decomposition,” in , ser. Mathematics and Visualization, H. Hege and
K. Polthier, Eds. Springer, 2002, pp. 113–134.[25] J. M.Greene, “Locating three-dimensional roots by a bisectionmethod,”
Journal of Computational Physics , vol. 98, no. 2, pp. 194–198, 1992.[26] S. Mann and A. P. Rockwood, “Computing singularities of 3Dvector fields with geometric algebra,” in
Proc. IEEE Visualization’02 , 2002, pp. 283–289. [27] W. Li, B. Vallet, N. Ray, and B. L´evy, “Representing higher-ordersingularities in vector fields on piecewise linear surfaces,”
IEEE
EEE TRANSACTIONS OF VISUALIZATION AND COMPUTER GRAPHICS, VOL. X, NO. X, NOVEMBER 2020 17
Transactions on Visualization and Computer Graphics , vol. 12, no. 5,pp. 1315–1322, 2006.[28] G. Chen, K. Mischaikow, R. S. Laramee, and E. Zhang, “EfficientMorse decompositions of vector fields,”
IEEE Transactions on Visu-alization and Computer Graphics , vol. 14, no. 4, pp. 848–862, 2008.[29] G. Reeb, “Sur les points singuliers d’une forme de Pfaffcompl`etement int´egrable ou d’une fonction num´erique,”
ComptesRendus Acad. Sciences , vol. 222, pp. 847–849, 1946.[30] H. Carr, J. Snoeyink, and U. Axen, “Computing contour trees inall dimensions,” in
Proc. 11th Annual ACM-SIAM Symposium onDiscrete Algorithms , 2000, pp. 918–926.[31] H. Theisel and H.-P. Seidel, “Feature flow fields,” in
Proc. Euro-Graphics/IEEE VGTC Symposium on Visualization ’03 , 2003, pp. 141–148.[32] T. Weinkauf, H. Theisel, A. V. Gelder, and A. T. Pang, “Stablefeature flow fields,”
IEEE Transactions on Visualization and ComputerGraphics , vol. 17, no. 6, pp. 770–780, 2011.[33] T. Klein and T. Ertl, “Scale-space tracking of critical points in 3Dvector fields,” in
Topology-based Methods in Visualization , ser. Math-ematics and Visualization, H. Hauser, H. Hagen, and H. Theisel,Eds. Springer, 2007, pp. 35–49.[34] J. Reininghaus, J. Kasten, T. Weinkauf, and I. Hotz, “Efficient com-putation of combinatorial feature flow fields,”
IEEE Transactions onVisualization and Computer Graphics , vol. 18, no. 9, pp. 1563–1573,2012.[35] B. Wang, P. Rosen, P. Skraba, H. Bhatia, and V. Pascucci, “Visualiz-ing robustness of critical points for 2D time-varying vector fields,”
Computer Graphics Forum , vol. 32, no. 3, pp. 221–230, 2013.[36] P. Skraba and B. Wang, “Interpreting feature tracking throughthe lens of robustness,” in
Topological Methods in Data Analysisand Visualization III, Theory, Algorithms, and Applications , P. Bremer,I. Hotz, V. Pascucci, and R. Peikert, Eds. Springer, 2014, pp. 19–37.[37] C. L. Phillips, T. Peterka, D. Karpeyev, and A. Glatz, “Detectingvortices in superconductors: Extracting one-dimensional topolog-ical singularities from a discretized complex scalar field,”
PhysicsReview E , vol. 91, no. 023311, pp. 1–12, 2015.[38] Y. Guo, X. Liu, C. Xiong, X. Xu, and C. Fu, “Towards high-quality visualization of superfluid vortices,”
IEEE Transactions onVisualization and Computer Graphics , vol. 24, no. 8, pp. 2440–2455,2018.[39] C. L. Phillips, H. Guo, T. Peterka, D. Karpeyev, and A. Glatz,“Tracking vortices in superconductors: Extracting singularitiesfrom a discretized complex scalar field evolving in time,”
PhysicalReview E , vol. 93, no. 023305, pp. 1–13, 2016.[40] M. Jiang, R. Machiraju, and D. Thompson, “Detection and visu-alization of vortices,” in
The Visualization Handbook , C. D. Hansenand C. R. Johnson, Eds., 2005.[41] J. Jeong and F. Hussain, “On the identification of a vortex,”
Journalof Fluid Mechanics , vol. 285, pp. 69–94, 1995.[42] N. J. Zabusky, O. N. Boratav, R. B. Pelz, M. Gao, D. Silver, andS. P. Cooper, “Emergence of coherent patterns of vortex stretchingduring reconnection: A scattering paradigm,”
Phys. Rev. Lett. ,vol. 67, no. 18, pp. 2469–2472, 1991.[43] D. Silver and X. Wang, “Tracking scalar features in unstructureddatasets,” in
Proc. IEEE Visualization ’98 , 1998, pp. 79–86.[44] ——, “Volume tracking,” in
Proc. IEEE Visualization ’96 , 1996, pp.157–164.[45] M. Ishii, M. Fernando, K. Saurabh, B. Khara, B. Ganapathysub-ramanian, and H. Sundar, “Solving PDEs in space-time: 4D tree-based adaptivity, mesh-free and matrix-free approaches,” in
Proc.SC ’19: International Conference for High Performance Computing,Networking, Storage, and Analysis , 2019, pp. 61:1–61:22.[46] S. Thite, “Adaptive spacetime meshing for discontinuous galerkinmethods,”
Computational Geometry , vol. 42, no. 1, pp. 20–44, 2009.[47] J. A. D. Loera, J. Rambau, and F. Santos,
Triangulations: Structuresfor Algorithms and Applications . Springer, 2010.[48] S. Shadden,
Time-dependent double gyre , ac-cessed November 14, 2020. [Online]. Available:https://shaddenlab.berkeley.edu/uploads/LCS-tutorial/examples.html[49] T. Weinkauf and H. Theisel, “Streak lines as tangent curves ofa derived vector field,”
IEEE Transactions on Visualization andComputer Graphics , vol. 16, no. 6, pp. 1225–1234, 2010. [50] C. S. Chang and S. Ku, “Spontaneous rotation sources in a quies- cent tokamak edge plasma,”
Physics of Plasmas , vol. 15, no. 062510,pp. 1–9, 2008. [51] K.-L. Ma,
TVDR: Time-Varying Volume Data Reposi-tory ∼ ma/ITR/tvdr.html[52] Z. Wang, J. Xu, Y. E. Kovach, B. T. Wolfe, E. Thomas, H. Guo, J. E.Foster, and H.-W. Shen, “Microparticle cloud imaging and trackingfor data-driven plasma science,” Physics of Plasmas , vol. 27, no. 3,p. 033703, 2020.[53] J. Zhang, H. Guo, F. Hong, X. Yuan, and T. Peterka, “Dynamicload balancing based on constrained K-D tree decomposition forparallel particle tracing,”
IEEE Transactions on Visualization andComputer Graphics , vol. 24, no. 1, pp. 954–963, 2018.[54] D. Pugmire, T. Peterka, and C. Garth, “Parallel integral curves,”in
High Performance Visualization: Enabling Extreme Scale ScientificInsight , E. W. Bethel, H. Childs, and C. Hansen, Eds. CRC Press,2012.[55] I. Sadovskyy, A. Koshelev, C. Phillips, D. Karpeyev, and A. Glatz,“Stable large-scale solver for Ginzburg-Landau equations for su-perconductors,”
Journal of Computational Physics , vol. 294, no. C,pp. 639–654, 2015.[56] A. Glatz, V. K. Vlasko-Vlasov, W. K. Kwok, and G. W. Crabtree,“Vortex cutting in superconductors,”
Physical Review B , vol. 94, no.064505, pp. 1–11, 2016.[57] V. Morozova, D. Radice, A. Burrows, and D. Vartanyan, “Thegravitational wave signal from core-collapse supernovae,”
TheAstrophysical Journal , vol. 861, no. 10, pp. 1–19, 2018.[58] G. Guennebaud and B. Jacob,
Eigen v3 , accessed November 14,2020. [Online]. Available: http://eigen.tuxfamily.org[59] C. Riccio,
OpenGL Mathematics (GLM) , accessed November 14,2020. [Online]. Available: https://github.com/g-truc/glm[60]
Boost C++ Libraries
IEEE Computer Graphics and Applications
IPDPS’09: Proc. IEEE International Symposium on Parallel andDistributed Processing , 2009, pp. 1–10.[65] B. Dudson, M. Umansky, X. Xu, P. Snyder, and H. R. Wilson,“BOUT++: A framework for parallel plasma fluid simulations,”
Computer Physics Communications , vol. 180, no. 9, pp. 1467–1480,2009.[66] I. Foster, M. Ainsworth, J. Bessac, F. Cappello, J. Choi, S. Di, A. M.Gok, H. Guo, K. A. Huck, C. Kelly, S. Klasky, K. Kleese van Dam,X. Liang, K. Mehta, M. Parashar, T. Peterka, L. Pouchard, T. Shu,H. van Dam, J. M. Wozniak, M. Wolf, W. Xu, I. Yakushin, S. Yoo,and T. Munson, “Online data analysis and reduction: An impor-tant co-design motif for extreme-scale computers,”
InternationalJournal of High-Performance Computing Applications , vol. in press,2020., vol. in press,2020.