Vis Comput
DOI 10.1007/s00371-017-1352-2

ORIGINAL ARTICLE

Approximations for the distribution of microflake normals
Nelson Max1,2· Tom Duff1 · Ben Mildenhall3 · Yajie Yan4

© Springer-Verlag Berlin Heidelberg 2017

Abstract Scenes in computer animation can have extreme
complexity, especially when high resolution objects areplaced in the distance and occupy only a few pixels. A usefultechnique for level of detail in these cases is to use a sparsevoxel octree containing both hard surfaces and a participating
medium consisting of microflakes. In this paper, we discuss
three different methods for approximating the distributionof normals of the microflakes, which is needed to computeextinction, inscattering of attenuated direct illumination, andmultiple scattering in the participating medium. Specifically,
we consider (a) k means approximation with k weighted representatives, (b) expansion in spherical harmonics, and (c) the
distribution of the normals of a specific ellipsoid. We compare their image quality, data size, and computation time.
Keywords Microflake · Sparse voxel octree · Volume
rendering · Distribution of normals

Electronic supplementary material The online version of thisarticle (doi:10.1007/s00371-017-1352-2) contains supplementarymaterial, which is available to authorized users.

B Nelson Max

max@cs.ucdavis.edu
Tom Duff
td@pixar.com
Ben Mildenhall
bmild@berkeley.edu
Yajie Yan
yajieyan@wustl.edu

1

Pixar Animation Studios, Emeryville, CA, USA

2

Present Address: University of California, Davis, CA, USA

3

University of California, Berkeley, CA, USA

4

Washington University in St. Louis, St. Louis, MO, USA

1 Introduction
Scene models for computer animation such as in Fig. 1 often
have extreme complexity if the original detailed models areused for objects in the distance. One method of dealing withthis data explosion is to have multiple levels of detail, sothat simplified models can be used for more distant objects.
Another method is to summarize into a volume data structure
the aspects of the model that affect rendering. We are using
the latter method, with a sparse voxel octree (SVO), of thesort described by [9]. The polygons (or micropolygons, see[2]) are scan converted into the leaves of the octree structure,
averaging the plane equation, plane normal variance, specular and diffuse color, and other shading quantities such as
surface roughness. These quantities are then averaged further
into the higher octree levels. When marching along a viewing
ray through the octree, cells are selected at the octree levelappropriate to the ray differential (randomly jittered to breakup any visible transitions in shading from different SVO levels), and the ray is intersected with the cell’s stored averagesurface plane, as clipped to the cell volume. The first pointintersected is shaded by combining the plane normal variance into the surface roughness. We refer to this as the “hardsurface” method.
One of our goals was to model distant vegetation, as inFig. 1. The hard surface method does not work well for this,because we are intersecting the ray with the average plane,
while the octree cell should actually represent many small
plant leaves inside its volume. In addition, even a single tinyleaf in the volume may cause an intersection with its extended
plane, making the apparent size of the leaf bloom larger.
To solve these problems, we treat the leaf surfaces in an
octree cell as a volume density of microflakes (infinitesimalpieces of flat surface) and use volume rendering methods tocomposite the inscattering and extinction along the ray as123

N. Max et al.

Fig. 1 A wide angle view of a forest of 218 trees of various species, path traced using the spherical harmonic microflake normal distribution. The
sparse voxel octree used 2,434,796,669 bytes, while a single pine tree polygon representation needed 1,816,829,593 bytes

it steps through the octree cells. Both the inscattering andthe extinction depend on the distribution of the microflake
surface normals. This paper compares several methods ofrepresenting this distribution by approximations specified bya small number of parameters. Ideally, such approximations
should give:
1. A compact representation,
2. Efficient averaging up the octree hierarchy,
3. Efficient computation of the extinction coefficient,
4. Efficient computation of the inscattering,
5. Efficient sampling of the scattered ray direction for global
illumination path tracing, and
6. Final renderings that well approximate converged ray
traced or path traced images of the original geometry.
The contributions of this paper are:
– Two new methods to approximate a distribution of
normals, by k means representatives and by sphericalharmonics (SH), and– Their comparison with each other and the SGGX distribution of Heitz et al. [8], which uses the distributionof normals of a specific ellipsoid.
Section 2 discusses related work, and Sect. 3 describes
the three approximations and how their fitting parameters
are computed. Section 4 discusses how the extinction andinscattering are computed from the distribution of microflake
normals. Section 5 shows how sample directions for pathtracing are computed, Sect. 6 discusses optimizations fordouble-sided microflakes, and Sect. 7 describes compressionof the data. Section 8 presents image accuracy, data size, andtiming results, which are discussed in Sect. 9. The appendix
has mathematical details for the SH approximation.

Octrees have long been used to describe sparse volumedata and speed up ray tracing. Neyret [23] used a hierarchicaloctree volume with the reflectance in each cell represented
by the distribution of normals of an ellipsoid. Heitz et al.
[8] did something similar in the context of the SVO of [9].
We have also used a version of this SVO, enhanced to store
the various approximations to the distribution of microflake
normals, and with some new compression mechanisms.
The basic principles of volume rendering participating
media, accounting for the optical effects of extinction andinscattering, are explained in [20]. Normally, the extinction
is isotropic, independent of the viewing direction. The scattering phase function, giving the distribution of scattering
directions ωo from an incoming direction ωi , is also normally isotropic, in the sense that it depends only on the angle
between ωi and ωo . However, in the case of an anisotropic
distribution of microflake normals, these isotropies are notpresent, and the full dependence of the extinction and inscattering on both ωo and ωi must be accounted for. This involves
the distribution of visible normals, which accounts for bothmasking and foreshortening for bumpy surfaces [6], and forforeshortening only in the case of microflakes, as discussed in[8], and briefly in Sect. 4 below. An early computer graphics
paper considering the optical effects of a microflake normal
distribution is [21], which simulated multiple scattering ina tree canopy. More recently, Jacob et al. [13] discussed theradiation transport equation for a general anisotropic distribution of microflake normals, and also the corresponding
diffusion approximation. Zhao et al. [32,33] modeled cloth asa collection of fibers extracted from micro-CT and modeled
the microfacet normal distribution of a fiber as a gaussian in
the dot product of the microflake normal and the fiber axis.

2 Related work
The fundamental paper on level of detail was Crow [4], which
first suggested designing several different polygonal models
for the same object. This idea was refined by Hoppe to create
the simplified models automatically [11], and continuouslytransition between them [12].

123

3 Representing the distribution of microflake
normals
Let g(x, ωn ) be the local probability distribution of the
microflake normals with direction ωn at position x, let ρ(x)be the density of total microflake area per unit volume, andApproximations for the distribution of microflake normals

let h(x, ωn ) = ρ(x)g(x, ωn ). Then, h(x, ωn )dωn is the total
area per unit volume of microflakes with normals within asolid angle dωn about ωn . We will usually mean this scaledh(x, ωn ) when we say the distribution of microflake normals
below, even though this is not a probability distribution.
We will compare three different methods of approximating
this distribution of microflake normals, by:(a) A finite number k of weighted δ-functions of unit direction,(b) An SH expansion of fixed maximum order L, or(c) the distribution of normals of an origin-centered ellipsoid.
If the degrees of freedom are all expressed as floats, thenmethod (a) requires 4k floats, method (b) requires (L + 1)2
floats, and method (c) requires only six floats at each SVO
node. Later we will discuss compression of these data.
When scan converting polygons into the SVO, we recursively split the polygons by the slicing planes defining the
SVO volume cells, until the fragments are the size of theleaf cells. For the hard surface contribution, we compute thearea inside each intersected leaf cell and add that area andalso the area-weighted diffuse color, specular color, unit normal, and other shading parameters like surface roughness, tothe accumulated values for that leaf cell. After all the polygons are scan converted, we add these weighted quantities
to the parent cells up the octree hierarchy and then divide
all these summed quantities by the sum of the weights. We
also compute the variance of the normals. An SVO cell cancontain both hard surfaces and microflakes, and the colorsof the microflakes are averaged separately into the SVO asabove, so that the plant leaves can have a different color thanthe hard surface trunks and branches. The calculations to fit
each of the three microflake normal distribution approximation methods into this SVO framework are described next.
The k means algorithm for method (a) converges to a local
minimum of its cost function, which depends on the initialcluster center normal representatives chosen. To get a good
spread in the initial representatives, a first pass through thepolygon fragments finds up to k representative normals peroctree cell. We take a new area-weighted representative normal whenever there are less than k and the current polygon
fragment normal differs from all existing representative normals by more than a threshold angle. Otherwise we add thearea-weighted fragment normal to the representative to which
is closest in direction. Then, in subsequent passes through thepolygons, we reassign each polygon fragment’s contribution
to the representative to which it is closest and compute newweighted average representative directions, thus implementing a weighted k means algorithm. When all polygons arehandled, we use the same procedure to go up the octree hierarchy, assigning child cell weighted representatives to parent

cell representatives, and computing the area-weighted average of the normals in each representative’s group. Currently,
we use only two passes through all the polygons in computing the k means for the leaf cells of the SVO, since suchpasses are expensive, but more passes to propagate the leaf
representatives into the k means of their ancestor cells up
the hierarchy, since this does not need the full polygon data.
The resulting approximation is an area-weighted sum of theδ-functions for the representative unit normals in each SVO
cell. We also save the variance of the normals in the cluster
assigned to each representative, estimated from the lengthof the weighted average normal using the method of [29],which does not require a separate pass through the polygon
data. This variance per representative is the fourth of the 4k
floats mentioned above. For greater accuracy, we could store
the covariance matrix of the distribution of the microfacet
normals assigned to each representative normal, as in [25],as done for method (c) below, but this would take 6 extra
floats per representative normal, instead of one.
Method (b) fits the distribution of normals using the real
spherical harmonic basis functions Ylm (θ, φ).
Since the SH basis functions are orthonormal, the best L 2
fit to a function h(θ, φ), of the formh(θ, φ) ≈L lh lm Ylm (θ, φ)(1)

l=0 m=−lhas coefficients given byh lm =

 2π  π
0

Ylm (θ, φ)h(θ, φ) sinθ dθ dφ.

(2)

0

So for an area-weighted δ-function representing the contribution of a polygon fragment with a constant unit normal toan SH basis function coefficient, this integral is just the area
times that basis function evaluated at that normal. These contributions are added into the leaf cells and then into parent
cells up the octree hierarchy. There is no need to divide by
the sum of the weights, because the goal is to approximate
an area-weighted density h(x, ωn ). However, since this is anarea per unit volume, the accumulated sum must be divided
by the volume of the octree cell. [This is also the case formethod (a).] We have experimented with spherical harmonicorders L up to 4.
For method (c), Heitz et al. [8] define the SGGX distribution of squared projected areas of an ellipsoid as a quadratic
form in the unit projection direction vector. They define it in
the coordinate system of the unit eigenvector directions (theellipsoid’s major, semi-major, and minor axis directions) bya diagonal matrix of the squared projected areas in these
directions and use the rotation matrix formed from these unit
eigenvector directions to transform this quadratic form into
the world coordinate system. They prove that the square root123

N. Max et al.

of this quadratic form evaluated at any other direction gives
the projected area of the ellipsoid in that direction. For apolygon fragment from a smooth surface, they use a “surfacelike” ellipsoid with unit projected area in the surface normal
direction, and very small projected area in the perpendicular
directions, while for a fragment of a curve (for fur, hair, or
pine needles), they use a “fiber-like” ellipsoid, with a small
projected area along the curve tangent direction, and unit
projected areas along perpendicular directions.
For the images in that paper, they just accumulate thesix independent elements of the positive definite 3 by 3
symmetric matrix for the quadratic form into the leaf cells
and then into the octree hierarchy, but they admit that thisdoes not properly accumulate the distributions of normalsof child cells into parents. So we use the two pass method
they describe in their Sect. 4.3, subsection “Parameter Estimation from Arbitrary Distributions.” In the first pass, we
accumulate the covariance matrices for the distributions. The
covariance matrix contribution to an octree leaf cell of aclipped polygon fragment with unit normal column vector ωis simply the fragment’s area times the outer product matrixωωT . These contributions are summed in each leaf cell andthen up the hierarchy. At the end of this pass, we computeand save the eigenvectors of the resulting symmetric covariance matrix in each octree cell. In a second pass throughthe polygons, we compute the projected areas of the polygon
fragments in these eigenvector directions, not just for the leaf
cells, but for all the octree cells encountered in the recursive
polygon slicing. Then, we use the square of these summed
projected areas to get the diagonal matrix for the quadratic
form in the eigenvector coordinate system and transform it
into the world coordinate system using the matrix of unit
eigenvectors.

4 Computing extinction and inscattering
The microflake area per unit volume near a position x, as
projected in the viewing direction ωo , is
σt (x, ωo ) =h(x, ωn )ωn , ωo dωn ,(3)where is the unit sphere, and the foreshortening projection
factor −, − represents the nonnegatively clamped dot product if one-sided flakes are being used, or the absolute valueof the dot product if double-sided flakes are being used, andωo is actually the negative of the viewing ray direction, as
in the usual convention for representing BRDFs. According
to the reasoning in [20], this σt (x, ωo ) is the extinction coefficient, or volumetric attenuation coefficient, expressing thefraction of flux intercepted by the flakes per (infinitesimal)
unit length. Figure 2 (left) shows a slab of front face area 1 and123

infinitesimal thickness ds, and thus with volume ds, a sample of microflakes from the distribution of normals h(x, ωn ),and their projections on the face of the slab closest to theviewer. The slab thickness ds is assumed to be so small thatthe chance that the microflake projections overlap each otherapproaches zero. Therefore, sum of their projected areas is
the slab’s opacity, expressed by multiplying Eq. (3) by theslab volume ds to convert area per unit volume to area in theslab. So the slab’s transparency is 1 − σt (x, ω0 )ds.
To compute the inscattering, let f (ωi , ωo ; ωn ) be the
BRDF of a microflake with normal ωn . Accounting for boththe foreshortening in the incident direction ωi (for irradiance)
and in the viewing direction ωo (for viewing visibility), theinscattering per unit length along the viewing ray isS(x, ωo ) =h(ωn ) f (ωi , ωo ; ωn )ωn , ωi ωn , ωo L(x, ωi )dωn dωi ,(4)where L(x, ωi ) is the incoming radiance at x coming fromdirection ωi . See the right half of Fig. 2, where the colors
represent the microflakes’ reflected radiance in direction ωo .
The quantities σt (x, ωo ) and S(x, ωo ) are used to compute the volume rendering integral along the viewing ray indirection −ωo , by integrating this inscattering contribution
along the ray, as attenuated by the transparency between eachslab and the viewpoint E:
L(E, ωo ) = D   sexp − (σt (x(t), ωo )dt S(x(s), ωo )ds,
0

0

(5)
where x(t) = E − tωo is the arclength parametrization of
the ray leaving the viewpoint E in direction −ωo , and D
is the parameter at which this ray hits a hard surface orexits the SVO data volume. See [20] for a detailed derivation. As shown by the equations above, the extinction andinscattering depend on the distribution of visible normals,
due to the foreshortening factor ωn , ωi . However, the fullgeometric attenuation factor accounting for shadowing andmasking when considering surface microfacets is not present
here, because these effects are instead included in the volume
attenuation calculations in Eq. (5) above, and in computing
similarly how much light from a distant source reaches x andcontributes to L(x, ωi ) in Eq. (4).
As is usual in volume rendering, we numerically computethe integral in equation (5) by marching along the viewing
ray through the octree, determining the successive segments
in which it intersects cells at the level locally appropriate
to its ray differential, and incrementally accumulating thetransparency

Approximations for the distribution of microflake normals
Fig. 2 Left extinction in a slabof area 1 and infinitesimal
thickness ds, Right inscatteringin that slabs

= S(ωo )exp(−σt (ωo )s)ds0= S(ωo )(1 − exp(−σt (ωo ) s))/σt (ωo )= (S(ωo )/σt (ωo ))(1 − Tseg ).

Fig. 3 Mutually occluding microflakes in a slab of finite thicknessT (s) = exp −

 s

s


(σt (x(t), ωo )dt(6)0

and also the ray segment’s contribution to the radiance arriving at E in Eq. (5).
We assume that h(x, ωn ), and therefore σ (x, ωo ) andS(x, ωo ), are constant in each cell, so we leave the location
x out of the equations below. If the ray segment has lengths, the update for T should not just be T = T × (1 − σt s),
accounting for the unoccluded fraction of the slab face in
Fig. 2, because s is not infinitesimal, and the projected
microflakes can overlap, as shown in Fig. 3. It is more accurate to compute the segment transparency integral Tseg = exp −s(σt (ωo )dt

0= exp(−σt (ωo ) s)(7)

which accounts for this overlap. Then, we use the update T =T × Tseg . Similarly, instead of incrementing L(E, ωo ) by
T ×S(x, ωo )× s as is sometimes done in volume rendering,
it is more accurate to compute the inscattering integral fromthe segment analytically, asL seg =
0

s


exp −

 s
0


(σt (ωo )dt S(ωo )ds(8)

The expression in the final line just above can be interpreted
in Fig. 3. Without accounting for overlap, S(ωo ) s is thesummed radiance times projected area of the microflakes,
and σt (ωo ) s is the summed projected area. So the fraction
S(ωo )/σt (ωo ) is the projected area-weighted average radiance of the microflakes. Since the depth order for occlusion
is random, this is also the average radiance of the colored
segments shown in Fig. 3 on the slab face. The area not covered by these segments is Tseg , so the area that is covered is
1 − Tseg . Thus, the inscattered radiance L seg in equation (8)is the average radiance times the area covered. The updateto L(E, ωo ) in equation (5) to include this segment’s contribution is L(E, ωo ) = L(E, ωo ) + T × L seg where T is
the accumulated transparency of all the previous segments,
before being updated to include the current one.
Now we discuss how to compute σt (x, ωo ) and S(x, ωo )for each of the three methods for approximating h(x, ω) in acell. For method (a), the distribution of normals is a sum of
weighted δ-functions, with weights c j , in directions ω j :kh(x, ω) =c j δ(ω − ω j ).
(9)
j=1

Substituting into equation (3), the extinction iskσt (x, ωo ) =c j ω j , ωo .

(10)

j=1

For a perfectly diffuse surface with albedo α and thus
diffuse BRDF α/π , and a small light source in direction ωiwhose radiance times subtended solid angle at x is I , equation(4) for the inscattering gives
Iα 
c j ω j , ωo ω j , ωi .
πkS(x, ωo ) =

(11)

j=1

123

N. Max et al.

Specular reflection can be handled similarly, using the variance of the polygon fragment normals clustered with each of
the k representative normals, calculated when constructing
the octree, to determine the width and maximum radiance ofeach specular peak.
For method (b), the integral for the extinction coefficientin equation (3) of the SH expansion for the microflake normal density times the clamped cosine can be computed as aweighted sum of the spherical harmonic basis functions evaluated at the viewing direction, according to Sect. 4 of [26][See Eq. (15) in the “Appendix” below]. If the two clamped
cosines in the integral for the inscattering in Eq. (4) are both
expanded in spherical harmonics as in [26], this integral canbe expressed as a weighed sum of triple product integralsof spherical harmonic basis functions [See the “Appendix”,
equation (16)]. Most of these triple product integrals are zerofor the higher orders, and we precomputed the rest and stored
them in a sparse array. Jacob et al. [13] also discussed theabove method of computing the extinction coefficient froma spherical harmonic representation, but not how to computethe inscattering integral.
For microfacet-generated mirror reflection, we need toevaluate the spherical harmonic visible microfacet normal
density approximation ω, ωo h(ω) at the half angle direction ω = ωh = (ωi + ωo )/|ωi + ωo | that specifies the
microfacet normal which reflects light from direction ωiinto direction ωo , as in [30] and [3]. We must multiply by
the change of variables factor dωh /dωo = 1/(4ωh , ωi ) as
given in [27] and explained in the introduction to Chap. 7 of
[16] and briefly in [31]. Since ωh , ωo  = ωh , ωi  these two
factors cancel.
For method (c), we use the routines in [7]. For the projected microflake area in the viewing direction in equation(3), we take the square root of the quadratic form of squared
projected areas, evaluated at the viewing direction. The perfect mirror microflake specular reflection can be found fromthe microflake normal distribution evaluated at the half angle
direction ωh as in method (b) above. This distribution is given
by equations (24), (11), and (12) of [8]. The diffuse reflection is more difficult, because, unlike methods (a) and (b), an
analytic expression for the integral in Eq. (4) is not known,even for a point or directional light source. So as proposed in[8], we instead sample the distribution of visible normals (seethe next section) and use the diffuse BRDF for a flat surface
with the resulting normal. This technique introduces noise,
however, and is only useful in a context where multiple ray
samples per pixel are already being used for other purposes
like anti-aliasing, depth of field, or global illumination. We
may also need to use it for the spherical harmonics method(b), if there are double-sided microflakes. (See Sect. 6 below.)123

5 Computing sample directions for global
illumination
To efficiently compute inscattered illumination from a wholeenvironment sphere in equation (4), or the next bounce
direction in path tracing methods for global illumination byMonte Carlo integration, once the Monte Carlo volume rendering algorithm determines that volume scattering shouldtake place (see [15]), we need to importance-sample the
distribution of scattering directions ωi for a viewing ray orcontinuing path direction −ωo . This involves first samplinga normal from the distribution of visible microfacet normals
ωo , ωn h(x, ωn ), or rather the normalized probability distribution proportional to it, and then sampling the BRDF fora flat microflake surface with this normal.
For method (a), we generate a random number uniformly
distributed in the interval [0, 1], use it with the inverse of
the cumulative distribution of the visible normal probabilitiesω j , ωo c j / ω j , ωo c j to chose surface j with probability
proportional to ω j , ωo c j , and then sample the BRDF for asurface with normal ω j .
For method (b), we show in equation (18) of the “Appendix”
that the importance of the incoming radiance from the next
diffuse bounce direction also has an SH representation, withcoefficients computed from ωo and the SH coefficients of the
microflake normal distribution. Thus, we can use the algorithm of [14] for sampling spherical harmonics distributions.
It divides the (φ, θ ) range of the unit sphere iteratively into
quadrants as in a quadtree, integrates the distribution overeach quadrant (using precomputed tables), and uses a random number for the inverse cumulative distribution method
to choose a quadrant at each stage of the iteration.
One problem with this sampling algorithm is that the best
spherical harmonic fit of a nonnegative function on the unit
sphere can have some negative values, due to ringing. In
our implementation, quadrants with negative integrals are
given probability zero, and the probabilities of the remaining quadrants are normalized to add up to one. However,when negative lobes of the SH function cause distortions in
the distribution of our samples by inappropriately decreasing
the probability of sampling the positive lobes in a quadrant,
we resort to rejection sampling, precomputing and savingin each octree cell the maximum of its spherical harmonic
fitting function. The rejection test rejects samples with negative function values and produces directions with probability
proportional to positive function values.
For method (c), we use the sampling algorithm in Sect. 1.6of [7].

Approximations for the distribution of microflake normals

6 Double-sided microflakes
In our application, the vegetation leaf polygons are shadedon both sides with the same color, as is the case for the
microflakes in [13]. For method (a), we can take advantageof this and use a smaller number k of representatives of the
microflake normals, by comparing both a new normal n andits negative −n with all the representatives when choosing
the closest. Then, when shading, we can use either the representative normal, or its negative, whichever is front-facing.
For method (b), we can conceptually use the δ-functions
with both n and −n when computing the spherical harmonic
representation, as in Sect. 3. The basis functions Ylm (n) areeven functions of the unit vector n, i.e., Ylm (−n) = Ylm (n),when l is even, and odd functions of n, i.e., Ylm (−n) =−Ylm (n), when l is odd. Thus, when summing the contributions from n and −n, all terms with odd l disappear, and only(L/2 + 1)(2L/2 + 1) terms remain, which for large L isabout half of the total (L +1)2 . For order L = 4, there remain
only 15 terms instead of 25. This saves both space and computation time. However, when Eq. (16) in the “Appendix” isapplied to the shading of a double-sided smooth surface we
get an effect like double-sided lighting (as if from direction−ωi as well as ωi , see Fig. 10), because the reversed normal also contributes to the shading. To avoid this, we canshade by sampling the visible normal distribution, reversing
the resulting normal if it has a negative dot product with ω0 ,and then applying the BRDF for a surface with this normal.
This normal sampling introduces noise, as with the shading for the ellipsoid method. An alternative is to keep all(L + 1)2 terms, and apply Eq. (16) twice, using the singlesided clamped version of the ωn , ωi  and ωn , ωo  factors,
and for the second time reversing both the viewing and lighting directions, which is equivalent to reversing the normal.
On the other hand, the double-sided lighting may actually beappropriate for vegetation leaves, if diffuse transmission isconsidered, as well as diffuse reflection.
Method (c) naturally represents double-sided microflakes,
since the ellipsoids are centrally symmetric and reflect fromall sides.

7 Compression
In order to reduce the size of our SVO, we compress the data
inside each cell after octree construction and dynamicallydecompress the requested values as rendering occurs. We use
the Best Fit Normal [19] method to compress normals, and adictionary-based method to compress other related floating
point numbers.
Best Fit Normal compresses a 3-float normal (12 bytes) to
3 bytes, achieving a compression ratio of 4. Essentially, thismethod assigns a normal to the best fit voxel’s index within avoxel space of size 2563 . The trick is how to find the best fitting voxel for each unit normal. To do this, instead of treating
the normal as a unit vector, we treat it as a ray with arbitrary
length. We start this ray at the center of the 2563 voxel space
and then walk along it. For every voxel intersected by theray, we compute the distance from that voxel center to the
ray and assign the index of the closest voxel to the normal
as its compressed form. The decompression of such a form
amounts to a subtraction and a normalization.
Our dictionary-based compression scheme takes in a set offloats, builds a dictionary, and assigns each float to its closest
entry in the dictionary. To build the dictionary, we use anapproximate k means method [18] for fast clustering, andextract N entries as specified by the user. For a single float,
the compression ratio is then 32/logN  bits. Decompression
is merely a lookup in the dictionary. We use this method forcolor components, and normal variance.
For compressing the 6 coefficients of the quadratic form inmethod (c), we use the suggestion of [8]. We take the squareroot of the absolute value, give it the sign of the original
coefficient, and then convert it to a 16 bit half float [17]. For
method (a), we specify each weighted representative normal
as a scalar weight multiplying a unit normal vector, and use
the Best Fit Normal method above for the unit normals. We
use the half float of the square root for the scalar weights
and the variances, and also for the spherical harmonics coefficients for method (b).

8 Results
We tested the methods described above on several treemodels, including a broad-leafed cottonwood tree and aneedle-leafed lodgepole pine. For the cottonwood “groundtruth” images, we ray traced the full polygonal model at 8 by
8 supersampling resolution. We also generated images by ray
tracing in the SVO, using the hard surface representation forthe leaves, as well as for the trunks and branches. Then, we
generated octrees and images using the hard surface representation for the trunk and branches, and ray tracing throughthe microflake volume density for the leaves, using the three
approximations methods discussed above for the distributionof microflake normals. The octrees had resolution 10243 . We
tested the spherical harmonic method for both order L = 2
with 9 basis functions and for order L = 4, with 25 basis
functions, and tested all methods at both 8 by 8 supersampling (with a correspondingly reduced ray differential), andwithout supersampling, where the summarization qualities of
the SVO are more important. Figure 4 shows the cottonwood
tree in 640 by 640 pixel resolution images with 8 by 8 supersampling, which used the lowest level leaf cells of the octree
when ray tracing. Note that the shading is washed out for the
L = 2 spherical harmonics, because the distribution of nor-123

N. Max et al.

Fig. 4 Cottonwood tree, at 640 by 640 resolution, and 8 by 8 supersampling. From left to right ground truth image from original polygons, hardsurface image from SVO, k means with k = 3, SH with L = 2, SH with L = 4, and ellipsoidal microflakes

Fig. 5 Cottonwood tree, at 200 by 200 resolution, and no supersampling except for the ground truth image which is 8 by 8 supersampled.
Top row from left to right ground truth image from original polygons,
hard surface image from SVO, k means with k = 3, SH with L = 2,

SH with L = 4, and ellipsoidal microflakes. Bottom row differences of
the rightmost five images in the top row and the leftmost ground truth
image, magnified by a factor of five, with neutral gray representing zeromals is too band limited. Figure 5 shows non-supersampled
images at 200 by 200 resolution, on the same octrees as Fig. 4,
where the much larger ray differentials cause the higher levelinternal nodes of the octree to be used. Both figures arewithout any shadows. The second row in Fig. 5 shows thedifferences between the SVO-based images in the first row
and the ground truth ray traced and supersampled polygonbased image of the same resolution, multiplied by five.
For the lodgepole pine, we expanded the vertices ofthe polylines defining the curved pine needles into regularhexagons, so that the needle segments were represented byapproximate hexagonal prisms. This means that the distribution of the normals to the six prism faces for a single
needle segment can be exactly represented by the k means
approximation with k = 3, when vectors of exactly opposite
directions are considered equivalent. The octrees had resolution 5123 . Figure 6 shows 200 by 400 resolution images. At
the left is a ray traced image of the polygonal model, at 40
by 40 supersampling, and the subsequent images are fromthe SVO, using the different microflake normal distribution
approximations, with 2 by 2 supersampling. The second row
again shows differences.

Tables 1, 2, and 3 show the performance statistics for generating the images in Figs. 4, 5, and 6, respectively. The SVO
sizes are in bytes, and the RMS errors are for the 3 floating
point color components in the range from 0 to 1, including thesky and ground regions with no error. Table 2 used the sameoctrees as Table 1. The timings were done on an HP Z800
workstation with 12 Intel Xeon X5660 2.88 GHz cores. The
octree construction was serial, but the ray tracing was done inparallel using 10 cores, with the other 2 for task management
and system calls. For the cottonwood tree, the rendering times
were per frame and do not include initialization and loading
the precomputed octree, while for the lodgepole pine, these
initialization and loading times are included.
Figure 7 shows the cottonwood tree with shadows at 400
by 400 resolution, with 4 by 4 supersampling, rendered froma 5123 SVO using the k means method (a) with k = 3. The
left-hand image used ray traced shadows and took 346.2 s
to render. For the right-hand image, we first propagated theillumination into the leaf cells of the SVO by tracing a grid of
2068 by 1462 rays from the light source through the octree,
depositing in each intersected cell the ray segment lengthand the length-weighted attenuated flux it carried. After all

123Approximations for the distribution of microflake normals

Fig. 6 Lodgepole pine tree, at 200 by 400 resolution. Top row fromleft to right ground truth image from original polygons; hard surface
image from SVO; k means microflakes with k = 3; spherical harmonics
microflakes with L = 4 with all 25 terms, and double-sided shading
as in Sect. 6; spherical harmonics microflakes with L = 4 with only

the 15 odd order terms, shaded by sampling from the microflake normal distribution; and ellipsoidal microflakes. Bottom row differences of
the rightmost five images in the top row and the leftmost ground truth
image, magnified by a factor of five

Table 1 Statistics for cottonwood tree images at 640 by 640 resolution,with 8 by 8 supersampling

Table 3 Statistics for lodgepole pine at 200 by 400 resolution

MethodSVO sizeRMS

SVO tmRenderHard surface

1,062,142,164

0.0449

7:24

17.87

MethodSVO sizeRMS

SVO tmRenderHard surface

123,734,707

0.067

13:54

0.48

k means

190,884,125

0.056

16:54

1.16
1.58

k means

2,468,438,261

0.0450

3:44

20.05

SH: 25 terms

201,329,829

0.050

15:49

SH: L = 2

2,012,479,768

0.0454

8:13

45.19

SH: 15 terms

180,358,269

0.062

16:04

1.42

SH: L = 4

2,906,925,160

0.0442

8:34

146.64

Ellipsoidal

178,261,368

0.057

34:29

1.97

Ellipsoidal

2,347,893,040

0.0451

9:57

30.07

The last two columns show SVO construction time, in minutes:seconds,
and rendering time, in seconds

The last two columns show SVO construction time, in minutes:seconds,
and rendering time, in seconds

Table 2 Statistics for cottonwood tree at 200 by 200 resolution
Method64 RMS

64 Render

1 RMS

1 RenderHard surface

0.04491.805

0.08570.097

k means0.04492.036

0.08580.096

SH: L = 2

0.04644.531

0.06400.133

SH: L = 4

0.046414.336

0.0633

0.313

Ellipsoidal

0.04643.150

0.06700.105

RMS errors with the ray traced polygonal model and rendering times
in seconds are shown for 64 samples per pixel, and 1 sample per pixel

Fig. 7 Cottonwood tree with shadows. Left By ray tracing from eachshaded point, and right precomputed shadows123

N. Max et al.

Fig. 8 Path traced beech tree. Left reflection plus transmission, righttransmission only

rays were propagated, we divided the accumulated flux ineach leaf cell by the accumulated length and averaged thisillumination up the octree hierarchy. This illumination precomputation took 13.6 s, and then the rendering took 60.1 s.
The flux propagation was done in parallel, by dividing thegrid of rays up into square subregions indexed by column
i and row j, and doing four passes, with pass k including
only the squares with (i mod 2) + 2( j mod 2) = k. Squares
were assigned to different threads, and since the squares ineach pass were well separated, there was no hazard from two
threads attempting to simultaneously add to the flux or length
in the same octree cell.
Figure 1 shows a forest of 218 trees of various species, rendered by path tracing. The next direction for each path waschosen by sampling a surface normal from the spherical harmonics microflake distribution using only the 15 even terms
up to L = 4, reversing the normal if it was not facing theviewer, and then sampling from the diffuse BRDF with that
normal. Please see the supplemental material videos showing
path traced 360 frame cycles with the camera rotating aroundthe center of this forest, using the various methods discussed
in this paper, and also a camera dolly showing smooth transition between octree levels.
Figure 8 shows a path traced 500 by 500 image of a beech
tree with double-sided microflakes from an SH representation with L = 4. The left-hand image includes both the
diffuse reflection and the diffuse transmission, while theright-hand image includes only the diffuse transmission, toshow its contribution. The main visible effect of the pathtracing is the illumination from the sky, particularly evident
in the shadows. The supplementary material also includes apath traced video with the camera rotating around this tree.
Figure 9 shows images of a head model with 60,000 hairs,
each represented by a polygonal quad mesh on a curved cylinder following a Bezier curve, with 25 segments along the
curve and 12 around the circular cross section. They are rendered at 250 by 250 resolution, with 6 by 6 supersamples
except for the left-hand ground truth image, in which theoriginal polygonal representation was ray traced at 20 by123

20 supersamples. The polygon mesh file was 1,355,638,498
bytes. The images in the top row were rendered from uncompressed octrees. Table 4 gives the octree sizes for theseimages, and the RMS errors to ground truth and from compression. The construction and rendering times shown in
Table 4 were on a 4 core machine. The construction times
include precomputing the illumination, and compression, andthe compressed rendering times in the last column include
decompression. The octree resolution was only 2563 sothat each octree cell intersected multiple hairs, stressing theboth the summarization ability of the different representation methods, and also their compression accuracy. Moon etal. [22] also used spherical harmonics for rendering hair, but
used them to store the incoming multiple scattered radianceat each cell, as computed from a Monte Carlo simulation,rather than to store the distribution of microflake normals.
Figure 10 shows a sphere, rendered as a hard surface byadjusting the extinction and inscattering computations, at thefirst nonempty cell the viewing ray encounters, to represent acompletely opaque surface with the probability distributionof normals given by the approximated h(x, ωn ). The hardsurface and k means method agree with the “ground truth”
image ray traced from the polygonal model. The sphericalharmonics method shows a main highlight that is too broad,and also a spurious highlight ring near the profile, caused bya positive ringing lobe in the SH fit, and also by the effect of“double-sided lighting” from the double-sided microflakes.
This spurious highlight is mostly removed by the precomputed shadows. The ellipsoidal microflake method is more
effective, but must still broaden the surface normal distribution somewhat, because an almost flat ellipsoid generates
floating point exceptions or numerical problems in the calculations and cannot be used. Thus, the shadow terminator
is softer than in the ground truth image.

9 Discussion
Bruneton and Neyret [1] survey methods of prefiltering
surface properties, specifically discussing distributions of
microfacet normals in their Sect. 3, and a more general classification of prefiltering methods into three classes in theirSect. 7. In terms of that classification, our method (a) using
the variance is a spanning set method, our method (b) is abasis function method, and our method (c) appears to be a
moment method, like those of Olano et al. [24] and [25].
However, the ellipsoids of [8] represent projected microflake
area, instead of the unprojected area-weighted normals in[25]. As mentioned above, this means that the coefficients
do not add linearly when the distributions of child cells are
combined into their parent cell. On the other hand, the 3D
covariance method of [25] gives the same information as ourspherical harmonic method with L = 2. In general, since theApproximations for the distribution of microflake normals

Fig. 9 Head with hair: top row from left to right ground truth image
ray traced from the polygonal model, hard surfaces, k means with k =3, SH with L = 4, and ellipsoidal microflakes. Middle row differences
between the image in the top row and the ground truth image on the left,scaled by a factor of 5, with zero as the neutral gray in the background.
Bottom row differences between the image rendered from a compressed
octree and the corresponding image in the top row, scaled by 5

Table 4 Statistics for compression
MethodSVO size

Compr. size

Polygons RMS

Compr. RMS

Constr. tm (s)Render tm (s)Hard surface
k means

Compr. rend (s)165,151,806

64,001,340

0.12930.0196

457

3.02

4.76

410,912,420

108,682,498

0.11260.0061

596

9.63

86.3

SH: 25 terms

432,539,169

166,440,279

0.09270.0467

568

21.1

147.9

Ellipsoidal

388,167,513

84,975,047

0.09760.0185

1088

21.1

36.1

Construction time includes SVO construction, precomputed illumination, and compression

Fig. 10 From left to right a sphere rendered from polygons, hard surfaces, k means with k = 3, spherical harmonics with L = 4 and only
even-ordered terms (showing the effects of double-sided microflakes),

spherical harmonics with L = 4 and only even-ordered terms with precomputed shadows to suppress the spurious highlights, and ellipsoidal
microflakes

123

N. Max et al.

spherical harmonics of order up to L are a basis for the polynomials in x, y, and z on the unit sphere of degree up to L,
our method (b) is equivalent to the moment method of [1].
For the cottonwood tree at one sample per pixel and 200 by
200 resolution, simulating a distant view, the L = 4 sphericalharmonics method has lowest error, although it requires mostprocessing time. However, this test puts the ellipsoid method
at a disadvantage, since a single sampled normal per pixel,
required because there is no analytic integral for the diffuse
reflectance from this distribution, introduces a lot of noise.
For the lodgepole pine tree, we used 2 by 2 supersampling,
and the noise in the ellipsoid method is less. For this model,
the SH method with all 25 terms up to L = 4 had the least
RMS error, but its SVO used the most storage.
For the top left image in Fig. 6, the optimized Embree
ray tracer, running in parallel on the same 12 core machine,
took 12.42 s to create its acceleration structure and ray trace
the 200 by 400 image at 40 by 40 supersamples per pixel.
Thus, the times quoted in Table 3 are not a significant speed
improvement. However, the acceleration structure for ray
tracing the scene in Fig. 1 with Embree could not fit in memory. The model for Fig. 1 was in fact created from rotated
and translated copies of the same four trees, but if every tree
was different, its data size would be unmanageable, while
even the uncompressed SVO sizes are practical for rendering, varying from 1,670,384,307 bytes for the hard surfaces
to 2,717,912,229 bytes for the 25 SH terms up to L = 4, at
10243 voxel resolution.
The bottom row of difference images in Fig. 9 and theRMS errors from Table 4 show that the compression does
not introduce much error. The top row of images shows
visible artifacts from the SVO methods, particularly in thehighlights. When many hairs are in an octree cell, the single
average normal used in the hard surface method cannot adequately represent them. The SH normals are better for diffuse
shading, but inadequate to represent the highlights, becauseof the band-limiting effect of a practical maximum order L.
The k means method gives discontinuous shading unless thestarting angular phase for the 12 facets around the cylinder
is chosen randomly per hair, as was done for Fig. 9, so that
there are no preferred representative directions, but then thehighlights are noisy. As expected, since the distribution of
normals of a long thin ellipsoid is a good approximation tothat of coherently oriented hairs, the ellipsoidal method performs best in this case.
The band-limiting effect of the SH approximation makes
the highlight too large in Fig. 10 as well as in Fig. 9. The SH
method performs better for wider, more general distributions,where concentrated highlights do not occur. The k means
method would require more representatives to handle general distributions. When the variance is taken into account, it
is similar to a gaussian mixture model. The iterative Expectation Maximization method of fitting a gaussian mixture

123

Table 5 The suitability of the three approximation methods for different sorts of distributions of microflake normals: + means suitable, −means unsuitable, a blank means suitable, but not the best, and e means
could work with enough terms
Distribution typek means

SH

Smooth surface+−

One narrow peak+−

One wide peak+

Ellipsoids++

Hair cylinderse+

Many peaks

e

e

−

Generalee−model might give better approximations for the same data
size, but it would take much longer than the k means method.
Table 5 summarizes the suitability of the three methods fordifferent sorts of applications.
Publications [9] and [10] consider microfacet height fieldsurface models where the color is correlated with the heightand can account for the masking effects related to this correlation during rendering. None of the volume microflake models
we consider can account for this sort of volume masking,because we assume that the summarized microflake properties are homogeneous within each octree cell.
Acknowledgements We thank Mark Meyer, Tony DeRose, Eric Heitz,
Wojciech Jarosz, Derek Nowrouzezahrai and Ted Kim for helpful discussions, and the SIGGRAPH, Pacific Graphics, and Visual Computer
reviewers for useful suggestions. Nelson Max thanks the University ofCalifornia, Davis for sabbatical salary.

Appendix
In this appendix, we give some formulas for the spherical harmonics calculations discussed above. Good references for themathematics of spherical harmonics as applied to computer
graphics are [28] and [5]. Sloan [28] gives simple formulas
for the Ylm (ω) in the form of polynomials of degree l in thex, y, and z coordinates of the unit vector ω, which we usedin our implementation.
To compute the extinction integral in Eq. (3) using anh(ω) of the form of Eq. (1), we expand the nonnegatively
clamped cosine of θ , A(θ ), into spherical harmonics. Since
this function is independent of φ, all the coefficients with m =0 are zero so we get just the “zonal” harmonics expansionA(θ, φ) ≈LAl0 Yl0 (θ, φ)

(12)

l=0

where, by Eq. (2),Al0 = 2π π/20Yl0 (θ, φ) cosθ sinθ dθ dφ.

(13)

Approximations for the distribution of microflake normals

Ramamoorthi and Hanrahan [26] (please see the correction of their Eq. (19) at the bottom of the web page http://cseweb.ucsd.edu/ravir/papers/invlamb/) give the formula forthe coefficients Al0 , which are zero for odd l > 1, anddecrease rapidly for increasing even l. We have used √only the√√first few values A00 = π/2, A10 = π/3, A20 = 5π /8,√
A30 = 0, and A40 = − π /16.
This zonal harmonic expansion is only useful for evaluating Eq. (3) if ωo is at the north pole of the unit sphere. For
other directions of ωo , one must rotate this zonal harmonic
expansion so that its axis of symmetry is along the direction ωo . This is simpler for zonal harmonics than for general
spherical harmonics. According to [26], the result of rotating
the expansion in equation (12) to make the axis of symmetrylie in the direction ωo = (θo , φo ) isR θo φo A(θ, φ)
≈

L 
l



Ylm (θo , φo )l=0 m=−l4π
Al0 Ylm (θ, φ).
2l + 1(14)

Thus, if Eq. (1) is the spherical harmonic expansion of
the microflake normal density h(ω), by equation (3) and theorthonormality of the SH basis functions,σt (x, θo , φo ) ≈



L 
l


Ylm (θo , φo )l=0 m=−l4πAl0 h lm . (15)2l + 1


L
L
l
l


4π
Iα  
h lm
Yl m (θo , φo )π2l + 1=l=0 m=−l

Al 0

L


l


l =0 m =−ll =0 m =−l

Yl m (θi , φi )4π
Al 0 Tlml m l m
2l + 1(16)

where Tlml m l m is the “triple product integral”
Tlml m l m 2π  π=
Ylm (θ, φ)Yl m (θ, φ)Yl m (θ, φ) sinθ dθ dφ.
0

0

(17)
For order L = 4, the potential number of triple productintegrals is 253 = 15,625, not accounting for symmetries in
the indices, but only 1,158 of them are nonzero, and theseare precomputed and stored in a sparse table. In fact, sinceA30 = 0, only 605 terms are actually included in the sum
in Eq. (16), and if the microflakes are double-sided, and weapproximate |cos(θ )| by A(θ, φ) + A(−θ, φ), all terms withl = 1 or l = 1 also disappear, and only 585 terms remain.
Note that if ωo = (θo , φo ) is fixed, grouping the terms
and factors from h(ωn )ωn , ωo  in the last form in Eq. (16)
into the large parentheses below gives a spherical harmonic
expansion in the direction variable ωi = (θi , φi ).
S(ωo )


l
l
l
L
L
L




Iα  
h lm
Yl m (θo , φo )π
l=0 m=−ll =0 m =−ll =0 m =−l



4π
4π
Al 0
Al 0 Tlml m l m Yl m (θi , φi )×
2l + 12l + 1=

The inscattering integral in Eq. (4) contains the clamped
cosine ωn , ωi  as well as ωn , ωo , so if ωn , ωi  is also
expanded in spherical harmonics as in Eq. (14), we get, forthe perfectly diffuse reflection of a small light source fromdirection ωi of radiance times solid angle equal to I ,S(ωo ) 
α=h(ωn ) ωn , ωi ωn , ωo L(ωi )dωn dωiπI α 2π π=h(θ, φ)R θi φi A(θ, φ)R θo φo A(θ, φ)π 00sinθ dθ dφLl
I α 2π π  ≈h lm Ylm (θ, φ)
π 00
l=0 m=−l

L


l




Yl m (θo , φo )l =0 m =−l
l

m =−l

Yl m (θi , φi )
4π
Al 0 Yl m (θ, φ)
2l + 1L

l =0

4π
Al 0 Yl m (θ, φ) sinθ dθ dφ
2l + 1(18)
Thus, it can be importance-sampled for path tracing using
the techniques of [14], or by rejection sampling. Note also
from this form that once the basis function coefficients inthe large parentheses have been computed, the inscatteringof illumination from multiple light sources can be computed with only (L + 1)2 multiplications and adds for each.
For a light source “at infinity” with constant ωi , a different grouping of the terms involving instead the factors
ωn , ωi ωn , ωo  can be computed once per viewing ray thathits microflakes, allowing only (L + 1)2 multiplications peroctree cell (9.2 s for the image in Fig. 9 instead of 21.1). For
finite distance light sources, using more storage, the inscattering and the extinction coefficient can be computed once
per voxel as first needed and saved for subsequent viewing rays crossing the same voxel (8.4 s for the image in
Fig. 9).

123

N. Max et al.

References
1. Bruneton, E., Neyret, F.: A survey of non-linear pre-filtering methods for efficient and accurate surface shading. IEEE Trans. Vis.
Comput. Gr. 18(3), 242–260 (2008)
2. Cook, R., Carpenter, L., Catmull, E.: The Reyes image renderingarchitecture. In: SIGGRAPH ’87: proceedings of the 14th annual
conference on computer graphics and interactive techniques, pp.
95–102. ACM (1987). doi:10.1145/37402.37414
3. Cook, R., Torrance, K.: A reflectance model for computer graphics.
ACM Trans. Gr. 1(1), 7–12 (1982). doi:10.1145/357290.357293
4. Crow, F.: A more flexible image generation environment. In: ACM
SIGGRAPH ’82: Proceedings of the 9th annual conference on computer graphics and interactive techniques, vol. 16, pp. 9–18 (1982).
doi:10.1145/800064.801253
5. Green, R.: Spherical harmonic lighting: the gritty details. In:
Archives of the Game Developers Conference, vol. 56 (2003)
6. Heitz, E.: Understanding the masking-shadowing function in
microfacet-based brdfs. J. Comput. Gr. Tech. 3(2), 32–91 (2014)
7. Heitz, E., Dupuy, J., Crassin, C., Dachsbacher, C.: The SGGX
microflake distribution: supplementary material. ACM Digital
Library (2015)
8. Heitz, E., Dupuy, J., Crassin, C., Dachsbacher, C.: The SGGX
microflake distribution. ACM Trans. Gr. 34(4), 48 (2015)
9. Heitz, E., Neyret, F.: Representing appearance and prefiltering subpixel data in sparse voxel octrees. In: Dachsbacher,C., Munkberg, J., Pantaleoni, J. (eds.) Eurographics/ ACM
SIGGRAPH symposium on high performance graphics. The
Eurographics Association (2012). doi:10.2312/EGGH/HPG12/
125-134
10. Heitz, E., Nowrouzezahrai, D., Poulin, P., Neyret, F.: Filtering color
mapped textures and surfaces. In: Proceedings of I3D, pp. 129–136.
ACM (2013)
11. Hoppe, H.: Progressive meshes. In: ACM SIGGRAPH ’96: proceedings of the 24th annual conference on computer graphics andinteractive techniques, pp. 99–108. ACM Press (1996). doi:10.
1145/237170.237216
12. Hoppe, H.: View-dependent refinement of progressive meshes. In:
ACM SIGGRAPH ’97: proceedings of the 24th annual conferenceon computer graphics and interactive techniques, pp. 189–198.
ACM Press/Addison-Wesley Publishing Co., New York (1997).
doi:10.1145/258734.258843
13. Jakob, W., Arbee, A., Moon, J.T., Bala, K., Marschner, S.: A radiative transfer framework for rendering materials with anisotropic
structure. ACM Transactions on Graphics (TOG)—proceedingsof ACM SIGGRAPH 2010 29(4) (2010). doi:10.1145/1778765.
1778790
14. Jarosz, W., Carr, N., Jensen, H.W.: Importance sampling sphericalharmonics. Comput. Gr. Forum 28(2), 577–586 (2009). doi:10.
1111/j.1467-8659.2009.01398.x
15. Jensen, H.W.: Realistic Image Synthesis Using Photon Mapping.
A K Peters/CRC Press, Boca Raton (2012)
16. Joy, K.I., Grant, C.W., Max, N.L., Hatfield, L.: Tutorial: computer
graphics: image synthesis. ACM Press, p 2 (1988)
17. Kainz, F., Bogart, R.: https://github.com/openexr/openexr/tree/
master/ilmbase/half (2002)
18. Kanungo, T., Mount, D.M., Netanyahu, N., Piatko, C., Silverman,
R., Wu, A.Y.: A local search approximation algorithm for k-meansclustering. Comput. Geom. Theory Appl. 28, 89–112 (2004)
19. Kaplanyan: Cryengine 3: reaching the speed of light. In: SIGGRAPH 2010: advances in Real-time Rendering course (2010)
20. Max, N.: Optical models for direct volume rendering. IEEE Trans.
Vis. Comput. Gr. 1, 99–108 (1995)
21. Max, N., Mobley, C., Wu, E.H.: Plane-parallel radiance transport
for global illumination in vegetation. In: Rendering techniques ’97:

123

proceedings of eurographics workshop on rendering, pp. 239–250.
Springer, Berlin (1997)
22. Moon, J., Walter, B., Marschner, S.: Efficient multiple scattering inhair using spherical harmonics. ACM Trans. Gr. 27(3), 1–7 (2008)
23. Neyret, F.: Modeling animating and rendering complex scenes
using volumetric textures. IEEE Trans. Vis. Comput. Gr. 4(1), 55–
70 (1998). doi:10.1109/2945.675652
24. Olano, M., Baker, D.: Lean mapping. In: Proceedings of I3D’10,
pp. 181–188. ACM (2010)
25. Olano, M., North, M.: Normal distribution mapping. Tech. Rep.
TR97-041, Department of Computer Science, University of North
Carolina, Chapel Hill (1997)
26. Ramamoorthi, R., Hanrahan, P.: On the relationship between radiance and irradiance: determining the illumination from images of aconvex Lambertian object. J. Opt. Soc. Am. A 18(10), 2448–2459(2001)
27. Rense, W.: Polarization studies of light diffusely reflected fromground and etched glass surfaces. J. Opt. Soc. Am. 40(1), 55–59(1950). doi:10.1364/JOSA.40.000055
28. Sloan, P.P.: Stupid spherical harmonics (sh) tricks. Microsoft
Corporation (2008). Expanded material from presentation at
the 2008 Game Developer Conference. http://www.ppsloan.org/
publications/StupidSH36.pdf
29. Toksvig, M.: Mipmapping normal maps. J. Gr. Tools 10(3), 65–71(2005)
30. Torrance, K., Sparrow, E.: Theory for off-specular reflection fromroughened surfaces. J. Opt. Soc. Am. 57(9), 1105–1114 (1967)
31. Walter, B., Marschner, S.R., Li, H., Torrance, K.E.: Microfacet models for refraction through rough surfaces. In: Kautz,J., Pattanaik, S. (eds.) Rendering techniques, pp. 105–205. The
Eurographics Association (2007). doi:10.2312/EGWR/EGSR07/
195-206
32. Zhao, S.: Modeling and rendering fabrics at micron resolution.
Ph.D. thesis, Cornell University (2014)
33. Zhao, S., Jacob, W., Marchner, S., Bala, K.: Building volumetricappearance of fabrics using micro CT imaging. ACM Trans. Gr.
30(4), 44:1–44:10 (2011)

Nelson Max received his Ph.D.
in Mathematics from Harvard
University in 1967 and is currently a Distinguished Professor of Computer Science at the
University of California, Davis.
His research interests are in
the areas of scientific visualization, computer animation, realistic computer graphics rendering, and multi-view stereo reconstruction. In visualization, heworks on molecular graphics,
and volume and flow visualization, particularly on irregular
finite element meshes. He has rendered realistic lighting effects inclouds, trees, and water waves and has produced numerous computer
animations, shown at the annual ACM SIGGRAPH conferences, andin OMNIMAX stereo at the Fujitu Pavilions at Expo ’85 in Tsukuba
Japan, and Expo ’90 in Osaka Japan. He is a fellow of the Association
for Computing Machinery, and the winner of its 2007 Steven A. Coonsaward.

Approximations for the distribution of microflake normals
Tom Duff was educated at
the universities of Waterloo andToronto and has worked at the
NYIT Computer Graphics Lab,
the Lucasfilm computer division,
Bell Labs and, currently, Pixar.
His main research interest is incomputer graphics and animation, but he has also publishedin operating systems, polymer
physics, wireless networking,mathematics and computer security. Past projects include writing(with Bill Reeves) the first version of Pixar’s animation system,
special effects in Star Trek II, Return of the Jedi and A Bugs Life, afull-scale radio-controlled Walle and a light-field movie camera. He hasa dozen patents and two Academy Awards for contributions to motion
picture science and technology. When not at work, he is an opera singer,
musical tinkerer and grandparent.

Yajie Yan is a Ph.D. student
at Washington University in St
Louis. His research at school
focuses on geometry processing, specifically robust computation and applications of medial
axis. He obtained his master
degree from Beihang University in 2012, working on feature extraction from 3D point
clouds. Besides geometry, he isalso interested in realistic rendering and 3D fusion applied in
VR/AR contexts.

Ben Mildenhall is a Ph.D. student at the University of California, Berkeley, where his work
is supported by a Fannie andJohn Hertz Foundation fellowship. His current research interests are 3D reconstruction for virtual reality, structured light capture, and light field reconstruction. He received his undergraduate degree from Stanford University in 2015, where his work
focused on applications of probabilistic programming to content
generation.

123