P_N-Method for Multiple Scattering in Participating Media

by   David Koerner, et al.

Rendering highly scattering participating media using brute force path tracing is a challenge. The diffusion approximation reduces the problem to solving a simple linear partial differential equation. Flux-limited diffusion introduces non-linearities to improve the accuracy of the solution, especially in low optical depth media, but introduces several ad-hoc assumptions. Both methods are based on a spherical harmonics expansion of the radiance field that is truncated after the first order. In this paper, we investigate the open question of whether going to higher spherical harmonic orders provides a viable improvement to these two approaches. Increasing the order introduces a set of complex coupled partial differential equations (the P_N-equations), whose growing number make them difficult to work with at higher orders. We thus use a computer algebra framework for representing and manipulating the underlying mathematical equations, and use it to derive the real-valued P_N-equations for arbitrary orders. We further present a staggered-grid P_N-solver and generate its stencil code directly from the expression tree of the P_N-equations. Finally, we discuss how our method compares to prior work for various standard problems.



page 1

page 2

page 7


Fast global spectral methods for three-dimensional partial differential equations

Global spectral methods offer the potential to compute solutions of part...

Sparse spectral methods for partial differential equations on spherical caps

In recent years, sparse spectral methods for solving partial differentia...

Partial differential equation solver based on optimization methods

The numerical solution methods for partial differential equation (PDE) s...

Symbolic Solutions of Simultaneous First-order PDEs in One Unknown

We propose and implement an algorithm for solving an overdetermined syst...

A Scalable Framework for Solving Fractional Diffusion Equations

The study of fractional order differential operators is receiving renewe...

FDTD: solving 1+1D delay PDE

We present a proof of concept for adapting the finite-difference time-do...

On time-domain NRBC for Maxwell's equations and its application in accurate simulation of electromagnetic invisibility cloaks

In this paper, we present analytic formulas of the temporal convolution ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Simulating light transport in participating media remains a challenging problem for image synthesis in computer graphics. Due to their ability to produce unbiased results and conceptual simplicity, Monte Carlo based techniques have become the standard approach [NGHJ18]. The main downside of these methods are their computational demands when rendering media with strong scattering or anisotropy.

Deterministic methods have enjoyed less popularity, because they suffer from discretization artifacts, produce biased results, cannot be coupled easily with surface rendering problems and are trickier to implement. However, their appeal lies in the fact that they produce a global solution across the whole domain and have better performance for certain problems [Bru02].

The work on path-guiding techniques from recent years [MGN17] has shown how approximate representations of the steady-state transport in a scene can be used to accelerate Monte Carlo integration techniques, such as path tracing. Instead of generating these approximate representations using Monte Carlo methods, deterministic methods may offer a viable alternative. Hybrid methods could combine the performance benefits of deterministic methods with accurate and unbiased Monte Carlo results. Deterministic methods also lend themselves to applications where fast approximate solutions are preferable over correct, but slowly converging results.

For these reasons, we suggest it is important for volume-rendering researchers to study deterministic methods and have a solid understanding of their characteristics and performance traits for typical rendering problems.

The -method is a deterministic method of solving the radiative transfer equation (RTE) which is used in other fields such as medical imaging and nuclear sciences, but has not found use in computer graphics thus far. The purpose and main contribution of our paper is to gain a solid understanding of its foundations and present a method for using it in the context of rendering. In particular, we present these theoretical and practical contributions:

  • We derive and present the time-independent real-valued -equations and write them down in a very concise and compact form which we have not found anywhere else in the literature.

  • We introduce a staggered-grid solver, for which we generate stencil code automatically from a computer algebra representation of the -equations. This allows us to deal with the increasingly complex equations which the -method produces for higher order. It further allows our solver to be used for any (potentially coupled) partial differential equations, which result in a system of linear equations after discretization.

  • Finally, we compare the -method for higher orders against flux-limited diffusion and ground truth Monte Carlo integration.

In the next section, we will discuss related work and its relation to our contribution. In Section 3 we revisit the deterministic approach to light transport simulation in participating media and outline the discretization using spherical harmonics. In Section 4 we introduce our computer algebra representation, which we required to derive the real-valued -equations, presented in Section 5. This representation is also a key component of our solver, which we present in Section 6. Section 7 discusses application of the solution in the context of rendering. We compare our -solver against flux-limited diffusion for a set of standard problems in Section 8. Finally, Section 9 concludes with a summary and review of future work.

2 Previous work

Light transport in participating media is governed by the RTE, first studied in the context of astrophysics by Chandrasekhar [Cha60] and later introduced to computer graphics by Kajiya [Kaj86]. In computer graphics today, this equation is typically solved using Monte Carlo methods [NGHJ18]. However in strongly scattering or highly anisotropic media these methods can become prohibitively expensive, for example in the case of a high albedo medium such as milk where tracing paths with a huge number of scattering events is necessary.

In contrast to path-tracing, the -method [Bru02] gives a solution by solving a system of linear equations. It is derived by discretizing the angular variable of the radiative transfer equation into spherical harmonics (SH). This gives rise to a set of coupled, complex-valued partial differential equations, called the -equations. The subscript refers to the spherical harmonics truncation order.

The -method has a long history in other fields and was never applied in graphics. Kajiya [KVH84] explained the theory, but did not give any details on implementation or how to solve it. In fact, as Max [Max95] pointed out, it is not clear if Kajiya succeeded at all at applying the method, as all of the results in his paper were produced with a simpler method. This is further strengthened by the fact that a straightforward finite difference discretization of the -equations produces unusable results, due to oscillation artifacts in the solution  [SF14]. We use a staggered-grid solver, motivated by the solver of Seibold et al. [SF14], that produces artifact-free solutions (see Figure 1).

Figure 1: Solving the 2D checkerboard problem using naive collocated grids produces oscillating artifacts (left). Our solver uses staggered grids and produces artifact free results (right).

Related to the -method, the classical diffusion approximation (CDA) is a deterministic method which arrives at a solution by solving a system of linear equations. It corresponds to the -equations when (truncation after the first SH order), which can be collapsed to a simple diffusion equation, giving the method its name. CDA has a long history in other domains, such as astrophysics and nuclear physics [Bru02], and was introduced to graphics by Stam [Sta95].

CDA suffers from severe energy loss close to regions with strong density gradients [KPS14]. The problem can be addressed by a modification known as the Variable Eddington factor (VEF) method [Bru02], which non-linearly adjusts the diffusion coefficient to improve the solution near density gradients and low-density regions. Flux-limited diffusion, developed in the context of astrophysics by Levermore et al. [LP81] and later introduced to graphics by Koerner et al. [KPS14]

, is the most prominent example. VEF is based on modifying the closure of the moment hierarchy in a way which more accurately than CDA models the transition from the diffusive regime (in which photons effectively undergo a random walk) to the transport regime (in which photons travel along straight lines) as the medium opacity decreases. VEF methods produce better results than CDA for isotropic media, but they do not provide a satisfactory treatment of anisotropic media.

Although the -method should provide an increasingly accurate solution to the full RTE including anisotropy as the truncation moment order is increased, at the expense of more computation time, it is still an open and unresolved question how high a truncation order is needed for a particular problem in order to produce significantly better results than first order non-linear diffusion methods (FLD and VEF). This question has also been raised in other domains [OAH00].

3 Discretized Radiative Transfer Equation

Our method derives from the radiative transfer equation (RTE), which expresses the change of the radiance field , with respect to an infinitesimal change of position in direction at point :

The left hand side (LHS) is the transport term, and we refer to the terms on the right hand side (RHS) as collision, scattering, and source term, respectively. The symbols , , , and refer to the extinction- and scattering coefficient, phase function and emission.

The RTE is often given in operator notation, where transport, collision, and scattering are expressed as operators , and , which are applied to the radiance field :


Deterministic methods are derived by discretizing the angular and spatial domain. This gives rise to a linear system of equations, which can be solved using standard methods. For the -method, the angular variable is first discretized, using a truncated spherical harmonics expansion. This results in the -equations, a system of coupled PDEs that still depend on a continuous spatial variable.

The number of equations grows with the truncation order . This is why the discretization is laborious and difficult to do without errors if done by hand. We therefore use a computer algebra representation to automate this process. After giving an outline of the general discretization in this section, we will present our computer algebra representation in the next section. The -equations that result from our automated discretization are given in Section 5.

Since the radiance field is real, we use the real-valued SH basis functions , which are defined in terms of the complex-valued SH basis functions as follows [JCJ09]:


We express the projection into the spherical harmonics basis functions with a projection operator :

The -equations are derived by first expressing all direction-dependent parameters in spherical harmonics. The radiance field in Equation 3 is therefore replaced by its SH reconstruction , introducing an error due to truncation at order :

After substitution, all angular parameters are expressed in terms of spherical harmonics, but they still depend on the continuous angular variable . As a next step, we project each term of the RTE into spherical harmonics, using the projection operator . This produces a single equation for each -pair. The -equations therefore can be written as:


Once the -equations have been found, the spatial variable is discretized using a finite difference (FD) voxel grid (using central differences for differential operators).

Following this discretization, the radiance field

, is represented as a set of SH coefficients per voxel. Flattening these over all voxels into a single vector, gives the solution vector

. The RHS vector

is produced similarly. The projected operators can be expressed as linear transformations, which can be collapsed into a single coefficient matrix

(see Figure 2):


, , are matrices, which result from the discretized transport, collision and scattering operators respectively.

Figure 2: Structure of coefficient matrix and solution vector after discretization of the -equations on a finite difference grid.

4 Computer Algebra Representation

So far, we have only given the -equations in high-level operator notation (Equation 3). Carrying out the full derivation creates large, unwieldy equations and requires a string of expansions and applications of identities. These are challenging to manipulate if done by hand. We therefore used a computer algebra representation, which allowed us to derive and discretize the -equation in a semi-automatic fashion (Figure 3)).

It represents the equations using a tree of mathematical expressions, which represent numbers, symbols and other expression types, such as integrals, derivatives, sums, products and functions. Further, manipulators can be executed on these expression trees to perform substitution, constant folding, reordering of nested integrals, application of identities and more complex operations. Finally, frontends allow rendering the expression tree into different forms, such as LaTeX and C++ source code. While we ended up implementing our own lightweight framework, off-the-shelf packages, such as SymPy (www.sympy.org), exist and would be equally suitable for our use case.

Figure 3: A computer algebra framework allows us to represent equations as mathematical expressions trees. It further provides a set of functions for manipulating the tree according to valid mathematical operations, such as the binomial expansion above. Frontends allow generation of source code from the expression tree.

Using the computer algebra representation, we perform the derivation steps required to arrive at the real-valued -equations (the derivation steps shown in the Appendix were almost all rendered from the expression tree). More importantly, we use the representation to perform the discretization and generate the stencil code used by our solver. This is detailed in section 6.1.

5 Real-valued -Equations

With the help of a computer algebra representation framework, we are able to easily derive and work with the large and unwieldy -equations. We present here the final real-valued -equations, for a general and in three dimensions, in a very compact form which we have not found elsewhere in the literature. The derivation of the real-valued equations is rather long and takes many pages to describe in detail. We therefore give only the final result in this section, and present the full derivation for reference in Appendix A.

Since the real-valued SH bases (Equation 2) have different definitions for , or , we get different projections for , depending on the sign of .

For we have


For (upper sign) and (lower sign) we have




In the next section we will present the solver that we use to solve the -equations.

Figure 4: Overview of our -solver. After generating the stencil source code from the expression trees representing the -equations, the linear system is built using RTE parameter fields and additional user input, such as grid resolution and type of boundary conditions. The resulting system is solved for , which is then used in our rendering application.

6 -Solver

The truncation order is the key input parameter to the solver. With higher values, manual implementation of the solver from the equations would be arduous, error-prone and time-consuming. We therefore decided to make use of the computer algebra representation and designed our solver around it.

The solver consists of two components. The first is a precomputation (Section 6.1), which is executed once for every single value of . This step runs a partial evaluation on the mathematical expression tree and applies the spatial discretization in a reference space we call stencil space. The precomputation step automatically generates source code from the resulting expression tree.

The generated stencil code is compiled with the runtime component (Section 6.2) of our solver. This component receives the actual problem as input, including grid resolution and RTE parameter fields. It then builds the linear system and solves it using standard methods. An overview of the solver is given in Figure 4.

6.1 Precomputation

The result of the precomputation is a stencil, which can be used during runtime to build the linear system for a given problem. The stencil is a pattern of indices into the solution vector, along with values. It expresses how the sum of the weighted solution vector components relate to the RHS for a given unknown in the system, and therefore contains most information required to fill the system matrix and RHS-vector row-by-row. Note that while the coefficients may change, the sparsity pattern of the stencil is identical for different rows.

Stencil generation entails discretizing a PDE at a hypothetical center voxel (assumed to mean the voxel center most of the time). Finite differences create weighted references to other voxels (e.g. ). After bringing the discretized equation into canonical form (a weighted sum of unknowns), one can write the stencil by reading off the weights and offsets (Figure 5). Voxel will only be known during runtime, when the stencil is executed for a particular unknown (row). Then the offsets can be used to find the component index into the solution-vector, and weights can be evaluated for concrete world space position. We refer to the space with the hypothetical voxel at the center as stencil space.

Figure 5: Creating the stencil code requires several steps, usually done by hand. We express the given problem in a computer algebra representation and use it to fully automate the process.

The spatial discretization is done by parsing the expression tree from the root. The discrete substitute for the continuous position variable is initialized with . Differential operator nodes are replaced by a subtree, which expresses the finite difference approximation (including voxel-size factor ). The subtree of the differential operator node is duplicated for different offsets to the discrete variable . Since this offset only applies to the subtree, a stack is maintained by the parser to manage scope. Whenever the parser encounters the continuous variable in the expression tree, its node in the tree is replaced by the discrete substitute, currently on top of the stack. Nested differential operators yield higher order finite difference stencils as expected.

Factorization into canonical form is done as a manipulation pass on the mathematical expression tree. The result allows us to implement the code generation in a straightforward fashion. For each term, the -offset is retrieved from the unknown. The factor-expression, which is multipled with the unknown, is extracted from the tree and used to generate source code for evaluating the factor-expression during runtime (including calls for evaluating RTE-parameter fields).

Figure 6: Staggering distributes the coefficients of the solution vector onto four disjoint grids, indicated by the symbols, in a way which guarantees second order accuracy [SF14].

Our solver supports placement of coefficients at arbitrary staggered grid locations (see Figure 6). This means that during the discretization step, the discrete location

(at which the unknown is meant to be evaluated) might not coincide with the location of the unknown. To solve this, depending on how the two are located relative to each other, the parser returns an expression which interpolates the coefficients at

from their defined locations. If those happen to coincide, the coefficient itself is returned. This is also done for RTE parameters, such as or , which are always located at the voxel center. We use the staggering scheme introduced by Seibold et al. [SF14] for their solver StaRMAP, which they proved ensures second order accuracy, i.e. a truncation error of , and prevents the growth of non-physical oscillations in the solution (see their paper for the full details of the staggered coefficient placement).

6.2 Runtime

The stencil code is generated once for every value of and compiled with the runtime component of our solver. The runtime executes the stencil for every voxel to populate the system matrix and RHS vector with numerical values.

The number of rows is determined by the number of voxels times the number of coefficients per voxel (see Figure 2) and can therefore become very large for high resolution and high truncation order. The matrix is square and extremely sparse, due to the local structure of the finite differences discretization. Unfortunately, it is non-symmetric due to the transport term and not diagonal dominant, which rules out many iterative methods for solving linear systems. Iterative methods are useful, as they allow balancing accuracy against performance by tweaking the convergence threshold. We address this by solving the normal form instead. This gives a symmetric and positive definite system matrix , albeit with a higher condition number. Investigation of other solution schemes (e.g. multigrid) would be an interesting avenue for future work. However, more importantly, in the presence of vacuum regions, the matrix becomes singular and the system cannot be solved at all. This requires the introduction of a minimum threshold for the extinction coefficient .

Our solver supports both Neumann and Dirichlet boundary conditions (BC). They are handled transparently by the code which generates the stencil. Whenever the stencil accumulates coefficients into a boundary location, the framework either ignores the write operation (Dirichlet BC) or accumulates into the row and column in of the coefficient in the closest voxel inside the domain (Neumann BC). This is done by changing the index of the written component.

7 Rendering

We use an approach similar to Koerner et al. [KPS14], where we separate the radiance field into single scattered light and multiple scattered light :


The single scattered light is folded into the emission term :


This means that our solver will solve for the multiple scattered light . The quantity is the “uncollided” light, which was emitted from the light source and attenuated by the volume without undergoing any scattering event. We compute it using a few light samples per voxel, which quickly converges to a useful result for Dirac delta light sources.

Running the solver gives solution vector . We then un-stagger the solution by interpolating all coefficients to voxel centers. The additional coefficients at boundary voxels are no longer needed. This operation is represented as a matrix that produces a three-dimensional voxel grid with SH coefficients for order at the center of each voxel.

For rendering, we use a simple forward path tracing approach, where we start tracing from the camera. At the first scattering event, we use next event estimation to account for

. Then we sample a new direction according to the phase function. Instead of continuing tracing into the new direction, we evaluate the in-scattering integral using . The SH coefficients at are found by trilinear interpolation from the voxel grid of SH coefficients.

8 Results

In this section, we present results for a set of standard problems in order to validate and evaluate our method. We also compare against classical diffusion (CDA) and flux-limited diffusion (FLD).

Our computer algebra framework and the precomputation has been implemented in Python. The runtime component of our solver has been implemented in C++. We use a naive implementation of a CG solver, which has been modified such that we do not need to explicitly compute the normal form of the linear system to solve. We use the sparse matrix representation and sparse matrix vector product from the Eigen linear algebra library (eigen.tuxfamily.org).

The solver for classical diffusion is based on the diffusion equation, which is derived by collapsing the -equations:


Since our solver can work with any PDE which results in a linear system of equations, we put Equation 9 into our computer algebra representation and provide it as an input to our solver, which generates the correct stencil code automatically.

Since FLD is based on a non-linear diffusion equation, we were not able to use our system in the same way. Our implementation closely follows the implementation in [KPS14] (though ours runs on CPU) and we refer to their paper for more details.

8.1 2D checkerboard

First we ran our solver on the 2D checkerboard, a very common test case in other fields. The problem has dimensions and is discretized with resolution . Unit size blocks are filled with purely absorbing medium in a checkerboard pattern. All other blocks are filled with a purely scattering medium with .

Solving the standard checkerboard problem allows us to validate our solver against the StaRMAP solver from Seibold et al. [SF14], which solves for the time-dependent and complex-valued -equations on a staggered grid, but in the 2D case only. The 2D case is derived by assuming that all SH coefficients, RTE parameters and boundary conditions are z-independent. This causes all SH coefficients and moment equations for which

is odd to vanish. Due to the time-dependency, their approach is to do explicit incremental steps in time. We run their solver for many timesteps to get a result close to steady state.

(a) vs. ground truth
(b) vs. CDA and FLD
Figure 7: Lineplot through the 3D solution of our solver for the point source problem for various order (left). Solution for compared against classical diffusion, flux-limited diffusion and analytical solution (right).
Figure 8: Comparison of the result (here the coefficient, or fluence) for the checkerboard test using StaRMAP’s time-stepping solver [SF14] (left) against our steady-state solver (right) with . Our results are in good agreement.

As can be seen in Figure 8, the results from our solver are in good agreement with the results from Seibold et al. [SF14] and verify the correctness of our implementation. Converging to a residual of takes for and for .

8.2 Point source problem

We also run our solver for the point source problem, a single point light in a homogeneous medium. This not only helps to validate our implementation for the 3D case, but also provides information on the accuracy of these methods. We use the Grosjean approximation, which was introduced by D’Eon et al.[dI11] as a very accurate approximation to the ground truth solution.

Figure 9: Convergence behaviour of our solver with for the nebulae dataset and for varying minimum thresholds of the extinction coefficient . Threshold values and an estimate for the condition number of (MATLAB’s condest function) are shown next to the plots. The convergence deteriorates as the threshold decreases. Once it reaches zero, the presence of pure vacuum makes the condition number infinite.

For our test case, we choose a FD resolution of , an extinction coefficient and albedo . We run the solver for different truncation values . In Figure 7 fig:pointsource_pn, we see that the solution becomes increasingly accurate for higher truncation order. The ground truth is underestimated when is odd, and overestimated when is even. We further see in  7 fig:pointsource_p5 that the solution exactly matches the results from CDA, which confirms that the latter is only a collapsed version of the former. The time to solve is significant with for and with . With these performance characteristics, our -solver is clearly not competitive in comparison with our CDA and FLD solver, which are much faster.

8.3 Nebulae

Finally, we run our solver on a procedural cloud dataset to get an idea of its performance in more practical applications. Figure -Method for Multiple Scattering in Participating Mediafig:nebulae_ours shows the result of for a procedurally generated heterogeneous cloud with an isotropic phase function, with path-traced ground truth in Figure -Method for Multiple Scattering in Participating Mediafig:teaser_gt. We see that at order , our method can capture indirect illumination similarly well as FLD (-Method for Multiple Scattering in Participating Mediafig:teaser_fld) and is significantly better than CDA (-Method for Multiple Scattering in Participating Mediafig:teaser_cda) as expected. The indirectly illuminated region at the bottom appears to be closer to the path-traced result as opposed to the solution from FLD which is very flat in that region. However, in many other areas, seems to still suffer a lot from the problem of energy loss near transitions from dense regions to vacuum. It appears that going higher order a few steps mitigates this only slowly at a high cost of compute and storage. The main characteristic of the nebulae dataset is the presence of vacuum. We found that having vacuum regions in the dataset will cause the condition number to become infinite and the solver practically does not converge. We therefore introduced a minimum threshold for the extinction coefficient . Every voxel with an extinction coefficient smaller than the threshold is set to the threshold value. In Figure 9 we show the effect of the minimum threshold on the convergence behaviour. As it increases, convergence improves.

9 Conclusion

In this paper we introduced the -method to the toolbox of deterministic methods for rendering participating media in computer graphics. We derived and presented the real-valued -equations, along with a staggered grid solver with numerical stencils constructed automatically from the equations via a computer algebra system. We showed how to use the results in a rendering system and ran our solver for various standard problems for validation.

We originally set out to understand how non-linear diffusion methods compare to the -method for increasing order. Our results indicate that although the lack of higher moments makes the FLD solution overly smooth, its energy conserving nature and comparably small resource footprint make it a better approach at present for most graphics applications compared to the -method, which becomes increasingly costly for higher values of .

The literature in other fields often states that the method—when solved in normal form as we do—is able to deal with vacuum regions. We found this misleading. The -method in normal form indeed does not break down completely in the presence of vacuum as diffusion based methods do (due to in the diffusion coefficient). However, in the presence of vacuum, the condition number of the system matrix becomes infinite and does not converge either. Therefore based methods also require minimum thresholding of the extinction coefficient and offer no benefit for vacuum regions.

Much more work needs to be done in order to make the -method competitive in performance to alternative solutions for volume rendering. We believe this can be made possible by employing a multigrid scheme for solving the linear system of equations. We implemented a multigrid solver, but did not find the expected performance improvements. This is possibly due to the coupling between coefficients within a voxel, which does not work well together with the smoothing step. We want to study this in the future.

Unique to our system is that it uses a computer algebra representation of the equation to solve as input. Discretization in angular and spatial domain is done using manipulation passes on the representation. The stencil code, which is used to populate the system of linear equations, is generated from the expression tree. This way, we can easily deal with coupled PDEs and avoid the time consuming and error prone process of writing stencil code by hand.

When researching the application of the -method in other fields, we came across a rich variety of variations, such as simplified ([McC10], filtered ([RARO13], diffusion-correction ([SFL11] and least-squares ([HPM14]. These variations have been introduced to address certain problems of the standard -method, such as ringing artifacts, dealing with vacuum regions and general convergence. Our solver can be applied to any (potentially coupled) PDE and therefore can generate stencil code for all these variations by simply expressing the respective PDEs in our computer algebra representation and providing this as an input to our solver. This would allow an exhaustive comparison of all these methods and we consider this as future work.

Finally, since the approach of our solver is very generic, we also would like to explore its application to other simulation problems in computer graphics.


We thank the anonymous reviewers for their valuable and encouraging comments and feedback. We also thank Martin Frank and Benjamin Seibold for the very valuable discussions and answers to our questions. We thank Bernd Eberhardt for feedback and support. This project has been partially funded by the MSC-BW project.


  • [Bru02] Brunner T. A.: Forms of Approximate Radiation Transport. Tech. Rep. SAND2002-1778, Sandia National Laboratories (2002).
  • [Cha60] Chandrasekhar S.: Radiative Transfer. Dover Publications, 1960.
  • [dI11] d’Eon E., Irving G.: A Quantized-Diffusion Model for Rendering Translucent Materials. ACM TOG (Proc. of SIGGRAPH) 30, 4 (July 2011), 56:1–56:14.
  • [HPM14] Hansen J., Peterson J., Morel J., Ragusa J., Wang Y.: A Least-Squares Transport Equation Compatible with Voids. Journal of Computational and Theoretical Transport 43, 1-7 (2014), 374–401.
  • [JCJ09] Jarosz, W., Carr, N. A., Jensen, H. W.: Importance Sampling Spherical Harmonics. Computer Graphics Forum (Proceedings of Eurographics) (2009), 577–586.
  • [Kaj86] Kajiya J. T.: The rendering equation. Computer Graphics (Proc. of SIGGRAPH) (1986), 143–150.
  • [KPS14] Koerner D., Portsmouth J., Sadlo F., Ertl T., Eberhardt B.: Flux-limited Diffusion for Multiple Scattering in Participating media. CoRR abs/1403.8105 (2014).
  • [KVH84] Kajiya J. T., Von Herzen B. P.: Ray tracing volume densities. Computer Graphics (Proc. of SIGGRAPH) 18, 3 (Jan. 1984), 165–174.
  • [LP81] Levermore C. D., Pomraning G. C.: A flux-limited diffusion theory. Astrophysical Journal 248 (1981), 321–334.
  • [Max95] Max N.: Efficient Light Propagation for Multiple Anisotropic Volume Scattering. Springer Berlin Heidelberg, Berlin, Heidelberg, 1995, pp. 87–104.
  • [McC10] McClarren R. G.: Theoretical Aspects of the Simplified Pn Equations. Transport Theory and Statistical Physics 39, 2-4 (2010), 73–109.
  • [MGN17] Muller, T., Gross, M., Novak, J.: Practical Path Guiding for Efficient Light Transport Simulation. Computer Graphics Forum (2017), 36:91–36:100
  • [NGHJ18] Novak, J., Georgiev, I., Hanika, J., Jarosz, W.: Monte Carlo Methods for Volumetric Light Transport Simulation. Computer Graphics Forum (Proceedings of Eurographics - State of the Art Reports) (2018), 37(2).
  • [OAH00] Olson G. L., Auer L. H., Hall M. L.: Diffusion, P1, and other approximate forms of radiation transport. Journal of Quantitative Spectroscopy and Radiative Transfer 64, 6 (2000), 619 – 634.
  • [RARO13] Radice D., Abdikamalov E., Rezzolla L., Ott C. D.: A New Spherical Harmonics Scheme for Multi-Dimensional Radiation Transport I. Static Matter Configurations. Journal of Computational Physics 242 (2013), 648 – 669.
  • [SF14] Seibold B., Frank M.: StaRMAP—A second order staggered grid method for spherical harmonics moment equations of radiative transfer. ACM Trans. Math. Softw. 41, 1 (Oct. 2014), 4:1–4:28.
  • [SFL11] Schäfer M., Frank M., Levermore C. D.: Diffusive Corrections to Pn Approximations. Multiscale Modeling and Simulation 9 (2011), 1–28.
  • [Sta95] Stam J.: Multiple scattering as a diffusion process. In Proc. of Eurographics Workshop on Rendering Techniques. Springer-Verlag, 1995, pp. 41–50.

Appendix A Full derivation of the -equations

2 Derivation of the complex-valued -equations

Deriving the -equations consists of two main steps. First, all angular dependent quantities in the RTE are expressed in terms of spherical harmonics (SH) basis functions. After this, the RTE still depends on the angular variable. Therefore, the second step projects each term of the RTE by multiplying with the complex conjugate of the SH basis functions, followed by integration over solid angle to integrate out the angular variable. This gives an equation for each spherical harmonics coefficient.

Spherical Harmonics are a set of very popular and well known functions on the sphere. The complex-valued SH basis functions are given by


where are the associated Legendre polynomials. The factor is called the Condon-Shortley phase and is not part of the associated Legendre Polynomial (unlike some other definitions).

2.1 Projecting Radiative Transfer Quantities

2.1.1 Radiance Field and Emission Field

Radiative transfer quantities, which depend on position and angle , are projected into spatially dependent SH coefficients for each SH basis function:

The function is completely reconstructed by using all SH basis functions up to infinite order. The -equations introduce a truncation error by only using SH basis functions up to order for the reconstruction and :


2.1.2 Phase Function

Throughout our article, we assume an isotropic phase function, which only depends on the angle between incident and outgoing vector and (note that in the graphics literature, these would often be called anisotropic). We will see later in section 2.2.3, that this allows us to fix the outgoing vector at the pole axis and compute the phase function SH coefficients by just varying the incident vector .

The expansion of the phase function can be further simplified because the phase function is rotationally symmetric around the pole axis . Consider the definition of the spherical harmonics basis function :

Now we apply a rotation of degrees around the pole axis. In spherical harmonics, this is expressed as:

If the phase function is rotationally symmetric around the pole axis, we have:

and in spherical harmonics this would be:

By equating coefficients we get:

Since for all only when , we can conclude that for all . This means that for a function which is rotationally symmetric around the pole axis, only the coefficients will be valid. Therefore, our phase function reconstruction for a fixed outgoing vector () only requires SH coefficients with :


2.2 Projecting Terms of the RTE

2.2.1 Transport Term

The transport term of the RTE is given as

Replacing with its expansion gives:

Next we multiply with and integrate over solid angle:

We can pull the spatial derivative out of the integral to get:


We apply the following recursive relation for the spherical harmonics basis functions:



Note that the signs for the - and - component depend on the handedness of the coordinate system in which the SH basis functions are defined. Using this in Equation 13 gives

Integrating the vector term over solid angle can be expressed as separate solid angle integrals over each component. These integrals over a sum of terms are split into separate integrals. We arrive at:

Applying the orthogonality property to the solid angle integrals will will select specific in each term:

Which gives the final moment equation for the transport term:

2.2.2 Collision Term

The collision term of the RTE is given as:

We first replace the radiance field with its spherical harmonics expansion:

Multiplying with and integrating over solid angle gives, after pulling some factors out of the integral:

2.2.3 Scattering Term

The scattering term in the RTE is given as:

The phase function used in isotropic scattering medium only depends on the angle between incident and outgoing direction and therefore is rotationally symmetric around the pole-defining axis. This property allows us to define a rotation , which rotates the phase function such that the pole axis aligns with direction vector . The rotated phase function is defined as:

where is the rotation operator, which can be implemented by applying the inverse rotation to the arguments of . With this rotated phase function, we now can express the integral of the scattering operator as a convolution denoted with the symbol :


As we evaluate the inner product integral of the convolution, the phase function rotates along with the argument .

We now use the spherical harmonics expansions of (Equation 2.1.1) and (Equation 23) in the definition for the inner product of our convolution (Equation 15):

Due to linearity of the inner product operator, we can pull out the non-angular dependent parts of the expansions:

and further:

The rotation of a function with frequency gives a function of frequency again. In addition the spherical harmonics basis functions are orthogonal. We therefore have:

which further simplifies our inner product integral to:

What remains to be resolved is the inner product. We use the fact that the spherical harmonics basis functions

are eigenfunctions of the inner product integral operator in the equation above: