A radiative transfer framework for rendering materials with anisotropic structure
Wenzel Jakob

Adam Arbree

Jonathan T. Moon
Cornell University

Kavita Bala

Steve Marschner(c) Detail (isotropic)(a) Isotropic scattering(b) Scattering by anisotropic micro-flakes(d) Detail (micro-flakes)

Figure 1: Renderings of a scarf represented as a high-resolution volume. Accounting for the anisotropic structure of the medium leads to asignificantly changed appearance, including realistic highlights and color variations.

Abstract

CR Categories: I.3.7 [Computing Methodologies]: Computer
Graphics—Three-Dimensional Graphics and Realism

The radiative transfer framework that underlies all current renderingof volumes is limited to scattering media whose properties are invariant to rotation. Many systems allow for “anisotropic scattering,”
in the sense that scattered intensity depends on the scattering angle,
but the standard equation assumes that the structure of the medium
is isotropic. This limitation impedes physics-based rendering of volume models of cloth, hair, skin, and other important volumetric ortranslucent materials that do have anisotropic structure. This paper
presents an end-to-end formulation of physics-based volume rendering of anisotropic scattering structures, allowing these materials to
become full participants in global illumination simulations.

Keywords: Subsurface scattering, anisotropy, light transport, diffusion theory, finite element method, dipole model, BSSRDF.

We begin with a generalized radiative transfer equation, derived
from scattering by oriented non-spherical particles. Within thisframework, we propose a new volume scattering model analogous to
the well-known family of microfacet surface reflection models; wederive an anisotropic diffusion approximation, including the weakform required for finite element solution and a way to compute thediffusion matrix from the parameters of the scattering model; andwe also derive a new anisotropic dipole BSSRDF for anisotropictranslucent materials. We demonstrate results from Monte Carlo,
finite element, and dipole simulations. All these contributions arereadily implemented in existing rendering systems for volumes andtranslucent materials, and they all reduce to the standard practice in
the isotropic case.

1

Introduction

Volume models are valuable in rendering materials with complex
surface or volumetric structure [Perlin and Hoffert 1989; Kajiya andKay 1989; Xu et al. 2001; Neyret 1998], and many very successfulprojects have made detailed visualizations and realistic renderings
using them. However, there is a fundamental gap in the technology
for rendering volumes, which impedes progress in improving therealism of volume renderings. On one hand, the tradition of volume
visualization, or “volume rendering,” achieves high detail and realistic shading effects using heuristic shading models that can describe
directional scattering from surfaces and fibers, but these models
lack a sound physical basis. On the other hand, physically based
rendering systems support volumes as “participating media” using
physically sound formulations originating from atmospheric optics;
unfortunately, these models are limited to isotropic media1 and are
inherently unable to describe the directional scattering needed toachieve realistic appearance in nontrivial solid materials.
This paper aims to bridge this gap by upgrading the underlying
1 We use the word “isotropic” to refer to a medium with properties thatare invariant to rotating the medium, in contrast to the alternative meaning of
scattering that is independent of scattered direction.

isotropic medium

anisotropic medium

Figure 2: The distinction between isotropic and anisotropic media.

framework of participating media to handle anisotropic effects, suchas specular reflections from surfaces or fibers. This allows themachinery of physics-based rendering to be brought to bear on awider range of volume models, thus providing the foundation for
a new level of realism in volume models. Currently, participating
media rendering algorithms model scattering media using a radiativetransfer model limited to isotropic media—media whose properties
are invariant to rotation. For instance, a cloud of steam filled withspherical water droplets fits within the standard model; or a cirrus
cloud filled with randomly oriented ice crystals: either medium looks
the same to light traveling in different directions. But a volume filledwith aligned ice crystals; or hair, wood, or cloth fibers; or fibrous
tissue like muscle or collagen; does not fit the usual model because
its properties depend on orientation. This is the type of medium
we call anisotropic in this paper (Figure 2). (This is a deviationfrom the common usage in graphics, where “anisotropic” refers to
scattering that depends on the scattering angle, such as the strong
forward scattering exhibited by spherical water drops.)
Anisotropic media have fundamentally different behavior fromisotropic media, but they can be handled without burdensomechanges to existing renderers. This paper outlines the changes
required to the theory and implementation of renderers for participating media and translucent materials, to make them handle
anisotropic media. Specifically, we provide:(i) an anisotropic version of the radiative transport equation (RTE)
to handle anisotropic media, derived from scattering by oriented non-spherical particles;(ii) a new model for oriented phase functions that generalizes
existing phase functions while also representing surface-likeand fiber-like scattering;(iii) an anisotropic generalization of the widely used diffusion approximation for rendering translucent materials; and(iv) methods for computing light transport solutions using Monte
Carlo to directly solve the RTE, using finite elements to solve
the diffusion approximation in a translucent material, or using
a new anisotropic version of the widely used dipole diffusion
solution.
These contributions establish the framework required for further
development of volume models as first-class objects in modern
physics-based rendering systems.

2

Prior work

Radiative transfer is used in many areas including astrophysics,
hydrology, neutron transport, remote sensing, acoustics, seismology,
and medical physics, and the standard form of the RTE comesfrom references in these fields [van de Hulst 1957; Chandrasekhar
1959; Preisendorfer 1976; Ishimaru 1978; Mishchenko et al. 2006].
The theory is usually stated for the prevailing case of sphericalor randomly oriented particles, disallowing anisotropic media (butallowing angle-dependent scattering), and this form has generally
been used in graphics.
Research in volume rendering has introduced a variety of methods
for mapping volume data into an image, including ray tracing [Levoy1988], resampling [Drebin et al. 1988; Lacroute and Levoy 1994],
and rasterization [Westover 1990] approaches, with many recent
systems relying on GPU ray casting [Crassin et al. 2009]. While
the standard ray integration equation comes from the RTE, thesesystems do not solve the full equation and instead use heuristic
shading methods derived from surface reflection models, often withthe density gradient playing the role of the surface normal.

The development of physics-based surface renderers has seen mostrendering algorithms extended to handle participating media, fromearly papers on ray tracing [Kajiya and Herzen 1984; Blinn 1982]
and radiosity [Rushmeier and Torrance 1987] to subsequent global
illumination systems [Lafortune and Willems 1996; Jensen andChristensen 1998; Pauly et al. 2000; Cerezo et al. 2005]. Most often,
participating media are used to add general fog or haze to scenes, butthey are also used to make realistic images of detailed atmosphericphenomena like smoke [Fedkiw et al. 2001], and clouds [Max 1994].
Many papers have investigated efficient real-time rendering methods
based on analytic solutions in special cases [Premoze et al. 2004;
Sun et al. 2005]. Because of the goal of physical accuracy and theneed to integrate with global illumination systems, all these methods
pay attention to correct implementation of the RTE, inheriting theimportant restriction to isotropic media.
Between these two camps lie techniques that use volume rendering
to make realistic images of materials, such as fur [Kajiya and Kay
1989], cloth [Xu et al. 2001], and other materials [Perlin and Hoffert1989], that are difficult to represent with surfaces. These models
must either stick with isotropic scattering or accept being incompatible with the physics-based-rendering toolbox—both undesirable
choices when realism is the goal. Some papers, notably the volume
modeling system of Neyret [1998], have recognized and discussed
many of the relevant effects, but have not made the connection tophysically based models of volume scattering that is required forthese models to be used in full physics-based simulations.
Many models for scattering from surfaces or in volumes have been
proposed, with microfacet models for surfaces [Cook and Torrance
1982; Ashikhmin et al. 2000; Walter et al. 2007] and the HenyeyGreenstein model for volumes [Cerezo et al. 2005] being the mostwidely used. The model proposed here has similarities to microfacet
models and to Neyret’s volumetric textures [1998], but differs in thatit provides a physically based model for the volume setting.
The diffusion approximation to radiative transfer is widely used toapproximate transport of light or other energy in highly scattering
media. It is derived from the RTE either using a first-order approximation to the directional radiance [Ishimaru 1978], as we do in thispaper, or by taking a limit of diminishing mean free path [Ryzhik
et al. 1996]. In graphics, diffusion is used to describe light transport
through translucent materials [Jensen et al. 2001]. Previous rendering work has always used isotropic diffusion, but other fields havedeveloped methods for anisotropic diffusion. Because many humantissues exhibit significant anisotropy, recent work in medical physics
has examined the effects of anisotropy [Kienle et al. 2004]. In particular, optical tomography [Arridge 1999; Gibson et al. 2005] uses
the diffusion approximation inversely to recover volume structurefrom external observations, and recent works [Heiskala et al. 2005;
Heino et al. 2003] have found anisotropic diffusion approximations
similar to ours by replacing individual components of the isotropic
RTE with more expressive counterparts. However, these approaches
do not account for the subtle interrelations between components of
the RTE and complications regarding reciprocity and normalizationin anisotropic media. Both models also lack expressiveness whenrendering directly using the RTE, being limited to isotropic scattering or unable to handle distributions of particles, respectively. Theseissues motivated the derivation of the anisotropic RTE presented inthis paper.
Dipole sources can provide simple approximate solutions when diffusion occurs in a homogeneous region with a planar boundary. Thisis the prevailing way to handle translucency in rendering [Jensen
et al. 2001], and similar formulations can be found in other fields. In
geoscience, dipole solutions have been used for transport of electriccharge in the earth, and in that context an anisotropic dipole was proposed [Li and Uren 1998], which inspired our anisotropic BSSRDF

model. In practice, many differences between the two models exist
due to the different underlying equations they solve.
In biophysics, recent work has led to both time-resolved [Dudko
and Weiss 2005] and steady-state [Johnson and Lagendijk 2009]
anisotropic multipole models for slab geometries. Both models use
absorptive boundary conditions that are inappropriate for graphics
applications, since they force the average illumination to zero on
the boundary. Another issue is that the supplied solution for theradiance on the slab boundaries only applies to the transmissionmode, whereas reflection is generally needed for rendering. Ourmodel provides an expression for this direction and, as will be seen
later, can also easily be extended to a multipole configuration.
Johnson and Lagendijk’s [2009] work presented an approach to relate micro- and macroscale anisotropy in forming their anisotropic
multipole model. Because they consider only the shape of the modeled particles, but not their scattering properties, the final diffusion
matrix derives its anisotropy solely from the particles’ directionallyvarying projected area. In comparison, we present a formulationthat accounts for anisotropy in the phase function, which leads to an
additional integral term in the resulting diffusion matrix.
The contribution of this paper is to upgrade the framework of participating media rendering so that it works for anisotropic media,
which makes it applicable to volume models of the sort that could
previously only be rendered using less-principled volume rendering methods. This involves generalizing the RTE, the diffusion
approximation, and the dipole BSSRDF model—each by employing
methods similar to methods found in the physics, medical imaging,and geosciences literature, but each requiring changes to become
suitable for use in rendering systems. Rendering anisotropic media
also requires expressive and detailed models for anisotropic scattering, which have not been developed in other fields, leading us topropose the new anisotropic scattering model in section 4.

3

Anisotropic Radiative Transfer

For example, consider a volume containing a distribution of aligned
needle-shaped particles with high albedo. Such a volume has much
lower σt (and σs ) for light traveling parallel to the needles thanperpendicular to them. If a volume of these particles is illuminatedfrom the end and viewed from the side, then illuminated from theside and viewed from the end, reciprocity demands that the total
scattering remain the same. But since the amount of light interactingwith the particles is very different, the phase function values must
likewise undergo large changes. Reciprocity still holds in the physical system—its manifestation in the equation is just different. With asomewhat different definition, the phase function can be made symmetrical, but at the expense of no longer being a probability distribution. Again, energy conservation still holds in the physical system,
but its manifestation in the equation is different. We have opted fora normalized fp because it is convenient in the context of Monte
Carlo rendering, and the equation is simpler and more familiar.
The isotropic RTE can be derived straightforwardly by analyzing
the behavior of a volume filled with spherical scattering particles(spherical because their properties are invariant to rotation). Here
we follow this derivation in the style of Ishimaru [1978] but with
non-spherical particles, leading to the anisotropic RTE. In both cases
it’s important to recognize that the resulting equation can also be
used successfully to model media that are not made up of particles—
biological tissue, cloth, hair, foam, etc.

3.1

Isolated non-spherical particle

Consider an isolated particle illuminated by incident radiance L(ω).
We can characterize the particle using three functions:

The isotropic radiative transfer equation (RTE), in the form customarily used in graphics, isZ(ω·∇)L(ω) = −σt L(ω)+σs
fp (ω 0 ·ω)L(ω 0 ) dω 0+Q(ω). (1)S2(In this equation all quantities depend on position, so this dependence is left implicit for compactness.) This equation relates the
directional derivative of the radiance distribution L to its value via
the attenuation coefficient σt and to an integral of L at the same
point via the scattering coefficient σs and the phase function fp .
Energy comes into the equation via the volume source term Q. Thisequation is for an isotropic2 medium because the interaction coefficients σs and σt are independent of direction, and because the phasefunction depends only on the angle between ω and ω 0 .
In order for this equation to describe an anisotropic medium, fp
must be a function of two directions independently, rather than of asingle angle: fp (ω → ω 0 ). Furthermore σs and σt also depend on
ω, leading to the anisotropic RTE:
(ω · ∇)L(ω) + σt (ω)L(ω) =Zσs (ω) fp (ω 0 → ω)L(ω 0 ) dω 0 + Q(ω).

before we can define a valid phase function model or derive a diffusion approximation. In the isotropic RTE, fp is a probability
distribution (when either argument is held fixed) and reciprocity
requires that it be invariant to exchanging its arguments. In theanisotropic case we must choose between keeping fp normalized asa probability distribution and keeping it symmetrical; both cannot
be achieved at once.

(2)

S2

Although this generalization is superficially straightforward, subtle
complications with normalization and reciprocity must be resolved
2 Again, we mean “invariant to rotation,” not “independent of direction.”(i) σ(ω) is the area of the particle’s projection onto ω ⊥ .
(ii) 0 ≤ α(ω) ≤ 1 is the albedo of the particle when illuminatedfrom direction ω.
(iii) p(ω → ω 0 ) is the phase function exhibited by the particle
when illuminated from direction ω and forms a probability
density in the outgoing direction ω 0 .
When the particle receives collimated light in direction ω, σ(ω) is
proportional to the probability of it hitting the particle; α(ω) is theprobability that light is scattered rather than absorbed, conditionedon having hit the particle; and p(ω → ω 0 ) is the probability density for scattering to direction ω 0 , conditioned on having scattered.
Note the distinction between the functions fp and p: whereas fp
describes the overall behavior of the statistical medium, p refers to
the properties of one of its constituent particles.
When the particle is illuminated by a radiance distribution L(ω), thescattered intensity is found by integrating the scattered light overincoming directions:Z
I(ω) =p(ω 0 → ω)α(ω 0 )σ(ω 0 )L(ω 0 ) dω 0 .
(3)
S2

Unlike in the case of a spherical particle, for which σ and α areconstant, p(ω → ω 0 ) 6= p(ω 0 → ω) for a non-spherical particle,
and p is generally not a probability density as a function of its firstargument. The appropriate reciprocity principle is that the scatteredintensity to a single direction under collimated illumination remains

invariant under an interchange of the source and receiver. This
causes p to obey the more involved reciprocity relationp(ω → ω 0 )α(ω)σ(ω) = p(ω 0 → ω)α(ω 0 )σ(ω 0 ).

(4)

It is the whole chain of events, not just the choice of scattereddirection, that has to occur with a reciprocal probability.

3.2

Distribution of particles in a volume

The properties of a volume containing many scattering particles
depend on the characteristics of the particles—both their individual
properties and the distribution of particle orientations that is present.
Up to now there has been no restriction on the shape of particles, but
at this point, because it suffices for the model proposed in the next
section, we assume that particles are rotationally symmetric aboutsome axis, so that their orientation can be entirely described by the
direction of the axis. For the same reason we further assume that theparticles are identical.3
Under these two assumptions we can characterize the volume in theneighborhood of a particular point by two quantities:(i) ρ is the density of particles per unit volume.

Figure 3: The effect of attenuation and scattering on radiance traveling through a volume of particles with a distribution of orientations.
where σt (ω) is the anisotropic volume attenuation coefficient:Z
σt (ω) = ρσ(m, ω)D(m) dm(6)
S2

and fa is an inscattering weighting function:
fa (ω 0 ↔ ω) = ρZp(m, ω 0 → ω)α(m, ω 0 )σ(m, ω 0 )D(m) dm.
S2

We also define σ(m, ω), α(m, ω), p(m, ω → ω 0 ), and I(m, ω) to
be the properties of a particle with orientation m.

To express fa in a more convenient form, we define the anisotropic
volume scattering coefficientσs (ω) as the inscattering due to uniR
form unit radiance, or fa (ω 0 ↔ ω) dω 0 , which works out to:Z
σs (ω) = ρα(m, ω)σ(m, ω)D(m) dm.
(7)

As in the isotropic case, the effect on light traveling along a ray canbe seen by thinking of a cylindrical beam of cross-sectional area A(Figure 3) along which radiance L propagates in the direction ω. An
infinitesimal segment ds of the beam contains n = A ds ρ particles,
each of projected area σ(m, ω).

Here we have made use of the reciprocity relation (4) and the normalization of p over the outgoing direction. We then divide fa by
σs , thereby normalizing it over the incoming direction, to obtain thephase function fp :(ii) D(m), a probability density on the sphere, gives the probability for a particle to be oriented in direction m.

Since they are small and well separated, they do not shadow one
another, and together they intercept an intensity equal to the radiancetimes their combined area:» Z–
dIatten (ω)=A ρσ(m, ω)D(m) dm L(ω).
ds
S2

S2

fa (ω 0 ↔ ω)
fp (ω 0 → ω) ==σs (ω)
Zρp(m, ω 0 → ω)α(m, ω 0 )σ(m, ω 0 )D(m) dm.
σs (ω) S 2(8)which inherits from fa the reciprocity property:

As in derivations of the isotropic RTE, the previous assumption
allows particles to be treated independently, and is presumed to holdin the limit of increasingly smaller cylindrical regions. The integral
over m accounts for the presence of particles with a distribution oforientations. At the same time, the particles scatter intensity into
the direction ω proportional to their number times the intensity (3)averaged over orientations:» Z–
dIinscat (ω)
I(m, ω)D(m) dm .
=A ρ
ds
S2

Together the expressions (6), (7), and (8) define the anisotropic
RTE based on the properties of the non-spherical particles that fill
the medium. In the next section we examine a particular type of
simple particle that can be used to define a phase function model
encompassing a range of scattering behaviors.

In the limit, these changes in intensity for our finite volume define thedirectional derivative of radiance (details in supplementary material):4(ω · ∇)L(ω) = − dIatten (ω)/(A ds) + dIinscat (ω)/(A ds)
which can be rearranged asZ

fa (ω 0 ↔ ω)L(ω 0 ) dω 0

(ω · ∇)L(ω) + σt (ω)L(ω) =(5)S23 Retaining the possibility of more asymmetric particles would require
generalizing the ensuing derivation to a distribution over unit quaternions
rather than over unit 3-vectors. A mixture of particle types can be accommodated easily by summing or integrating over the particles’ properties.

σs (ω)fp (ω 0 → ω) = σs (ω 0 )fp (ω → ω 0 ).

(9)

Substituting σs (ω)fp (ω → ω 0 ) for fa (ω ↔ ω 0 ) in (5), and adding
the source term Q, results in the desired form (2).

Micro-flake model

We consider a model based on two-sided specularly reflecting flakes
oriented according to a known directional distribution. We callthis a micro-flake model. There are obvious parallels to the popular family of microfacet models for surface reflection [Cook andTorrance 1982]: as with microfacet models, the flakes are not individually resolved but instead aggregated statistically to capture
the overall optical behavior of the medium. Another similarity is
the use of a half-vector formulation, which follows naturally fromthe description of a specular surface. The freedom of choosing the
micro-flake distribution will later allow the description of various
types of anisotropic effects in the context of volume scattering.

We assume that the flakes are planar, and that the albedo of a single
flake depends only on the inclination of incident illumination. Using
the framework of section 3 to find the phase function for this type ofparticle entails defining the appropriate functions σ, α and p.

(a)First, due to their planarity, the projected area of a flake with normalm onto ω ⊥ is equal toσ(m, ω) = a |ω · m|
where a is the surface area of one side. The above restriction` on the´albedo of a single particle translates into α(m, ω) = α |ω · m| .
That is, α is an arbitrary function of the angle between ω and m that
does not distinguish between sides of a flake. The flakes are ideal
two-sided reflectors, causing the phase function of a single flake to
be a sum of two delta functions associated with the measure dω:(b)p(m, ω 0 → ω) = δω (m, h(ω, −ω 0 ))+δω (−m, h(ω, −ω 0 )) (10)Here, h(ωi , ωo ) := (ωi + ωo ) /kωi + ωo k denotes the half-direction
vector. Note the sign on ω 0 —this direction needs to be flipped dueto the common convention that the incident direction of a phasefunction points inward. Normalization follows since only one of thedelta functions becomes nonzero for fixed incoming ω 0 and m.
To express the delta functions in (10) with respect to particle orientations m and their associated solid angle measure, a correction factor
must be applied [Walter et al. 2007]:p(m, ω 0 → ω) =δm (m, h(ω, −ω 0 )) + δm (m, −h(ω, −ω 0 )).
4|ω · m|and this function satisfies the reciprocity criterion (4), as light can
only be scattered when m · ω = m · ω 0 . Applying the tools from
Section 3 then leads to the following expressions for σt , σs and fp(details in supplementary material):Z
σt (ω) = aρ|ω · m|D(m) dm
(11)
2
ZS
`
´
σs (ω) = aρ
α |ω · m| |ω · m|D(m) dm(12)S2

fp (ω 0 → ω) =aρ
4 σs (ω)`´
D(h(ω, −ω 0 )) + D(−h(ω, −ω 0 ))`´
α |ω 0 · h(ω, −ω 0 )|(13)

In the special case of a uniform flake distribution D(m) = 1/4π,
this model reduces to traditional “anisotropic” scattering, where
the phase function depends only on the scattering angle: here, thescattering coefficient σs reduces to a constant related to the albedo
Z π2
σs = aρα(cos θ) cos θ sin θ dθ,0the extinction coefficient turns into σt = aρ, and fp becomes a2function of cos θ = ω · −ω 0„«aρθfp (cos θ) =
α cos
.
(14)
8πσs
2
which can be solved for the α required to archive a particular fp :α(cos θ) =´´
8πσs `σs `fp cos(2θ) = 4π
fp cos(2θ)
aρ
σt

(15)

For instance, to represent the Henyey-Greenstein or Mie phase functions, equation (15) provides a way of encoding their properties into
an albedo function α. Any isotropic medium can be expressed inthis fashion, even when the material itself is not physically realizible

Figure 4: Two sample flake normal distributions and reflection profiles corresponding to the indicated incident direction ωi :(a) D(ω) = cos20 (v, ω)+ : Reflection similar to a rough surface,where v acts as surface normal. (b) D(ω) = sin20 (v, ω): Reflection similar to a rough fiber, where v acts as the fiber’s tangent
direction. In both cases α is constant.
using flakes. For instance, by the previous equation, a medium without absorption and strong forward scattering will produce α > 1 insome directions. Although the physical interpretation is lost in suchcases, the framework remains internally consistent and all results
hold regardless. This situation is similar to microfacet models, which
can often still provide convenient abstractions when the underlying
material breaks the statistical heightfield assumption.
Various types of anisotropic scattering can be modeled dependingon the choice of flake distribution. For instance, concentrating
D(m) around a fixed direction (Figure 4a) produces behavior that
resembles surface reflection, and in fact the form of the model is notunlike the commonly used Blinn-Phong BRDF, with the term σt (ω)
playing the role of the usual cosine factor for surface illumination.
Similarly, it is also possible to produce fiber-like scattering by choosing a D(m) that is concentrated near the plane perpendicular to thefiber’s tangent direction (Figure 4b). In this case, the resulting phasefunction resembles the reflection model of Kajiya and Kay [1989],
with σt (ω) playing the role of the usual sine factor for fiber illumination. Some differences exist; for instance, the microflake model
produces a reflection cone that narrows towards forward scattering,
resulting in less forward scattering than Kajiya-Kay if α(ω) is constant. However, the micro-flake approach does not share some ofthat model’s shortcomings, particularly its lack of energy conservation: especially for illumination near grazing angles, the Kajiya-Kaymodel is far from normalized. Attempts to remedy this problem
using a parameterized normalization constant invariably lead to a
non-reciprocal phase function. The scattering model presented here,on the other hand, is always energy conserving, and it satisfies
the appropriate reciprocity constraint (9) in conjunction with thepreviously defined directionally varying scattering coefficients.
Section 6 will present three solution techniques for rendering withthis model, including a Monte Carlo approach and two approximate
methods that involve the diffusion equation.

5

Anisotropic diffusion

To reason about diffusion in the presence of anisotropy, we must
propagate the earlier changes to the RTE through the steps of thisapproximation. Classic derivations of the diffusion equation (e.g.
Ishimaru [1978]) set the stage by considering solutions to the RTE
under smoothness constraints.
The key assumption is that because each scattering event effectively
blurs the incident illumination, high frequencies disappear fromthe angular radiance distribution as the number of scattering events
increases. Motivated by this observation, the radiance function
is restricted to a simple two-term expansion in terms of angularmoments. Formally, for a function f : S 2 → R, the n-th momenton the unit sphere may be defined asZ(µn [f ])i,j,k,... :=ωi ωj ωk · · · f (ω) dω.
(16)
{z }
S2 |

V is a 3-dimensional vector and M a 3 × 3 matrix, which both
capture information about the low-frequency scattering behavior. In
the case of isotropy, these equations together appropriately reduce to(20), and M −1 = I/(3σt0 ) takes the role of the diffusion constant.
This approximation may be applied to any reflection model expressed in the framework of Section 3, and will be accurate to theextent that radiance does become uniform in the region of interest.
In Section 5.3 we will further specialize to the micro-flake model.

5.2

Boundary conditions

We can proceed to derive an anisotropic generalization of the commonly used Robin boundary conditions: these specify that the net
inward flux Γin (x) matches the fraction Fdr of the outward flux thatis reflected by interface effects Γref (x), i.e.
Γin (x) = Γref (x)n factors

For non-negative f , the 0th moment for instance gives the function’s
integral over the sphere, the 1st moment can be interpreted as a“center of mass” 3-vector, and the 2nd moment is a positive definite
3×3-matrix. The above restriction then translates into settingL(x, ω) :=3
1
φ(x) +
ω · E(x)
4π
4π

where
Γin (x) =Γref (x) =Fdr

The traditional derivation then proceeds by substituting this approximation into the RTE and enforcing equality amongst the 0th and1st -order projections, which leads to the coupled set of equations∇ · E(x) = Q0 (x) − σa φ(x)(18)

∇φ(x) = 3Q1 (x) − 3σt0 E(x)(19)

where Qi (x) := µi [Q(x, ·)], σt0 := σa + σRs (1 − g) and the constantg specifies the mean cosine defined as g := 4π fp (ω·ω 0 )(ω·ω 0 ) dω 0 .
Eliminating E in (19) produces the isotropic diffusion equation
λ∇ φ(x) = σa φ(x) − Q0 (x) + 3λ∇ · Q1 (x)L(x, ω) |ω · n| dω(24)

L(x, ω) |ω · n| dω

(25)

2−

1+Fdr
where we have defined A := 1−F. For mismatched boundaries,
dr
we make use of the following rational approximation to measured
data due to Egan and Hilgeman [1979]:(20)

where λ := 3σ1 0 is the diffusion constant.

ZH

H 2− indicates integration over the inward hemisphere and the normal n is assumed to point outwards. Substituting the two-term
expansion into (23) leads to
»
–
1
1
1
1
φ(x) − n · E(x) = Fdr
φ(x) + n · E(x)
4
2
4
2
1
A
⇔
φ(x) − n · E(x) = 0
(26)
4
2

Fdr = −

2

Z

H 2+

(17)

where the quantities φ(x) = µ0 [L(x, ·)] and E(x) = µ1 [L(x, ·)] are
known as fluence and vector irradiance, respectively.

(23)

1.440
0.710
+
+ 0.668 + 0.0636η
η2
η

(27)

where η is the relative refractive index at the boundary.

t

5.1

Anisotropic diffusion—general form

The same set of ideas can also be applied to the anisotropic RTE.
Here, the angular dependence of the scattering coefficients andthe more general form of fp causes the appearance of additionalterms in the final set of equations. Since performing the two projections amounts to a mechanical application of rules also used in
the isotropic derivation, we will only state the final result; the fullderivation may be found in the supplementary material to this paper.
The first and second-order projections of equation (2) are given by
1∇ · E(x) = Q0 (x) −µ0 [σa (x, ·)] φ(x) − V (x) · E(x) (21)4π
3µ1 [σa (x, ·)] φ(x) − M (x) E(x) (22)∇φ(x) = 3Q1 (x) −
4π
where
»
–
Z
3
µ1 [σt (x, ·)] −σs (x, ω)µ1 [fp (· → ω)] dω
4π
S2
»
–
Z
9
M (x) :=
µ2 [σt (x, ·)] −ω σs (x, ω) µ1 [fp (· → ω)]T dω
4π
S2
V (x) :=

5.3

Anisotropic diffusion by micro-flakes

To apply the diffusion approximation to the micro-flake model, thecoefficients in (21) and (22) must be computed. In this section, we
derive these coefficients, discovering some important simplifications
for micro-flakes compared to the general case.
Beginning with the coefficient of φ in (21), we findµ0 [σa (x, ·)] = 2aρπ(1 − α̂)
where we have reversed the order of integration and defined
Z π
`
´
α̂ :=
α | cos θ| |cos θ| sin θ dθ.

(28)

(29)

0

This constant may be interpreted as the overall reflection from aflake under constant illumination. When the albedo does not varywith respect to inclination, α = α̂.
An immediate observation is that all first-order moments of thescattering and absorption coefficients vanish due to their mirrorsymmetry (e.g. σa (x, ω) = σa (x, −ω)). This eliminates the firstterm in the expression for V. The remaining vector-valued integralexpression in V (x) also vanishes: after a change of variables tohalf-directions, we have
Z
Z
`
´
3aρ
V (x) =
D(x, m) ω 0 α |ω 0 · m| |ω 0 · m| dω 0 dm = 0
4π S 2
S2

Finally, with the V (x) term gone, the equation also takes on a morefamiliar form: solving (22) for E(x) and substituting into (21) as isdone in the isotropic variant yields

which follows from the symmetry of the inner integrand (i.e.
f (ω 0 ) = −f (−ω 0 )).

(35)

ˆ
˜
∇ · M (x)−1 (∇φ(x) − 3Q1 (x))
1µ0 [σa (x, ·)] φ(x) − Q0 (x).
=4π

By the same argument used for the first term of V, the coefficient of
φ in (22) vanishes.

This step can be repeated once more to eliminate E(x) from theboundary conditions (26):Turning finally to the matrix coefficient M, the first term in its
definition takes on a simple form relating it to the low-frequency
properties of D(m):1A
φ(x) + nT M −1 (∇φ(x) − 3Q1 (x)) = 0.
42aρπµ2 [σt (x, ·)] =(I + µ2 [D(x, ·)]) .
2(30)

The second term can be handled by a change of variables to halfdirections and some simplifications, and then M is found to havea surprisingly simple closed form as a function of the second flake
distribution moment:
M (x) =´
9aρ `(1+2c1 +2c2 )µ2 [D(x, ·)] + (1−2c1 )I .
8(31)

where the constants
Z π
c1 :=
α(| cos θ|) |cos θ| sin3 θ dθ0Z πc2 := 2α(| cos θ|) |cos θ|3 sin θ dθ.

(32)
(33)

0

only depend on the choice of the albedo falloff curve. A property thatwill later be required is the positive definiteness of M (x), which directly follows from c1 , c2 ≥ 0, c1 ≤ 12 and the positive definitenessof µ2 [D(x, ·)].
When α is constant, c1 = α/2, c2 = α and (31) further reduces to
M (x) =´9aρ `(1 + 3α)µ2 [D(x, ·)] + (1 − α)I .
8(34)

The simplicity of (31) is of particular convenience when the flake distribution D(m) is available in the form of a real spherical harmonics
expansion: suppose that
D(x, ω) =

∞
X

l
X

where
8
0
>
(ω),
if m = 0<Y√lmyl (ω) =2 Re Ylm (ω),
if m > 0>:√2 Im Y −m (ω), if m < 0l
and Ylm are the usual complex-valued orthonormal spherical harmonic basis functions. Then µ2 [D(x, ·)] is readily obtained in termsof the expansion coefficients:µ2 [D(x, ·)] =
2q

2

π
15

6
6
6
6
4

a0
5 a0 − √2 + a2
2
3 0
3
−2
2

−2
2
q
a0
5 a0 − √2 − a2
2
3 0
3

−a1
2

−1
−a
2

a

−a1
2

a

−1
−a
2

q

The two diffusion-based solution techniques in Section 6 are built
on equation (35) and (36) and they only require that the problem
be in this form. This means that as long as M is a positive definite
matrix and µ0 [σa (x, ·)] a nonnegative constant, they can be appliedirrespective of what model was used to generate them.
To summarize: we have derived a diffusion equation and suitableboundary conditions based on first- and second-order projectionsof our anisotropic RTE, which resulted in the introduction of anadditional vector term V and the matrix M . For the micro-flakemodel, the V term reduced to zero and M took on a simple closed
form solution (31) relating it to the low-frequency characteristics of
the underlying flake distribution. With all of these pieces at hand,
we are now able to compute approximate diffusion-based solutions
to the anisotropic RTE when micro-flakes are used as the underlying
phase function model.

6

Solution techniques

In this section, we present three different approaches for computing
light transport in anisotropic media, spanning a range of performance and accuracy characteristics. First, we describe the necessary changes to support the flake model in a Monte Carlo renderer.
Next, we adapt the finite element method to the anisotropic diffusion equation found in Section 5. And finally, we propose a new
dipole BSSRDF for anisotropic translucent materials. Both diffusion
approaches are general in the sense that they only depend on thecancellation of the V (x) and µ1 [σa (x, ·)] terms in (21) and (22),which means that they may also be applied to any future scatteringmodels sharing these properties.

6.1

Monte Carlo rendering

To use micro-flake distributions in a Monte Carlo rendering context,
an efficient way of evaluating the coefficients σs , σt , and the phasefunction fp will be necessary. Furthermore, good importance sampling support for fp is crucial to achieve reasonably fast convergence,in particular when working with scatter-dominated media.

maml yl (ω)

l=0 m=−l

r

(36)

3
7
7
7
7
0 5

a
5 a0 + 2 √2
3 0
3

For simple families of flake distributions, it may be feasible to findclosed-form importance sampling procedures and expressions for σs
and σt . For general D(m), however, we face a prohibitive numerical
2D integration for each scattering coefficient evaluation, since (11)and (12) are expressed in terms of spherical convolutions.
Both problems are greatly alleviated by moving to the frequency domain: for this, we represent D(m) as a truncated spherical harmonics expansion using enough coefficients to ensure that the introduced
error is negligible (< 0.1%). The convolutions in the expressions
for σs and σt are then found to involve
invariant` the azimuthally´
kernel functions |ω · m| and |ω · m| α ω · m , making it easy tofind the coefficients of their respective expansions by an applicationof the Funk-Hecke theorem [Groemer 1996].

To support importance sampling of the phase function (13), we perform a 2D discretization over the outgoing direction: for each fixedoutgoing ω, we tabulate the coefficients for a truncated sphericalharmonics expansion of fp (· → ω). During rendering, an expansion
over the incident direction is obtained using bilinear interpolationfrom the nearest four outgoing direction samples. In all our examples we use between 7 and 10 coefficient bands and a 50 × 100
discretization over the (θω , φω )-space.
With the spherical harmonics coefficients at hand, it is possible toto apply the sampling technique proposed by Jarosz et al. [2009],
which we found to work well in practice. Any procedure of thisform of course introduces truncation and interpolation errors, andwe therefore use the interpolated phase function only to guide theplacement of samples, with (13) taken as the actual phase function
during evaluation.

6.2

Anisotropic Diffusive Source Boundary Condition

The ADSBC combines the basic Robin boundary condition with a
new term describing illumination incident on the boundary due tolight sources outside of the medium. This is often referred to as the
reduced intensity source. This illumination is approximated by adiffusive flux term equal to the net refracted irradiance attenuatedby the expected absorption before the first scattering event. Thischanges (23) to Γin (x) = Γref (x) + Γs (x) whereZ
σ (x,n)− aΓs (x) = e σs (x,n)
Ft (ω)L(x, −ω)(n · ω) dω
(37)
H 2+

Following the mathematics presented in Section 5.2 ((24)–(36))
results in the anisotropic diffusive source boundary condition:
1
A
1
φ(x) + nT M −1 ∇φ(x) =Γs (x)
42
Fdt

Anisotropic Finite Element Diffusion Equation

Next, this boundary condition is used to find the finite element
solution of (22). This process has two steps: the derivation of the
weak form and the discretization of that form into a matrix equation.
Since these steps are essentially the same as in [Arbree et al. 2009],
with (22) replacing the original isotropic equation, we present only
the final matrix equation here:(R + S + T )a = (q + g)

gi =

2
AFdt

Z
Γs βi ds∂ΩThe βi are basis functions defined at the vertices of the P
mesh and adenotes the solution’s coefficient vector satisfying φ = n−1i=0 ai βi .
Since M −1 is always positive definite (Section 5.3), this equation is
guaranteed to have a solution.
The final step is the derivation of the query function that converts
the solution fluence of the matrix equation into the outgoing surface
radiance (see [Jakob et al. 2010]). For completeness, we replicate
this equation, which turns out to be the same as the isotropic case.
L(x, ω) =Ft (ω)
4π

»„
1+

1
A

«
φ(x) −

4
Fdt A–
Γs (x)

6.3

BSSRDF model

In this section, we will first solve the anisotropic diffusion equation for a point source located in an infinite homogeneous medium.
Following this, we use our solution to extend the BSSRDF dipole
model proposed by Jensen et al. [2001] to anisotropic media.
6.3.1

Monopole

To arrive at the monopole solution, we model a point source locatedat the origin radiating the total amount of power Φ uniformly into
all directions:Φδ(x).
Q(x) =
4π
For this type of source, Q1 (x) = 0 which further simplifies theanisotropic diffusion equation (35). The mixed derivatives on theleft-hand side of this equation can be eliminated by switching toa set of linearly transformed coordinates ξ = P x. The particularmatrix P which accomplishes this is given by1

P := QD 2 QT

∇2 φ(ξ) =1
φ(ξ) µ0 [σa (·)] − Q0 (ξ)
4π∇βi · (M −1 ∇βj ) dxΩ

After a similar change of variables in the source term, its 0th -orderprojection is found to be Q0 (ξ) = Φ |det P | δ(ξ). Finally, theGreen’s function describing the fluence due to this type of source isfound for example in Ishimaru [1978]:Zqi =Q0 βi dx
Ω

(42)

(39)

where
Rij =

(41)

where QDQT is an eigendecomposition of M . This intuitively
means that the anisotropic diffusion equation reduces to the isotropic
variant when expressed using the appropriate coordinates:!
r1
Φ det Pφ(x) =exp −µ0 [σa (·)] kP xk .
4πkP xk
4πZ(40)

Mathematically, this final solution is nearly identical to the isotropic
case. The two key differences are the new definitions of the matrices
R and S in (39).

(38)

where Fdt := 1 − Fdr denotes the overall Fresnel transmittance.
6.2.2

Z
1
µ0 [σa ] βi βj dx
4π Ω
Z
1
Tij =
βi βj ds
2A ∂Ω

Sij =

Finite Element Solution

The finite element method (FEM) provides a faster solution technique for the anisotropic diffusion equation. Our finite element
solution has three components: first, we modify the Robin boundary
condition (23) to produce a so-called anisotropic diffusive sourceboundary condition (ADSBC). Second, we use this new boundary
condition to compute a matrix solution to the anisotropic diffusion
equation (35). Finally, we present a query function that converts
the computed solutions into radiance. Since our derivation closely
follows the isotropic case [Arbree et al. 2009], we highlight only thedifferences here. More details can also be found in the supplementary material [Jakob et al. 2010].
6.2.1

and

which is expressed back in terms of the original coordinates.

(43)

n

φ(x) ≈ 0

P xv

nT M −1 n
kM −1 nk2

2AkM −1 nk2

2AnT M −1 nxv

p1 × p2

Figure 5: Illustration of the approximate boundary conditions usedby the anisotropic BSSRDF model.
xr

6.3.2

Dipole boundary conditions

We now switch to a half-space geometry with the interface locatedat z = 0 and a surface normal given by n = (0, 0, 1)T . To describe
this setup, the dipole model [Farrell and Patterson 1992] relies on
an approximation to the Robin boundary conditions presented earlier: in the isotropic setting, Moulton [1990] has shown that (36)
is approximately satisfied by requiring the fluence to be zero on an
extrapolated boundary at distance 2Aλ (where λ is the diffusion
constant).
To find the corresponding boundary condition at the interface of ananisotropic medium, we can interpret (36) as a directional derivativeof φ(x) along the direction (M −1 n)/kM −1 nk, i.e.
„«
M −1 n·∇φ(x) = −φ(x).
2AkM −1 nk
kM −1 nk
When this approximation is good, we expect φ to vanish at a distanceof 2AkM −1 nk in this direction. Since x was an arbitrary pointon the surface, this property will in fact hold on a whole surface
extrapolated at a distance of(b) Transformed space

Figure 6: Geometric intuition of the anisotropic dipole: under alinear transformation, the diffusion becomes isotropic, which allows
the straightforward construction of a dipole by mirroring the realsource across the extrapolated surface (dashed line).
q1
where β :=µ [σ (·)] and dr and dv denote distances com4π 0 aputed in a space transformed by P , i.e. dr := kP (x − xr ) k anddv := kP (x − xv ) k.
The final BSSRDF is obtained by computing the radiant exitancein the normal direction, n · E(x). Applying these steps leads to thefamiliar-looking"#
ᾱ(n) det P
e−β dre−β dvdMo (x) = dΦC1
+ C2
, (47)
4π
d3r
d3v
where

dp = 2A nT M −1 nwhich can be obtained by projecting the vector M −1 n/kM −1 nk
onto n (Figure 5).
6.3.3(a) Original space

P xr

Dipole solution˛˛˛˛
C1 := ˛nT xr ˛ (βdr + 1) and C2 := ˛nT xv ˛ (βdv + 1)

(48)

As in Jensen et al. [2001], we have added an albedo term ᾱ(n) =σs (n)/σt (n) to account for the first scattering event in the earlier
collimated beam approximation.

Following Farrell et al. [1992], an approximate half-space solution
is now found based on a dipole source configuration. In this model,
a negative ‘virtual’ light source is introduced to satisfy the zero
fluence constraint on the extrapolated boundary.

It is possible to further extend this model to multipole source configurations similar to Donner and Jensen [2005]. This is straightforwardand entails simply adding more sources to cancel out the fluence on
two extrapolated boundaries in the skewed solution space.

In the isotropic case, placement of the sources is simple: the realsource is positioned one mean free path below the surface (zr =−1/σt0 ) to simulate the first scattering event due to a collimatedbeam at normal incidence. The zero fluence constraint at z = 2Aλis satisfied by mirroring the real source along this plane, which leads
to a virtual source position of zv = 1/σt0 + 4Aλ.

7

The natural extension to anisotropic media places the real source at
zr = −1/σt (n) below the surface, which matches the mean freepath at normal incidence. Determination of the virtual source position is slightly more complicated since the monopole fluence (43)is now a function of kP (x − xp )k. Under the linear transformation
P , the extrapolated surface z = dp turns into a plane with a normal
pointing into direction p1 × p2 , where pi are the columns of P (Figure 6). Mirroring the light source along this plane and transforming
the obtained position back into the original space results in
xr := −σt (n)−1 n,`´
xv := xr + 2 σt (n)−1 + dp

(44)
−1

P (p1 × p2 )nT (P −1 (p1 × p2 ))The fluence due to both sources is then equal to»–Φ det P e−β dr
e−β dvφ(x) =
−
4π
dr
dv

(45)

(46)

Results

Our final rendering pipeline, written in C++ and Java, ran on an
2.93Ghz Intel Xeon X5570 workstation with 8 cores. We implemented all three solution techniques, and this section presents an
overview of our experimental results.
The scarf model was kindly provided by theauthors of the recent work on yarn-level knit simulation [Kaldor
et al. 2008]. To create volume densities based on their spline representation, we used a technique similar to the Lumislice algorithm[Xu et al. 2001], but took advantage of the flexibility of a generalvolume to build 2-ply yarns (yarns made by twisting multiple smalleryarns called “plies”). The fibers were well organized, so they should
produce specular highlights running along the twisted plies. The
final volume had an effective resolution of 1408 × 256 × 1152 voxels, which required about 3 gigabytes of storage using a sparse grid
data structure. Every occupied cell in the volume consisted of 4
single-precision values specifying density and local fiber orientation,and the latter was used both as an input to the Kajiya-Kay model,
and to rotate flake distributions into the desired coordinate frame.
Scarf renderings:

Three Monte Carlo simulations with full multiple scattering were
performed at 512 samples per pixel using an isotropic phase function

(a) Micro-flakes (SS)(b) Kajiya-Kay (SS)(c) Kajiya-Kay (MS)(d) Lost energy (MS)Figure 7: (a), (b): Micro-flake and Kajiya-Kay renderings with exponents tuned for similar highlight sharpness in single-scattering images—
already here, the lack of energy conservation in the Kajiya-Kay model is visible. (c): Kajiya-Kay rendering with full multiple scattering. (d):difference image showing the energy loss incurred in (c) as compared to the micro-flake rendering in Figure 1b.

(a) Illuminated from the left(b) Illuminated from the right(c) Collimated illumination

Figure 8: (a), (b): Using spatially varying flake distributions oriented according to the direction maps captured by Marschner et al. [2005],
we are able to reproduce the shifting reflection patterns (see the accompanying video) exhibited by their model, while furthermore accounting
for multiple scattering. (c): Projecting collimated illumination onto the wood reveals interesting diffusion effects along the grain direction.
(Figure 1a), micro-flakes (Figure 1b) and the Kajiya-Kay model(Figure 7), respectively. The isotropic rendering lacks any speculareffects, leading to dull and somewhat blurry appearance, while theanisotropic renderings show highlights delineating the plies. In theanisotropic case, multiple scattering naturally leads to near-whitehighlights (caused by single scattering), and in all cases it leads
to increasing color saturation in the body of the material and theninto the shadows, due to the many scattering events required forenergy to reach these regions. The rendering using the Kajiya-Kay
phase function shows similar local shading to the micro-flake model,
since the two phase functions are similar. An important differenceis that the micro-flake model is exactly normalized, whereas the
Kajiya-Kay model loses energy towards grazing angles. A fractionof the total light energy is lost at every bounce, and this causes
blurring in multiple scattering renderings, as well as a noticabledecrease in the overall image brightness (Figure 7d). As theseexamples demonstrate, accounting for oriented structures leads to asignificant gain in visual realism, and the new model achieves thiswhile retaining the full advantages of physics based rendering.
Changing from isotropic scattering to the micro-flake model whenrendering heterogeneous materials does incur additional computational costs. Most of the extra time is spent evaluating σt (ω) at
every step of ray marching between scattering events. In our code
this currently involves evaluating a spherical harmonics expansion,leaving considerable opportunities for future optimization.
As an example of a solid material withanisotropic structure, we used the measurements of finished wood
released by Marschner et. al. [2005], which contain both diffuse
Wood example:

color and fiber direction textures for several types of wooden materials. To use these textures in a volume rendering context, thesingle scattering albedo α needs to be known, which we computedby numerically inverting the dipole BRDF approximation [Jensen
et al. 2001] for each pixel of the diffuse color texture.
We then projected the curly maple textures into a slab filled with aflake medium and chose a distribution that simulates scattering fromrough fibers (D = sin20 θ). In the resulting MC renderings (Figure 8), the flake medium is able to reproduce the shifting anisotropic
reflections modeled by the original paper, and the ability to accountfor multiple scattering eliminates the need for an ad-hoc diffuse
component. By illuminating the wood with a a grid of collimatedbeam sources, we can observe anisotropic light diffusion occurring
along the grain direction.
To validate the finite element computation against Monte Carlo, we used a simple scene to test the
spreading of light from an illuminated spot. Because collimated illumination and half-spaces cannot be represented in the FEM system,
we projected a 1mm diffuse spot on a 4×4cm cube with material
parameters aρ = 5.2482/cm and α = 0.998. Except for the spot’s
immediate neighborhood, where the simulations are not expected
to match because of the large contribution of low-order scattering,
the solutions were found to be in good agreement (Figure 10a–c).
Because of the structured sources in our examples, we modified themesh refinement algorithm to increase resolution along edges in the
projected source image. Since the finite element solution does not
require any particular geometry or a homogeneous material, it isable to model a wide range of scenes, such as Figure 9a.
Finite element simulations:(a) FEM solution to heterogeneous anisotropy(b) Isotropic BSSRDF model by Jensen et al.

(c) Our new anisotropic BSSRDF model

Figure 9: (a): Finite element solution to the Stanford Dragon illuminated by a grid of circular spots. The model is filled by an anisotropic
micro-flake medium with its flake distribution aligned to a procedurally generated vector field, causing illumination to spread along curved
paths. (b), (c): Comparison of the isotropic and anisotropic BSSRDF: a light source projects colored illumination on a translucent slab. Bothuse the same albedo and average density, but the anisotropic slab uses a sin20 flake distribution about a diagonal axis parallel to the surface.
Scene
Scarf model (Kajiya-Kay)
Scarf model (Micro-flakes)
Heterogeneous dragon
Isotropic slab (1 sample/pixel)
Anisotropic slab (1 sample/pixel)

Algorithm Rendering time
MC
MC
FEM
DipoleDipole5 hr 24 min
22 hr 11 min
6 min 26 sec
3.4 sec
9.8 sec

Table 1: Timing results for some of our example scenes
Anisotropic dipole: We have added support for anisotropy to
our existing implementation of the dipole BSSRDF [2001], which
required only minimal changes to the code. Figure 9b–c shows
the result of projecting colored illumination onto an anisotropic
slab modeled in this way. In any dipole model renderer, a heuristic
for transforming points into the half-space geometry is required,and in the anisotropic case a simple solution based on distance isinsufficient, since dMo is now a function of two parameters insteadof just the radial distance. To keep our implementation simple,
we simply drop one coordinate, and this can be seen to introduce
artifacts on the sides of the slab. We believe that a heuristic withoverall smoother behavior could be used to remove such artifacts.

Similarly to the FEM test case, we constructed a simple scene tovalidate the dipole model, and we ran Monte Carlo simulations ofa collimated beam impinging on an anisotropic half-space with thesame parameters as before. Figure 10d–f shows a result of thiscomparison; as in the isotropic case, the dipole is only found to beaccurate in the far field. A curious observation in our simulations
was that the lower-order scattering caused by the very concentrated
beam source tends to propagate perpendicularly to the main directionof diffusion (Figure 10e), which was also noticed by Kienle et al.
[2004]. The dipole does not attempt to model this effect, which
accounts for some of the differences. We do not see this as an
obstacle for graphics applications, since the lower-order scattering
effect quickly disappears as one moves to finite beam widths.
Another noteworthy difference is that the anisotropic BSSRDF may
not achieve its maximum value at the point of incidence, whichhappens when increased diffusion occurs along an axis that is tilted
with respect to the surface plane. This is not neccessarily a failureof the dipole, but a materialization of the disagreement between
the embedded point source and collimated beam solutions in suchmedia. When this effect is not desired, the dipole solution can be
re-centered to remove the inconsistency. A more extensive overview
of the Monte Carlo comparisons can be found in the accompanying
technical report [Jakob et al. 2010].

8

Conclusion

This paper has provided the theory and practice required to integrate
anisotropic media into physically based rendering systems. An
important piece of this is the introduction of a new model for media,
based on micro-flake distributions, which generalizes previous phasefunction models and provides a physically sound alternative to themore heuristic reflectance models that are normally used in volume
models of materials with directional scattering properties.
We also have shown how the more general anisotropic radiativetransfer equation still can be approximated by diffusion (thougha diffusion matrix is required where a scalar diffusion constant
was formerly used), which provides fast, approximate solutions
to multiple scattering in highly scattering media, by using a finite
element procedure or by using a more approximate but extremely
fast anisotropic dipole solution. The diffusion approximation shares
the limitations familiar from the isotropic case: it applies only whenradiance becomes uniform in direction, which happens only in highalbedo media and only far from the source—and these are satisfied
more quickly the lower the degree of anisotropy.
The framework introduced in this paper provides a basis for future
work in several areas. First of all, new volume models for complex
geometry at any scale, from cloth viewed at arm’s length to a forest
viewed from across a valley, can now be built based on physicallyrealistic scattering models, with multiple scattering naturally adding
richness to the appearance, rather than being built as volume aggregations of surfaces shaded by heuristic single-scattering models.
Volume models are also a natural setting for multiresolution approaches [Neyret 1998]. Aggregating geometry to larger scales canlead to anisotropic behavior even when the volume is isotropic at thesmall scale, so our anisotropic framework provides an ideal platform
for building volume level-of-detail in a physically based setting.
The framework we have presented does have limitations. One iscomputational cost: volume models consume significant memory,and the more so when they are augmented with spatially varying vectors or other parameters for anisotropy; also, evaluating anisotropic
models consumes more computing time than looking up scalar parameters. The radiative transfer formulation itself has limitations,and may not correctly model short-distance effects below the correlation length of the random medium being simulated. The diffusion
approximation has further limitations; it only can be expected tobe accurate when the radiance distribution becomes directionallysmooth, which happens less easily with highly anisotropic or sharply
heterogeneous media, and the dipole solution, while it gives reasonable estimates in the far field for half-space geometry, is (just as in(b) MC(c) FEM(e) MC(f) Dipole

1

0.1
0.01

0.1

0.001

0.01

10

4

0.001

10

5

10

6

2

1

0

(a) D(ω) = 1 − (x · ω)2

1

2

cm

10

4

4

2

0

(d) D(ω) = (x · ω)4

2

4

cm

Figure 10: Validation of the finite element and dipole solutions against Monte Carlo references for two different distributions D(ω) centeredaround the x-axis. The blue and red graphs represent horizontal and vertical slices through the origin, respectively. (a)–(c): The finite element
diffusion solution agrees well with the Monte Carlo simulation away from the source. (d)–(f): The dipole solution is more approximate but
provides a useful match far from the source. The Monte Carlo simulation shows a characteristic feature of anisotropic media: low-orderscattering produces horizontal spread near the center and high-order scattering produces vertical spread in the diffusion region.
the isotropic case) only a very rough model when applied to more
complex surfaces.
Also, though the micro-flake model has considerable flexibility andcan capture the basic phenomena of surface-like, fiber-like, andpoint-like scattering, not all volume scattering functions can be
represented by micro-flakes. Future work will need to discover thepractical significance of its limits in representing realistic materials,
developing extensions where needed. Aggregating given geometry
into volume models (another form of volume level of detail) byselecting appropriate micro-flake distributions is also an important
potential application.
The most important implication of our work for the future of volume models is that they can now act as full-fledged participants inglobal illumination systems. This means volume models, whichhave proven so successful for fur, cloth, trees, and other complex
structures, can take full advantage of the richness and realism that isnaturally provided by complete light transport simulations.

9

Acknowledgements

This work was supported by the National Science Foundation (grants
CCF-0347303 and CCF-0541105) and by Unilever Corporation. The
authors also thank Jonathan Kaldor, Manuel Vargas and Manolis
Savva for providing the scarf dataset.

References
A RBREE , A., WALTER , B., AND BALA , K. 2009. Heterogeneous
subsurface scattering using the finite element method. To appear
in IEEE Transactions on Visualization and Computer Graphics.
A RRIDGE , S. R. 1999. Optical tomography in medical imaging.
Inverse Problems 15, 2, R41–R93.

A SHIKHMIN , M., P REMOZE , S., AND S HIRLEY, P. S. 2000. A
microfacet-based BRDF generator. In Proceedings of ACM SIGGRAPH 2000, 65–74.
B LINN , J. F. 1982. Light reflection functions for simulation ofclouds and dusty surfaces. In Computer Graphics (Proceedingsof SIGGRAPH 82), 21–29.
C EREZO , E., P ÉREZ , F., P UEYO , X., S ERON , F. J., AND S ILLION ,F. X. 2005. A survey on participating media rendering techniques.
The Visual Computer 21, 5, 303–328.
C HANDRASEKHAR , S. 1959. Radiative Transfer. Oxford University
Press.
C OOK , R. L., AND T ORRANCE , K. E. 1982. A reflectance model
for computer graphics. ACM Trans. Graph. 1, 1 (Jan.), 7–24.
C RASSIN , C., N EYRET, F., L EFEBVRE , S., AND E ISEMANN , E.
2009. Gigavoxels : Ray-guided streaming for efficient and detailed voxel rendering. In ACM SIGGRAPH Symposium on Interactive 3D Graphics and Games (I3D).
D ONNER , C., AND J ENSEN , H. W. 2005. Light diffusion in
multi-layered translucent materials. In SIGGRAPH ’05: ACM
SIGGRAPH 2005 Papers, ACM, New York, NY, USA, 1032–
1039.
D REBIN , R. A., C ARPENTER , L., AND H ANRAHAN , P. 1988.
Volume rendering. In Computer Graphics (Proceedings of SIGGRAPH 88), 65–74.
D UDKO , O. K., AND W EISS , G. H. 2005. Estimation of anisotropic
optical parameters of tissue in a slab geometry. BiophysicalJournal 88, 5, 3205 – 3211.
E GAN , W., AND H ILGEMAN , T. 1979. Optical Properties of
Inhomogeneous Materials. Academic Press New York.
FARRELL , T. J., AND PATTERSON , M. S. 1992. A diffusion theorymodel of spatially resolved, steady-state diffuse reflectance forthe noninvasive determination of tissue optical properties in vivo.
In Medical Physics, Volume 19, Issue 4, 879–888.
F EDKIW, R., S TAM , J., AND J ENSEN , H. W. 2001. Visual simulation of smoke. In Proceedings of ACM SIGGRAPH 2001,
15–22.
G IBSON , A. P., H EBDEN , J. C., AND A RRIDGE , S. R. 2005.
Recent advances in diffuse optical imaging. Physics in Medicine
and Biology 50, R1–R43.
G ROEMER , H. 1996. Geometric applications of Fourier series andspherical harmonics. Cambridge Univ Press.
H EINO , J., A RRIDGE , S., S IKORA , J., AND S OMERSALO , E. 2003.
Anisotropic effects in highly scattering media. Phys. Rev. E 68, 3.
H EISKALA , J., N ISSILA , I., N EUVONEN , T., JARVENPAA , S., ANDS OMERSALO , E. 2005. Modeling anisotropic light propagation
in a realistic model of the human head. Applied Optics 44, 11,2049–2057.
I SHIMARU , A. 1978. Wave Propagation and Scattering in RandomMedia. Academic Press, New York, USA.
JAKOB , W., A RBREE , A., M OON , J. T., BALA , K., AND
M ARSCHNER , S. 2010. A radiative transfer framework for rendering materials with anisotropic structure. Tech. rep., Cornell University. (Expanded version), http://hdl.handle.net/1813/14982.
JAROSZ , W., C ARR , N. A., AND J ENSEN , H. W. 2009. ImportanceSampling Spherical Harmonics. Computer Graphics Forum (Proc.
Eurographics EG’09) 28, 2 (4), 577–586.
J ENSEN , H. W., AND C HRISTENSEN , P. H. 1998. Efficient simulation of light transport in scenes with participating media using
photon maps. In Proceedings of SIGGRAPH 98, 311–320.

L I , P., AND U REN , N. F. 1998. Analytical solution for the electric
potential due to a point source in an arbitrarily anisotropic halfspace. Journal of Engineering Mathematics 33, 129–140.
M ARSCHNER , S. R., W ESTIN , S. H., A RBREE , A., AND M OON ,J. T. 2005. Measuring and modeling the appearance of finishedwood. ACM Trans. Graph. 24 (July), 727–734.
M AX , N. L. 1994. Efficient light propagation for multiple
anisotropic volume scattering. In Fifth Eurographics Workshop
on Rendering, 87–104.
M ISHCHENKO , M. I., T RAVIS , L. D., AND L ACIS , A. A. 2006.
Multiple Scattering of Light by Particles. Cambridge U. Press.
M OULTON , J. 1990. Diffusion modeling of picosecond laser pulse
propagation in turbid media. Master’s thesis, McMaster University.
N EYRET, F. 1998. Modeling, animating, and rendering complex
scenes using volumetric textures. IEEE Transactions on Visualization and Computer Graphics 4, 1 (Jan./Mar.), 55–70.
PAULY, M., KOLLIG , T., AND K ELLER , A. 2000. Metropolis lighttransport for participating media. In Rendering Techniques 2000:
11th Eurographics Workshop on Rendering, 11–22.
P ERLIN , K., AND H OFFERT, E. M. 1989. Hypertexture. In
Computer Graphics (Proceedings of SIGGRAPH 89), 253–262.
P REISENDORFER , R. 1976. Hydrologic Optics. US Dept Commerce.
P REMOZE , S., A SHIKHMIN , M., T ESSENDORF, J., R AMAMOOR THI , R., AND NAYAR , S. 2004. Practical rendering of multiple
scattering effects in participating media. In Rendering Techniques
2004: 15th Eurographics Workshop on Rendering, 363–374.

J ENSEN , H. W., M ARSCHNER , S. R., L EVOY, M., AND H ANRA HAN , P. 2001. A practical model for subsurface light transport.
In Proceedings of SIGGRAPH 01, 511–518.

RUSHMEIER , H. E., AND T ORRANCE , K. E. 1987. The zonal
method for calculating light intensities in the presence of a participating medium. In Computer Graphics (Proceedings of SIGGRAPH 87), 293–302.

J OHNSON , P. M., AND L AGENDIJK , A. 2009. Optical anisotropic
diffusion: new model systems and theoretical modeling. Journalof Biomedical Optics 14, 5.

RYZHIK , L., PAPANICOLAOU , G., AND K ELLER , J. B. 1996.
Transport equations for elastic and other waves in random media.
Wave Motion 24, 1–44.

K AJIYA , J. T., AND H ERZEN , B. P. V. 1984. Ray tracing volume
densities. In Computer Graphics (Proceedings of SIGGRAPH84), 165–174.

S UN , B., R AMAMOORTHI , R., NARASIMHAN , S. G., AND NAYAR ,S. K. 2005. A practical analytic single scattering model for realtime rendering. ACM Trans. Graph. 24, 3 (Aug.), 1040–1049.

K AJIYA , J. T., AND K AY, T. L. 1989. Rendering fur with three
dimensional textures. In Computer Graphics (Proceedings of
SIGGRAPH 89), 271–280.
K ALDOR , J. M., JAMES , D. L., AND M ARSCHNER , S. 2008.
Simulating knitted cloth at the yarn level. In SIGGRAPH ’08:
ACM SIGGRAPH 2008 papers, 1–9.
K IENLE , A., F ORSTER , F., AND H IBST, R. 2004. Anisotropyof light propagation in biological tissue. Optics letters 29, 22,2617–2619.
L ACROUTE , P., AND L EVOY, M. 1994. Fast volume rendering
using a shear-warp factorization of the viewing transformation.
In Proceedings of SIGGRAPH 94, 451–458.
L AFORTUNE , E. P., AND W ILLEMS , Y. D. 1996. Rendering participating media with bidirectional path tracing. In Eurographics
Rendering Workshop 1996, 91–100.
L EVOY, M. 1988. Display of surfaces from volume data. IEEE
Computer Graphics & Applications 8, 3 (May), 29–37.

VAN DE H ULST, H. C. 1957. Light Scattering by Small Particles.

John Wiley & Sons.
WALTER , B., M ARSCHNER , S. R., L I , H., AND T ORRANCE , K. E.
2007. Microfacet models for refraction through rough surfaces.
In Eurographics Workshop on Rendering 2007, 195–206.
W ESTOVER , L. 1990. Footprint evaluation for volume rendering. In
Computer Graphics (Proceedings of SIGGRAPH 90), 367–376.
X U , Y.-Q., C HEN , Y., L IN , S., Z HONG , H., W U , E., G UO , B.,
AND S HUM , H.-Y. 2001. Photo-realistic rendering of knitwear
using the lumislice. In Proc. ACM SIGGRAPH 2001, 391–398.