- Email: [email protected]
Visualization of wavelet compressed mantle convection data
Physics of the Earth and Planetary Interiors 163 (2007) 292–298
Visualization of wavelet compressed mantle convection data Benjamin J. Kadlec a,∗ , Daniel E. Goldstein b , David A. Yuen c , Alexei Vezolainen d a
b
Department of Computer Science, University of Colorado at Boulder, Boulder, CO 80309, United States NorthWest Research Associates, Inc., Colorado Research Associates Division, Boulder, CO 80301, United States c Department of Geology and Geophysics and Minnesota Supercomputing Institute, University of Minnesota, Minneapolis, MN 55455, United States d Department of Mechanical Engineering, University of Colorado at Boulder, Boulder, CO 80309, United States Received 5 February 2007; accepted 20 March 2007
Abstract Wavelet compression methodologies promise significant computational and storage cost savings through the benefit of producing results on adapted grids that significantly reduce storage and data manipulation costs. Yet, with these storage methodologies come new challenges in the visualization of temporally adaptive datasets. Wavelet compressed three-dimensional mantle convection simulations are visualized as a challenging case study. We employ adaptive 3D visualization using cache friendly octal tree volume visualization as an effective visualization approach. To explore the localized multi-scale structures in the datasets, the wavelet coefficients are also visualized allowing visualization of energy contained in local physical regions as well as in local wavenumber space. Our technique allows for exploring localized multi-scale structures using multi-resolution visualization from the compressed state of large geophysical simulations. © 2007 Elsevier B.V. All rights reserved. Keywords: 3D visualization; Wavelets; Mantle convection; Volume visualization
1. Introduction The race towards petascale computing has created an explosion in data from large-scale computations of nonlinear phenomena in many fields, including the geo- and environmental sciences. Today, the geosciences are facing an onslaught of information coming from terabyte-scale numerical simulations of nonlinear phenomena and higher resolution data from satellite missions and precise laboratory experiments. In addi∗
Corresponding author. E-mail addresses: [email protected] (B.J. Kadlec), [email protected] (D.E. Goldstein), [email protected] (D.A. Yuen), [email protected] (A. Vezolainen). 0031-9201/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.pepi.2007.03.009
tion, great torrents of data are being produced by new earthquake models coupled with sensing information, which can easily overwhelm researchers due to their high dynamic range and multi-scale character in space and time. As scientists in the geo and environmental disciplines continue their march towards petascale computing, the cost of large-scale data analysis will be a dire problem requiring novel tools to process and represent data more efficiently. The efficient storage and subsequent visualization of these large datasets is a trade off in storage costs versus data quality. The last decade and a half has witnessed the development of wavelets or wavelet analysis. Wavelets have a wide range of application in areas such as feature extraction, data compression, numerical simulation (Goldstein
B.J. Kadlec et al. / Physics of the Earth and Planetary Interiors 163 (2007) 292–298
et al., 2005; Goldstein and Vasilyev, 2004; Schneider et al., 2005), and visualization (Mallat, 1999). Wavelets are an effective tool for efficiently representing data with localized multi-scale structures due to the inherent localized support of the wavelet basis functions in both physical and wavenumber space. This dual-space localization property of wavelets allows for wavelet compressed datasets to represent localized multi-scale structures with a high level of compression and minimal loss of important features. The ability of wavelets to efficiently represent data containing localized multi-scale structure makes them an ideal tool for the simulation of geophysical systems and the subsequent storage and visualization of simulation and experimental results (Yuen et al., 2004; Vecsey et al., 2003; Vezolainen et al., 2005) In this work our emphasis will be on the issues and possible dramatic advantages to be gained by representing existing geo-physical datasets in a wavelet compressed form for storage and visualization. This paper is organized as follows: we will first describe the basic properties of second generation wavelets. Then we will discuss the de-noising properties through wavelet thresholding, which allow decomposition of localized multi-scale coherent structures from an incoherent background. Next we will motivate the need for visualization in our work and describe the visualization techniques we derive. Following this we will present, as an example case, the visualization of a mantle convection simulation dataset at different levels of wavelet compression. We will conclude with some discussion of future improvements that are being studied for the visualization of wavelet compressed datasets.
293
2. Wavelets Wavelets are basis functions which are localized in both physical space (due to their finite support) and wavenumber space (due to their vanishing moments), e.g. Fig. 1. For comparison, the classical Fourier transform is based on functions (sines and cosines) that are well localized in wavenumber space, but do not provide localization in physical space due to their global support. Because of this space/scale localization, the wavelet transform provides both spatial and scale (frequency) information while the Fourier transform on the other hand only provides frequency information. A scalar field f(x) can be represented in terms of wavelet basis functions as f (x) =
+∞ 2 −1 n
cl0 φl0 (x) +
μ,j
μ,j
dk ψk (x), (1)
j=0 μ=1 k ∈ κμ,j
l ∈ L0 μ,j
where φl0 (x) and ψk are respectively n-dimensional scaling functions and wavelets of different families (μ) and levels of resolution (j). One may think of a wavelet decomposition as a multi-level or multi-resolution representation of a function, where each level of resolution j j (except the coarsest one) consists of wavelets ψk or μ,j family of wavelets ψk having the same scale but located at different positions. Scaling function coefficients represent the averaged values of the field, while the wavelet coefficients represent the details of the field at different scales. The wavelet functions have a zero mean, while the scaling functions do not. Note that in n-dimensions there are 2n − 1 distinctive n-dimensional wavelets (Daubechies, 1992). Also note that due to the
Fig. 1. Lifted interpolating wavelet ψ, of order 6 (a) and its Fourier transform (b).
294
B.J. Kadlec et al. / Physics of the Earth and Planetary Interiors 163 (2007) 292–298
local support of both scaling functions and wavelets, there is a one-to-one correspondence between the location of each scaling function or wavelet with a grid point. As a result each scaling function coefficient cl0 and each μ,j wavelet coefficient dk is uniquely associated with a single grid point with the indices l and k, respectively. Traditionally, one-dimensional first generation j wavelets ψk are defined as translates and dilates of j one basic wavelet ψ, i.e. ψk (x) = ψ(2j x − k). Second generation wavelets (Sweldens, 1996, 1998) are a generalization of first generation wavelets that supplies additional freedom to deal with arbitrary boundary conditions, and irregular sampling intervals. Second generation wavelets form a Riesz basis for L2 space, with the wavelets being local in both space and frequency and often having many vanishing polynomial moments, but without the translation and dilation invariance of their first generation cousins. Despite the loss of these two fundamental properties of wavelet bases, second generation wavelets retain many of the useful features of first generation wavelets, including a fast O(N) transform. The construction of second generation wavelets is based on the lifting scheme that is discussed in detail by Sweldens (1996, 1998). For this study we use a set of second generation wavelets known in the literature as lifted interpolating wavelets (Vasilyev and Bowman, 2000; Sweldens, 1996). In particular a lifted interpolating wavelet of order 6, which is shown in Fig. 1 along with its Fourier transform is used in this work. For a more in-depth discussion on the construction of these wavelets the reader is referred to the papers by Sweldens (1996, 1998), and Vasilyev and Bowman (2000). For a more general discussion on wavelets we refer the reader to the books of Daubechies (1992), Mallat (1999), and FoufoulaGeorgiou and Kumar (1994). 2.1. Wavelet filters Wavelet filtering is performed in wavelet space using wavelet coefficient thresholding, which can be considered as a nonlinear filter that depends on each flow realization. The wavelet thresholding filter is defined by f¯ >ε (x) = cl0 φl0 (x) l ∈ L0 +∞ 2 −1 n
+
j=0 μ=1
k ∈ κμ,j
μ,j
|dk |> ∈ ||f ||WTF
μ,j
μ,j
dk ψk (x),
(2)
where f(x) is a scalar field, ε > 0 stands for the nondimensional (relative) threshold parameter, and ||·||WTF being the Wavelet Threshold Filtering (WTF) norm that provides the (absolute) dimensional scaling for filtered variable f. For instance, in the case of velocity, the (absolute) dimensional scaling can be specified as the L2 norm (||ui ||WTF = ||ui ||2 ) or the L∞ norm (||ui ||WTF = ||ui ||∞ ). Note that once the WTF-norm ||·||WTF is specified, the wavelet thresholding filter (2) is uniquely defined by the non-dimensional threshold parameter, ε. The reconstruction error due to wavelet filtering with non-dimensional threshold parameter ε can be shown to be (Vasilyev, 2003; Donoho, 1992): ||f (x) − f¯ >ε (x)||2 ≤ Cε||f ||WTF ,
(3)
for a sufficiently smooth function f(x), where C is of order unity. 2.2. Wavelet compression and wavelet de-noising The major strength of wavelet filtering decomposition (2), is the ability to compress signals. For functions that contain isolated small scales on a large-scale background, most wavelet coefficients are small, thus, we can retain good approximation even after discarding a large number of wavelets with small coefficients. Intuitively, μ,j the coefficient dk will be small unless u(x) has variation on the scale of j in the immediate vicinity of wavelet μ,j ψk (x). Another important property of wavelet analysis used in this work is the ability of wavelets to de-noise signals. The wavelet de-noising procedure, also called waveletshrinkage, was introduced by Donoho (1993, 1994) based on orthogonal wavelet decompositions. It can be described as follows: given a function that consists of a smooth function with superimposed noise, one performs a forward wavelet transform and sets to zero “noisy” wavelet coefficients (i.e. those wavelet coefficients whose modulus squared is less than the noise variance σ 2 ), otherwise the wavelet coefficient is kept. This procedure is known as hard thresholding. Donoho (1993) demonstrated that hard thresholding is optimal for de-noising signals in the presence of Gaussian white noise. In the case of the plume data discussed in this work the “noise” is actually the areas of low grade, less coherent, temperature between the coherent plume regions. 3. Need for localized multi-scale visualization The multitude of data being generated from highresolution simulations of multi-scale phenomena, such
B.J. Kadlec et al. / Physics of the Earth and Planetary Interiors 163 (2007) 292–298
as convection in the Earth’s mantle, requires novel techniques for efficiently analyzing numerical results. Interactive visualization is the most feasible and natural tool for allowing human observation to extract scientific insight from such complex data. Advances in high-performance graphics processing and large highresolution LCD displays have allowed interactive 3D visualization to be a realizable dream that allows gigasized datasets to be investigated on commodity graphics workstations. As the range of scales and depth of resolution used within simulations become more vast and precise, new innovations must be developed to maintain and increase the interaction available with many of today’s simulations. Harnessing the multi-scale capabilities of wavelet analysis and coupling it with direct visualization is a novel solution to this problem. Wavelet analysis allows for the compression of datasets while retaining the multi-scale physical features of the dataset. In order for this to be visually realized, proper choice of visualization techniques and parameters must be made. In the case of our mantle convection data, we employ volume visualization to capture the structure and behavior of plume upwellings and downwellings. Representing these plumes as they were intended to be modeled requires data that includes connectivity between points and appropriate opacity-values that enhance plume features. In addition, three-dimensional interaction is necessary in order to determine the best visual representation of a dataset by looking at different perspectives, while recognizing key physical features Fig. 2. The powerful result of our work is being able to visualize large three-dimensional datasets with only a compressed subset of the full volume, allowing for larger datasets to be interactively explored. As the size of datasets continue to increase faster than the rate of memory and storage technology, this research will be timely to everyone computing volumes of data on the order of
Fig. 2. Volume rendered mantle convection simulation with Rayleigh number 10ˆ 8 shows the localized multi-scale plume structures to be wavelet compressed.
295
terabytes and larger. The end goal is allowing researchers to conduct scientific visualization with a highly reduced number of bytes that are still able to represent the key physical information in a dataset. 4. Visualization technique In this work we use direct wavelet space visualization and unstructured volume visualization in real space to study mantle convection data. This unit problem represents the complex visualization and data storage and management issues seen with larger more complex geophysical datasets from experiment and simulation. Direct wavelet space visualization is a technique for low cost qualitative data analysis. The unstructured volume visualization technique used falls between direct wavelet space visualization and structured volume visualization techniques in computation and memory cost, yet the quality of this technique is on par with full structured volume visualization techniques. The direct visualization of wavelet coefficient magnitudes allows for a quick low cost view of localized features in a dataset. By interactively thresholding the range of wavelet coefficient magnitudes from below, a small subset of the data can be used to highlight the location and general physical character of data features. This can be done easily by representing point centered data as spherical glyphs scaled and colored by wavelet coefficient magnitude. The cost of rendering glyphs is significantly lower than volume visualization allowing for exploration of even the largest datasets at interactive speeds for physically significant features. This technique can be used to determine areas of interest that can then be explored further with more costly visualization techniques. Volume visualization is an ideal technique for capturing structure and behavior of physical features in a dataset. Direct volume visualization of wavelet compressed datasets is a way to significantly reduce the cost of volume visualization in comparison with volume visualization on a full regular mesh. We note that, although state of the art regular mesh volume visualization algorithms are faster on a per cell basis, significant computational saving, in terms of the full data domain, can be realized with the unstructured volume visualization applied to the wavelet compressed data due to the significant levels of compression realized. In this work we use an octal tree-based ray-tracing volume rendering technique developed at the MultiScale Modeling and Simulation Laboratory at the University of Colorado-Boulder. First an octal tree-based cell structure must be created from wavelet compressed
296
B.J. Kadlec et al. / Physics of the Earth and Planetary Interiors 163 (2007) 292–298
points that lie on a subset of the original data’s regular mesh. This cell creation is done recursively based on wavelet level. The multi-level nature of the wavelet coefficients allows for a multi-resolution algorithm that is initiated on the scaling function points that are all retained in the wavelet thresholding algorithm. The points associated with increasing wavelet level are then iterated on at increasing level to create an efficient mesh representation. The technique is implemented in a cache-friendly fashion such that requests are minimized for information stored in main memory outside of the cache. 5. Results Mantle convection simulation data at 2563 with constant viscosity and Rayleigh number of 108 has been wavelet filtered at increasing threshold (ε) values resulting in dataset wavelet compression rates from 0% to 99.27%. Table 1 displays statistics on the range of epsilon values in terms of significant wavelet coefficients, number of points used in rendering, and fraction of data used. We consider the significant wavelet coefficients those that were above the threshold value thus were retained. Fig. 3 shows the wavelet space coefficients of the wavelet transformed mantle convection simulation data. The spherical glyphs are scaled and colored by coef-
Table 1 2563 Mantle convection wavelet compression Compression (%)
Points
Epsilon
0.00 16.90 66.20 94.96 98.30 99.27
16,777,216 13,939,031 5,667,187 845,846 284,642 121,955
0.00000 0.00050 0.00500 0.05000 0.10000 0.20000
ficient magnitude. In this figure the scaling function coefficients, that represent the mean of the dataset, are not shown. Only wavelet coefficients with magnitude greater than 0.05 are shown to highlight the significant plume structures. This is also a clear visual representation of the collocation of the retained data points after wavelet threshold filtering and areas of localized structures, in this case plumes, in the dataset. Fig. 4 shows three-dimensional volume visualization directly on unstructured meshes generated from the wavelet compressed datasets. Three-dimensional visualization was completed using ray-tracing with 5002 rays and an average of 170 points per ray. We show rendering results at six compressions (0%, 16.9%, 66.2%, 94.96%, 98.3%, and 99.27%). Threshold values of approximately 0.05 (94.96% compression) and lower generate render-
Fig. 3. Wavelet coefficients from the mantle convection simulation plotted as spherical glyphs that are scaled and colored according to wavelet coefficient magnitude.
B.J. Kadlec et al. / Physics of the Earth and Planetary Interiors 163 (2007) 292–298
297
Fig. 4. Mantle convection simulation at 0%, 16.9%, 66.2%, 94.96%, 98.3%, and 99.27% wavelet compression (left to right, top to bottom).
ing results that appear unchanged from the raw data (top four images). The lower two images show that at higher wavelet compression rates critical feature information is lost. The middle right plot corresponds to the
wavelet coefficients shown in Fig. 3. We can see from this that visualizing the wavelet coefficients themselves highlights the regions of plume activity accurately. Clearly the direct volume visualization does a better job of cap-
298
B.J. Kadlec et al. / Physics of the Earth and Planetary Interiors 163 (2007) 292–298
turing the structure and behavior of the plume upwellings and downwellings. 6. Conclusion Mantle convection data has been visualized both in wavelet space and real space in a significantly wavelet compressed form and is compared to regular uncompressed visualization. Wavelet coefficient visualizations have been shown that highlight the joint space/wavenumber representation of the wavelet coefficients. The significant plume structures in the dataset are clearly seen in the low overhead wavelet coefficient plots (Fig. 3). Direct volume rendering on the wavelet compressed datasets has been shown to compare well with the full dataset volume visualization up to a compression of 94.96% (Fig. 4). The techniques developed in this and ongoing work along with dynamically adaptive wavelet simulation techniques being developed in related work, will allow for simulation on a supercomputer, transfer of results to a visualization workstation and subsequent visualization to be all done at a fraction of the cost of using un-adapted Cartesian or zonal mesh representations. Future work is to develop more specialized algorithms for cell generation from adaptive wavelet collocation datasets and volume rendering techniques based on higher order interpolations using pre-calculated derivatives. These will allow for increased fidelity volume rendering on the wavelet collocation mesh. This is of highest importance in regions where the wavelet collocation points are extremely sparse causing a ghosting effect when cell-based linear interpolation methods are used. Acknowledgements We would like to thank Fabien W. Dubuffet for his mantle convection solution and Oleg Vasilyev and Gordon Erlebacher for stimulating discussions on wavelets and visualization.
References Daubechies, L., 1992. Ten Lectures on Wavelets. No. 61 in CBMS-NSF Series in Applied Mathematics. SIAM, Philadelphia. Donoho, D., 1993. Unconditional bases are optimal bases for data compression and for statistical estimation. Appl. Comput. Harmon. Anal. 1, 100–115. Donoho, D.L., 1992. Interpolating wavelet transforms. Tech. Rep. 408. Department of Statistics, Stanford University. http://playfair. Stanford.EDU/reports/donoho/interpol.ps.Z. Donoho, D.L., 1994. De-noising by soft-thresholding. IEEE Trans. Inf. Theory 41 (3), 613–627. Foufoula-Georgiou, E., Kumar, P. (Eds.), 1994. Wavelets in Geophysics. Academic Press. Goldstein, D., Vasilyev, O., Kevlahan, N.-R., 2005. CVS and SCALES simulations of 3D isotropic turbulence. J. Turbul. 06. (37). Goldstein, D.A., Vasilyev, O.V., 2004. Stochastic coherent adaptive large eddy simulation method. Phys. Fluids 16 (7), 2497– 2513. Mallat, S.G., 1999. A Wavelet Tour of Signal Processing. Academic Press, Paris. Schneider, K., Farge, M., Pellegrino, G., Rogers, M., 2005. Coherent vortex simulation of three-dimensional turbulent mixing layers using orthogonal wavelets. J. Fluid Mech. 534, 39– 66. Sweldens, W., 1996. The lifting scheme: a custom-design construction of biorthogonal wavelets. Appl. Comput. Harmon. Anal. 3 (2), 186–200. Sweldens, W., 1998. The lifting scheme: a construction of second generation wavelets. SIAM J. Math. Anal. 29 (2), 511– 546. Vasilyev, O.V., 2003. Solving multi-dimensional evolution problems with localized structures using second generation wavelets. Int. J. Comp. Fluid Dyn. 17 (2), 151–168 (Special issue on Highresolution Methods in Computational Fluid Dynamics). Vasilyev, O.V., Bowman, C., 2000. Second generation wavelet collocation method for the solution of partial differential equations. J. Comp. Phys. 165, 660–693. Vecsey, L., Hier-Majumder, C.A., Yuen, D.A., 2003. Multiresolution tectonic features over the earth inferred from a wavelet transformed geoid. Visual Geosci. 8, 26–44. Vezolainen, A., Erlebacher, G., Vasilyev, O., Yuen, D., 2005. Volumetric rendering of geophysical data on adaptive wavelet grid. Eos Trans. AGU 86, IN43B-0335. Yuen, D.A., Erlebacher, G., Vasilyev, O.V., Goldstein, D.E., Fuentes, M., 2004. Role of wavelets in the physical and statistical modelling of complex geological processes. Pure Appl. Geophys. 161, 2231–2244.
Visualization of wavelet compressed mantle convection data Benjamin J. Kadlec a,∗ , Daniel E. Goldstein b , David A. Yuen c , Alexei Vezolainen d a
b
Department of Computer Science, University of Colorado at Boulder, Boulder, CO 80309, United States NorthWest Research Associates, Inc., Colorado Research Associates Division, Boulder, CO 80301, United States c Department of Geology and Geophysics and Minnesota Supercomputing Institute, University of Minnesota, Minneapolis, MN 55455, United States d Department of Mechanical Engineering, University of Colorado at Boulder, Boulder, CO 80309, United States Received 5 February 2007; accepted 20 March 2007
Abstract Wavelet compression methodologies promise significant computational and storage cost savings through the benefit of producing results on adapted grids that significantly reduce storage and data manipulation costs. Yet, with these storage methodologies come new challenges in the visualization of temporally adaptive datasets. Wavelet compressed three-dimensional mantle convection simulations are visualized as a challenging case study. We employ adaptive 3D visualization using cache friendly octal tree volume visualization as an effective visualization approach. To explore the localized multi-scale structures in the datasets, the wavelet coefficients are also visualized allowing visualization of energy contained in local physical regions as well as in local wavenumber space. Our technique allows for exploring localized multi-scale structures using multi-resolution visualization from the compressed state of large geophysical simulations. © 2007 Elsevier B.V. All rights reserved. Keywords: 3D visualization; Wavelets; Mantle convection; Volume visualization
1. Introduction The race towards petascale computing has created an explosion in data from large-scale computations of nonlinear phenomena in many fields, including the geo- and environmental sciences. Today, the geosciences are facing an onslaught of information coming from terabyte-scale numerical simulations of nonlinear phenomena and higher resolution data from satellite missions and precise laboratory experiments. In addi∗
Corresponding author. E-mail addresses: [email protected] (B.J. Kadlec), [email protected] (D.E. Goldstein), [email protected] (D.A. Yuen), [email protected] (A. Vezolainen). 0031-9201/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.pepi.2007.03.009
tion, great torrents of data are being produced by new earthquake models coupled with sensing information, which can easily overwhelm researchers due to their high dynamic range and multi-scale character in space and time. As scientists in the geo and environmental disciplines continue their march towards petascale computing, the cost of large-scale data analysis will be a dire problem requiring novel tools to process and represent data more efficiently. The efficient storage and subsequent visualization of these large datasets is a trade off in storage costs versus data quality. The last decade and a half has witnessed the development of wavelets or wavelet analysis. Wavelets have a wide range of application in areas such as feature extraction, data compression, numerical simulation (Goldstein
B.J. Kadlec et al. / Physics of the Earth and Planetary Interiors 163 (2007) 292–298
et al., 2005; Goldstein and Vasilyev, 2004; Schneider et al., 2005), and visualization (Mallat, 1999). Wavelets are an effective tool for efficiently representing data with localized multi-scale structures due to the inherent localized support of the wavelet basis functions in both physical and wavenumber space. This dual-space localization property of wavelets allows for wavelet compressed datasets to represent localized multi-scale structures with a high level of compression and minimal loss of important features. The ability of wavelets to efficiently represent data containing localized multi-scale structure makes them an ideal tool for the simulation of geophysical systems and the subsequent storage and visualization of simulation and experimental results (Yuen et al., 2004; Vecsey et al., 2003; Vezolainen et al., 2005) In this work our emphasis will be on the issues and possible dramatic advantages to be gained by representing existing geo-physical datasets in a wavelet compressed form for storage and visualization. This paper is organized as follows: we will first describe the basic properties of second generation wavelets. Then we will discuss the de-noising properties through wavelet thresholding, which allow decomposition of localized multi-scale coherent structures from an incoherent background. Next we will motivate the need for visualization in our work and describe the visualization techniques we derive. Following this we will present, as an example case, the visualization of a mantle convection simulation dataset at different levels of wavelet compression. We will conclude with some discussion of future improvements that are being studied for the visualization of wavelet compressed datasets.
293
2. Wavelets Wavelets are basis functions which are localized in both physical space (due to their finite support) and wavenumber space (due to their vanishing moments), e.g. Fig. 1. For comparison, the classical Fourier transform is based on functions (sines and cosines) that are well localized in wavenumber space, but do not provide localization in physical space due to their global support. Because of this space/scale localization, the wavelet transform provides both spatial and scale (frequency) information while the Fourier transform on the other hand only provides frequency information. A scalar field f(x) can be represented in terms of wavelet basis functions as f (x) =
+∞ 2 −1 n
cl0 φl0 (x) +
μ,j
μ,j
dk ψk (x), (1)
j=0 μ=1 k ∈ κμ,j
l ∈ L0 μ,j
where φl0 (x) and ψk are respectively n-dimensional scaling functions and wavelets of different families (μ) and levels of resolution (j). One may think of a wavelet decomposition as a multi-level or multi-resolution representation of a function, where each level of resolution j j (except the coarsest one) consists of wavelets ψk or μ,j family of wavelets ψk having the same scale but located at different positions. Scaling function coefficients represent the averaged values of the field, while the wavelet coefficients represent the details of the field at different scales. The wavelet functions have a zero mean, while the scaling functions do not. Note that in n-dimensions there are 2n − 1 distinctive n-dimensional wavelets (Daubechies, 1992). Also note that due to the
Fig. 1. Lifted interpolating wavelet ψ, of order 6 (a) and its Fourier transform (b).
294
B.J. Kadlec et al. / Physics of the Earth and Planetary Interiors 163 (2007) 292–298
local support of both scaling functions and wavelets, there is a one-to-one correspondence between the location of each scaling function or wavelet with a grid point. As a result each scaling function coefficient cl0 and each μ,j wavelet coefficient dk is uniquely associated with a single grid point with the indices l and k, respectively. Traditionally, one-dimensional first generation j wavelets ψk are defined as translates and dilates of j one basic wavelet ψ, i.e. ψk (x) = ψ(2j x − k). Second generation wavelets (Sweldens, 1996, 1998) are a generalization of first generation wavelets that supplies additional freedom to deal with arbitrary boundary conditions, and irregular sampling intervals. Second generation wavelets form a Riesz basis for L2 space, with the wavelets being local in both space and frequency and often having many vanishing polynomial moments, but without the translation and dilation invariance of their first generation cousins. Despite the loss of these two fundamental properties of wavelet bases, second generation wavelets retain many of the useful features of first generation wavelets, including a fast O(N) transform. The construction of second generation wavelets is based on the lifting scheme that is discussed in detail by Sweldens (1996, 1998). For this study we use a set of second generation wavelets known in the literature as lifted interpolating wavelets (Vasilyev and Bowman, 2000; Sweldens, 1996). In particular a lifted interpolating wavelet of order 6, which is shown in Fig. 1 along with its Fourier transform is used in this work. For a more in-depth discussion on the construction of these wavelets the reader is referred to the papers by Sweldens (1996, 1998), and Vasilyev and Bowman (2000). For a more general discussion on wavelets we refer the reader to the books of Daubechies (1992), Mallat (1999), and FoufoulaGeorgiou and Kumar (1994). 2.1. Wavelet filters Wavelet filtering is performed in wavelet space using wavelet coefficient thresholding, which can be considered as a nonlinear filter that depends on each flow realization. The wavelet thresholding filter is defined by f¯ >ε (x) = cl0 φl0 (x) l ∈ L0 +∞ 2 −1 n
+
j=0 μ=1
k ∈ κμ,j
μ,j
|dk |> ∈ ||f ||WTF
μ,j
μ,j
dk ψk (x),
(2)
where f(x) is a scalar field, ε > 0 stands for the nondimensional (relative) threshold parameter, and ||·||WTF being the Wavelet Threshold Filtering (WTF) norm that provides the (absolute) dimensional scaling for filtered variable f. For instance, in the case of velocity, the (absolute) dimensional scaling can be specified as the L2 norm (||ui ||WTF = ||ui ||2 ) or the L∞ norm (||ui ||WTF = ||ui ||∞ ). Note that once the WTF-norm ||·||WTF is specified, the wavelet thresholding filter (2) is uniquely defined by the non-dimensional threshold parameter, ε. The reconstruction error due to wavelet filtering with non-dimensional threshold parameter ε can be shown to be (Vasilyev, 2003; Donoho, 1992): ||f (x) − f¯ >ε (x)||2 ≤ Cε||f ||WTF ,
(3)
for a sufficiently smooth function f(x), where C is of order unity. 2.2. Wavelet compression and wavelet de-noising The major strength of wavelet filtering decomposition (2), is the ability to compress signals. For functions that contain isolated small scales on a large-scale background, most wavelet coefficients are small, thus, we can retain good approximation even after discarding a large number of wavelets with small coefficients. Intuitively, μ,j the coefficient dk will be small unless u(x) has variation on the scale of j in the immediate vicinity of wavelet μ,j ψk (x). Another important property of wavelet analysis used in this work is the ability of wavelets to de-noise signals. The wavelet de-noising procedure, also called waveletshrinkage, was introduced by Donoho (1993, 1994) based on orthogonal wavelet decompositions. It can be described as follows: given a function that consists of a smooth function with superimposed noise, one performs a forward wavelet transform and sets to zero “noisy” wavelet coefficients (i.e. those wavelet coefficients whose modulus squared is less than the noise variance σ 2 ), otherwise the wavelet coefficient is kept. This procedure is known as hard thresholding. Donoho (1993) demonstrated that hard thresholding is optimal for de-noising signals in the presence of Gaussian white noise. In the case of the plume data discussed in this work the “noise” is actually the areas of low grade, less coherent, temperature between the coherent plume regions. 3. Need for localized multi-scale visualization The multitude of data being generated from highresolution simulations of multi-scale phenomena, such
B.J. Kadlec et al. / Physics of the Earth and Planetary Interiors 163 (2007) 292–298
as convection in the Earth’s mantle, requires novel techniques for efficiently analyzing numerical results. Interactive visualization is the most feasible and natural tool for allowing human observation to extract scientific insight from such complex data. Advances in high-performance graphics processing and large highresolution LCD displays have allowed interactive 3D visualization to be a realizable dream that allows gigasized datasets to be investigated on commodity graphics workstations. As the range of scales and depth of resolution used within simulations become more vast and precise, new innovations must be developed to maintain and increase the interaction available with many of today’s simulations. Harnessing the multi-scale capabilities of wavelet analysis and coupling it with direct visualization is a novel solution to this problem. Wavelet analysis allows for the compression of datasets while retaining the multi-scale physical features of the dataset. In order for this to be visually realized, proper choice of visualization techniques and parameters must be made. In the case of our mantle convection data, we employ volume visualization to capture the structure and behavior of plume upwellings and downwellings. Representing these plumes as they were intended to be modeled requires data that includes connectivity between points and appropriate opacity-values that enhance plume features. In addition, three-dimensional interaction is necessary in order to determine the best visual representation of a dataset by looking at different perspectives, while recognizing key physical features Fig. 2. The powerful result of our work is being able to visualize large three-dimensional datasets with only a compressed subset of the full volume, allowing for larger datasets to be interactively explored. As the size of datasets continue to increase faster than the rate of memory and storage technology, this research will be timely to everyone computing volumes of data on the order of
Fig. 2. Volume rendered mantle convection simulation with Rayleigh number 10ˆ 8 shows the localized multi-scale plume structures to be wavelet compressed.
295
terabytes and larger. The end goal is allowing researchers to conduct scientific visualization with a highly reduced number of bytes that are still able to represent the key physical information in a dataset. 4. Visualization technique In this work we use direct wavelet space visualization and unstructured volume visualization in real space to study mantle convection data. This unit problem represents the complex visualization and data storage and management issues seen with larger more complex geophysical datasets from experiment and simulation. Direct wavelet space visualization is a technique for low cost qualitative data analysis. The unstructured volume visualization technique used falls between direct wavelet space visualization and structured volume visualization techniques in computation and memory cost, yet the quality of this technique is on par with full structured volume visualization techniques. The direct visualization of wavelet coefficient magnitudes allows for a quick low cost view of localized features in a dataset. By interactively thresholding the range of wavelet coefficient magnitudes from below, a small subset of the data can be used to highlight the location and general physical character of data features. This can be done easily by representing point centered data as spherical glyphs scaled and colored by wavelet coefficient magnitude. The cost of rendering glyphs is significantly lower than volume visualization allowing for exploration of even the largest datasets at interactive speeds for physically significant features. This technique can be used to determine areas of interest that can then be explored further with more costly visualization techniques. Volume visualization is an ideal technique for capturing structure and behavior of physical features in a dataset. Direct volume visualization of wavelet compressed datasets is a way to significantly reduce the cost of volume visualization in comparison with volume visualization on a full regular mesh. We note that, although state of the art regular mesh volume visualization algorithms are faster on a per cell basis, significant computational saving, in terms of the full data domain, can be realized with the unstructured volume visualization applied to the wavelet compressed data due to the significant levels of compression realized. In this work we use an octal tree-based ray-tracing volume rendering technique developed at the MultiScale Modeling and Simulation Laboratory at the University of Colorado-Boulder. First an octal tree-based cell structure must be created from wavelet compressed
296
B.J. Kadlec et al. / Physics of the Earth and Planetary Interiors 163 (2007) 292–298
points that lie on a subset of the original data’s regular mesh. This cell creation is done recursively based on wavelet level. The multi-level nature of the wavelet coefficients allows for a multi-resolution algorithm that is initiated on the scaling function points that are all retained in the wavelet thresholding algorithm. The points associated with increasing wavelet level are then iterated on at increasing level to create an efficient mesh representation. The technique is implemented in a cache-friendly fashion such that requests are minimized for information stored in main memory outside of the cache. 5. Results Mantle convection simulation data at 2563 with constant viscosity and Rayleigh number of 108 has been wavelet filtered at increasing threshold (ε) values resulting in dataset wavelet compression rates from 0% to 99.27%. Table 1 displays statistics on the range of epsilon values in terms of significant wavelet coefficients, number of points used in rendering, and fraction of data used. We consider the significant wavelet coefficients those that were above the threshold value thus were retained. Fig. 3 shows the wavelet space coefficients of the wavelet transformed mantle convection simulation data. The spherical glyphs are scaled and colored by coef-
Table 1 2563 Mantle convection wavelet compression Compression (%)
Points
Epsilon
0.00 16.90 66.20 94.96 98.30 99.27
16,777,216 13,939,031 5,667,187 845,846 284,642 121,955
0.00000 0.00050 0.00500 0.05000 0.10000 0.20000
ficient magnitude. In this figure the scaling function coefficients, that represent the mean of the dataset, are not shown. Only wavelet coefficients with magnitude greater than 0.05 are shown to highlight the significant plume structures. This is also a clear visual representation of the collocation of the retained data points after wavelet threshold filtering and areas of localized structures, in this case plumes, in the dataset. Fig. 4 shows three-dimensional volume visualization directly on unstructured meshes generated from the wavelet compressed datasets. Three-dimensional visualization was completed using ray-tracing with 5002 rays and an average of 170 points per ray. We show rendering results at six compressions (0%, 16.9%, 66.2%, 94.96%, 98.3%, and 99.27%). Threshold values of approximately 0.05 (94.96% compression) and lower generate render-
Fig. 3. Wavelet coefficients from the mantle convection simulation plotted as spherical glyphs that are scaled and colored according to wavelet coefficient magnitude.
B.J. Kadlec et al. / Physics of the Earth and Planetary Interiors 163 (2007) 292–298
297
Fig. 4. Mantle convection simulation at 0%, 16.9%, 66.2%, 94.96%, 98.3%, and 99.27% wavelet compression (left to right, top to bottom).
ing results that appear unchanged from the raw data (top four images). The lower two images show that at higher wavelet compression rates critical feature information is lost. The middle right plot corresponds to the
wavelet coefficients shown in Fig. 3. We can see from this that visualizing the wavelet coefficients themselves highlights the regions of plume activity accurately. Clearly the direct volume visualization does a better job of cap-
298
B.J. Kadlec et al. / Physics of the Earth and Planetary Interiors 163 (2007) 292–298
turing the structure and behavior of the plume upwellings and downwellings. 6. Conclusion Mantle convection data has been visualized both in wavelet space and real space in a significantly wavelet compressed form and is compared to regular uncompressed visualization. Wavelet coefficient visualizations have been shown that highlight the joint space/wavenumber representation of the wavelet coefficients. The significant plume structures in the dataset are clearly seen in the low overhead wavelet coefficient plots (Fig. 3). Direct volume rendering on the wavelet compressed datasets has been shown to compare well with the full dataset volume visualization up to a compression of 94.96% (Fig. 4). The techniques developed in this and ongoing work along with dynamically adaptive wavelet simulation techniques being developed in related work, will allow for simulation on a supercomputer, transfer of results to a visualization workstation and subsequent visualization to be all done at a fraction of the cost of using un-adapted Cartesian or zonal mesh representations. Future work is to develop more specialized algorithms for cell generation from adaptive wavelet collocation datasets and volume rendering techniques based on higher order interpolations using pre-calculated derivatives. These will allow for increased fidelity volume rendering on the wavelet collocation mesh. This is of highest importance in regions where the wavelet collocation points are extremely sparse causing a ghosting effect when cell-based linear interpolation methods are used. Acknowledgements We would like to thank Fabien W. Dubuffet for his mantle convection solution and Oleg Vasilyev and Gordon Erlebacher for stimulating discussions on wavelets and visualization.
References Daubechies, L., 1992. Ten Lectures on Wavelets. No. 61 in CBMS-NSF Series in Applied Mathematics. SIAM, Philadelphia. Donoho, D., 1993. Unconditional bases are optimal bases for data compression and for statistical estimation. Appl. Comput. Harmon. Anal. 1, 100–115. Donoho, D.L., 1992. Interpolating wavelet transforms. Tech. Rep. 408. Department of Statistics, Stanford University. http://playfair. Stanford.EDU/reports/donoho/interpol.ps.Z. Donoho, D.L., 1994. De-noising by soft-thresholding. IEEE Trans. Inf. Theory 41 (3), 613–627. Foufoula-Georgiou, E., Kumar, P. (Eds.), 1994. Wavelets in Geophysics. Academic Press. Goldstein, D., Vasilyev, O., Kevlahan, N.-R., 2005. CVS and SCALES simulations of 3D isotropic turbulence. J. Turbul. 06. (37). Goldstein, D.A., Vasilyev, O.V., 2004. Stochastic coherent adaptive large eddy simulation method. Phys. Fluids 16 (7), 2497– 2513. Mallat, S.G., 1999. A Wavelet Tour of Signal Processing. Academic Press, Paris. Schneider, K., Farge, M., Pellegrino, G., Rogers, M., 2005. Coherent vortex simulation of three-dimensional turbulent mixing layers using orthogonal wavelets. J. Fluid Mech. 534, 39– 66. Sweldens, W., 1996. The lifting scheme: a custom-design construction of biorthogonal wavelets. Appl. Comput. Harmon. Anal. 3 (2), 186–200. Sweldens, W., 1998. The lifting scheme: a construction of second generation wavelets. SIAM J. Math. Anal. 29 (2), 511– 546. Vasilyev, O.V., 2003. Solving multi-dimensional evolution problems with localized structures using second generation wavelets. Int. J. Comp. Fluid Dyn. 17 (2), 151–168 (Special issue on Highresolution Methods in Computational Fluid Dynamics). Vasilyev, O.V., Bowman, C., 2000. Second generation wavelet collocation method for the solution of partial differential equations. J. Comp. Phys. 165, 660–693. Vecsey, L., Hier-Majumder, C.A., Yuen, D.A., 2003. Multiresolution tectonic features over the earth inferred from a wavelet transformed geoid. Visual Geosci. 8, 26–44. Vezolainen, A., Erlebacher, G., Vasilyev, O., Yuen, D., 2005. Volumetric rendering of geophysical data on adaptive wavelet grid. Eos Trans. AGU 86, IN43B-0335. Yuen, D.A., Erlebacher, G., Vasilyev, O.V., Goldstein, D.E., Fuentes, M., 2004. Role of wavelets in the physical and statistical modelling of complex geological processes. Pure Appl. Geophys. 161, 2231–2244.