Volumetric Rendering on Wavelet-Based Adaptive Grid

: Numerical modeling of physical phenomena frequently involves processes across a wide range of spatial and temporal scales. In the last two decades, the advancements in wavelet-based numerical methodologies to solve partial differential equations, combined with the unique properties of wavelet analysis to resolve localized structures of the solution on dynamically adaptive computational meshes, make it feasible to perform large-scale numerical simulations of a variety of physical systems on a dynamically adaptive computational mesh that changes both in space and time. Volumetric visualization of the solution is an essential part of scientiﬁc computing, yet the existing volumetric visualization techniques do not take full advantage of multi-resolution wavelet analysis and are not fully tailored for visualization of a compressed solution on the wavelet-based adaptive computational mesh. Our objective is to explore the alternatives for the visualization of time-dependent data on space-time varying adaptive mesh using volume rendering while capitalizing on the available sparse data representation. Two alternative formulations are explored. The ﬁrst one is based on volumetric ray casting of multi-scale datasets in wavelet space. Rather than working with the wavelets at the ﬁnest possible resolution, a partial inverse wavelet transform is performed as a preprocessing step to obtain scaling functions on a uniform grid at a user-prescribed resolution. As a result, a solution in physical space is represented by a superposition of scaling functions on a coarse regular grid and wavelets on an adaptive mesh. An efﬁcient and accurate ray casting algorithm is based just on these coarse scaling functions. Additional details are added during the ray tracing by taking an appropriate number of wavelets into account based on support overlap with the interpolation point, wavelet coefﬁcient magnitude, and other characteristics, such as opacity accumulation (front to back ordering) and deviation from frontal viewing direction. The second approach is based on complementing of wavelet-based adaptive mesh to the traditional Adaptive Mesh Reﬁnement (AMR) mesh. Both algorithms are illustrated and compared to the existing volume visualization software for Rayleigh-Benard thermal convection and electron density data sets in terms of rendering time and visual quality for different data compression of both wavelet-based and AMR adaptive meshes.


Introduction
Numerical modeling of physical phenomena frequently involves processes across a wide range of spatial and temporal scales.To ensure better capturing of the flow physics on a near optimal adaptive computational mesh with substantially smaller computational cost, while resolving the localized structures of the solution, a number of wavelet-based adaptive computational methods have been developed recently (e.g., [1][2][3][4]).
Similarly to adaptive mesh refinement (AMR) methods (e.g., [5][6][7]), adaptive waveletbased methods [3,8] provide a convenient and efficient approach to resolving multi-scale processes on dynamically adaptive computational mesh, which would be hard-pressed to be solved by conventional finite-difference techniques on a uniform non-adaptive grid.Wavelet-based methods take advantage of the wavelet compression properties, as a result, functions with localized regions of sharp transition are well compressed using wavelet decomposition briefly discussed in Section 2. The adaptation is achieved by retaining only those wavelets, whose coefficients are greater than a given wavelet threshold (see Equation (2)).Thus, high-resolution computations are carried out only in those regions, where sharp transitions occur.In addition to controlling global L 2 approximation error, wavelet based adaptive grid provides a compression factor of 10-10 2 relative to the uniform grid of the same effective resolution.
Since all the computations are performed in wavelet space on the adaptive grid, it makes sense to utilize the available data compression and perform data visualization in wavelet space on the same adaptive grid.Despite the wide use of AMR methods in high performance computing (e.g., [9,10]) and fast growing development of adaptive waveletbased methods (e.g., [3,4,[11][12][13][14]) visualizing data on adaptive grids has only limited level of support in widely used visualization packages such as VTK [15], ParaView [16], VisIt [17], and ChomboVis [18].Rendering such data on adaptive meshes interactively with high visual quality, i.e., without holes and artifacts, still constitutes a challenge [19].
One of the first adaptive mesh volume rendering systems based on cell projection volume rendering was introduced by Ma and Crockett [20].Over the years, a number of approaches for volume rendering of data on AMR meshes have been developed.For example, a slice-based volume rendering of 3D textures [21] and its multi-pass [22] and single-pass [23] extensions were developed for adaptive mesh overlapping coarse and fine-resolution AMR levels.An efficient method for rendering AMR data-based interblock interpolation was developed in [24].Examples of multi-resolution volume rendering methods for adaptive meshes based on the submission of lower-resolution hexahedral cells and local interpolants can be found in [25].
Despite the development of many successful volume-rendering methods for AMR data and the constantly growing need for widely available AMR data visualization tools, these methods have not been implemented in existing visualization packages.We present here several techniques for volume rendering of multi-scale data sets while capitalizing on the available sparse data representation.The first visualization technique is based on the direct summation of wavelet coefficients during volumetric ray casting.The second approach is to complement our adaptive wavelet grid with the appropriate number of nodes in order to produce a traditional AMR grid.Volume rendering on an AMR grid is performed either with the existing or with our own software.Our work has some parallels with a recently developed non-wavelet-based volume rendering method for representing and traversing adaptive octree data [19], where, similarly to the present paper, it was also demonstrated that the developed approach could be adopted for production use in ParaView [16].

Wavelet Based Solution of PDEs
The use of second generation wavelets [26] for the solution of partial differential equations has been extensively discussed in [4, 27,28].Partial differential equations are discretized at collocation points thus leading to a system of algebraic equations for the unknown wavelet coefficients at these collocation points.Grid adaptation is obtained by keeping the points with wavelet coefficients (from a previous time step) larger than a predefined threshold limit.
As a result of the wavelet transform, a function in physical space f (x) is represented by a superposition of scaling functions φ 0 l (x) (l ∈ L 0 ), at the coarsest level of resolution, and wavelets ψ µ,j k (x) (k ∈ K µ,j ) of different families (µ = 1, . . ., 2 n − 1) and levels of resolution (j = 1, . . ., j max ) on an adaptive mesh: where n is the dimensionality of the space, the bold subscripts l and k denote n-dimensional indices, while L 0 and K µ,j stand for the associated index sets, and J is the maximum level of resolution that is present or allowed in the wavelet approximation.The scaling functions φ 0 l and the wavelets ψ µ,j k are constructed on a set of nested tensorial meshes with one-to-one correspondence between grid points and functions.The coefficients c 0 l and d µ,j k represent, respectively, the averaged values and the local variation details of the field f (x) at different scales.The value of f (x) in physical space can be obtained either through an inverse wavelet transform to a regular grid of required resolution or through a direct application of Equation ( 1) at the required point of physical space.Scaling and wavelet functions normally have local support (see one-dimensional scaling function φ(ξ) and wavelet ψ(ξ) Figure 1), which significantly limits the number of coefficients to sum in the Equation (1).
Wavelet-compressed solution f > (x) on an adaptive mesh is obtained in the wavelet space by using wavelet coefficient thresholding.Namely, the wavelet filtered function is defined by where > 0 stands for the non-dimensional (relative) threshold value, f being the (absolute) dimensional scale of f .For the details of the adaptive wavelet collocation method, the reader is referred to Refs.[4, 27,28].Note that a substantially larger threshold value can be used for visualization purposes compared to the one used in computations.
Thus, an additional level of compression is possible if the data are further compressed with larger , which would further improve the efficiency of the visualization.

Volume Rendering in Compression Domain
Compression properties of wavelet transform have been widely used for the visualization of large data sets [29][30][31][32][33]. Regularly sampled data have been projected into a wavelet basis.Depending on the required image quality, a large number of wavelet coefficients were neglected, thus leading to a significant reduction of memory usage during a local rendering (or during a network transmission for the rendering on a remote client).Interactive time rendering has been implemented for X-ray-like images via wavelet splatting and Fourier projection slice theorem [34,35].Another interactive rendering technique, a walk through the large data set by Guthe et al. [32] has been performed with the help of a block-wise wavelet transform [36,37].The relevant blocks of voxels have been decompressed on the fly and rendered using hardware texture mapping [32].
Ray tracing in wavelet domain has been performed in non-interactive time.A function value f (x) has been reconstructed from its wavelet decomposition and used for the computing of rendering integral [31].Alternatively, separate wavelet decompositions have been used for the reconstruction of opacity and color values used in the computing of rendering integral [38].
Interactive frame rates have been obtained in a parallel ray casting implementation by using block-wise wavelet transform on a multiprocessor system [36].Also, interactive rendering of large data sets (fastest along axes directions) has been obtained by using "thick rays" in a shear-warped "foveated volume" [39].This approach resembles volume rendering on adaptive mesh refinement (AMR) grids.
Adaptive mesh refinement has been extensively used for the simulations of multi-scale physical phenomena.A number of visualization techniques have been developed directly for the adaptive meshes, without resampling the data on a uniform grid (e.g., [40][41][42]).Extrapolation of wavelet-based grid data into some additional nodes will generate a correspondent AMR grid suitable for the application of the existing adaptive mesh visualization algorithms and software.

Numerical Results and Discussion
A number of open-source software for volume rendering is currently available.Mostly, the rendering is performed on a regular grid, the size of which is limited by the available machine memory, e.g., VolPack library [43,44].Unstructured grid software is more appropriate for the purposes of wavelet-based data visualization.VTK, the visualization toolkit, provides a number of libraries capable of dealing with unstructured grids [45].ParaView [16], a VTK-based parallel visualization application, can display data on an unstructured grid [16].Additionally, VisIt [17] and ChomboVis [18] are capable of displaying volume data on unstructured adaptive mesh.

Available Volume Rendering Software
A number of open-source software for volume rendering is currently available.Mostly, the rendering is performed on a regular grid, the size of which is limited by the available machine memory, e.g., VolPack library [43,44].Unstructured grid software is more appropriate for the purposes of wavelet-based data visualization.VTK, the visualization toolkit, provides a number of libraries capable of dealing with unstructured grids [45].ParaView, a VTK-based parallel visualization application, can display data on an unstructured grid [16].Additionally, VisIt [17] and ChomboVis [18] are capable of displaying volume data on unstructured adaptive mesh.

The Proposed Rendering Techniques
Several techniques for volume rendering of multi-scale data sets are presented.The first approach is based on the direct summation of wavelet coefficients during volumetric ray casting.

Direct Summation of Wavelets
First of all, rather than working with a wavelet scaling function of the largest possible support at the initially defined coarse grid, we perform a partial inverse wavelet transform as a preprocessing step to obtain scaling functions on a coarse grid at a user-prescribed resolution (64 3 seems to be an optimal one for the current implementation).This preprocessing step is performed in a fraction of a second, although it may require some additional memory to store new scaling coefficients.
Interactive time ray casting algorithm is based just on these scaling coefficients.Additional detail is added during the ray tracing by taking an appropriate number of wavelets into account based on support overlap with the interpolation point and wavelet amplitude.Other characteristics, such as opacity accumulation (front to back ordering) and deviation from frontal viewing direction could also be taken into account.

Ray Casting
The main goal of ray casting is to evaluate the solution of radiative transfer equation: where color and opacity are determined by the function f (x) in physical space, C = C( f (x)) and α = α( f (x)), respectively.In the current implementation we use the following approximation of the rendering integral: where i corresponds to the index of ray point, d i is the current step along the ray, and n(α) is the total number of points along the ray.
In our current implementation of the algorithm, the step along the ray d i is constant, depending on the number of fine levels on the adaptive mesh.A coarse image (of 64 3 scaling coefficients) is obtained in an interactive time by linear interpolations from the nearest eight coarse mesh scaling coefficients into the correspondent points on a ray.Computing of additional details requires the summation of fine level wavelet coefficients, which is currently non-interactive.It should also be noted that no hardware acceleration is used in the current implementation.

Data Structure
A data structure has been built for the rendering process (Figure 2).A linear array of coarse mesh nodes has been created for the storage of the scaling coefficients (an array of size of 64 3 in the current implementation) (Figure 2).Each coarse mesh node keeps the values of 8 nearest scaling coefficients.It should be noted that such redundancy provides 30% or so acceleration versus a non-redundant array of the scaling coefficients due to a better cache coherence of the redundant structure.Each coarse mesh node has a pointer to a dynamically allocated adaptive array of wavelet coefficients located inside the correspondent coarse mesh block ("Wavelet Nodes" at Figure 2).The length of the array of wavelet coefficients is determined during the preprocessing step and is a constant for the rendering of time-independent data.The length is also stored inside the coarse mesh node.While computing the function f (x) value during ray tracing by using the Equation (1), one cache miss happens during access to the scaling coefficients and at least one cache miss happens during the summation of the wavelet coefficients inside each coarse mesh block.To increase the cache coherence of the data structure, the wavelet coefficients in the adaptive array are sorted by the wavelet levels and types.Each wavelet coefficient is stored together with its coordinates and wavelet type.Though storing of wavelet type is redundant, it saves computation time during the expensive operation of summation of wavelet coefficients according to Equation (1) thus leading to an additional speedup of 10% or so.
In order to perform the summation of wavelet coefficients even faster, we compute the values of scaling and wavelet functions φ 0 i (x) and ψ j i (x) through a table lookup.Additionally, we artificially decrease the support of scaling and wavelet functions by introducing "cutoff" limits ξ 0 such that φ(ξ), ψ(ξ) ≡ 0, for ξ > ξ 0 (Figure 1).The current implementation uses "cutoff" limits 1 and 2 for scaling and wavelet functions, respectively.

Adaptive Wavelet Grid and AMR
The second approach to volume rendering of wavelet-based data is to complement the adaptive wavelet grid with the required number of nodes in order to produce a traditional AMR grid.Though some additional storage is required for the supplementary nodes, several software products have been already developed for volume rendering of AMR grid data.

Completing Wavelet Based Grid to AMR
As a preprocessing step, first, we create a node-centered AMR grid.Then an inverse wavelet transform is performed on all the nodes of the obtained AMR grid.It should be noted that contrary to the previously considered approach of direct summation in wavelet space, in this approach, we deal with the actual function values on a traditional nodecentered AMR grid.The preprocessing step requires some additional memory and may take a significant time, depending on the complexity of the wavelet grid.Although, for a time-dependent problem the timing of the preprocessing step is expected to be of the same order of magnitude as the computational time for a single time step.
The process of AMR creation is illustrated for a two-dimensional grid at Figure 3. Any wavelet node inside a cell will induce that cell "bisection" up to the level of that node.A wavelet node located at the cell's face or edge will induce "bisection" of all of the adjacent cells.

AMR Volume Rendering
To perform volume rendering, we write the obtained AMR grid data in the form of an unstructured (tetrahedral) mesh VTK file [VTK, 2000].This format is readable by Paraview, VisIt, and by VTK itself.A different input file format, HDF5, is required by ChomboVis.Although ChomboVis provides some volume rendering for cell-centered AMR data [Ligocki, 2000] and the opacity function can be tabulated, ChomboVis uses gray colormap for its volume rendering.Additionally, the quality of volume rendering is low and seems to be inappropriate for our purposes of visualization.Thus, AMR grid data are to be written in an unstructured mesh VTK format.
Each cubic cell of the AMR grid is considered to be a combination of six tetrahedra (Figure 4).If we number all the nodes in the cell according to the numbering of VTK_VOXEL cell type, six correspondent tetrahedra, or VTK_TETRA cells, will be: 0125, 1325, 0245, 2645, 3725, 7625 [VTK, 2000].It should be noted that such a triangulation is not compatible across the neighboring cells (i.e., T-junctions may appear) and therefore is potentially subject to various rendering artifacts, including gap (empty space) appearance during isosurface generation [46].Note when a size of the input file is a concern one may save AMR data in a VTK_VOXEL cell format instead of VTK_TETRA cell format [VTK, 2000].It might be convenient for slice or isosurface generation.As for the volume rendering, VTK seems to perform it on an unstructured tetrahedral mesh only; hence the volume is still to be tetrahedralized (and the additional memory for the tetrahedra to be requested).

Results
We have timed our volume rendering implementations for several data sets on a single processor system: Pentium 4, 3.2 GHz, 2G RAM, with Nvidia GeForce Fx5200-TD128 graphics card.Preprocessing time was not taken into account.The direct summation algorithm has been implemented through a C++ code and the rendering time was the time to compute the values of all the pixels of 500 × 500 image matrix.VTK-based rendering has been performed through a Tcl/Tk script and the rendering time has been computed as the time of execution of vtkRenderWindow class command Render.We note that the rendering studies were performed on old generation system and the reported timings would be much shorter for systems with more modern processors.

Convection Data Sets
Constant viscosity 256 3 thermal convection data with Rayleigh number of 10 10 has been processed to simulate the output of a wavelet-based PDE solver.A 4-th order second generation wavelet transform [26] has been used for the data processing.The data has been assumed to be periodic in order not to deal with the distortions of wavelet functions near the boundaries.
The summary of the data set properties is shown in Table 1.Different thresholds for wavelet coefficients produce different data compressions (relative to the regular 256 3 grid).Completing the original wavelet-based grid to the AMR grid may significantly increase the number of grid nodes for the data set, especially if the original wavelet grid compression is relatively high.Table 1.Data set properties and timings: is the threshold, σ is the compression for the original wavelet based grid, σ amr is the compression for the correspondent AMR grid, F and F vtk are binary file sizes (in megabytes) to store the original wavelet and tetrahedralized AMR grid in VTK format, respectively.T wlt , T PT , T RC , and T ZS are the rendering times (in seconds) for the direct wavelet summation algorithm and for the VTK library unstructured grid renderers: ProjectedTetrahedra, RayCast, and ZSweep, respectively.A hyphen (-) implies the lack of RAM for the rendering.

Hydrogen 3D Orbital Data Sets
An analytical function, electron density for 3D hydrogen orbital, has been used to generate the data set on a regular 256 3 grid.A 4-th order interpolating wavelet transform has been used to compress the data to simulate the output of a wavelet-based PDE solver.The summary of the data set properties is shown in Table 1.

Rendering Results
Rendering times for both data sets are shown in the Table 1 for the direct wavelet summation approach and for three VTK library unstructured grid renderers with VTK's default settings.The projected tetrahedra renderer seems to provide fast volume rendering with adequate quality.Although the best rendering quality is provided by the much slower ZSweep renderer.It should be noted that VTK rendering with the default settings requires a lot of memory (Table 1) and therefore is not appropriate for large data sets.Apparently, ∼5 × 10 7 is the maximum number of tetrahedra to be rendered on our system by VTK library.
VTK rendering through projected tetrahedra gives the rendering time linear to the number of tetrahedra in the unstructured mesh (Figure 5).In the direct wavelet summation approach, the rendering time is linear relative to the number of rays, the average number of points on a ray, and the number of wavelet coefficients in the sum of Equation (1) (Figure 5).Contrary to VTK, the direct summation approach does not require much memory and works successfully for all the data sets (Table 1).The rendering quality of direct wavelet summation on adaptive mesh is comparable to that of VTK's renderers as well as to the quality of Amira [47], volume rendering on a regular 256 3 grid for both convection and hydrogen data sets (Figures 6 and 7).Noticeable color discrepancies in the rendering results are attributed to different ray integration approaches in the different renderers.Some rendering artifacts appear in the images produced by some VTK's renderers (Figure 7).
Increasing the threshold value leads to a smaller number of wavelet coefficients being taken into account during the rendering.The error in the determining the function value f (x) through the truncated version of Equation ( 1) with all wavelet coefficients greater than the threshold is O( ) (e.g., [27]).For both data sets considered, the rendering results appear not to be changed for all the threshold values below 0.01.Therefore, we recommend to use that threshold value as an optimal one for a preprocessing step, before direct summation or tetrahedralization.The estimated errors in the rendering image, ∼1%, seem to be admissible for visualization purposes.Additionally, a larger threshold value strongly decreases the number of tetrahedra in the correspondent AMR grid hence making memory-consuming VTK-based visualization approach feasible.

Conclusions and Future Work
Two visualization techniques have been presented for volume rendering on waveletbased adaptive grids: direct summation of wavelet coefficients during volumetric ray casting, and complementing adaptive wavelet grid to a correspondent AMR grid for the subsequent VTK-based rendering.Both approaches provide relatively good image quality and reasonable timing.The direct summation approach provides a better image quality.In addition, it does not subject to memory restrictions, i.e., if the data set has been computed at some system using a wavelet-based PDE solver, the volume rendering can be performed on the same system.Although, "completing to an AMR" approach seems to be more general in the sense that it moves the problem of volume rendering on a wavelet-based grid into the more developed area of AMR grid visualization.Additionally, VTK provides a large number of tools to make the visualization process, if memory permits, more general.For memory saving, the number of nodes in a wavelet-based grid and therefore the number of tetrahedra in the corresponding AMR grid should be decreased.The obvious approach is to change thresholding.It is also possible to consider a different tetrahedralization technique than the one presented.
As for the summation of wavelet coefficients approach, alternative ways of ray integration could be considered.Particularly, the step along the ray is to be adaptive and the number of rays should not be fixed from the beginning of the rendering.Additionally, direct summation through Equation (1) could be changed to an exact wavelet transform into the corners of the block (of appropriate resolution) around the point of interest on a ray.
Hardware acceleration could be used for wavelet coefficient summation as well as for the coarse mesh interactive rendering.Additionally, coarse mesh interactive rendering could be performed through a faster algorithm, e.g., through shear-warp approach [43] It would allow interactive rendering on a regular grid larger than the current 64 3 grid, which could filter out several levels of wavelet coefficients to sum during the non-interactive phase of the rendering.
These alternative approaches to the volume rendering on wavelet-based adaptive grids are currently under consideration.

Conflicts of Interest:
The authors declare no conflict of interest.

Figure 2 .
Figure 2. Data structure used for volume rendering.Scaling coefficients are stored as N 3 array of "Coarse Mesh Nodes" (N = 64 in the current implementation).Wavelet coefficients inside a coarse mesh block are stored in a dynamically allocated array of "Wavelet Nodes".

Figure 3 .
Figure 3. Completing wavelet based grid (filled circles) to AMR in two-dimensional case by adding the required nodes (white circles).Total number of cells in the final AMR grid is 23.

Figure 5 .Figure 6 .
Figure 5. Direct summation rendering time T wlt vs. data compression σ for the convection (a) and electron density (c) data sets.The coarse mesh is 64 3 , the number of rays is 500 2 , the average number of points per ray is 170.The cutoff limits for the scaling and wavelet functions are 1 and 2, respectively.VTK projected tetrahedra rendering time T PT vs. the number of tetrahedra in the correspondent AMR grid for for the Rayleigh-Benard convection (b) and electron density (d) data sets.

Funding:
D.A.Y. was partially supported by the DOE under grant DE-SC0019759 and by the NSF under grant EAR-1918126.