A practical path guiding method for participating media

Rendering translucent materials is costly: light transport algorithms need to simulate a large number of scattering events inside the material before reaching convergence. The cost is especially high for materials with a large albedo or a small mean-free-path, where higher-order scattering effects dominate. In simple terms, the paths get lost in the medium. Path guiding has been proposed for surface rendering to make convergence faster by guiding the sampling process. In this paper, we introduce a path guiding solution for translucent materials. We learn an adaptive approximate representation of the radiance distribution in the volume and use it to sample the scattering direction, combining it with phase function sampling by resampled importance sampling. The proposed method significantly improves the performance of light transport simulation in participating media, especially for small lights and media with refractive boundaries. Our method can handle any homogeneous participating medium, with high or low scattering, with high or low absorption, and from isotropic to highly anisotropic.


Introduction
Computing illumination simulation in scenes with participating media is still a costly process, even with algorithms specifically designed to handle them.To render high-albedo materials, such as milk or skin, we need to simulate a large number of scattering events before convergence.Path tracing, as a simple and powerful method, has been used to simulate these materials in the industry recently.However, it takes a very long time to converge.It's even worse for small light sources.The main reason is that sampling the out-scatter direction only depends on the importance sampling of the phase function, and has no knowledge of the place of the light.This results in a lot of low contribution path in the end.
Path guiding has been recently introduced for surface rendering [6,23,28,31].The common goal is to find an "optimal" distribution that can approximate the actual path integral and make convergence faster.In these prior works, the incoming radiance distribution of some samples is learned and further used by combining with Bidirectional Reflection Distribution Function with multiple importance sampling [23] or product importance sampling [6].Besides path guiding methods in path space, several works have focused on the primary sample space [5,24,32].Some of these works used machine learning for path guiding.Path guiding makes convergence faster compared to the original path tracing.
Zero-variance methods [20,22] have been proposed for participating media, but only for isotropic homogeneous media.Recently, it is further improved in a concurrent work [7], which used path guiding for arbitrary homogeneous media based on zero-variance random walk theory.
In this paper, we present a different path guiding method for arbitrary homogeneous participating media.Similar to Müller et al. [23], we use a SD-Tree to represent the incoming radiance distribution in the scene.However, the distribution is sampled inside the volumes rather than on the surfaces, and also learned in a bidirectional manner, considering difficult paths due to the participating media boundaries.We propose a resampling method, combining this incoming radiance distribution with the phase function to better sample the product of two high frequency functions.Our resampling method only requires a couple of samples to approximate the outgoing radiance distribution.For more efficiency, we introduce a selective sampling strategy, choosing the sampling method based on the phase function.We use path guiding to select direction in the medium.Our method results in faster convergence than original path tracing, greatly reducing noise in the pictures without impact on accuracy.On top of our method, denoising methods can further remove the noise.Our method produces a better input for denoising compared to classical path tracing.To summarize, we make the following technical contributions: • a modified learning method for SD-Tree to represent incoming radiance distribution in the volume, • a resampled importance sampling method to sample the product of incoming radiance distribution function and phase function, • a selective sampling strategy to choose the sampling method based on the phase function anisotropy of the media.In the next Section, we review some of the previous work on rendering participating media.In Section 4, we present path guiding in participating media.We present our results, compare with previous work and analyze performances in Section 5, and conclude in Section 6.

Participating Media Rendering
Path tracing: Simulating light transport in participating media requires simulating multiple scattering events in the medium.Early approaches used Monte-Carlo integration to compute illumination [21,26,27].
In recent works, Georgiev et al. [4] used the product of geometry and scattering terms for importance sampling in rendering participating media.Compared to ours, this method is limited to low-order scattering, while we handle both high-order and low-order scattering.
Photon Mapping: Chandrasekhar [3] introduced the radiative transfer equation, describing radiation transport in participating media.
Jensen and Christensen [16] apply this equation to light transport, and presented an algorithm based on Photon Mapping.Jarosz et al. [12][13][14] extend this algorithm using beams instead of photons for faster computations, with less noise.Křivánek et al. [19] automatically selects between beams, points and paths in participating media using multiple importance sampling.Bitterli and Jarosz [2] further extend the idea by tracing photon planes and volume.
Photon-mapping-based methods can provide high quality simulations of light transport in participating media, but they are usually biased.In contrast, our method provides unbiased result.
Radiance Caching: Jarosz et al. [11] introduced radiance caching in participating media and decreased the rendering cost, by storing radiance at sample places in an octree.Marco et al. [17] improved this work by considering the second-order occlusion.Compared with these methods, our method also improves the performance, but in a different way.Furthermore, our method is unbiased, while theirs are biased.

Path Guiding
Several path guiding methods have been proposed to drive the importance sampling of outgoing ray directions based on learned information.Most of them focus on surface-tosurface interactions and one work presented a path tracing to participating media.
Path Space: Jensen [15] proposed to place regular bins to collect photons traced from light and reconstructs histograms from these bins to guide sampling directions when tracing paths from camera.Vorba et al. [31] learned the incoming radiance distribution of sampled points in the scene by shooting photons and use the trained distribution to guide direction sampling in path tracing.The learned incoming radiance distribution is represented with Gaussian Mixture Model.Müller et al. [23] use an adaptive binary tree in the spatial domain and a quad tree in the angular domain to store the incoming radiance distribution.They use multiple importance sampling of the Bidirectional Reflection Distribution Function (BRDF) and the learned distribution.Herholz et al. [6] used the product of the trained incoming light distribution and BRDF to sample the outgoing direction, resulting in higher sampling quality, at the cost of the product operator.
Simon et al. [28] proposed to use path guiding only for paths with high variance or difficult paths, and use regular path tracing for other paths, as path guiding is more expensive compared to the regular path tracing.
Deep learning has also been used in path guiding.Müller et al. [24] proposed a deep neural network model to present the probability density function of samples.The learned model was leveraged for sampling the ray outgoing direction.The advantage of this work is independence from scenes and local point parameters (e.g.including textured BRDF).However, this method is biased and the sampling is very expensive.
Primary Sample Space: Other researchers have worked on path guiding in Primary Sample Space (PSS) instead of path space.PSS is the space of random numbers that are used to  generate paths.PSS has been widely used in Metropolis Light transport algorithms [30].Guo et al. [5] proposed a path guiding method in PSS.The random numbers are sampled from an trained structure, kd tree, to generate the paths.Müller et al. [24] and Zheng and Zwicker [32] provided PSS based methods for path sampling using deep learning.
One issue of the PSS based methods is that they can not handle long path very well, as the random number count depends on the path length.We did not use them for participating media as these materials result in very long paths, due to many scattering events.

Zero-Variance Based Methods
Contrary to path space and primary sample space path guiding methods who use data-driven models or neural networks to guide sampling, zero-variance random walks rely on analytical models for sampling.
Křivánek and d'Eon [20] used zero-variance random walk in dense, isotropic and highly scattering media, by fitting a slab to the bounding geometry, which helps to guide the random walks back to the surface.Later, Meng et al. [22] improved it, considering thin geometry or back-lit cases.
Herholz et al. [7] used zero-variance path sampling theory and exploited the vMF model to represent the involved radiance distributions to guide all the sampling strategy in volumes including distance, direction, russian roulette and splitting.
Compared to their work, we use different representation and only focus on direction sampling.

Radiative Transfer Equation
We consider a scene containing objects with translucent material.Each of these is assumed to be made of an homogeneous material, with index of refraction η, scattering coefficient σ s , absorption coefficient σ a and phase function p( ω, ω t ).We note the mean-free path inside the material (mfp), with 1/ = σ t = σ s + σ a .
Radiance leaving point x in the direction ω is the sum of exitant radiance from the nearest surface along this direction and in-scattered radiance from the medium among the whole length of the ray [3] (see Figure 1): where T r is the transmittance, defined as: s is the distance through the medium to the nearest surface at x s = x − s ω, and x t = x − t ω with t ∈ (0, s).L(x s , ω) is the exit radiance from the nearest surface, which is governed by the rendering equation [18].L i (x s , ω) is the in-scattering radiance at x t from all direction ω t over the sphere of directions Ω 4π using the phase function p, defined as:

SD-Tree Path Guiding
Our proposed method is based on the practical path guiding method proposed by Müller et al. [23], thus we give a detailed review.They used a binary tree to represent the spatial distribution and quad trees to represent the angular distribution of incoming radiance, so we call it SD-Tree Path Guiding.
They learned an approximate representation of the scene's spatio-directional radiance field in an unbiased and iterative manner.To that end, they propose an adaptive spatiodirectional hybrid data structure, referred to as SD-Tree, for storing and sampling incident radiance.The SD-tree consists of an upper part -a binary tree that partitions the 3D spatial domain of the light field and a lower part-a quadtree that partitions the 2D directional domain.
They use an iterative scheme similar to the one proposed by Vorba et al. [31].They trained the incoming radiance sequence L 1 , L 2 , ..., L M , where L 1 was estimated with just BSDF sampling, and for all k > 1, L k was estimated by combining samples of L k−1 and the BSDF via multiple importance sampling.Two SD Trees are kept at the same time, while one is used to represent L k−1 and the other is used for updating the incoming radiance distribution in the k th iteration.
For a given shading point, two steps are performed to sample the outgoing direction of reflected ray: first, they descend spatially through the binary tree to find the leaf node containing the shading point position; next, they sample the outgoing direction from the quadtree contained in the spatial leaf node via hierarchical sample warping.
The SD-Tree method [23] can not be used for participating media directly, for two reasons.First, the SD-Tree is learned from the camera side, this leads to very expensive learing costs for participating media with boundaries.In figure 5, we use the unmodified SD-Tree method [23] on the Candle scene.It produces a significantly darker result than the reference, due to the difficulty in finding the path to the light source from the camera side.Second, multiple importance sampling is not efficient to sample the product of two high-frequency functions (the incoming light distribution and the phase function).Figure 3 illustrates this issue: multiple importance sampling produces a darker image than resampling.This issue also appears without participating media, although it is less visible.Figure 4 shows a comparison between multiple importance sampling and resampling with high-frequency BRDFs and illumination: resampling reduces the noise in the picture compared with multiple importance sampling.
In this paper, we extend this method to participating media and solve these issues.

Resampled Importance Sampling
Monte Carlo integration of a function f ( ω) is defined: , where M is the sample count and q is the probability density function.
Ideally, the optimal sampling method is sampling q with the same distribution as f , which results in zero variance and error.However, this is impractical.q( ω) may not be sampled.Talbot et al [29] introduced resampled importance sampling (RIS) to approximately sample q( ω) by firstly  We define some notations at first.q( ω) is called target function, and s( ω) is called source function, which can be arbitrary function which can be sampled.
And the entire RIS algorithm runs as follows (see Figure 5): 1. Sample N candidates (N ≥ 1) from s( ω i ).
2. Calculate weight for each candidate: w( 3. Choose the final sample from the candidates with probability proportional to the candidate weights w( ω i o ).
A weighting term W is added to the standard Monte Carlo estimator for unbiasedness: Combining the weighting term W, we get the RIS Monte Carlo estimator: In Figure 6, we compare importance resampling and multiple importance sampling to sample two function product.Given the analytic formula of two functions, we know the target function.By comparison, we can see RIS provides closer curve to the target function, thus RIS produces higher sampling quality than MIS.

Our Algorithm
We presents a path guiding method for participating media.We use a SD-tree structure to represent the incident light field (Section 4.1).Our method has two steps: in the training step, the distribution is learned by shooting light rays; in the rendering step, we sample the out-scattered direction (Section 4.2) using importance resampling of the trained distribution and the phase function (Section 4.3).To further improve the efficiency of our method, we propose a method to select the optimal sampling approach between multiple importance sampling and importance resampling, depending on the anisotropy of the phase function (Section 4.4).

Volume Representation for Radiance Distribution
We propose a representation for incoming radiance distribution in participating media with SD trees: a spatial tree for discrete samples in the volume and a directional tree in each leaf node for angular distribution.
For each media, the spatial tree is built starting form the axis aligned bounding box (AABB) of the scene.Nodes are divided into two child nodes when the contained sample count is larger than a threshold.The division dimension is chosen iteratively.The threshold depends on the iteration Fig. 7 Our algorithm runs in two steps: training step and rendering step.In the training step, we train the incoming radiance distribution iteratively.In each iteration, the light rays are shot from light sources, scattered in the the participating media, and stored in a spatial-directional hierarchy.In the rendering step, the camera rays are scattered in the media, while the outgoing direction sampling method are chosen by considering the anisotropy of the phase function.For high anisotropic media, the outgoing directions are sampled via the resampling product of the learned incoming radiance and the phase function.For low anisotropic media, the outgoing directions are sampled via multiple importance sampling of the learned incoming radiance and the phase function.
(see Section 4.2).The total count of STree node is bounded by the memory budget, or precisely the total STree node count (N).
When arriving the leaf node, a directional node is built starting from a sphere, and subdivided into four child nodes.Each node is subdivided when its contained energy is larger than 1% of the entire energy in the entire directional tree.

Learning and Rendering
Our method runs in two steps: a learning step and rendering step.Figure 7 shows an overview of our method.
In the learning step, we learn the incoming radiance distribution of sampled position in the volume from the light side.Participating media usually has refractive boundaries.If the path starts from the camera, it's difficult to reach the light sources, especially when the light sources are small.On the other side, it's easier in a bidirectional manner.Further more, the learned distribution from the light side is view independent, so it can be reused when only camera is modified.
We train the incoming radiance distribution iteratively, use c × 2048 photon count in the first iteration and increase photon count as a factor of 2. In each iteration, we shoot rays from the light sources by sampling the light sources.The rays are traced in the scene until getting refracted into the media and then start the random walk in the media.For each scattering event, the outgoing direction is got by importance sampling the phase function.For the STree, the leaf node is subdivided, when the leaf node contains more than c × √ 2 k .c is a constant value related to the resolution of the DTree, and set as 12000, similar to [23].k represents the iteration count.The DTree is also updated according to the rule described in Section 4.1.The maximum iteration count is set as 3-6 for all our test scenes.
In the rendering step, following the classical path tracing, we shoot rays from the camera, and trace them in the scene.When the rays arrive at the volume, for each scattering event, we first traverse the spatial tree to find the leaf node that contains the scattering event.Then we use the directional tree in the leaf node for direction sampling.The outgoing direction is computed by a combination of importance resampling and multiple importance sampling of the learned incoming direction and the phase function.We will show the details about our importance resampling method in Section 4.3.

Resampling of Incoming Radiance Distribution and Phase Function
In the rendering step, we have two known functions: the phase function p( ω i , ω o ) and the learned incoming radiance distribution function Li (x, ω o ) (we use ω o from the point view of camera) at position x.We use Li ( ω o ) for short.Our goal is to importance sample the direction ω o considering the product of the two functions p( ω i , ω o ) × Li (x, ω o ), called the target function.
We use importance resampling approach [29] to approximate the product sampling of the two functions (see Figure 6): 4. evaluate the scattering event for the chosen direction (Equation 5).
The weight of each sample is defined: where k represents the k th candidate, q Li represents the pdf of incoming radiance distribution function, q p represents the pdf of phase function and λ is the weight of the phase function (set as 0.9 in all of our test scenes).
In our implementation, with only 6 candidates, our method can provide convincing results.
After choosing the out-scattering direction ω o , we evaluate the scatter event: where t is the weight computed as:

Selective Importance Sampling
On top of the importance resampling method, we propose a selective importance sampling method, combining multiple importance sampling (MIS) and resampled importance sampling (RIS), based the the anisotropy of the phase function.We use MIS for low anisotropic media (g smaller than 0.5) and RIS for high anisotropic media (g larger than 0.5).
For MIS, we sample the phase function with probability µ and the incoming radiance distribution function with probability 1 − µ, where µ is set as 0.5 for all the test scenes.The details of RIS can be found in Section 4.3.
The main reason for the selective importance sampling is that MIS method can achieve faster convergence than RIS for isotropic media or very low anisotropic media but suffers from high variance for high anisotropic media.

Results and Discussion
We have implemented our algorithm inside the Mitsuba Renderer [10].We compared our algorithm against original path tracing and UPBP (Unified Points, Beams and Paths) [19] for quality validation.We use UPBP as the reference in most cases and Path tracing for the Bathroom Scene and Candle Scene with rough boundaries.All timings in this section are measured on a 2.20GHz Intel i7 (20 cores, 40 threads) with 32 GB of main memory.
We measure the convergence rate with Mean Square Error (MSE).

Qualitative Validation
We first compare our method with the reference solution and also classical Path Tracing with equal time in Figure 8, 9, 11, 12, 13 and 14.To valid the ability to handle high frequency effects, we also provide a comparison between our method and path tracing for Bumpysphere Scene with single scattering only in Figure 10.In all these test scenes, our method provides better quality (smaller MSE) than path tracing with equal time, thanks to our path guiding.We also compare the denoised results (using Intel Open Image Denoise [9]) with our method and equal time path tracing as the inputs.By both visual and quantitative comparison, our method provides a better input for denoising.The cost of denoising is negligible.
Among all these scenes, we provide both high anisotropic media, such Wax in Figure 9 and Bumpyshpere media in Figure 11 and also isotropic media, such as Marble in Figure 12 and Fog in Figure 14.

Performance and Timings
For all test scenes, we report the computation time of our method, Path tracing with equal time and the reference in Table 2.In our method, we also provide the cost of training step.
Comparing our method and classical path tracing with equal time, our method has less errors (MSE), which confirms the higher quality of our method.With equal time, our method has less sample count than path tracing, due to the extra cost for the resampled importance sampling in the path guiding.
Our method has an extra training cost to learn the incoming radiance distribution.The training costs depend on the scene complexity, media types and the learning iterations.For the equal time comparison, the training time is included in the total time.

Selective Sampling Method Validation
In Figure 15, we show the results of two sampling methods: multiple importance sampling (MIS) and importance resampling (RIS).For isotropic media (g = 0.0), MIS provides higher quality than RIS with equal rendering time.The cost of MIS is less expensive than RIS, so it can have more sample count with same time.For high anisotropic media (g = 0.8), RIS has higher quality              than MIS.Using high anisotropic media and a small area light source, both the incoming light distribution and the phase function have high frequency.Multiple importance sampling between these two functions might result in very low probability density function in one of the functions.This is shown in Figure 3.The same problem also exists in surface rendering, but it becomes even more obvious in participating media due to the high dimension of paths.On the other hand, resampling tries to sample the product of the two functions, though only uses 6 candidates, and still provides higher sampling quality.Figure 16 displays the impact of the media anisotropy on path tracing and our two different sampling methods: MIS and RIS.With low anisotropy (from 0.0 to about 0.5), MIS provides better quality than RIS.After increasing the anisotropy further, RIS provides better quality than MIS, thanks to the higher sampling quality.Thus, in our implementation, we use 0.5 as a threshold to choose the sampling methods.

Parameter Analysis
Figure 17 shows convergence (error) of our method and path tracing with varying sample count.Increasing the sampling count, both methods have decreasing error.Our method has less error than path tracing consistently.
Figure 18 and 19 show the impact of several parameters on our algorithm for the Candle Scene (Figure 9):  So we use 3-6 iterations for all of our test scenes.Figure 20 and 21 show the impact of surface and media parameters on our algorithm for the Candle Scene (Figure 9): • Figure 20 shows the impact of surface roughness on our algorithm.We compare our method and path tracing with varying surface roughness (0.01, 0.1 and 0.5).With low roughness (0.01 or 0.1), our method has higher quality than path tracing with equal time.
With high roughness (0.5), path tracing has higher quality than ours.When the surface is rough, the rays for the camera have higher chance to reach the light source, thus our method loses the benefit of guiding the camera rays to the light sources.With the extra cost for guiding, our method results in slower convergence for roughness surfaces.• Figure 21 displays the impact of media anisotropy and mean free path on our algorithm.From isotropic media   to high anisotropic media or from small mean free path to large mean free path, our method produces better quality results than the equal-time path tracing consistently.The rendered images can be found in the supplemental material.

Discussion and limitation
In Figure 2, we compare our method with Müller et al. [23] on the Candle Scene.Müller et al. produces darker result, as their method has difficulties to find paths to light from the camera side due to the refractive boundaries and the multiple importance sampling used in their method leads the path to exit the participating media too early.Thus, we believe Müller et al. [23] can not be used for participating media.
Our method shows advantage for media with refractive boundaries.As shown in Figure 20, when the surface roughness increases, our method loses its benefits due to the expensive sampling method.An automatic switching between the proposed sampling methods and the simple phase sampling will be helpful.
We use path tracing as our comparison method, as our method is based on path tracing.We would like to keep the comparison fair.
Our method does not rely on denoising method.We provided the denoised results to show that our method can produce better inputs than classical path tracing.
We only consider the out scattering direction sampling in our paper, but distance sampling is also very important for media simulation.Including distance sampling will further improve the convergence speed.

Conclusion
We have presented a practical path guiding method for participating media.We use a SD-tree to represent the incoming radiance distribution in the volume.We propose a novel product importance sampling method, combining the incoming radiance distribution and the phase function, with only a couple of samples.The proposed method provides faster convergence compared to classical path tracing.
We also did a full analysis of the light and phase function distributions, and provided a selective importance sampling strategy, to avoid expensive product importance sampling where it does not required.In the end, we got even faster convergence.
Our algorithm is suitable for any media, from low frequency effect where multiple scattering dominant to high frequency effects where single scattering dominant.
In the future, we will consider introduce our method to heterogeneous participation media and we also want to combine the volume path guiding representation and the surface path guiding representation.As for nonexponential correlated media [1], changing exponential distance sampling to non-exponential one will stimulate deeper thoughts in volumetric path guiding domain.

Fig. 1
Fig. 1 Light transport in participating media: symbols used in equation 1.

Fig. 2 Fig. 3
Fig.2Comparison between SD-Tree path guiding[23] and our method, with equal time, on the Candle scene.The SD-Tree method has difficulty finding paths connecting to the light source and produces a darker image.

Fig. 4 Fig. 5
Fig.4Comparison between multiple importance sampling and resampling on the Kitchen scene, with highly glossy surfaces.Resampling reduces the noise even without participating media.

Fig. 6
Fig. 6 Comparison between multiple importance sampling and RIS of two function f 1 and f 2 product.The RIS produces result much closer to the target function (product of two functions) than MIS.

1 ., ω k o 2 .
use weighted addition of phase function and incoming radiance distribution function as the source function, sample of the source function and get n outgoing direction candidatesevaluate the weight of each candidate (Equation4),3.compute the Cumulative Distribution Function (CDF)of the weights of the candidate out-scattering direction and use it to importance sample the final out-scattered direction.

Fig. 8
Fig. 8 Comparison between our method and path tracing with equal time on the Oil Scene.

Fig. 9 4 Fig. 10
Fig. 9 Comparison between our method and path tracing with equal time on the Candle Scene.

Fig. 11
Fig. 11 Comparison between our method and path tracing with equal time on the Bumpy Sphere Scene (Full Solution).

Fig. 12
Fig. 12 Comparison between our method and path tracing with equal time on the Athena Scene.

4 Fig. 13 Fig. 14
Fig. 13 Comparison between our method and path tracing with equal time on the Stilllife Scene.

Fig. 15
Fig.15Comparison between our method with multiple importance sampling (MIS) and our method with importance resampling (RIS) on the Candle Scene.Material: g = 0.0 and g = 0.8.

Fig. 16
Fig.16Error (MSE) for our method with multiple importance sample, our method with resampled importance sampling and path tracing with equal time (about 20 minutes), as a function of media anisotropy, for the Candle Scene.

Fig. 17
Fig. 17Error (MSE) for our method and path tracing, as a function of the sample count per pixel, for the Candle Scene.
Fig. 17Error (MSE) for our method and path tracing, as a function of the sample count per pixel, for the Candle Scene.

Fig. 18
Fig. 18 Error (MSE) as a function of resampling sample count with equal samples per pixel (spp: 1029) (left) and equal time (30 minutes) (right) for the Candle Scene.

Fig. 19
Fig. 19 Error (MSE) as a function of training iterations with different training time and same rendering time (30 minutes) for the Candle Scene.

Fig. 20
Fig. 20 Comparison between our method and path tracing with equal time on the Candle Scene with varying surface roughness.

Fig. 21
Fig. 21 Error (MSE) for our method and path tracing with equal time (about 20 minutes), as a function of media anisotropy and the media mean free path, for the Candle Scene.
Parameters for the materials used in this paper.
Our MethodPath Tracing (Equal Time) Reference scene #iter.time (train) memory spp (render) time (render) time (total) MSE #sample time Computation time, memory costs and error for our test scenes.#iter.represents the training iterations.