A Practical Stochastic Algorithm for Rendering Mirror-Like Flakes
Asen Atanasov Vladimir Koylazov
Chaos Group∗

Figure 1: A) A metal plate with our flakes material. B) A Beetle toy with diffuse base layer, stochastic flakes middle layer, and reflective coat.
Flakes use triplanar mapping and have colors with random hue in HSL color space. C) Sparkling snow with sub-surface scattering base layerand stochastic flakes coat material. The model parameters (N, αGGX , γ) for A), B), and C) are (16 × 107 , 0.04, 1◦ ), (8 × 107 , 0.09, 0.3◦ ),and (108 , 0.25, 0.4◦ ) respectively. B) and C) have ray-traced depth of field and slight post-process lens effects.

Abstract
Materials, such as snow, sand, metallic paints, rough plastics, andmetals, often exhibit small-scale phenomena observed as brightsparkling or glittering surface features. These features become
more pronounced under narrow-angle illumination and vary basedon the orientation of the surface with respect to the viewer andlight sources. Microfacet-based surface models, composed of alarge finite number of microscopic mirror-like flakes, can mimic
this effect. An associated microfacet BRDF and a memory-efficientstochastic algorithm are explored in [Jakob et al. 2014]. We present
a new stochastic algorithm that inherits the good properties of theoriginal algorithm, but does not require any precomputation; implements optimal importance sampling which is extended to efficiently
sample wide and heavy-tailed microfacet distributions (i.e. GGX),and offers better overall performance. In addition, a triplanar mapping technique is employed to handle geometry without texture coordinates. The algorithm is both practical and easier to implement.
Keywords: Stochastic algorithm, microfacet BRDF, importancesampling, sparkling snow, metallic paintConcepts: •Computing methodologies → Rendering; Reflectance modeling;

1Introduction

Microfacet theory idealizes surfaces as composed of infinitely many
microscopic facets, described by two statistical measures - themicrofacet distribution D(x, h) and the shadowing-masking function G(i, o, h). The related microfacet BRDF is fr (x, i, o) =1
F (i · h)D(x, h)G(i, o, h), where F is the Fresnel term,
4(i·n)(o·n)
x - the shading point, i - direction from which light is incident, o ∗ e-mail: {asen.atanasov,vlado}@chaosgroup.com
Permission to make digital or hard copies of part or all of this work forpersonal or classroom use is granted without fee provided that copies are
not made or distributed for profit or commercial advantage and that copies
bear this notice and the full citation on the first page. Copyrights for thirdparty components of this work must be honored. For all other uses, contact
the owner/author(s). c 2016 Copyright held by the owner/author(s).
SIGGRAPH ’16, July 24-28, 2016, Anaheim, CA,
ISBN: 978-1-4503-4282-7/16/07
DOI: http://dx.doi.org/10.1145/2897839.2927391

i+o, and n - the surface
direction in which light is reflected, h = ||i+o||
normal [Walter et al. 2007]. Converged renders of such surfaces
appear perfectly smooth. However, the visible microstructure, especially observed in motion or change of illumination, gives morerealism. Usually, the rendering of small specular surface features isachieved through normals maps, but as the features become muchsmaller than a pixel, resolving the correct pixel value using Monte
Carlo sampling becomes inefficient.

2Previous work

The memory-efficient rendering of large finite number of flakes isinvestigated in [Jakob et al. 2014]. N flakes are uniformly distributed in the unit texture space and their normals follow a microfacet distribution on the unit hemisphere. The multiscale BRDF
is defined as the microfacet BRDF, averaged over a finite surface
area A, containing the shading point x and a finite solid angle Ωaround the incident direction i. In the implementation, A is a parallelogram approximation of the pixel footprint, defined by the raydifferentials [Idehy 1999], and Ω is a cone of radius γ, centered
at i. The continuous microfacet distribution D is substituted by a4b where Db is the fraction of(i · h)D,discrete counterpart a(A)σ(Ω)
flakes which are contained in A and reflect o in Ω. The expresb is an approximation to D̃ = 1 P i · mj , where mjsion (i · h)DN
are the normals of this fraction of flakes. A hierarchical traversalof the texture-direction space, starting from the unit texture squareb A key advanand the unit hemisphere, is performed to compute D.
tage of the algorithm is that the flakes are not stored in memory, but
their counts are reproduced by a deterministically seeded stochastic
process during traversal. The result converges to a standard microfacet BRDF as N increases. This method has practical issues. To
build the data structure efficiently on average, it relies on statistics
collected during a full animation precomputation step. The importance sampling strategy does not sample the discrete integrand, but
samples from a different smooth distribution. Thus, it practically
fails to sample wide and heavy-tailed microfacet distributions.

3Our approach

We present a new stochastic algorithm, based on [Jakob et al.
2014], which inherits the good properties and solves the issues
listed above. For a given shading point with a footprint A, we per-form one texture space query to locate the number of flakes in it
during which we explicitly generate the normals of the flakes andcollect data to calculate the lighting. A comparison between thetwo algorithms is given in Table 1.

Table 1: A comparison between Jakob et al. and our method for
N = 2 × 107 and distribution width α. Variance-based image
sampler with a fixed threshold is used to render 200 × 300 regionof the metal plate on Intel Core i7-980X, 3.33GHz, 6-core machine.

Texture query. We query the texture space with the footprint parallelogram A, using a quad-tree, stack-based, depth-first traversal. The root of the tree is the unit texture square which has
the flakes count N and seed s used to initialize a pseudorandom
generator for this query. Each square node is split into 4 equal
nodes. Node-footprint intersections are performed to filter out thenodes which do not overlap with the footprint. On each split operation we distribute the node flakes count to its children using
normal approximation to multinomial distribution with probabilities p = ( 41 , 14 , 14 , 41 ) [Jakob et al. 2014]. Our query generator
combines four xorshift-based generators which share a 64-bit state.
During the traversal push operation, we use a different generator foreach child node and the state is stored on the stack. On the pop operation, the generator state is restored. A node is a leaf if it has lessthan K flakes, or the stack depth T is reached. For our implementation N < 231 and (K, T ) = (16, 15) proved adequate. If a leaf
is entirely inside the footprint, then all of its flakes are processed.
Otherwise we generate random positions of the flakes inside theleaf and process these flakes which belong to the footprint.

αJakob et al. (Beckmann, γ = 1◦ )Jakob et al. (Beckmann, γ = 5◦ )Our method (Beckmann, γ = 1◦ )Our method (Beckmann, γ = 5◦ )Our method (GGX, γ = 1◦ )Our method (GGX, γ = 5◦ )Reflection cache. During the query, we accumulate n footprint
flakes and generate their random normals mj , according to the desired microfacet distribution, using the query generator. Then o isreflected from each normal mj to get per-flake directions rj . These
directions represent n cones of radius γ from which light can be reflected along the outgoing direction o. Furthermore, we define perflake weights wj = mj · rj . We call the set of all n pairs (rj , wj )reflection cache. For the microfacet
P distribution, we compute anaccurate approximation D̃ ≈ N1i·rj ≥cos γ wj by checking how
many cached cones contain the incident light direction i.

31 10 2oAmjrj

1
K=4AFigure 2: Four non-empty leaves (red) overlapping the footprint(left). Partially overlapping leaves generate flakes positions andcheck if they belong to the footprint (middle). Generated normals
mj for the three flakes and reflection cache directions rj (right).
Optimal importance sampling. Along with the reflection cache,
we construct a cumulative table with n + 1 elements c0 = 0,
cj = cj−1 + wj , j = 1, n. We pick a random cache directionr, proportional to wj , by generating a uniform random numberξ ∈ (0, cn ) and binary searching the cumulative table with it. Thenwe generate the BRDF direction i as a uniform random direction ina cone of radius γ around r. The corresponding geometric prob1
ability is π(1−cos, where π(1 − cos γ) is the solid angle of theγ)kcone. Finally, the probability of generating i is p(i) = π(1−cos,γ)cnP
k = i·rj ≥cos γ wj is the sum of the weights of the cached cones
which contain i and could have generated it proportionally to theirweights. We compute k as a by-product of the BRDF evaluation.
Distance blending. As the number of flakes within the footprint increases, the material appearance conforms to the smooth microfacet
model and the calculations of our algorithm increase. Therefore,0.01
19s
18s
4s
6s
5s
6s

0.1
631s
49s
8s
8s
9s
9s

0.4
1945s
170s
10s
10s
9s
9s

we implement blending based on the expected number of flakes
per pixel footprint E[n] = N · a(A). We define blending range[bmin , bmax ]. For all footprints with expected number of flakes
E[n] ≤ bmin we compute only the discrete microfacet distribution;
for E[n] ≥ bmax we compute only the smooth microfacet distribution; for the case of bmin ≤ E[n] ≤ bmax we compute both andlinearly interpolate between them, using E[n] as an interpolation
parameter. Our default blending range is [500, 2000].
Colored flakes. Optionally, we generate a color table with a fixed
number of entries which is either spanning the hue parameter in
HSL space for a fixed saturation and lightness, or it is sampled froma texture. Per-flake color is assigned by generating a random index
in this table. We use the same blending strategy, but we multiply
the smooth BRDF by the average color of the table. For importancesampling we pick random cache direction proportional to the product of the weights of the flakes wj and their color intensities. In ourimplementation the table size is 64.
Triplanar mapping. In practice, texture coordinates are not always
available or have bad quality. Alternatively, our method works withtriplanar mapping. We transform the surface normal and the footprint, defined in world coordinates, into the local coordinates of theshaded object. Then we project the footprint along the axis corresponding to the maximal absolute component of the local normal.
Then we scale the resulting 2D footprint depending on the shadedobject size. We find the square with integer coordinates and size
1, to which the footprint belongs, and perform the texture query in
it. In order to avoid tiling, we generate the seed s of the root node
based on the square coordinates. If the footprint is at the boundaryof such squares, we query all of them.

Acknowledgements
We thank Stephen Parker, Rossen Dimov, Yordan Madzhunkov,
Simeon Balabanov, Kalina Panteleeva, Yavor Stoikov, Emanuele
Lecchi, Christopher Nichols and Rhys Dippie.

References
I DEHY, H. 1999. Tracing ray differentials. In Proceedings ofSIGGRAPH 99, ACM Press/Addison-Wesley Publishing Co.
JAKOB , W., H AAN , M., YAN , L.-Q., L AWRENCE , J., R A MAMOORTHI , R., AND M ARSCHNER , S. 2014. Discrete
stochastic microfacet models. ACM Transactions on Graphics(Proceedings of SIGGRAPH 2014) 33, 4.
WALTER , B., M ARSCHNER , S., AND T ORRANCE , K. 2007.
Microfacet models for refraction through rough surafaces. In
Proceedings of the 18th Eurographics Conference on RenderingTechniques, 195–206.