{'fraction_non_alphanumeric': 0.026990696522692703, 'fraction_numerical': 0.006438147794403764, 'mean_word_length': 4.887026239067056, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Level Set Restricted Voronoi Tessellation for Large scale Spatial Statistical Analysis

Spatial statistical analysis of multivariate volumetric data can be challenging due to scale, complexity, and occlusion. Advances in topological segmentation, feature extraction, and statistical summarization have helped overcome the challenges. This work introduces a new spatial statistical decomposition method based on level sets, connected components, and a novel variation of the restricted centroidal Voronoi tessellation that is better suited for spatial statistical decomposition and parallel efficiency. The resulting data structures organize features into a coherent nested hierarchy to support flexible and efficient out-of-core region-of-interest extraction. Next, we provide an efficient parallel implementation. Finally, an interactive visualization system based on this approach is designed and then applied to turbulent combustion data. The combined approach enables an interactive spatial statistical analysis workflow for large-scale data with a top-down approach through multiple-levels-of-detail that links phase space statistics with spatial features.


INTRODUCTION
Simulation scientists continue to push computational limits for the modeling of complex spatial systems. These simulations can generate extreme-scale data with high spatial and temporal fidelity, but I/O bottlenecks and storage constraints place limitations on how much of this data can be saved and analyzed post hoc. This problem motivates the development of in situ data analysis and processing workflows with the goal of reducing I/O while retaining crucial information.
One way to approach this problem is spatial statistical aggregation. Such approaches usually involve segmentation or decomposition in the spatial domain, and then statistical aggregation within the resulting segments. Several methods have been developed for spatial segmentation of statistical aggregation, such as blockwise decomposition, Voronoi decomposition, and voxel clustering. These methods have demonstrated value; however, it remains a challenge to organize such a tessellation around dynamic flow features, typically defined by surfaces.
To address this problem, we segment the spatial domain dynamically using level set restricted centroidal Voronoi tessellation (CVT). The CVT may also be weighted to vary the spatial fidelity based on a spatial saliency function. Since the decomposition is aligned to the level sets, this also provides a hierarchical segmentation, from isbands (voxels between pairs of isosurfaces or contours) to connected components within isobands, to the centroidal tessellations within the connected components. The statistical aggregations may also be merged over Voronoi cells and throughout the hierarchy. For example, one may view statistics of the Voronoi cells individually, or locally merge Voronoi cells to obtain aggregations over multi-segment volumetric patches that are aligned to the surfaces of interest, or merge all Voronoi cells within a connected component or entire isoband.
For practical purposes, the implementation must be scalable and efficient. Efficient algorithms have been developed for computing CVTs; however, the focus of these methods has been on different application areas, such as re-meshing and construction of computation grids. These other applications have requirements that Voronoi must conform to, and are typically defined through precise geometry. For statistical aggregation, however, the requirements may be relaxed, since heterogeneity in region size but conflict with the goal of generating a conforming mesh in complex geometries. Meanwhile, the benefit of the discrete tessellation is that it can be done efficiently in parallel on GPUs and it fits the analysis workflow whereby voxel elements are sub-selected discretely for visualization. Since our requirements differ from the typical use case, and because we are aiming for high performance, we implemented our own discrete method for a non-conforming level set restricted tessellation that favors volumetric and geodesic regularity.
Additionally, interacting with the data in this form is difficult without a supporting visualization system. As scientists put great preparation into planning an expensive, large-scale simulation run, it is important to also adequately prepare the configuration of the in situ workflow so that it can be relied on without concern. For example, with our approach, one must choose level sets, overall Voronoi site density, and the optional saliency function for concentrating site density. For statistical aggregation, there are several possible models for distribution functions, each with different possible parameters that should be tuned appropriately. Among models and their parameters, there are also tradeoffs affecting performance, storage costs, fidelity, and other measures of quality. Thus a working visualization system is a valuable tool to have before one attempts to apply the in situ methods. Our visualization system is able to help with this, since it can be used on existing data where the methods can be tested and the aggregates can be compared with the raw data to get an idea what information may be lost or misrepresented by the aggregates. Figure 1 shows the system interface.
Our primary contributions are • a parallel level set restricted CVT algorithm and implementation, • a supporting visualization system, and • an overall approach to visualizing large scale multivariate spatial statistical information. We apply our work to simulation data of turbulent combustion and wall-bounded turbulence, and we evaluate the performance and scalability on the Summit supercomputer at the Oak Ridge National Laboratory.

RELATED WORK
This work is inspired by past success using spatial statistical aggregation to reduce and analyze large simulation data. Chadhuri et al. introduced scalable distributed algorithms for computing histogram-based spatial summaries of large 3D mesh data [13]. Lu et al. introduced compact histograms representations based on bit-strings and space-filling curves [35]. Histograms have been used to summarize space phase distributions of statistical super-particles in particle-in-cell simulation data [28,39]. For simulation monitoring and anomaly detection in transonic jet engine simulations, distribution functions were used to quickly and compactly represent flow conditions, so that they could be monitored in real-time for human in the loop stall detection [19]. Ye et al. decomposed the domain of combustion simulations into blocks, then computed spatial statistics for querying Lagrangian particles based on histogram and coupled data structure [54]. Many more related methods have been introduced. Gaussian mixture models (GMMs) and copula functions have provided small footprint models of joint statistics [20,23]. Additional work has focused on the recovery of the spatial distributions [48]. In terms of spatial segmentation for statistical aggregation, further approaches have been applied, including kd-trees, and simple linear iterative clustering [20].
The aforementioned work [19,39,54], each included interactive visualization systems for local spatial statistical analysis with distribution functions computed over local regions, but were each designed with different goals or for application domains. Neuroth et al. introduced a design which relied on a 2D decomposition for Tokamak data [39]; where symmetry justifies radial integration to project the data onto the 2D plane. They used a static Voronoi tessellation. Like our system, their system also supported selecting and merging the statistics in different spatial regions, which could then be analyzed through linked plots. However, they only had static plot parameters. Ye et al. [54] utilized a Cartesian spatial decomposition of a 3D domain for combustion data. This has an efficiency advantage but cannot support feature-aligned spatial statistical analysis. Dutta et al. probed spatial statistics around the surface of engine blades [19]. Their system provided human-inthe-loop anomaly detection utilizing a streaming timeline visualization. While these systems may be considered state of the art in their respective ways, we are not aware of any approach and supporting system which can provide the type of interactive statistical spatial analysis around features that our system supports.
Besides tailored visualization systems for large-scale spatial statistical analysis, there are state-of-the-art tools for accommodating many general visualization tasks, both in situ and at scale. Paraview [4] is one such platform, and Visit [14] is another such tool. These platforms are indispensable for supporting large-scale visualization and analysis across an extremely wide range of dataset formats and scientific domains. However, while they can also be heavily customized and used in integrated workflows, they do not have the capability by default to replicate the interactivity and features of our system for linked spatial statistical analysis. Nevertheless, alternate or similar workflows could be used, perhaps with an integrated approach incorporating multiple tools. The visualization toolkit (VTK) [43] pipelines would work well to support the data processing pipeline of the Voronoi processing and statistical aggregation. Our algorithms, which were implemented in CUDA [40] (NVIDIA's GPU computing API) and Thrust [8] (a high-level interface for data-parallel computing), could also be implemented using VTK-m [38] (a toolkit of parallel scientific visualization algorithms), which is similar in its model to Thrust and CUDA.
Besides approaches based on statistical modeling, there has been a great deal of work done on sub-sampling data from large complex domains while preserving surface and topological features in flows [2,6,51,55]. Topological methods for visualization, in general, are too numerous to mention. Heine et al. provided an in-depth survey [24]. Carr et al. introduced efficient parallel algorithms to compute contour trees [11,12]. Bremer et al. proposed using merge tree as a way to semgement combustion data into hierarchical features that can be summarzied through statistics and explored efficiently through their graphs [10]. Another great idea that has pushed us forward is the work on summarizing large time-varying simulations through the use of dynamic nested tracking graphs [37] in a coupled interface that links the topological features to in situ rendered images within a cinema database [5]. In addition, recent work has offered a generalization of isosurfaces to multi-fields, which inspires an entirely new approach towards level set and topological based analysis in multivariate data [25]. These are just a few examples, as the area of topological segmentation is a rich and active area of research. Our method utilizing a restricted centroidal Voronoi tessellation (CVT) is designed to compliment these methods, by offering an additional capability to decompose the topological segments further without losing the benefits.
Voronoi tessellations, also called Voronoi diagrams or Voronoi decompositions, are among the oldest and most fundamental methods in computational geometry and related fields [7,17]. The CVT is also useful for glyph placement and stippling [18,21]. The centroidal Voronoi tessellation, where each Voronoi site is at the centroid of its region, is commonly computed using Lloyd's Algorithm [36]. A centroidal Voronoi generates a segmentation with more homogeneous region sizes and shapes, which is more favorable for spatial statistical analysis. Our contribution is a CVT algorithm that is designed for spatial statistical analysis and summarization of flow features, which are complex in shape and encapsulated by level sets.
CVTs with restrictions/constraints have been introduced [47], however this has been made possible in complex shapes by concentrating more Voronoi sites around the sharp/complex features, often with the goal of geometric modelling and surface simplification or reconstruction. For example, Amenta et al. introduced PowerCrust, which is able to generate a 'watertight' boundary 3D polyhedral solid, with theoretical guarantees on quality. VoroCrust expanded on this line of work with additional provable guarantees under more general conditions [1]. Once a 'conforming' decomposition of the boundary has been constructed, the interior could then be decomposed into a centroidal tessellation using existing algorithms. Voronoi tessellations on manifolds and within non-Euclidean spaces based on geodesic distances have been introduced as well [41,49], where the geodesic represents the shortest path between two points on the manifold. The manifold could be surfaces in 3D space or higher dimensional, e.g. the 3+1D hyperbolic space in general relativity. Note the geodesics in our method are not geodesics on a manifold, but rather graph theoretic geodesics representing shortest paths obeying the level set restrictions.

Motivation, Research Challenge and Novelty
Restrictions and Homogeneity. In Section 2, we discussed the state of the art in creating Voronoi tessellations around complex restrictions. The state of the art methods work by building a conforming crust around the boundaries. But in building the crust so that it is conforming, we must concentrate more sites around sharp features, and thus we will sacrifice homogeneity of the Voronoi regions on the boundaries. But the regions on the boundaries are those of primary interest for spatial statistical analysis. This is why we relax conditions required for surface reconstruction, and instead focus on homogeneity of the region sizes. However, this creates challenges for a restricted CVT algorithm since individual regions will wrap around the flow features, and so the normal centroid is no longer valid since it may not even be within the region, and the centroidal updates could carry a Voronoi site through a boundary/restriction. Therefor, we create a novel restricted CVT algorithm based on graph geodesic distances, where the graph is constructed within the volume for each voxel to have a shortest path around the restriction boundaries to the graph geodesically nearest Voronoi sites. Besides CVTs, homogeneity guided decompositions for spatial statistics have been approach using simple linear iterative clustering (SLIC) [20], however, to our knowledge, an algorithm for restricted SLIC is not available, and since the spatial component of SLIC is based on Euclidean distance, a restricted version also could not reliably meet complex boundary restrictions while maintaining spatial homogeneity of region sizes.
Scalability and Computational Performance. Since the domain is evolving, and one application of the decomposition is for in situ processing, it is important that our algorithm is efficient, and exploits the heterogeneous parallelism of modern super computers, including distributed data parallelism, and GPGPU parallelism. Comparisons between our algorithm and others are difficult since to our knowledge, no other algorithm computes the same form of CVT as ours. However, algorithms have been proposed for standard CVTs on 2D planes and surfaces which exploit GPU parallelism [42], and in 3D parallel clipped CVTs [52]. A general data distributed parallelism for CVTs was intruduced by Starinshak et al. [45], which also handles restrictions only by clipping. However, clipping will leave randomly fragmented regions along the boundaries depending on the data.

METHODS
First, we introduce a local statistical aggregation approach based on a discrete geodesically based level set restricted CVTs. It is geodesically based in a graph theoretical sense, with voxels as nodes and the geodesics being the shortest restricted paths from each voxel to its nearest Voronoi site. We then present a prototype system that utilizes the method for hierarchical spatial statistical analysis.

Level Set Restricted Voronoi Tessellation
Level Sets, Isobands, and Connected Components. We define a set of disjoint discrete isobands where each isobands B a,b is the set of voxels which have values between scalar values a, b, i.e.
is adjacent to any voxel in C j a,b for any pair (i, j). Each of these components is a bounded region which is to be tessellated by a CVT.
Voronoi Sites and Regions. We define the Voronoi site s is a defined by a point x s , and its Voronoi region restricted to the component C as where S C is the set of Voronoi sites in component C, and gd(x v , x u ) is the graph geodesic distance from voxel v to Voronoi site u, where each voxel may be node in the graph, and no edges in the graph cross any boundaries, so that the distance gd(x v , x u ) is the shortest restricted path from the voxel to the Voronoi site. Centroids or Centers of Mass. A standard discrete CVT in a volume is a Voronoi tessellation where each Voronoi site point is at the centroid or center of mass of the region it encapsulates, where center of mass is defined as 1 m(R s ) ∑ v∈R s m v x v , where R s is the region of the site s, v ∈ R s are the encapsulated voxels, m v is the weighting function evaluated at volume element v, x v is the position vector of v, and m(R s ) is the net mass for all v ∈ R s . In the unweighted version, However, by this definition, the center of mass of a Voronoi region could be outside the boundaries, which would make the centroidal distribution impossible, create disconnected Voronoi regions, and break Loyde's algorithm (which repeatedly moves each site to the center of mass). To avoid escaping its component, we can move the site towards the centroid only until it runs into the boundary, but the sites can still remain stuck moving into the obstacles instead of around them. Instead we compute a weighted centroid of only line-of-site (LOS) voxels. The LOS voxels are those for which a strait line between them and their Voronoi site does not intersect a boundary. To account for voxels which are not LOS, we propagate their mass to the line of site nodes in their geodesic path back to their Voronoi site. This leads us to a new definition of the centroid x c , where instead of

Parallel Algorithm and Implementation
Our algorithm is based on Loyde's algorithm [16], which repeatedly computes the Voronoi tessellation (VT) and then updates each Voronoi site to the centroid of its region as shown in Algorithm 1. The difference is we use the geodesic distances defined in Equation 1, and illustrated in Figure 3, to compute the VT, and we update the site positions based on x c as defined in Equation 2. The full process is illustrated in Figure 4. For distributed parallelism, each rank computes a LSRCVT for one block of data independently. The result is that the CVT is aligned to both the level set restrictions and the block boundaries, as can be seen in the large scale tessellation in Figure 1. The benefits and limitations are discussed in Section 7. Within each block, GPU parallelism is used to compute the CVT. The code is written in C++ using MPI for distributed parallelism, and Thrust [8] and CUDA [40] for GPU parallelism. Due to space limitations, the full implementation details are included in the supplemental material.

Stratified Random Initial Site Distribution.
The initial Voronoi sites are generated to cover the set of connected components with a target number of sites n c in each component chosen based on its net mass and a site-density parameter α using the formula Where γ is a parameter, the attraction, that reduces or increases the effect of the weighting, m(C, γ) = ∑ v∈C m(v) γ , and m(v) is the voxel's weight. We also use a block-wise spatial stratification, where each component is segmented into pieces and then the n c sites are distributed fairly over those pieces. This is illustrated in Figure 4 row 1 (r1), columns 3,4, and 5, where the components (c3) are split into blocks (c4), the net mass of each component block is computed (c5), and then sites are placed based on the net mass 3 (r2,c1).
The process is done in parallel on the GPU. First we compute the net component masses m(C, γ) and the net component block masses m(b i , γ). Next, the n c sites are distributed by first seeding n c sites in each block, and then placing the remaining sites to the most underrepresented blocks based on the equation,

Discrete Restricted Geodesic Voronoi Tessellation.
The Voronoi tessellation, defined by Equation 1, is computed using a GPU parallel region growing approach emanating from the Voronoi sites. The voxels communicate with their neighbors to determine if there is a possible shorter path, and the process continues until all voxels have paths and no voxel can find a shorter one. The paths P v are represented recursively through a map src(v) : V −→ V , and can be constructed by the formula Equation 2) is the node in the path which is in line-of-site of s. Each voxel also has an 1-byte bit-field, state v . One bit indicates los(v, site(v)) (whether v is line-of-site of its site), another indicates active(v) (whether v has a current path), and another indicates node(v) (whether v is a valid node).
In the initial step, the nearest volume element to each site is classified and activated. Then, repeatedly until convergence, for each voxel v in parallel, we ask each other adjacent voxel w if it is active, and if it is, then their are two possible paths for v back to s: The former path (2) might intersect a boundary, so a ray must be cast from v to src(w) to check for a collision. If v finds a shorter path based on the best of these possible choices, then it records the distance and new src(v), but the update is delayed until a next kernel invocation to avoid race conditions. In order to determine which voxels are LOS to their Voronoi sites, in a initial stage of the above procedure, we only allow s to be active nodes so that all classified voxels are line-of-site. Then we activate the rest of the voxels as nodes and perform the procedure repeatedly until voxels have paths and cease to find a better path.  . This is illustrated in Figure 4 row 2, column 4, and Figure 3. In a GPU kernel parallelized over voxels, we sum the moments m v x v of each line-of-sight voxel in each Voronoi region. The Cuda atomic add operator is used for the summation. The geodesic weighting comes into play as the line of sight voxels which are path nodes that connect voxels which are around corners are weighted by the amount of voxels that they are sources of.
Next, a GPU kernel parallelized over the Voronoi sites is invoked, and the resulting site's sum of moments from the previous stage is normalized by the net mass to give the vector x c from Equation 2, which is the new target position for the x s . We then cast a ray from x s to x c and move x s only as far as possible before crossing the boundary.

Data Reduction and Structured Data Layout.
Oftentimes, flow features of interest will comprise only a small portion of the full data. In this case, it makes sense to summarize and/or extract only the data within the level set restrictions of our tessellation as a way to reduce the data. Our method can optionally extract the raw data and couple it with the statistical summarization in a structured data layout before saving on disk. The data structure is hierarchically ordered by level set and connected components as shown in Figure 5. This ensures the components and layers are directly indexable and contiguous in storage, so one can directly load only the feature of interest for an efficient out of core approach. The statistical summaries, which are lighter weight, can be loaded first, explored, and then other data for selected features of interest can be loaded selectively.
The reduction from discarding the regions beyond the isobands can be approximated as (n − r) × (m + 1) + n l + n c where n is the number of volume elements in the full domain, r is the size of the subset within the levels, and m is the number of variables per volume element. n l and n c are the number of layers and components, respectively (for storing the offsets) which should be very small relative to r. Thus, when (n − r) is large and or m is large, the additional cost of storing coordinates becomes nearly negligible. In extreme-scale combustion simulations, m can be greater than 100, in which case the coordinates can be as little as 1% of the total footprint of the reduced data.

Role and Choice of Parameters
Tessellation. The tessellation is both dynamic and data dependent, and also depends on a few important parameters. First, one must have the scalar field to use for the level set restrictions. For example, in our combustion application we may choose temperature isotherms, or we may use the distance function from the flame surface. The set of isovalues must then also be chosen. If the method is used in situ, one must know in advance what kinds of value ranges the simulation will produce, and what isobands will be of interest for representing the flow features. In addition, if the tessellation is weighted based on a second scalar field, then this adds another aspect to what must be known in advance to appropriately. The parameter γ controls the attraction force is equivalent to applying a non-linear re-scaling of the second scalar field used for weighting, and the effect it has depends on the data. The site density parameter α represents the choice of the number ratio of Voronoi sites to voxels. So this parameter can reasonably be decided in advanced with limited knowledge.
Aggregation. When aggregating data for spatial statistical summarization, information is lost if the raw data is not also output. Additionally, statistical aggregation is usually not a competitive form of pure compression, so if the primary goal is a reduced form of data to reconstruct the original data from, then this is probably not the best approach. Instead, we consider the statistical summarization as an explorable set of analysis results, which is useful for summarizing and analyzing certain aspects of the data, especially the dynamic spatial multivariate correlation structure. In this paper, we assume existing methods for statistical summarization of the segments in the tessellation will be used (e.g. histograms, Gaussian's, Gaussian mixtures, etc.). It also depends on the method and application what the value of the summarization will be, how much space it will take, and how much of the raw data if any can be replaced with the summary. Our view is that the summaries have multiple opportunities for added value. First, they are lightweight and, when coupled appropriately with a visualization system, can support scalable interactive visualization of a larger domain than would otherwise be possible from using the raw data directly. Second, even if the summary is not used to replace the raw data, it represents an analysis result that would be useful anyways in addition to the raw data. And since simulations are already commonly unable to output full raw data at each time step due to I/O and storage limitations, outputting the more lightweight summaries at higher temporal fidelity would add value to the simulation result without reducing any information at only a nominal extra cost. The lighter weight summaries can be possible to stream through a network for in situ monitoring. Finally, considering the hierarchically structured data layout that our tessellation provides, the summaries corresponding to flow features can be extracted out of core, and if coupled with raw data, if it exists, can also be a basis for drilling down into the full details from the raw data.

INTERACTIVE VISUALIZATION SYSTEM
To showcase our feature-based approach, we designed and implemented a prototype visualization system. The system represents one promising way that the dynamic spatial decomposition approach could be leveraged. The system demonstrates the qualitative value of the approach, and we hope it can inspire more adaptations and ideas for similar work.
The visualization workflow follows the natural hierarchical pattern for sub-selection that is implicit in spatial decomposition. As the user drills down into details, based on layers, connected components, Voronoi cells, smaller subsets of volume elements are brought into play for spatial statistical analysis. The layers, components, Voronoi sites, and volume elements are interacted on and visualized in 3 associated sets of views as shown in Figure. 1. Since the system supports interactive multivariate statistical analysis, data from many fields may need to be active simultaneously. The hierarchical sub-selection becomes useful in this case, not only as a means to narrow down the features of interest, but also to improve performance since sub-selections at each level will significantly reduce the amount of active data.

Motivation and Design Goals
This system in its first iteration was designed based on discussions between visualization researchers and combustion scientists. Both the visualization researchers and scientists are co-authors of this paper. Since, to our knowledge, no other systems exist that could explore the same data in sufficiently comparable way, and due to the limited length and scope of the paper, we are not able to do a thorough presentation of the system and robust evaluation at this stage. Instead we view this system as a starting point, and proof of concept for how to utilize the tessellation result, which may be built on, optimized, and evaluated in depth in future work. We provide below the set of visualization and analysis goals used to inform our design. In Section 7 we further lay out our vision for improvement and evaluation of the system.
Design Goals: 1. Be able to use the system on both raw and aggregated data for testing and comparison. 2. Be able to link and visualize the multivariate statistical plots and the features in the spatial domain. 3. Support the following common statistical plots: 1/2D probability density functions (PDF's), cumulative distribution functions (CDF's), 1/2D conditional plots, and scatter plots. 4. Be able to compare different candidate models for the probability distribution functions. 5. Be able to inspect the higher order joint statistical moments and how well the PDF models preserve them. 6. Be able to drill down into the details with a focus on the hierarchical structure of the data that is decomposed around flow features.

Isobands and Connected Components
The first set of views that the user will interact with are for selecting isobands and their connected components. A tree-based visualization is used to show the nested hierarchy as depicted in Figure. 5, and another linked view represents the components based on a 2D t-distributed stochastic neighbor embedding (t-SNE) projection, as can be seen in Figure. 6 (Left). The input to the projection is a featurization of the components, for example based on the bin values of histograms, parameters of parametric statistical models, or statistical moments. t-SNE will try to group the points with similar statistical features together. The points can be selected by the user through a lasso selection tool in the projection or from the tree view. The lasso selection operation has several modes besides new selection (•), based on set combination operations with the existing selection (∩, ∪, or \).
In the projection, the points are color coded by isoband, but in the tree view, that information is implicitly encoded by level of the node. Thus, for tree nodes, color encodes other information. For example, in Figure. 5 color is used to encode µ 4 (H 2 O, p), which is the cokurtosis of hydrogen mixture fraction and pressure. Since very small components have too few voxels to get quality statistics, we color encode them in gray instead. This can also be helpful information since those small components may be particularly meaningful (e.g. in combustion they may indicate the incipient development of an ignition event).

Voronoi Cell Views
The Voronoi regions are also projected into t-SNE projection views. Figure. 6 (lower left) shows such a projection. After selecting a set of components, the Voronoi cells in those components are highlighted in the Voronoi region projections and become selectable using the same lasso tool with the same set operator modes.
In addition to t-SNE projection based on statistical features, we found that a projection based on spatial similarity is a useful option for one or more of the linked Voronoi region projections, since it can group points based on spatial proximity. In Figure 1 (right) the bottom projection is based on such a 2D spatial projection. For this projection, we use a pre-computed distance/similarity matrix based on a custom pairwise distance measure between Voronoi sites. Since the tessellation is restricted by level sets, it makes sense to utilize the geodesic/path distance obeying the restrictions, as illustrated in Figure. 5 (A). The blue line shows the Euclidean distance, while the white line shows geodesic distance. An example where this useful in capturing the condition where regions around a folded surfaces (e.g. flame) are geodesically distant but yet close together in Euclidean distance. In combustion, for example, this may preclude or represent interesting interactions of the flame as the Euclidean distance approaches the laminar flame thickness while the geodesic distance remains larger. One example distance function that can be used is, where P i, j is the graph geodesic of adjacent Voronoi sites, gd is the path's distance, c ≥ 1 is a parameter where a higher value results in more separation of sites outside of the same component. To find P i, j , we compute the connectivity of the Voronoi sites based on region adjacency, and then Johnson's all-pairs shortest path algorithm (which internally uses Dijkstra's algorithm) is used. Since the neighborhood is also restricted to within connected components the computation is parallelized by letting different threads operate on different components. The Boost graph library's implementation is used, which has time complexity O(|V | · |E| · log(|V |)), where V is the set of nodes, and E is the set of edges. An alternative to Johnson's/Dijkstra's is to use the Floyd-Warshall all pairs shortest path algorithm which has time complexity O(|V | 3 ). Since our graph is only locally connected it is sparse enough that O(|V | · |E| · log(|V |)) is better than O(|V | 3 ).

Interactive Joint Plots
Once a selection of Voronoi regions is made, statistical plots are generated, either from their pre-aggregated statistics, or from the raw data loaded out of core based on the aforementioned structured data layout.
Based on the experience of the domain scientists, the following statistical plots are chosen: joint PDF's and their marginals, joint CDFs, conditional plots (which estimate one variable conditioned on 1 or 2 others), and scatter plots. Since each of these plots share the same axis with the same scaling, we can switch between them, or combine layers on a single interactive generalized joint plot.
The joint plot includes information from both the background layer which is generated from the set of selected components, and foreground layer which is generated from the set of Voronoi regions within the components. The background layer in the joint plot is a static gray scatter plot, while the foreground layer is plotted over the background and depends on the selected plot mode. There are also marginals for both the background layer which are shown as gray histograms drawn as bar charts, and marginals of the foreground data as heatmaps beneath the bar charts as seen in Figure 6. The user can zoom and pan by scrolling or dragging the mouse when within either the joint plot (to zoom/pan both axis), or over one of the marginals (to zoom/pan only one axis). Highlighting on the background marginals, and arrows are used as a form of mini-map to depict the zoom level.
The joint plot of the foreground layer can be either scatter plot, histogram, CDF, kernel density estimation (KDE), or Gaussian mixture models (GMM's). The KDE plots are useful for capturing high-level detail and smoothness in the potentially complex nonlinear distributions. Histograms and GMM's are models which are commonly used for in situ statistical summarization. By using our system on existing data, before employing in situ, the choices of models and their parameters can be tested and the user can get an idea what information might be lost or misrepresented using the candidate methods. As the user interacts with the data, the PDF's may be interactively recomputed as the data inputs and the data range changes. Since the KDE and GMM plots are often not able to be recomputed immediately at interactive rates, the plot will automatically switch from KDE or GMM to a histogram, for example, while zooming/panning is occurring.
To evaluate the PDF's, the system estimates how well they preserve joint moments, which are used in the analysis of turbulent flows [9,44]. The moments (i.e. mean, covariance, coskewness, and cokurtosis) estimated from the PDF's are compared against the moments computed from the raw data. Each plot has a button (in the upper right corner, labeled µ k ) that can be pressed to open/close an attached table that is auto-generated as a LaTeX document as shown in Figure. 1. This functionality can also aid featurization of high-order moments-based anomaly detection algorithms to [3,27].
If voxel data is available, the voxels deriving the foreground can be selected using the lasso tool and aforementioned set combination operators. The selections are linked across views and plotted as blue points. This layer can be toggled on/off per-plot to manage occlusion. These selections are also linked in the 3D view where their bounding surfaces are rendered. Conditional plots can be toggled on/off as well, and are rendered on the top most layer as a red curve estimating mean y as a function of x, with blue lines above and below the curve at ±1 standard deviation. 2D conditional plots, showing estimated mean z as a function of (x, y), rendered as a heatmap, is also supported. For these plots, the forground layer can also be toggled to show a scatter plot of (x, y) as shown in Figure 6 (second from right).

PERFORMANCE
The data used to obtain these results is the hydrogen/ammonia/nitrogen flame data described in Section 6. Five isobands of temperature were used. The bands are highly contorted as can be seen in Figure 7. This makes the data a good stress test for the algorithm as the Voronoi sites need to find their way around complex obstacles.  Fig. 8. Computation time per-node for convergence to within a mean d s of 0.5 (unit voxel length). Each line represents a different block size for distributed computation. Since the computation is independent for each block of data in a distributed setting, the scalability can be assumed to depend primarily on the number of GPU nodes available.
First, we evaluate the computation time with respect to site density and number of voxels on a node; shown in Figure 8. The performance drops as the the number of voxels on a GPU node increases. The primary bottleneck is the Voronoi classification described in Section 3.2.2. This computation is parallelized on the GPU using one thread per voxel without communication. Thus we can expect the GPU parallelism to scale well with the number of GPU threads. However, the performance decreases approximately linearly with respect to site density, as lower site density means longer paths traced out per voxel/thread. The algorithm converges similarly in each case to low mean d s (less than 0.5 voxel widths), very rapidly. The continues to improve slightly over time with diminishing returns on the benefit for the cost.
The second result, in Figure 9, shows how quickly the tessellation converged in terms of steps in Algorithm 1. Diminishing returns for the computation cost are seen after only a handful iterations. Since the computation is independent over each node/block of the block decomposition domain, we were able to use a single user machine to perform these first two tests; an Arch-Linux machine with an Intel(R) Core(TM) i7-5930K CPU @ 3.50GHz, 64 GB of main memory, and an NVIDIA GeForce GTX TITAN X which has 12 GB of memory. The number of Cuda threads per-block was set at 128, which gave us the best performance. The Cuda code was compiled with compiler flags, 'nvcc -std=c++11 -O3 -w -m 64 -rdc=false -ftz=true -use fast math -Xptxas -O4, -Xcompiler -O4 -arch=sm 52'. Last, we performed a full distributed parallel test for the combustion data on the Summit supercomputer. The volume has a size of 3456 × 1280 × 2560 voxels. The site density was 2700 voxels per Voronoi site. We ran the system until it reached the point where mean d s dropped below 0.25 unit voxel widths. The results are shown in Table 1. Each of the nodes on Summit has 6 NVIDIA V100 Tesla accelerators, 512 GB of DDR4 memory, and two IBM POWER9 CPUs. We ran our code with 1 GPU and 1 CPU per rank, and 6 ranks per node. We testing the performance using 48,96,192,288, and 360 GPUs. Since Summit has 6 GPUs per-node, we mapped this to 8,16,32,48, and 60 nodes respectively. The run command for example was 'jsrun -n48 -a1 -c4 -g1 ./t8' for the case of 48 GPUs. The CUDA code was compiled with NVCC, with command-line arguments: -std=c++11 -O3 -w -m 64 -rdc=true -Xptxas -O3. The CUDA version was 11.4.2. The C++ program that links the GPU code was compiled with GCC 10.2.0 (the GNU Compiler Collection) with arguments: -std=c++11 -O3 -fpic -fopenmp. Since each node took a different amount of time to complete, we report the worst case, which represents the overall time for the entire data to process. The CUDA kernels used a block size of 128.

APPLICATIONS
We tested our algorithms on two types of flow data, turbulent combustion and wall bounded channel flow. The flame surface in the turbulent combustion data is more connected, but still highly complex, while the channel flow data has many small vortical structures. The combustion data features strong coupling between the flame wrinkling and the underlying turbulent strain, where the flame response occurs in a composition phase space (dependent locally on species concentrations and temperature which control the burning rate and are modulated by the turbulent strain), while turbulent strain evolves in physical space. Hence, there is an inherent need to interactively toggle between physical and composition space to drill down and understand causality. Depictions of the two datasets and tessellations made using them are shown in Figure 7 and Figure 10, respectively.

Hydrogen/Ammonia/Nitrogen Turbulent Flame
We apply the nested hierarchical approach presented here to turbulent combustion data from Direct Numerical Simulation (DNS) of a premixed hydrogen/ammonia/nitrogen flame [50]. The DNS study aims to gain an understanding of the combustion behavior of hydrogen/ammonia/nitrogen as a carbon-free replacement for conventional fuel (i.e., natural gas). The simulation was run on 900 nodes of the Oak Ridge National Laboratory OLCF Summit cluster requiring approximately 110k node-hours. The aggregate amount of data produced comprises 400 checkpoint files, where each checkpoint file is 2TB.
Due to the nature of turbulent combustion and the need to resolve all flow and chemical scales, turbulent combustion DNS data is characterized by a large separation of scales, with portions of the data often not relevant for statistical analysis or over-resolved on the numerical grid that is designed to capture the smallest scales (e.g., regions of hot products feature small gradients/large flow and chemical features). The Voronoi tessellation effectively reduces the amount of data to be processed for the purpose of statistical analysis and allows for efficient and fast data exploration. While conventional statistical analysis (e.g., (joint) statistical distributions, statistical moments, correlations between species and temperature and between reaction rates and strain rate) is important in turbulent combustion, the analysis of spatially and temporally intermittent features is often of interest (e.g., flame-flame interactions, reactant pocket formation enriched with minor species that are chemically crucial, and formation of localized flame extinction holes). The hierarchical selection of connected components together with distance-based metrics, and with the linked views between spatial and statistical features, allows scientists to quickly identify and select specific flame features and directly extract necessary statistics to understand causality between the 'turbulence-chemistry' interactions and to steer further downstream analysis, and adaptive labeling of 'feature' data for searchable training and validation for machine learning.
The mass fraction of H 2 O, normalized by the values obtained in burned and unburned gas of a laminar, unperturbed flame at the same conditions as the turbulent flame, acts as a so-called progress variable, which monotonically increases from 0 to 1 across the flame brush. This progress variable is usually used to identify specific zones of a premixed flame and specific values can be selected to obtain iso-surfaces that identify the various zones in the flame front (i.e. preheat zone, reaction zone, oxidation zone in a premixed flame). Here, we choose a specific value that represents the reaction zone, i.e., the location where the heat released by the flame is at its peak. We refer to the iso-surface corresponding to this value as the flame surface. Of particular interest is the analysis of the value of other quantities on the flame surface. For example, a large value of the OH radical species represents a region that is burning particularly strongly and in the present dataset, this corresponds to a flame element that experiences preferential diffusion effects (i.e., hydrogen reactant diffusing faster than heat), where the burning is amplified locally by a large influx of fuel. In other regions of the flame, OH is low, which represents flame elements that are locally extinguishing, i.e., strain and stretch-induced by turbulence disrupt the reactions to a large extent. Distance-based metrics also provide a conditioning variable in physical space indicating the response of different zones in the flame, upstream and downstream of the flame isosurface, to the strain rate. Distance function iso-contours represent the locus of all points that are equidistant from the flame surface. Fig. 10. Channel flow decomposition using the λ 2 field to encapsulate vortical structures, and the distance function to define iso-levels. The dark grey surface represents the vortical structures.

Turbulent Channel Flow at High Reynolds Number
Here we examine is the incompressible turbulent channel flows at a high Reynolds number, Re (Figure 10). High Re wall-bounded turbulence has great importance in many engineering applications with moving objects at high speed. However, the fundamental dynamics of wallbounded flows are not yet fully understood. The turbulent channel flow is one of the simplest canonical flows to study the high Re wall-bounded turbulence. Still, the DNS of high Re turbulent channel flows is quite expensive due to the scale separation similar to the turbulent combustion as explained above. The simulation performed with 33K nodes in Mira (Bluegene/Q, ALCF) used 242 billion DOFs, and it remains the largest DNS simulation of wall-bounded turbulence [29][30][31].
In this study, we use the λ 2 field, which is calculated from a subset of the aforementioned DNS. The interaction of coherent structures in turbulent flows is a long-standing problem. Naturally, the methods to identify vortex structures are well developed. λ 2 is one of the widely used vortex identification methods. It uses the eigenvalues of a 3 × 3 matrix constructed with a velocity gradient tensor [26]. Even though the vortex identification method is mature itself, using it to study the dynamics of high Re wall-bounded turbulence is still a challenging task because the identified vortex structures are too numerous to analyze readily from the raw data [34]. Note that other types of high Re turbulence also encounter the same flow complexity. The dynamic nested hierarchical statistical decomposition approach can efficiently identify the vortical structures and analyze their interactions with other state variables, such as velocity and pressure fields.

DISCUSSION
Block Based Distributed Parallelism. For large-scale domains, a distributed approach is needed. There are two possibilities; (1) perform the CVT globally over the whole domain, which will be computationally intensive due to the communication requirements across nodes, or (2) compute the CVT separately on large blocks of the domain on different ranks. We have only implemented the second option. The result is that the boundaries of the blocks will act as further restrictions so that the Voronoi sites will be aligned to those boundaries in addition to the isobands. Such alignments are already implicit at the boundaries of the full domain, and the alignment to the block boundaries does not break the level set based hierarchy, or alignment of the CVT to the level set surfaces. Additionally, it offers one advantage over a global CVT; whole blocks can be loaded without clipping the tessellation. This is practical because large datasets are often analyzed post-hoc through block-wise sub-selection due to the large scale. It also makes a more efficient implementation possible and is trivial to scale up.
CVT Algorithm and Optimality. In general, it is difficult to prove guarantees for CVTs [15]. Liu et al. provide theory on CVT convergence [33]. Despite theoretical challenges, CVT algorithms adaptively minimize energy and are argued to have practical expectations of convergence. For our algorithm it depends on the definition of x c , the new target site position in the update step for Loyde's algorithm. As site density increases and/or the restriction complexity decreases, x c approaches the Euclidean centroid. While the tessellation will evolve from less optimal to more optimal distributions, we are not able to prove theoretical guarantees. A more optimal definition is also possible, but it is unclear what that should be in order to balance efficiency and quality. Our implementation can be extended or improved in the future by redefining x c and modifying the update step. As an alternative to Loyde's algorithm, one can also try an approach based on limited memory BFGS method for large scale optimization (L-BFGS) [32].
Need to Choose Parameters A Priori and Loss of Information. As discussed in Section 3.3, there is a significant limitation in the need to specify parameters before hand, and the choice of parameters depends on dynamic data. This becomes a bigger concern when one is less clear what kind of data they can expect the simulation to produce. When one has some existing data which they expect to have similar characteristics to the data the simulation will produce, then they can apply the methods on the raw data and use our prototype system to explore the results while being able to compare what information may be lost from the raw data given a parameterization of the methods. As future work, one may explore robust dynamic methods to adapt the statistical aggregation/modelling approach as the data evolves based on the uncertainty or error. For example, when a model is no longer able to represent the data sufficiently, one can trigger more output of the raw data, or employ a more expensive yet superior algorithm or model.
Applicability and Extensibility. Our implementation works on Cartesian grids and rectilinear grids. For rectilinear grids, one can perform the tessellation as if it were a Cartesian grid in which case the tessellation will be stretched to match the grid. In some cases, this approach is sensible because it results in more homogeneity in region sample sizes (volume elements per region), which may yield greater statistical consistency. However, one can also use the stretching function for the tessellation. The algorithm can be extended to work in most cases where a discrete binning of volume elements is acceptable.
Impact of the Visualization Approach While visualization tools for 3D rendering and surface analysis of the large scale data have matured over the years, there is still a lack of scalable tools for exploring large scale multivariate data and spatial statistical correlation. Our work contributes new ideas and software to help push the state of the art in this area, and can be expanded on by future research.
Evaluation and Refinement of Visualization Design. The visualization system and overall approach has clear motivation and has been applied to help combustion scientists better understand their data. The design is still in an early stage and represents only one for approach visualizing spatial statistics using the hierarchical segmentation. We have determined several ways it can be improve. First, a direct selection of connected components in the 3D view would be useful. This could be supported either based on an explorable image with depth buffers [46], through 3D ray intersection testing. Second, we would like to region growing of the selections based on spatial or statistical proximity. Lastly, the visualization system needs better support for temporal visualization. The statistical summaries of components can be efficiently visualized through time utilizing a dynamic tracking graph [37] and the structured data layout from Section 3.2.4. Once the visualization system has matured, a more robustly evaluation can be done. Evaluation of methods designed for use by qualified scientists for exploratory visualization and analysis can be challenging, because few are qualified to participate in studies, and appropriate tasks are difficult to define for visualizations when the result of the visualization is qualitative insight. If a user study is overly simplified or narrow in scope, it can give a false impression of the usefulness and potential of the design [22]. Feedback from the participating scientists on the value of the methods and tools is included in Section 6. However, the qualified scientists we have worked with are co-authors and helped design the system around their interests, so there is a level of potential bias that may be expected in their direct feedback, and the potential use cases for the methodology spans domains we are not experts in. Our vision for the future evaluation and refinement of the system involves iterative refinement based on lessons learned as we apply it in our domain [53], as well as through as much feedback as can be gained by other users outside our domain.

CONCLUSION
We believe that the workflow presented in this paper can help drive future scientific researchers closer to realizing the promises of ever increasing large and complex scientific data.
In the near term, we envision semi-automated in situ workflows where the scientist needs to keep up with and play an active role in the process from simulation to analysis. The ability to identify statistical patterns and their characteristics manifested in both physical and high dimensional phase space coordinates will enable 'feature or event' based parameterization. As part of a holistic in situ dynamic workflow, and the detection of these features could trigger further statistical spatiotemporal analysis, 'event' labeling and storage for data driven reduced-order surrogate modeling, adaptive mesh refinement, or adaptive in situ manifold learning as anomalous events occur.
In addition, our visualization software utilizing this methodology can be expanded to support a workflow for interactive statistical featurization, and event modelling, and then can automatically export the models into the target language and workflow model/system, to be rapidly deployed for into their scientific pipelines. Ultimately, this would act as a coupling of the discovered spatial statistical pattern with a reproducible workflow that can be automated. This could then be applied in situ, e.g. within a dynamic workflow for feature detection and subroutine triggering.