1 Background
1.1 Fiber Surfaces
We briefly summarize the fiber surface extraction technique proposed by Carr et al. [9] for bivariate scalar fields. Broadly, for a multivariate function , the spatial domain is mapped to an attribute space . The attribute space encodes scientific observations or simulations comprising scalar, vector, and tensor fields. In this paper, we consider bivariate scalar fields with and . Let be the attribute space with two attributes or scalar functions and . A fiber is then defined as the inverse image of a point [49, 9, 32], and it corresponds to an intersection of isosurfaces [37] for isovalues and .
The fibers of interest can be specified with FSCP, which we denote as a trait, . Each point on a trait boundary has a corresponding fiber, and such fibers carve out or represent a fiber surface. Consider the fiber surface example illustrated for the ethandiol dataset in Figure 1, which is similar to the one by Carr et al. [9]. Figure 1a visualizes a continuous scatterplot [5] of the grid vertex data with the electron density (Rho) as attribute on the horizontal axis and the reduced gradient (s) [28] as attribute on the vertical axis. The continuous scatterplot in all our examples is computed with the topology toolkit [60]. The four traits are selected as colored polygons with arbitrary shapes in an attribute space. The respective colored fiber surfaces (i.e., the set of fibers corresponding to points along the polygon boundaries) are visualized in Figure 1b.
The process of fiber surface extraction is similar to the process of isosurface extraction, which uses the marching cubes algorithm [37] but with a few tweaks [9]. Initially, each grid vertex of domain
is classified as
interior or exterior depending on whether the function values lie inside (interior) or outside (exterior) FSCP. Mathematically, a vertex is classified as interior if; otherwise, it is classified as exterior. Such classification determines the marching cubes topology case within each grid cell for fiber surface reconstruction, and inverse linear interpolation may be applied in the attribute space
to estimate fiber surface positions on grid edges (please refer to Algorithm 1 of [9]). Feature level sets [25] generalize the idea of fiber surfaces and extract features based on distance from trait . For , the feature corresponds to the fiber surface itself. For , features present a spatial evolution of fiber surfaces in the vicinity of trait .1.2 Fiber Uncertainty for Rectangular FSCP
The uncertainty in multivariate data can result in ambiguity as to whether a vertex should be classified as interior or exterior. The problem of uncertain vertex classification was recently addressed by Zheng and Sadlo [69] and Sane et al. [52]. Zheng and Sadlo proposed a statistical framework to compute the probability of data at point being in the interior of a rectangular FSCP. The computed probabilities were then visualized via direct volume rendering and referred to as uncertain rangefibers. Their analysis assumed the independent/correlated Gaussian (parametric) noise model and a rectangular shape of trait .
We illustrate the probability computation process proposed by Zheng and Sadlo in Figure 2
. Let random variable
denote uncertain bivariate data at position , where and are the two random variables that indicate uncertain data ranges for two attributes. Let and denote the probability distributions of random variables and , respectively. Let the rectangular trait be defined by ranges and for attributes and , respectively. Let denote the probability of uncertain data at point being in the interior of a trait, which we refer to as interior probability. The interior probability (i.e., Pr( = )) can be computed in two steps: (1) the joint probability distribution of and (i.e., ) is computed, and then (2) the joint probability distribution is integrated over trait . Mathematically, the probability of point being in the interior of a fiber surface can be expressed as the following double integral over a rectangular trait:(1) 
Zheng and Sadlo [69] proposed a closedform computation of the double integral in Equation 1 when
is Gaussian distributed. Sane et al.
[52]proposed confidencefeature level sets for Gaussian distributed data and rectangular FSCP for bivariate fields. In their technique, they considered confidence intervals of uncertain ranges at a position
(i.e., random variables and ) to generate respective confidence visualizations of fiber surfaces.2 Related Work
Our proposed methods mainly relate to two important research areas of scientific visualization: multivariate data visualization and uncertainty visualization. Next, we discuss prior work in each of these areas that is relevant to fiber surface visualization.
2.1 Multivariate Data Visualization
Understanding correlations among variables is one of the main challenges in multivariate data analysis. Scatterplots [15] and parallel coordinate plots [23] are fundamental tools that facilitate effective exploration of attribute spaces of multivariate data. Several research efforts have focused on improving challenges associated with scatterplot and parallelcoordinate plot visualization. Continuous versions of scatterplots [5] and parallel coordinate plots [21] have been proposed to avoid sampling artifacts of their traditional discrete versions. Quadri and Rosen [46] used topological data analysis to identify and explore clusters in scatterplots. Sauber et al. [54] proposed multifield graphs to study correlations among variables and their correlation strength.
Data variables can be encoded into glyphs [63] to visualize a correlation in a physical domain but can suffer from occlusion and cluttering issues. Other effective alternatives for multivariate visualization include brushing of scatterplots or parallel coordinate plots and linking them with the spatial domain for exploration of correlations in a physical space. For example, in Figure 1, each trait indicates brushing in the attribute space, and their respective visualization in physical space corresponds to a fiber surface. Hauser et al. [19] introduced angular brushes for exploration of parallel coordinate plots and scatterplots. Kniss et al. [33] proposed multidimensional transfer functions as a brushing technique for volume rendering of multidimensional data. Jänicke et al. [24] presented an approach for the projection of multivariate data to a 2D space (referred to as attribute clouds) followed by clustering and brushing for visualization.
We provide a brief overview of several variants and applications of fiber surfaces for visualizing multivariate data. Klacansky et al. [32] extended the fiber surface extraction technique to tetrahedral meshes and proposed a parallel computational framework for interactive exploration of fiber surfaces. Wu et al. [65] presented a direct volume rendering framework for visualization of fiber surfaces without explicit extraction of fiber surfaces. Sakurai et al. [51] extended the idea of flexible isosurfaces [10] to flexible fiber surfaces for the visualization of fiber surface components of interest without occlusion. Raith et al. [47] applied fiber surface extraction to tensorfield invariant spaces for the visualization of tensor fields. Tierny and Carr [59] derived fiber surfaces around the Jacobi sets of bivariate fields (a conceptual equivalent of critical points of univariate fields) for computing fiber surface topological descriptors known as Reeb spaces [11]. Our literature review covered a small subset of multivariate visualizations relevant to fiber surfaces. For a comprehensive overview of multivariate data analysis and visualization, we encourage readers to refer to the survey papers by Fuchs et al. [16] and Kehrer et al. [31].
2.2 Uncertainty Visualization
To date, most of the research in uncertainty visualization has analyzed noise propagation in univariate data and the associated visualization algorithms. Specifically, these algorithms, including level sets [48, 64, 45, 4, 1], direct volume rendering [38, 14, 35, 50, 2], and topologybased visualizations [66, 17, 12, 68, 3], have been extensively studied in the context of uncertain univariate data. A few studies have investigated the uncertainty in visualizations of vectorfield [36, 40, 55, 13, 18] and tensorfield [29, 26, 56] data.
There have been a few developments in uncertainty visualization of multivariate data. Zheng and Sadlo [69] proposed a framework for visualizing continuous scatterplots and uncertain fibers when the uncertainty in bivariate data is modeled with Gaussian distributions. Hazarika et al. [20] characterized multivariate data distributions with copulabased models and proposed a sampling strategy for copulabased models to quantify the uncertainty of features such as level sets and vortices. Xie et al. [67] encoded data quality into scatterplots and parallel coordinate plots for multivariate data in the context of information visualization. Nagaraj et al. [39] studied the consistency and inconsistency of gradient fields derived across attributes of multivariate data. In this paper, we take a step toward uncertainty quantification for multivariate data by investigating uncertainty in fiber visualizations of bivariate data for parametric and nonparametric noise distributions and arbitrary FSCP shapes.
3 Fiber Uncertainty for Independent Parametric Noise Models
In this section, we discuss two approaches, the closedform integration and Monte Carlo integration, for fiber uncertainty computations. Although we present a closedform fiber uncertainty quantification framework for integrable kernels (uniform, Epanechnikov, and Gaussian), a more generic Monte Carlo approach can be utilized to estimate fiber uncertainties for integrable as well as nonintegrable kernels.
3.1 ClosedForm Integration
We address the challenge of uncertainty quantification of fibers, similar to Zheng and Sadlo [69], for arbitrary shapes of FSCP. Specifically, we present a closedform framework for computing interior probabilities (i.e., Pr( = )) per grid vertex when data noise is modeled with the independent uniform, Epanechnikov, and Gaussian distributions, which are integrable. The techniques detailed in this section are then leveraged as building blocks of our derivations for nonparametric noise models in subsection 4.1. When the FSCP shape is not rectangular, the double integral in Equation 1 cannot be utilized because it considers the ranges of two attributes. For FSCP with an arbitrary shape, we apply Green’s theorem to integrate the joint probability distribution . Our approach is depicted in Figure 3.
In Figure 3, the blue region indicates the joint probability distribution of random variables and (i.e., ). The polygon indicated by the orange dotted lines in Figure 3 results from the intersection of the uncertain region (blue) with trait . Let () denote the coordinates of the intersection polygon. We compute the interior probability (i.e., Pr( = )) by integrating the joint probability distribution over the intersection polygon (depicted by the orange dotted lines) using Green’s theorem.
To use Green’s theorem, we compute the new polynomials, and . Assuming that and are independent, the = . Thus, the polynomials and can be rewritten as follows:
(2) 
The integral over a line segment with end coordinates () and (), denoted by , can be computed as follows:
(3) 
Substituting the expressions for and presented in subsection 3.1 into Equation 3, we have the following:
(4) 
For each vertex () of a polygon, the line integral (subsection 3.1) is computed for the line segment that connects coordinates () and () in a counterclockwise direction (refer to the arrow heads of the polygon with dotted orange lines in Figure 3), where indicates the number of vertices of the intersection polygon. The expressions of integration of the uniform, Epanechnikov, and Gaussian distribution functions (computed with Wolfram alpha [22]) are provided in Table 1. Finally, the interior probability, Pr( = ( ), can be computed using Green’s theorem (i.e., by summing the line integrals [subsection 3.1] computed for each edge of the intersection polygon). Please refer to Sect. 1 of the supplementary material for a numerical integration experiment validating our derivation in subsection 3.1.
3.2 Monte Carlo Integration
In subsection 3.1, we derived interior probabilities in closed form for parametric distributions with the uniform, Epanechnikov, and Gaussian kernels. Such analytical derivations are possible since a closedform integration exists for each of the three kernels (see Table 1). For parametric distributions that are difficult to integrate or do not have a closedform integration, a more generic Monte Carlo integration can be employed to obtain an approximate solution. In the Monte Carlo approach for the independent noise assumption, samples are independently drawn from distributions and at a vertex. If samples lie in the interior of trait , we quantify or estimate the interior probability as Pr( ) = . This estimation technique is similar to the probabilistic marching cubes [43] technique, which estimates the spatial probability for level sets via Monte Carlo sampling of probability distributions. The accuracy of Monte Carlo solutions increases and converges with an increase in the number of samples.
Kernel 


Uniform  
Epanechnikov  
Gaussian  

4 Fiber Uncertainty for Nonparametric Noise Models
4.1 Independent Noise Assumption
We will first discuss a closedform integration approach for fiber uncertainty quantification. For our closedform approach, we utilize the analytical derivations of independent parametric noise models proposed in subsection 3.1 as building blocks for deriving polynomial integration of nonparametric noise models. In nonparametric density models, the probability distribution of each uncertain variable per vertex can be estimated by deriving a histogram or by using KDE [41]. Let denote independently drawn samples from an unknown probability distribution of a random variable . The can then be estimated from samples with KDE as follows:
(5) 
The in Equation 5 denotes a nonnegative bandwidth of a kernel associated with a sample. represents a scaled kernel, where . The functions in Table 1 denote kernels with a unit bandwidth centered at . We estimate the bandwidth of a kernel from samples using Silverman’s rule of thumb [57]. Similarly, for samples independently drawn from an unknown probability distribution , the probability distribution for bandwidth can be estimated as follows:
(6) 
Substituting Equation 5 and Equation 6 in subsection 3.1, the line integral for the independent nonparametric density estimation can be written as
(7) 
Rearranging the order of summations in subsection 4.1 gives us
(8) 
subsection 4.1 takes a form similar to the one for the parametric statistics presented in subsection 3.1, except that it loops through all pairs of kernels for random variables and . Kernel in subsection 4.1 can be replaced with the uniform, Epanechnikov, or Gaussian functions presented in Table 1.
Histograms could be considered as a special case of KDE with nonoverlapping uniform kernels associated with fixedbin centers (precomputed from noise samples and userspecified number of bins) and with a bandwidth equivalent to the bin width. Thus, for histograms, looping through all pairs of bins of histograms of random variables and is required for computing the line integral in subsection 4.1. Finally, the interior probability (i.e., Pr( = )) can be computed using Green’s theorem by summing line integrals for all edges of the intersection polygon, as illustrated in Figure 3.
When kernel of a nonparametric distribution is difficult to integrate or nonintegrable, a Monte Carlo approach similar to the one for independent parametric distributions (subsection 3.2) can be employed. In other words, for samples independently drawn from nonparametric distributions and with the base kernel , the interior probability can be estimated as Pr( ) = if samples lie in the interior of trait .
4.2 Correlated Noise Assumption
When random variables and are correlated, is not equal to the product . In the correlated noise case, 2D histograms or bivariate KDE can be used to capture the correlation between variables and and the multimodality of their probability distributions. Bivariate KDE extends the idea of univariate KDE in Equation 5 to two variables [62]. An example of 2D histograms and bivariate KDE is provided in Sect. 2 of the supplementary material.
For a 2D histogram, the interior probability (i.e., Pr( = )) can be computed in closed form. Specifically, we loop through each 2D bin of a 2D histogram and quantify the amount of bin overlap with trait by taking into account the bin extent. We then scale the quantified overlap for a 2D bin with the weight of a bin computed from a 2D histogram. Our approach here is similar to the histogram approach for the independent nonparametric case (subsection 4.1), except that the weight of each bin is not a product , but the weight computed based on a 2D histogram.
In the case of a bivariate KDE, we propose a numerical integration approach similar to the one presented in Sect. 1 of the supplementary material. Specifically, we discretize the uncertain ranges of random variables and with a uniform grid. We then estimate the density at each grid vertex using a bivariate KDE. We use the gaussian_kde function from the Python SciPy package [61] for bivariate KDE. We then perform a pointinpolygon test for userspecified FSCP at each grid vertex and sum the weights of all vertices that lie inside FSCP or trait to compute the interior probability (i.e., Pr( = )).
5 Memory and Computational Complexity
Here, we discuss the memory and computational complexity of our proposed parametric and nonparametric statistical models. Let the size of uncertain input data be , where denotes the dimension of a uniform grid that represents domain , indicates the number of variables ( for bivariate data), and is the number of ensemble members or noise samples that represent uncertainty in the data. For closedform derivations with the independent parametric noise models in subsection 3.1, we summarize
noise samples with two parameters per grid vertex. For example, we store the mean and width per vertex for a uniform distribution and the mean and standard deviation for a Gaussian distribution. Thus, for parametric noise models, the memory consumed is
, which amounts to data reduction by a factor of . To compute the interior probability volume for parametric noise models, one polynomial integration is performed per trait per grid vertex if uncertain data at a vertex intersect trait . Thus, the worstcase computational complexity is proportional to the grid size (i.e., ).For closedform derivations with the independent nonparametric statistical approach (subsection 4.1), the probability density per vertex can be estimated with KDE or histograms. For KDE, if a kernel is associated with each of the noise samples per vertex and per variable, then the memory requirement is the same as the original data (i.e., ). If nonparametric density is estimated with a histogram with bins per vertex and per variable, then memory consumption is reduced to , which amounts to data reduction by a factor of . To compute the interior probability volume, the number of computations per vertex increases quadratically with the number of kernels or histogram bins owing to the nested summation in subsection 4.1 (see the experiment in Figure 6).
For the more generic Monte Carlo approach described in subsection 3.2 and subsection 4.1, the number of computations grows in proportion to the grid size and number of samples drawn per grid vertex (i.e., ). In the case of the correlated nonparametric statistical approach (subsection 4.2), the space complexity is again equivalent to that of the independent nonparametric approach. For 2D histograms, the number of computations grows quadratically with the number of bins because we loop through each 2D bin. For bivariate KDE, we perform numerical integration by discretizing the uncertain data range at a vertex with a uniform grid of resolution, . Thus, the worstcase time complexity is proportional to per trait if the uncertain data range at each grid vertex intersects trait .
6 Visualization of Uncertain Fibers
We visualize the positional uncertainty of fibers, i.e., the interior probability volumes derived using our proposed parametric (section 3) and nonparametric (section 4) statistics, via most probable fiber surface extraction, direct volume rendering, and probabilistic segmentation. Athawale et al. proposed a vertexbased classification framework [4] for extracting the most probable isosurface from uncertain scalar fields. We extend their idea for extracting the most probable fiber surface from uncertain bivariate data. Specifically, we derive the interior probability Pr( = ) per vertex, as described in section 3 and section 4. We then probabilistically predict the vertex sign as if Pr( ; otherwise, we predict the vertex sign to be . The isosurface extracted from a grid with this sign classification indicates the most probable fiber surface topology.
We visualize the positional uncertainty of fibers via direct volume rendering of interior probability volumes derived using parametric and nonparametric statistics, similar to the visualizations by Zheng and Sadlo [69]. We render the positions with relatively high interior probability using higher opacity and the ones with relatively low interior probability using lower opacity. We visualize the probabilistic segmentation by thresholding computed interior probabilities. Specifically, for threshold , we visualize the positions with Pr( . Such visualizations provide probabilistic insight into fiber positions.
7 Results and Discussions
We validate our techniques and present their effectiveness via experiments on synthetic and simulation datasets.
7.1 Synthetic Data
7.1.1 Independent Noise Models
We present the results of our proposed closedform independent parametric (section 3) and nonparametric (subsection 4.1) statistical models for the synthetic data in Figure 4. For our experiment, we sample the synthetic sphere and tangle [34] functions on a 3D grid with resolution . Figure 4a visualizes a 2D continuous scatterplot of sphere and tangle values plotted on the horizontal and vertical axes, respectively, for each grid vertex. The cyan polygon in Figure 4a denotes FSCP or trait , and the corresponding fiber surface is visualized in Figure 4b, which we treat as the reference.
We mix the synthetic data sampled on a grid with noise to generate an ensemble of
members. Specifically, we draw samples from a bimodal probability distribution, in which the mode with 80% cumulative probability density is centered around the ground truth, and the mode with 20% cumulative probability density is situated relatively far away from the ground truth. Thus, the noise samples from the mode with 20% cumulative probability density denote the outlier samples. We study the uncertainty in fiber positions arising from the ensemble.
Figure 4c visualizes the fiber surface for the meanfield. The meanfield fiber surface is colormapped with respect to its signed distance [6] from the ground truth fiber surface. Figs. 4d–4f visualize the results for the independent Gaussian noise assumption. Figure 4d visualizes the interior probability volume (i.e., Pr( )) per grid vertex via direct volume rendering, in which yellow with higher opacity denotes positions with relatively high interior probabilities, and white with lower opacity denotes positions with relatively low interior probabilities. Figure 4e visualizes the most probable fiber surface computed using the vertexbased classification (section 6) colormapped with the signed distance from the ground truth fiber surface. Figure 4f visualizes the probabilistic segmentation (section 6) extracted for isovalues , , and , and mapped to high, moderate, and low opacity, respectively. The probabilistic segmentation provides insight into the probability of fiber positions being in the interior of a fiber surface for different segments. Figs. 4g–4i visualize results similar to Figs. 4d–4f for the independent nonparametric KDE with the Gaussian base kernel.
The most probable fiber surface for the nonparametric noise assumption (Figure 4h) is spatially closer to the ground truth fiber surface (Figure 4b) when compared with the meanfield fiber surface (Figure 4c) and the most probable fiber surface for the parametric assumption (Figure 4e), as is evident by the colormapped signed distance. The result implies a greater accuracy of interior probability computations per grid vertex with the nonparametric approach than with the meanfield or parametric approach. This result is expected because the nonparametric models can capture the bimodality of the underlying noise distribution more accurately than the parametric models, which results in nonparametric models that are less sensitive to outliers than the parametric models.
We further confirm the enhanced accuracy of nonparametric models compared to meanfield and parametric models using a quantitative error analysis, as shown in Figure 5. Specifically, we compute and plot a Euclidean 2norm of the difference between the interior probability volume (i.e., Pr( ) per vertex) computed using different statistical techniques and the interior probability volume for the ground truth data, which we refer to as interior probability error. In the case of the ground truth data, the interior probability is when data at a vertex lie inside trait ; otherwise, it is . Figure 5a visualizes the results for parametric noise models with uniform (red), Epanechnikov (green), and Gaussian (blue) assumptions. Figure 5b visualizes the results for nonparametric KDE with uniform (red), Epanechnikov (green), and Gaussian (blue) base kernels. The dotted magenta line in Figure 5a represents the error for the meanfield. In Figure 5, we observe that the error computed for the nonparametric approach (Figure 5b) is consistently lower than the error computed for the parametric density assumption and meanfield (Figure 5a). The Gaussian kernel for the parametric and nonparametric statistics has the highest computational accuracy when compared with the uniform and Epanechnikov kernels.
We compare the results of the Monte Carlo sampling technique (subsection 3.2 and subsection 4.1) with our closedform solutions in Figure 5. In Figure 5, the dashed lines denote the fixed error corresponding to our proposed closedform solutions, and the solid curves plot the error as a function of the number of Monte Carlo samples. The error computed for the Monte Carlo solutions (solid curves) converges to the one computed for our proposed analytical solutions (dashed lines) as we increase the number of Monte Carlo samples. Such convergence confirms the correctness of our derivations (Equations 3.1 and 4.1) and implementations.
Figure 6 visualizes the accuracy and time curves as a result of histogram rebinning for the independent nonparametric noise assumption. At each grid vertex of the domain, we perform nonparametric density estimation by computing a histogram of ensemble members followed by the computation of interior probability. We then compute and plot the interior probability error (i.e., the Euclidean 2norm of the difference between the interior probability volumes for the histogram noise model and the ground truth data). The number of bins determines the memory and computational requirements of the proposed nonparametric approach (section 5) and impacts the accuracy. As illustrated in Figure 6a, the error reduces as the number of histogram bins increases. Note the sharp decline in errors for a relatively small number of bins. On the other hand, the computational complexity appears to grow quadratically, as we had anticipated. Since the interior probability at each grid vertex can be computed independently, the computations can be parallelized. As illustrated in Figure 6b, the red and magenta curves indicate the speedup achieved with the OpenMP C++ parallel implementation with four and threads, respectively. The GPU CUDA version of our code achieved 23X average speedup with the NVIDIA V100 graphics card (the green curve in Figure 6b) compared to the Power9 CPU with openMP threads. The computing resources are courtesy of the Summit Supercomputer at the Oak Ridge National Laboratory.
7.1.2 Correlated Noise Models
We study the correlated noise models (subsection 4.2) via a synthetic experiment. The experimental settings are similar to the ones for the independent noise, except the bimodal noise now comprises two correlated Gaussians (e.g., Fig. 2a in the supplementary material) with the correlation specified by covariance matrices. A correlated Gaussian with an 80% probability concentration is situated close to the ground truth, and one with 20% probability concentration is situated relatively far away from the ground truth (denoting the outliers). The noise samples drawn from such bimodal distributions are mixed with the ground truth to generate an ensemble.
We observe that the ability of bivariate KDE or 2D histograms to reliably capture the correlation between random variables and at each grid vertex is influenced by the number of ensemble members and by the grid resolution used for numerical integration of bivariate KDE. When the number of noise samples is relatively large (e.g., samples for bivariate KDE estimation in Fig. 2c of the supplementary material), the correlation is captured reliably. For fewer samples, the correlation might not be captured well and might lead to relatively large errors. As shown in Figure 7a, the interior probability error (blue curve) reduces with an increase in the number of ensemble members. The interior probabilities computed with the independent nonparametric noise assumption (the green dotted line in Figure 7a) appear to be more accurate than the correlated noise assumption when the number of ensemble members is .
The error for nonparametric models with the correlated or independent noise assumptions (Figure 7a) is again smaller than the errors for the meanfield and parametric noise models owing to the higher robustness of nonparametric models to outliers. The meanfield for ensemble members resulted in an interior probability error equal to . The parametric Gaussian noise assumption resulted in an interior probability error equal to . The error for the meanfield and parametric noise assumption does not fluctuate much with an increase in the number of ensemble members. Figure 7b visualizes the error as a function of a grid resolution used for numerical integration in the case of bivariate KDE. The error reduces slowly with an increase in grid resolution.
7.2 Rectangular Vs. Arbitrary Shape of Polygonal Traits
The uncertainty visualization framework proposed by Zheng and Sadlo [69] is limited to a rectangular trait (subsection 1.2). Our proposed framework extends their work to a polygonal trait with arbitrary shape. Figure 8 demonstrates the advantage of an arbitrary polygon over rectangle for trait selection in the contexts of the original as well as uncertain data. Figure 8a visualizes a continuous scatterplot of the electron density (Rho) and reduced gradient (s) attributes of the ethanediol molecule. The main separating axis or diagonal of the continuous scatterplot (with strong red hue in Figure 8a) indicates regions with no chemical interactions. Thus, regions away from this axis are important for chemists in analyzing chemical interactions. The two traits and are selected in the continuous scatterplot space to enclose a nondiagonal feature. As observed in Figure 8a, the arbitrary shape of a polygon allows us to select a trait ( in magenta) that is closer to a nondiagonal feature than the rectangular shape of a trait ( in cyan).
The fiber surfaces corresponding to traits and are visualized for the original data in Figure 8b. The fiber surface for the nonrectangular trait (the left image in Figure 8b) segments out carbon atoms more accurately than the fiber surface for the rectangular trait (the right image in Figure 8b). Next, we reduce the original data with the Gaussiandistributed hixel representation [58]. Specifically, we partition the original data into blocks of size and represent each block as a Gaussian distribution with mean and standard deviation. The fiber surfaces corresponding to traits and are visualized for the hixel representation in Figure 8c. The carbon atoms are again segmented well for trait with the high probability (yellow) regions (the left image of Figure 8c). The two newly formed blobs in the left image of Figure 8c that do not exist in the left image of Figure 8b are due to the uncertainty in data arising from the hixel representation, but they have a very low probability (purple) of existence. As observed in the right image of Figure 8c, the rectangular trait results in an overestimation of features, i.e., a multiple high probability (yellow) regions, that do not quite correctly represent carbon atoms.
7.3 Analysis of Simulation Datasets with Parametric Vs. Nonparametric Models of Uncertainty
We demonstrate an application of our proposed parametric and nonparametric statistical techniques (section 3 and section 4) in the context of oceanology for the analysis of ocean eddies. The knowledge of eddy positions is important for oceanologists to understand the transport of energy and biogeochemical particles in oceans. Figure 9 visualizes the results of our probabilistic analysis of fibers that represent vortical features of the Red Sea dataset [53]. The dataset was downloaded from the IEEE SciVis Contest 2020 website (https://kaustvislab.github.io/SciVis2020/) and comprises ensembles of multiple variables pertinent to the oceanology simulated for time steps. In this experiment, we visualize the positional uncertainty of fibers corresponding to eddy features of the Red Sea by analyzing ensemble members of the velocity field for a time step of over the Gulf of Aden sampled on a grid of resolution .
Initially, we derive the vorticity field for each ensemble member by computing the curl of the velocity field. We then compute the vorticity magnitude and Z (vertical) component of a vorticity vector per grid vertex for each member. Figure 9a visualizes a continuous scatterplot with the mean vorticity magnitude plotted on the vertical axis and the mean Z component of vorticity vector plotted on the horizontal axis. Trait (cyan polygon) is selected to understand the fiber positions with relatively high vorticity and a negative value of the Z component of a vorticity vector which denote anticyclonic eddy positions.
Figure 9b visualizes the meanfield fiber surface corresponding to trait . The magenta box in the meanfield visualization indicates a rare occurrence of vortical features. Figure 9c visualizes the results for the independent Gaussian (parametric) noise model (the same assumption as in the paper by Zheng and Sadlo [69]), but for nonrectangular polygonal shape (our contributions). The volume rendering for the Gaussian noise model (top image in Figure 9c) indicates the presence of a relatively large number of vortical features inside the region enclosed by the magenta box because there are multiple regions with high interior probability (mapped to yellow). The probabilistic segmentation of the interior probability volume with an isosurface for isovalue (the bottom image in Figure 9c) clearly shows the presence of multiple vortical features inside the region enclosed by the magenta box. Figure 9d visualizes results similar to Figure 9c, but for the independent histrogram noise model with bins. Both volume rendering and probabilistic segmentation for the nonparametric noise model indicate a relatively low probability of vortical features inside the region enclosed by the magenta box. Such inconsistency regarding the presence of vortical features inside the magenta box across the three statistical models indicates the need for further eddy analysis in the same region.
We parallelized the implementation of the independent nonparametric histogram models with OpenMP C++ and CUDA using an approach similar to the one for the histogram results in Figure 6. The serial computation of interior probabilities required seconds, whereas computations with four and threads resulted in a significant speedup with and seconds, respectively, of computational time. The GPU CUDA computations required only seconds.
Figure 10 visualizes the results of our proposed statistical techniques for the asteroid impact dataset [42]. For our experiment, we analyze time step of a simulation, in which an asteroid m in diameter bursts at an elevation of 5 km above sea level and impacts the sea at an angle of 45 (simulation series YB31 on the IEEE SciVis 2018 contest website (https://sciviscontest2018.org/). We treat the dataset with resolution as the reference. Figure 10a visualizes a continuous scatterplot of the temperature and volume fraction of water variables for the reference dataset. Figure 10b visualizes the ground truth fiber surface extracted from the reference volume for trait , as indicated by the dark blue polygon in Figure 10a. Trait denotes the positions with a relatively high volume fraction of water and temperature. Next, we partition the data into bricks and summarize each brick with a probability distribution similar to the hixel technique proposed by Thompson et al. [58].
Figure 10c visualizes the meanfield, which results in a significant loss of features caused by data variations in each brick. Figs. 10d10f visualize the interior probabilities derived assuming the independent Epanechnikov (parametric) distribution (section 3), independent histograms with two bins per variable (subsection 4.1), and correlated histograms with two bins per variable (subsection 4.2), respectively. Note that for the correlated 2D histograms in Figure 10f, we compute interior probabilities in closed form (without requiring numerical integration), as described in subsection 4.2. The probabilistic feature recovery with our proposed statistical frameworks is quite remarkable compared to the meanfield features, as illustrated by the boxes in Figure 10, with Figure 10b as the reference. The nonparametric statistics in Figs. 10e–10f show a slightly improved recovery of probabilistic features than the parametric statistics in Fig. 10d, as anticipated based on our synthetic data experiments (e.g., Figure 5). A relatively small number of noise samples (here, eight per vertex) may not effectively capture noise correlation in the case of correlated nonparamertic statistics (e.g., Figure 7), thus yielding quite similar results in Figs. 10e–10f.
8 Conclusion and Future Work
We extend the prior work by Zheng and Sadlo [69] to study the positional uncertainty of fibers for uncertain bivariate data. Specifically, we present a closedform statistical framework for spatial uncertainty quantification of fibers when data noise is characterized by independent parametric and nonparametric probability distributions. We perform our analysis for the uniform, Epanechnikov, and Gaussian kernels. Additionally, we present a numerical integration approach for uncertainty quantification of fibers when noise in bivariate data is correlated. We present our statistical frameworks for arbitrary shapes of FSCP as opposed to being limited to rectangular FSCP [69, 52]. We show that the nonparametric statistics improve the accuracy of uncertainty quantification and visualization when compared to the parametric and meanfield statistics via experiments on synthetic and simulation datasets.
In the future, we would like to extend our research to multivariate data with more than two variables. For bivariate data with correlated noise assumption, we would like to derive closedform solutions (instead of numerical integration) similar to the independent noise assumption, if computable. Finally, we are interested in investigating the effect of uncertainty in multivariate data on the topology of fiber surfaces.
Acknowledgements.
This work was partially supported by the Scientific Discovery through Advanced Computing (SciDAC) program in the U.S. Department of Energy, the Intel Graphics and Visualization Institutes of XeLLENCE, the Intel OneAPI CoE, the NIH under award number R24 GM136986, the DOE under grant number DEFE0031880, and the Utah Office of Energy Development. We wish to thank Dr. Jieyang Chen at the Oak Ridge National Laboratory for helping us with the GPU code implementation of fiber uncertainty quantification. We would also like to thank the reviewers of this article for their valuable feedback.References
 [1] T. M. Athawale and C. R. Johnson. Probabilistic asymptotic decider for topological ambiguity resolution in levelset extraction for uncertain 2D data. IEEE Transactions on Visualization and Computer Graphics, 25(1):1163–1172, Jan. 2019. doi: 10 . 1109/TVCG . 2018 . 2864505
 [2] T. M. Athawale, B. Ma, E. Sakhaee, C. R. Johnson, and A. Entezari. Direct volume rendering with nonparametric models of uncertainty. IEEE Transactions on Visualization and Computer Graphics, 27(2):1797–1807, Feb. 2021. doi: 10 . 1109/TVCG . 2020 . 3030394
 [3] T. M. Athawale, D. Maljovec, L. Yan, C. R. Johnson, V. Pascucci, and B. Wang. Uncertainty visualization of 2D Morse complex ensembles using statistical summary maps. IEEE Transactions on Visualization and Computer Graphics, 28(4):1955–1966, Apr. 2022. doi: 10 . 1109/TVCG . 2020 . 3022359
 [4] T. M. Athawale, E. Sakhaee, and A. Entezari. Isosurface visualization of data with nonparametric models for uncertainty. IEEE Transactions on Visualization and Computer Graphics, 19(12):2723–2732, Jan. 2016. doi: 10 . 1109/TVCG . 2015 . 2467958
 [5] S. Bachthaler and D. Weiskopf. Continuous scatterplots. IEEE Transactions on Visualization and Computer Graphics, 14(6):1428–1435, Nov.Dec. 2008. doi: 10 . 1109/TVCG . 2008 . 119
 [6] J. Baerentzen and H. Aanaes. Signed distance computation using the angle weighted pseudonormal. IEEE Transactions on Visualization and Computer Graphics, 11(3):243–253, MayJune 2005. doi: 10 . 1109/TVCG . 2005 . 49
 [7] G. Bonneau, H. Hege, C. R. Johnson, M. Oliveira, K. Potter, P. Rheingans, and T. Schultz. Overview and stateoftheart of uncertainty visualization. In M. Chen, H. Hagen, C. Hansen, C. R. Johnson, and A. Kauffman, eds., Scientific Visualization: Uncertainty, Multifield, Biomedical, and Scalable Visualization, pp. 3–27. Springer London, 2014. doi: 10 . 1007/9781447164975_1
 [8] K. Brodlie, R. A. Osorio, and A. Lopes. A review of uncertainty in data visualization. In J. Dill, R. Earnshaw, D. Kasik, J. Vince, and P. C. Wong, eds., Expanding the Frontiers of Visual Analytics and Visualization, pp. 81–109. Springer Verlag London, 2012.
 [9] H. Carr, Z. Geng, J. Tierny, A. Chattopadhyay, and A. Knoll. Fiber surfaces: Generalizing isosurfaces to bivariate data. Computer Graphics Forum, 34(3):241–250, July 2015. doi: 10 . 1111/cgf . 12636
 [10] H. Carr, J. Snoeyink, and M. van de Panne. Flexible isosurfaces: Simplifying and displaying scalar topology using the contour tree. Computational Geometry, 43(1):42–58, Jan. 2010. Special Issue on the 14th Annual Fall Workshop. doi: 10 . 1016/j . comgeo . 2006 . 05 . 009
 [11] H. Edelsbrunner, J. Harer, and A. K. Patel. Reeb spaces of piecewise linear mappings. In Proceedings of the TwentyFourth Annual Symposium on Computational Geometry, SCG ’08, p. 242–250. Association for Computing Machinery, New York, NY, USA, June 2008. doi: 10 . 1145/1377676 . 1377720
 [12] G. Favelier, N. Faraj, B. Summa, and J. Tierny. Persistence atlas for critical point variability in ensembles. IEEE Transactions on Visualization and Computer Graphics, 25(1):1152–1162, Jan. 2019. doi: 10 . 1109/TVCG . 2018 . 2864432
 [13] F. Ferstl, K. Bürger, and R. Westermann. Streamline variability plots for characterizing the uncertainty in vector field ensembles. IEEE Transactions on Visualization and Computer Graphics, 22(1):767–776, Jan. 2016. doi: 10 . 1109/TVCG . 2015 . 2467204
 [14] N. Fout and K.L. Ma. Fuzzy volume rendering. IEEE Transactions on Visualization and Computer Graphics, 18(12):2335–2344, Dec. 2012. doi: 10 . 1109/TVCG . 2012 . 227
 [15] M. Friendly and D. Denis. The early origins and development of the scatterplot. Journal of the History of the Behavioral Sciences, 41(2):103–130, Spring 2005. doi: 10 . 1002/jhbs . 20078
 [16] R. Fuchs and H. Hauser. Visualization of multivariate scientific data. Computer Graphics Forum, 28(6):1670–1690, Sept. 2009. doi: 10 . 1111/j . 14678659 . 2009 . 01429 . x
 [17] D. Günther, J. Salmon, and J. Tierny. Mandatory critical points of 2D uncertain scalar fields. Computer Graphics Forum, 33(3):31–40, July 2014. doi: 10 . 1111/cgf . 12359
 [18] H. Guo, W. He, T. Peterka, H.W. Shen, S. M. Collis, and J. J. Helmus. Finitetime Lyapunov exponents and Lagrangian coherent structures in uncertain unsteady flows. IEEE Transactions on Visualization and Computer Graphics, 22(6):1672–1682, June 2016. doi: 10 . 1109/TVCG . 2016 . 2534560
 [19] H. Hauser, F. Ledermann, and H. Doleisch. Angular brushing of extended parallel coordinates. In IEEE Symposium on Information Visualization, 2002. INFOVIS 2002., pp. 127–130, Oct. 2002. doi: 10 . 1109/INFVIS . 2002 . 1173157
 [20] S. Hazarika, A. Biswas, and H.W. Shen. Uncertainty visualization using copulabased analysis in mixed distribution models. IEEE Transactions on Visualization and Computer Graphics, 24(1):934–943, Jan. 2018. doi: 10 . 1109/TVCG . 2017 . 2744099
 [21] J. Heinrich and D. Weiskopf. Continuous parallel coordinates. IEEE Transactions on Visualization and Computer Graphics, 15(6):1531–1538, Nov.Dec. 2009. doi: 10 . 1109/TVCG . 2009 . 131
 [22] W. R. Inc. Mathematica, Version 13.0.0. Champaign, IL, 2021.
 [23] A. Inselberg and B. Dimsdale. Parallel coordinates. In A. Klinger, ed., HumanMachine Interactive Systems, pp. 199–233. Springer US, Boston, MA, 1991. doi: 10 . 1007/9781468458831_9
 [24] H. Jänicke, M. Böttinger, and G. Scheuermann. Brushing of attribute clouds for the visualization of multivariate data. IEEE Transactions on Visualization and Computer Graphics, 14(6):1459–1466, Nov.Dec. 2008. doi: 10 . 1109/TVCG . 2008 . 116
 [25] J. Jankowai and I. Hotz. Feature levelsets: Generalizing isosurfaces to multivariate data. IEEE Transactions on Visualization and Computer Graphics, 26(2):1308–1319, Feb. 2020. doi: 10 . 1109/TVCG . 2018 . 2867488
 [26] F. Jiao, J. M. Phillips, Y. Gur, and C. R. Johnson. Uncertainty visualization in HARDI based on ensembles of ODFs. In 2012 IEEE Pacific Visualization Symposium, pp. 193–200, Feb.Mar. 2012. doi: 10 . 1109/PacificVis . 2012 . 6183591
 [27] C. R. Johnson and A. R. Sanderson. A next step: Visualizing errors and uncertainty. IEEE Computer Graphics and Applications, 23(5):6–10, Sept.Oct. 2003. doi: 10 . 1109/MCG . 2003 . 1231171
 [28] E. R. Johnson, S. Keinan, P. MoriSánchez, J. ContrerasGarcía, A. J. Cohen, and W. Yang. Revealing noncovalent interactions. Journal of the American Chemical Society, 132(18):6498–6506, Apr. 2010. PMID: 20394428. doi: 10 . 1021/ja100936w
 [29] D. K. Jones. Determining and visualizing uncertainty in estimates of fiber orientation from diffusion tensor MRI. Magnetic Resonance in Medicine, 49(1):7–12, Dec. 2002. doi: 10 . 1002/mrm . 10331
 [30] A. Kamal, P. Dhakal, A. Javaid, V. K. Devabhaktuni, D. Kaur, J. Zaientz, and R. Marinier. Recent advances and challenges in uncertainty visualization: a survey. Journal of Visualization, 24(5):861–890, May 2021. doi: 10 . 1007/s12650021007551
 [31] J. Kehrer and H. Hauser. Visualization and visual analysis of multifaceted scientific data: A survey. IEEE Transactions on Visualization and Computer Graphics, 19(3):495–513, Mar. 2013. doi: 10 . 1109/TVCG . 2012 . 110
 [32] P. Klacansky, J. Tierny, H. Carr, and Z. Geng. Fast and exact fiber surfaces for tetrahedral meshes. IEEE Transactions on Visualization and Computer Graphics, 23(7):1782–1795, July 2017. doi: 10 . 1109/TVCG . 2016 . 2570215
 [33] J. Kniss, G. Kindlmann, and C. Hansen. Multidimensional transfer functions for interactive volume rendering. IEEE Transactions on Visualization and Computer Graphics, 8(3):270–285, Jul.Sept. 2002. doi: 10 . 1109/TVCG . 2002 . 1021579
 [34] A. Knoll, Y. Hijazi, A. Kensler, M. Schott, C. Hansen, and H. Hagen. Fast ray tracing of arbitrary implicit surfaces with interval and affine arithmetic. Computer Graphics Forum, 28(1):26–40, Feb. 2009. doi: 10 . 1111/j . 14678659 . 2008 . 01189 . x
 [35] S. Liu, J. A. Levine, P.T. Bremer, and V. Pascucci. Gaussian mixture model based volume visualization. In IEEE Symposium on Large Data Analysis and Visualization (LDAV), pp. 73–77, Oct. 2012. doi: 10 . 1109/LDAV . 2012 . 6378978
 [36] S. Lodha, A. Pang, R. Sheehan, and C. Wittenbrink. Uflow: visualizing uncertainty in fluid flow. In Proceedings of Seventh Annual IEEE Visualization ’96, pp. 249–254, Oct.Nov. 1996. doi: 10 . 1109/VISUAL . 1996 . 568116
 [37] W. E. Lorensen and H. E. Cline. Marching cubes: A high resolution 3D surface construction algorithm. SIGGRAPH Computer Graphics, 21(4):163–169, Aug. 1987. doi: 10 . 1145/37402 . 37422
 [38] C. Lundström, P. Ljung, A. Persson, and A. Ynnerman. Uncertainty visualization in medical volume rendering using probabilistic animation. IEEE Transactions on Visualization and Computer Graphics, 13(6):1648–1655, Nov.Dec. 2007. doi: 10 . 1109/TVCG . 2007 . 70518
 [39] S. Nagaraj, V. Natarajan, and R. S. Nanjundiah. A gradientbased comparison measure for visual analysis of multifield data. In Proceedings of the 13th Eurographics / IEEE  VGTC Conference on Visualization, EuroVis’11, p. 1101–1110. The Eurographs Association & John Wiley & Sons, Ltd., Chichester, GBR, June 2011. doi: 10 . 1111/j . 14678659 . 2011 . 01959 . x
 [40] M. Otto, T. Germer, and H. Theisel. Uncertain topology of 3D vector fields. In 2011 IEEE Pacific Visualization Symposium, pp. 67–74, Mar. 2011. doi: 10 . 1109/PACIFICVIS . 2011 . 5742374
 [41] E. Parzen. On estimation of a probability density function and mode. The Annals of Mathematical Statistics, 33(3):1065–1076, Sept. 1962. doi: 10 . 1214/aoms/1177704472
 [42] J. M. Patchett, F. Samsel, K. Tsai, G. R. Gisler, D. H. Rogers, G. Abram, T. L. Turton, and L. Alamos. Visualization and analysis of threats from asteroid ocean impacts. Technical report, Nov. 2016.
 [43] K. Pöthkow, B. Weber, and H.C. Hege. Probabilistic marching cubes. Computer Graphics Forum, 30(3):931–940, June 2011. doi: 10 . 1111/j . 14678659 . 2011 . 01942 . x
 [44] K. Potter, P. Rosen, and C. R. Johnson. From quantification to visualization: A taxonomy of uncertainty visualization approaches. In A. M. Dienstfrey and R. F. Boisvert, eds., Uncertainty Quantification in Scientific Computing, pp. 226–249. Springer Berlin Heidelberg, Berlin, Heidelberg, 2012. doi: 10 . 1007/9783642326776_15
 [45] K. Pöthkow and H.C. Hege. Nonparametric models for uncertainty visualization. Computer Graphics Forum, 32(3pt2):131–140, July 2013. doi: 10 . 1111/cgf . 12100
 [46] G. Quadri and P. Rosen. Modeling the influence of visual density on cluster perception in scatterplots using topology. IEEE Transactions on Visualization & Computer Graphics, 27(02):1829–1839, Feb. 2021. doi: 10 . 1109/TVCG . 2020 . 3030365
 [47] F. Raith, C. Blecha, T. Nagel, F. Parisio, O. Kolditz, F. Günther, M. Stommel, and G. Scheuermann. Tensor field visualization using fiber surfaces of invariant space. IEEE Transactions on Visualization and Computer Graphics, 25(1):1122–1131, Jan. 2019. doi: 10 . 1109/TVCG . 2018 . 2864846
 [48] P. J. Rhodes, R. S. Laramee, R. D. Bergeron, and T. M. Sparr. Uncertainty Visualization Methods in Isosurface Rendering. In Eurographics 2003  Short Presentations. Eurographics Association, 2003. doi: 10 . 2312/egs . 20031054
 [49] O. Saeki. Topology of singular fibers of differentiable maps, vol. 1854 of Lecture Notes in Mathematics. Springer Heidelberg, Germany, 2004. doi: 10 . 1007/b100393
 [50] E. Sakhaee and A. Entezari. A statistical direct volume rendering framework for visualization of uncertain data. IEEE Transactions on Visualization and Computer Graphics, 23(12):2509–2520, Dec. 2017. doi: 10 . 1109/TVCG . 2016 . 2637333
 [51] D. Sakurai, K. Ono, H. Carr, J. Nonaka, and T. Kawanabe. Flexible Fiber Surfaces: A ReebFree Approach, pp. 187–201. Mathematics and Visualization. Springer Science and Business Media Deutschland GmbH, Germany, Dec. 2020. doi: 10 . 1007/9783030430368_12
 [52] S. Sane, T. M. Athawale, and C. R. Johnson. Visualization of Uncertain Multivariate Data via Feature Confidence LevelSets. In M. Agus, C. Garth, and A. Kerren, eds., Proc. EuroVis  Short Papers. The Eurographics Association, 2021. doi: 10 . 2312/evs . 20211053
 [53] S. Sanikommu, H. Toye, P. Zhan, S. Langodan, G. Krokos, O. Knio, and I. Hoteit. Impact of atmospheric and model physics perturbations on a highresolution ensemble data assimilation system of the red sea. Journal of Geophysical Research: Oceans, 125(8):e2019JC015611, July 2020. doi: 10 . 1029/2019JC015611
 [54] N. Sauber, H. Theisel, and H.p. Seidel. Multifieldgraphs: An approach to visualizing correlations in multifield scalar data. IEEE Transactions on Visualization and Computer Graphics, 12(5):917–924, Sept.Oct. 2006. doi: 10 . 1109/TVCG . 2006 . 165
 [55] D. Schneider, J. Fuhrmann, W. Reich, and G. Scheuermann. A variance based FTLElike method for unsteady uncertain vector fields. In R. Peikert, H. Hauser, H. Carr, and R. Fuchs, eds., Topological Methods in Data Analysis and Visualization II: Theory, Algorithms, and Applications, pp. 255–268. Springer Berlin Heidelberg, Berlin, Heidelberg, 2012. doi: 10 . 1007/9783642231759_17
 [56] F. Siddiqui, T. Höllt, and A. Vilanova. A progressive approach for uncertainty visualization in diffusion tensor imaging. Computer Graphics Forum, 40(3):411–422, June 2021. doi: 10 . 1111/cgf . 14317
 [57] B. W. Silverman. Density Estimation for Statistics and Data Analysis. Taylor & Francis, New York, Oct. 1998. doi: 10 . 1201/9781315140919
 [58] D. Thompson, J. A. Levine, J. C. Bennett, P.T. Bremer, A. Gyulassy, V. Pascucci, and P. P. Pébay. Analysis of largescale scalar data using hixels. In 2011 IEEE Symposium on Large Data Analysis and Visualization, pp. 23–30, Oct. 2011. doi: 10 . 1109/LDAV . 2011 . 6092313
 [59] J. Tierny and H. Carr. Jacobi fiber surfaces for bivariate Reeb space computation. IEEE Transactions on Visualization and Computer Graphics, 23(1):960–969, Jan. 2017. doi: 10 . 1109/TVCG . 2016 . 2599017
 [60] J. Tierny, G. Favelier, J. A. Levine, C. Gueunet, and M. Michaux. The topology toolkit. IEEE Transactions on Visualization and Computer Graphics, 24(1):832 – 842, Jan. 2018. doi: 10 . 1109/TVCG . 2017 . 2743938
 [61] P. e. a. Virtanen. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261–272, Feb. 2020. doi: 10 . 1038/s4159201906862
 [62] M. P. Wand and M. C. Jones. Comparison of smoothing parameterizations in bivariate kernel density estimation. Journal of the American Statistical Association, 88(422):520–528, June 1993. doi: 10 . 1080/01621459 . 1993 . 10476303
 [63] M. O. Ward. Multivariate Data Glyphs: Principles and Practice, pp. 179–198. Springer Berlin Heidelberg, Berlin, Heidelberg, 2008. doi: 10 . 1007/9783540330370_8
 [64] R. T. Whitaker, M. Mirzargar, and R. M. Kirby. Contour boxplots: A method for characterizing uncertainty in feature sets from simulation ensembles. IEEE Transactions on Visualization and Computer Graphics, 19(12):2713–2722, Dec. 2013. doi: 10 . 1109/TVCG . 2013 . 143
 [65] K. Wu, A. Knoll, B. J. Isaac, H. Carr, and V. Pascucci. Direct multifield volume ray casting of fiber surfaces. IEEE Transactions on Visualization and Computer Graphics, 23(1):941–949, Jan. 2017. doi: 10 . 1109/TVCG . 2016 . 2599040
 [66] K. Wu and S. Zhang. A contour tree based visualization for exploring data with uncertainty. International Journal for Uncertainty Quantification, 3:203–223, 2012. doi: 10 . 1615/Int . J . UncertaintyQuantification . 2012003956
 [67] Z. Xie, S. Huang, M. O. Ward, and E. A. Rundensteiner. Exploratory visualization of multivariate data with variable quality. In 2006 IEEE Symposium On Visual Analytics Science And Technology, pp. 183–190, Oct.Nov. 2006. doi: 10 . 1109/VAST . 2006 . 261424
 [68] L. Yan, Y. Wang, E. Munch, E. Gasparovic, and B. Wang. A structural average of labeled merge trees for uncertainty visualization. IEEE Transactions on Visualization and Computer Graphics, 26(1):832–842, Jan. 2020. doi: 10 . 1109/TVCG . 2019 . 2934242
 [69] B. Zheng and F. Sadlo. Uncertainty in continuous scatterplots, continuous parallel coordinates, and fibers. IEEE Transactions on Visualization and Computer Graphics, 27(2):1819–1828, Feb. 2021. doi: 10 . 1109/TVCG . 2020 . 3030466
 [70] T. Zuk and S. Carpendale. Visualization of uncertainty and reasoning. In Proceedings of the 8th international symposium on Smart Graphics, pp. 164–177, June 2007. doi: 10 . 1007/9783540732143_15