PN-Method for Multiple Scattering in Participating Media

P _N -Method for Multiple Scattering in Participating Media

Supplemental Material

David Koerner

Jamie Portsmouth

Wenzel Jakob

1University of Stuttgart ²Solid Angle ³´Ecole polytechnique f´ed´erale de Lausanne (EPFL)

This document contains the detailed derivations of the complex-valued and real-valued PN- equations, referred to by the article for better readability.

Contents

1 Isotropic Radiative Transfer Equation 2

2 Derivation of the complex-valued P_N-equations 2

2.1 Projecting Radiative Transfer Quantities . . . 3

2.1.1 Radiance FieldLand Emission FieldQ . . . 3

2.1.2 Phase Function . . . 3

2.2 Projecting Terms of the RTE . . . 4

2.2.1 Transport Term . . . 4

2.2.2 Collision Term . . . 5

2.2.3 Scattering Term . . . 6

2.2.4 Emission Term. . . 8

2.3 Final Equation . . . 8

3 Derivation of the real-valued P_N-equations 8 3.1 Projecting Radiative Transfer Quantities . . . 8

3.1.1 Radiance FieldLand Emission FieldQ . . . 8

3.1.2 Phase Function . . . 9

3.2 Projecting Terms of the RTE . . . 10

3.2.1 Transport Term . . . 10

3.2.2 Collision Term . . . 21

3.2.3 Scattering Term . . . 22

3.2.4 Emission Term. . . 24

3.3 Final equation . . . 25

1 Isotropic Radiative Transfer Equation

ThePN-equations are derived from the radiative transfer equation (RTE), which expresses the change of the radiance fieldL, with respect to an infinitesimal change of position into directionω at point~x:

(∇·ω)L(~x, ω) =−σt(~x)L(~x, ω) +σ_s(~x)

Z

Ω

p(ω⁰·ω)L(~x, ω) dω⁰ +Q(~x, ω)

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

The derivation is done in two steps: First, the directional-dependent quantities are replaced by their SH- projected counterparts. 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ω. This is why in a second step, the RTE is projected into spherical harmonics, which is done by multiplying each term with the conjugate complex of the SH basis functions.

The SH basis functions are complex, which produces complex coefficients and complexPN-equations. How- ever, there are also real SH basis functions, which are defined in terms of the complex SH basis functions and which produce real coefficients and reconstructions. Since the radiance fieldLis real, the real SH basis functions is used in our article.

In the next section, the complexP_N-equations are derived. In order to give the derivation a better structure, the two steps mentioned above are applied to each term individually in a seperate subsection. The section concludes by putting all derived terms together. In Section3we derive the realP_N-equations which are used in the article.

2 Derivation of the complex-valued P

-equations

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

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

Y^l,m

C (ω) =Y^l,m

C (θ, φ) =







(−1)^mq

2l+1 4π

(l−m)!

(l+m)!e^imφP^l,m(cos (θ)), form≥0 (−1)^mY^l|m|

C (θ, φ), form <0

(1)

whereP^l,mare the associated Legendre-Polynomials. The(−1)^mfactor is called the Condon-Shortley phase and is not part of the associated Legendre Polynomial like in other definitions.

2.1 Projecting Radiative Transfer Quantities

2.1.1 Radiance FieldL and Emission FieldQ

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

L^l,m(~x) = Z

Ω

L(~x, ω)Y^l,m

C dω Q^l,m(~x) =

Z

Ω

Q(~x, ω)Y^l,m

C dω

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

L(~x, ω)≈Lˆ(~x, ω) =

N

X

l=0 l

X

m=−l

L^l,m(~x)Y^l,m

C (ω) =X

l,m

L^l,m(~x)Y^l,m

C (ω) Q(~x, ω)≈Qˆ(~x, ω) =

N

X

l=0 l

X

m=−l

Q^l,m(~x)Y^l,m

C (ω) =X

l,m

Q^l,m(~x)Y^l,m

C (ω) (2)

2.1.2 Phase Function

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

p^l,m= Z

Ω

p(ω_i·~e₃)Y^l,m

C (ω_i) dω_i

The expansion of the phase function can be further simplified, because the phase function is rotationally sym- metric around the pole axis~e3. Consider the definition of the spherical harmonics basis functionY_C^l,m:

Y^l,m

C (θ, φ) =C^l,me^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_C^l,m) =e^−imαY_C^l,m

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

ρ_R(α)(p) =p and in spherical harmonics this would be:

X

l,m

e^−imαp^l,mY^l,m

C (ω_i) =X

l,m

p^l,mY^l,m

C (ω_i) By equating coefficients we get:

p^l,m=p^l,me^−imα

Since e^−imα = 1 for all αonly when m = 0, we can conclude that p^l,m = 0 for all m 6= 0. This means that for a function, which is rotationally symmetric around the pole axis, only them= 0coefficients will be valid. Therefore, our phase function reconstruction for a fixed outgoing vector (ωo =~e3) only requires SH coefficients withm= 0:

p(ωi) =X

l

p^l0Y_C^l0(ωi) (3)

2.2 Projecting Terms of the RTE

2.2.1 Transport Term

The transport term of the RTE is given as

(ω·∇)L(~x, ω) ReplacingLwith its expansion gives:

(ω·∇)



 X

l,m

L^l,m(~x)Y^l,m

C (ω)





Next we multiply withY^l0m0

C and integrate over solid angle:

Z

Ω

Y^l0m0(ω·∇)X

l,m

L^l,m(~x)Y^l,m

C (ω) dω We can pull the spatial derivative out of the integral to get:

∇·Z

Ω

ωY^l0m0

C

X

l,m

L^l,m(~x)Y^l,m

C (ω) dω (4)

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

ωY^l,m

C = 1 2









c^l−1,m−1Y^l−1,m−1

C −d^l+1,m−1Y^l+1,m−1

C −e^l−1,m+1Y^l−1,m+1

C +f^l+1,m+1Y^l+1,m+1

C

i

−c^l−1,m−1Y^l−1,m−1

C +d^l+1,m−1Y^l+1,m−1

C −e^l−1,m+1Y^l−1,m+1

C +f^l+1,m+1Y^l+1,m+1

C

2

a^l−1,mY^l−1,m

C +b^l+1,mY^l+1,m

C







 (5)

with a^l,m=

s

(l−m+ 1) (l+m+ 1)

(2l+ 1) (2l−1) b^l,m= s

(l−m) (l+m)

(2l+ 1) (2l−1) c^l,m= s

(l+m+ 1) (l+m+ 2) (2l+ 3) (2l+ 1) d^l,m=

s

(l−m) (l−m−1)

(2l+ 1) (2l−1) e^l,m= s

(l−m+ 1) (l−m+ 2)

(2l+ 3) (2l+ 1) f^l,m= s

(l+m) (l+m−1) (2l+ 1) (2l−1) Note that the signs for thex- andy- component depend on the handedness of the coordinate system in which the SH basis functions are defined. Using this in equation4gives





1 2∂x i 2∂y

∂z



·Z

Ω







c^{l0−1,m0−1}Y^{l0−1,m0−1}

C −d^l0+1,m0−1Y^l0+1,m0−1

C −e^l0−1,m0+1Y^l0−1,m0+1

C +f^l0+1,m0+1Y^l0+1,m0+1

C

−c^{l0−1,m0−1}Y^{l0−1,m0−1}

C +d^l0+1,m0−1Y^l0+1,m0−1

C −e^l0−1,m0+1Y^l0−1,m0+1

C +f^l0+1,m0+1Y^l0+1,m0+1

C

a^l0−1,m0Y^l0−1,m0

C +b^l0+1,m0Y^l0+1,m0

C





 X

l,m

L^l,m(~x)Y^l,m

C (ω) dω

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





1 2∂x i 2∂_y

∂_z



·









c^{l0−1,m0−1}P

l,mL^l,m(~x)R

ΩY^{l0−1,m0−1}

C (ω)Y^l,m

C (ω) dω − ...

−c^{l0−1,m0−1}P

l,mL^l,m(~x)R

ΩY^{l0−1,m0−1}

C (ω)Y^l,m

C (ω) dω + ...

a^l0−1,m0P

l,mL^l,m(~x)R

ΩY^l0−1,m0

C (ω)Y^l,m

C (ω) dω + ...









Applying the orthogonality property to the solid angle integrals will will select specificl, min each term:





1 2∂_x

i 2∂_y

∂_z



·





c^l−1,m−1L^l−1,m−1−d^l+1,m−1L^l+1,m−1−e^l−1,m+1L^l−1,m+1+f^l+1,m+1L^l+1,m+1

−c^l−1,m−1L^l−1,m−1+d^l+1,m−1L^l+1,m−1−e^l−1,m+1L^l−1,m+1+f^l+1,m+1L^l+1,m+1 a^l−1,mL^l−1,m+b^l+1,mL^l+1,m





Which gives the final moment equation for the transport term:

=1

2∂_x c^l−1,m−1L^l−1,m−1−d^l+1,m−1L^l+1,m−1−e^l−1,m+1L^l−1,m+1+f^l+1,m+1L^l+1,m+1 + i

2∂y −c^l−1,m−1L^l−1,m−1+d^l+1,m−1L^l+1,m−1−e^l−1,m+1L^l−1,m+1+f^l+1,m+1L^l+1,m+1 +

∂_z a^l−1,mL^l−1,m+b^l+1,mL^l+1,m

=1

2c^l−1,m−1∂xL^l−1,m−1−1

2d^l+1,m−1∂xL^l+1,m−1−1

2e^l−1,m+1∂xL^l−1,m+1+1

2f^l+1,m+1∂xL^l+1,m+1+

− i

2c^l−1,m−1∂_yL^l−1,m−1+i

2d^l+1,m−1∂_yL^l+1,m−1− i

2e^l−1,m+1∂_yL^l−1,m+1+ i

2f^l+1,m+1∂_yL^l+1,m+1+ a^l−1,m∂zL^l−1,m+b^l+1,m∂zL^l+1,m

2.2.2 Collision Term

The collision term of the RTE is given as:

−σ_t(~x)L(~x, ω)

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

−σt(~x)X

l,m

L^l,m(~x)Y_C^l,m(ω)

Multiplying with Y^l0m0

C and integrating over solid angle gives after pulling some factors out of the inte- gral:

−σ_t(~x)X

l,m

L^l,m(~x) Z

Ω

Y^l0m0

C (ω)Y^l,m

C (ω) dω

=−σt(~x)X

l,m

L^l,m(~x)δll⁰δmm⁰

=−σt(~x)L^l,m(~x) 2.2.3 Scattering Term

The scattering term in the RTE is given as:

σs(~x) Z

Ω

p(~x, ω⁰·ω)L(~x, ω⁰) dω⁰

The phase function used in isotropic scattering medium only depends on the angle between incident and outgoing direction and therefore is rotationally symmetric around the pole defining axis. This property allows us to define a rotationR(ω), which rotates the phase function, such 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 ofp. With this rotated phase function, we now can express the integral of the scattering operator as a convolution denoted with the symbol◦:

Z

Ω

p(~x, ω⁰·ω)L(~x, ω⁰) dω⁰=L◦ρ_R(ω)(p)

= Z

Ω⁰

L(~x, ω⁰)ρ_R(ω)(p)(ω⁰) dω⁰

=hL, ρR(ω)(p)i (6) As we evaluate the inner product integral of the convolution, the phase function rotates along with the argumentω.

We now use the spherical harmonics expansions ofL (equation (2)) andp(equation (14)) in the definition for the inner product of our convolution (equation (6)):

hL, ρR(ω)(p)i=

* X

l,m

L^l,m(~x)Y_C^l,m, ρR(ω)

X

l

p^l0Y_C^l0

!+

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

hL, ρ_R(ω)(p)i=X

l,m

L^l,m(~x)

* Y^l,m

C , ρ_R(ω) X

l

p^l0Y_C^l0

!+

and further:

hL, ρR(ω)(p)i=X

l⁰

X

l,m

p^l00L^l,m(~x)D Y^l,m

C , ρ_R(ω) Y_C^l00E

The rotation ρ_R(ω) of a function with frequency l gives a function of frequency l again. In addition the spherical harmonics basis functionsY^l,m

C are orthogonal. We therefore have:

D Y^l,m

C , ρ_R(ω)

Y_C^l0m0E

= 0 for all l6=l⁰ which further simplifies our inner product integral to:

hL, ρ_R(ω)(p)i=X

l,m

p^l0L^l,m(~x)D Y^l,m

C , ρ_R(ω) Y_C^l0E

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 of the inner product integral operator in the equation above:

DY^l,m

C , ρ_R(ω) Y_C^l0E

=λ_lY^l,m

C

with

λ_l= r 4π

2l+ 1 Replacing the inner product gives:

hL, ρR(ω)(p)i=X

l,m

λlp^l0L^l,m(~x)Y_C^l,m

This allows us to express the scattering term using SH expansions of phase function p and radiance field L:

σs(~x) Z

Ω

p(~x, ω⁰·ω)L(~x, ω⁰) dω⁰=σs(~x)hL, ρR(ω)(p)i

=σs(~x)X

l,m

λlp^l0L^l,m(~x)Y^l,m

C

However, we haven’t done a spherical harmonics expansion of the term itsself. It is still a scalar function which depends on directionω. We project the scattering term into spherical harmonics, by multiplying with Y^l0m0

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 the complex-valuedPN-equations:

Z

Ω

Y_C^l0m0(ω)σs(~x)X

l,m

λlp^l0L^l,m(~x)Y_C^l,m(ω) dω

=λ_lσ_s(~x)p^l0L^l,m(~x)X

l,m

Z

Ω

Y^l0m0

C (ω)Y^l,m

C (ω) dω

=λlσs(~x)p^l0L^l,m(~x)X

l,m

δll⁰δmm⁰

=λ_lσ_s(~x)p^l0L^l,m(~x) (7)

2.2.4 Emission Term

The emission term of the RTE is given as:

Q(~x, ω) (8)

The derivation of the SH projected term is equal to the derivation of the projected collision term. After replacing the emission field with its SH projection and multiplying the term with the conjugate complex of Y_Cresults after applying the orthogonality property in:

Q^l,m(~x, ω) (9)

2.3 Final Equation

We arrive at the complex-valuedP_N-equations after putting all the projected terms together:

1

2c^l−1,m−1∂xL^l−1,m−1−1

2d^l+1,m−1∂xL^l+1,m−1−1

2e^l−1,m+1∂xL^l−1,m+1+1

2f^l+1,m+1∂xL^l+1,m+1+

−i

2c^l−1,m−1∂_yL^l−1,m−1+ i

2d^l+1,m−1∂_yL^l+1,m−1− i

2e^l−1,m+1∂_yL^l−1,m+1+ i

2f^l+1,m+1∂_yL^l+1,m+1+ a^l−1,m∂zL^l−1,m+b^l+1,m∂zL^l+1,m=−σt(~x)L^l,m(~x) +λlσs(~x)p^l0L^l,m(~x) +Q^l,m(~x, ω)

3 Derivation of the real-valued P

-equations

The real-valuedP_N-equations are derived similar to their complex-valued counterpart. The key difference is, that the real-valued SH basis functions Y_R are used, instead of the the complex-valued SH basis functions.

The real-valued SH basis functions are defined in terms of complex-valued SH basis functions:

Y^l,m

R =











√i 2

Y_C^l,m−(−1)^mY_C^l,−m

, form <0

Y_C^l,m, form= 0

√1 2

Y^l,−m

C −(−1)^mY^l,m

C

, form >0

(10)

We use the subscriptR andCdo distinguish between real- and complex-valued SH basis functions respec- tively.

3.1 Projecting Radiative Transfer Quantities

3.1.1 Radiance FieldL and Emission FieldQ

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(~x) = Z

Ω

L(~x, ω)Y_R^l,mdω Q^l,m(~x) =

Z

Ω

Q(~x, ω)Y_R^l,mdω

The reconstruction Lˆ andQˆ, is found by a truncated linear combination of SH basis functions weighted by their respective coefficients:

Lˆ(~x, ω) =X

l,m

L^l,m(~x)Y^l,m

R (ω) (11)

Qˆ(~x, ω) =X

l,m

Q^l,m(~x)Y^l,m

R (ω) (12)

We later will have to apply identities and properties for the complex-valued SH basis functions and therefore need to expand the real-valued basis function in Lˆ. The real-valued basis function is different depending on mand therefore gives different expansions for the sign ofm:

Lˆ(~x, ω) =









 P

l,mL^l,m(~x)^√i

2

Y^l,m

C −(−1)^mY^l,−m

C

, form <0 P

l,mL^l,m(~x)Y_C^l,m, form= 0 P

l,mL^l,m(~x)^√1

2

Y^l,−m

C −(−1)^mY^l,m

C

, form >0

=









 P

l

P−1

m=−lL^l,m(~x)^√i

2

Y^l,m

C −(−1)^mY^l,−m

C

, form <0 L^l,0(~x)Y^l,0

C , form= 0

P

l

Pl

m=1L^l,m(~x)^√1

2

Y^l,−m

C −(−1)^mY^l,m

C

, form >0

=

N

X

l=0

−1

X

m=−l

L^l,m(~x) i

√2Y^l,m

C (ω)− i

√2(−1)^mY^l,−m

C (ω)

+L^l,0(~x)Y^l,0

C (ω) +

l

X

m=1

L^l,m(~x) 1

√2Y_C^l,−m(ω) + 1

√2(−1)^mY_C^l,m(ω) !

= i

√2

N

X

l=0

−1

X

m=−l

L^l,m(~x)Y^l,m

C (ω)

!

− i

√2

N

X

l=0

−1

X

m=−l

L^l,m(~x) (−1)^mY^l,−m

C (ω)

!

+

N

X

l=0

L^l,0(~x)Y^l,0

C (ω) + 1

√2

N

X

l=0 l

X

m=1

L^l,m(~x)Y^l,−m

C (ω)

! + 1

√2

N

X

l=0 l

X

m=1

L^l,m(~x) (−1)^mY^l,m

C (ω)

!

(13)

3.1.2 Phase Function

The real-valued spherical harmonics expansion of the phase function follows the the same derivation as the complex-valued expansion from section 2.1.2. We first fix the outgoing direction vector ωo to always align with thez-axis (ωo=~e3). We compute the spherical harmonics projection over incident direction vectorωi, using the real-valued spherical harmonics basis functions:

p^l,m= Z

Ω

p(ωi·~e3)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 rotationR 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:

X

l,m

e^−imαp^l,mY^l,m

R (ω_i) =X

l,m

p^l,mY^l,m

R (ω_i) By comparing coefficients we get

p^l,m=p^l,me^−imα

From this we can infer thatp^l,m = 0 for allm 6= 0, if the phase function is rotationally symmetric around the pole axis. We therefore have the same property as we have with complex-valued expansions of functions, which are symmetric about the pole axis: Only the m = 0 coefficients are needed for reconstruction. We therefore have

p(ω_i) =X

l

p^l0Y_R^l0(ω_i) (14)

3.2 Projecting Terms of the RTE

3.2.1 Transport Term

The transport term of the RTE is given as

(ω·∇)L(~x, ω) =ωx∂xL(~x, ω) +ωy∂yL(~x, ω) +ωz∂zL(~x, ω) (15) To improve readability, we first project the term into SH by multiplying with the conjugate complex of the Sh basis, and replace L by its Sh expansion afterwards. This order was reversed, when we derived the complex-valuedPN-equation in section2.2.1.

We now multiply equation15with the real-valued SH basis and integrate over solid angle. However, the SH basis is different form⁰<0,m⁰ = 0andm⁰>0, and therefore will give us differentPN-equations depending onm⁰. We will go through the derivation in detail for them⁰<0case and give the results for the other cases at the end.

Multiplying the expanded transport term with the SH basis for m⁰ < 0 and integrating over solid angle gives:

Z −i

√2Y^l0,m0(ω)− −i

√2(−1)^m0Y^l0,−m0(ω)

(ω_x∂_xL(~x, ω) +ω_y∂_yL(~x, ω) +ω_z∂_zL(~x, ω)) dω

= Z

− i

√2Y^l0,m0(ω)ωx∂xL(~x, ω)− i

√2Y^l0,m0(ω)ωy∂yL(~x, ω)− i

√2Y^l0,m0(ω)ωz∂zL(~x, ω) + i

√

2(−1)^m0Y^l0,−m0(ω)ωx∂xL(~x, ω) + i

√

2(−1)^m0Y^l0,−m0(ω)ωy∂yL(~x, ω) + i

√2(−1)^m0Y^l0,−m0(ω)ω_z∂_zL(~x, ω) dω

After expanding the integrand and splitting the integral, we apply the recursive relation from equation5 to get:

√i 2

1

2c^{l0−1,m0−1} Z

∂_xL(~x, ω)Y^{l0−1,m0−1}(ω) dω− i

√2 1

2d^l0+1,m0−1 Z

∂_xL(~x, ω)Y^l0+1,m0−1(ω) dω

− i

√2 1

2e^l0−1,m0+1 Z

∂_xL(~x, ω)Y^l0−1,m0+1(ω) dω+ i

√2 1

2f^l0+1,m0+1 Z

∂_xL(~x, ω)Y^l0+1,m0+1(ω) dω

− i

√ 2

i

2c^{l0−1,m0−1} Z

∂yL(~x, ω)Y^{l0−1,m0−1}(ω) dω+ i

√ 2

i

2d^l0+1,m0−1 Z

∂yL(~x, ω)Y^l0+1,m0−1(ω) dω

− i

√2 i

2e^l0−1,m0+1 Z

∂yL(~x, ω)Y^l0−1,m0+1(ω) dω+ i

√2 i

2f^l0+1,m0+1 Z

∂yL(~x, ω)Y^l0+1,m0+1(ω) dω

− i

√2a^l0−1,m0 Z

∂zL(~x, ω)Y^l0−1,m0(ω) dω− i

√2b^l0+1,m0 Z

∂zL(~x, ω)Y^l0+1,m0(ω) dω

− i

√2(−1)^m0 1

2c^{l0−1,−m0−1} Z

∂_xL(~x, ω)Y^{l0−1,−m0−1}(ω) dω + i

√2(−1)^m0 1

2d^{l0+1,−m0−1} Z

∂_xL(~x, ω)Y^{l0+1,−m0−1}(ω) dω + i

√

2(−1)^m0 1

2e^{l0−1,−m0+1} Z

∂xL(~x, ω)Y^{l0−1,−m0+1}(ω) dω

− i

√2(−1)^m0 1

2f^l0+1,−m0+1 Z

∂xL(~x, ω)Y^l0+1,−m0+1(ω) dω + i

√2(−1)^m0 i

2c^{l0−1,−m0−1} Z

∂yL(~x, ω)Y^{l0−1,−m0−1}(ω) dω

− i

√2(−1)^m0 i

2d^{l0+1,−m0−1} Z

∂_yL(~x, ω)Y^{l0+1,−m0−1}(ω) dω + i

√2(−1)^m0 i

2e^{l0−1,−m0+1} Z

∂yL(~x, ω)Y^{l0−1,−m0+1}(ω) dω

− i

√

2(−1)^m0 i

2f^l0+1,−m0+1 Z

∂yL(~x, ω)Y^l0+1,−m0+1(ω) dω + i

√2(−1)^m0a^l0−1,−m0 Z

∂zL(~x, ω)Y^l0−1,−m0(ω) dω+ i

√2(−1)^m0b^l0+1,−m0 Z

∂zL(~x, ω)Y^l0+1,−m0(ω) dω Before we further expand the radiance fieldLinto 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 (16)

This allows us to rewrite the equation above as:

−iα_c Z

∂_yL(~x, ω)Y^{l0−1,m0−1}(ω) dω+ (−1)^m0iα_c Z

∂_yL(~x, ω)Y^{l0−1,−m0+1}(ω) dω +iα_d

Z

∂_yL(~x, ω)Y^l0+1,m0−1(ω) dω−(−1)^m0iα_d Z

∂_yL(~x, ω)Y^l0+1,−m0+1(ω) dω

−iαe

Z

∂yL(~x, ω)Y^l0−1,m0+1(ω) dω+ (−1)^m0iαe

Z

∂yL(~x, ω)Y^{l0−1,−m0−1}(ω) dω +iαf

Z

∂yL(~x, ω)Y^l0+1,m0+1(ω) dω−(−1)^m0iαf

Z

∂yL(~x, ω)Y^{l0+1,−m0−1}(ω) dω +αc

Z

∂xL(~x, ω)Y^{l0−1,m0−1}(ω) dω+ (−1)^m0αc

Z

∂xL(~x, ω)Y^{l0−1,−m0+1}(ω) dω

−α_e Z

∂_xL(~x, ω)Y^l0−1,m0+1(ω) dω−(−1)^m0α_e Z

∂_xL(~x, ω)Y^{l0−1,−m0−1}(ω) dω +α_f

Z

∂_xL(~x, ω)Y^l0+1,m0+1(ω) dω+ (−1)^m0α_f Z

∂_xL(~x, ω)Y^{l0+1,−m0−1}(ω) dω

−αd

Z

∂xL(~x, ω)Y^l0+1,m0−1(ω) dω−(−1)^m0αd

Z

∂xL(~x, ω)Y^l0+1,−m0+1(ω) dω

−αa

Z

∂zL(~x, ω)Y^l0−1,m0(ω) dω+ (−1)^m0αa

Z

∂zL(~x, ω)Y^l0−1,−m0(ω) dω

−αb

Z

∂zL(~x, ω)Y^l0+1,m0(ω) dω+ (−1)^m0αb

Z

∂zL(~x, ω)Y^l0+1,−m0(ω) dω

with

αc= i

√2 1

2c^{l0−1,m0−1}, αe= i

√2 1

2e^l0−1,m0+1, αd= i

√2 1

2d^l0+1,m0−1 αf = i

√2 1

2f^l0+1,m0+1, αa= i

√2a^l0−1,m0, αb= i

√2b^l0+1,m0

In the next step, we substitute the radiance field functionLwith its spherical harmonics expansion and arrive