Beyond ExaBricks: GPU Volume Path Tracing of AMR Data

Adaptive Mesh Refinement (AMR) data is currently emerging as one of the most prevalent and important types of data that any visualization-focused renderer needs to be able to handle. Unfortunately, from a rendering point of view AMR data comes with several challenges: first, AMR data can come in many forms, from octree-AMR [BWG11] to various forms of block-structured AMR [CGL^∗00]. Second, for rendering AMR data one requires a sample reconstruction method to be defined over the AMR cells’ scalar values; and how this is defined can have a huge influence on the images that the renderer produces. Third, most AMR data is highly non-uniform in nature, requiring the renderer to be able to concentrate its sampling effort where it matters most. Fourth, rendering AMR data is expensive partly because of the aforementioned issue of non-uniformity of the signal, and partly because of the hierarchical nature of the data. This means that each sample typically requires some form of hierarchical data structure traversal. Fifth, the recent advances in real-time ray- and path-tracing have raised the stakes for sci-vis rendering, too. While arguably, effects such as indirect illumination or soft shadows and ambient lighting are most important for photorealism, demand for high-quality Monte Carlo rendering is gradually becoming the de facto for sci-vis rendering as well. Particularly, this change in paradigm has been initiated by industry-quality tools like OSPRay [WJA^∗17] or OpenVKL [KJM21], which are readily available to sci-vis users through ParaView or VisIt plugins [WUP^∗18].

We explore the problem of high-quality path traced rendering of non-trivial AMR data sets (cf. Beyond ExaBricks: GPU Volume Path Tracing of AMR Data). For that we start with the ExaBrick AMR rendering framework by Wald et al. [WZU^∗21], and look at what is required to extend this to volumetric path tracing. ExaBrick is both a framework and an acceleration data structure, similar to what kd-trees and BVHs are for surface ray tracing, and while fundamentally, there is no reason this acceleration data structure should not work with volumetric path tracing, it was not originally optimized for that. It is an open research question if this data structure optimized for sci-vis style volume ray marching can be easily adapted to support the fundamentally different operation of computing transmittance estimates as performed by a volumetric path tracer.

In contrast to ray marching, the dominating operation in volume path tracing is computing free flight distances, which require an upper bound estimate for the extinction of the volume density, called majorant. Though the method is correct irrespective of how much that upper bound over-shoots the real value, too high a difference between conservative majorant and actual extinction that varies in space causes unnecessary samples to be taken, so that key to accelerating volume path tracing is finding majorants that vary in space themselves and then form an acceleration structure.

In sci-vis, however–and in particular, for the type of data we are looking at in this paper–computing good majorants is hard: one factor is that actual extinction at any point in space depends on a transfer function that can map any data value to any extinction value, and which can—and typically will—undergo radical changes as the user explores the volume. A second factor is that the original AMR structure wants to adapt to the data, so that even if we were hypothetically allowed to relax the real-time requirement, finding good alpha-mapped majorants is a hard problem in itself.

Building upon these challenges in AMR visualization, the key contributions of this paper are the following:

• •

  We extend the AMR framework ExaBrick by replacing alpha-composited ray marching with volumetric path tracing, giving rise to significantly improved time to image for sci-vis rendering, plus the possibility to interactively create images with much better visual fidelity by using full global illumination.

• •

  We investigate the question how efficient the existing ExaBrick data structure is by systematically identifying and implementing all of the most obvious alternative variations of that approach.

• •

  We thoroughly analyze and compare the different alternatives and discuss their technical implications. This includes a comparison to a hypothetical (and due to its construction times impractical) reference data structure that is allowed to perform near unlimited preprocessing.

Volume rendering is important for both scientific visualization (sci-vis) and production rendering, but opposing goals regarding visual fidelity on the one hand, and interactivity on the other hand, have led the two fields to take largely orthogonal approaches. Even early production rendering systems used scattering and virtual point lights [LSK^∗07] and have since transitioned to full global illumination with multi scattering and hundreds of bounces [FHP^∗18]. Scientific visualization has traditionally focused on interactive exploration [HKRs^∗06] based on ray marching [Lev88] with absorption and emission [Max95]. Rendering with this lighting model can still be found in common scientific visualization systems like VisIt [CBW^∗12] or ParaView [AGL05].

Exposure renderer [KPB12] is one of the first examples to use volumetric path tracing for sci-vis. It supports multi-scattering and can produce highly realistic images [MB17, IZM17]. However, volumetric path tracing only started to gain more attention from the community recently when optimized software and hardware ray tracing frameworks became widely available [WJA^∗17, IGMM22]. Sci-vis renderers use simple isotropic phase functions covering the whole medium, RGB albedo from color lookup tables as opposed to full spectral rendering, and simple lighting models. Common techniques to compute free-flight distance estimates—the stochastic distance a photon travels without colliding with a particle from the density—are unbiased, which is attractive also for scientific visualization [MSG^∗23].

A number of techniques exist to compute free flight distances and transmittance estimates using stochastic sampling. Woodcock (or delta-) tracking [WMPT65] is one of the oldest methods and virtually “homogenizes” the volume by introducing fictitious particles that produce null collisions where the ray is extended without a change in direction. Other estimators include decomposition tracking [KHLN17], residual ratio tracking [NSJ14], or unbiased ray marching using power series expansion [KDPN21].

Combined with spatially varying local majorants [NSJ14], unbiased transmittance and radiance estimates can be found much faster than with quadrature techniques, which is attractive for interactive exploration. In this sense, local majorants allow for the equivalent of empty space skipping and adaptive sampling, which are popular acceleration techniques in scientific volume rendering [ZSL19, ZSL21, WZM21]. In contrast to approaches that approximate the local frequency using the number of local cells and their sizes [MUWP19], majorants directly consider the actual frequency of the, which is represented by the local extinction.

Commonly used data structures that allow for traversal with local majorants are kd-trees [YIC^∗10] as well as uniform macrocell grids that are traversed using the 3D digital differential analyzer (DDA) algorithm [SKTM11]. Bounding volume hierarchies (BVH) have received less attention as direct volume rendering accelerators, as the predominant techniques in sci-vis traditionally required ordered traversal. Morrical et al. [MSG^∗23] used BVHs to accelerate an interactive volume path tracer, but concluded that a simple DDA grid implementation was more efficient. Approaches using BVHs as empty space skipping accelerators for volume rendering with front-to-back alpha compositing required the bounding boxes at the leaf nodes to not overlap [ZHL19]. The OpenVDB data structure that is popular in production rendering uses a hierarchical grid and DDA traversal with a fixed hierarchy size [Mus13]. Hofmann et al. [HE21] also provide a recent summarization of multi-level DDA traversal for volumetric path tracing on GPUs.

The term adaptive mesh refinement (AMR) was originally coined by Berger and Collela [BO84, BC89] at NASA. It refers to simulation codes that hierarchically subdivide a (often uniform) grid, both in space and in time, in regions where the simulation domain is more interesting than in others. Subdivision schemes include octrees [BWG11], or block-structured AMR [CGL^∗00] where the rectangular regions are allowed to be irregular.

Virtually all existing AMR simulation codes (e.g., FLASH [DFG^∗08], Lava [KBH^∗14], etc.) output data that is cell-centric and results in T-junctions at odd level boundaries. Work by Wald et al. [WBUK17] and by Wang et al. [WWW^∗19] has proposed interpolators that are continuous even at level boundaries, which is important for artifact-free visualization. These interpolators use sampling routines that perform neighbor queries by traversing tree data structures. If not carefully designed they require traversing those data structures per sample taken, resulting in non-trivial reconstruction costs infeasible for GPUs. ExaBrick by Wald et al. [WZU^∗21], which is explained in more detail in Section 3.2, provides an acceleration structure to accomplish fast AMR cell location on GPUs.

For rendering octree-based AMR specifically, Wang et al. [WMU^∗20] developed high-quality reconstruction filters for direct ray tracing. Other authors have also focused on large-scale out-of-core AMR rendering, such as the interactive streaming and caching framework proposed by Wu et al. [WDM22] targeted at CPU rendering, or the framework by Zellmann et al. [ZWS^∗22], who proposed a similar method for AMR streaming and rendering on GPUs, but focused on time-varying data.

GPU path tracing native AMR data, to our knowledge, has not been proposed by the literature yet. Closest to this are OpenVKL [KJM21] that can path trace on the CPU but does not use any of the optimizations necessary for GPUs, and the framework proposed by Zellmann et al. [ZWMW23] that renders the AMR data by converting it to an unstructured mesh with voxels.

In this section we discuss how to integrate volumetric path tracing (VPT) into a framework like Wald et al.’s ExaBrick [WZU^∗21], which can render AMR data natively. A key requirement of our method is that VPT is compatible with RGB$\alpha$ transfer functions that contribute albedo and extinction coefficients.

We implement volumetric path tracing (VPT) with OptiX 7. We support three visualization modes with different quality settings: absorption plus emission, realized with Woodcock tracking; single scattering, where an additional bounce to a light source chosen at random is performed to compute volumetric shadows; and multi scattering using a global phase function that allows the light to bounce arbitrarily. The different quality settings can have a tremendous impact on visual fidelity, as can be seen in Fig. 3.

Inside the OptiX ray generation program, we generate cameras ray with random jitter to compute screen space samples. For the tracking estimator implementation, the rays traverse an acceleration structure (as previously described); at each sampling position they compute the extinction coefficient $\mu(x)$ by first sampling into the volume (using an AMR sample acceleration structure), and then post-classifying the sample using the RGB$\alpha$ transfer function.

Given those two operations, which are described and evaluated in detail below—segment traversal and volume sampling—free flight distances are computed that are required by all three visualization modes. For the absorption and emission model, the sample color is proportional to the albedo from the transfer function at the position where the collision occurred. The single scattering estimator casts a shadow feeler towards a randomly chose light source, and the albedo gets weighted by the transmission between sample position and light source that we compute using ratio tracking [NSJ14]. In the case of multi scattering we evaluate a Henyey Greenstein phase function and perform bounces in a while loop inside the ray generation program. Based on the throughput, we terminate paths at random using Russian roulette.

Since we interpret the RGB component of the RGB$\alpha$ transfer function as albedo, the path throughput is an RGB tuple as well—and so is the final screen space sample. We accumulate the screen space samples in a device-side accumulation buffer using a weighted average. Using an interactive prototypical viewer, whenever the camera changes, we reset accumulation. The accumulation buffer content is periodically presented as sRGB color tone mapped and drawn into an OpenGL pixelbuffer object mapped using CUDA/GL interop.

We describe how to implement the two core routines for VPT—ray segment traversal and AMR data sampling—in the following sections. The most obvious way to implement those operations is by using the existing ABRs (cf. Section 3.2). In the following sections, we explain the computation of ABR majorants and introduce alternative methods for spatial and object subdivision, which facilitate traversal and sampling of the AMR data structure.

Volume acceleration data structures partition the ray into segments $[t_{0i},t_{1i}]$ that each have their own majorant $\bar{\mu}_{i}$. Inside the segments, the rays are sampled. In this section we describe different strategies to implement those two operations. A restriction of hardware ray tracing is that when a ray traverses a BVH and the intersection or hit program is called, tracing rays against another hardware BVH is not permitted. It is thus not valid to use one BVH to traverse segments, and while integrating a segment traversing another BVH for sample reconstruction.

Alternatives to using two different hardware accelerated BVHs are using software acceleration structures for either of the operations; using the same acceleration structure for traversal and sampling; or halting the (outer) ray segment traversal for traversal, and then restarting it by setting the ray’s $t_{min}$ parameter to the previous segment’s $t_{max}$. Each strategy we propose implements one of those alternatives.

For the evaluation we concentrate on runtime and memory performance of the different combinations of acceleration structures. Our goal is to compare the performance of two segment traversal structures: one adaptable to transfer function changes (sections Section 5.1.1, Section 5.1.2, Section 5.1.3), and one optimized for a single transfer function (Section 5.1.4). Another question we explore is if spatial acceleration structures (grids, kd-trees) have an advantage over object structures (BVHs)—or the other way around, when being used for ray segment traversal; we are also interested in if, for this kind of volume data it is beneficial to use the same acceleration structures for traversal and sampling. And finally, we want to explore if a universal combination of acceleration structures exists, or if the choice is predominantly dependent on the data set and transfer function.

We compare the various combinations of sampling methods (ABR BVH vs. extended brick BVH), spatial subdivion/majorant traversal data structures (extended bricks, brick bounds, macrocell grid), and implementation using ray tracing hardware (RTX BVH) vs. software DDA, against the baseline of using ABRs for both traversal and sampling. We use the data sets from Table 1, which also shows majorants for the acceleration structures.

Our target system for the performance evaluation is a GPU workstation with an NVIDIA GeForce RTX 4090 GPU (24 GB GDDR memory) and an Intel i7-6800K CPU and 128 GB main memory on the host. Our evaluation concentrates on memory consumption and rendering performance. We compute full global illumination with multi-scattering (isotropic phase function, russian roulette path termination, ambient light surrounding the whole scene). We use two transfer functions—one that generates a high-opacity visualization (cf. first column in Table 1, “spiky”), and a second that results in a more translucent appearance (“foggy”). For the grid benchmark, we first empirically determined optimal grid dimensions, the result of which is presented in Fig. 7. Based on that we chose $128^{3}$ (cloud), $128\times 128\times 64$ (meteor-20k, meteor-46k), $512\times 512\times 256$ (gear), and $512\times 256\times 384$ (exajet) for the remaining experiments. The performance results in frames/sec. can be found in Table 2. We also report peak GPU memory consumption and the total memory consumption (i.e., during rendering, when acceleration structures were built and temporary memory freed) in Table 3.

Our results allow for multiple observations. One such observation is that data structures with overlap (ABRs, bricks, brick domains) can suffer from the aforementioned problem of single cells infecting large regions with their contribution. This can also be seen in Table 1, where the majorant regions for ABRs or bricks sometimes appear excessively coarse. This observation suggests that spatial data structures are to be preferred over object data structures as this problem is easier to control. However, it is well-known that grids suffer from the “teapot in a stadium” problem. To illustrate this we include the landing gear data set whose cells are contained within a large boundary of coarser cells that do not contribute important features (cf. Fig. 8). Apart from this extreme case, we however note that the proposed traversal and sampling methods come close to and often even outperform the majorant of a laborously pre-computed kd-tree.

We also observed a fixed cost associated with using software traversal with DDA, which we suspect is due to register pressure (OptiX however does not allow us to exactly determine per-kernel register usage). This becomes obvious in the comparison between DDA and RTX grid traversal for small grid sizes (e.g., $\leq 10^{3}$ macrocells, cf. Fig. 7). We however also note that the memory consumption for this option can become quite high; in the special case of a transfer function without empty space, we have to build a full BVH over the grid, which limits the number of grid cells we can use effectively.

We have shown how to take a visualization framework targeted at large-scale AMR data and modify it to support interactive volumetric path tracing. The benefits of this are increased visual fidelity and improved time to first image. Volumetric path tracing seamlessly integrates complex lighting models. It also supports traditional models if the user desires not to focus on realism. Even in these cases, we still find path tracing preferable, as it enables smarter sample placement compared to ray marching, resulting in fewer samples needed for volume-rendered images. We contribute the first framework plus acceleration structures optimized for volume path tracing native AMR data on GPUs with ray tracing hardware.

The base operation performed by a volumetric path tracer is free flight distance tracking, which requires routines to obtain extinction and majorant estimates for local regions. The local majorant extinctions to traverse the volume with the minimum number of samples usually have an irregular spatial arrangement. AMR hierarchies already focus their less uniform structures in regions where the entropy of the data is high and consequently the extinction varies a lot. One might therefore intuitively assume that the AMR hierarchy would be the best possible majorant traversal data structure. However, our experiments show that the task of finding majorants for this kind of data is more complicated.

We find it hard to make a general recommendation as of the best combination of acceleration structures for traversal and sampling, yet we observe a tendency for object order data structures to do a bad job in this case because their overlap allows small regions from neighboring structures to “infect” other structures that are otherwise homogeneous or even empty. Notably, both combinations that use the same acceleration structure for traversal and sampling generally do not outperform the other alternatives.

In fact, grids often present a viable alternative; in this case the outer loop traversal over ray segments is implemented in software and the inner loop uses hardware ray tracing. We find that using hardware ray tracing to traverse the grid’s macrocells pays off for large grids, though this is quite memory intensive, restricting the maximum grid size one can use.

The solutions we focus on allow us to update RGB$\alpha$ transfer functions interactively, even for our largest data sets. The alternative to this is using pre-classified data structures, which due to the nature of how extinction is computed via alpha transfer function lookups often even results in inferior acceleration structures, and can take hours to build even when using optimizations such as priority queues and binning.

Another observation we made is that transfer functions that have high frequencies—as is often the case when extracting ISO-surface like features—can result in extreme performance degredations. Adjusting the slope of “alpha peaks” can lead to extremely “spiky” transfer functions, causing high sample counts and significant frame rate reductions.

Another point important to discuss is that our solution accounts for the (pre-tabulated) albedo by multiplying it to the result of the evaluated phase function. Because of that we cannot use residual tracking methods [KHLN17, NSJ14] as the control variates cannot vary in space, and reconstructing the correct albedo for a control variate sample is impossible. To this end, the solution employed in production rendering would be to have individual majorants plus control volumes for each R,G, and B channel of the volume. In our case, this would limit memory availability even further, so we opt to use larger grids, or deeper hierarchies with more leaves, in favor of residual tracking algorithms, which we believe gives us better overall performance. A more thorough investigation of this would however be interesting future work.

In conclusion, we have presented how to extend a sci-vis volume renderer with non-trivial sampling structures to support volumetric path tracing. Volumetric path tracing can significantly improve visual fidelity, but even with standard sci-vis lighting models has several advantages, such as better sample placement—and resulting from that faster time to image—as well as unbiasedness and fewer sampling artifacts. We have evaluated multiple alternatives to implement this by proposing extensions to the ExaBrick data structure; our experiments show that computing majorant extinctions for this kind of data is not trivial—not even in the case where a single hand-picked transfer function is used, and even less so when transfer functions are allowed to change interactively. Our results indicate a slight performance advantage for spatial over object data structures. Overall we have presented numerous different strategies to accelerate traversal and sampling that allow for interactive transfer function updates and demonstrated that these do not perform significantly worse—and sometimes even outperform—the offline reference.

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—grant no. 456842964. We also express our gratitude to NVIDIA, who kindly provided us with hardware we used for the evaluation.