P N -Method for Multiple Scattering in Participating Media
EEurographics Symposium on Rendering - Experimental Ideas & Implementations (2018)T. Hachisuka and W. Jakob (Editors) P N -Method for Multiple Scattering in Participating Media David Koerner Jamie Portsmouth Wenzel Jakob University of Stuttgart Solid Angle École polytechnique fédérale de Lausanne (EPFL) (a)
Classical diffusion (b) P (ours) (c) Flux-limited diffusion (d)
Brute force path tracing
Figure 1:
Flux-limited diffusion (c) is an extension to the classical diffusion approximation (a). It improves accuracy, but is based on ad-hocassumptions about volumetric transport. We add the P N -method (b), which allows solution of the RTE with configurable accuracy, to thetoolbox of methods in volume rendering and investigate its benefits and trade-offs against flux-limited diffusion. Abstract
Rendering highly scattering participating media using brute force path tracing is a challenge. The diffusion approximationreduces the problem to solving a simple linear partial differential equation. Flux-limited diffusion introduces non-linearities toimprove the accuracy of the solution, especially in low optical depth media, but introduces several ad-hoc assumptions. Bothmethods 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 thesetwo 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 forrepresenting and manipulating the underlying mathematical equations, and use it to derive the real-valued P N -equations forarbitrary orders. We further present a staggered-grid P N -solver and generate its stencil code directly from the expression treeof the P N -equations. Finally, we discuss how our method compares to prior work for various standard problems.
1. Introduction
Simulating light transport in participating media remains a chal-lenging problem for image synthesis in computer graphics. Due totheir ability to produce unbiased results and conceptual simplic-ity, Monte Carlo based techniques have become the standard ap-proach [NGHJ18]. The main downside of these methods are theircomputational demands when rendering media with strong scatter-ing or anisotropy.Deterministic methods have enjoyed less popularity, becausethey suffer from discretization artifacts, produce biased results,cannot be coupled easily with surface rendering problems and aretrickier to implement. However, their appeal lies in the fact that they produce a global solution across the whole domain and havebetter performance for certain problems [Bru02].The work on path-guiding techniques from recentyears [MGN17] has shown how approximate representationsof the steady-state transport in a scene can be used to accelerateMonte Carlo integration techniques, such as path tracing. Insteadof generating these approximate representations using Monte Carlomethods, deterministic methods may offer a viable alternative.Hybrid methods could combine the performance benefits ofdeterministic methods with accurate and unbiased Monte Carloresults. Deterministic methods also lend themselves to applicationswhere fast approximate solutions are preferable over correct, butslowly converging results. c (cid:13) (cid:13) a r X i v : . [ c s . G R ] J u l . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media For these reasons, we suggest it is important for volume-rendering researchers to study deterministic methods and have asolid understanding of their characteristics and performance traitsfor typical rendering problems.The P N -method is a deterministic method of solving the radiativetransfer equation (RTE) which is used in other fields such as medi-cal imaging and nuclear sciences, but has not found use in computergraphics thus far. The purpose and main contribution of our paperis to gain a solid understanding of its foundations and present amethod for using it in the context of rendering. In particular, wepresent these theoretical and practical contributions: • We derive and present the time-independent real-valued P N -equations and write them down in a very concise and compactform which we have not found anywhere else in the literature. • We introduce a staggered-grid solver, for which we generatestencil code automatically from a computer algebra represen-tation of the P N -equations. This allows us to deal with the in-creasingly complex equations which the P N -method produces forhigher order. It further allows our solver to be used for any (po-tentially coupled) partial differential equations, which result in asystem of linear equations after discretization. • Finally, we compare the P N -method for higher orders againstflux-limited diffusion and ground truth Monte Carlo integration.In the next section, we will discuss related work and its relationto our contribution. In Section 3 we revisit the deterministic ap-proach to light transport simulation in participating media and out-line the discretization using spherical harmonics. In Section 4 weintroduce our computer algebra representation, which we requiredto derive the real-valued P N -equations, presented in Section 5. Thisrepresentation is also a key component of our solver, which wepresent in Section 6. Section 7 discusses application of the solu-tion in the context of rendering. We compare our P N -solver againstflux-limited diffusion for a set of standard problems in Section 8.Finally, Section 9 concludes with a summary and review of futurework.
2. Previous work
Light transport in participating media is governed by theRTE, first studied in the context of astrophysics by Chan-drasekhar [Cha60] and later introduced to computer graphics byKajiya [Kaj86]. In computer graphics today, this equation is typ-ically solved using Monte Carlo methods [NGHJ18]. However instrongly scattering or highly anisotropic media these methods canbecome prohibitively expensive, for example in the case of a highalbedo medium such as milk where tracing paths with a huge num-ber of scattering events is necessary.In contrast to path-tracing, the P N -method [Bru02] gives a so-lution by solving a system of linear equations. It is derived bydiscretizing the angular variable of the radiative transfer equationinto spherical harmonics (SH). This gives rise to a set of cou-pled, complex-valued partial differential equations, called the P N -equations. The subscript N refers to the spherical harmonics trun-cation order.The P N -method has a long history in other fields and was never Figure 2:
Solving the 2D checkerboard problem using naive col-located grids produces oscillating artifacts (left). Our solver usesstaggered grids and produces artifact free results (right). applied in graphics. Kajiya [KVH84] explained the theory, but didnot give any details on implementation or how to solve it. In fact,as Max [Max95] pointed out, it is not clear if Kajiya succeeded atall at applying the method, as all of the results in his paper wereproduced with a simpler method. This is further strengthened bythe fact that a straightforward finite difference discretization of the P N -equations produces unusable results, due to oscillation artifactsin the solution [SF14]. We use a staggered-grid solver, motivatedby the solver of Seibold et al. [SF14], that produces artifact-freesolutions (see Figure 2).Related to the P N -method, the classical diffusion approxima-tion (CDA) is a deterministic method which arrives at a solutionby solving a system of linear equations. It corresponds to the P N -equations when N = ∗ ∗ P N -method should provide an increasingly accu-rate solution to the full RTE including anisotropy as the truncationmoment order is increased, at the expense of more computationtime, it is still an open and unresolved question how high a trunca-tion order is needed for a particular problem in order to produce sig-nificantly better results than first order non-linear diffusion meth-ods (FLD and VEF). This question has also been raised in otherdomains [OAH00]. c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media
3. Discretized Radiative Transfer Equation
Our method derives from the radiative transfer equation (RTE),which expresses the change of the radiance field L , with respect toan infinitesimal change of position in direction ω at point (cid:126) x : ( ∇ · ω ) L ( (cid:126) x , ω ) = − σ t ( (cid:126) x ) L ( (cid:126) x , ω )+ σ s ( (cid:126) x ) (cid:90) Ω p (cid:0) ω (cid:48) · ω (cid:1) L (cid:0) (cid:126) x , ω (cid:48) (cid:1) d ω (cid:48) + Q ( (cid:126) x , ω ) . The left hand side (LHS) is the transport term, and we refer tothe terms on the right hand side (RHS) as collision, scattering, andsource term, respectively. The symbols σ t , σ s , p , and Q refer to theextinction- and scattering coefficient, phase function and emission.The RTE is often given in operator notation, where transport,collision, and scattering are expressed as operators T , C and S ,which are applied to the radiance field L : T ( L ) = −C ( L ) + S ( L ) + Q . (1)Deterministic methods are derived by discretizing the angularand spatial domain. This gives rise to a linear system of equations,which can be solved using standard methods. For the P N -method,the angular variable is first discretized, using a truncated sphericalharmonics expansion. This results in the P N -equations, a system ofcoupled PDEs that still depend on a continuous spatial variable.The number of equations grows with the truncation order N . Thisis why the discretization is laborious and difficult to do without er-rors if done by hand. We therefore use a computer algebra repre-sentation to automate this process. After giving an outline of thegeneral discretization in this section, we will present our computeralgebra representation in the next section. The P N -equations thatresult from our automated discretization are given in Section 5.Since the radiance field L is real, we use the real-valued SH basisfunctions Y l , m R , which are defined in terms of the complex-valuedSH basis functions Y l , m C as follows [JCJ09]: Y l , m R = i √ (cid:16) Y l , m C − ( − ) m Y l , − m C (cid:17) , for m < Y l , m C , for m = √ (cid:16) Y l , − m C − ( − ) m Y l , m C (cid:17) for m > P : P l , m ( f ) = (cid:90) Ω f ( (cid:126) x , ω ) Y l , m R ( ω ) d ω = f l , m ( (cid:126) x ) . The P N -equations are derived by first expressing all direction-dependent parameters in spherical harmonics. The radiance field L in Equation 1 is therefore replaced by its SH reconstruction (cid:98) L ,introducing an error due to truncation at order N : (cid:98) L ( (cid:126) x , ω ) = N ∑ l = l ∑ m = − l L l , m ( (cid:126) x ) Y l , m R ( ω ) ≈ L ( (cid:126) x , ω ) . After substitution, all angular parameters are expressed in terms
SH coefficients per voxelSpatial discretization(FD grid) Y × X Y × X Y × X Figure 3:
Structure of coefficient matrix A and solution vector (cid:126) uafter discretization of the P N -equations on a finite difference grid. of spherical harmonics, but they still depend on the continuous an-gular variable ω . As a next step, we project each term of the RTEinto spherical harmonics, using the projection operator P . This pro-duces a single equation for each l , m -pair. The P N -equations there-fore can be written as: P l , m T (cid:16)(cid:98) L (cid:17) = −P l , m C (cid:16)(cid:98) L (cid:17) + P l , m S (cid:16)(cid:98) L (cid:17) + P l , m ( Q ) . (3)Once the P N -equations have been found, the spatial variable (cid:126) x isdiscretized using a finite difference (FD) voxel grid (using centraldifferences for differential operators).Following this discretization, the radiance field L , is representedas a set of SH coefficients per voxel. Flattening these over all vox-els into a single vector, gives the solution vector (cid:126) u . The RHS vector (cid:126) Q is produced similarly. The projected operators can be expressedas linear transformations, which can be collapsed into a single co-efficient matrix A (see Figure 3): ( T + C − S ) (cid:126) u = A (cid:126) u = (cid:126) Q . (4) T , C , S are matrices, which result from the discretized transport,collision and scattering operators respectively.
4. Computer Algebra Representation
So far, we have only given the P N -equations in high-level oper-ator notation (Equation 3). Carrying out the full derivation createslarge, unwieldy equations and requires a string of expansions andapplications of identities. These are challenging to manipulate ifdone by hand. We therefore used a computer algebra representa-tion, which allowed us to derive and discretize the P N -equation ina semi-automatic fashion (Figure 4)).It represents the equations using a tree of mathematical ex-pressions, which represent numbers, symbols and other expressiontypes, such as integrals, derivatives, sums, products and functions.Further, manipulators can be executed on these expression trees toperform substitution, constant folding, reordering of nested inte-grals, application of identities and more complex operations. Fi-nally, frontends allow rendering the expression tree into differentforms, such as L A P N -equations c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media ) b + a ( y x + ba representation expansionrendering ∗ + b y x a + y x ∗ a b Figure 4:
A computer algebra framework allows us to representequations as mathematical expressions trees. It further provides aset of functions for manipulating the tree according to valid math-ematical operations, such as the binomial expansion above. Fron-tends allow generation of source code from the expression tree. (the derivation steps shown in the Appendix were almost all ren-dered from the expression tree). More importantly, we use the rep-resentation to perform the discretization and generate the stencilcode used by our solver. This is detailed in section 6.1.
5. Real-valued P N -Equations With the help of a computer algebra representation framework,we are able to easily derive and work with the large and unwieldy P N -equations. We present here the final real-valued P N -equations,for a general N and in three dimensions, in a very compact formwhich we have not found elsewhere in the literature. The derivationof the real-valued equations is rather long and takes many pagesto describe in detail. We therefore give only the final result in thissection, and present the full derivation for reference in Appendix A.Since the real-valued SH bases (Equation 2) have different defi-nitions for m < m = m >
0, we get different projections for S l , m , depending on the sign of m .For m = √ c l − , − ∂ x L l − , − √ d l + , − ∂ x L l + , √ c l − , − ∂ y L l − , − − √ d l + , − ∂ y L l + , − + a l − , ∂ z L l − , + b l + , ∂ z L l + , + σ t L l , m − σ s λ l p l , L l , m = Q l , m . (5)For m < m > c l − , ± m − ∂ x L l − , m ∓ − d l + , ± m − ∂ x L l + , m ∓ − β m x e l − , m ± ∂ x L l − , m ± + β m x f l + , ± m + ∂ x L l + , m ± ∓ c l − , ± m − ∂ y L l − , − m ± ± d l + , ± m − ∂ y L l + , − m ± ∓ β m y e l − , ± m + ∂ y L l − , − m ∓ ± β m y f l + , ± m + ∂ y L l + , − m ∓ + a l − , ± m ∂ z L l − , ∓ m + b l + , ± m ∂ z L l + , ∓ m + σ t L l , m − σ s λ l p l , L l , m = Q l , m , (6) with β mx = , for m = − √ , for m (cid:54) = , otherwise , β my = √ , for m = − , for m (cid:54) = , otherwiseand a l , m = (cid:115) ( l − m + ) ( l + m + )( l + ) ( l − ) b l , m = (cid:115) ( l − m ) ( l + m )( l + ) ( l − ) c l , m = (cid:115) ( l + m + ) ( l + m + )( l + ) ( l + ) d l , m = (cid:115) ( l − m ) ( l − m − )( l + ) ( l − ) e l , m = (cid:115) ( l − m + ) ( l − m + )( l + ) ( l + ) f l , m = (cid:115) ( l + m ) ( l + m − )( l + ) ( l − ) λ l = (cid:114) π l + . In the next section we will present the solver that we use to solvethe P N -equations. P N -Solver The truncation order N is the key input parameter to the solver.With higher values, manual implementation of the solver from theequations would be arduous, error-prone and time-consuming. Wetherefore decided to make use of the computer algebra representa-tion and designed our solver around it.The solver consists of two components. The first is a precompu-tation (Section 6.1), which is executed once for every single valueof N . This step runs a partial evaluation on the mathematical ex-pression tree and applies the spatial discretization in a referencespace we call stencil space. The precomputation step automaticallygenerates source code from the resulting expression tree.The generated stencil code is compiled with the runtime com-ponent (Section ) of our solver. This component receives the ac-tual problem as input, including grid resolution and RTE parameterfields. It then builds the linear system and solves it using standardmethods. An overview of the solver is given in Figure 5. The result of the precomputation is a stencil, which can be usedduring runtime to build the linear system for a given problem. Thestencil is a pattern of indices into the solution vector, along withvalues. It expresses how the sum of the weighted solution vectorcomponents relate to the RHS for a given unknown in the system,and therefore contains most information required to fill the systemmatrix A and RHS-vector (cid:126) Q row-by-row. Note that while the coef-ficients may change, the sparsity pattern of the stencil is identicalfor different rows. c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media stencilcode code generation ·· φ ∗ + ... h ,j + i Precomputation ) (cid:25)x ( q ) = (cid:25)x ( φ ∆ Runtime renderingbuild solve
Q(cid:25) = A(cid:25)u (cid:25)u
Figure 5:
Overview of our P N -solver. After generating the stencil source code from the expression trees representing the P N -equations, thelinear system A (cid:126) u = (cid:126) Q 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 (cid:126) u, which is then used in our rendering application.
Stencil generation entails discretizing a PDE at a hypotheticalcenter voxel ( i , j , k ) (assumed to mean the voxel center most of thetime). Finite differences create weighted references to other voxels(e.g. i + , j , k ). After bringing the discretized equation into canon-ical form (a weighted sum of unknowns), one can write the stencilby reading off the weights and offsets (Figure 6). Voxel ( i , j , k ) willonly be known during runtime, when the stencil is executed for aparticular unknown (row). Then the offsets can be used to find thecomponent index into the solution-vector, and weights can be eval-uated for concrete world space position. We refer to the space withthe hypothetical voxel ( i , j , k ) at the center as stencil space. ) (cid:31)x ( q = ... ) + (cid:31)x ( φ ∂x∂∂x∂ h φ i +1 − h φ i + 1 h φ i − + ... = q i,j φ i +1 − φ i h − φ i − φ i − h h + ... = q i,j code generationcanonical formdiscretization Figure 6:
Creating the stencil code requires several steps, usuallydone by hand. We express the given problem in a computer algebrarepresentation and use it to fully automate the process. x ∂ ·· − φh ∗ (cid:28)x ,j − i,j − i − ·· φh ∗ φ discretization The spatial dis-cretization is doneby parsing the ex-pression tree fromthe root. The dis-crete substitute forthe continuous po-sition variable (cid:126) x isinitialized with i jk . Differential operator nodes are replaced by asubtree, which expresses the finite difference approximation (in-cluding voxel-size factor h ). The subtree of the differential opera-tor node is duplicated for different offsets to the discrete variable ( i , j , k ) . Since this offset only applies to the subtree, a stack is main-tained by the parser to manage scope. Whenever the parser encoun-ters the continuous variable (cid:126) x in the expression tree, its node in thetree is replaced by the discrete substitute, currently on top of thestack. Nested differential operators yield higher order finite differ-ence stencils as expected. ·· φ ∗ code generation + ... h ,j + i Factorizationinto canonical formis done as a ma-nipulation pass onthe mathematicalexpression tree.The result allowsus to implement thecode generation ina straightforward fashion. For each term, the i jk -offset is retrievedfrom the unknown. The factor-expression, which is multipledwith the unknown, is extracted from the tree and used to generatesource code for evaluating the factor-expression during runtime(including calls for evaluating RTE-parameter fields). staggered grid interpolation ) i, j, k ( Our solver sup-ports placement ofcoefficients at arbi-trary staggered gridlocations (see Fig-ure 7). This meansthat during the dis-cretization step, the discrete location ( i , j , k ) (at which the unknownis meant to be evaluated) might not coincide with the location ofthe unknown. To solve this, depending on how the two are locatedrelative to each other, the parser returns an expression which in- c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media active unknowns witha row in Ax=b voxel centervoxelface in xvoxelface in yvoxel intersectionboundary voxelinterior voxelghost coefficients forboundary conditions Figure 7:
Staggering distributes the coefficients of the solution vec-tor (cid:126) u onto four disjoint grids, indicated by the symbols, in a waywhich guarantees second order accuracy [SF14]. terpolates the coefficients at ( i , j , k ) from their defined locations. Ifthose happen to coincide, the coefficient itself is returned. This isalso done for RTE parameters, such as σ t or p l , m , which are alwayslocated at the voxel center. We use the staggering scheme intro-duced by Seibold et al. [SF14] for their solver StaRMAP , whichthey proved ensures second order accuracy, i.e. a truncation errorof O ( h ) , and prevents the growth of non-physical oscillations inthe solution (see their paper for the full details of the staggered co-efficient placement). The stencil code is generated once for every value of N and com-piled with the runtime component of our solver. The runtime exe-cutes the stencil for every voxel to populate the system matrix A and RHS vector (cid:126) Q with numerical values.The number of rows is determined by the number of voxels timesthe number of coefficients per voxel (see Figure 3) and can there-fore become very large for high resolution and high truncation or-der. The matrix A is square and extremely sparse, due to the localstructure of the finite differences discretization. Unfortunately, itis non-symmetric due to the transport term and not diagonal dom-inant, which rules out many iterative methods for solving linearsystems. Iterative methods are useful, as they allow balancing accu-racy against performance by tweaking the convergence threshold.We address this by solving the normal form A T A (cid:126) u = A T (cid:126) Q instead.This gives a symmetric and positive definite system matrix A T A ,albeit with a higher condition number. Investigation of other solu-tion schemes (e.g. multigrid) would be an interesting avenue for fu-ture work. However, more importantly, in the presence of vacuumregions, the matrix A becomes singular and the system cannot besolved at all. This requires the introduction of a minimum thresh-old for the extinction coefficient σ t .Our solver supports both Neumann and Dirichlet boundary con-ditions (BC). They are handled transparently by the code whichgenerates the stencil. Whenever the stencil accumulates coefficientsinto a boundary location, the framework either ignores the write op-eration (Dirichlet BC) or accumulates into the row and column in A of the coefficient in the closest voxel inside the domain (NeumannBC). This is done by changing the index of the written component.
7. Rendering
We use an approach similar to Koerner et al. [KPS ∗ L ss andmultiple scattered light L ms : L ( (cid:126) x , ω ) = L ss ( (cid:126) x , ω ) + L ms ( (cid:126) x , ω ) . (7)The single scattered light is folded into the emission term Q : Q ( (cid:126) x , ω ) = L ss ( (cid:126) x , ω ) = σ s ( (cid:126) x ) (cid:90) Ω p (cid:0) ω (cid:48) → ω (cid:1) L u (cid:0) (cid:126) x , ω (cid:48) (cid:1) d ω (cid:48) . (8)This means that our solver will solve for the multiple scatteredlight L ms . The quantity L u is the “uncollided” light, which was emit-ted from the light source and attenuated by the volume without un-dergoing any scattering event. We compute it using a few light sam-ples per voxel, which quickly converges to a useful result for Diracdelta light sources.Running the solver gives solution vector (cid:126) u . We then un-staggerthe solution by interpolating all coefficients to voxel centers. Theadditional 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 N at the cen-ter 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 L ss . Then we sample anew direction according to the phase function. Instead of contin-uing tracing into the new direction, we evaluate the in-scatteringintegral using (cid:98) L ms . The SH coefficients at (cid:126) x are found by trilinearinterpolation from the voxel grid of SH coefficients.
8. Results
In this section, we present results for a set of standard prob-lems in order to validate and evaluate our method. We also compareagainst classical diffusion (CDA) and flux-limited diffusion (FLD).Our computer algebra framework and the precomputation hasbeen implemented in Python. The runtime component of our solverhas been implemented in C++. We use a naive implementation ofa CG solver, which has been modified such that we do not need toexplicitly compute the normal form of the linear system to solve.We use the sparse matrix representation and sparse matrix vectorproduct from the Eigen linear algebra library (eigen.tuxfamily.org).The solver for classical diffusion is based on the diffusion equa-tion, which is derived by collapsing the P -equations: ∇ (cid:18) σ t ∇ L , (cid:19) = − Q , . (9)Since our solver can work with any PDE which results in a linearsystem of equations, we put Equation 9 into our computer algebrarepresentation and provide it as an input to our solver, which gen-erates the correct stencil code automatically.Since FLD is based on a non-linear diffusion equation, we werenot able to use our system in the same way. Our implementationclosely follows the implementation in [KPS ∗
14] (though ours runson CPU) and we refer to their paper for more details.
First we ran our solver on the 2D checkerboard, a very common c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media (a) P N vs. ground truth (b) P vs. CDA and FLD Figure 8:
Lineplot through the 3D solution of our solver for thepoint source problem for various order N (left). Solution for P compared against classical diffusion, flux-limited diffusion and an-alytical solution (right). test case in other fields. The problem has dimensions 7 × ×
71. Unit size blocks are filled withpurely absorbing medium σ a =
10 in a checkerboard pattern. Allother blocks are filled with a purely scattering medium with σ s = StaRMAP solver from Seibold et al. [SF14],which solves for the time-dependent and complex-valued P N -equations on a staggered grid, but in the 2D case only. The 2Dcase is derived by assuming that all SH coefficients, RTE parame-ters and boundary conditions are z-independent. This causes all SHcoefficients and moment equations for which l + m is odd to vanish.Due to the time-dependency, their approach is to do explicit incre-mental steps in time. We run their solver for many timesteps to geta result close to steady state. Figure 9:
Comparison of the result (here the L , coefficient, or fluence ) for the checkerboard test using StaRMAP ’s time-steppingsolver [SF14] (left) against our steady-state solver (right) with N = . Our results are in good agreement. As can be seen in Figure 9, the results from our solver are in goodagreement with the results from Seibold et al. [SF14] and verifythe correctness of our implementation. Converging to a residual of10 e − takes 0 . s for P and 25 s for P . We also run our solver for the point source problem, a single
Figure 10:
Convergence behaviour of our solver with N = forthe nebulae dataset and for varying minimum thresholds of the ex-tinction coefficient σ t . Threshold values and an estimate for thecondition number of A (MATLAB’s condest function) are shownnext to the plots. The convergence deteriorates as the threshold de-creases. Once it reaches zero, the presence of pure vacuum makesthe condition number infinite. point light in a homogeneous medium. This not only helps to val-idate our implementation for the 3D case, but also provides infor-mation on the accuracy of these methods. We use the Grosjean ap-proximation, which was introduced by D’Eon et al. [dI11] as a veryaccurate approximation to the ground truth solution.For our test case, we choose a FD resolution of 80 × × σ t = . α = .
9. We run thesolver for different truncation values N . In Figure 8 a, we see thatthe solution becomes increasingly accurate for higher truncationorder. The ground truth is underestimated when N is odd, and over-estimated when N is even. We further see in 8 b that the P solutionexactly matches the results from CDA, which confirms that the lat-ter is only a collapsed version of the former. The time to solve issignificant with 10 m for P and 45 m with P . With these perfor-mance characteristics, our P N -solver is clearly not competitive incomparison with our CDA and FLD solver, which are much faster. Finally, we run our solver on a procedural cloud dataset to getan idea of its performance in more practical applications. Figure 1bshows the result of P for a procedurally generated heterogeneouscloud with an isotropic phase function, with path-traced groundtruth in Figure 1d. We see that at order N =
5, our method cancapture indirect illumination similarly well as FLD (1c) and is sig-nificantly better than CDA (1a) as expected. The indirectly illumi-nated region at the bottom appears to be closer to the path-tracedresult as opposed to the solution from FLD which is very flat in thatregion. However, in many other areas, P seems to still suffer a lotfrom the problem of energy loss near transitions from dense regionsto vacuum. It appears that going higher order a few steps mitigatesthis only slowly at a high cost of compute and storage. The maincharacteristic of the nebulae dataset is the presence of vacuum. Wefound that having vacuum regions in the dataset will cause the con-dition number to become infinite and the solver practically does notconverge. We therefore introduced a minimum threshold for the ex-tinction coefficient σ t . Every voxel with an extinction coefficientsmaller than the threshold is set to the threshold value. In Figure 10we show the effect of the minimum threshold on the convergencebehaviour. As it increases, convergence improves. c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media
9. Conclusion
In this paper we introduced the P N -method to the toolbox of de-terministic methods for rendering participating media in computergraphics. We derived and presented the real-valued P N -equations,along with a staggered grid solver with numerical stencils con-structed automatically from the equations via a computer algebrasystem. We showed how to use the results in a rendering systemand ran our solver for various standard problems for validation.We originally set out to understand how non-linear diffusionmethods compare to the P N -method for increasing order. Our re-sults indicate that although the lack of higher moments makes theFLD solution overly smooth, its energy conserving nature and com-parably small resource footprint make it a better approach at presentfor most graphics applications compared to the P N -method, whichbecomes increasingly costly for higher values of N .The literature in other fields often states that the P N method—when solved in normal form as we do—is able to deal with vacuumregions. We found this misleading. The P N -method in normal formindeed does not break down completely in the presence of vacuumas diffusion based methods do (due to σ − t in the diffusion coeffi-cient). However, in the presence of vacuum, the condition numberof the system matrix becomes infinite and does not converge either.Therefore P N based methods also require minimum thresholding ofthe extinction coefficient and offer no benefit for vacuum regions.Much more work needs to be done in order to make the P N -method competitive in performance to alternative solutions for vol-ume rendering. We believe this can be made possible by employ-ing a multigrid scheme for solving the linear system of equations.We implemented a multigrid solver, but did not find the expectedperformance improvements. This is possibly due to the couplingbetween coefficients within a voxel, which does not work well to-gether with the smoothing step. We want to study this in the future.Unique to our system is that it uses a computer algebra represen-tation of the equation to solve as input. Discretization in angularand spatial domain is done using manipulation passes on the repre-sentation. The stencil code, which is used to populate the system oflinear equations, is generated from the expression tree. This way,we can easily deal with coupled PDEs and avoid the time consum-ing and error prone process of writing stencil code by hand.When researching the application of the P N -method in otherfields, we came across a rich variety of variations, such assimplified P N ( SP N ) [McC10], filtered P N ( FP N ) [RARO13],diffusion-correction P N ( DP N ) [SFL11] and least-squares P N ( LSP N ) [HPM ∗ P N -method, such as ringingartifacts, dealing with vacuum regions and general convergence.Our solver can be applied to any (potentially coupled) PDE andtherefore can generate stencil code for all these variations by sim-ply expressing the respective PDEs in our computer algebra repre-sentation and providing this as an input to our solver. This wouldallow an exhaustive comparison of all these methods and we con-sider this as future work.Finally, since the approach of our solver is very generic, we alsowould like to explore its application to other simulation problemsin computer graphics. Acknowledgements
We thank the anonymous reviewers for their valuable and encour-aging comments and feedback. We also thank Martin Frank andBenjamin Seibold for the very valuable discussions and answers toour questions. We thank Bernd Eberhardt for feedback and support.This project has been partially funded by the MSC-BW project.
References [Bru02] B
RUNNER
T. A.: Forms of Approximate Radiation Transport.
Tech. Rep. SAND2002-1778, Sandia National Laboratories (2002). 1, 2[Cha60] C
HANDRASEKHAR
S.:
Radiative Transfer . Dover Publications,1960. 2[dI11] D ’E ON E., I
RVING
G.: A Quantized-Diffusion Model for Render-ing Translucent Materials.
ACM TOG (Proc. of SIGGRAPH) 30 , 4 (July2011), 56:1–56:14. 7[HPM ∗
14] H
ANSEN
J., P
ETERSON
J., M
OREL
J., R
AGUSA
J., W
ANG
Y.: A Least-Squares Transport Equation Compatible with Voids.
Journalof Computational and Theoretical Transport 43 , 1-7 (2014), 374–401. 8[JCJ09] J
AROSZ , W., C
ARR , N. A., J
ENSEN , H. W.: Importance Sam-pling Spherical Harmonics.
Computer Graphics Forum (Proceedings ofEurographics) (2009), 577–586. 3[Kaj86] K
AJIYA
J. T.: The rendering equation.
Computer Graphics(Proc. of SIGGRAPH) (1986), 143–150. 2[KPS ∗
14] K
OERNER
D., P
ORTSMOUTH
J., S
ADLO
F., E
RTL
T., E
BER - HARDT
B.: Flux-limited Diffusion for Multiple Scattering in Participat-ing media.
CoRR abs/1403.8105 (2014). 2, 6[KVH84] K
AJIYA
J. T., V ON H ERZEN
B. P.: Ray tracing volume densi-ties.
Computer Graphics (Proc. of SIGGRAPH) 18 , 3 (Jan. 1984), 165–174. 2[LP81] L
EVERMORE
C. D., P
OMRANING
G. C.: A flux-limited diffu-sion theory.
Astrophysical Journal 248 (1981), 321–334. 2[Max95] M AX N.:
Efficient Light Propagation for Multiple AnisotropicVolume Scattering . Springer Berlin Heidelberg, Berlin, Heidelberg,1995, pp. 87–104. 2[McC10] M C C LARREN
R. G.: Theoretical Aspects of the Simplified PnEquations.
Transport Theory and Statistical Physics 39 , 2-4 (2010), 73–109. 8[MGN17] M
ULLER , T., G
ROSS , M., N
OVAK , J.: Practical Path Guid-ing for Efficient Light Transport Simulation.
Computer Graphics Forum (2017), 36:91–36:100 1[NGHJ18] N
OVAK , J., G
EORGIEV , I., H
ANIKA , J., J
AROSZ , W.: MonteCarlo Methods for Volumetric Light Transport Simulation.
ComputerGraphics Forum (Proceedings of Eurographics - State of the Art Reports) (2018), 37(2). 1, 2[OAH00] O
LSON
G. L., A
UER
L. H., H
ALL
M. L.: Diffusion, P1, andother approximate forms of radiation transport.
Journal of QuantitativeSpectroscopy and Radiative Transfer 64 , 6 (2000), 619 – 634. 2[RARO13] R
ADICE
D., A
BDIKAMALOV
E., R
EZZOLLA
L., O TT C. D.:A New Spherical Harmonics Scheme for Multi-Dimensional RadiationTransport I. Static Matter Configurations.
Journal of ComputationalPhysics 242 (2013), 648 – 669. 8[SF14] S
EIBOLD
B., F
RANK
M.: StaRMAP—A second order stag-gered grid method for spherical harmonics moment equations of radia-tive transfer.
ACM Trans. Math. Softw. 41 , 1 (Oct. 2014), 4:1–4:28. 2, 6,7[SFL11] S
CHÄFER
M., F
RANK
M., L
EVERMORE
C. D.: Diffusive Cor-rections to Pn Approximations.
Multiscale Modeling and Simulation 9 (2011), 1–28. 8[Sta95] S
TAM
J.: Multiple scattering as a diffusion process. In
Proc.of Eurographics Workshop on Rendering Techniques . Springer-Verlag,1995, pp. 41–50. 2 c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media A. Full derivation of the P N -equations This appendix presents the detailed derivations of the complex-valued and real-valued P N -equations. This was actually performed semi-automatically using our computer algebra representation of the equations, guaranteeing correctness, and here we just report the result of eachstage in the derivation.The starting point is the radiative transfer equation (RTE), which expresses the change of the radiance field L , with respect to an infinites-imal change of position into direction ω at point (cid:126) x : ( ∇ · ω ) L ( (cid:126) x , ω ) = − σ t ( (cid:126) x ) L ( (cid:126) x , ω )+ σ s ( (cid:126) x ) (cid:90) Ω p (cid:0) ω (cid:48) · ω (cid:1) L (cid:0) (cid:126) x , ω (cid:48) (cid:1) d ω (cid:48) . + Q ( (cid:126) x , ω ) where 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 sourceterm, respectively. The symbols σ t , σ s , p , and Q refer to the extinction coefficient, scattering coefficient, phase function and emission term.The derivation of the P N -equations is then done in two steps. First, the directional-dependent quantities are replaced by their SH-projectedcounterparts. For example the radiance field L is replaced by its SH projection. This way the quantities are expressed in spherical harmonics,but still depend on direction ω . In the second step, the RTE is projected into spherical harmonics, which is done by multiplying each termwith the complex conjugate of the SH basis functions.The SH basis functions are complex, which produces complex coefficients and complex P N -equations. However, there are also real SHbasis functions, which are defined in terms of the complex SH basis functions and which produce real coefficients and reconstructions. Sincethe radiance field L is real, it is more convenient to work with the real SH basis functions.In the next section, the complex P N -equations are derived. In order to give the derivation a clearer structure, the two steps mentioned aboveare applied to each term individually in a separate subsection. The section concludes by putting all derived terms together. In Section 11 theanalogous derivation is followed to obtain the real P N -equations which are used in the article. Contents
10 Derivation of the complex-valued P N -equations10.1 Projecting Radiative Transfer Quantities10.1.1 Radiance Field L and Emission Field Q P N -equations11.1 Projecting Radiative Transfer Quantities11.1.1 Radiance Field L and Emission Field Q c (cid:13) (cid:13)(cid:13)
10 Derivation of the complex-valued P N -equations10.1 Projecting Radiative Transfer Quantities10.1.1 Radiance Field L and Emission Field Q P N -equations11.1 Projecting Radiative Transfer Quantities11.1.1 Radiance Field L and Emission Field Q c (cid:13) (cid:13)(cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media
10. Derivation of the complex-valued P N -equations Deriving the P N -equations consists of two main steps. First, all angular dependent quantities in the RTE are expressed in terms of sphericalharmonics (SH) basis functions. After this, the RTE still depends on the angular variable. Therefore, the second step projects each term ofthe RTE by multiplying with the complex conjugate of the SH basis functions, followed by integration over solid angle to integrate out theangular 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 givenby Y l , m C ( ω ) = Y l , m C ( θ , φ ) = ( − ) m (cid:113) l + π ( l − m ) ! ( l + m ) ! e im φ P l , m ( cos ( θ )) , for m ≥ ( − ) m Y l | m | C ( θ , φ ) , for m < (10) where P l , m are the associated Legendre polynomials. The ( − ) m factor is called the Condon-Shortley phase and is not part of the associatedLegendre Polynomial (unlike some other definitions). L and Emission Field Q Radiative transfer quantities, which depend on position (cid:126) x and angle ω , are projected into spatially dependent SH coefficients for each SHbasis function: L l , m ( (cid:126) x ) = (cid:90) Ω L ( (cid:126) x , ω ) Y l , m C d ω . Q l , m ( (cid:126) x ) = (cid:90) Ω Q ( (cid:126) x , ω ) Y l , m C d ω . The function is completely reconstructed by using all SH basis functions up to infinite order. The P N -equations introduce a truncation errorby only using SH basis functions up to order N for the reconstruction ˆ L and ˆ Q : L ( (cid:126) x , ω ) ≈ ˆ L ( (cid:126) x , ω ) = N ∑ l = l ∑ m = − l L l , m ( (cid:126) x ) Y l , m C ( ω ) = ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) . Q ( (cid:126) x , ω ) ≈ ˆ Q ( (cid:126) x , ω ) = N ∑ l = l ∑ m = − l Q l , m ( (cid:126) x ) Y l , m C ( ω ) = ∑ l , m Q l , m ( (cid:126) x ) Y l , m C ( ω ) . (11) Throughout our article, we assume an isotropic phase function, which only depends on the angle between incident and outgoing vector ω i and ω o (note that in the graphics literature, these would often be called anisotropic). We will see later in section 10.2.3, that this allows us tofix the outgoing vector ω o at the pole axis (cid:126) e and compute the phase function SH coefficients by just varying the incident vector ω i . p l , m = (cid:90) Ω p ( ω i · (cid:126) e ) Y l , m C ( ω i ) d ω i . The expansion of the phase function can be further simplified because the phase function is rotationally symmetric around the pole axis (cid:126) e . Consider the definition of the spherical harmonics basis function Y l , m C : Y l , m C ( θ , φ ) = C l , m e im φ P l , m ( cos ( θ )) . Now we apply a rotation R ( α ) of α degrees around the pole axis. In spherical harmonics, this is expressed as: ρ R ( α ) ( Y l , m C ) = e − im α Y l , m C . If the phase function is rotationally symmetric around the pole axis, we have: ρ R ( α ) ( p ) = p . and in spherical harmonics this would be: ∑ l , m e − im α p l , m Y l , m C ( ω i ) = ∑ l , m p l , m Y l , m C ( ω i ) . By equating coefficients we get: p l , m = p l , m e − im α . c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media Since e − im α = α only when m =
0, we can conclude that p l , m = m (cid:54) =
0. This means that for a function which is rotationallysymmetric around the pole axis, only the m = ω o = (cid:126) e ) only requires SH coefficients with m = p ( ω i ) = ∑ l p l Y l C ( ω i ) . (12) The transport term of the RTE is given as ( ω · ∇ ) L ( (cid:126) x , ω ) Replacing L with its expansion gives: ( ω · ∇ ) (cid:32) ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) (cid:33) Next we multiply with Y l (cid:48) m (cid:48) C ( ω ) and integrate over solid angle: (cid:90) Ω Y l (cid:48) m (cid:48) ( ω )( ω · ∇ ) ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) d ω We can pull the spatial derivative out of the integral to get: ∇ · (cid:90) Ω ω Y l (cid:48) m (cid:48) C ( ω ) ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) d ω (13) We apply the following recursive relation for the spherical harmonics basis functions: ω Y l , m C = c l − , m − Y l − , m − C − d l + , m − Y l + , m − C − e l − , m + Y l − , m + C + f l + , m + Y l + , m + C i (cid:16) − c l − , m − Y l − , m − C + d l + , m − Y l + , m − C − e l − , m + Y l − , m + C + f l + , m + Y l + , m + C (cid:17) (cid:16) a l − , m Y l − , m C + b l + , m Y l + , m C (cid:17) (14) with a l , m = (cid:115) ( l − m + ) ( l + m + )( l + ) ( l − ) , b l , m = (cid:115) ( l − m ) ( l + m )( l + ) ( l − ) , c l , m = (cid:115) ( l + m + ) ( l + m + )( l + ) ( l + ) d l , m = (cid:115) ( l − m ) ( l − m − )( l + ) ( l − ) , e l , m = (cid:115) ( l − m + ) ( l − m + )( l + ) ( l + ) , f l , m = (cid:115) ( l + m ) ( l + m − )( l + ) ( l − ) Note that the signs for the x - and y - component depend on the handedness of the coordinate system in which the SH basis functions aredefined. Using this in Equation 13 gives ∂ xi ∂ y ∂ z · (cid:90) Ω c l (cid:48) − , m (cid:48) − Y l (cid:48) − , m (cid:48) − C − d l (cid:48) + , m (cid:48) − Y l (cid:48) + , m (cid:48) − C − e l (cid:48) − , m (cid:48) + Y l (cid:48) − , m (cid:48) + C + f l (cid:48) + , m (cid:48) + Y l (cid:48) + , m (cid:48) + C − c l (cid:48) − , m (cid:48) − Y l (cid:48) − , m (cid:48) − C + d l (cid:48) + , m (cid:48) − Y l (cid:48) + , m (cid:48) − C − e l (cid:48) − , m (cid:48) + Y l (cid:48) − , m (cid:48) + C + f l (cid:48) + , m (cid:48) + Y l (cid:48) + , m (cid:48) + C a l (cid:48) − , m (cid:48) Y l (cid:48) − , m (cid:48) C + b l (cid:48) + , m (cid:48) Y l (cid:48) + , m (cid:48) C ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) d ω Integrating the vector term over solid angle can be expressed as separate solid angle integrals over each component. These integrals overa sum of terms are split into separate integrals. We arrive at: ∂ xi ∂ y ∂ z · c l (cid:48) − , m (cid:48) − ∑ l , m L l , m ( (cid:126) x ) (cid:82) Ω Y l (cid:48) − , m (cid:48) − C ( ω ) Y l , m C ( ω ) d ω − ... − c l (cid:48) − , m (cid:48) − ∑ l , m L l , m ( (cid:126) x ) (cid:82) Ω Y l (cid:48) − , m (cid:48) − C ( ω ) Y l , m C ( ω ) d ω + ... a l (cid:48) − , m (cid:48) ∑ l , m L l , m ( (cid:126) x ) (cid:82) Ω Y l (cid:48) − , m (cid:48) C ( ω ) Y l , m C ( ω ) d ω + ... c (cid:13) (cid:13)(cid:13)
0. This means that for a function which is rotationallysymmetric around the pole axis, only the m = ω o = (cid:126) e ) only requires SH coefficients with m = p ( ω i ) = ∑ l p l Y l C ( ω i ) . (12) The transport term of the RTE is given as ( ω · ∇ ) L ( (cid:126) x , ω ) Replacing L with its expansion gives: ( ω · ∇ ) (cid:32) ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) (cid:33) Next we multiply with Y l (cid:48) m (cid:48) C ( ω ) and integrate over solid angle: (cid:90) Ω Y l (cid:48) m (cid:48) ( ω )( ω · ∇ ) ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) d ω We can pull the spatial derivative out of the integral to get: ∇ · (cid:90) Ω ω Y l (cid:48) m (cid:48) C ( ω ) ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) d ω (13) We apply the following recursive relation for the spherical harmonics basis functions: ω Y l , m C = c l − , m − Y l − , m − C − d l + , m − Y l + , m − C − e l − , m + Y l − , m + C + f l + , m + Y l + , m + C i (cid:16) − c l − , m − Y l − , m − C + d l + , m − Y l + , m − C − e l − , m + Y l − , m + C + f l + , m + Y l + , m + C (cid:17) (cid:16) a l − , m Y l − , m C + b l + , m Y l + , m C (cid:17) (14) with a l , m = (cid:115) ( l − m + ) ( l + m + )( l + ) ( l − ) , b l , m = (cid:115) ( l − m ) ( l + m )( l + ) ( l − ) , c l , m = (cid:115) ( l + m + ) ( l + m + )( l + ) ( l + ) d l , m = (cid:115) ( l − m ) ( l − m − )( l + ) ( l − ) , e l , m = (cid:115) ( l − m + ) ( l − m + )( l + ) ( l + ) , f l , m = (cid:115) ( l + m ) ( l + m − )( l + ) ( l − ) Note that the signs for the x - and y - component depend on the handedness of the coordinate system in which the SH basis functions aredefined. Using this in Equation 13 gives ∂ xi ∂ y ∂ z · (cid:90) Ω c l (cid:48) − , m (cid:48) − Y l (cid:48) − , m (cid:48) − C − d l (cid:48) + , m (cid:48) − Y l (cid:48) + , m (cid:48) − C − e l (cid:48) − , m (cid:48) + Y l (cid:48) − , m (cid:48) + C + f l (cid:48) + , m (cid:48) + Y l (cid:48) + , m (cid:48) + C − c l (cid:48) − , m (cid:48) − Y l (cid:48) − , m (cid:48) − C + d l (cid:48) + , m (cid:48) − Y l (cid:48) + , m (cid:48) − C − e l (cid:48) − , m (cid:48) + Y l (cid:48) − , m (cid:48) + C + f l (cid:48) + , m (cid:48) + Y l (cid:48) + , m (cid:48) + C a l (cid:48) − , m (cid:48) Y l (cid:48) − , m (cid:48) C + b l (cid:48) + , m (cid:48) Y l (cid:48) + , m (cid:48) C ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) d ω Integrating the vector term over solid angle can be expressed as separate solid angle integrals over each component. These integrals overa sum of terms are split into separate integrals. We arrive at: ∂ xi ∂ y ∂ z · c l (cid:48) − , m (cid:48) − ∑ l , m L l , m ( (cid:126) x ) (cid:82) Ω Y l (cid:48) − , m (cid:48) − C ( ω ) Y l , m C ( ω ) d ω − ... − c l (cid:48) − , m (cid:48) − ∑ l , m L l , m ( (cid:126) x ) (cid:82) Ω Y l (cid:48) − , m (cid:48) − C ( ω ) Y l , m C ( ω ) d ω + ... a l (cid:48) − , m (cid:48) ∑ l , m L l , m ( (cid:126) x ) (cid:82) Ω Y l (cid:48) − , m (cid:48) C ( ω ) Y l , m C ( ω ) d ω + ... c (cid:13) (cid:13)(cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media Applying the orthogonality property to the solid angle integrals will will select specific l , m in each term: ∂ xi ∂ y ∂ z · c l − , m − L l − , m − − d l + , m − L l + , m − − e l − , m + L l − , m + + f l + , m + L l + , m + − c l − , m − L l − , m − + d l + , m − L l + , m − − e l − , m + L l − , m + + f l + , m + L l + , m + a l − , m L l − , m + b l + , m L l + , m Which gives the final moment equation for the transport term: = ∂ x (cid:16) c l − , m − L l − , m − − d l + , m − L l + , m − − e l − , m + L l − , m + + f l + , m + L l + , m + (cid:17) + i ∂ y (cid:16) − c l − , m − L l − , m − + d l + , m − L l + , m − − e l − , m + L l − , m + + f l + , m + L l + , m + (cid:17) + ∂ z (cid:16) a l − , m L l − , m + b l + , m L l + , m (cid:17) = c l − , m − ∂ x L l − , m − − d l + , m − ∂ x L l + , m − − e l − , m + ∂ x L l − , m + + f l + , m + ∂ x L l + , m + + − i c l − , m − ∂ y L l − , m − + i d l + , m − ∂ y L l + , m − − i e l − , m + ∂ y L l − , m + + i f l + , m + ∂ y L l + , m + + a l − , m ∂ z L l − , m + b l + , m ∂ z L l + , m . The collision term of the RTE is given as: − σ t ( (cid:126) x ) L ( (cid:126) x , ω ) We first replace the radiance field L with its spherical harmonics expansion: − σ t ( (cid:126) x ) ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) Multiplying with Y l (cid:48) m (cid:48) C and integrating over solid angle gives, after pulling some factors out of the integral: − σ t ( (cid:126) x ) ∑ l , m L l , m ( (cid:126) x ) (cid:90) Ω Y l (cid:48) m (cid:48) C ( ω ) Y l , m C ( ω ) d ω = − σ t ( (cid:126) x ) ∑ l , m L l , m ( (cid:126) x ) δ ll (cid:48) δ mm (cid:48) = − σ t ( (cid:126) x ) L l , m ( (cid:126) x ) . The scattering term in the RTE is given as: σ s ( (cid:126) x ) (cid:90) Ω p ( (cid:126) x , ω (cid:48) · ω ) L ( (cid:126) x , ω (cid:48) ) d ω (cid:48) The phase function used in isotropic scattering medium only depends on the angle between incident and outgoing direction and thereforeis rotationally symmetric around the pole-defining axis. This property allows us to define a rotation R ( ω ) , which rotates the phase functionsuch that the pole axis aligns with direction vector ω . The rotated phase function is defined as: ρ R ( ω ) ( p ) where ρ is the rotation operator, which can be implemented by applying the inverse rotation R ( ω ) − to the arguments of p . With this rotatedphase function, we now can express the integral of the scattering operator as a convolution denoted with the symbol ◦ : (cid:90) Ω p ( (cid:126) x , ω (cid:48) · ω ) L ( (cid:126) x , ω (cid:48) ) d ω (cid:48) = L ◦ ρ R ( ω ) ( p )= (cid:90) Ω (cid:48) L ( (cid:126) x , ω (cid:48) ) ρ R ( ω ) ( p )( ω (cid:48) ) d ω (cid:48) = (cid:104) L , ρ R ( ω ) ( p ) (cid:105) . (15) As we evaluate the inner product integral of the convolution, the phase function rotates along with the argument ω . c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media We now use the spherical harmonics expansions of L (Equation 11) and p (Equation 23) in the definition for the inner product of ourconvolution (Equation 15): (cid:104) L , ρ R ( ω ) ( p ) (cid:105) = (cid:42) ∑ l , m L l , m ( (cid:126) x ) Y l , m C , ρ R ( ω ) (cid:32) ∑ l p l Y l C (cid:33)(cid:43) . Due to linearity of the inner product operator, we can pull out the non-angular dependent parts of the expansions: (cid:104) L , ρ R ( ω ) ( p ) (cid:105) = ∑ l , m L l , m ( (cid:126) x ) (cid:42) Y l , m C , ρ R ( ω ) (cid:32) ∑ l p l Y l C (cid:33)(cid:43) . and further: (cid:104) L , ρ R ( ω ) ( p ) (cid:105) = ∑ l (cid:48) ∑ l , m p l (cid:48) L l , m ( (cid:126) x ) (cid:68) Y l , m C , ρ R ( ω ) (cid:16) Y l (cid:48) C (cid:17)(cid:69) . The rotation ρ R ( ω ) of a function with frequency l gives a function of frequency l again. In addition the spherical harmonics basis functions Y l , m C are orthogonal. We therefore have: (cid:68) Y l , m C , ρ R ( ω ) (cid:16) Y l (cid:48) m (cid:48) C (cid:17)(cid:69) = l (cid:54) = l (cid:48) which further simplifies our inner product integral to: (cid:104) L , ρ R ( ω ) ( p ) (cid:105) = ∑ l , m p l L l , m ( (cid:126) x ) (cid:68) Y l , m C , ρ R ( ω ) (cid:16) Y l C (cid:17)(cid:69) . What remains to be resolved is the inner product. We use the fact that the spherical harmonics basis functions Y l , m C are eigenfunctions ofthe inner product integral operator in the equation above: (cid:68) Y l , m C , ρ R ( ω ) (cid:16) Y l C (cid:17)(cid:69) = λ l Y l , m C with λ l = (cid:114) π l + . Replacing the inner product gives: (cid:104) L , ρ R ( ω ) ( p ) (cid:105) = ∑ l , m λ l p l L l , m ( (cid:126) x ) Y l , m C . This allows us to express the scattering term using SH expansions of phase function p and radiance field L : σ s ( (cid:126) x ) (cid:90) Ω p ( (cid:126) x , ω (cid:48) · ω ) L ( (cid:126) x , ω (cid:48) ) d ω (cid:48) = σ s ( (cid:126) x ) (cid:104) L , ρ R ( ω ) ( p ) (cid:105) = σ s ( (cid:126) x ) ∑ l , m λ l p l L l , m ( (cid:126) x ) Y l , m C . However, we haven’t done a spherical harmonics expansion of the scattering term itself. It is still a scalar function which depends ondirection ω . We thus project the scattering term into spherical harmonics by multiplying with Y l (cid:48) m (cid:48) C and integrating over solid angle ω .We further pull out all factors which do not depend on ω , and apply the SH orthogonality property to arrive at the scattering term of thecomplex-valued P N -equations: (cid:90) Ω Y l (cid:48) m (cid:48) C ( ω ) σ s ( (cid:126) x ) ∑ l , m λ l p l L l , m ( (cid:126) x ) Y l , m C ( ω ) d ω = λ l σ s ( (cid:126) x ) p l L l , m ( (cid:126) x ) ∑ l , m (cid:90) Ω Y l (cid:48) m (cid:48) C ( ω ) Y l , m C ( ω ) d ω = λ l σ s ( (cid:126) x ) p l L l , m ( (cid:126) x ) ∑ l , m δ ll (cid:48) δ mm (cid:48) = λ l σ s ( (cid:126) x ) p l L l , m ( (cid:126) x ) . (16) c (cid:13) (cid:13)(cid:13)
0. This means that for a function which is rotationallysymmetric around the pole axis, only the m = ω o = (cid:126) e ) only requires SH coefficients with m = p ( ω i ) = ∑ l p l Y l C ( ω i ) . (12) The transport term of the RTE is given as ( ω · ∇ ) L ( (cid:126) x , ω ) Replacing L with its expansion gives: ( ω · ∇ ) (cid:32) ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) (cid:33) Next we multiply with Y l (cid:48) m (cid:48) C ( ω ) and integrate over solid angle: (cid:90) Ω Y l (cid:48) m (cid:48) ( ω )( ω · ∇ ) ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) d ω We can pull the spatial derivative out of the integral to get: ∇ · (cid:90) Ω ω Y l (cid:48) m (cid:48) C ( ω ) ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) d ω (13) We apply the following recursive relation for the spherical harmonics basis functions: ω Y l , m C = c l − , m − Y l − , m − C − d l + , m − Y l + , m − C − e l − , m + Y l − , m + C + f l + , m + Y l + , m + C i (cid:16) − c l − , m − Y l − , m − C + d l + , m − Y l + , m − C − e l − , m + Y l − , m + C + f l + , m + Y l + , m + C (cid:17) (cid:16) a l − , m Y l − , m C + b l + , m Y l + , m C (cid:17) (14) with a l , m = (cid:115) ( l − m + ) ( l + m + )( l + ) ( l − ) , b l , m = (cid:115) ( l − m ) ( l + m )( l + ) ( l − ) , c l , m = (cid:115) ( l + m + ) ( l + m + )( l + ) ( l + ) d l , m = (cid:115) ( l − m ) ( l − m − )( l + ) ( l − ) , e l , m = (cid:115) ( l − m + ) ( l − m + )( l + ) ( l + ) , f l , m = (cid:115) ( l + m ) ( l + m − )( l + ) ( l − ) Note that the signs for the x - and y - component depend on the handedness of the coordinate system in which the SH basis functions aredefined. Using this in Equation 13 gives ∂ xi ∂ y ∂ z · (cid:90) Ω c l (cid:48) − , m (cid:48) − Y l (cid:48) − , m (cid:48) − C − d l (cid:48) + , m (cid:48) − Y l (cid:48) + , m (cid:48) − C − e l (cid:48) − , m (cid:48) + Y l (cid:48) − , m (cid:48) + C + f l (cid:48) + , m (cid:48) + Y l (cid:48) + , m (cid:48) + C − c l (cid:48) − , m (cid:48) − Y l (cid:48) − , m (cid:48) − C + d l (cid:48) + , m (cid:48) − Y l (cid:48) + , m (cid:48) − C − e l (cid:48) − , m (cid:48) + Y l (cid:48) − , m (cid:48) + C + f l (cid:48) + , m (cid:48) + Y l (cid:48) + , m (cid:48) + C a l (cid:48) − , m (cid:48) Y l (cid:48) − , m (cid:48) C + b l (cid:48) + , m (cid:48) Y l (cid:48) + , m (cid:48) C ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) d ω Integrating the vector term over solid angle can be expressed as separate solid angle integrals over each component. These integrals overa sum of terms are split into separate integrals. We arrive at: ∂ xi ∂ y ∂ z · c l (cid:48) − , m (cid:48) − ∑ l , m L l , m ( (cid:126) x ) (cid:82) Ω Y l (cid:48) − , m (cid:48) − C ( ω ) Y l , m C ( ω ) d ω − ... − c l (cid:48) − , m (cid:48) − ∑ l , m L l , m ( (cid:126) x ) (cid:82) Ω Y l (cid:48) − , m (cid:48) − C ( ω ) Y l , m C ( ω ) d ω + ... a l (cid:48) − , m (cid:48) ∑ l , m L l , m ( (cid:126) x ) (cid:82) Ω Y l (cid:48) − , m (cid:48) C ( ω ) Y l , m C ( ω ) d ω + ... c (cid:13) (cid:13)(cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media Applying the orthogonality property to the solid angle integrals will will select specific l , m in each term: ∂ xi ∂ y ∂ z · c l − , m − L l − , m − − d l + , m − L l + , m − − e l − , m + L l − , m + + f l + , m + L l + , m + − c l − , m − L l − , m − + d l + , m − L l + , m − − e l − , m + L l − , m + + f l + , m + L l + , m + a l − , m L l − , m + b l + , m L l + , m Which gives the final moment equation for the transport term: = ∂ x (cid:16) c l − , m − L l − , m − − d l + , m − L l + , m − − e l − , m + L l − , m + + f l + , m + L l + , m + (cid:17) + i ∂ y (cid:16) − c l − , m − L l − , m − + d l + , m − L l + , m − − e l − , m + L l − , m + + f l + , m + L l + , m + (cid:17) + ∂ z (cid:16) a l − , m L l − , m + b l + , m L l + , m (cid:17) = c l − , m − ∂ x L l − , m − − d l + , m − ∂ x L l + , m − − e l − , m + ∂ x L l − , m + + f l + , m + ∂ x L l + , m + + − i c l − , m − ∂ y L l − , m − + i d l + , m − ∂ y L l + , m − − i e l − , m + ∂ y L l − , m + + i f l + , m + ∂ y L l + , m + + a l − , m ∂ z L l − , m + b l + , m ∂ z L l + , m . The collision term of the RTE is given as: − σ t ( (cid:126) x ) L ( (cid:126) x , ω ) We first replace the radiance field L with its spherical harmonics expansion: − σ t ( (cid:126) x ) ∑ l , m L l , m ( (cid:126) x ) Y l , m C ( ω ) Multiplying with Y l (cid:48) m (cid:48) C and integrating over solid angle gives, after pulling some factors out of the integral: − σ t ( (cid:126) x ) ∑ l , m L l , m ( (cid:126) x ) (cid:90) Ω Y l (cid:48) m (cid:48) C ( ω ) Y l , m C ( ω ) d ω = − σ t ( (cid:126) x ) ∑ l , m L l , m ( (cid:126) x ) δ ll (cid:48) δ mm (cid:48) = − σ t ( (cid:126) x ) L l , m ( (cid:126) x ) . The scattering term in the RTE is given as: σ s ( (cid:126) x ) (cid:90) Ω p ( (cid:126) x , ω (cid:48) · ω ) L ( (cid:126) x , ω (cid:48) ) d ω (cid:48) The phase function used in isotropic scattering medium only depends on the angle between incident and outgoing direction and thereforeis rotationally symmetric around the pole-defining axis. This property allows us to define a rotation R ( ω ) , which rotates the phase functionsuch that the pole axis aligns with direction vector ω . The rotated phase function is defined as: ρ R ( ω ) ( p ) where ρ is the rotation operator, which can be implemented by applying the inverse rotation R ( ω ) − to the arguments of p . With this rotatedphase function, we now can express the integral of the scattering operator as a convolution denoted with the symbol ◦ : (cid:90) Ω p ( (cid:126) x , ω (cid:48) · ω ) L ( (cid:126) x , ω (cid:48) ) d ω (cid:48) = L ◦ ρ R ( ω ) ( p )= (cid:90) Ω (cid:48) L ( (cid:126) x , ω (cid:48) ) ρ R ( ω ) ( p )( ω (cid:48) ) d ω (cid:48) = (cid:104) L , ρ R ( ω ) ( p ) (cid:105) . (15) As we evaluate the inner product integral of the convolution, the phase function rotates along with the argument ω . c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media We now use the spherical harmonics expansions of L (Equation 11) and p (Equation 23) in the definition for the inner product of ourconvolution (Equation 15): (cid:104) L , ρ R ( ω ) ( p ) (cid:105) = (cid:42) ∑ l , m L l , m ( (cid:126) x ) Y l , m C , ρ R ( ω ) (cid:32) ∑ l p l Y l C (cid:33)(cid:43) . Due to linearity of the inner product operator, we can pull out the non-angular dependent parts of the expansions: (cid:104) L , ρ R ( ω ) ( p ) (cid:105) = ∑ l , m L l , m ( (cid:126) x ) (cid:42) Y l , m C , ρ R ( ω ) (cid:32) ∑ l p l Y l C (cid:33)(cid:43) . and further: (cid:104) L , ρ R ( ω ) ( p ) (cid:105) = ∑ l (cid:48) ∑ l , m p l (cid:48) L l , m ( (cid:126) x ) (cid:68) Y l , m C , ρ R ( ω ) (cid:16) Y l (cid:48) C (cid:17)(cid:69) . The rotation ρ R ( ω ) of a function with frequency l gives a function of frequency l again. In addition the spherical harmonics basis functions Y l , m C are orthogonal. We therefore have: (cid:68) Y l , m C , ρ R ( ω ) (cid:16) Y l (cid:48) m (cid:48) C (cid:17)(cid:69) = l (cid:54) = l (cid:48) which further simplifies our inner product integral to: (cid:104) L , ρ R ( ω ) ( p ) (cid:105) = ∑ l , m p l L l , m ( (cid:126) x ) (cid:68) Y l , m C , ρ R ( ω ) (cid:16) Y l C (cid:17)(cid:69) . What remains to be resolved is the inner product. We use the fact that the spherical harmonics basis functions Y l , m C are eigenfunctions ofthe inner product integral operator in the equation above: (cid:68) Y l , m C , ρ R ( ω ) (cid:16) Y l C (cid:17)(cid:69) = λ l Y l , m C with λ l = (cid:114) π l + . Replacing the inner product gives: (cid:104) L , ρ R ( ω ) ( p ) (cid:105) = ∑ l , m λ l p l L l , m ( (cid:126) x ) Y l , m C . This allows us to express the scattering term using SH expansions of phase function p and radiance field L : σ s ( (cid:126) x ) (cid:90) Ω p ( (cid:126) x , ω (cid:48) · ω ) L ( (cid:126) x , ω (cid:48) ) d ω (cid:48) = σ s ( (cid:126) x ) (cid:104) L , ρ R ( ω ) ( p ) (cid:105) = σ s ( (cid:126) x ) ∑ l , m λ l p l L l , m ( (cid:126) x ) Y l , m C . However, we haven’t done a spherical harmonics expansion of the scattering term itself. It is still a scalar function which depends ondirection ω . We thus project the scattering term into spherical harmonics by multiplying with Y l (cid:48) m (cid:48) C and integrating over solid angle ω .We further pull out all factors which do not depend on ω , and apply the SH orthogonality property to arrive at the scattering term of thecomplex-valued P N -equations: (cid:90) Ω Y l (cid:48) m (cid:48) C ( ω ) σ s ( (cid:126) x ) ∑ l , m λ l p l L l , m ( (cid:126) x ) Y l , m C ( ω ) d ω = λ l σ s ( (cid:126) x ) p l L l , m ( (cid:126) x ) ∑ l , m (cid:90) Ω Y l (cid:48) m (cid:48) C ( ω ) Y l , m C ( ω ) d ω = λ l σ s ( (cid:126) x ) p l L l , m ( (cid:126) x ) ∑ l , m δ ll (cid:48) δ mm (cid:48) = λ l σ s ( (cid:126) x ) p l L l , m ( (cid:126) x ) . (16) c (cid:13) (cid:13)(cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media The emission term of the RTE is given as: Q ( (cid:126) x , ω ) (17) The derivation of the SH projected term is equivalent to the derivation of the projected collision term. Replacing the emission field with itsSH projection and multiplying the term with the conjugate complex of Y C results, after applying the orthogonality property, in: Q l , m ( (cid:126) x , ω ) (18) We arrive at the complex-valued P N -equations after putting all the projected terms together:12 c l − , m − ∂ x L l − , m − − d l + , m − ∂ x L l + , m − − e l − , m + ∂ x L l − , m + + f l + , m + ∂ x L l + , m + + − i c l − , m − ∂ y L l − , m − + i d l + , m − ∂ y L l + , m − − i e l − , m + ∂ y L l − , m + + i f l + , m + ∂ y L l + , m + + a l − , m ∂ z L l − , m + b l + , m ∂ z L l + , m = − σ t ( (cid:126) x ) L l , m ( (cid:126) x ) + λ l σ s ( (cid:126) x ) p l L l , m ( (cid:126) x ) + Q l , m ( (cid:126) x , ω ) .
11. Derivation of the real-valued P N -equations The real-valued P N -equations are derived similar to their complex-valued counterpart, except that the real-valued SH basis functions Y R areused instead of the complex-valued SH basis functions. The real-valued SH basis functions are defined in terms of complex-valued SH basisfunctions as follows: Y l , m R = i √ (cid:16) Y l , m C − ( − ) m Y l , − m C (cid:17) , for m < Y l , m C , for m = √ (cid:16) Y l , − m C + ( − ) m Y l , m C (cid:17) , for m > (19) Note we use the subscript R and C do distinguish between real- and complex-valued SH basis functions respectively. L and Emission Field Q As with the complex-valued case, the angular dependent quantities are projected into SH coefficients. Here, those coefficients will be real-valued, since we use the real-valued SH basis. L l , m ( (cid:126) x ) = (cid:90) Ω L ( (cid:126) x , ω ) Y l , m R d ω . Q l , m ( (cid:126) x ) = (cid:90) Ω Q ( (cid:126) x , ω ) Y l , m R d ω . The reconstruction ˆ L and ˆ Q , is found by a truncated linear combination of SH basis functions weighted by their respective coefficients:ˆ L ( (cid:126) x , ω ) = ∑ l , m L l , m ( (cid:126) x ) Y l , m R ( ω ) . (20) ˆ Q ( (cid:126) x , ω ) = ∑ l , m Q l , m ( (cid:126) x ) Y l , m R ( ω ) . (21) We later will have to apply identities and properties for the complex-valued SH basis functions and therefore need to expand the real-valued c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media basis function in ˆ L . The real-valued basis function is different depending on m and therefore gives different expansions for the sign of m :ˆ L ( (cid:126) x , ω ) = ∑ l , m L l , m ( (cid:126) x ) i √ (cid:16) Y l , m C − ( − ) m Y l , − m C (cid:17) , for m < ∑ l , m L l , m ( (cid:126) x ) Y l , m C , for m = ∑ l , m L l , m ( (cid:126) x ) √ (cid:16) Y l , − m C − ( − ) m Y l , m C (cid:17) , for m > = ∑ l ∑ − m = − l L l , m ( (cid:126) x ) i √ (cid:16) Y l , m C − ( − ) m Y l , − m C (cid:17) , for m < L l , ( (cid:126) x ) Y l , C , for m = ∑ l ∑ lm = L l , m ( (cid:126) x ) √ (cid:16) Y l , − m C − ( − ) m Y l , m C (cid:17) , for m > = N ∑ l = (cid:32) − ∑ m = − l L l , m ( (cid:126) x ) (cid:18) i √ Y l , m C ( ω ) − i √ ( − ) m Y l , − m C ( ω ) (cid:19) + L l , ( (cid:126) x ) Y l , C ( ω )+ l ∑ m = L l , m ( (cid:126) x ) (cid:18) √ Y l , − m C ( ω ) + √ ( − ) m Y l , m C ( ω ) (cid:19)(cid:33) = i √ (cid:32) N ∑ l = − ∑ m = − l L l , m ( (cid:126) x ) Y l , m C ( ω ) (cid:33) − i √ (cid:32) N ∑ l = − ∑ m = − l L l , m ( (cid:126) x ) ( − ) m Y l , − m C ( ω ) (cid:33) + N ∑ l = L l , ( (cid:126) x ) Y l , C ( ω )+ √ (cid:32) N ∑ l = l ∑ m = L l , m ( (cid:126) x ) Y l , − m C ( ω ) (cid:33) + √ (cid:32) N ∑ l = l ∑ m = L l , m ( (cid:126) x ) ( − ) m Y l , m C ( ω ) (cid:33) . (22) The real-valued spherical harmonics expansion of the phase function follows the same derivation as the complex-valued expansion fromsection 10.1.2. We first fix the outgoing direction vector ω o to always align with the z -axis ( ω o = (cid:126) e ). We compute the spherical harmonicsprojection over incident direction vector ω i , using the real-valued spherical harmonics basis functions: p l , m = (cid:90) Ω p ( ω i · (cid:126) e ) Y l , m R ( ω i ) d ω i . Phase functions which only depend on the angle between incident and outgoing vectors are rotationally symmetric around the pole axis.Like with the complex-values spherical harmonics basis functions, such a rotation R of angle α around the pole axis is given by: ρ R ( α ) ( Y l , m R ) = e − im α Y l , m R . We formulate symmetry around the pole axis with the following constraint: ∑ l , m e − im α p l , m Y l , m R ( ω i ) = ∑ l , m p l , m Y l , m R ( ω i ) . By comparing coefficients we get p l , m = p l , m e − im α . From this we can infer that p l , m = m (cid:54) =
0, if the phase function is rotationally symmetric around the pole axis. We thereforehave the same property as we have with complex-valued expansions of functions, which are symmetric about the pole axis: only the m = p ( ω i ) = ∑ l p l Y l R ( ω i ) . (23) The transport term of the RTE is given as ( ω · ∇ ) L ( (cid:126) x , ω ) = ω x ∂ x L ( (cid:126) x , ω ) + ω y ∂ y L ( (cid:126) x , ω ) + ω z ∂ z L ( (cid:126) x , ω ) . (24) c (cid:13) (cid:13)(cid:13)
0, if the phase function is rotationally symmetric around the pole axis. We thereforehave the same property as we have with complex-valued expansions of functions, which are symmetric about the pole axis: only the m = p ( ω i ) = ∑ l p l Y l R ( ω i ) . (23) The transport term of the RTE is given as ( ω · ∇ ) L ( (cid:126) x , ω ) = ω x ∂ x L ( (cid:126) x , ω ) + ω y ∂ y L ( (cid:126) x , ω ) + ω z ∂ z L ( (cid:126) x , ω ) . (24) c (cid:13) (cid:13)(cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media To improve readability, we first project the term into SH by multiplying with the conjugate complex of the SH basis, and replace L by itsSH expansion afterwards. This order was reversed, when we derived the complex-valued P N -equation in section 10.2.1.We now multiply Equation 24 with the real-valued SH basis and integrate over solid angle. However, the SH basis is different for m (cid:48) < m (cid:48) = m (cid:48) >
0, and therefore will give us different P N -equations depending on m (cid:48) . We will go through the derivation in detail for the m (cid:48) < m (cid:48) < (cid:90) (cid:18) − i √ Y l (cid:48) , m (cid:48) ( ω ) − − i √ ( − ) m (cid:48) Y l (cid:48) , − m (cid:48) ( ω ) (cid:19) ( ω x ∂ x L ( (cid:126) x , ω ) + ω y ∂ y L ( (cid:126) x , ω ) + ω z ∂ z L ( (cid:126) x , ω )) d ω = (cid:90) − i √ Y l (cid:48) , m (cid:48) ( ω ) ω x ∂ x L ( (cid:126) x , ω ) − i √ Y l (cid:48) , m (cid:48) ( ω ) ω y ∂ y L ( (cid:126) x , ω ) − i √ Y l (cid:48) , m (cid:48) ( ω ) ω z ∂ z L ( (cid:126) x , ω )+ i √ ( − ) m (cid:48) Y l (cid:48) , − m (cid:48) ( ω ) ω x ∂ x L ( (cid:126) x , ω ) + i √ ( − ) m (cid:48) Y l (cid:48) , − m (cid:48) ( ω ) ω y ∂ y L ( (cid:126) x , ω )+ i √ ( − ) m (cid:48) Y l (cid:48) , − m (cid:48) ( ω ) ω z ∂ z L ( (cid:126) x , ω ) d ω . After expanding the integrand and splitting the integral, we apply the recursive relation from Equation 14 to get: i √ c l (cid:48) − , m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − ( ω ) d ω − i √ d l (cid:48) + , m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − ( ω ) d ω − i √ e l (cid:48) − , m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + ( ω ) d ω + i √ f l (cid:48) + , m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + ( ω ) d ω − i √ i c l (cid:48) − , m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − ( ω ) d ω + i √ i d l (cid:48) + , m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − ( ω ) d ω − i √ i e l (cid:48) − , m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + ( ω ) d ω + i √ i f l (cid:48) + , m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + ( ω ) d ω − i √ a l (cid:48) − , m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) ( ω ) d ω − i √ b l (cid:48) + , m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) ( ω ) d ω − i √ ( − ) m (cid:48) c l (cid:48) − , − m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − ( ω ) d ω + i √ ( − ) m (cid:48) d l (cid:48) + , − m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − ( ω ) d ω + i √ ( − ) m (cid:48) e l (cid:48) − , − m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + ( ω ) d ω − i √ ( − ) m (cid:48) f l (cid:48) + , − m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + ( ω ) d ω + i √ ( − ) m (cid:48) i c l (cid:48) − , − m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − ( ω ) d ω − i √ ( − ) m (cid:48) i d l (cid:48) + , − m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − ( ω ) d ω + i √ ( − ) m (cid:48) i e l (cid:48) − , − m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + ( ω ) d ω − i √ ( − ) m (cid:48) i f l (cid:48) + , − m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + ( ω ) d ω + i √ ( − ) m (cid:48) a l (cid:48) − , − m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) ( ω ) d ω + i √ ( − ) m (cid:48) b l (cid:48) + , − m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) ( ω ) d ω . Before we further expand the radiance field L into its SH expansion, we will simplify coefficients by using the following relations: a l , m = a l , − m , b l , m = b l , − m , c l , m = e l , − m , d l , m = f l , − m . (25) c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media This allows us to rewrite the equation above as: − i α c (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − ( ω ) d ω + ( − ) m (cid:48) i α c (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + ( ω ) d ω + i α d (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − ( ω ) d ω − ( − ) m (cid:48) i α d (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + ( ω ) d ω − i α e (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + ( ω ) d ω + ( − ) m (cid:48) i α e (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − ( ω ) d ω + i α f (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + ( ω ) d ω − ( − ) m (cid:48) i α f (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − ( ω ) d ω + α c (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − ( ω ) d ω + ( − ) m (cid:48) α c (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + ( ω ) d ω − α e (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + ( ω ) d ω − ( − ) m (cid:48) α e (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − ( ω ) d ω + α f (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + ( ω ) d ω + ( − ) m (cid:48) α f (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − ( ω ) d ω − α d (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − ( ω ) d ω − ( − ) m (cid:48) α d (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + ( ω ) d ω − α a (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) ( ω ) d ω + ( − ) m (cid:48) α a (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) ( ω ) d ω − α b (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) ( ω ) d ω + ( − ) m (cid:48) α b (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) ( ω ) d ω with α c = i √ c l (cid:48) − , m (cid:48) − , α e = i √ e l (cid:48) − , m (cid:48) + , α d = i √ d l (cid:48) + , m (cid:48) − α f = i √ f l (cid:48) + , m (cid:48) + , α a = i √ a l (cid:48) − , m (cid:48) , α b = i √ b l (cid:48) + , m (cid:48) . In the next step, we substitute the radiance field function L with its spherical harmonics expansion and arrive at the following expressionafter further expansions and transformations: − i α c ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + ( − ) m (cid:48) i α c ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω (26) + i α d ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) − C ( ω ) d ω − ( − ) m (cid:48) i α d ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω (27) − i α e ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) + C ( ω ) d ω + ( − ) m (cid:48) i α e ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω (28) + i α f ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) + C ( ω ) d ω − ( − ) m (cid:48) i α f ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω (29) + α c ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + ( − ) m (cid:48) α c ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω (30) − α e ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) + C ( ω ) d ω − ( − ) m (cid:48) α e ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω (31) + α f ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) + C ( ω ) d ω + ( − ) m (cid:48) α f ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω (32) − α d ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) − C ( ω ) d ω − ( − ) m (cid:48) α d ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω (33) − α a ∂ z ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) C ( ω ) d ω + ( − ) m (cid:48) α a ∂ z ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) C ( ω ) d ω (34) − α b ∂ z ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) C ( ω ) d ω + ( − ) m (cid:48) α b ∂ z ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) C ( ω ) d ω (35) The real-valued P N -equation have an intricate structure which causes many terms to cancel out. We take the first two terms (Equation 26) c (cid:13) (cid:13)(cid:13)
0, and therefore will give us different P N -equations depending on m (cid:48) . We will go through the derivation in detail for the m (cid:48) < m (cid:48) < (cid:90) (cid:18) − i √ Y l (cid:48) , m (cid:48) ( ω ) − − i √ ( − ) m (cid:48) Y l (cid:48) , − m (cid:48) ( ω ) (cid:19) ( ω x ∂ x L ( (cid:126) x , ω ) + ω y ∂ y L ( (cid:126) x , ω ) + ω z ∂ z L ( (cid:126) x , ω )) d ω = (cid:90) − i √ Y l (cid:48) , m (cid:48) ( ω ) ω x ∂ x L ( (cid:126) x , ω ) − i √ Y l (cid:48) , m (cid:48) ( ω ) ω y ∂ y L ( (cid:126) x , ω ) − i √ Y l (cid:48) , m (cid:48) ( ω ) ω z ∂ z L ( (cid:126) x , ω )+ i √ ( − ) m (cid:48) Y l (cid:48) , − m (cid:48) ( ω ) ω x ∂ x L ( (cid:126) x , ω ) + i √ ( − ) m (cid:48) Y l (cid:48) , − m (cid:48) ( ω ) ω y ∂ y L ( (cid:126) x , ω )+ i √ ( − ) m (cid:48) Y l (cid:48) , − m (cid:48) ( ω ) ω z ∂ z L ( (cid:126) x , ω ) d ω . After expanding the integrand and splitting the integral, we apply the recursive relation from Equation 14 to get: i √ c l (cid:48) − , m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − ( ω ) d ω − i √ d l (cid:48) + , m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − ( ω ) d ω − i √ e l (cid:48) − , m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + ( ω ) d ω + i √ f l (cid:48) + , m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + ( ω ) d ω − i √ i c l (cid:48) − , m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − ( ω ) d ω + i √ i d l (cid:48) + , m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − ( ω ) d ω − i √ i e l (cid:48) − , m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + ( ω ) d ω + i √ i f l (cid:48) + , m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + ( ω ) d ω − i √ a l (cid:48) − , m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) ( ω ) d ω − i √ b l (cid:48) + , m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) ( ω ) d ω − i √ ( − ) m (cid:48) c l (cid:48) − , − m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − ( ω ) d ω + i √ ( − ) m (cid:48) d l (cid:48) + , − m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − ( ω ) d ω + i √ ( − ) m (cid:48) e l (cid:48) − , − m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + ( ω ) d ω − i √ ( − ) m (cid:48) f l (cid:48) + , − m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + ( ω ) d ω + i √ ( − ) m (cid:48) i c l (cid:48) − , − m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − ( ω ) d ω − i √ ( − ) m (cid:48) i d l (cid:48) + , − m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − ( ω ) d ω + i √ ( − ) m (cid:48) i e l (cid:48) − , − m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + ( ω ) d ω − i √ ( − ) m (cid:48) i f l (cid:48) + , − m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + ( ω ) d ω + i √ ( − ) m (cid:48) a l (cid:48) − , − m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) ( ω ) d ω + i √ ( − ) m (cid:48) b l (cid:48) + , − m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) ( ω ) d ω . Before we further expand the radiance field L into its SH expansion, we will simplify coefficients by using the following relations: a l , m = a l , − m , b l , m = b l , − m , c l , m = e l , − m , d l , m = f l , − m . (25) c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media This allows us to rewrite the equation above as: − i α c (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − ( ω ) d ω + ( − ) m (cid:48) i α c (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + ( ω ) d ω + i α d (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − ( ω ) d ω − ( − ) m (cid:48) i α d (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + ( ω ) d ω − i α e (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + ( ω ) d ω + ( − ) m (cid:48) i α e (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − ( ω ) d ω + i α f (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + ( ω ) d ω − ( − ) m (cid:48) i α f (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − ( ω ) d ω + α c (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − ( ω ) d ω + ( − ) m (cid:48) α c (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + ( ω ) d ω − α e (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + ( ω ) d ω − ( − ) m (cid:48) α e (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − ( ω ) d ω + α f (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + ( ω ) d ω + ( − ) m (cid:48) α f (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − ( ω ) d ω − α d (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − ( ω ) d ω − ( − ) m (cid:48) α d (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + ( ω ) d ω − α a (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) ( ω ) d ω + ( − ) m (cid:48) α a (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) ( ω ) d ω − α b (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) ( ω ) d ω + ( − ) m (cid:48) α b (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) ( ω ) d ω with α c = i √ c l (cid:48) − , m (cid:48) − , α e = i √ e l (cid:48) − , m (cid:48) + , α d = i √ d l (cid:48) + , m (cid:48) − α f = i √ f l (cid:48) + , m (cid:48) + , α a = i √ a l (cid:48) − , m (cid:48) , α b = i √ b l (cid:48) + , m (cid:48) . In the next step, we substitute the radiance field function L with its spherical harmonics expansion and arrive at the following expressionafter further expansions and transformations: − i α c ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + ( − ) m (cid:48) i α c ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω (26) + i α d ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) − C ( ω ) d ω − ( − ) m (cid:48) i α d ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω (27) − i α e ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) + C ( ω ) d ω + ( − ) m (cid:48) i α e ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω (28) + i α f ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) + C ( ω ) d ω − ( − ) m (cid:48) i α f ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω (29) + α c ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + ( − ) m (cid:48) α c ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω (30) − α e ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) + C ( ω ) d ω − ( − ) m (cid:48) α e ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω (31) + α f ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) + C ( ω ) d ω + ( − ) m (cid:48) α f ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω (32) − α d ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) − C ( ω ) d ω − ( − ) m (cid:48) α d ∂ x ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω (33) − α a ∂ z ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) C ( ω ) d ω + ( − ) m (cid:48) α a ∂ z ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) C ( ω ) d ω (34) − α b ∂ z ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) C ( ω ) d ω + ( − ) m (cid:48) α b ∂ z ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) C ( ω ) d ω (35) The real-valued P N -equation have an intricate structure which causes many terms to cancel out. We take the first two terms (Equation 26) c (cid:13) (cid:13)(cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media of the P N -equations and apply the following orthogonality property of SH: (cid:90) Ω Y l , m R Y l , m C d ω = i √ (cid:16) δ l = l m = m − ( − ) m δ l = l m = − m (cid:17) , for m < δ l = l m = m , for m = √ (cid:16) δ l = l m = − m + ( − ) m δ l = l m = m (cid:17) , for m > . (36) This way we get for the first two terms: − i α c ∂ y N ∑ l = − ∑ m = − l L l , m ( (cid:126) x ) √ i δ l = l (cid:48) − m = m (cid:48) − + i α c ∂ y N ∑ l = − ∑ m = − l L l , m ( (cid:126) x ) √ i ( − ) m δ l = l (cid:48) − m = − m (cid:48) + − i α c ∂ y N ∑ l = L l , ( (cid:126) x ) √ δ l = l (cid:48) − = m (cid:48) − − i α c ∂ y N ∑ l = l ∑ m = L l , m ( (cid:126) x ) √ δ l = l (cid:48) − m = − m (cid:48) + − i α c ∂ y N ∑ l = l ∑ m = L l , m ( (cid:126) x ) √ ( − ) m δ l = l (cid:48) − m = m (cid:48) − + ( − ) m (cid:48) i α c ∂ y N ∑ l = − ∑ m = − l L l , m ( (cid:126) x ) √ i δ l = l (cid:48) − m = − m (cid:48) + − ( − ) m (cid:48) i α c ∂ y N ∑ l = − ∑ m = − l L l , m ( (cid:126) x ) √ i ( − ) m δ l = l (cid:48) − m = m (cid:48) − + ( − ) m (cid:48) i α c ∂ y N ∑ l = L l , ( (cid:126) x ) √ δ l = l (cid:48) − = − m (cid:48) + + ( − ) m (cid:48) i α c ∂ y N ∑ l = l ∑ m = L l , m ( (cid:126) x ) √ δ l = l (cid:48) − m = m (cid:48) − + ( − ) m (cid:48) i α c ∂ y N ∑ l = l ∑ m = L l , m ( (cid:126) x ) √ ( − ) m δ l = l (cid:48) − m = − m (cid:48) + . We apply the delta function for the sums which run over the variable l : N ∑ l = b ∑ m = a L l , m δ l = xm = y = b ∑ m = a L x , m δ m = y . (37) We get for the first two terms of the transport term of the P N -equation (Equation 26): − i α c ∂ y − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) √ i δ m = m (cid:48) − + i α c ∂ y − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) √ i ( − ) m δ m = − m (cid:48) + − i α c ∂ y L l (cid:48) − , ( (cid:126) x ) √ δ m (cid:48) = − i α c ∂ y l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) √ δ m = − m (cid:48) + − i α c ∂ y l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) √ ( − ) m δ m = m (cid:48) − + ( − ) m (cid:48) i α c ∂ y − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) √ i δ m = − m (cid:48) + − ( − ) m (cid:48) i α c ∂ y − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) √ i ( − ) m δ m = m (cid:48) − + ( − ) m (cid:48) i α c ∂ y L l (cid:48) − , ( (cid:126) x ) √ δ m (cid:48) = + ( − ) m (cid:48) i α c ∂ y l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) √ δ m = m (cid:48) − + ( − ) m (cid:48) i α c ∂ y l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) √ ( − ) m δ m = − m (cid:48) + The variables l (cid:48) and m (cid:48) specify a particular equation within the given set of P N -equations. We remember that m (cid:48) originated from multiplyingthe transport term with the real-valued SH basis function Y R for the projection. The real-valued basis function is different for the sign of m (cid:48) and we derived the transport term of the P N -equations under the assumption of m (cid:48) < m (cid:48) = m (cid:48) > m (cid:48) and that m (cid:48) < m , up to − m > m < m (cid:48) < δ m = m (cid:48) − or δ m = − m (cid:48) + , an even m is selected if m (cid:48) is odd and viceversa. Therefore, we have ( − ) m ( − ) m (cid:48) = −
1. This causes term one and seven (red) to vanish and term four and ten (black) to collapse intoone term. c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media Therefore, the first two terms in the expansion (Equation 26), simplify to: − i α c ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + ( − ) m (cid:48) i α c ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω = − √ i α c ∂ y L l (cid:48) − , − m (cid:48) + ( (cid:126) x )= − √ i i √ c l (cid:48) − , m (cid:48) − ∂ y L l (cid:48) − , − m (cid:48) + ( (cid:126) x )= c l (cid:48) − , m (cid:48) − ∂ y L l (cid:48) − , − m (cid:48) + ( (cid:126) x ) . The terms in Equation 27 are derived in the same way with the difference, that the signs are reversed and that we have l (cid:48) + l (cid:48) −
1. However, this does not affect the simplification: i α d ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) − C ( ω ) d ω − ( − ) m (cid:48) i α c ∂ y ∑ l , m L l , m ( (cid:126) x ) (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω = − √ i α d ∂ y L l (cid:48) − , − m (cid:48) + ( (cid:126) x )= − √ i i √ d l (cid:48) + , m (cid:48) − ∂ y L l (cid:48) + , − m (cid:48) + ( (cid:126) x )= d l (cid:48) + , m (cid:48) − ∂ y L l (cid:48) + , − m (cid:48) + ( (cid:126) x ) . Carrying out the same simplifications for the remaining terms, results in the following real-valued P N -equations for m (cid:48) < − c l (cid:48) − , m (cid:48) − ∂ y L l (cid:48) − , − m (cid:48) + + d l (cid:48) + , m (cid:48) − ∂ y L l (cid:48) + , − m (cid:48) + − β m (cid:48) e l (cid:48) − , m (cid:48) + ∂ y L l (cid:48) − , − m (cid:48) − + β m (cid:48) f l (cid:48) + , m (cid:48) + ∂ y L l (cid:48) + , − m (cid:48) − + c l (cid:48) − , m (cid:48) − ∂ x L l (cid:48) − , m (cid:48) − − δ m (cid:48) (cid:54) = − e l (cid:48) − , m (cid:48) + ∂ x L l (cid:48) − , m (cid:48) + + δ m (cid:48) (cid:54) = − f l (cid:48) + , m (cid:48) + ∂ x L l (cid:48) + , m (cid:48) + − d l (cid:48) + , m (cid:48) − ∂ x L l (cid:48) + , m (cid:48) − + a l (cid:48) − , m (cid:48) ∂ z L l (cid:48) − , m (cid:48) + b l (cid:48) + , m (cid:48) ∂ z L l (cid:48) + , m (cid:48) . with β x = (cid:40) √ , for | x | = , for | x | (cid:54) = . (38) We now carry out the same derivation for the assumption of m (cid:48) >
0. We multiply Equation 10.2.1 with the definition of the real-valued SHbasis for m (cid:48) > (cid:90) (cid:18) √ Y l (cid:48) , − m (cid:48) C ( ω ) + √ ( − ) m (cid:48) Y l (cid:48) , m (cid:48) C ( ω ) (cid:19) ( ω x ∂ x L ( (cid:126) x , ω ) + ω y ∂ y L ( (cid:126) x , ω ) + ω z ∂ z L ( (cid:126) x , ω )) d ω c (cid:13) (cid:13)(cid:13)
0. We multiply Equation 10.2.1 with the definition of the real-valued SHbasis for m (cid:48) > (cid:90) (cid:18) √ Y l (cid:48) , − m (cid:48) C ( ω ) + √ ( − ) m (cid:48) Y l (cid:48) , m (cid:48) C ( ω ) (cid:19) ( ω x ∂ x L ( (cid:126) x , ω ) + ω y ∂ y L ( (cid:126) x , ω ) + ω z ∂ z L ( (cid:126) x , ω )) d ω c (cid:13) (cid:13)(cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media We expand the integrand and split the integral. Then we apply the recursive relation from Equation 14 and get:1 √ c l (cid:48) − , − m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω − √ d l (cid:48) + , − m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω − √ e l (cid:48) − , − m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω + √ f l (cid:48) + , − m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω − √ i c l (cid:48) − , − m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω + √ i d l (cid:48) + , − m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω − √ i e l (cid:48) − , − m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω + √ i f l (cid:48) + , − m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω + √ a l (cid:48) − , − m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) C ( ω ) d ω + √ b l (cid:48) + , − m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) C ( ω ) d ω + √ ( − ) m (cid:48) c l (cid:48) − , m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − C ( ω ) d ω − √ ( − ) m (cid:48) d l (cid:48) + , m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − C ( ω ) d ω − √ ( − ) m (cid:48) e l (cid:48) − , m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + C ( ω ) d ω + √ ( − ) m (cid:48) f l (cid:48) + , m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + C ( ω ) d ω − √ ( − ) m (cid:48) i c l (cid:48) − , m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + √ ( − ) m (cid:48) i d l (cid:48) + , m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − C ( ω ) d ω − √ ( − ) m (cid:48) i e l (cid:48) − , m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + C ( ω ) d ω + √ ( − ) m (cid:48) i f l (cid:48) + , m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + C ( ω ) d ω + √ ( − ) m (cid:48) a l (cid:48) − , m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) C ( ω ) d ω + √ ( − ) m (cid:48) b l (cid:48) + , m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) C ( ω ) d ω We simplify these using the identities from Equation 25: α c (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω − ( − ) m (cid:48) α c (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + C ( ω ) d ω − α d (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω + ( − ) m (cid:48) α d (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + C ( ω ) d ω − α e (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω + ( − ) m (cid:48) α e (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + α f (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω − ( − ) m (cid:48) α f (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − C ( ω ) d ω − i α c (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω − ( − ) m (cid:48) i α c (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + C ( ω ) d ω + i α d (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω + ( − ) m (cid:48) i α d (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + C ( ω ) d ω − i α e (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω − ( − ) m (cid:48) i α e (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + i α f (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω + ( − ) m (cid:48) i α f (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − C ( ω ) d ω + α a (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) C ( ω ) d ω + ( − ) m (cid:48) α a (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) C ( ω ) d ω + α b (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) C ( ω ) d ω + ( − ) m (cid:48) α b (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) C ( ω ) d ω c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media with α c = √ c l (cid:48) − , − m (cid:48) − , α e = √ e l (cid:48) − , − m (cid:48) + , α d = √ d l (cid:48) + , − m (cid:48) − α f = √ f l (cid:48) + , − m (cid:48) + , α a = √ a l (cid:48) − , − m (cid:48) , α b = √ b l (cid:48) + , − m (cid:48) . We substitute the radiance field function L with its spherical harmonics expansion and arrive at the following expression after furtherexpansions and transformations: α c ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω − ( − ) m (cid:48) α c ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) + C ( ω ) d ω − α d ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω + ( − ) m (cid:48) α d ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) + C ( ω ) d ω − α e ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω + ( − ) m (cid:48) α e ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + α f ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω − ( − ) m (cid:48) α f ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) − C ( ω ) d ω − i α c ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω − ( − ) m (cid:48) i α c ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) + C ( ω ) d ω + i α d ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω + ( − ) m (cid:48) i α d ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) + C ( ω ) d ω − i α e ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω − ( − ) m (cid:48) i α e ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + i α f ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω + ( − ) m (cid:48) i α f ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) − C ( ω ) d ω + α a ∂ z ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) C ( ω ) d ω + ( − ) m (cid:48) α a ∂ z ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) C ( ω ) d ω + α b ∂ z ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) C ( ω ) d ω + ( − ) m (cid:48) α b ∂ z ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) C ( ω ) d ω Again we apply the identity given in Equation 36. For the first two terms we for example get: α c ∂ x − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) i √ δ m = − m (cid:48) − − α c ∂ x − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) i √ ( − ) m δ m = m (cid:48) + + α c ∂ x L l (cid:48) − , ( (cid:126) x ) √ δ − m (cid:48) = + α c ∂ x l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) √ δ m = m (cid:48) + + α c ∂ x l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) √ ( − ) m δ m = − m (cid:48) − − ( − ) − m (cid:48) α c ∂ x − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) i √ δ m = m (cid:48) + + ( − ) − m (cid:48) α c ∂ x − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) i √ ( − ) m δ m = − m (cid:48) − − ( − ) − m (cid:48) α c ∂ x L l (cid:48) − , ( (cid:126) x ) √ δ − m (cid:48) = − ( − ) − m (cid:48) α c ∂ x l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) √ δ m = − m (cid:48) − − ( − ) − m (cid:48) α c ∂ x l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) √ ( − ) m δ m = m (cid:48) + As with the m (cid:48) < c (cid:13) (cid:13)(cid:13)
0. We multiply Equation 10.2.1 with the definition of the real-valued SHbasis for m (cid:48) > (cid:90) (cid:18) √ Y l (cid:48) , − m (cid:48) C ( ω ) + √ ( − ) m (cid:48) Y l (cid:48) , m (cid:48) C ( ω ) (cid:19) ( ω x ∂ x L ( (cid:126) x , ω ) + ω y ∂ y L ( (cid:126) x , ω ) + ω z ∂ z L ( (cid:126) x , ω )) d ω c (cid:13) (cid:13)(cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media We expand the integrand and split the integral. Then we apply the recursive relation from Equation 14 and get:1 √ c l (cid:48) − , − m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω − √ d l (cid:48) + , − m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω − √ e l (cid:48) − , − m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω + √ f l (cid:48) + , − m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω − √ i c l (cid:48) − , − m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω + √ i d l (cid:48) + , − m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω − √ i e l (cid:48) − , − m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω + √ i f l (cid:48) + , − m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω + √ a l (cid:48) − , − m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) C ( ω ) d ω + √ b l (cid:48) + , − m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) C ( ω ) d ω + √ ( − ) m (cid:48) c l (cid:48) − , m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − C ( ω ) d ω − √ ( − ) m (cid:48) d l (cid:48) + , m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − C ( ω ) d ω − √ ( − ) m (cid:48) e l (cid:48) − , m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + C ( ω ) d ω + √ ( − ) m (cid:48) f l (cid:48) + , m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + C ( ω ) d ω − √ ( − ) m (cid:48) i c l (cid:48) − , m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + √ ( − ) m (cid:48) i d l (cid:48) + , m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − C ( ω ) d ω − √ ( − ) m (cid:48) i e l (cid:48) − , m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + C ( ω ) d ω + √ ( − ) m (cid:48) i f l (cid:48) + , m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + C ( ω ) d ω + √ ( − ) m (cid:48) a l (cid:48) − , m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) C ( ω ) d ω + √ ( − ) m (cid:48) b l (cid:48) + , m (cid:48) (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) C ( ω ) d ω We simplify these using the identities from Equation 25: α c (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω − ( − ) m (cid:48) α c (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + C ( ω ) d ω − α d (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω + ( − ) m (cid:48) α d (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + C ( ω ) d ω − α e (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω + ( − ) m (cid:48) α e (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + α f (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω − ( − ) m (cid:48) α f (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − C ( ω ) d ω − i α c (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω − ( − ) m (cid:48) i α c (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + C ( ω ) d ω + i α d (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω + ( − ) m (cid:48) i α d (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + C ( ω ) d ω − i α e (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω − ( − ) m (cid:48) i α e (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + i α f (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω + ( − ) m (cid:48) i α f (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − C ( ω ) d ω + α a (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , − m (cid:48) C ( ω ) d ω + ( − ) m (cid:48) α a (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) C ( ω ) d ω + α b (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , − m (cid:48) C ( ω ) d ω + ( − ) m (cid:48) α b (cid:90) ∂ z L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) C ( ω ) d ω c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media with α c = √ c l (cid:48) − , − m (cid:48) − , α e = √ e l (cid:48) − , − m (cid:48) + , α d = √ d l (cid:48) + , − m (cid:48) − α f = √ f l (cid:48) + , − m (cid:48) + , α a = √ a l (cid:48) − , − m (cid:48) , α b = √ b l (cid:48) + , − m (cid:48) . We substitute the radiance field function L with its spherical harmonics expansion and arrive at the following expression after furtherexpansions and transformations: α c ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω − ( − ) m (cid:48) α c ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) + C ( ω ) d ω − α d ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω + ( − ) m (cid:48) α d ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) + C ( ω ) d ω − α e ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω + ( − ) m (cid:48) α e ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + α f ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω − ( − ) m (cid:48) α f ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) − C ( ω ) d ω − i α c ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) − C ( ω ) d ω − ( − ) m (cid:48) i α c ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) + C ( ω ) d ω + i α d ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) − C ( ω ) d ω + ( − ) m (cid:48) i α d ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) + C ( ω ) d ω − i α e ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) + C ( ω ) d ω − ( − ) m (cid:48) i α e ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) − C ( ω ) d ω + i α f ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) + C ( ω ) d ω + ( − ) m (cid:48) i α f ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) − C ( ω ) d ω + α a ∂ z ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , − m (cid:48) C ( ω ) d ω + ( − ) m (cid:48) α a ∂ z ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) C ( ω ) d ω + α b ∂ z ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , − m (cid:48) C ( ω ) d ω + ( − ) m (cid:48) α b ∂ z ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) C ( ω ) d ω Again we apply the identity given in Equation 36. For the first two terms we for example get: α c ∂ x − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) i √ δ m = − m (cid:48) − − α c ∂ x − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) i √ ( − ) m δ m = m (cid:48) + + α c ∂ x L l (cid:48) − , ( (cid:126) x ) √ δ − m (cid:48) = + α c ∂ x l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) √ δ m = m (cid:48) + + α c ∂ x l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) √ ( − ) m δ m = − m (cid:48) − − ( − ) − m (cid:48) α c ∂ x − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) i √ δ m = m (cid:48) + + ( − ) − m (cid:48) α c ∂ x − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) i √ ( − ) m δ m = − m (cid:48) − − ( − ) − m (cid:48) α c ∂ x L l (cid:48) − , ( (cid:126) x ) √ δ − m (cid:48) = − ( − ) − m (cid:48) α c ∂ x l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) √ δ m = − m (cid:48) − − ( − ) − m (cid:48) α c ∂ x l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) √ ( − ) m δ m = m (cid:48) + As with the m (cid:48) < c (cid:13) (cid:13)(cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media P N -equations for the transport term therefore are: α c ∂ x l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) √ δ m = m (cid:48) + − ( − ) − m (cid:48) α c ∂ x l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) √ ( − ) m δ m = m (cid:48) + = √ α c ∂ x L l (cid:48) − , m (cid:48) + ( (cid:126) x ) = √ (cid:18) √ c l (cid:48) − , − m (cid:48) − (cid:19) ∂ x L l (cid:48) − , m (cid:48) + ( (cid:126) x )= c l (cid:48) − , − m (cid:48) − L l (cid:48) − , m (cid:48) + ( (cid:126) x ) . Following this through for the remaining terms gives us the real-valued P N -equations for m (cid:48) > c l (cid:48) − , − m (cid:48) − ∂ x L l (cid:48) − , m (cid:48) + ( (cid:126) x ) − d l (cid:48) + , − m (cid:48) − ∂ x L l (cid:48) + , m (cid:48) + ( (cid:126) x ) − β m (cid:48) e l (cid:48) − , m (cid:48) − ∂ x L l (cid:48) − , m (cid:48) − ( (cid:126) x ) β m (cid:48) f l (cid:48) + , − m (cid:48) + ∂ x L l (cid:48) + , m (cid:48) − ( (cid:126) x ) c l (cid:48) − , − m (cid:48) − ∂ y L l (cid:48) − , − m (cid:48) − ( (cid:126) x ) − d l (cid:48) + , − m (cid:48) − ∂ y L l (cid:48) + , − m (cid:48) − ( (cid:126) x ) δ m (cid:48) (cid:54) = e l (cid:48) − , − m (cid:48) + ∂ y L l (cid:48) − , − m (cid:48) + ( (cid:126) x ) − δ m (cid:48) (cid:54) = f l (cid:48) + , − m (cid:48) + ∂ y L l (cid:48) + , − m (cid:48) + ( (cid:126) x ) a l (cid:48) − , − m (cid:48) ∂ z L l (cid:48) − , m (cid:48) ( (cid:126) x ) b l (cid:48) + , − m (cid:48) ∂ z L l (cid:48) + , m (cid:48) ( (cid:126) x ) Finally the m (cid:48) = P N -equations as in this case, thereal-valued SH basis function is identical to the complex-valued SH basis function. We multiply Equation 10.2.1 with the definition of thereal-valued SH basis for m (cid:48) = (cid:90) Y l (cid:48) , m (cid:48) C ( ω ) ( ω x ∂ x L ( (cid:126) x , ω ) + ω y ∂ y L ( (cid:126) x , ω ) + ω z ∂ z L ( (cid:126) x , ω )) d ω Expanding the integrand and applying the recursion relation (Equation 25) produces the following set of terms:12 c l (cid:48) − , m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − C ( ω ) d ω − e l (cid:48) − , m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + C ( ω ) d ω − d l (cid:48) + , m (cid:48) − (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − C ( ω ) d ω + f l (cid:48) + , m (cid:48) + (cid:90) ∂ x L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + C ( ω ) d ω − i c l (cid:48) − , m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) − C ( ω ) d ω − i e l (cid:48) − , m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) − , m (cid:48) + C ( ω ) d ω + i d l (cid:48) + , m (cid:48) − (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) − C ( ω ) d ω + i f l (cid:48) + , m (cid:48) + (cid:90) ∂ y L ( (cid:126) x , ω ) Y l (cid:48) + , m (cid:48) + C ( ω ) d ω + a l (cid:48) − , m (cid:48) (cid:90) Y l (cid:48) − , m (cid:48) C ( ω ) ∂ z L ( (cid:126) x , ω ) d ω + b l (cid:48) + , m (cid:48) (cid:90) Y l (cid:48) + , m (cid:48) C ( ω ) ∂ z L ( (cid:126) x , ω ) d ω Again we will replace the radiance field L with its real-valued SH projection and get:12 c l (cid:48) − , m (cid:48) − ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) − C ( ω ) d ω − e l (cid:48) − , m (cid:48) + ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) + C ( ω ) d ω − d l (cid:48) + , m (cid:48) − ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) − C ( ω ) d ω + f l (cid:48) + , m (cid:48) + ∂ x ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) + C ( ω ) d ω − i c l (cid:48) − , m (cid:48) − ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) − C ( ω ) d ω − i e l (cid:48) − , m (cid:48) + ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) − , m (cid:48) + C ( ω ) d ω + i d l (cid:48) + , m (cid:48) − ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) − C ( ω ) d ω + i f l (cid:48) + , m (cid:48) + ∂ y ∑ l , m L l , m (cid:90) Y l , m R Y l (cid:48) + , m (cid:48) + C ( ω ) d ω + a l (cid:48) − , m (cid:48) ∂ z ∑ l , m L l , m (cid:90) Y l (cid:48) − , m (cid:48) C ( ω ) Y l , m R d ω + b l (cid:48) + , m (cid:48) ∂ z ∑ l , m L l , m (cid:90) Y l (cid:48) + , m (cid:48) C ( ω ) Y l , m R d ω These terms also have an intricate structure where many terms cancel out and simplify. This is seen once we apply the SH orthogonality c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media property (Equation 36) and further consider that m (cid:48) =
0. We show this for the first two terms, which expand to:12 c l (cid:48) − , − i √ ∂ x (cid:32) − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) δ − , m (cid:33) − c l (cid:48) − , − i √ ∂ x (cid:32) − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) ( − ) m δ − , − m (cid:33) c l (cid:48) − , − ∂ x (cid:16) L l (cid:48) − , ( (cid:126) x ) δ − , (cid:17) c l (cid:48) − , − √ ∂ x (cid:32) l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) δ − , − m (cid:33) c l (cid:48) − , − √ ∂ x (cid:32) l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) ( − ) m δ − , m (cid:33) − e l (cid:48) − , i √ ∂ x (cid:32) − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) δ , m (cid:33) e l (cid:48) − , i √ ∂ x (cid:32) − ∑ m = − l (cid:48) + L l (cid:48) − , m ( (cid:126) x ) ( − ) m δ , − m (cid:33) − e l (cid:48) − , ∂ x (cid:16) L l (cid:48) − , ( (cid:126) x ) δ , (cid:17) − e l (cid:48) − , √ ∂ x (cid:32) l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) δ , − m (cid:33) − e l (cid:48) − , √ ∂ x (cid:32) l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) ( − ) m δ , m (cid:33) Again the blue terms vanish since the delta functions will never be non-zero under the sums. The red terms cancel each other out since c l , − = e l , and − m = − m = −
1. The terms in black simplify to:12 c l (cid:48) − , − √ ∂ x (cid:32) l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) δ − , − m (cid:33) − e l (cid:48) − , √ ∂ x (cid:32) l (cid:48) − ∑ m = L l (cid:48) − , m ( (cid:126) x ) ( − ) m δ , m (cid:33) = c l (cid:48) − , − √ ∂ x L l (cid:48) − , ( (cid:126) x ) − c l (cid:48) − , − √ ∂ x L l (cid:48) − , ( (cid:126) x ) ( − ) = √ c l (cid:48) − , − ∂ x L l (cid:48) − , ( (cid:126) x ) . Similar simplifications apply to the remaining terms of the SH expansion of the transport term for m =
0, resulting in the final expression:1 √ c l (cid:48) − , − ∂ x L l (cid:48) − , ( (cid:126) x ) − √ d l (cid:48) + , − ∂ x L l (cid:48) + , ( (cid:126) x ) √ c l (cid:48) − , − ∂ y L l (cid:48) − , − ( (cid:126) x ) − √ d l (cid:48) + , − ∂ y L l (cid:48) + , − ( (cid:126) x ) a l (cid:48) − , ∂ z L l (cid:48) − , ( (cid:126) x ) + b l (cid:48) + , ∂ z L l (cid:48) + , ( (cid:126) x ) The collision term of the RTE is given as: − σ t ( (cid:126) x ) L ( (cid:126) x , ω ) We first replace the radiance field L with its real-valued SH expansion: − σ t ( (cid:126) x ) ∑ l , m L l , m ( (cid:126) x ) Y l , m R ( ω ) In order to project the term into SH, we have to multiply with the real-valued SH basis function and integrate over solid angle. Since thebasis function is different depending on m (cid:48) < m = m >
0, we have to derive separate P N -equations for each case.We first derive the SH projection of the collision term for the case m (cid:48) <
0. Multiplying with the SH basis and integrating over solid angle c (cid:13) (cid:13)(cid:13)
0. Multiplying with the SH basis and integrating over solid angle c (cid:13) (cid:13)(cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media gives, after some further transformations and application of the SH orthogonality property: − i √ σ t ( (cid:126) x ) i √ − ∑ m = − l (cid:48) L l (cid:48) , m ( (cid:126) x ) δ m (cid:48) , m + i √ σ t ( (cid:126) x ) i √ − ∑ m = − l (cid:48) L l (cid:48) , m ( (cid:126) x ) ( − ) m δ m (cid:48) , − m − i √ σ t ( (cid:126) x ) L l (cid:48) , ( (cid:126) x ) δ m (cid:48) , − i √ σ t ( (cid:126) x ) √ l (cid:48) ∑ m = L l (cid:48) , m ( (cid:126) x ) δ m (cid:48) , − m − i √ σ t ( (cid:126) x ) √ l (cid:48) ∑ m = L l (cid:48) , m ( (cid:126) x ) ( − ) m δ m (cid:48) , m + i √ ( − ) m (cid:48) σ t ( (cid:126) x ) i √ − ∑ m = − l (cid:48) L l (cid:48) , m ( (cid:126) x ) δ − m (cid:48) , m − i √ ( − ) m (cid:48) σ t ( (cid:126) x ) i √ − ∑ m = − l (cid:48) L l (cid:48) , m ( (cid:126) x ) ( − ) m δ − m (cid:48) , − m + i √ ( − ) m (cid:48) σ t ( (cid:126) x ) L l (cid:48) , ( (cid:126) x ) δ − m (cid:48) , + i √ ( − ) m (cid:48) σ t ( (cid:126) x ) √ l (cid:48) ∑ m = L l (cid:48) , m ( (cid:126) x ) δ − m (cid:48) , − m + i √ ( − ) m (cid:48) σ t ( (cid:126) x ) √ l (cid:48) ∑ m = L l (cid:48) , m ( (cid:126) x ) ( − ) m δ − m (cid:48) , m As for the transport term derivation, the blue terms vanish due to the delta function being always zero under the sum. The red terms canceleach other out. The remaining term (black) determines the SH projection of the collision term for m (cid:48) < σ t L l (cid:48) , m (cid:48) (39) The derivation for the SH projection of the collision term for m > σ t L l (cid:48) , m (cid:48) (40) The real-values SH projection of the collision term for m = σ t L l (cid:48) , m (cid:48) (41) The scattering term is given as a convolution of the radiance field L with the phase function p using a rotation R ω : σ s ( (cid:126) x ) (cid:90) Ω (cid:48) p ( (cid:126) x , ω (cid:48) · ω ) L ( (cid:126) x , ω (cid:48) ) d ω (cid:48) = σ s ( (cid:126) x )( L ◦ ρ R ( ω ) ( p ))( ω ) where the convolution can be also expressed as a inner product integral: ( L ◦ ρ R ( ω ) ( p )) = (cid:90) Ω (cid:48) L ( (cid:126) x , ω (cid:48) ) ρ R ( ω ) ( p )( ω (cid:48) ) d ω (cid:48) = (cid:104) L , ρ R ( ω ) ( p ) (cid:105) . We substitute L with its real-valued SH-expansion (Equation 20) in the inner product integral and perform some further factorizations toget: i √ N ∑ l = N ∑ l (cid:48) − ∑ m = − l L l , m ( (cid:126) x ) f l (cid:48) (cid:68) Y l , m C ( ω ) , ρ R ( ω ) (cid:16) Y l (cid:48) C (cid:17)(cid:69) − i √ N ∑ l = N ∑ l (cid:48) − ∑ m = − l ( − ) m L l , m ( (cid:126) x ) f l (cid:48) (cid:68) Y l , − m C ( ω ) , ρ R ( ω ) (cid:16) Y l (cid:48) C (cid:17)(cid:69) + N ∑ l = N ∑ l (cid:48) L l , ( (cid:126) x ) f l (cid:48) (cid:68) Y l , C ( ω ) , ρ R ( ω ) (cid:16) Y l (cid:48) C (cid:17)(cid:69) + √ N ∑ l = N ∑ l (cid:48) l ∑ m = L l , m ( (cid:126) x ) f l (cid:48) (cid:68) Y l , − m C ( ω ) , ρ R ( ω ) (cid:16) Y l (cid:48) C (cid:17)(cid:69) + √ N ∑ l = N ∑ l (cid:48) l ∑ m = ( − ) m L l , m ( (cid:126) x ) f l (cid:48) (cid:68) Y l , m C ( ω ) , ρ R ( ω ) (cid:16) Y l (cid:48) C (cid:17)(cid:69) The spherical harmonics basis functions Y lm C are orthogonal. We therefore have (cid:68) Y lm , ρ R ( ω ) (cid:16) Y l (cid:48) m (cid:48) (cid:17)(cid:69) =
0, for all l (cid:54) = l (cid:48) , which further c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media simplifies our scattering operator to i √ N ∑ l = − ∑ m = − l L l , m ( (cid:126) x ) f l (cid:68) Y l , m C ( ω ) , ρ R ( ω ) (cid:16) Y l C (cid:17)(cid:69) (42) − i √ N ∑ l = − ∑ m = − l ( − ) m L l , m ( (cid:126) x ) f l (cid:68) Y l , − m C ( ω ) , ρ R ( ω ) (cid:16) Y l (cid:17)(cid:69) (43) + N ∑ l = L l , ( (cid:126) x ) f l (cid:68) Y l , C ( ω ) , ρ R ( ω ) (cid:16) Y l C (cid:17)(cid:69) (44) + √ N ∑ l = l ∑ m = L l , m ( (cid:126) x ) f l (cid:68) Y l , − m C ( ω ) , ρ R ( ω ) (cid:16) Y l C (cid:17)(cid:69) (45) + √ N ∑ l = l ∑ m = ( − ) m L l , m ( (cid:126) x ) f l (cid:68) Y l , m C ( ω ) , ρ R ( ω ) (cid:16) Y l C (cid:17)(cid:69) (46) What remains to be resolved are the inner products. We use the fact that the spherical harmonics basis functions Y lm C are eigenfunctions ofthe inner product integral operator in the equation above, i.e. (cid:68) Y lm C , ρ R ( ω ) (cid:16) Y l C (cid:17)(cid:69) = λ l Y lm C (47) which results in: i √ N ∑ l = − ∑ m = − l L l , m ( (cid:126) x ) f l λ l Y l , m C ( ω ) (48) − i √ N ∑ l = − ∑ m = − l ( − ) m L l , m ( (cid:126) x ) f l λ l Y l , − m C ( ω ) (49) + N ∑ l = L l , ( (cid:126) x ) f l λ l Y l , C ( ω ) (50) + √ N ∑ l = l ∑ m = L l , m ( (cid:126) x ) f l λ l Y l , − m C ( ω ) (51) + √ N ∑ l = l ∑ m = ( − ) m L l , m ( (cid:126) x ) f l λ l Y l , m C ( ω ) . (52) The next step is to project the scattering term into real-valued SH. Again we will have to use different terms for m < m = m > m < i √ σ s ( (cid:126) x ) i √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) − ∑ m = − l (cid:48) L l (cid:48) , m ( (cid:126) x ) δ m (cid:48) , m − i √ σ s ( (cid:126) x ) i √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) − ∑ m = − l (cid:48) ( − ) m L l (cid:48) , m ( (cid:126) x ) δ m (cid:48) , − m + i √ σ s ( (cid:126) x ) L l (cid:48) , ( (cid:126) x ) p l (cid:48) , ( (cid:126) x ) λ l (cid:48) δ m (cid:48) , + i √ σ s ( (cid:126) x ) √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) l (cid:48) ∑ m = L l (cid:48) , m ( (cid:126) x ) δ m (cid:48) , − m + i √ σ s ( (cid:126) x ) √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) l (cid:48) ∑ m = ( − ) m L l (cid:48) , m ( (cid:126) x ) δ m (cid:48) , m − i √ ( − ) m (cid:48) σ s ( (cid:126) x ) i √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) − ∑ m = − l (cid:48) L l (cid:48) , m ( (cid:126) x ) δ − m (cid:48) , m + i √ ( − ) m (cid:48) σ s ( (cid:126) x ) i √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) − ∑ m = − l (cid:48) ( − ) m L l (cid:48) , m ( (cid:126) x ) δ − m (cid:48) , − m − i √ ( − ) m (cid:48) σ s ( (cid:126) x ) L l (cid:48) , ( (cid:126) x ) p l (cid:48) , ( (cid:126) x ) λ l (cid:48) δ − m (cid:48) , − i √ ( − ) m (cid:48) σ s ( (cid:126) x ) √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) l (cid:48) ∑ m = L l (cid:48) , m ( (cid:126) x ) δ − m (cid:48) , − m − i √ ( − ) m (cid:48) σ s ( (cid:126) x ) √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) l (cid:48) ∑ m = ( − ) m L l (cid:48) , m ( (cid:126) x ) δ − m (cid:48) , m Again, the blue terms vanish, because the delta functions will always be zero under the sum. The red terms cancel each other out. The blackterms reduce to: − σ s ( (cid:126) x ) λ l (cid:48) p l (cid:48) , ( (cid:126) x ) L l (cid:48) , m (cid:48) ( (cid:126) x ) The same happens for the derivation for m > m = c (cid:13) (cid:13)(cid:13)
0, for all l (cid:54) = l (cid:48) , which further c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media simplifies our scattering operator to i √ N ∑ l = − ∑ m = − l L l , m ( (cid:126) x ) f l (cid:68) Y l , m C ( ω ) , ρ R ( ω ) (cid:16) Y l C (cid:17)(cid:69) (42) − i √ N ∑ l = − ∑ m = − l ( − ) m L l , m ( (cid:126) x ) f l (cid:68) Y l , − m C ( ω ) , ρ R ( ω ) (cid:16) Y l (cid:17)(cid:69) (43) + N ∑ l = L l , ( (cid:126) x ) f l (cid:68) Y l , C ( ω ) , ρ R ( ω ) (cid:16) Y l C (cid:17)(cid:69) (44) + √ N ∑ l = l ∑ m = L l , m ( (cid:126) x ) f l (cid:68) Y l , − m C ( ω ) , ρ R ( ω ) (cid:16) Y l C (cid:17)(cid:69) (45) + √ N ∑ l = l ∑ m = ( − ) m L l , m ( (cid:126) x ) f l (cid:68) Y l , m C ( ω ) , ρ R ( ω ) (cid:16) Y l C (cid:17)(cid:69) (46) What remains to be resolved are the inner products. We use the fact that the spherical harmonics basis functions Y lm C are eigenfunctions ofthe inner product integral operator in the equation above, i.e. (cid:68) Y lm C , ρ R ( ω ) (cid:16) Y l C (cid:17)(cid:69) = λ l Y lm C (47) which results in: i √ N ∑ l = − ∑ m = − l L l , m ( (cid:126) x ) f l λ l Y l , m C ( ω ) (48) − i √ N ∑ l = − ∑ m = − l ( − ) m L l , m ( (cid:126) x ) f l λ l Y l , − m C ( ω ) (49) + N ∑ l = L l , ( (cid:126) x ) f l λ l Y l , C ( ω ) (50) + √ N ∑ l = l ∑ m = L l , m ( (cid:126) x ) f l λ l Y l , − m C ( ω ) (51) + √ N ∑ l = l ∑ m = ( − ) m L l , m ( (cid:126) x ) f l λ l Y l , m C ( ω ) . (52) The next step is to project the scattering term into real-valued SH. Again we will have to use different terms for m < m = m > m < i √ σ s ( (cid:126) x ) i √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) − ∑ m = − l (cid:48) L l (cid:48) , m ( (cid:126) x ) δ m (cid:48) , m − i √ σ s ( (cid:126) x ) i √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) − ∑ m = − l (cid:48) ( − ) m L l (cid:48) , m ( (cid:126) x ) δ m (cid:48) , − m + i √ σ s ( (cid:126) x ) L l (cid:48) , ( (cid:126) x ) p l (cid:48) , ( (cid:126) x ) λ l (cid:48) δ m (cid:48) , + i √ σ s ( (cid:126) x ) √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) l (cid:48) ∑ m = L l (cid:48) , m ( (cid:126) x ) δ m (cid:48) , − m + i √ σ s ( (cid:126) x ) √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) l (cid:48) ∑ m = ( − ) m L l (cid:48) , m ( (cid:126) x ) δ m (cid:48) , m − i √ ( − ) m (cid:48) σ s ( (cid:126) x ) i √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) − ∑ m = − l (cid:48) L l (cid:48) , m ( (cid:126) x ) δ − m (cid:48) , m + i √ ( − ) m (cid:48) σ s ( (cid:126) x ) i √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) − ∑ m = − l (cid:48) ( − ) m L l (cid:48) , m ( (cid:126) x ) δ − m (cid:48) , − m − i √ ( − ) m (cid:48) σ s ( (cid:126) x ) L l (cid:48) , ( (cid:126) x ) p l (cid:48) , ( (cid:126) x ) λ l (cid:48) δ − m (cid:48) , − i √ ( − ) m (cid:48) σ s ( (cid:126) x ) √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) l (cid:48) ∑ m = L l (cid:48) , m ( (cid:126) x ) δ − m (cid:48) , − m − i √ ( − ) m (cid:48) σ s ( (cid:126) x ) √ p l (cid:48) , ( (cid:126) x ) λ l (cid:48) l (cid:48) ∑ m = ( − ) m L l (cid:48) , m ( (cid:126) x ) δ − m (cid:48) , m Again, the blue terms vanish, because the delta functions will always be zero under the sum. The red terms cancel each other out. The blackterms reduce to: − σ s ( (cid:126) x ) λ l (cid:48) p l (cid:48) , ( (cid:126) x ) L l (cid:48) , m (cid:48) ( (cid:126) x ) The same happens for the derivation for m > m = c (cid:13) (cid:13)(cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media The derivation of the real-valued SH projection of the emission term is exactly the same as for the collision and scattering term. Afterreplacing the emission term Q with its real-valued SH expansion, we multiply by the real-valued SH basis function for m < − i √ i √ − ∑ m = − l (cid:48) Q l (cid:48) , m ( (cid:126) x ) δ m (cid:48) , m + i √ i √ − ∑ m = − l (cid:48) Q l (cid:48) , m ( (cid:126) x ) ( − ) m δ m (cid:48) , − m − i √ δ m (cid:48) , Q l (cid:48) , ( (cid:126) x ) − i √ √ l (cid:48) ∑ m = Q l (cid:48) , m ( (cid:126) x ) δ m (cid:48) , − m − i √ √ l (cid:48) ∑ m = Q l (cid:48) , m ( (cid:126) x ) ( − ) m δ m (cid:48) , m + i √ ( − ) m (cid:48) i √ − ∑ m = − l (cid:48) Q l (cid:48) , m ( (cid:126) x ) δ − m (cid:48) , m − i √ ( − ) m (cid:48) i √ − ∑ m = − l (cid:48) Q l (cid:48) , m ( (cid:126) x ) ( − ) m δ − m (cid:48) , − m + i √ ( − ) m (cid:48) δ − m (cid:48) , Q l (cid:48) , ( (cid:126) x ) + i √ ( − ) m (cid:48) √ l (cid:48) ∑ m = Q l (cid:48) , m ( (cid:126) x ) δ − m (cid:48) , − m + i √ ( − ) m (cid:48) √ l (cid:48) ∑ m = Q l (cid:48) , m ( (cid:126) x ) ( − ) m δ − m (cid:48) , m Again, the blue terms vanish and the red terms cancel each other out. The black terms collapses to: Q l (cid:48) , m (cid:48) for m < m = m > Putting all projected terms from previous subsections together, we get for m = √ c l − , − ∂ x L l − , − √ d l + , − ∂ x L l + , √ c l − , − ∂ y L l − , − − √ d l + , − ∂ y L l + , − a l − , ∂ z L l − , + b l + , ∂ z L l + , + σ t L l , m − σ s λ l p l , L l , m = Q l , m . (53) For m < − c l − , m − ∂ y L l − , − m + + d l + , m − ∂ y L l + , − m + − β m e l − , m + ∂ y L l − , − m − + β m f l + , m + ∂ y L l + , − m − + c l − , m − ∂ x L l − , m − − δ m (cid:54) = − e l − , m + ∂ x L l − , m + + δ m (cid:54) = − f l + , m + ∂ x L l + , m + − d l + , m − ∂ x L l + , m − + a l − , m ∂ z L l − , m + b l + , m ∂ z L l + , m + σ t L l , m − σ s λ l p l , L l , m = Q l , m . (54) And for m >
0: 12 c l − , − m − ∂ x − d l + , − m − ∂ x L l + , m + − β m e l − , m − ∂ x L l − , m − + β m f l + , − m + ∂ x L l + , m − + c l − , − m − ∂ y L l − , − m − − d l + , − m − ∂ y L l + , − m − + δ m (cid:54) = e l − , − m + ∂ y L l − , − m + − δ m (cid:54) = f l + , − m + ∂ y L l + , − m + + a l − , − m ∂ z L l − , m + b l + , − m ∂ z L l + , m + σ t L l , m − σ s λ l p l , L l , m = Q l , m . (55) Where we defined: β x = (cid:40) √ , for | x | = , for | x | (cid:54) = , δ x (cid:54) = y = (cid:26) , for x (cid:54) = y , for x = y . c (cid:13) (cid:13) . Koerner & J. Portsmouth & W. Jakob / P N -Method for Multiple Scattering in Participating Media and a l , m = (cid:115) ( l − m + ) ( l + m + )( l + ) ( l − ) , b l , m = (cid:115) ( l − m ) ( l + m )( l + ) ( l − ) c l , m = (cid:115) ( l + m + ) ( l + m + )( l + ) ( l + ) , d l , m = (cid:115) ( l − m ) ( l − m − )( l + ) ( l − ) e l , m = (cid:115) ( l − m + ) ( l − m + )( l + ) ( l + ) , f l , m = (cid:115) ( l + m ) ( l + m − )( l + ) ( l − ) λ l = (cid:114) π l + . These equations can be written in a more compact form by using ± and ∓ to write the equations for m < m > c (cid:13) (cid:13)(cid:13)