Implementing the microscopic self-shadowing microflake model

(1)

EUROGRAPHICS 2018 / D. Gutierrez and A. Sheffer (Guest Editors)

(2018),

Implementing the microscopic self-shadowing microflake model

Supplemental material of

"A new microflake model with microscopic self-shadowing for accurate volume downsampling"

Guillaume Loubet

INRIA, Univ. Grenoble Alpes/LJK, CNRS/LJK guillaume.loubet@gmx.fr

In the main paper, we have introduced a new participating medium model that generalizes the standard microflake model [JAM^∗10]. It has new parameters that characterize self-shadowing effects at the microscopic scale. This model is more complex than the standard model, and its implementation is difficult unless normal distributions and self-shadowing functions are carefully chosen, so that closed-form expressions and sampling procedures can be derived. In this document, we present all the details, proofs and derivations for implementing our self- shadowing model. We first recall the general expressions of this model in Sec.2and of the simplified version in which self-shadowing is isotropic (Sec.2.5). Then, we discuss the implementation of the model in Sec.4, and the implementation of the simplified model based on SGGX distribution [HDCD15] in Sec.5. We briefly recall the theory of rejection sampling in Sec.3because we use it several times in this document. We provide code that implements our model along with this document. We tested our sampling procedures using chi-square tests, and verified that our phase functions are normalized using numerical integration.

Contents 1 Notations

2 The microscopic self-shadowing model 2.1 Anisotropic RTE

2.2 Helmholtz’s law of reciprocity 2.3 Self-shadowing function 2.4 Self-shadowing model

2.5 Simplified model with isotropic self-shadowing 3 Rejection sampling

4 Implementing the self-shadowing model with trigonometric lobes 4.1 Choice of a self-shadowing function

4.2 Trigonometric lobes:Dcos(ω,ξ,n)andDsin(ω,ξ,n) 4.3 Normalization factors

4.4 Proof that distributionsDsin(ω,ξ,n)can be written as a sum of distributionsDcos(ω,ξ,n) 4.5 Closed-form expression of the attenuation coefficientσt(ω)

4.6 Single scattering coeffσssfor cosⁿand sinⁿmicroflake distributions 4.7 Closed-form expression for^Rσtand^Rσss

4.8 Sampling the distribution of visible normals for a cosine lobeDcos(m,ξ_D,n) 4.9 Sampling the distribution of visible normals for a cosine lobeD_sin(m,ξ_D,n)

c

2018 The Author(s)

(2)

Guillaume Loubet / Implementing the self-shadowing model

4.10 Sampling single scattering phase functionfssforDcos(m,ξ_D,n) 4.11 Sampling single scattering phase functionfssforDsin(m,ξD,n) 4.12 Sampling the multiple scattering phase function fms

4.13 Sampling a normal fromDcos

4.14 Sampling a normal fromDsin

4.15 Sampling a cosine lobe

4.16 Case of isotropic distributionD(ω) =_4π¹ 4.17 Combining trigonometric lobes

5 Implementing the simplified self-shadowing model with the SGGX distribution 6 Evaluation of the estimation of multiple scattering albedos

References

1. Notations

In this document, we use often use spherical coordinates and we callφ∈(0,2π)the azimuth angle andθ∈(0,π)the polar angle, so that a normalized vectorωwrites

ω=





cos(φ)sin(θ) sin(φ)sin(θ)

cos(θ)



.

Other commonly used symbols can be found in the following table:

Symbol Meaning Unit

D(ω) Microflake normal distribution st⁻¹

A(ω) Microscopic self-shadowing function 1

ρ Area of microflake per unit volume m⁻¹

σt Attenuation coefficient m⁻¹

σs,σss,σms Scattering coefficients (wavelength dependent) m⁻¹ f,fss,fms Normalized phase functions st⁻¹ α,αss,αms Albedos (wavelength dependent) 1

h · i Clamped dot product -

χ+(x) Indicator function of positive reals -

2. The microscopic self-shadowing model

We recall here the main expressions of our self-shadowing model. Motivations for this model can be found in the main paper.

2.1. Anisotropic RTE

The anisotropic radiative transfer equation (RTE) introduced by Jakob et al. [JAM^∗10] writes:

(ω· ∇)L(ω) +σt(ω)L(ω) =σs(ω) Z

S²

f(ω→ω⁰)L(ω⁰)dω⁰+Q(ω) (2.1)

withσt(ω)the anisotropic attenuation coefficient,σs(ω)the anisotropic scattering coefficient (which is usually wavelength dependent), and f the anisotropic phase function, in the sense that it depends onωandω⁰, and not just on the angle betweenωandω⁰. In this document, phase function are considered normalized in the second parameter:

Z

S²

f(ω→ω⁰)dω⁰=1, ∀ω. (2.2)

(3)

2.2. Helmholtz’s law of reciprocity

The Helmholtz’s law of reciprocity for anisotropic media writes

σs(ω)f(ω→ω⁰) =σs(ω⁰)f(ω⁰→ω). (2.3)

This law means that the amount of scattered energy from directionωin directionω⁰, given by the scattering coefficientσs(ω)and the phase function, should be equal to the amount of energy scattered from directionω⁰in directionω.

2.3. Self-shadowing function

Before introducing our self-shadowing model, we introduce a directional self-shadowing functionAwith the following properties:

∀ω∈S², 0<A(ω)≤1 and A(ω) =A(−ω). (2.4)

This function represents the probability of self-shadowing at a microscopic scale, depending on the directionω.

2.4. Self-shadowing model

The anisotropic attenuation coefficient in our self-shadowing model writes σt(ω) =A(ω)ρ

Z

D(m)hω·midm. (2.5)

This is the attenuation coefficient in the standard model times the factorA(ω)which takes into account self-shadowing effects. The anisotropic single scattering coefficient is given by

σss(ω) =ραss(λ)A(ω) Z

A(ω⁰)D(m)hω·midm (2.6)

whereω⁰ =2m(m·ω)−ωis the reflected direction, given an input directionωand a microflake normal m. Using the Jacobian of the transformation from normalsω_hto reflected directions provided by Walter et al. [WMLT07],

∂ωh

∂ω⁰

= 1

4|ωh·ω⁰|, we can also write

σss(ω) =ραss(λ)A(ω) 4

Z

A(ω⁰)D(ω_h)dω⁰ (2.7)

withωh= ω+ω⁰

kω+ω⁰kthe half vector. The single scattering phase functionfssis the standard phase function attenuated in some directions due to microscopic shadowing, and re-normalized:

fss(ω→ω⁰) = A(ω⁰)D(ω_h) Z

A(ω⁰⁰)D

ω+ω⁰⁰ kω+ω⁰⁰k

dω⁰⁰

(2.8)

= ραss(λ)A(ω)A(ω⁰)D(ω_h)

4σss(ω) (2.9)

For energy conservation, we introduce a local multiple scattering coefficient:

σms(ω) =αms(λ)ρA(ω) Z

1−A(ω⁰)

D(m)hω·midm (2.10)

=αms(λ)

σt(ω)−σss(ω) αss(λ)

(2.11) and its associated phase function

fms(ω→ω⁰) =fms(ω⁰) = σms(ω⁰)

Rσms(ω⁰⁰)dω⁰⁰= σms(ω⁰) αms(λ)^Rσt(ω⁰⁰)dω⁰⁰−αms(λ)

αss(λ)

Rσss(ω⁰⁰)dω⁰⁰

(2.12)

c

2018 The Author(s)

(4)

2.5. Simplified model with isotropic self-shadowing

Our model greatly simplifies when the self-shadowing function is isotropic, meaning thatA(ω) =A, ∀ω. In the case of isotropic self- shadowing, expressions reduce to

σt(ω) =Aρ Z

D(m)hω·midm (2.13)

σss(ω) =αss(λ)Aσt(ω) (2.14)

fss(ω→ω⁰) =ραss(λ)A²D(ω_h)

4σss(ω) (2.15)

σms(ω) =αms(λ)σt(ω) (1−A) (2.16)

fms(ω⁰) = σms(ω⁰)

Rσms(ω⁰⁰)dω⁰⁰ = σt(ω⁰)

Rσt(ω⁰⁰)dω⁰⁰ (2.17)

Note thatσt(ω)is equal to the projected area of microflakes timesA, andfssis exactly the specular phase function of the standard microflake model.

3. Rejection sampling

Probability distribution functions (pdf) are ubiquitous in physically-based rendering, and procedures that generate samples from such distribution are often needed. In particular, generating samples from phase functions tends to decrease noise in volume path tracing.

Unfortunately, it is not always possible to find an efficient sampling procedure for given apdf. The most common method, known as inverse transform sampling, requires an expression of the inverseF⁻¹of the cumulative distribution function of thepdf:

F(u) =

u Z

x=−∞

pdf(x)dx,

but it is not always possible to find closed-form expressions forF⁻¹. Some authors have found sampling procedures using the fact that theirpdfcomes from a geometric problem – for instance, sampling visible normals of an ellipsoid viewed from a given direction [HDCD15, DWMG15]. This strategy cannot be applied to allpdf.

Rejection sampling is a very general method that can generate samples from any normalizedpdf. For sampling the distributionpdf₁, it is possible to generate samples from any other distributionpdf₂and to accept or reject these samples with a probability such that accepted samples follow exactly the distributionpdf₁. Rejection sampling is very powerful in the sense that it can be used for sampling arbitrary distributions, but the algorithm is only efficient ifpdf₂isclose enoughtopdf₁. The morepdf₂is different frompdf₁, the more samples will be rejected by the algorithm, leading to a very inefficient sampling procedure. For samplingpdf₁usingpdf₂, rejection sampling needs a scalar Msuch that, for allx,

pdf₁(x)≤M pdf₂(x) Then, the sampling procedure is the following:

1. Generate a samplexfrompdf₂ 2. Generate a random numberUin(0,1).

IfU≤ pdf₁(x)

M pdf₂(x), accept the samplexand the procedure is finished.

Else, go back to step 1.

The valueMis the average number of iterations needed before accepting one sample, and 1/Mis the average probability of accepting one sample frompdf₂. It is thus possible to predict how efficient is a sampling procedure – in average.

4. Implementing the self-shadowing model with trigonometric lobes

In this section, we derive closed-form expressions and sampling procedures for implementing our self-shadowing model, using an anisotropic self-shadowing distribution (Sec.4.1) and microflake normal distributions based on trigonometric lobes (Sec.4.2). The derivations include:

• Normalization factors for trigonometric lobes (Sec.4.3).

• Expressions for the attenuation coefficientσt(Sec.4.5).

• Expressions for the single scattering coefficientσss(Sec.4.6).

• Expressions for^Rσtand^Rσss, involved in the multiple scattering albedoσms(Sec.4.7).

(5)

• Sampling procedures for the single scattering phase functionfss(Sec.4.10and4.11).

• Sampling procedures for the multiple scattering phase functionfms(Sec.4.12).

In the next section (Sec.5), we discuss an implementation of our simplified self-shadowing model (when self-shadowing is anisotropic) based on the SGGX distribution.

4.1. Choice of a self-shadowing function

We introduce an anisotropic self-shadowing functionAfor which we can find closed-form expressions for the spherical convolutions and integrals involved in our self-shadowing model (Eq.2.6, 2.12). This anisotropic self-shadowing function writes

A(ω) =a1(ξ_A₁·ω)²+a2(ξ_A₂·ω)²+a3(ξ_A₃·ω)². (4.1) It can also be written as the square projected area of an ellipsoid in directionω[HDCD15] with axesξAi:

A(ω) =ω^T



ξA₁ ξA₂ ξA₃









a₁ 0 0

0 a2 0

0 0 a₃









 ξ^T_A₁

ξ^T_A₂

ξ^T_A₃







ω=ω^TSAω. (4.2)

4.2. Trigonometric lobes:Dcos(ω,ξ,n)andDsin(ω,ξ,n) We use distributions based on trigonometric lobes:

Dcos(ω,ξ,n) =cos²ⁿ(ω,ξ)

Ncos(n) = (ω·ξ)²ⁿ Ncos(n)

Dsin(ω,ξ,n) =sin²ⁿ(ω,ξ)

Nsin(n) = (1−(ω·ξ)²)ⁿ Nsin(n) D_iso(ω) = 1

4π withn∈NandNcos(n)andNsin(n)the normalization factors.

4.3. Normalization factors

For the cosine lobe, the normalization factor writes

Ncos(n) = Z

(m·ξD)²ⁿ dm= 4π 2n+1. The normalization factor for the sine lobe writes

N_sin(n) = Z

1−(m·ξ_D)²n

dm=2π³²Γ(1+n) Γ

3 2+n . In our implementation, we pre-computed the normalization factors for the sine lobe forn∈[1..20].

4.4. Proof that distributionsD_sin(ω,ξ,n)can be written as a sum of distributionsDcos(ω,ξ,n) Sine lobesD_sin(ω,ξ,n)have an interesting property: they can be written as a sum of 2ncosine lobes:

D_sin(ω,ξ,n) = 1 2n

2n−1 p=0

∑

Dcos(ω,ξp,n) (4.3)

withξp being 2nevenly spaced axes in the plane orthogonal toξ. For example, aD_sin(ω,ξ,5)distribution is equivalent to a sum of 10 Dcos(ω,ξp,5)distributions:

= + + + + + + + + +

c

2018 The Author(s)

(6)

Without loss of generality, we derive the proof in spherical coordinates in the caseξ^sin_D = [0,0,1]^T. The sum of cosine lobes writes

1 2n

2n+1 4π

2n−1

∑

p=0











 cos

pπ 2n

sinpπ 2n

0







·





cos(θ)











2n

=2n+1 8nπ

2n−1

∑

p=0

cospπ

2n

cos(φ)sin(θ) +sinpπ 2n

sin(φ)sin(θ)2n

(4.4)

=2n+1 8nπ sin(θ)²ⁿ

2n−1

∑

p=0

cospπ

2n

cos(φ) +sinpπ 2n

sin(φ)2n

(4.5)

2n−1

∑

p=0

cospπ 2n−φ2n

. (4.6)

We can expand cos(x)²ⁿ[Wei17] as follows:

cos(x)²ⁿ= 1 2²ⁿ

2n n

!

+ 1

2²ⁿ⁻¹

n−1

∑

k=0

2n k

!

cos(2(n−k)x) so we can write

2n−1

∑

p=0

cospπ 2n−φ2n

=

2n−1

∑

p=0

1 2²ⁿ

2n n

!

+ 1

2²ⁿ⁻¹

n−1

∑

k=0

2n k

! cos

2(n−k)pπ 2n−φ

! .

We can check that terms involvingφcancel to 0:

2n−1 p=0

∑

n−1

∑

k=0

2n k

! cos

2(n−k)pπ 2n−φ

=

n−1

∑

k=0

2n k

!_2n−1

p=0

∑

cos

2(n−k)pπ 2n−φ

(4.7)

=

n−1

∑

k=0

2n k

!

<

2n−1

∑

p=0

exp

2i(n−k)pπ 2n−φ

!

(4.8)

=

n−1

∑

k=0

2n k

!

< exp(−2i(n−k)φ)

2n−1

∑

p=0

exp

i(n−k)2π 2n

p!

(4.9)

=

n−1

∑

k=0

2n k

!

<







exp(−2i(n−k)φ) 1−exp(i(n−k)2π) 1−exp

i(n−k)2π 2n







(4.10)

=0. (4.11)

Finally, we get

2n+1 8nπ

2n−1

∑

p=0











 cospπ

2n

sin pπ

2n

0







·





cos(θ)











2n

2n−1

∑

p=0

1 2²ⁿ

2n n

!

(4.12)

=sin(θ)²ⁿ2n+1 8nπ

2n 2²ⁿ

2n n

!

(4.13)

=sin(θ)²ⁿ2n+1 π4ⁿ⁺¹

(2n)!

2n! (4.14)

=sin(θ)²ⁿ

√π(2n+2)!

4ⁿ⁺¹(n+1)!

(n+1)

(2n+2)π^3/2n! (4.15)

=sin(θ)²ⁿΓ(n+3/2)

2π^3/2n! (4.16)

=D_sin(ω,ξ^sinD,n). (4.17)

(7)

4.5. Closed-form expression of the attenuation coefficientσt(ω) 4.5.1. Derivation of general closed-from expressions forσt(ω)

A general closed-form expression can be found forσt(ω)but this expression is not convenient for evaluation because it involves costly evaluations of Gauss hypergeometric functions₂F1. For information, we provide a derivation of this expression for the cosine lobe (Eq.4.41).

For efficient implementation, we have pre-computed simpler closed-form expressions for eachn∈[1..20]as explained in Section4.5.2. Note that a general closed-form expression could also be derived for the sine lobe, using the fact that a sine lobe can be written as a sum of cosine lobes (Sec.4.4).

In the case of a cos²ⁿdistribution,σt(ω)is given by

σt(ω) =ρA(ω) Z

Dcos(m,ξD,n)hm·ωidm=ρ A(ω) Ncos(n)

Z

(m·ξD)²ⁿhm·ωidm.

Without loss of generality, we choose a convenient spherical parametrization such thatω= [cos(φ_i),sin(φ_i),0]andξD= [1,0,0]. We have

σt(ω) =ρ A(ω) Ncos(n)

Z









cos(θ)



·



 1 0 0









2n





cos(θ)



·



 cos(φi) sin(φ_i)

0





sin(θ)dφdθ (4.18)

=ρ A(ω) Ncos(n)

π/2−φi

Z

φ=−π/2−φi

π Z

θ=0

(cos(φ)sin(θ))²ⁿ(cos(φ)sin(θ)cos(φi) +sin(φ)sin(θ)sin(φi))sin(θ)dφdθ (4.19)

=ρ A(ω) Ncos(n)

π/2−φi

Z

φ=−π/2−φi

π Z

θ=0

cos(φ)²ⁿ(cos(φ)cos(φi) +sin(φ)sin(φi))sin(θ)²ⁿ⁺²dφdθ (4.20)

=ρ A(ω) Ncos(n)

√ π Γ

₃

2+n Γ(2+n)

π/2−φi

Z

φ=−π/2−φi

cos(φ)²ⁿ(cos(φ)cos(φi) +sin(φ)sin(φi))dφ (4.21)

=ρ A(ω) Ncos(n)

√ π Γ

3 2+n

Γ(2+n)





cos(φ_i)

π/2−φi

Z

φ=−π/2−φi

cos(φ)²ⁿ⁺¹dφ+sin(φ_i)

π/2−φi

Z

φ=−π/2−φi

cos(φ)²ⁿsin(φ)dφ





 (4.22)

=ρ A(ω) Ncos(n)

√ π Γ

3 2+n Γ(2+n)





cos(φi)

π/2−φi

Z

φ=−π/2−φi

cos(φ)²ⁿ⁺¹dφ−2 sin(φ_i)²ⁿ⁺² 2n+1





. (4.23)

Now integral^Rcos(x)^mdxcan be written as a finite sum using Z

cos(x)^mdx= 1

mcos(x)^m−1sin(x) +m+1 m

Z

cos(x)^m−2dx.

It is also possible to write it with the Gauss hypergeometric function. Given that 0≤φi≤π/2, on the interval(−^π₂−φi,−^π₂)we perform a change of variable withX=cos(φ), and get

−π/2 Z

φ=−π/2−φi

cos(φ)²ⁿ⁺¹dφ=

0 Z

X=cos(−π/2−φi)

X²ⁿ⁺¹

√1−X² dX (4.24)

=

0 Z

X=−sin(φi)

X²ⁿ⁺¹

√1−X² dX. (4.25)

c

2018 The Author(s)

(8)

Using the change of variableY= X²

−sin(φ_i), with

∂X

∂Y

= sin(φ_i) 2√

Y we get

0 Z

X=−sin(φi)

X²ⁿ⁺¹

√1−X² dX=

1 Z

0

Yⁿ⁺¹²sin(φ_i)²ⁿ⁺¹ p1−Ysin(φ_i)²

sin(φ_i) 2√

Y dY (4.26)

=−sin(φ_i)²ⁿ⁺² 2

0 Z

1

Yⁿ

p1−Ysin(φ_i)² dY. (4.27)

Using

Γ(b)Γ(c−b)

Γ(c) ²F1(a,b,c,z) = Z1

0

t^b−1(1−t)^c−b−1 (1−tz)^a dt we can write

−sin(φ_i)²ⁿ⁺² 2

0 Z

1

Yⁿ⁺¹²

p1−Ysin(φ_i)² dY=−sin(φ_i)²ⁿ⁺² 2

Γ(b)Γ(c−b)

Γ(c) ²F₁(a,b,c,z) (4.28)

witha=1

2,b=n+1,c=n+2 andz=sin(φi)². We have

−π/2 Z

φ=−π/2−φi

cos(φ)²ⁿ⁺¹dφ=−sin(φi)²ⁿ⁺² 2

Γ(n+1)Γ(1) Γ(n+2) ²F1

1

2,n+1,n+2,sin(φ_i)²

(4.29)

=−sin(φi)²ⁿ⁺² 2n+2 ²F1

1

2,n+1,n+2,sin(φi)²

. (4.30)

Now, we integrate between−π/2 and 0:

0 Z

φ=−π/2

cos(φ)²ⁿ⁺¹dφ= Z 1

0

X²ⁿ⁺¹

√1−X² dX (4.31)

= Z 1

0

Yⁿ⁺¹²

√1−Y 1 2√

Y dY (4.32)

=1 2

Z ₁ 0

Yⁿ

√1−Y dY (4.33)

=Γ(b)Γ(c−b)

2Γ(c) ²F1(a,b,c,1) (4.34)

(4.35) witha=¹₂,b=n+1,c=n+2. We have

0 Z

φ=−π/2

cos(φ)²ⁿ⁺¹dφ=1 2B

1 2,n+1

. (4.36)

Finally, using the same approach,

π/2−φi

Z

φ=0

cos(φ)²ⁿ⁺¹dφ=

π/2 Z

φ=0

cos(φ)²ⁿ⁺¹dφ+

π/2−φi

Z

φ=π/2

cos(φ)²ⁿ⁺¹dφ (4.37)

=

0 Z

1

X²ⁿ⁺¹

−√

1−X² dX+

cos(π/2−φi) Z

0

X²ⁿ⁺¹

−√

1−X² dX (4.38)

= 1

2n+2B 1

2,n+1

−sin(φi)²ⁿ⁺² 2n+2 ²F1

1

2,n+1,n+2,sin(φ_i)²

. (4.39)

(9)

We can add the three terms and get

π/2−φi

Z

φ=−π/2−φi

cos(φ)²ⁿ⁺¹dφ=B 1

2,n+1

−sin(φi)²ⁿ⁺² n+1 ²F1

1

2,n+1,n+2,sin(φi)²

. (4.40)

We proved that

σt(ω) =ρ A(ω) Ncos(n)

√ π Γ

3 2+n

Γ(2+n) cos(φi) B 1

2,n+1

−sin(φi)²ⁿ⁺² n+1 ²F1

1

2,n+1,n+2,sin(φi)² !

−2 sin(φi)²ⁿ⁺² 2n+1

!

. (4.41) As explained before, we don’t use this expression in our implementation because more efficient expressions can be found for each particular value ofn, whennis small.

4.5.2. Pre-computing simple closed-form expressions for particular valuesn

Using symbolic integration [MGH^∗05], we have found that for each particularn,σ^cost (ω)can be written

σ^cost (ω) =ρA(ω)





 1 (ω·ξD)² (ω·ξD)⁴ . . . . (ω·ξD)²ⁿ







·C1

withC₁an+1 vector of coefficients. The evaluation cost is linear innand this expression is very efficient whennis small. In practice we pre-computed coefficient up ton=20 (corresponding to a cos⁴⁰lobe). We have found similar expressions for sine lobes. These coefficients can be found in our code (filesmatrixcoslobe.handmatrixsinlobe.h).

4.6. Single scattering coeffσssforcosⁿandsinⁿmicroflake distributions

General closed-form expressions forσssare even more complex than expressions forσt(Eq.4.41) and have no practical interest for rendering.

Like we did forσt, we pre-computed simpler and exact closed-form solutions for eachnup ton=20. We have found using symbolic integration that

Z

(ξ_A·ω⁰)²Dcos(m,ξD,n)hm·ωidm=





 1 (ω·ξD)² (ω·ξD)⁴ . . . . (ω·ξD)²ⁿ







·C2·



 (ω·ξA)² (ξ_A·ξD)²

1



+ (ξ_A·ξD)(ξ_A·ω)







(ω·ξD) (ω·ξD)³ . . . . . (ω·ξD)²ⁿ⁻¹







·C3 (4.42)

withω⁰=−ω+2m(ω·m),C2a(n+1)×3 matrix of coefficients andC3a vector of coefficients of lengthn. Given Equation4.1, we sum contributions of the three shadowing axes and get

Z

A(ω⁰)Dcos(m,ξD,n)hm·ωidm=





 1 (ω·ξD)² (ω·ξD)⁴ . . . . (ω·ξD)²ⁿ







·C2·





a1(ω·ξ_A1)²+a2(ω·ξA2)²+a3(ω·ξA3)² a1(ξA1·ξD)²+a2(ξA2·ξD)²+a3(ξA3·ξD)²

a1+a2+a3



 (4.43)

+ (a1(ξA1·ξD)(ξA1·ω) +a2(ξA2·ξD)(ξA2·ω) +a3(ξA3·ξD)(ξA3·ω))







(ω·ξD) (ω·ξD)³ . . . . . (ω·ξD)²ⁿ⁻¹







·C3 (4.44)

We can write

a1(ω·ξA1)²+a2(ω·ξA2)²+a3(ω·ξA3)²=ω^TSAω a1(ξA1·ξD)²+a2(ξA2·ξD)²+a3(ξA3·ξD)²=ξ^TDSAξD

a₁(ξ_A1·ξ_D)(ξ_A1·ω) +a₂(ξ_A2·ξD)(ξ_A2·ω) +a₃(ξ_A3·ξD)(ξ_A3·ω) =ω^TS_AξD c

2018 The Author(s)

(10)

Guillaume Loubet / Implementing the self-shadowing model a₁+a₂+a₃=tr(A)

Finally we get, for a cosine lobeDcos,

σss(ω) =αss(λ)ρA(ω)





 1 (ω·ξD)² (ω·ξD)⁴ . . . . (ω·ξD)²ⁿ







·C₂·



 ω^TSAω ξ^T_DS_AξD

trace(A)



+αss(λ)ρA(ω) ω^TS_AξD







(ω·ξD) (ω·ξD)³ . . . . . (ω·ξD)²ⁿ⁻¹







·C₃ (4.45)

Expressions in the case of sine lobesD_sinare similar, only coefficients inC₂andC₃are different. All coefficients can be found in our code in filesmatrixcoslobe.handmatrixsinlobe.h.

4.7. Closed-form expression for^Rσtand^Rσss

Integrals^Rσt(ω)dωand^Rσss(ω)dωare involved in our multiple scattering phase function (Eq.2.12). We provide here derivations of their closed-from expressions.

4.7.1. Closed-form expression of^Rσt(ω)for cosine lobesDcos(m,ξ_D,n) We want to find a closed-form expression for a cosine distribution:

Z

σt(ω)dω=ρ Z

A(ω) Z

Dcos(m,ξD,n)hm·ωidmdω. (4.46)

We first invert integrals:

Z A(ω)

Z

Dcos(m,ξ_D,n)hm·ωidmdω= Z

Dcos(m,ξD,n) Z

A(ω)hm·ωidωdm. (4.47) We can decompose our shadowing function (Eq.4.1):

Z

A(ω)hm·ωidω=

∑

^aⁱ

Z

(ω·ξ_A_i)²hm·ωidω. (4.48)

We choose a convenient parametrization: we takem= [1,0,0]^T andξ_A_i= [sin(θ_A),0,cos(θ_A)]^T, so that we can write Z

(ω·ξAi)²hm·ωidω=

2π Z

φ=0 π/2 Z

θ=0

(sin(θ_A)cos(φ)sin(θ) +cos(θ_A)cos(θ))²cos(θ)sin(θ)dφdθ=π 4

cos(θ_A)²+1

= π 4

(ξ_A_i·m)²+1 . (4.49) We can now substitute this expression in Eq.4.47and get

Z

Dcos(m,ξ_D,n)π 4

(ξ_A_i·m)²+1

dm=

Z 2n+1

4π (ξ_D·m)²ⁿπ 4

(ξ_A_i·m)²+1

dm. (4.50)

Choosing a parametrization such thatξD= [0,0,1]^T andξAi= [sin(θ_A),0,cos(θ_A)]^T, this writes

2π Z

φ=0 Zπ

θ=0

2n+1

4π cos(θ)²ⁿπ 4

(sin(θ_A)cos(φ)sin(θ) +cos(θ_A)cos(θ))²+1

sin(θ)dm= π

n(ξ_D·ξAi)²+n+2

4n+6 . (4.51)

Finally, we can sum and factor the formula for the three axes, and we get Z

σt(ω)dω=ρa1

π

n(ξ_D·ξA₁)²+n+2

4n+6 +ρa2

π

n(ξ_D·ξA₂)²+n+2

4n+6 +ρa3

π

n(ξ_D·ξA₃)²+n+2

4n+6 . (4.52)

Using Eq.4.2, we can write, for a distributionDcos(m,ξ_D,n):

Z

σt(ω)dω=ρ Z

A(ω) Z

Dcos(m,ξD,n)hm·ωidmdω=ρ π

n(ω^TSAω) +tr(SA)(n+2)

4n+6 (4.53)

(11)

4.7.2. Closed-form expression of^Rσt(ω)for sine lobeD_sin(m,ξ_D,n)

The derivation is almost the same as for theDcos(m,ξ_D,n)in the previous subsection. We need to compute Z

Dsin(m,ξ_D,n)π 4

(ξ_A_i·m)²+1

dm=

Z 2π³²Γ(1+n) Γ

3 2+n

1−(ξ_D·m)²nπ 4

(ξ_A_i·m)²+1

dm. (4.54)

Choosing the same convenient frame in whichξD= [0,0,1]^TandξAi= [sin(θA),0,cos(θA)]^T, this writes 2×

2π Z

φ=0 π/2 Z

θ=0

π 4

(sin(θ_A)cos(φ)sin(θ) +cos(θ_A)cos(θ))²+1

sin(θ)²ⁿ⁺¹dm= π

3n+4−n(ξD·ξ_A_i)²

8n+12 . (4.55)

Again, we have to sum this for each shadowing axis. For aDsin(m,ξ_D,n)distribution, we get:

Z

σt(ω)dω=ρ Z

A(ω) Z

Dsin(m,ξ_D,n)hm·ωidmdω=ρ π

(3n+4)tr(S_A)−n(ω^TSAω)

8n+12 (4.56)

4.7.3. Closed-form expression of^Rσss(ω)for cosine lobeDcos(m,ξ_D,n)

In the case of a distributionDcos(m,ξD,n), we look for a closed-form expression of the integral Z

σss(ω)dω=ραss(λ) Z

A(ω) Z

A(ω⁰)Dcos(m,ξD,n)hω·midmdω, (4.57)

withω⁰=−ω+2m(ω·m). We invert integrals and get Z

Dcos(m,ξD,n) Z

A(ω)A(ω⁰)hω·midωdm. (4.58)

We expand each shadowing term using Equation4.1and get 9 terms of the form:

Z

Dcos(m,ξD,n) Z

A(ω)A(ω⁰)hω·midωdm=

∑

i,j∈{1,2,3}

Z

Dcos(m,ξD,n)aiaj Z

(ω·ξAi)²(ω⁰·ξA j)²hω·midω. (4.59) We first consider the casei= j, and then the casei6= j.

4.7.3.1. Casei= j

Without loss of generality, we work in the reference frame in whichm= [1,0,0]andξAi= [cos(φ_Ai),sin(φ_Ai),0]. We find Z

(ω·ξAi)²((−ω+2m(ω·m))·ξAi)²hω·midω (4.60)

=

π/2 Z

φ=−π/2 Zπ

θ=0

(cos(φ_Ai)cos(φ) +sin(φ_Ai)sin(φ))²(cos(φ_Ai)cos(φ)−sin(φ_Ai)sin(φ))²cos(φ)sin(θ)⁶dωdθ (4.61)

= 1 24π

15(m·ξ_Ai)⁴−10(m·ξAi)²+3

. (4.62)

We integrate this result withDcos(m,ξD,n)(recall that we want a closed-form expression for Equation4.57) and find Z

Dcos(m,ξ_D,n)1 24π

15(m·ξAi)⁴−10(m·ξAi)²+3

dm=π15n(n−1) (ξD·ξAi)⁴−10n(n−2) (ξD·ξAi)²+3n²+7n+10

24n²+96n+90 . (4.63)

4.7.3.2. Casei6= j

We choose a reference frame such thatm= [0,0,1]^T,Ai= [sin(θAi),0,cos(θAi)]^TandAj= [cos(φA j)sin(θA j),sin(φA j)sin(θA j),cos(θA j)]^T. We have

ξAi⊥ξA j⇒sin(θ_Ai)cos(φ_{A j})sin(θ_{A j}) +cos(θ_{A j})cos(θ_Ai) =0.

We develop and integrate as before, and find Z

(ω·ξ_Ai)²((−ω+2m(ω·m))·ξ_{A j})²hω·midω= π 24

15(m·ξ_Ai)²(m·ξ_{A j})²+ (m·ξ_Ai)²+ (m·ξ_{A j})²+1

. (4.64)

c

2018 The Author(s)

(12)

Guillaume Loubet / Implementing the self-shadowing model Finally, we get

Z

Dcos(m,ξ_D,n)π 24

15(m·ξAi)²(m·ξA j)²+ (m·ξAi)²+ (m·ξA j)²+1

dm=π15C1C2n(n−1) + (n²+10n)(C1+C2) +n²+5n+10 24n²+96n+90

(4.65) withC1= (ξ_Ai·ξD)²andC2= (ξ_{A j}·ξD)².

4.7.3.3. Sum of all terms

Given Equation4.59,4.63and4.65, we can write the closed-form expression of^Rσss(ω)in the case of a cosine lobeDcos(m,ξD,n). After simplification, we get:

Z

σss(ω)dω= πραss(λ) 24n²+96n+90

15n(n−1) ξ^TDS_AξD

2

−10n(n−2)ξ^T_DSASAξD+2n(n+10)

ξ^TDSAξDtr(S_A)−ξ^TDSASAξD

+ (3n²+7n+10)tr(SASA) + (n²+5n+10)

tr(SA)²−tr(SASA)

(4.66)

4.7.4. Closed-form expression of^Rσss(ω)for sine lobeDsin(m,ξ_D,n) The derivation is almost the same as for cosine lobes.

4.7.4.1. Casei= j

We have shown (Sec.4.7.3.1) that Z

(ω·ξ_A_i)²((−ω+2m(ω·m))·ξ_Ai)²hω·midω= π 24

15(m·ξ_Ai)⁴−10(m·ξ_Ai)²+3

(4.67) Then

Z

Dsin(m,ξD,n)π 24

15(m·ξAi)⁴−10(m·ξAi)²+3

dm=π(45n²−45n)(ξD·ξAi)⁴+ (−50n²+10n)(ξD·ξAi)²+29n²+91n+80 192n²+768n+720

(4.68) 4.7.4.2. Casei6= j

We have found previously (Sec.4.7.3.2) that Z

(ω·ξ_Ai)²((−ω+2m(ω·m))·ξ_{A j})²hω·midω= 1 24π

(4.69) We get

Z

D_sin(m,ξD,n)π 24

dm=π45C₁C₂n(n−1)−(C₁+C₂)n(19n+25) +31n²+105n+80 192n²+768n+720

(4.70) WithC1= (ξ_A1·ξD)²andC2= (ξ_A2·ξD)²

4.7.4.3. Sum of all terms

Again, we sum terms from casesi= jand casesi6= j, and after simplification we get the closed-form expression of^Rσssfor sine lobes Dsin(m,ξD,n):

Z

σss(ω)dω= πραss(λ) 192n²+768n+720

45n(n−1) ξ^TDSAξD

2

+10n(1−5n)ξ^TDSASAξD

−2n(19n+25)

ξ^T_DS_AξDtr(S_A)−ξ^T_DS_AS_AξD

+ (29n²+91n+80)tr(S_AS_A) + (31n²+105n+80)

tr(S_A)²−tr(S_AS_A)

(4.71)

(13)

4.8. Sampling the distribution of visible normals for a cosine lobeDcos(m,ξ_D,n)

We address here the problem of sampling the distribution of visible microflake normals [HDCD15] given that the microflake distribution is a cosine lobeDcos(m,ξD,n). This sampling procedure is useful when using trigonometric lobes in the standard microflake model or in the simplified self-shadowing model. Indeed, sampling visible normal distributions is actually the same problem as sampling the specular phase function. We want to sample the probability density function

Dωi(ωo) = Dcos(ω_h,ξD,n)

4^RDcos(m,ξD,n)hωi·midm (4.72)

withωh= ωi+ωo

kω_i+ωok. It is more convenient to work in the space of microflake normals, and to sample normals instead of outgoing directions.

Using the Jacobian of the transformation from microflake normals to outgoing directions [WMLT07], thepdffor sampling microflake normals writes

Dωi(m) = Dcos(m,ξ_D,n)hωi·mi

RDcos(m,ξ_D,n)hωi·midm. (4.73)

We could not find a procedure for sampling directly this function. We propose a sampling procedure based on rejection sampling (Sec.3).

4.8.1. Using rejection sampling for sampling visible normals ofDcos(m,ξ_D,n) We have found that it is possible to sample the function

Dcos(m,ξD,n)(ωi·m)² RDcos(m,ξD,n)(ωi·m)²dm

which is relatively similar toDωi(m)(Eq.4.73) on the hemisphere for whichDωi(m)is not null. We propose a rejection sampling method based on this function. We introduce the normalizedpdf

D^s_ω_i(m) = Dcos(m,ξD,n)((ω1·m)²+1/4) RDcos(m,ξD,n)((ω₁·m)²+1/4)dm because we have

|ω₁·m| ≤(ω₁·m)²+1/4 so that we can write

Dcos(m,ξ_D,n)|ωi·m|

RDcos(m,ξ_D,n)|ω_i·m|dm≤MD^s_ω_i(m) with

M=

RDcos(m,ξ_D,n)((ω₁·m)²+1/4)dm RDcos(m,ξ_D,n)|m·ωi|dm . Using rejection sampling, we can generate samples from thepdf

Dcos(m,ξD,n)|ωi·m|

RDcos(m,ξD,n)|ωi·m|dm by generating samples fromD^s_ω_i(m)and accepting samples with probability

Dcos(m,ξ_D,n)|ωi·m|

RDcos(m,ξ_D,n)|ω_i·m|dm 1

MD^s_ω_i(m)= |ω1·m|

(ω₁·m)²+1/4. Once we get a samplem, we invert its direction ifωi·m<0 and get a sample forD_ωi.

4.8.2. SamplingD^s_ωi(m)

We show here how to sample thepdf

D^s_ω_i(m) = Dcos(m,ξD,n)((ω₁·m)²+1/4) RDcos(m,ξD,n)((ω₁·m)²+1/4)dm. We choose a reference frame such thatξD= [0,0,1]^T andω= [sin(ti),0,cos(ti)]^T, and find

D^s_ωi(φ,θ) =cos(θ)²ⁿ((cos(t_i)cos(φ)sin(θ) +sin(t_i)cos(θ))²+1/4) (2n+1)(2n+3) π(8ncos(ti)²+2n+7).

c

2018 The Author(s)

(14)

The marginalpdfforθis given by f_θ(θ) =

Z

D^s_ω_i(φ,θ)sin(θ)dφ=(2n+3)(2n+1)cos(θ)²ⁿsin(θ) 8 cos(ti)²n+2n+7

cos(θ)²

3 cos(ti)²−1

−cos(ti)²+3 2

.

For sampling this function, we use inverse transform sampling, meaning that we sample a random numberU∈(0,1)and solve forθthe equation

U= Z_θ

0 f_θ(x)dx.

We find that

Z_θ

0 f_θ(x)dx= −cos(θ)²ⁿ⁺¹ 8 cos(ti)²n+2n+7

(2n+1)

3 cos(t_i)²−1

cos(θ)²+ (2n+3)

−cos(t_i)²+3 2

.

In our implementation, we solve this function numerically using Brent’s method, although other methods could also be used. Onceθhas been sampled,φmust be sampled from thepdf

D^s_ωi(φ,θ) f_θ(θ) . This leads to the following equation for samplingφ:

V= Z_φ

0

D^s_ω_i(x,θ)

f_θ(θ) dx=2 sin(t_i)²sin(θ)²(cos(φ)sin(φ) +φ) +8 sin(θ)cos(θ)cos(t_i)sin(t_i)sin(φ) +4 cos(θ)²cos(t_i)²φ+φ

2π(6 cos(θ)²cos(t_i)²−2 cos(θ)²−2 cos(t_i)²+3)) , (4.74) V being another random number in(0,1). Again, we invert this equation numerically with Brent’s method, and we obtain a samplem=





cos(θ)



.

4.9. Sampling the distribution of visible normals for a cosine lobeDsin(m,ξD,n)

For sine lobes, we use the fact thatDsin(m,ξ_D,n)can be written as a sum of cosine lobesDcos(m,ξ_D,n)(Sec.4.4). We first sample one of the cosine lobes proportionally to their projected area in directionωi, and then sample a visible normal from this lobe using the procedure derived in the previous subsection (Sec.4.8).

4.10. Sampling single scattering phase functionfssforDcos(m,ξD,n)

In subsections4.8and4.9, we have derived sampling procedures for the distribution of visible normals. Now, we derive similar sampling procedures for the single scattering phase function fss. The single scattering phase function is given by

fss(ωi→ωo) =ραss(λ)A(ωo)A(ω_i)D(ω_h) 4σs(ωi) . We useD^s_ω_i(Eq.4.8.1) in the space of outgoing directions

D^s_ωi(ωo) = D(ωh)((ωh·ωi)²+1/4) 4|ω_h·ωi|^RD(ω_h)((ω_h·ωi)²+1/4)dm and write

fss(ω_i→ωo) =ραss(λ)A(ωo)A(ωi)D(ωh)

4σs(ω_i) (4.75)

=D^s_ω_i(ωo)ραss(λ)A(ωo)A(ω_i) σs(ω_i)

|ωh·ωi|^RD(ω_h)((ω_h·ωi)²+1/4)dm

((ω_h·ωi)²+1/4) (4.76)

≤D^s_ωi(ωo)ραss(λ)AmaxA(ωi)^RD(ωh)((ωh·ωi)²+1/4)dm

σs(ω_i) (4.77)

≤D^s_ωi(ωo)M_A (4.78)

withAmax=max

ω

A(ω)and

MA= AmaxR

D(m)((m·ωi)²+1/4)dm RA(−ωi+2m(m·ωi))D(m)|m·ωi|dm.

(15)

Using this relation betweenD^s_ω_i and fss, with can generate samples from fss by generating samples fromD^s_ω_i and accepting them with probability

fss(ω_i→ωo)

MD^s_ω_i(ω) = A(ωo)|ω_h·ωi|

Amax((ω_h·ωi)²+1/4). (4.79)

4.11. Sampling single scattering phase functionfssforDsin(m,ξD,n)

In the case of a sine distributionDsin(m,ξ_D,n), we rely on the sampling procedure for visible normals proposed in subsection4.9. This method generates samples from thepdf

D_ωi(ωo) = D(ω_h) 4^RD(m)|m·ω_i|dm. We can write

fss(ω_i→ωo) =ραss(λ)A(ω_i)A(ωo)D(ω_h)

4σs(ω_i) (4.80)

≤Dωi(ωo)ραss(λ)AmaxA(ωi)^RD(m)|m·ωi|dm

σs(ω_i) (4.81)

(4.82) so we can sample fssfor a sine distributionD_sin(m,ξD,n)by generating samplesωofromDωi(ωo)and accepting them with probability

pdf₁(ωo)

Mpdf₂(ωo)=ραss(λ)A(ω_i)A(ωo)D(ω_h) 4σs(ωi)

4^RD(m)|m·ωi|dm D(ω_h)

σs(ω_i)

ραss(λ)AmaxA(ωi)^RD(m)|m·ωi|dm (4.83)

=A(ωo) Amax

. (4.84)

4.12. Sampling the multiple scattering phase function fms

Our multiple scattering phase function is given by

fms(ωo) = σt(ωo)−σs(ωo)/αss(λ)

Rσt(ω)−^Rσs(ω)/αss(λ). (4.85)

We observe that in the case of isotropic self-shadowing, this phase function reduces to thepdf A(1−A)^RD(m)hm·ωoidm

RA^RD(m)hm·ωoidmdω−^RA^RAD(m)hm·ωoidmdω=

RD(m)hm·ωoidm R RD(m)hm·ωoidmdω

which is proportional the projected area of microflakes in directionωo. Given a single microflake with normalm, we can sample a direction ωproportionally to the projection area of the microflake|m·ω|by sampling a cosine lobe around directionm. Similarly, we can sample the projected area of a cloud of microflakes by sampling first a normal from the distribution of normalsD, and then by sampling a cosine lobe around this normal. In subsection4.13and4.14, we show how to sample normals forDcosandDsin.

We derive a rejection sampling procedure forfmsbased on thispdf:

σt(ωo)−σs(ωo)/αss(λ)

Rσt(ω)−^Rσs(ω)/αss(λ)= A(ωo)^R(1−A(−ωo+2m(m·ωo)))D(m)hm·ωoidm

RA(ω)^R(1−A(−ω+2m(m·ω)))D(m)hm·ωidmdω (4.86)

<

RD(m)hm·ωoidm R RD(m)hm·ωoidmdωo

Amax(1−Amin)^{R R}D(m)hm·ωoidmdωo

RA(ω)^R(1−A(−ω+2m(m·ω)))D(m)hm·ωidmdω. (4.87) We find the following probability of acceptance:

A(ω)^R(1−A(−ω+2m(m·ω)))D(m)hm·ωidm

Amax(1−Amin)^RD(m)hm·ωoidm . (4.88) This probability tends to one when self-shadowing is almost anisotropic, and tends to be smaller whenAminandAmax are very different, meaning that this sampling procedure is less efficient is such case.

4.13. Sampling a normal fromDcos

The distribution

Dcos(ω,ξ,n) = (ω·ξ)²ⁿ2n+1 4π

c

2018 The Author(s)

(16)

can be sampled easily in the reference frame in whichξ= [0,0,1], using inverse transform sampling. We first sample the azimuthφuniformly in(0,2π). Then, the marginal distribution forθis

2π Z

0

cos(θ)²ⁿ2n+1

4π sin(θ)dφ= 2n+1

2 cos(θ)²ⁿsin(θ).

The cumulative distribution function writes

x Z

0

2n+1

2 cos(θ)²ⁿsin(θ)dθ= 1 2

1−cos(x)²ⁿ⁺¹ .

This function can be inverted:

U=1 2

1−cos(x)²ⁿ⁺¹

⇒cos(x)²ⁿ⁺¹=1−2U, IfU≥1/2:

x=arccos

(1−2U)²ⁿ⁺¹¹

. (4.89)

IfU≤1/2:

x=π−arccos

(2U−1)²ⁿ⁺¹¹

. (4.90)

Therefore,θcan be sampled by samplingUin(0,1)and then using Eq.4.89or4.90. Finally, the sample m=





cos(θ)





must be rotated back in the original reference frame ofDcos(ω,ξ,n).

4.14. Sampling a normal fromDsin

Again, we use the fact thatD_sin(m,ξD,n)can be written as a sum of cosine lobesDcos(m,ξD,n)(Sec.4.4). We first sample one of the cosine lobes randomly and follow the procedure derived in the previous subsection (Sec.4.13).

4.15. Sampling a cosine lobe

Sampling simple cosine lobes cos(θ)is very common in path tracing, we recall here the procedure. Given a lobe axism, sampling a cosine lobe around this direction can be done easily in the frame in whichm= [0,0,1]. In this frame, the cosine lobe on the upper hemisphere in spherical coordinates writes

cos(θ)

π 1θ∈(0,^π₂)(θ).

We first sampleφuniformly in(0,2π). The marginalpdfforθwrites

2π Z

0

cos(θ)

π sin(θ)dφ=2 cos(θ)sin(θ).

Then the cumulative distribution function writes

x Z

0

2 cos(θ)sin(θ)dθ=sin(x)².

Soθcan be sampled using a random sampleUin(0,1)and the inverse of the cumulative function:

θ=arcsin√ U

.