{'fraction_non_alphanumeric': 0.030853345182666487, 'fraction_numerical': 0.006916404242576457, 'mean_word_length': 4.927851580393953, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning Multiple-Scattering Solutions for Sphere-Tracing of Volumetric Subsurface Effects

Accurate subsurface scattering solutions require the integration of optical material properties along many complicated light paths. We present a method that learns a simple geometric approximation of random paths in a homogeneous volume of translucent material. The generated representation allows determining the absorption along the path as well as a direct lighting contribution, which is representative of all scattering events along the path. A sequence of conditional variational auto-encoders (CVAEs) is trained to model the statistical distribution of the photon paths inside a spherical region in presence of multiple scattering events. A first CVAE learns to sample the number of scattering events, occurring on a ray path inside the sphere, which effectively determines the probability of the ray being absorbed. Conditioned on this, a second model predicts the exit position and direction of the light particle. Finally, a third model generates a representative sample of photon position and direction along the path, which is used to approximate the contribution of direct illumination due to in-scattering. To accelerate the tracing of the light path through the volumetric medium toward the solid boundary, we employ a sphere-tracing strategy that considers the light absorption and is able to perform statistically accurate next-event estimation. We demonstrate efficient learning using shallow networks of only three layers and no more than 16 nodes. In combination with a GPU shader that evaluates the CVAEs' predictions, performance gains can be demonstrated for a variety of different scenarios. A quality evaluation analyzes the approximation error that is introduced by the data-driven scattering simulation and sheds light on the major sources of error in the accelerated path tracing process.


INTRODUCTION
The realistic rendering of translucent solid materials such as wax or skin requires simulating the scattering of light in the interior of the body. Subsurface scattering refers to the mechanism where light that penetrates the boundary of a translucent object is scattered internally until it exits the body at a different location on the boundary. Due to the internal scattering of light and the resulting complicated light paths that have to be considered in the simulation, an accurate subsurface scattering solution becomes computationally very expensive.
To reduce the computational requirements, approximation methods exploit special shape and material configurations, such as planar surface geometry, isotropic scattering, or separability of the bidirectional scattering-surface reflectance distribution function (BSSRDF) into location and direction dependent terms. While it has been shown that efficient analytic transport solutions can be developed for some special cases, rendering quality decreases when the resulting algorithms are applied to scene configurations beyond the scope of the restrictive assumptions. To improve the generalizability of rendering arXiv:2011.03082v1 [cs.GR] 5 Nov 2020 acceleration solutions, there has been growing interest in data-driven approaches that can infer approximate multiple-scattering solutions at render time.
One of the earliest methods precomputes long-distance light transport through a homogeneous medium using transition probabilities at varying spatial scales [33]. This so-called shell tracing approach enables to move rays through the volume in large steps without considering individual scattering events. It builds upon the concept of sphere-tracing for surface rendering [15], to adapt the step size to how deep the ray is in the solid body. Recently, for rendering volumetric media without an explicit boundary representation, deep-scattering [21] proposed to fed geometric information into a network by progressively evaluating a hierarchical volumetric geometry descriptor. Inference is based on an expressive multi-layer neural network. To address specifically the issue of boundary effects, a method for learning a shape-aware BSSRDF model from ground truth volumetric transport simulations using path tracing has been proposed by Vicini et al [47]. The method avoids many limiting assumptions of prior analytic models regarding the planarity of surfaces and isotropy of volumetric transport, and it generalizes the built-in notion of spatio-directional separability. It assumes a local polynomial surface expansion around a shading point, and uses this expansion to re-project predicted off-surface samples back onto the surface.
Our work builds upon prior works in data-driven subsurface scattering simulation, by leveraging the core idea of "bypassing a potentially lengthy internal scattering process" [47] through a network that has learned to sample outgoing light locations on an object's surface from an incident light location. In contrast to existing work in subsurface scattering simulation, our method avoids an explicit encoding of the local surface geometry in the stochastic modelling of the scattering process, and it also infers the outgoing light direction from a reference distribution generated via volumetric path tracing, instead of using importance sampling. To achieve this, we combine shell tracing with network-based learning of long-distance light transport between a point and a spherical surface centered at this point. A sequence of neural-network-based conditional variational auto-encoders (CVAEs) is trained to model the statistical distribution of light paths and absorption along them inside a spherical region in presence of multiple scattering events. A first CVAE learns to sample the number of scattering events, i.e., the absorption, occurring along a path inside the sphere. Conditioned on this, a second model predicts exit position and continuation direction of this path. By using shell tracing, light paths can be traced efficiently through the volumetric medium toward the solid boundary without an explicit encoding of the boundary geometry.
Furthermore, instead of modelling purely the one-to-one light transport between two surface points, our method supports extensions of volumetric path tracing, including direct in-scattering to achieve improved convergence. For this, a third network model generates a representative sample of photon position and direction along the path, which is used to approximate the contribution of direct illumination due to in-scattering. In combination with existing approaches for simulating the direct interior illumination, the path-tracing solution converges quickly even when none of the paths arrives directly at a light source.
We jointly train all three networks on a spherical geometry with unitary radius and varying interior density (different radii can be considered by adjusting the density appropriately). In the training and validation phase, ground truth transport solutions using Monte Carlo path tracing are used. An interesting result is that all proposed networks are light-weight, comprised of only up to three layers with no more than 16 nodes each. Despite the resulting compact internal latent-space representations, very realistic reconstructions of the modeled distributions are obtained. Due to this compactness, all inferencing steps can be implemented efficiently using short shader programs and matrix multiplication instructions on the GPU. We analyze the performance and accuracy of our approach on different geometry and material properties, and we also investigate the inferencing skills of each of our network models individually. Our method achieves quality en par with that of a reference Monte Carlo path tracing solution, at significantly improved performance.

RELATED WORK
Compared to traditional ray tracing algorithms that can assume a free traverse of the light on the vacuum between surfaces, ray tracing participating media represents a superior level in computational complexity. For some medium such as milk and clouds, the high value of the scattering albedo (i.e. the ratio of energy remaining after scattering) entails a large number of light bounce simulations before the ray is absorbed, leaves the volume or hits a surface. Due to the plethora of possible material interactions, the rendering equation for realistic scene settings does not possess a closed-form solution in presence of participating media. Thus, for synthesizing realistic images, approximate solutions are required.
Traditional approximation schemes include the exclusion of higherorder scattering effects, leading to the single-scattering approximation [46] -which under certain circumstances can admit a closed-form solution [40] -or approaches based on luminance diffusion [9,20]. Especially approximations of the latter type, however, typically impose severe assumptions on homogeneity of the considered material and surface geometry, which effectively disables application of these approaches in more complex settings. Even increasing model capacity, for example by considering dipole diffusion [19], does not resolve the shape restrictions. A more thorough review of traditional rendering techniques for participating media can be found in [3] Monte Carlo methods More flexible approaches arise from the application of Monte Carlo path tracing algorithms [28] for numerically solving the rendering equation for participating media [14,20,44]. Though being computationally expensive, Monte Carlo path tracing in a converged setting yields physically accurate images which can serve as a "ground-truth" for validating alternative approaches. The recent survey by Novak et al. [36] reviews the latest advances in Monte Carlo path tracing methods for solving the light transport in participating media.
In path tracing and as an alternative to ray-marching with constant step size, delta tracking is often used to determine free path lengths according to the optical depth in the participating medium. As an importance sampling of the corresponding probability density function, Woodcock tracking [51] adjusts the sampling distances so that dense regions in a volume can be sampled appropriately, and it has been adapted to achieve unbiased sampling of sparse inhomogeneous media [53]. Free path sampling with probabilities not necessarily proportional to the volume transmittance has been realized by means of weighted delta tracking approaches, e.g. [25,38,43]. All these techniques reduce noise in the estimated light paths for participating media, and they enable an automatic selection of the step size according to the underlying material distribution.
A different approach to move rays through the volume in adaptive steps even without considering individual scattering events is shell tracing [33]. It builds upon a statistical model that is based on a discretized array of transition probabilities, limiting the expressiveness of the simulated light transport paths. In the spirit of sphere-tracing for surface rendering [15], it performs larger step in the interior of a body and successively decreases the step size towards the object boundary.
Diffusion methods Despite the large generality of Monte Carlo Rendering, a variety of methods have been developed, which are based on diffusion approximations. Especially in the field of subsurface scattering, computation times can explode when highlydense materials require translucent renderings. Though the light does not enter the bodies deeply, Monte-Carlo approaches may run into lengthy sequences of frequent scattering or reflection events right underneath the surface. Diffusion models, which approximate the light transport by virtue of bidirectional scattering-surface reflectance distribution functions (BSSRDFs) may in these cases lead to speed-ups. While early approaches, like [19], can hardly keep up with the quality of ray-traced images, more recent methods like [7,10,52] can generate visually appealing results at reduced rendering time, however, at the cost of reduced physical realism. Notably, also hybrid methods exist, which weight out the advantages and disadvantages of Monte Carlo-based and diffusion-based rendering against each other [12].
Machine Learning for Rendering As Monte Carlo-based rendering can intrinsically be described as a statistical problem, the use of machine learning methods for rendering has recently attracted great interest. In the context of path tracing the goal is typically to use machine learning models to model sampling distributions or create guidance for efficient sample selection. Popular modeling approaches include Gaussian mixture models [17,48], and neural network methods.
Neural network-based methods In the context of path tracing, deep learning methods can be classified into image-based adaptive sampling and denoising algorithms on the one hand, and multiplescattering approximations on the other hand.
The recent approach by Kuznetsov et al. [26] facilitates neural networks for learning adaptivity in Monte-Carlo path-tracing and denoising of the final image. A first network learns to adapt the number of additional paths from an initial image at the target resolution, which is generated via one path per pixel. A second denoising network learns to model the relationship between an image with increased variance in the color samples to the ground truth rendering [31,42]. Recently, Weiss et al. [50] used neural networks to learn the positions of sample locations in image space from a low resolution image, using a differentiable sampling stage as well as a differentiable image reconstruction stage that can work on a sparse set of samples. Related to denoising approaches for path traced images is the screen space shading approach by Nalbach et al. [35], where a network is trained to infer missing shading information from image information in a view-dependent G-buffer.
For surface graphics, Zheng and Zwicker have used neural networks to model the relationships between scene parameters and light paths [55]. Recently, deep-scattering [21] employed a radiancepredicting neural network for simulating scattering events in clouds. Geometry information is fed into the network by progressively evaluating a hierarchical volumetric geometry descriptor, and inference is based on an expressive multi-layer neural network, optimized with recent design features. Nevertheless, the method was developed for cloud rendering and does, as such, not explicitly account for boundary effects, which may play a major role in sub-surface scattering. Adrressing specifically this issue, a method for learning a shape-aware BSSRDF model from ground truth volumetric transport simulations using path tracing has been proposed [47]. Exploiting the expressiveness of neural network models in a similar way, this method avoids many limiting assumptions of prior analytic models regarding the planarity of surfaces and isotropy of volumetric transport, and it generalizes the built-in notion of spatio-directional separability. It assumes a locally expandable surface geometry, which is explicitly encoded using a trivariate polynomial around a shading point. In this way, an approximate signed-distance function to the surface geometry can be used to model the outgoing on-surface light distribution.
In [47], a learned approximation of all the subsurface scattering allows them to overcome the limitation of the planar approximation, fitting the geometry with a cubic polynomial. Nevertheless, their assumption of the final outgoing direction with a uniform distribution and a single polynomial approximation to the whole geometry biases the final result.
Another idea, considering the light scattering in a medium as a random process (similar to the walk-on-sphere statistic technique [34] [27]) is to have a short-cut of all possible scattering inside a sphere with homogeneous density. This model can be used to simulate all those "events" as a single step in a Montecarlo technique. This strategy represents an acceleration rather than an approximation, and removes any bias regarding the geometry, limitation of other techniques. In [33] they propose a method named "shell-tracing" dealing with the rendering of discrete random media using precomputed scattering distributions in form of matrix factorization. With this approach, the number of variables considered is reduced to three due to the increment of the matrix complexity if a higher dimension is considered.

BACKGROUND AND METHODS
Our key idea is to train generative statistical models to bypass multiple scattering events occurring inside a spherical volume of material with constant density and potentially anisotropic scattering characteristics, and to use the learned representations in the rendering process to generate paths through the volume using large steps (Fig. 2a). Whenever the ray arrives at a sampling point, it tests for the largest step it can make without leaving the volume, and then evaluates a sequence of statistical models to infer the end position and direction when performing this step. The ray proceeds using delta tracking until a next collision event is determined, and recurrently performs learning-based path tracing (Fig. 2b). To determine the length of Figure 2: Method overview. a) Given optical properties of the volumetric medium ( , , ), the path (gray dotted lines) starting at 1 = (0, 0, 0) from the incoming direction 1 = in is bypassed by letting a network infer directly (blue) the random outgoing position and direction out in an unitary radius sphere. b) After free-travel events (black arrows), from a signed distance field the maximum sphere not intersecting the geometry is computed, and the radius is used to control the adaptive density and the length of the next step along a light path. This process is repeatedly applied until the path leaves the volume.
the step, a 3D signed distance function is computed for the current geometry in a pre-process. For this, the object's bounding volume is discretized using a 3D voxel grid, and for each voxel a conservative shortest Chebyshev distance to the object's surface is computed [6]. This representation enables to quickly find the radius of the sphere that can be safely traversed without intersecting the geometry. As the ray approaches the edges, it is forced to take smaller and smaller steps and eventually the tracking will pass across the boundary surface. To efficiently trace against a polygonal boundary surface, we employ a GPU ray-tracer using NVIDIA's RTX ray tracing interface [39] through the DirectX Raytracing API (DXR).

Scattering model and data generation
Learning is conducted in an offline process. We assume homogeneous extinction t and scattering albedo (with t = a + s considering the rate of photons getting absorbed a/ t or scattered into other directions = s/ t due to light-particle interaction). Anisotropy in scattering is modelled by using the Henyey-Greenstein [16] phase function, with cos(·, ·) indicating the cosine between the respective directions and anisotropy parameter , which equals the expected cosine between ray direction before and after scattering. For instance, > 8 results in a strong forward lobe of scattering into the incoming light direction. In this setting, the radiative transport [4] can be described in terms of the volume rendering equation [20], ( , ) describes the incoming radiance at position from direction direction , e ( , ), s ( , ) and d ( , ) denote emitted, scattered and direct radiance, and ( , ) is the transmittance along the ray (Beer-Lambert law [29]), with ( , ) = exp(− t ) in the case of homogeneous medium. The impact of the phase function is considered in s ( , ), which can be rewritten as where i ( , ′ ) denotes the directional incoming radiance at , and the integral runs over the full solid angle Ω. For a detailed discussion of the derivation of the rendering equation, we refer the reader to [41].
For training our model, we assume that rays start in the center of a sphere of unitary radius, and terminate when they cross the boundary of the sphere for the first time, after experiencing a sequence of scattering events in between. By assumption, the incident direction of the rays in the sphere center is set to CVAE in = (0, 0, 1) for all rays, and the first scattering event occurs directly in the center of the sphere. For each ray, Monte Carlo path tracing is used as a ground truth renderer to simulate radiative transport through scattering media [28,36]. The distance between subsequent scattering events is estimated as = − log(1 − )/ t , wherein is a uniformly distributed random number between 0 and 1. I.e., the mean free path is equal to the reciprocal of the extinction. At every scattering event, it is decided whether the path continues or terminates due to absorption using importance-based Russian roulette [8]. In case of continuation, the scattering direction is computed by importance sampling the Henyey-Greenstein phase function.
Within our model, we ignore the details of the full path trace. Instead, we assume that -in the simplest case of a standard path tracer -the only states which affect the outcome of the rendering are the number of scattering events affecting a ray inside the sphere, determining the likelihood of absorption (ray terminates due to absorption, = 1, or not, = 0), as well as the exit position on the shell of the spherical volume and the outgoing direction of the ray. To further support the use of the proposed method in the more general context of path tracing, we include the opportunity to sample additional ray statistics, which we summarize as Σ.
Following up these considerations, the likelihood of sampling a particular ray path with summary statistics ( , , , , Σ) is determined by a conditional probability density function (PDF) of the form ( , , , , Σ| , , ), wherein quantities to the right of the vertical bar affect the PDF only as conditions.

Model architectures and inference
Empowered by fast progress in machine learning, in particular deep learning, a variety of approaches have been proposed to efficiently learn complex statistical distributions. These include variational autoencoders [22], generative adversarial networks [11] or invertible neural networks [1], which all can be extended to learn conditional distributions, see e.g. [2,32,45]. Within this study, we propose to train conditional variational auto-encoders (CVAEs) to learn a representation of ( , , , , Σ| , , ), which offer a preferable trade-off between computational model complexity, model flexibility and training stability.
In general, a CVAE defines a directed graphical model, which learns to emulate complex conditional data distributions, X ( | ), with data ∈ X and condition ∈ C, by establishing a mapping between the original data space X and a latent feature space Z. In this feature space, the data samples follow a prior distribution, Z ( | ), ∈ Z, from which samples can be drawn more easily, such as a multi-variate Gaussian distribution N (0, I Z ) with zeromean and diagonal unit variance. The distributions are related via which can be combined to yield a self-supervised learning model, similar to the standard auto-encoder. Therein Enc and Dec denote probabilistic encoding and decoding mappings between data and latent space. Exploiting these relations, the CVAE defines a pair of neural network models, the encoder Enc( , ) and the decoder Dec( , ), which learn to parameterize the encoding and decoding PDFs Enc and Dec . While the data-space posterior Dec is typically considered to reflect the true mapping between Z and X, the latent-space posterior Enc can only be treated approximately, since the marginalization required to derive the exact latent-space posterior is computationally unfeasible. Within our study, we assume a conditional Gaussian shape for both distributions, i.e., wherein Enc ( | , ) refers to the parametric approximation of Enc ( | , ). Furthermore, Enc ( , ) and Σ Enc ( , ), as well as Dec ( , ) and Σ Dec ( , ) are outputs of the encoder model Enc( , ), and decoder model Dec( , ), respectively, with Enc ( , ) ∈ Z, Dec ( , ) ∈ X. Σ Enc and Σ Dec are diagonal covariance matrices, consistent with the dimensions of Z and X, respectively. At training time, encoder and decoder models are trained jointly to maximize the so-called evidence lower bound which constitutes a lower bound for the observed data log-likelihood [22]. The first term refers to Kullback-Leibler divergence [24] between the learned latent-space posterior and the latent-space prior, and penalizes a mismatch between both distributions. The second term denotes the expected log-likelihood of the learned data posterior, where the expectation is taken with respect to the latent-space posterior, and penalizes the reconstruction error of the model. At inference time, the encoder model is discarded and a latent-space sample is drawn from the latent-space prior distribution N (0, I Z ). Conditioned on this sample, the decoder model is evaluated and a data-space sample is obtained as where is another Gaussian random variable, drawn from N (0, I X ). This so-called reparameterization trick [22] assures that ∼ N ( Dec , Σ Dec ).
A sketch of the training and inference processes in the current setting is shown in Figure 3a.
In our work, we realize CVAE models via pairs of neural networks, encoder and decoder, respectively. During rendering, only the decoder is evaluated and, thus, needs to be kept as simple as possible.  , the latent space dimension , the depth of the networks , and the number of nodes per layer . Encoder and decoder share the same settings of , and , and output estimates for mean Enc/Dec and diagonal covariance matrix Σ Enc/Dec of the respectively modeled distributions.
All our models follow the scheme shown in Figure 3b and are determined through the number of input parameters in , the number of output parameters, out , the latent space dimension , the number of layers, i.e., the depth of the networks , and the number of nodes per layer . Both encoder and decoder networks use Softplus activation functions [54], between hidden layers. The outputs of the final layers, mean Enc/Dec and diagonal covariance matrix Σ Enc/Dec for the respective parameters, are not transformed by non-linear activation functions. Note however, that for standard deviations we found it useful to predict log(Σ Enc/Dec ) instead of Σ Enc/Dec directly. Both encoder and decoder networks share the same settings for and . For a given set of input and target parameters, we examined different settings for , and , and empirically selected those parameters resulting in the most successful models. Parameter selection was guided by considerations with respect to computational efficiency. Since matrix multiplications on the GPU are carried out in matrix sections of size 4 × 4, we adapted , and so that especially the decoder model is able to efficiently utilize the parallelization capabilities on the GPU. In particular, we considered multiples of four for -values, to speed-up matrix multiplications between hidden layers, and selected so that the input size of the decoder, in + , is a multiple of four. Concerning , we considered small values, as deeper models require longer computation time due to the sequential evaluation of subsequent layers. The models were trained using an AdamW optimizer [30], implemented in the PyTorch interface in Python.
Once the training was completed, the models were re-implemented in a shader-viable form, using mainly matrix multiplications and the nonlinear activation functions.
Recall now that the likelihood of sampling a particular ray path with summary statistics ( , , , , Σ) is determined by a conditional PDF of the form ( , , , , Σ| , , ). To enable sequential model evaluation and minimize the need for modeling complex correlation structures between random variables, we decompose the PDF according to and decompose the full model into three separate CVAEs, which we will call LENGTHGEN, PATHGEN and EVENTGEN (Figure 4), one for each of the distributions , , and Σ . At inference time, the models are evaluated sequentially. The first model then infers the number of scattering events on the path. The PDF , thereby depends only on t and , because we suppose in this step that a full path is observed, which starts in the center of the sphere and terminates when crossing the sphere boundary without being absorbed in between. This approach simplifies the generation of training data, especially in scenarios with significant absorption. Conditioned on and , the probability of absorption along the bypassed path follows a Bernoulli distribution with ( , ) = 1 − . This admits simple sampling, so that training of a CVAE is not required. At this stage, a path can be rejected, if sampling yields = 1. In this case, no further models have to be evaluated, increasing computational efficiency.
If the path is accepted, the second CVAE infers the position on the sphere surface where the ray exits the volume as well as the direction into which it continues to be propagated (Figure 2a). Note here that , is independent of , which can be assumed since the absorption process is treated separately, before and are sampled, and and are conditioned on the outcome of the absorption request . If applicable, the third model finally resolves the summary statistics Σ. Since in general no dependencies can be excluded, Σ may depend on all previously sampled variables, as well as on all material parameters. An application of these statistics will be discussed in Section 4.3, where Σ = ( , ) refers to a representative sample of all scattering events along the bypassed path inside the sphere. In this setting, denotes the position of the representative scattering event and the in-going direction of the photon path at .

LEARNING MULTIPLE SCATTERING IN SPHERICAL VOLUMES
While training of the CVAEs can be performed in normalized conditions, i.e. in a unit sphere with all rays coming from the same direction, the parameters have to be adapted to the local frame of reference of the current ray at render time. Also, the radius of the sphere has to be adapted to the geometry of the scene setting. To account for this, the output parameters of the models, i.e., , , and , have to be parameterized in a way that is invariant to rotations of the spatial coordinates and changes of the sphere radius. To account for this, we propose a parameterization scheme as shown in After sampling and , the outgoing position on a sphere shell of radius sphere is parameterized in spherical coordinates through a tuple ( sphere , , ), where cos( ) = cos( in , ). Due to invariance of the distribution , with respect to rotations around in , can be sampled as a random number from a uniform distribution between 0 and 2 . This invariance is inherited from the rotation invariance of the Henyey-Greenstein phase function ( , ′ ) and can be used to restrict the learning-space of the CVAE to only the -coordinate. The rotation invariance furthermore simplifies the matching between realspace coordinates and CVAE reference system, since the orientation of the basis vectors orthogonal to in can be chosen arbitrarily, as long as the CVAE-reference direction CVAE in = (0, 0, 1) is mapped appropriately to the true in in real-space coordinates. Given a suitable rotation matrix CVAE , the direction can be sampled as where˜= (sin( ), 0, cos( )), and (0,0,1) denotes the rotation matrix corresponding to the random rotation of angle around the CVAE-reference direction CVAE in = (0, 0, 1). Given , a local frame of reference can be obtained in by considering the coordinate system n =ˆ,ˆb = in ×ˆ,ˆt =ˆb ×ˆn, withˆ= / sphere , andˆn,ˆb andˆt referring to normal, binormal and tangent vectors of the frame of reference around . The operator × denotes the cross product in R 3 . Within this system, the out-going direction can be parameterized through a tuple ( , ), wherein = cos( ,ˆb) and = cos( ,ˆt). Again exploiting rotation invariance, can be sampled in the CVAE-reference system and can be transformed to real-space by applying a transformation similar to that in Equation 3.
The treatment of summary statistics Σ has to be considered separately and depending on which quantities are required. For the case of Σ = ( , ), with and describing position and direction of a representative scattering event, a symmetry-guided restriction of sample space is not admissible since the statistics may depend on and in a complex fashion. Therefore, and are sampled as real-valued 3D vectors in the CVAE system of reference with being normalized subsequently. Both quantities are then transformed to real-space coordinates by applying a suitable rotation matrix and scaling.
To account for varying radii sphere in the model setting, the extinction coefficient of the medium is rescaled according to t → t sphere before applying the models. This is justified because of homogeneity of the medium and independence of scattering events therein. All other variables are unaffected by the transformation between model-and real-space coordinates.

Model details and shader implementation
Building upon the above considerations, the statistical model described above can be decomposed into a small number of shader functions, which accept input parameters according to the probabilities defined in Equation (2) and sample output parameters as discussed in Section 4.1. Three of the shader functions involve CVAE-based sub-steps and can be summarized follows: • Sampling of the number of scattering events: LENGTHGEN( , , ( L ) , (1) ) → .
Therein, ( ) ∼ N (0, I ) denotes the Gaussian random variables used for sampling the CVAE priors, respectively, and ( o ) ∼ N (0, I out ) denote Gaussian random variables used for sampling the decoder posterior, as discussed in Section 3.2. The settings of hyper parameters , and used for setting up the CVAEs are summarized in Table 1.

In-scattered direct illumination
In the following, the inclusion of next-event estimation for in-scattering from direct illumination [5,36] is described as a concrete application of the third model, EVENTGEN. In standard path tracing, this information is either neglected during the internal scattering process and considered only at the position where the path exits the volume, or the in-scattering is considered at every internal sampling point, e.g., by using pre-computed photon maps [37] or direct illumination using next-event estimation approaches [5,13,23]. Since at render time, the samples at which internal scattering occurs are not available any more, we set in this case Σ = ( , ), wherein indicates the position of the scattering event in the interior of the sphere and the direction of the ray arriving to . At training time, ( , ) is sampled from the set of all scattering events on the path through the sphere with weights according to the contribution of the in-scattering events to the luminosity of the ray.
For example, the direct light contribution from a directional light source at position L with power Φ along a path inside a translucent medium can be considered at every internal sampling point of the ray path. Characterizing the ray through the occurring scattering events, ( 1 , 1 ), . . . , ( , ), the total in-scattering due to direct illumination along the path is obtained as Here, L → indicates the direction of the incoming direct illumination at position and ( L → ) denotes the extinction of the light when traveling from the light source through the volume to the scattering position. Since in path tracing typically the expectation value ⟨ ⟩ over many paths is considered, our approach is to approximate in Monte Carlo fashion by choosing one single representative sample ( , ) per path and considering the in-scattering from the direct light source only for this scattering event. A sample ( , ), thereby, is chosen to be the -th scattering event in a path {( , )} =1 with probability /Λ ( ) , where Λ ( ) = =0 . This approximation does not alter ⟨ ⟩ if a sufficient number of rays is considered in the average. The crucial step in applying this approximation lies in efficiently computing the incoming direction of the direct light, L → , as well as the attenuation coefficient ( L → ) for all positions . In absence of refractive boundaries, the expressions simplify to with ∥·∥ indicating the standard Euclidean vector norm and ( L → ) referring to the distance along the direct path between L and , during which the light travels through the medium. In the more general case of refractive medium boundaries being present, however, alternative computation methods have to be employed, which are able to account for refraction appropriately. A number of methods have been proposed to tackle the problem, including [18,23,49]. However, the choice of a specific method does not affect applicability and performance gain of the path-tracing acceleration method proposed here.

RESULTS
To demonstrate the validity of our approach, we conducted numerical experiments regarding fitting accuracy of the model predictions with respect to the real distributions, visual accuracy of the rendering results including and excluding the in-scattered direct illumination, as well as performance analysis of the method in several scenarios. Note here that for all scenarios the same model configuration was used, i.e., no models were retrained to specialize for particular tasks. Training of the models was conducted on 1.6 × 10 6 ray samples with material parameters sampled from the ranges t ∈ [0, 200] and ∈ [−1, 1] to guarantee an appropriate coverage of material parameters relevant for our study. For all our experiments we assume model height to be normalized to 1m. Figure 6 shows a direct comparison between renderings computed with our proposed sphere-tracing approach and with a brute force path tracer. In-scattered direct illumination is not yet included at this stage, so only LENGTHGEN and PATHGEN are employed for rendering. To illustrate the capabilities and limitations of our method, we compare the same geometric model in two density configurations (high density and low density) and under different lighting conditions (light coming from a single light source in front of and behind the object). Since both approaches perform ray tracing and, thus, compute the same quantities, a fair comparison is achieved by fixing the number of rays traced per pixel (5000 in the current setting) and comparing the results and rendering times. Figure 6 reveals hardly any perceptual differences between our approach and brute force path tracing. Quantitatively, our renderings achieve a root mean square error to the path-traced ground truth between 0.01 and 0.02 and, thus, come very close to the ground truth. Especially the low-density rendering under backside lighting demonstrates that our approach is well able to generate translucency effects in thin material parts, which closely resemble the effects observed in the brute-force rendering.

Performance analysis
In the low-density setting in Figure 6, our approach achieves a speedup of roughly 1.5x (73 s for sphere tracing vs. 114 s for path tracing) compared to brute force path tracing. The speed-up increases with material density, so that in the high-density case our approach is already 3.1x faster than path tracing (123 s for sphere tracing vs. 381 s for path tracing). This trend continues when the optical density of the material is further increased. As seen in the teaser (Figure 1), speed-ups of roughly 7.7x or more can be achieved for sufficiently dense materials.
Given the increase in relative performance gain of our method with respect to standard path tracing, we further shed light on the different factors causing this behaviour. Therefore, we rendered different geometric objects with varying parameter settings and examined the number of scattering events per path as well as the rendering performance depending on object geometry and optical material density. Figure 7 visually encodes the number of scattering events occurring per path for different geometries. The corresponding rendering times of our approach and brute force path tracing are shown in Figure 8.
The Stanford bunny in the top row of Figure 7 is characterized by a compact shape, with many convex parts. As a result, light can penetrate deeply into the body, yielding long light paths with a large number of scattering events. As to be expected, the number of scattering events per path increases strongly with increasing optical density in brute force path tracing (left). Our method, in contrast, can take advantage of the voluminous convex shape to construct large spheres, which efficiently propagate the rays through the medium. Even though also in our approach the number of scattering events increases with increasing optical density, the rate is much lower. While in the case of the highest density the brute force path tracer is required to simulate up to 3000 scattering events per path, our method propagates the rays in less than 50 steps. The need for sequential path evaluations, which cannot be parallelized on the GPU, is thus reduced by a factor of up to 60x, and drastically reduces the amount of time that is spent waiting for long sequences of scattering events to converge. The rendering times in Figure 8 confirm these findings. Notably, compared to our approach the time increases at a higher rate when using brute force path tracing. In the case of highest optical density, a rendering time for a per-pixel path-trace of about 30 ms stands against a rendering time of more than 500 ms, yielding a speed-up by a factor of more than 16x. The difference between the reduction in scattering events, 60x, and the final speed-up can be attributed to the increasing cost of evaluating the neural network models, compared to the few operations required in a standard path tracing.
When rendering the statue model in the bottom row of Figure 7, however, some limitations of our approach become paramount. Since the statue possesses many thin features, which do not form a voluminous material body, the sphere tracing approach is unable to build large spheres and propagate rays efficiently. Especially in the case of low-densities, our approach requires only slightly less model evaluations than the brute force method. In combination with comparatively costly model evaluations in our approach, the performance gain becomes very low, yielding a speed-up of less than 2x (see Figure 8). Nevertheless, for larger material densities the number of scattering events per path increases much stronger in the brute force method   than in our approach, so that again a significant performance gain can be achieved, i.e. roughly 5x at the highest density.
From these considerations, we conclude that our method can be employed at high efficiency if standard path tracing leads to exceedingly long scattering paths. Thus, another factor limiting efficiency of our approach is absorption, which has not been considered strongly in previous experiments. Figure 9 explores this effect in detail. For this experiment, we rendered a geometric object with constant material parameters and under fixed lighting conditions, while varying the absorption rate of the medium (1 − ). From the figure it is obvious that higher absorption leads to shorter ray paths and thus shortens the rendering time. At higher absorption, the advantage of our approach decreases.  path tracing to generate a distribution of ground-truth rays, starting in the center of a unit sphere in direction in = (0, 0, 1) and evaluate the ray statistics ( , , , , , ) for each of the rays. The resulting distribution is then compared to the distribution, which is obtained from sampling the respective CVAE models. For all experiments, we draw 10000 samples per parameter setting ( t , ), for both ground truth and CVAE statitics.. The distribution of samples drawn from the LENGTHGEN model is compared to the respective ground truth in Figure 10. Ground truth statistics are summarized in a log-scale histogram, shown in blue, and overlaid by an orange curve reflecting the corresponding statistics obtained from sampling LENGTHGEN. It can be seen that for small values of t (top row), a large fraction of rays leaves the sphere already after the first scattering event. This results in a bipartite structure of the histogram, consisting of a pronounced peak of the histogram counts around = 1 and a smooth tail indicating rays that are scattered more than once. LENGTHGEN learns to compensate for this and is able to reproduce the peak accurately, especially in the case of strongly anisotropic forward scattering, = 0.9. Only in the case of backward scattering, = −0.7, the ground truth is underestimated significantly. Nevertheless, the distribution of scatter counts in the case of > 1 is reproduced to reasonable accuracy in all cases. For larger values, the significance of the peak at = 1 decreases and the ground-truth statistics appear to converge to a unimodal shape, similar to a log-Gaussian. For larger values of t , notably t = 20, 100, the matching of the distributions is good, and also the dependence on is captured appropriately. The highest inaccuracy is observed at intermediate scattering coefficients, t = 4, where the multi-scattering contributions are dominant, but single scattering still has a notable contribution. There, it seems difficult for the model to predict the balance between both effects, especially as a function of . These observations can be used, however, to bias the statistics of t and . By increasing the number of samples in these difficult regions a more accurate fit can be achieved.

Quality of the learned distributions
For the PATHGEN model, three parameters need to be validated against ground truth. The sampling statistics of cos( ) are shown in Figure 11a. Again, the blue histogram refers to ground truth statistics, whereas the orange line reflects the PATHGEN predictions. Notably, strongly anisotropic scattering is observed only in cases where t is small or is close to extreme values −1 or 1. There is an overall good fit between learned and ground-truth distributions, only for strongly anisotropic cases, = 0.9 and = −0.7, the model has minor difficulties with correctly handling the importance of single-scattering rays.
To validate the distribution of and , we generate a contour map, which shows the distribution function of and as a function of cos( ). Two-dimensional histograms of the distribution of (cos( ), ) and (cos( ), ) are obtained and normalized by the number of observations. The resulting bin probability values are finally    smoothed with a Gaussian filter to reduce visual clutter. The groundtruth distribution is shown as blue-filled contours, the learned distribution is overlaied as orange isocontours at probability iso-values matching those of the ground truth. The comparison is shown in Figures 11b and 11c. Similar to the distribution of cos( ), the greatest anisotropy is observed for highly anisotropic cases, i.e., = 0.9 or = −0.7, and small t . Nevertheless, the learned distributions fit the ground truth to reasonable accuracy.
We conclude that the major challenges for all models lie in learning a reasonable PDF at small values of t and in correctly resolving the dependence of the PDF on . From a practical perspective, learning results can be improved, however, by specializing the model at training time towards those parameter regions, which are particularly relevant during upcoming rendering processes. By over-weighting ray samples with particular ( t , )-settings at training time, these regions of parameter spaces are assigned a larger weight in the loss function and thus affect the outcome of the training process more strongly. If, for example, only media with largely forwarddirected scattering are considered during rendering, there is no need to include negative values of at training time. For our purposes, however, we did not make use of such optimizations.

In-scattered direct illumination
To also illustrate the use of the third CVAE model, EVENTGEN, we compare the result of rendering a geometric object with the full sequence of three models, LENGTHGEN, PATHGEN and EVENT-GEN. For comparison purposes, we render the same image also with a standard path tracer, as well as with a brute-force path tracer with next-event estimation using direct illumination. To obtain a fair comparison, the image was first rendered with the standard path tracer, using 1000 ray samples per pixels and rendering time was measured to be 129 s. The remaining renderers finally were allowed to render the same geometry for the same amount of time. In the current implementation, both the extended brute-force path tracer as well as the sphere tracing model draw a shadow ray to a single light source in the scene setting to account for direct illumination. As such, they do not handle refractive object boundaries accurately and yield biased results in comparison with the standard path tracer. Nevertheless, Figure 12 illustrates nicely that our methods yields a lower variance than both other algorithms. This is a promising result and suggests that CVAEs can be employed successfully for learning summary statistics of the bypassed path sections. Also, the performance gain achieved is orthogonal to the problem of accurately treating refractive boundaries, since potential fixes, enabling accurate inclusion of boundary effects, would affect our model and the reference path tracer in identical ways.

CONCLUSION
In this work we have introduced and analyzed a network pipeline comprised of CVAEs that learns to bypass multiple-scattering paths in anisotropic volumetric material. We have demonstrated that this pipeline can be embedded into a sphere-tracing approach to speed up the performance of Monte Carlo path tracing. For different parameter settings and geometries, we have demonstrated numerical accuracy of the inferred results. We are particular intrigued about the visual quality of the rendering results compared to ground truth Monte Carlo path tracing, as well as the performance gain that can be achieved for objects with large volumes of material.
On the other hand, we have seen that for thin objects where sphere tracing cannot exploit its full potential, the performance gains are only marginal. To minimize the overhead in computation time in low-density or thin object parts, we will consider a case-sensitive implementation of the sphere-tracing algorithm in the future. For shell tracing [33], a fall-back strategy to standard path tracing was already used to avoid sphere tracing in the limit of very small radii. Efficiently parallelizing CVAE evaluation and standard path tracing on GPU architectures concurrently, however, appears challenging. Furthermore, the proposed combination of deep learning and sphere tracing opens the possibility to efficiently render volumetric objects like clouds. In this context, we envision, in particular, the integration of direct in-scattering to result in significantly improved performance.

PT
PT-DI ST-DI 129 s 129 s 129 s Figure 12: Convergence improvement using next-event estimation. All images are rendered in 129 s. Left: Brute-force path tracing (1000 rays per pixel), excluding next event estimation. Middle: Brute-force path tracing (400 rays per pixel), including next-event estimation using in-scattered direct illumination at every internal sample. Right: Our approach (1465 rays per pixel), including nextevent estimation using one representative sample ( , ) per sphere tracing step.