Importance Sampling for Physically-Based Hair Fiber Models
Eugene d’Eon1
Steve Marschner2Johannes Hanika1
1 Weta Digital
2 Cornell University

Figure 1: Our new importance sampling strategy allows easy inclusion of Marschner and related hair reflectance functions in physicallybased Monte Carlo renderers. Here we show hair volumes illuminated by environment maps and area lights with unbiased global illumination(computed using a forward path-tracer with multiple importance sampling). Our sampling strategy requires no precomputation, so it is easy
to vary the absorption along the fiber (second image), and to add noise to the index of refraction, roughness, and scale tilt to create subtle
heterogeneity along each fiber. Each image is 1024 samples/pixel.

Abstract
We present a new strategy for importance sampling hair reflectancemodels. To combine hair reflectance models with increasingly popular physically-based rendering algorithms, an efficient sampling
scheme is required to select scattered rays that lead to lower variance and noise. Our new strategy, which is tied closely to thederivation of physically-based fiber functions, works well for bothsmooth and rough fibers based on the Marschner et al. model andalso for Lambertian fibers. It should be directly usable with future
hair reflectance models that allow for more general cross-sections
and more complex surface properties, provided the lobes are derived in a similar, separable fashion. Our strategy includes lobe
selection and can efficiently sample complex lobe shapes like the
Marschner TRT function. The scheme is easy to implement andrequires no precomputation, allowing fully heterogeneous variationof all fiber parameters.
CR Categories: I.3.7 [Computer Graphics]: Three-Dimensional
Graphics and Realism—Color, shading, shadowing, and texture;
I.3.7 [Computer Graphics]: Three-Dimensional Graphics andRealism—Raytracing;
Keywords: Hair scattering, BCSDF, importance sampling

1

Introduction

Believable computer generated characters require believable hair.
For rendered hair to appear realistic the Bidirectional Curve Scattering Distribution Function (BCSDF) used to simulate the lightαa Longitudinal Analysis ΘFigure 2: Refracted pathways through a hair fiber. Single
Reflection (R), double transmission (TT), transmission-reflectiontransmission (TRT) and so on.
interaction with the individual hair fibers must be accurate—all of
the visibly relevant behaviours seen in plausible illuminations oftypical hair volumes must be well approximated.
Previously in graphics, the fiber model of Kajiya and Kay [1989]
remained predominant until Marschner et al. [2003] introduced thefactored-lobe analytic BCSDF that remains the basis of most highaccuracy, parameteric hair rendering today. Factored-lobe BCSDFs
decompose the reflectance into separate modes of propagation—
direct reflection (R), double transmission, (TT) and paths with one
or more internal reflections (TRT, TRRT, . . . ) (Figure 2). The total reflectance function S is the sum of all such lobes S p , indexedby p, the number of internal path segments traversed by light rays
contributing to that mode (Figure 2),∞S(θi , θo , φ ) = ∑ S p (θi , θo , φ ).

(1)

p=0

Many Monte Carlo rendering algorithms, such as path tracing, require efficient schemes for randomly selecting samples (directions
and weights) with a statistical distribution that is close to the reflectance function of the material being simulated. We derive apractical, analytic, easy-to-implement importance-sampling strategy for several recent physically-based factored-lobe BCSDFs, allowing their efficient combination with these rendering algorithms.

2

Related work

Our new sampling strategy works for fiber models including thoseof Zinke and Weber [2007], who introduced the formalism of BCSDFs. They also derived a near-field model, important for close
rendering of hair fibers. We use the energy-conserving longitudinalterm of d’Eon et al. [2011] and Zinke and Weber’s [2007] methodof treating azimuthal roughness—these remove the requirement forroot-solvers and careful treatment of caustics required in the original Marschner model.
The Marschner, Zinke and d’Eon models are the most physicallyaccurate parametric models, but no analytic sampling scheme has
been previously presented for these full models. Hery and Ramamoorthi [2012] presented importance sampling for the R lobe of
Marschner et al. [2003] and Ou et al. [2012] presented importancesampling for the full model of Sadeghi et al. [2010]. The Sadeghimodel allows more intuitive user control of the hair appearance,
but isn’t as physically accurate as other models. Our sampling of
the NR function closely corresponds to both of these related works,
but our approach also handles the complex T T and T RT lobes thatclosely match the transport of light within dielectric cylinders. In
addition, we handle lobe selection automatically in a fashion thatworks well regardless of the angle of incidence or absorption level
and we present the first longitudinal sampling scheme that accounts
for all angular variables.
Our sampling scheme can be used for direct illumination sampling(where indirect illumination may be approximated by Dual Scattering [Zinke et al. 2008]) and also for indirect illumination, where
sampling at each vertex of a random path is used to create photon
maps within hair volumes [Moon et al. 2008] or within a varietyof unbiased path-tracing algorithms [Veach 1997]. For an extensive
review of hair shading, animation and modeling techniques we refer
the reader to the survey by Ward et al. [2007].

3

Hair BCSDF Importance Sampling

We begin with some notation. Following the convention of
Marschner et al. [2003], directions are measured in a spherical coordinate system centered on the fiber axis, with θi and θo measuring
inclinations to the normal plane of the fiber and φ = φo − φi measuring the azimuthal difference between the incident and reflected
directions. Each lobe is factored into a longitudinal scattering function Mp and an azimuthal scattering function Np :S p (θi , θo , φ ) = Mp (θi , θo )Np (θi , θo , φ ).

(2)

Note that our notation differs from previous work by including all
inclination-dependent factors in Mp , such as the cos2 θd term of
the original Marschner model. Notation is simplified by referringto the longitudinal difference angle θd = (θo − θi )/2. The relative index of refraction η of the hair to the surrounding medium istypically fixed at 1.55. Given θd , η , and the offset h ∈ [−1, 1] ofan incoming ray from the fiber axis (Figure 2), a Bravais analysis
gives the azimuthal distributions
compactly using γi = arcsin(h),√γt = arcsin( ηh0 ), η 0 =η 2 −sin(θd )cos(θd )to predict the relative change inazimuthΦ(p, h) = 2 p γt − 2 γi + pπ

(3)for each mode p being considered. For rough fibers we use theazimuthal functions Np of d’Eon et al. [2011] evaluated using a 70point Gaussian quadratureZ
1 1
Np (φ ) =dh A(p, h)D p (v p,N , φ − Φ(p, h)).
(4)2 −1Here, the attenutation terms A(p, h) and the wrapped Gaussian D pof variance v p,N are those given by d’Eon et al. [2011]. The ef-

fects due to tilted scales are simulated approximately by perturbing
the specular cone angles slightly, θcone,R = −θi + 2α , θcone,T T =−θi − α , θcone,T RT = −θi − 4α , where α is the scale tilt angle.
We now describe our importance sampling strategy, using the hairreflectance model of d’Eon et al. [2011] as an example. Alternative longitudinal functions (and sampling for those functions) canbe substituted and the remaining elements of the approach remainunchanged.
3.1

Strategy Outline

A good importance sampling method for a BCSDF starts with afixed incoming direction (θi and φi ) and produces random outgoing directions (θo and φo ) with a distribution that is close toS(θi , θo , φ ) cos2 θo (one cosine for the spherical coordinate system
and another for the geometric projection to the fiber). We begin
with lobe selection, by using uniform random sampling of a position h along the cross section of the hair fiber. Given this offset h,
we first use the attenuations through a smooth fiber at this offset h
for selecting a lobe, then use the same h to choose θo and φ . Because
the longitudinal functions Mp in Equation 2 are independent of azimuth, we sample θo first, and then sample the azimuthal direction.
For the former procedure, we derive perfect importance sampling
for the energy-conserving longitudinal term of d’Eon et al. [2011].
Then we show how to use Equation 3 to importance sample the azimuthal functions Np for smooth (v p,N = 0) hairs. To sample roughhairs, we use an additional Gaussian random variable applied to theexitant azimuth predicted by the smooth hair sampling. The result
is an analytic, easy to implement, and efficient sampling scheme
requiring no precomputation.
3.2

Importance Sampling the Longitudinal M function

We use the longitudinal scattering function M(v, θcone , θo ) of d’Eon
et al. [2011] for all lobescsch(1/v) sin θc sin θocos θc cos θo1
M(v, θc , θo ) =e v I0. (5)2vv
cos θoFor additional numerical recipes for evaluating this function,see [d’Eon 2013]. The derivation of this M function was to convolve a Dirac circle δ (θ − θcone ) on the surface of a sphere by aspherical Gaussian and then take the longitudinal distribution ofthat result to be the spreading in θ due to roughness. To sample
this with a single random number seems intractable. However, with
two uniform random variables ξ1 and ξ2 each within [0, 1) we canimportance sample M cos2 θo perfectly. We start with the importance sampling of a spherical Gaussian of variance v (being carefulof numerical issues for low variance [Jakob 2013])


1
(6)
u(ξ1 ) = v log e1/v − 2ξ1 sinhv
and then take the final inclination θo of a point chosen uniformly on
the Dirac specular-cone circle parameterized by φ = 2πξ2 and tilted
away from the north pole by the angle θ sampled from the spherical Gaussian, cos θ = u(ξ1 ). Putting this all together, sampling θogiven variance v and the specular cone angle θcone (θi , φ , α ) isθo (ξ1 , ξ2 , v, θcone ) =arcsin(u(ξ1 ) cos θ 0 +0p
1 − u(ξ1 )2 cos(2πξ2 ) sin θ 0 )(7)

where θ = 2 − θcone . This is perfect importance sampling, so thesample weight is 1 and the PDF is M cos2 θo .
π

We note that this M function is not reciprocal due to the 1/ cos θiterm. Lack of reciprocity is common in other forms of energyconserving M functions that renormalize Gaussians over the finite

domain of inclinations based on the incoming angle θi . Also,
analagous to shading normals for surfaces, the tilted scales treatment generally yields non-reciprocal scattering functions in any
case. We haven’t found any issues with ignoring these details in
our implementation, however.
3.3

Importance Sampling Azimuthal N functions

The azimuthal functions Np of the Marschner model describe alarge family of lobe shapes that model the distributions of light
exiting a dielectric fiber when illuminated from a single direction.
However, despite their complexity and variety, they can be sampled
analytically and efficiently. The key insight is to exploit the original
derivation of Marschner et al. [2003]—the fiber function is a uniform integral over the fiber cross-section h ∈ [−1, 1] (Equation 4)
and the laws relating an incident azimuth to an exitant azimuth
are compact expressions. In fact, the derivation of the original
Marschner model was precisely the calculation of the outgoing densities that arise from this sampling scheme.
We first select a random cross-section offset h = 2ξh − 1 (whereξh ∈ [0, 1) is a uniformly distributed random number). For smoothfibers, the exitant normal-plane inclination is θo = θcone p (θi ) fromwhich θd is then known. In the case of no azimuthal roughness, theset of directions Φ(p, 2ξh − 1), using Equation 3, each with weightA(p, 2ξh − 1) ≤ 1, is an efficient scheme for importance sampling
the functions Np derived by Marschner et al. [2003]. In the casep = 0, this is identical to prior sampling schemes for NR [Hery andRamamoorthi 2012; Ou et al. 2012].
In the case of azimuthal roughness the functions we wish to sample
are the convolution of Φ(p, h) by a wrapped Gaussian D p (v p,N )with variance v p,N (Equation 4). These are easily sampled by selecting relative azimuths√
φ = Φ(p, 2ξh − 1) + g v p,N(8)
where g is a Gaussian distributed random variable with standard
deviation 1 (sampled with a lookup table or Box-Muller). Again,
the sample weight is A(p, h) ≤ 1 and the PDF is Equation 4 withA(p, h) replaced by 1. We mention that sampling the Np functions
is quite cheap—much cheaper than evaluating them or computing
their PDF, due to the quadrature.
3.4

Lobe Selection

We now return to lobe selection. The key to lobe selection is todelay the choice until after randomly selecting h. Given h, theattenuations through a smooth, ideally specular fiber, Aspec (p, h),
are known—they are the product of Fresnel terms and absorption terms for the single set of paths predicted for the smoothfiber (assuming M is a delta function at the ideal specular coneangle when computing θd ). We use these attenuations for thesmooth fiber to select a lobe p to sample, even for rough fibers.
When the inclination is high and Fresnel dominates, or absorption
is high, we tend to select the R lobe, for example. A random
number ξ p ∈ [0, 1) can be used to choose a lobe p in proportion
to w p = Aspec (p, h)/ ∑ p Aspec (p, h). Selecting p proportional tow p incurs a selection weight of 1/w p . The attenuation throughthe actual path (which deviates from the specular cone path after sampling M) has a throughput of A(p, h), so the final sample
weight is w = A(p, h)/w p , which is ∑ p Aspec (p, h) times the ratioA(p, h)/Aspec (p, h). Provided the roughness is reasonably low, theFresnel and absorption terms for slightly deflected paths will notdiffer significantly from the ideal specular case, keeping the final
sample weights low. However, for extreme inclinations and roughnesses, we clamp the upper bound of the final sample weights to

2.0 (and have noted no measurable energy loss in a variety of hairrenders).
3.5

Sampling Summary and PDF Computation

To summarize, our importance sampling strategy for the full fiber
function with rough surfaces is:• An incoming direction θi is given• We uniformly choose a random offset h ∈ [−1, 1] along the fibercross-section• We compute attenuations Aspec (h, p) for each lobe p assuming
no deflection away from the specular cone angles• We select a lobe in proportion to the specular attenuations• We now know which Mp function (of variance v p ) to sample
using Equation 7, giving θo and θd .
• We sample a random Gaussian variable g and use Equation 8 tocompute the relative outgoing azimuth.
• We return a sample weight of w = A(p, h)/w p
The PDFs of this sampling scheme are exactly analogous to evaluation of the model S, but with the attenuations A replaced by theselection weights w p , which are easily evaluated for R and T T bysolving for the single root h p . For the more complicated TRT lobe,
the quadrature evaluation of NT RT with the selection weights w p inplace of A is required, and this is of comparable cost to evaluationof TRT itself. The total PDF is the sum of these individual lobe pdf
terms.

4

Lambertian Fibers

We briefly mention two new (to our knowledge) results regarding
diffuse fibers. Zinke and Weber [2007] derived the BCSDF for
an opaque, cylindrical fiber with Lambertian reflectance properties.
They expressed the exact solution as an integral and gave an approximate expansion of it. We note that the exact solution has asimple form (taking care that this holds only for φ ∈ [0, 2π ]):SLambert (θi , θo , φ ) = kd

(π − φ ) cos(φ ) + sin(φ )4π(9)where kd is the diffuse albedo. We also note that importance sampling this BCSDF is done easily with three uniform random variables by, again, selecting h ∈ [−1, 1] uniformly
and then sampling a√
Lambertian BRDF about normal (h, 0, 1 − h2 ) (in the coordinate
space where we have chosen the fiber tangent to be aligned withthe y axis and the incoming direction to be (0, sin θi , cos θi ). After
sampling an outgoing direction in this coordinate space, a rotationback to world space gives the desired direction, and a sample weightof kd is returned. We mention this approach because the CDF forthe azimuthal angle seems to have no analytic inversion.

5

Results

We have implemented our new importance sampling for the d’Eon
et al. [2011] hair model in a path tracer. Russian roulette is used together with next event estimation and multiple importance samplingat each vertex. The geometric term includes a sine between the hairfiber axis and the ray direction instead of the cosine term applied
to surfaces. Timings for a single call to our sampling function arecomparable to the time of a single BCSDF evaluation (the mostexpensive calculations are the quadrature for evaluating the TRT
function and its PDF).
Figures 1 and 3 show several results under IBL or area light illumination. We easily support a variety of heterogeneous fiber6Conclusions

We have presented new methods for importance sampling
physically-based fiber reflectance functions. The analytic formof our sampling scheme makes it easy to adopt in a variety ofMonte Carlo renderers. We demonstrated the utility of our new
methods for rendering physically-based images in a path tracer.

7

Acknowledgements

Joe Letteri, Wenzel Jakob, Guillaume Francois, and Luca Fascione
provided helpful feedback, Nicholas Wilson modeled the hair curls
and Jean-Marie Aubry kindly provided the hair dynamics simulations.

References
D ’E ON , E., F RANCOIS , G., H ILL , M., L ETTERI , J., AND AUBRY,

J. 2011. An energy-conserving hair reflectance model. In Computer Graphics Forum, vol. 30, 1181–1187.
D ’E ON , E. 2013. publons.com/p/2867/. Publons 2867.

H ERY, C., AND R AMAMOORTHI , R. 2012. Importance samplingof reflection from hair fibers. Journal of Computer Graphics
Techniques (JCGT) 1, 1 (June), 1–17.

Figure 3: Various IBL illuminations of hair.

JAKOB , W. 2013. Light transport on path-space manifolds. PhD
thesis, Cornell University.
K AJIYA , J., AND K AY, T. 1989. Rendering fur with three dimensional textures. In ACM Siggraph Computer Graphics, vol. 23,
ACM, 271–280.
M ARSCHNER , S. R., J ENSEN , H. W., C AMMARANO , M., W OR LEY, S., AND H ANRAHAN , P. 2003. Light scattering fromhuman hair fibers. In SIGGRAPH ’03: ACM SIGGRAPH 2003
Papers, ACM, New York, NY, USA, 780–791.

(a) 0

(b) 1

(c) 2

(d) 3

(e) 7

(f) 12(g) 29Figure 4: Limiting the number of indirect bounces (from left toright: 0, 1, 2, ..., 29) clearly shows how much the appearance of
hair depends on multiple scattering. The concentration of melaninin the hair is varied from top to bottom to simulate a variety of hairtypes in one dataset.
properties—a challenge for sampling strategies that use precomputation, especially for smooth fibers, which have strongly peakedreflectances and would require very high resolution datasets to store
their sampling tables. With our new importance sampling, the images converge at about 1024 samples or less and are fast enough toallow reasonable iteration time for achieving a desired look, evenwith very low absorption hairs where most of the reflectance is due
to high-order scattering (Figure 4). The requirement for such a highsample count is reduced with lower variance illumination, but sunlight illumination of light, shiny hair remains a challenge to renderboth efficiently and accurately.

M OON , J. T., WALTER , B., AND M ARSCHNER , S. 2008. Efficient
multiple scattering in hair using spherical harmonics. In ACM
Transactions on Graphics (TOG), vol. 27, ACM, 31.
O U , J., X IE , F., K RISHNAMACHARI , P., AND P ELLACINI , F.
2012. ISHair: importance sampling for hair scattering. In Computer Graphics Forum, vol. 31, 1537–1545.
S ADEGHI , I., P RITCHETT, H., J ENSEN , H., AND TAMSTORF, R.
2010. An artist friendly hair shading system. ACM Transactions
on Graphics (TOG) 29, 4, 1–10.
V EACH , E. 1997. Robust Monte Carlo Methods for Light Transport
Simulation. PhD thesis, Stanford University.
WARD , K., B ERTAILS , F., K IM , T., M ARSCHNER , S., C ANI , M.,
AND L IN , M. 2007. A survey on hair modeling: Styling, simulation, and rendering. IEEE Transactions on Visualization andComputer Graphics, 213–234.
Z INKE , A., AND W EBER , A. 2007. Light scattering from filaments. IEEE Transactions on Visualization and Computer
Graphics, 342–356.
Z INKE , A., Y UKSEL , C., W EBER , A., AND K EYSER , J. 2008.
Dual scattering approximation for fast multiple scattering in hair.
ACM Transactions on Graphics (TOG) 27, 3, 1–10.