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^{n}and sin^{n}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)

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π}^{1}
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^{−1}

A(ω) Microscopic self-shadowing function 1

ρ Area of microflake per unit volume m^{−1}

σt Attenuation coefficient m^{−1}

σs,σss,σms Scattering coefficients (wavelength dependent) m^{−1}
f,fss,fms Normalized phase functions st^{−1}
α,α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^{2}

f(ω→ω^{0})L(ω^{0})dω^{0}+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ω^{0}, and not just on the angle betweenωandω^{0}. In this document,
phase function are considered normalized in the second parameter:

Z

S^{2}

f(ω→ω^{0})dω^{0}=1, ∀ω. (2.2)

2.2. Helmholtz’s law of reciprocity

The Helmholtz’s law of reciprocity for anisotropic media writes

σs(ω)f(ω→ω^{0}) =σs(ω^{0})f(ω^{0}→ω). (2.3)

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

2.3. Self-shadowing function

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

∀ω∈S^{2}, 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(ω^{0})D(m)hω·midm (2.6)

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

∂ωh

∂ω^{0}

= 1

4|ωh·ω^{0}|,
we can also write

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

Z

A(ω^{0})D(ω_{h})dω^{0} (2.7)

withωh= ω+ω^{0}

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

fss(ω→ω^{0}) = A(ω^{0})D(ω_{h})
Z

A(ω^{00})D

ω+ω^{00}
kω+ω^{00}k

dω^{00}

(2.8)

= ραss(λ)A(ω)A(ω^{0})D(ω_{h})

4σss(ω) (2.9)

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

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

1−A(ω^{0})

D(m)hω·midm (2.10)

=αms(λ)

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

(2.11) and its associated phase function

fms(ω→ω^{0}) =fms(ω^{0}) = σms(ω^{0})

Rσms(ω^{00})dω^{00}= σms(ω^{0})
αms(λ)^{R}σt(ω^{00})dω^{00}−αms(λ)

αss(λ)

Rσss(ω^{00})dω^{00}

(2.12)

c

2018 The Author(s)

Guillaume Loubet / Implementing the self-shadowing model

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(ω→ω^{0}) =ραss(λ)A^{2}D(ω_{h})

4σss(ω) (2.15)

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

fms(ω^{0}) = σms(ω^{0})

Rσms(ω^{00})dω^{00} = σt(ω^{0})

Rσt(ω^{00})dω^{00} (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 distri- bution 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^{−1}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^{−1}. 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_{1}, it
is possible to generate samples from any other distributionpdf_{2}and to accept or reject these samples with a probability such that accepted
samples follow exactly the distributionpdf_{1}. Rejection sampling is very powerful in the sense that it can be used for sampling arbitrary
distributions, but the algorithm is only efficient ifpdf_{2}isclose enoughtopdf_{1}. The morepdf_{2}is different frompdf_{1}, the more samples will be
rejected by the algorithm, leading to a very inefficient sampling procedure. For samplingpdf_{1}usingpdf_{2}, rejection sampling needs a scalar
Msuch that, for allx,

pdf_{1}(x)≤M pdf_{2}(x)
Then, the sampling procedure is the following:

1. Generate a samplexfrompdf_{2}
2. Generate a random numberUin(0,1).

IfU≤ pdf_{1}(x)

M pdf_{2}(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_{2}. 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).

• 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}_{1}·ω)^{2}+a2(ξ_{A}_{2}·ω)^{2}+a3(ξ_{A}_{3}·ω)^{2}. (4.1)
It can also be written as the square projected area of an ellipsoid in directionω[HDCD15] with axesξAi:

A(ω) =ω^{T}

ξA_{1} ξA_{2} ξA_{3}

a_{1} 0 0

0 a2 0

0 0 a_{3}

ξ^{T}_{A}_{1}

ξ^{T}_{A}_{2}

ξ^{T}_{A}_{3}

ω=ω^{T}SAω. (4.2)

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

Dcos(ω,ξ,n) =cos^{2n}(ω,ξ)

Ncos(n) = (ω·ξ)^{2n}
Ncos(n)

Dsin(ω,ξ,n) =sin^{2n}(ω,ξ)

Nsin(n) = (1−(ω·ξ)^{2})^{n}
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)^{2n} dm= 4π
2n+1.
The normalization factor for the sine lobe writes

N_{sin}(n) =
Z

1−(m·ξ_{D})^{2}n

dm=2π^{3}^{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)

Guillaume Loubet / Implementing the self-shadowing model

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(φ)sin(θ) sin(φ)sin(θ)

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}

2n−1

### ∑

p=0

cospπ

2n

cos(φ) +sinpπ 2n

sin(φ)2n

(4.5)

=2n+1
8nπ sin(θ)^{2n}

2n−1

### ∑

p=0

cospπ 2n−φ2n

. (4.6)

We can expand cos(x)^{2n}[Wei17] as follows:

cos(x)^{2n}= 1
2^{2n}

2n n

!

+ 1

2^{2n−1}

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}

2n n

!

+ 1

2^{2n−1}

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(φ)sin(θ) sin(φ)sin(θ)

cos(θ)

2n

=2n+1
8nπ sin(θ)^{2n}

2n−1

### ∑

p=0

1
2^{2n}

2n n

!

(4.12)

=sin(θ)^{2n}2n+1
8nπ

2n
2^{2n}

2n n

!

(4.13)

=sin(θ)^{2n}2n+1
π4^{n+1}

(2n)!

2n! (4.14)

=sin(θ)^{2n}

√π(2n+2)!

4^{n+1}(n+1)!

(n+1)

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

=sin(θ)^{2n}Γ(n+3/2)

2π^{3/2}n! (4.16)

=D_{sin}(ω,ξ^{sin}D,n). (4.17)

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_{2}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^{2n}distribution,σt(ω)is given by

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

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

Z

(m·ξD)^{2n}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(φ)sin(θ) sin(φ)sin(θ)

cos(θ)

·

1 0 0

2n

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

cos(θ)

·

cos(φi)
sin(φ_{i})

0

sin(θ)dφdθ (4.18)

=ρ A(ω) Ncos(n)

π/2−φi

Z

φ=−π/2−φi

π Z

θ=0

(cos(φ)sin(θ))^{2n}(cos(φ)sin(θ)cos(φi) +sin(φ)sin(θ)sin(φi))sin(θ)dφdθ (4.19)

=ρ A(ω) Ncos(n)

π/2−φi

Z

φ=−π/2−φi

π Z

θ=0

cos(φ)^{2n}(cos(φ)cos(φi) +sin(φ)sin(φi))sin(θ)^{2n+2}dφdθ (4.20)

=ρ A(ω) Ncos(n)

√ π Γ

_{3}

2+n Γ(2+n)

π/2−φi

Z

φ=−π/2−φi

cos(φ)^{2n}(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(φ)^{2n+1}dφ+sin(φ_{i})

π/2−φi

Z

φ=−π/2−φi

cos(φ)^{2n}sin(φ)dφ

(4.22)

=ρ A(ω) Ncos(n)

√ π Γ

3 2+n Γ(2+n)

cos(φi)

π/2−φi

Z

φ=−π/2−φi

cos(φ)^{2n+1}dφ−2 sin(φ_{i})^{2n+2}
2n+1

. (4.23)

Now integral^{R}cos(x)^{m}dxcan be written as a finite sum using
Z

cos(x)^{m}dx= 1

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

Z

cos(x)^{m−2}dx.

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

−π/2 Z

φ=−π/2−φi

cos(φ)^{2n+1}dφ=

0 Z

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

X^{2n+1}

√1−X^{2} dX (4.24)

=

0 Z

X=−sin(φi)

X^{2n+1}

√1−X^{2} dX. (4.25)

c

2018 The Author(s)

Guillaume Loubet / Implementing the self-shadowing model

Using the change of variableY= X^{2}

−sin(φ_{i}), with

∂X

∂Y

= sin(φ_{i})
2√

Y we get

0 Z

X=−sin(φi)

X^{2n+1}

√1−X^{2} dX=

1 Z

0

Y^{n+}^{1}^{2}sin(φ_{i})^{2n+1}
p1−Ysin(φ_{i})^{2}

sin(φ_{i})
2√

Y dY (4.26)

=−sin(φ_{i})^{2n+2}
2

0 Z

1

Y^{n}

p1−Ysin(φ_{i})^{2} dY. (4.27)

Using

Γ(b)Γ(c−b)

Γ(c) ^{2}F1(a,b,c,z) =
Z1

0

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

−sin(φ_{i})^{2n+2}
2

0 Z

1

Y^{n+}^{1}^{2}

p1−Ysin(φ_{i})^{2} dY=−sin(φ_{i})^{2n+2}
2

Γ(b)Γ(c−b)

Γ(c) ^{2}F_{1}(a,b,c,z) (4.28)

witha=1

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

−π/2 Z

φ=−π/2−φi

cos(φ)^{2n+1}dφ=−sin(φi)^{2n+2}
2

Γ(n+1)Γ(1)
Γ(n+2) ^{2}F1

1

2,n+1,n+2,sin(φ_{i})^{2}

(4.29)

=−sin(φi)^{2n+2}
2n+2 ^{2}F1

1

2,n+1,n+2,sin(φi)^{2}

. (4.30)

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

0 Z

φ=−π/2

cos(φ)^{2n+1}dφ=
Z 1

0

X^{2n+1}

√1−X^{2} dX (4.31)

= Z 1

0

Y^{n+}^{1}^{2}

√1−Y 1 2√

Y dY (4.32)

=1 2

Z _{1}
0

Y^{n}

√1−Y dY (4.33)

=Γ(b)Γ(c−b)

2Γ(c) ^{2}F1(a,b,c,1) (4.34)

(4.35)
witha=^{1}_{2},b=n+1,c=n+2. We have

0 Z

φ=−π/2

cos(φ)^{2n+1}dφ=1
2B

1 2,n+1

. (4.36)

Finally, using the same approach,

π/2−φi

Z

φ=0

cos(φ)^{2n+1}dφ=

π/2 Z

φ=0

cos(φ)^{2n+1}dφ+

π/2−φi

Z

φ=π/2

cos(φ)^{2n+1}dφ (4.37)

=

0 Z

1

X^{2n+1}

−√

1−X^{2} dX+

cos(π/2−φi) Z

0

X^{2n+1}

−√

1−X^{2} dX (4.38)

= 1

2n+2B 1

2,n+1

−sin(φi)^{2n+2}
2n+2 ^{2}F1

1

2,n+1,n+2,sin(φ_{i})^{2}

. (4.39)

We can add the three terms and get

π/2−φi

Z

φ=−π/2−φi

cos(φ)^{2n+1}dφ=B
1

2,n+1

−sin(φi)^{2n+2}
n+1 ^{2}F1

1

2,n+1,n+2,sin(φi)^{2}

. (4.40)

We proved that

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

√ π Γ

3 2+n

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

2,n+1

−sin(φi)^{2n+2}
n+1 ^{2}F1

1

2,n+1,n+2,sin(φi)^{2}
!

−2 sin(φi)^{2n+2}
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,σ^{cos}t (ω)can be written

σ^{cos}t (ω) =ρA(ω)

1
(ω·ξD)^{2}
(ω·ξD)^{4}
. . . .
(ω·ξD)^{2n}

·C1

withC_{1}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^{40}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^{n}andsin^{n}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}·ω^{0})^{2}Dcos(m,ξD,n)hm·ωidm=

1
(ω·ξD)^{2}
(ω·ξD)^{4}
. . . .
(ω·ξD)^{2n}

·C2·

(ω·ξA)^{2}
(ξ_{A}·ξD)^{2}

1

+ (ξ_{A}·ξD)(ξ_{A}·ω)

(ω·ξD)
(ω·ξD)^{3}
. . . . .
(ω·ξD)^{2n−1}

·C3 (4.42)

withω^{0}=−ω+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(ω^{0})Dcos(m,ξD,n)hm·ωidm=

1
(ω·ξD)^{2}
(ω·ξD)^{4}
. . . .
(ω·ξD)^{2n}

·C2·

a1(ω·ξ_{A1})^{2}+a2(ω·ξA2)^{2}+a3(ω·ξA3)^{2}
a1(ξA1·ξD)^{2}+a2(ξA2·ξD)^{2}+a3(ξA3·ξD)^{2}

a1+a2+a3

(4.43)

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

(ω·ξD)
(ω·ξD)^{3}
. . . . .
(ω·ξD)^{2n−1}

·C3 (4.44)

We can write

a1(ω·ξA1)^{2}+a2(ω·ξA2)^{2}+a3(ω·ξA3)^{2}=ω^{T}SAω
a1(ξA1·ξD)^{2}+a2(ξA2·ξD)^{2}+a3(ξA3·ξD)^{2}=ξ^{T}DSAξD

a_{1}(ξ_{A1}·ξ_{D})(ξ_{A1}·ω) +a_{2}(ξ_{A2}·ξD)(ξ_{A2}·ω) +a_{3}(ξ_{A3}·ξD)(ξ_{A3}·ω) =ω^{T}S_{A}ξD
c

2018 The Author(s)

Guillaume Loubet / Implementing the self-shadowing model
a_{1}+a_{2}+a_{3}=tr(A)

Finally we get, for a cosine lobeDcos,

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

1
(ω·ξD)^{2}
(ω·ξD)^{4}
. . . .
(ω·ξD)^{2n}

·C_{2}·

ω^{T}SAω
ξ^{T}_{D}S_{A}ξD

trace(A)

+αss(λ)ρA(ω)
ω^{T}S_{A}ξD

(ω·ξD)
(ω·ξD)^{3}
. . . . .
(ω·ξD)^{2n−1}

·C_{3} (4.45)

Expressions in the case of sine lobesD_{sin}are similar, only coefficients inC_{2}andC_{3}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}

^{i}

Z

(ω·ξ_{A}_{i})^{2}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)^{2}hm·ωidω=

2π Z

φ=0 π/2 Z

θ=0

(sin(θ_{A})cos(φ)sin(θ) +cos(θ_{A})cos(θ))^{2}cos(θ)sin(θ)dφdθ=π
4

cos(θ_{A})^{2}+1

= π 4

(ξ_{A}_{i}·m)^{2}+1
.
(4.49)
We can now substitute this expression in Eq.4.47and get

Z

Dcos(m,ξ_{D},n)π
4

(ξ_{A}_{i}·m)^{2}+1

dm=

Z 2n+1

4π (ξ_{D}·m)^{2n}π
4

(ξ_{A}_{i}·m)^{2}+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(θ)^{2n}π
4

(sin(θ_{A})cos(φ)sin(θ) +cos(θ_{A})cos(θ))^{2}+1

sin(θ)dm= π

n(ξ_{D}·ξAi)^{2}+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_{1})^{2}+n+2

4n+6 +ρa2

π

n(ξ_{D}·ξA_{2})^{2}+n+2

4n+6 +ρa3

π

n(ξ_{D}·ξA_{3})^{2}+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(ω^{T}SAω) +tr(SA)(n+2)

4n+6 (4.53)

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)^{2}+1

dm=

Z 2π^{3}^{2}Γ(1+n)
Γ

3 2+n

1−(ξ_{D}·m)^{2}nπ
4

(ξ_{A}_{i}·m)^{2}+1

dm. (4.54)

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

2π Z

φ=0 π/2 Z

θ=0

π 4

(sin(θ_{A})cos(φ)sin(θ) +cos(θ_{A})cos(θ))^{2}+1

sin(θ)^{2n+1}dm=
π

3n+4−n(ξD·ξ_{A}_{i})^{2}

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(ω^{T}SAω)

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(ω^{0})Dcos(m,ξD,n)hω·midmdω, (4.57)

withω^{0}=−ω+2m(ω·m). We invert integrals and get
Z

Dcos(m,ξD,n) Z

A(ω)A(ω^{0})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(ω^{0})hω·midωdm=

### ∑

i,j∈{1,2,3}

Z

Dcos(m,ξD,n)aiaj Z

(ω·ξAi)^{2}(ω^{0}·ξA j)^{2}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)^{2}((−ω+2m(ω·m))·ξAi)^{2}hω·midω (4.60)

=

π/2 Z

φ=−π/2 Zπ

θ=0

(cos(φ_{Ai})cos(φ) +sin(φ_{Ai})sin(φ))^{2}(cos(φ_{Ai})cos(φ)−sin(φ_{Ai})sin(φ))^{2}cos(φ)sin(θ)^{6}dωdθ (4.61)

= 1 24π

15(m·ξ_{Ai})^{4}−10(m·ξAi)^{2}+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)^{4}−10(m·ξAi)^{2}+3

dm=π15n(n−1) (ξD·ξAi)^{4}−10n(n−2) (ξD·ξAi)^{2}+3n^{2}+7n+10

24n^{2}+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)]^{T}andAj= [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})^{2}((−ω+2m(ω·m))·ξ_{A j})^{2}hω·midω= π
24

15(m·ξ_{Ai})^{2}(m·ξ_{A j})^{2}+ (m·ξ_{Ai})^{2}+ (m·ξ_{A j})^{2}+1

. (4.64)

c

2018 The Author(s)

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

Z

Dcos(m,ξ_{D},n)π
24

15(m·ξAi)^{2}(m·ξA j)^{2}+ (m·ξAi)^{2}+ (m·ξA j)^{2}+1

dm=π15C1C2n(n−1) + (n^{2}+10n)(C1+C2) +n^{2}+5n+10
24n^{2}+96n+90

(4.65)
withC1= (ξ_{Ai}·ξD)^{2}andC2= (ξ_{A j}·ξD)^{2}.

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^{2}+96n+90

15n(n−1)
ξ^{T}DS_{A}ξD

2

−10n(n−2)ξ^{T}_{D}SASAξD+2n(n+10)

ξ^{T}DSAξDtr(S_{A})−ξ^{T}DSASAξD

+ (3n^{2}+7n+10)tr(SASA) + (n^{2}+5n+10)

tr(SA)^{2}−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})^{2}((−ω+2m(ω·m))·ξ_{Ai})^{2}hω·midω= π
24

15(m·ξ_{Ai})^{4}−10(m·ξ_{Ai})^{2}+3

(4.67) Then

Z

Dsin(m,ξD,n)π 24

15(m·ξAi)^{4}−10(m·ξAi)^{2}+3

dm=π(45n^{2}−45n)(ξD·ξAi)^{4}+ (−50n^{2}+10n)(ξD·ξAi)^{2}+29n^{2}+91n+80
192n^{2}+768n+720

(4.68) 4.7.4.2. Casei6= j

We have found previously (Sec.4.7.3.2) that Z

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

15(m·ξ_{Ai})^{2}(m·ξ_{A j})^{2}+ (m·ξ_{Ai})^{2}+ (m·ξ_{A j})^{2}+1

(4.69) We get

Z

D_{sin}(m,ξD,n)π
24

15(m·ξ_{Ai})^{2}(m·ξ_{A j})^{2}+ (m·ξ_{Ai})^{2}+ (m·ξ_{A j})^{2}+1

dm=π45C_{1}C_{2}n(n−1)−(C_{1}+C_{2})n(19n+25) +31n^{2}+105n+80
192n^{2}+768n+720

(4.70)
WithC1= (ξ_{A1}·ξD)^{2}andC2= (ξ_{A2}·ξD)^{2}

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^{2}+768n+720

45n(n−1)
ξ^{T}DSAξD

2

+10n(1−5n)ξ^{T}DSASAξD

−2n(19n+25)

ξ^{T}_{D}S_{A}ξDtr(S_{A})−ξ^{T}_{D}S_{A}S_{A}ξD

+ (29n^{2}+91n+80)tr(S_{A}S_{A}) + (31n^{2}+105n+80)

tr(S_{A})^{2}−tr(S_{A}S_{A})

(4.71)

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^{R}Dcos(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)^{2}
RDcos(m,ξD,n)(ωi·m)^{2}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)^{2}+1/4)
RDcos(m,ξD,n)((ω_{1}·m)^{2}+1/4)dm
because we have

|ω_{1}·m| ≤(ω_{1}·m)^{2}+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)((ω_{1}·m)^{2}+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|

(ω_{1}·m)^{2}+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)((ω_{1}·m)^{2}+1/4)
RDcos(m,ξD,n)((ω_{1}·m)^{2}+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(θ)^{2n}((cos(t_{i})cos(φ)sin(θ) +sin(t_{i})cos(θ))^{2}+1/4) (2n+1)(2n+3)
π(8ncos(ti)^{2}+2n+7).

c

2018 The Author(s)

Guillaume Loubet / Implementing the self-shadowing model

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

Z

D^{s}_{ω}_{i}(φ,θ)sin(θ)dφ=(2n+3)(2n+1)cos(θ)^{2n}sin(θ)
8 cos(ti)^{2}n+2n+7

cos(θ)^{2}

3 cos(ti)^{2}−1

−cos(ti)^{2}+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(θ)^{2n+1}
8 cos(ti)^{2}n+2n+7

(2n+1)

3 cos(t_{i})^{2}−1

cos(θ)^{2}+ (2n+3)

−cos(t_{i})^{2}+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})^{2}sin(θ)^{2}(cos(φ)sin(φ) +φ) +8 sin(θ)cos(θ)cos(t_{i})sin(t_{i})sin(φ) +4 cos(θ)^{2}cos(t_{i})^{2}φ+φ

2π(6 cos(θ)^{2}cos(t_{i})^{2}−2 cos(θ)^{2}−2 cos(t_{i})^{2}+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(φ)sin(θ) sin(φ)sin(θ)

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)^{2}+1/4)
4|ω_{h}·ωi|^{R}D(ω_{h})((ω_{h}·ωi)^{2}+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|^{R}D(ω_{h})((ω_{h}·ωi)^{2}+1/4)dm

((ω_{h}·ωi)^{2}+1/4) (4.76)

≤D^{s}_{ω}i(ωo)ραss(λ)AmaxA(ωi)^{R}D(ωh)((ωh·ωi)^{2}+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)^{2}+1/4)dm
RA(−ωi+2m(m·ωi))D(m)|m·ωi|dm.

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)^{2}+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^{R}D(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)^{R}D(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_{1}(ωo)

Mpdf_{2}(ωo)=ραss(λ)A(ω_{i})A(ωo)D(ω_{h})
4σs(ωi)

4^{R}D(m)|m·ωi|dm
D(ω_{h})

σs(ω_{i})

ραss(λ)AmaxA(ωi)^{R}D(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)^{R}D(m)hm·ωoidm

RA^{R}D(m)hm·ωoidmdω−^{R}A^{R}AD(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)^{R}D(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}2n+1
4π

c

2018 The Author(s)

Guillaume Loubet / Implementing the self-shadowing model

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}2n+1

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

2 cos(θ)^{2n}sin(θ).

The cumulative distribution function writes

x Z

0

2n+1

2 cos(θ)^{2n}sin(θ)dθ= 1
2

1−cos(x)^{2n+1}
.

This function can be inverted:

U=1 2

1−cos(x)^{2n+1}

⇒cos(x)^{2n+1}=1−2U,
IfU≥1/2:

x=arccos

(1−2U)^{2n+1}^{1}

. (4.89)

IfU≤1/2:

x=π−arccos

(2U−1)^{2n+1}^{1}

. (4.90)

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

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

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,^{π}_{2})(θ).

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)^{2}.

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

θ=arcsin√ U

.