To appear in the SIGGRAPH conference proceedings

A Practical Model for Subsurface Light Transport
Henrik Wann Jensen

Stephen R. Marschner

Marc Levoy

Pat Hanrahan

Stanford University
Abstract

1.1This paper introduces a simple model for subsurface light transport
in translucent materials. The model enables efficient simulation ofeffects that BRDF models cannot capture, such as color bleeding
within materials and diffusion of light across shadow boundaries.
The technique is efficient even for anisotropic, highly scattering
media that are expensive to simulate using existing methods. The
model combines an exact solution for single scattering with a dipole
point source diffusion approximation for multiple scattering. We
also have designed a new, rapid image-based measurement technique for determining the optical properties of translucent materials. We validate the model by comparing predicted and measured
values and show how the technique can be used to recover the optical properties of a variety of materials, including milk, marble, andskin. Finally, we describe sampling techniques that allow the model
to be used within a conventional ray tracer.

Almost all BRDF models are derived exclusively from surface scattering, with any subsurface scattering approximated by a Lambertian component. An exception is the model by Hanrahan andKrueger [10] which includes an analytic expression for single scattering in a homogeneous, uniformly lit slab. However, all BRDF
models ultimately assume that light scatters at one surface point andthey do not model subsurface transport from one point to another.
Subsurface transport can be simulated accurately but slowly bysolving the full radiative transfer equation [1]. Only a few papers ingraphics have taken this approach to subsurface scattering. Dorsey
et al. [5] simulated full subsurface scattering using photon mappingto capture the appearance of weathering in stone. Pharr and Hanrahan [15] used scattering functions to simulate subsurface scattering. These approaches, while capable of simulating all of theeffects of subsurface scattering, are computationally very expensive compared to the simulation of opaque materials. Techniques
based on path sampling are particularly inefficient for highly scattering materials, such as milk and skin, in which light scatters multiple (often several hundred) times before exiting the material. For
highly scattering media Stam [17] introduced the use of diffusion
theory. He solved a diffusion equation approximation using a multigrid method, and used this method to render clouds with multiple
scattering.
Subsurface scattering is also important in medical physics, wheremodels have been developed to describe the scattering of laser lightin human tissue [6, 8]. In that context, diffusion theory is often used
to predict as well as to measure the optical properties of highly scattering materials. We have extended this theory for use in computer
graphics by adding exact single scattering, support for arbitrary geometry, and a practical sampling technique for rendering.
In measurements of appearance for computer graphics, subsurface scattering has rarely been considered. Debevec et al. [3] measured light reflection from human faces, which included contributions from subsurface scattering, but they did not relate the data to
the physical properties of the material. Again building on medical
physics research [8, 9], we have extended a methodology developed
for measuring biological tissues into a rapid image-based appearance measurement technique for translucent materials. This method
examines the radial reflectance profile resulting from a beam illuminating the sample material. By fitting an expression derived fromdiffusion theory it is possible to estimate the absorption and scattering properties of the material.

Keywords: Subsurface scattering, BSSRDF, reflection models,
light transport, diffusion theory, realistic image synthesis

1

Introduction

Accurately modeling the scattering of light by materials is fundamental for realistic image synthesis. Even the most sophisticatedlight transport algorithms fail to produce convincing results if thelocal scattering models are too simple. Therefore a great deal ofresearch has gone into describing the scattering of light from materials.
Previous research has focused on developing models for thebidirectional reflectance distribution function (BRDF). The BRDF
was introduced by Nicodemus [14] as a simplification of themore general bidirectional surface scattering distribution function(BSSRDF). The BSSRDF can describe light transport between any
two rays that hit a surface, whereas the BRDF assumes that light entering a material leaves the material at the same position (Figure 1).
This approximation is valid for metals, but it fails for translucent
materials, which exhibit significant transport below the surface.
Even for many materials that do not seem very translucent, using
the BRDF creates a hard, distinctly computer-generated appearance
because it does not locally blend surface features such as color andgeometry. Only methods that consider subsurface scattering cancapture the true appearance of translucent materials, such as marble, cloth, paper, skin, milk, cheese, bread, meat, fruits, plants, fish,
ocean water, snow, etc.

2

Previous Work

TheoryThe BSSRDF, S, relates the outgoing radiance, Lo (xo , ω~ o ) at the
point xo in direction ω~ o , to the incident flux, Φi (xi , ω~ i ) at the pointxi from direction ω~ i [14]:dLo (xo , ω~ o ) = S(xi , ω~ i ; xo , ω~ o ) dΦi (xi , ω~ i ).
The BRDF is an approximation of the BSSRDF for which it
is assumed that light enters and leaves at the same point (i.e.,
xo = xi ). Given a BSSRDF, the outgoing radiance is computed

To appear in the SIGGRAPH conference proceedings

(a)

S
Rd
Fr
Ft
Fdr
~
E
φ
σa
σs
σt
σt0
σtr
D
α
p
η
g
Q
Q0
~1
Q

(b)

Figure 1: Scattering of light in (a) a BRDF, and (b) a BSSRDF.
by integrating the incident radiance over incoming directions andarea, A:Z Z

S(xi , ω~ i ; xo , ω~ o ) Li (xi , ω~ i ) (~n·~ωi ) dωi dA(xi ).

Lo (xo , ω~ o) =

A

2π

Light propagation in a participating medium is described by theradiative transport equation, often referred to in computer graphics
as the volume rendering equation:

BSSRDF
Diffuse BSSRDF
Fresnel reflectanceFresnel transmittance
Diffuse Fresnel reflectance
Vector irradiance
Radiant fluence
Absorption coefficient
Scattering coefficient
Extinction coefficient
Reduced extinction coefficient
Effective extinction coefficient
Diffusion constant
AlbedoPhase function
Relative index of refraction
Mean cosine of the scattering angle
Volume source distribution
0th-order source distribution
1st-order source distribution
Figure 2: Selected symbols.

Z~(~ω ·∇)L(x,ω~ ) = −σt L(x, ω~ )+σs p(~
ω, ω~ 0 )L(x, ω~ 0 ) dω 0 +Q(x, ω
~ ).
4π

2.1

In this equation, the properties of the medium are described by
the absorption coefficient σa , the scattering coefficient σs , and thephase function p(~
ω, ω~ 0 ). The extinction coefficient σt is defined
as,
σ
=
σ
+
σ
.
We
assume the phase function is normalized,
a
s
R t
0
0
p(~
ω
,
ω
~
)
dω
=
1
and
is a function only of the phase angle,
4π
p(~
ω, ω
~ 0 ) = p(~
ω·ω~ 0 ). The mean cosine, g, of the scattering angle
isZ

The diffusion approximation is based on the observation thatthe light distribution in highly scattering media tends to become
isotropic. This is true even if the initial light source distribution
and the phase function are highly anisotropic. Each scattering event
blurs the light distribution, and as a result the light distribution tends
toward uniformity as the number of scattering events increases.
In this situation, the radiance may be approximated by a twoterm expansion involving the radiant fluence and the vector irradiance:
1
3~
L(x, ω
~) =φ(x) +
ω~ · E(x).
4π4πThe constants are determined by the definitions of fluence and vector irradiance.
The diffusion equation follows from this approximation. Specifically, we substitute this two-term expansion of the radiance into
the radiative transport equation and then integrate over ω~ ; for the
algebraic details consult Ishimaru [12]. The result is(~ω·ω~ 0 )p(~ω·ω~ 0 ) dω 0 .

g=

4π

If g is positive, the phase function is predominantly forward scattering; if g is negative, backward scattering dominates. A constantphase function results in isotropic scattering (g = 0).
For an infinitesimal beam entering a homogeneous medium,the incoming radiance will decrease exponentially with distance s.
This is referred to as the reduced intensity:
Lri (xi + s~
ωi , ω
~ i ) = e−σt s Li (xi , ω~ i ).

~~~ 1 (x).
∇φ(x)= −3σt0 E(x)+Q

The first-order scattering of the reduced intensity, Lri , may betreated as a volumetric source:
Q(x, ω~ ) = σsZ

given byσt0 = σs0 + σa where σs0 = σs (1 − g) .

4π

The reduced scattering coefficient σs0 scales the original scattering
coefficient by a factor of (1 − g). Intuitively, once light becomesisotropic, only backward scattering terms change the net flux; forward scattering is indistinguishable from no scattering.
In the case where there are no sources, or where the sources are~ 1 vanishes from Equation 2. Then the vector irradianceisotropic, Q
is the gradient of the scalar fluence,

To gain insight into the volumetric behavior of light propagation, it is useful to integrate the radiative transport equation over all
directions ω~ at a point x which yields(1)This equation R relates the scalar irradiance, or fluence,
φ(x)=L(x, ω~ ) dω, and the vector irradiance,
4π
R
~
E(x)
=
L(x, ω~ )~ω dω. In the absence of loss due to ab-

~
~
E(x)
= −D∇φ(x).

4πsorption or gain due to a volumetric light source (Q0 = 0), thedivergence of the vector irradiance equals zero. In this equation,
we introduce a 0th-order source term, Q0 , and later we will need
~ 1 , where
the 1st-order source term, Q
Q0 (x) =Z
4πQ(x, ω~ ) dω,~ 1 (x) =
Q

Z

(2)

Here we have used the reduced extinction coefficient, σt0 , which isp(~ω0 , ω~ )Lri (x, ω~ 0 ) dω 0 .

~ · E(x)~∇= −σa φ(x) + Q0 (x).

The Diffusion Approximation

Here D = 3σ1 0 is the diffusion constant. This equation makes pretcise the intuitive notion that there is net energy flow (i.e., non-zero
vector irradiance) from regions of high energy density (high fluence) to regions of low energy density.
Finally, substituting Equation 2 into Equation 1, we arrive at theclassic diffusion equation

Q(x, ω
~ )~
ω dω.

4π

~ ·Q
~ 1 (x).
D∇2 φ(x) = σa φ(x) − Q0 (x) + 3D∇2

To appear in the SIGGRAPH conference proceedings
Finally, the boundary condition allows us to compute the diffuse
BSSRDF, Rd . Rd is equal to the radiant exitance divided by the~ s ))
incident flux. The radiant exitance leaving the surface (~n · E(xis equal to the gradient of the fluence at the surface

The diffusion equation has a simple solution in the case of a single isotropic point light source in an infinite medium.
φ(x) =

Φ e−σtr r(x),4π D r(x)Rd (r) = −D

where Φ is the power of the point light source,
p r is the distance tothe location of the point source, and σtr = 3σa σt0 is the effective
transport coefficient. The point source results in an energy density
in the volume with an exponential falloff.
In the case of a scattering medium in a finite region of space, thediffusion equation must be solved subject to the appropriate boundary conditions. The boundary condition is that the net inward diffuse flux is zero at each point, xs , on the surface

Z

where r = ||xs − xi ||.
In the case of finite media, the diffusion equation does not ingeneral have an analytical solution. In this paper we are interested
in subsurface reflection, which is often modeled as a semi-infinite
plane-parallel medium. Several authors have analyzed the planeparallel problem for simple source geometries, in particular, approximations of a cylindrical beam entering the media. Exact formulas exist, but they involve an infinite sum of Bessel functions[9, 16]. We seek a simple formula suitable for modeling subsurface
reflection that does not involve infinite sums or numerical solutionof a partial differential equation.
Eason [6] and Farrell et al. [8] have developed a method forapproximating the volumetric source distribution using two point
sources; that is, a dipole. Eason introduced this idea and derived
explicit formulae for the dipoles for various source geometries,
such as a cylindrical beam, by expanding the source distributionsin terms of their moments. Farrell et al. proposed using a single
dipole to represent the incident source distribution. They found asingle dipole to be as accurate as, or, in some cases, more accurate
than using the diffusion approximation with the true source distribution.
The dipole method consists of positioning two point sources near
the surface in such a way as to satisfy the required boundary condition [6] (see Figure 3). One point source, the positive real lightsource, is located at the distance zr beneath the surface, and theother, the negative virtual light source, is located above the surface
at a distance zv = zr + 4AD. The resulting fluence isL(xs , ω~ )(~ω · ~n(xs )) dω = 0.

2π −Here, 2π − denotes integration over the hemisphere of inward directions. Using the two-term expansion, the boundary condition is~φ(xs ) − 2D(~n · ∇)φ(xs ) = 0.

(3)

The minus sign in the second term results from the convention thatthe surface normal points outward, whereas the integral is over inward directions.
Equation 3 covers the case where the two layers have matchingindices of refraction, but another important case is where these indices differ. When an interface exists between media with different
refractive indices, there is a reflection at the interface. Assuming Fr
is the Fresnel formula for the reflectance at a dielectric interface, theaverage diffuse Fresnel reflectance is
Fdr =

Z

~
(~n · ∇φ)(xs),
dΦi (xi )Fr (η, ~n · ω~ 0 )(~n · ω~ 0 ) dω 0 ,2π

where η is the relative index of refraction of the medium with thereflected ray to the other medium. Fdr may be computed analytically from the Fresnel formula [13]. However, we will use a rational
approximation of the measured diffuse reflectance [7]:Fdr = −

φ(x) =

0.710
1.440
+
+ 0.668 + 0.0636η.
η2
η

L(x, ω~ )(~ω · ~n− )dω = Fdr

2π −

Z

L(x, ω~ )(~ω · ~n+ ) dω.

2π +Rd (r) = −D

Here the + and − subscript means outward and inward directions
respectively. This yields=~~
φ(xs ) − 2D(~n · ∇)φ(xn · ∇)φ(xs ) = Fdr φ(xs ) + 2D(~s) .

 −σ d
e

tr r

dr

−

e−σtr dv
dv

,

where dr = ||x − xr || is the distance from x to the real source, anddv = ||x − xv || is the distance from x to the virtual source. Farrell et al. [8] proposed positioning the real light source at distancezr = 1/σt0 , or one mean free path, below the surface. They only
considered light parallel to the normal. For other light directions
reciprocity can be enforced by still placing the light source 1/σt0
straight below xi .
The diffuse reflectance due to the dipole source can now be computed.

The resulting boundary condition between two media with different
indices of refraction is

Z

Φ
4π D



~
(~n · ∇φ(x
s ))
dΦi



e−σtr dv
α0
e−σtr dr+ zv (σtr dv + 1)
(σtr dr + 1)
4π
σt0 d3r
σt0 d3v

.

(4)Note that the difference in signs between the two sides of this equation occurs because one integral is over outward directions and theother is over inward directions. Rearranging terms,

Lastly, we need to take into account the Fresnel reflection at theboundary for both the incoming light and the outgoing radiance.

~φ(xs ) − 2AD(~n · ∇)φ(xs ) = 0.

Sd (xi , ω~ i ; xo , ω~ o) =

This boundary condition is the same as when the indices of refraction match (Equation 3); the only difference is that 2D is replaced
by 2AD, where
1 + FdrA=.
1 − Fdr1
Ft (η, ω~ i )Rd (||xi − xo ||)Ft (η, ω~ o ) (5)πwhere Sd is the diffusion term of the BSSRDF. This term represents
multiple scattering (one scattering event is already included in theconversion to a point source). The next section explains how tocompute the contribution due to single scattering.

3

To appear in the SIGGRAPH conference proceedings

zv


xi

r

xo


zr

Figure 3: An incoming ray is transformed into a dipole source forthe diffusion approximation.

2.2

Figure 4: Single scattering occurs only when the refracted incoming
and outgoing rays intersect, and is computed as an integral over path
length s along the refracted outgoing ray.

Single Scattering Term

Hanrahan and Krueger [10] have derived a BRDF model for subsurface reflection that analytically computes the total first-order scattering from a flat, uniformly lit, homogeneous slab. In this section,
we show how their BRDF can be extended to a BSSRDF in orderto account for local variations in lighting over the surface.
(1)
The total outgoing radiance, Lo , due to single scattering iscomputed by integrating the incident radiance along the refracted
outgoing ray (see Figure 4):(1)

Lo (xo , ω~ o ) = σs (xo )

Z

0

0

F p(~
ωi · ω
~ o)

2π

=

Z Z
A

S

(1)

Z ∞

e

−σtc sLi (xi , ω~ i ) ds d~ωi

2.4

BRDF Approximation

We can approximate the BSSRDF with a BRDF by assuming thatthe incident illumination is uniform. This assumption makes it possible to integrate the BSSRDF over the surface. By integrating the
diffusion term we find the total diffuse reflectance Rd of the material as:
Z ∞

 √
√
4
0
0
α0
Rd = 2π
Rd (r) r dr =1 + e− 3 A 3(1−α ) e− 3(1−α ) .
2
0

(6)

0

Notice how the diffuse reflectance only depends on the reduced
albedo and the internal reflection parameter A.
The integration of the single scattering term results in the model
presented in [10]. For a semi-infinite medium this gives:(xi , ω~ i ; xo , ω~ o ) Li (xi , ω~ i ) (~n·ω~ i ) dωi dA(xi ).

2πHere F = Ft (η, ω~ o )Ft (η, ω~ i ) is the product of the two Fresnel
transmission terms, and ω~ i0 and ω~ o0 are the refracted incoming andoutgoing directions. The combined extinction coefficient σtc is
given by σtc = σt (xo ) + Gσt (xi ), where G is a geometry fac|~
n ·~
ω0 |
tor; for a flat surface G = |~nii ·~ωo0 | . The single scattering BSSRDF,

fr(1) (x, ω~ i, ω~ o ) = αFp(~ωi0 · ω~ o0 ).
0|~n · ω~ i | + |~n · ω~ o0 |

The complete BRDF model is the sum of the diffuse reflectancescaled by the Fresnel term and the single scattering approximation:i

S (1) , is defined implicitly by the second line of this equation. Notethat there is a change of variables between the first line, which integrates only over the configurations where the two refracted rays
intersect, and the second line, which integrates over all incoming
and outgoing rays. This implies that the distribution S (1) contains
a delta function.

2.3s

fr (x, ω~ i, ω~ o ) = fr(1) (x, ω~ i, ω~ o) + F

Rd.
πThis model has the same parameters as the BSSRDF. It is similar
to the BRDF model presented in [10], but with the important difference that the amount of diffusely reflected light is computed fromthe intrinsic material parameters. The BRDF approximation is useful for opaque materials, which have a very short mean free path.

The BSSRDF Model

The complete BSSRDF model is a sum of the diffusion approximation and the single scattering term:

3

S(xi , ω~ i ; xo , ω~ o ) = Sd (xi , ω~ i ; xo , ω~ o ) + S (1) (xi , ω~ i ; xo , ω~ o)Measuring the BSSRDF

To verify our BSSRDF model, and to determine appropriate parameters for rendering different kinds of materials, we used the
diffusion theory of Section 2 to make measurements of subsurface
scattering in several media. Our measurement approach applies totranslucent materials for which σa  σs , implying that far enoughaway from the point of illumination, we may neglect single scattering and use the diffusion term to relate measurements to material
parameters.
When multiple scattering dominates, Equation 4 predicts the radiant exitance per unit incident flux that will be observed due to anarrow incident beam, as a function of distance from the point of
illumination. To make the corresponding measurement, we illuminate the surface of a sample with a tightly focused beam of whitelight and take a photograph using a 3-CCD video camera to observe
the radiant exitance across the entire surface. We keep our observations at constant angles so that the Fresnel term remains constantfor all the measurements. Figure 5(a) illustrates our measurement
setup.

Here Sd is evaluated using Equation 5 and S (1) is evaluated using Equation 6. The parameters for the BSSRDF are: σa , σs0 , η,and possibly a phase function (without a phase function the scattering can be modeled as isotropic). This model accounts for lighttransport between different locations on the surface, and it simulates both the directional component (due to single scattering) as
well as the diffuse component (due to multiple scattering).
Finally, note the distances involved in both the single scatteringterm and the diffusion approximations. The average exit point isapproximately one mean free path from the entry point. However,these two mean free paths have quite different length scales. In
the single scattering case, the mean free path equals 1/σt ; in the
diffusion case, the mean free path equals 1/σtr . For translucent
materials where σa  σs0 and consequently σtr  σt , the single
scattering term decreases much faster than the diffusion term as the
distance to xo increases.

4

To appear in the SIGGRAPH conference proceedings
SourceMaterial

90cm

A

Camera


4˚

f /8
50mm

50˚

SampleApple
Chicken1
Chicken2
Cream
Ketchup
MarblePotatoSkimmilk
Skin1
Skin2
Spectralon
Wholemilk

σs0 [mm−1 ]RG
B
2.29 2.39 1.97
0.15 0.21 0.38
0.19 0.25 0.32
7.38 5.47 3.15
0.18 0.07 0.03
2.19 2.62 3.00
0.68 0.70 0.55
0.70 1.22 1.90
0.74 0.88 1.01
1.09 1.59 1.79
11.6 20.4 14.9
2.55 3.21 3.77

R
0.00300.015
0.018
0.00020.061
0.00210.00240.0014
0.032
0.013
0.00
0.0011(a)

σa [mm−1 ]GB
0.0034 0.046
0.077
0.19
0.088
0.20
0.0028 0.0163
0.97
1.45
0.0041 0.0071
0.0090 0.12
0.0025 0.0142
0.17
0.48
0.070
0.145
0.00
0.00
0.0024 0.014

Diffuse ReflectanceR
G
B
0.85 0.84 0.53
0.31 0.15 0.10
0.32 0.16 0.10
0.98 0.90 0.73
0.16 0.01 0.00
0.83 0.79 0.75
0.77 0.62 0.21
0.81 0.81 0.69
0.44 0.22 0.13
0.63 0.44 0.34
1.00 1.00 1.00
0.91 0.88 0.76

η
1.3
1.3
1.3
1.3
1.3
1.5
1.3
1.3
1.3
1.3
1.3
1.3

(b)
Figure 5: (a) Measurement apparatus, (b) measured parameters for several materials.

100

Because the signal falls off exponentially away from the point of
illumination, the measurement must span a wide dynamic range. To
this end we used a series of different exposure times, ranging from 1
millisecond to 4 seconds, and assembled a high-dynamic-range image using a modified version of Debevec and Malik’s technique [4].
To reduce the effects of stray light and fixed-pattern CCD noise, we
subtracted a dark image, taken with the illumination beam blocked
just before the focusing lens (point A in Figure 5(a)), from eachmeasurement and reference image. The resulting images had a dynamic range of around 105 (the small amount of total energy in theimage reduces the effects of lens and camera flare, allowing higher
dynamic range than might otherwise be possible).
To interpret the measurements, we examined only a 1D sliceof each measurement image, corresponding to a line on the surface through the illumination point and perpendicular to the camera’s view direction. Under the assumption that light exits diffusely1 , the pixel values pi in this slice (see Figure 6 for an example) are measurements of radiant exitance as a function of distance on the surface. Since Rd gives the ratio of this quantity to Φ,
pi = KΦRd (ri ), where K is an unknown constant. To eliminate
the scale factor, we also took a reference image with the sample
replaced by a white ideal diffuse reflector (Labsphere Spectralon,
reflectance > 0.99). By summing all the pixels in this image, we
can integrate the radiant exitance to get the total flux exiting the surface, which for this special material is equal to the incident flux Φ.
With the same constant K as above, this sum is KΦ/A, where A
is the (known) area on the sample’s surface subtended by one pixel.
The measured value for Rd (ri ) can then be computed as pi /(KΦ).
In principle, σa and σs0 can be determined by fitting the relative reflectance curve with Equation 4 over a range of distances
far enough from the illumination point to allow the use of diffusion theory [8]. However, we found this fitting problem to be illconditioned enough that the uncertainty in the resulting parameters
led to too much uncertainty in the appearance of the material, especially the total diffuse reflectance.
We remove this ill-conditioning by measuring the total diffuse
reflectance R (which is the sum of the measurement image divided
by the sum of the reference Rimage) and computing the least-squares
fit subject to the constraint Rd dA = R.
Figure 6 shows how these measurements confirm the diffusion
theory for a sample of white marble (only the camera’s green chan-







10-1

data not used in fit
data used in fit
diffusion theory using fitted parameters
Monte Carlo simulation using fitted parameters

Rd (mm –2 )

10-2

10-3



10-4

10-5

10-6
0

2

4



6



8


10
r (mm)

12

14

16

18

20

Figure 6: Measurements for marble (green wavelength band) plotted with fit to diffusion theory and confirming Monte Carlo simulation.
nel is shown). Fitting the theory (solid line) to the data (points) led
to the parameters σa = 0.0041/mm, σs0 = 2.6/mm. The reflectancecomputed by a Monte Carlo simulation using these values (dashedline) confirms the correctness of the computed parameters. Fitted
values for several other materials appear in the table in Figure 5(b).
Note, that we used empirical values for the index of refraction for
most of the materials. Also note that the diffusion theory is assuming that σs  σa , and as such the parameters for the relatively
opaque materials (such as the blue wavelength in ketchup) may beless accurate.

4

Rendering Using the BSSRDF

The BSSRDF model derived in the theory section only applies tosemi-infinite homogeneous media. A similar derivation is not possible in the presence of arbitrary geometry and texture variation.
However, we can use some of the intuition behind the theory to extend it to a practical model for computer graphics. Specifically, weneed to consider:

1 We verified this assumption for marble by examining the reflectance fordifferent outgoing angles, and it closely resembled a Lambertian material
scaled by a Fresnel transmission term.

5

To appear in the SIGGRAPH conference proceedings

(a)

(b)

(a)

(b)

Figure 8: Scattering of laser light in a marble block. The marble
block is 40mm. wide and has a significant amount of subsurface
scattering. The picture on the left is a photograph of the marble
block, and the picture on the right is a synthetic rendering of a similarly sized cube using the BSSRDF model and the measured scattering properties of the marble. Note how the appearance of the two
images is very similar.

Figure 7: (a) Sampling a BRDF (traditional sampling), (b) samplinga BSSRDF (the sample points are distributed both over the surface
as well as the light).
• Efficient integration of the BSSRDF including importancesampling• Single scattering evaluation for arbitrary geometry

Diffusion approximation for arbitrary geometry: An important component of the diffusion approximation is the use of thedipole source. If the geometry is locally flat we can get a very
good approximation by using a similar dipole source configuration
as that for flat materials (i.e., we always place the light source 1/σt0
straight below xi ). Special care must be taken in the presence ofhighly curved surfaces; we handle this case by always evaluating
the diffusion term with a minimum distance of 1/σt0 . In this way
we eliminate singularities at sharp edges where the source can beplaced arbitrarily close to xo . We found this approach to work very
well in our experiments.
Texture: We approximate textured materials by making a few
small changes to the usage of the BSSRDF. We only consider texture variation at the surface — effects due to volumetric texture
variation would require a full participating media simulation. For
the diffusion approximation we always use the material parameters
at xi , which ensures a natural local blending of the texture properties. For the single scattering term we use σs (xo ) and σt (xo ) along
the refracted outgoing ray, and σt (xi ) along the refracted incident
ray. This variation is included in Equation 6.

• Diffusion approximation for arbitrary geometry• Texture (spatial variation on the object surface).
In this section we explain how to do this in a ray-tracing context.
Integrating the BSSRDF: At each ray-object intersection traditional lighting models (based on BRDFs) need just a point anda normal to compute the outgoing radiance (Figure 7(a)). For the
BSSRDF it is necessary to integrate the incoming lighting over an
area of the surface (Figure 7(b)). We do this by stochastically sampling the location of both endpoints of the shadow ray — this canbe seen as an extension of the classical distribution ray tracing technique for sampling area light sources [2]. To efficiently sample
locations on the surface we exploit the exponential falloff in the
diffusion term and the single scattering term. We sample the twoterms of the BSSRDF separately, since the single scattering sample locations must be along the refracted outgoing ray whereas the
diffusion samples should be distributed around xo .
More specifically, for the diffusion term, we use standard Monte
Carlo techniques to randomly sample the surface with density(σtr e−σtr d ) at some distance d from xo .
Single scattering is reparameterized since the incoming ray andthe outgoing ray must intersect. Our technique is explained in thefollowing section.
Single scattering evaluation for arbitrary geometry: Single scattering is evaluated using Monte Carlo integration along
the refracted outgoing ray. We pick a random distance, s0o =log(ξ)/σt (xo ), along the refracted outgoing ray. Here ξ ∈ ]0, 1]
is a uniformly distributed random number. For this sample location
we compute the outscattered radiance as:
L(1)~ o )=o (xo , ω

5

We have implemented the BSSRDF model in a Monte Carlo ray
tracer, and in this section we will present a number of experimental results obtained with this implementation. All simulations havebeen done on a dual 800MHz Pentium III PC running Linux and theimages have been rendered with 4 samples per pixel and a width of1024 pixels.
Our first simulation is shown in Figure 8, which compares aside photograph of a marble cube illuminated from above with asynthetic rendering. The synthetic image is rendered using the
BSSRDF model and the measured parameters for marble (from thetable in Figure 5). We only used a simple cube to approximate therounded marble block, so there are natural visible differences along
the edges. Nonetheless, the BSSRDF model faithfully renders theappearance including the scattered light exiting from the side of themarble cube.
Figure 9 shows several different simulations of subsurface scattering in a marble bust (1.3 million triangles) illuminated from behind. The BSSRDF simulation mostly matches the appearance ofthe full Monte Carlo simulation, yet is significantly faster (5 minutes vs. 1250 minutes). The hair at the back of the head is slightly
darker in the BSSRDF simulation; we believe this is due to theforced 1/σt0 distance in the diffusion approximation. A similar rendering was done using photon mapping in [5] in roughly 12 minutes (scaled to the speed of our computer). However, the photon
mapping method requires a full 3D-description of the material, it
requires memory to store the photons, and it becomes costly forσs (xo )F p(~
ωi · ω~ o ) −s0i σt (xi ) −s0o σt (xo )eeLi (xi , ω~ i ).
σtc

Here s0i is the distance that the sample ray moves through the material. Optimizing this equation to sample direct illumination (withshadow rays) is difficult for arbitrary geometry since it requires
finding the point at the surface where the shadow ray is refracted.
However, in practice a good approximation can be found by usinga shadow ray that does not refract at the surface — this assumes
that the light source is far away compared to the mean free path ofthe medium. We can use Snell’s law to estimate the true refracted
distance through the medium of the incoming ray:s0i = si q
1−|~
ωi · ~ni |


1 2

η

Results.

(1 − |~
ωi · ~n(xi )|2 )Here si is the observed distance and s0i is the refracted distance.

6

To appear in the SIGGRAPH conference proceedings

(a)

(b)

(d)

(c)

(e)

(f)

Figure 9: A simulation of subsurface scattering in a marble bust. The marble bust is illuminated from behind and rendered using: (a) the
BRDF approximation (in 2 minutes), (b) the BSSRDF approximation (in 5 minutes), and (c) a full Monte Carlo simulation (in 1250 minutes).
Notice how the BSSRDF model matches the appearance of the Monte Carlo simulation, yet is significantly faster. The images in (d–f) show
the different components of the BSSRDF: (d) single scattering term, (e) diffusion term, and (f) Fresnel term.

highly scattering materials (such as milk and skin).
A particularly interesting aspect of the BSSRDF simulation isthat it is able to capture the smooth appearance of the marble surface. In comparison the BRDF simulation gives a very hard appearance where even tiny bumps on the surface are visible (this isa classic problem in realistic image synthesis where objects often
look hard and unreal).
For the marble we used synthetic scattering and absorption coefficients, since we wanted to test the difficult case when the average scattering albedo is 0.5 (here the contribution from diffusionand single scattering is approximately the same). Figure 9 demonstrates how the sum of both single scattering and the diffusion term
is necessary to match the Monte Carlo simulation.
Figure 10 contains three renderings of milk. The first rendering uses a diffuse reflection model; the others use the BSSRDF
model and our measurements for skim milk and whole milk. Notice
how the diffuse milk looks unreal and too opaque compared to the
BSSRDF images, even though multiple scattering dominates andthe radiant exitance due to subsurface scattering is very diffuse. It
is interesting that the BSSRDF simulations are capable of capturing
the subtle details in the appearance of milk, making the milk look
more bluish at the front and more reddish at the back. This is due
to Rayleigh scattering that causes shorter wavelengths of light to be
scattered more than longer wavelengths.
Skin is a material that is particularly difficult to render using
methods that simulate subsurface scattering by sampling ray paths
through the material. This is due to the fact that skin is highly
scattering (typical albedo is 0.95) and also very anisotropic (typical average cosine of the scattering angle is 0.85). Both of theseproperties mean that the average number of scattering events of aphoton is very high (often more than 100). In addition skin is verytranslucent, and it cannot be rendered correctly using a BRDF (seeFigure 11). A complete skin model requires multiple layers, but areasonable approximation can be obtained using just one layer. In
Figure 11 we have rendered a simple face model using the BSSRDF
and our measured values for skin (skin1). Here we also used theHenyey-Greenstein phase function [11] with g = 0.85 as the estimated mean cosine of the scattering angle. The skin measurements
are from an arm (which is likely more translucent than skin on theface), but the overall appearance is still realistic considering the lackof spatial variation (texture). The BSSRDF gives the skin a soft appearance, and it renders the color bleeding in the shadow region
below the nose. Here, the absorption by blood is particularly noticeable as the light that scatters deep in the skin is redder. For thissimulation the diffusion term is much larger than the single scattering term. This means that skin reflects light fairly diffusely, butalso that internal color bleeding is an important factor. The BRDF
image was rendered in 7 minutes, the BSSRDF image was renderedin 17 minutes.

6

Conclusion and Future WorkIn this paper we have presented a new practical BSSRDF model
for computer graphics. The model combines a dipole diffusion approximation with an accurate single scattering computation. We
have shown how the model can be used to measure the scattering
properties of translucent materials, and how the measured values
can be used to reproduce the results of the measurements as well
as synthetic renderings. We evaluate the BSSRDF by sampling theincoming light over the surface, and we demonstrate how this technique is capable of capturing the soft and smooth appearance oftranslucent materials.
In the future we plan to extend the model to multiple layers aswell as include support for efficient global illumination.

7

To appear in the SIGGRAPH conference proceedings

(a)

(b)

(c)

Figure 10: A glass of milk: (a) diffuse (BRDF), (b) skim (BSSRDF) and (c) whole (BSSRDF). (b) and (c) are using our measured values.
The rendering times are 2 minutes for (a), and 4 minutes for (b) and (c); this includes caustics and global illumination on the marble table anda depth-of-field simulation.

7

Acknowledgements

Special thanks to Steven Stahlberg for modeling the face. Thanks
to the SIGGRAPH reviewers and to Maryann Simmons and Heidi
Marschner for helpful comments on the manuscript. This research
was funded in part by the National Science Foundation Information
Technology Research grant (IIS-0085864). The first author was
also supported by DARPA (DABT63-95-C-0085), and the second
author was also supported by Honda North America, Inc.

References[1] S. Chandrasekhar. Radiative Transfer. Oxford Univ. Press, 1960.
[2] R. L. Cook, T. Porter, and L. Carpenter. Distributed ray tracing. In ACM Computer Graphics (SIGGRAPH’84), volume 18, pages 137–145, July 1984.
[3] P. Debevec, T. Hawkins, C. Tchou, H. Duiker, W. Sarokin, and M. Sagar. Acquiring the reflectance field of a human face. In Computer Graphics Proceedings,
Annual Conference Series, 2000, pages 145–156, July 2000.
[4] P. E. Debevec and J. Malik. Recovering high dynamic range radiance maps fromphotographs. In Computer Graphics Proceedings, Annual Conference Series,
1997, pages 369–378, August 1997.

BRDF[5] J. Dorsey, A. Edelman, H. W. Jensen, J. Legakis, and H. K. Pedersen. Modelingand rendering of weathered stone. In Computer Graphics Proceedings, Annual
Conference Series, 1999, pages 225–234, August 1999.
[6] G. Eason, A. Veitch, R. Nisbet, and F. Turnbull. The theory of the backscatteringof light by blood. J. Physics, 11:1463–1479, 1978.
[7] W. G. Egan and T. W. Hilgeman. Optical Properties of Inhomogeneous Materials. Academic Press, New York, 1979.
[8] T. J. Farell, M. S. Patterson, and B. Wilson. A diffusion theory model of spatiallyresolved, steady-state diffuse reflectance for the noninvasive determination oftissue optical properties in vivo. Med. Phys., 19:879–888, 1992.
[9] R. A. Groenhuis, H. A. Ferwerda, and J. J. Ten Bosch. Scattering and absorptionof turbid materials determined from reflection measurements. 1: Theory. Applied
Optics, 22:2456–2462, 1983.
[10] P. Hanrahan and W. Krueger. Reflection from layered surfaces due to subsurface scattering. In ACM Computer Graphics (SIGGRAPH’93), pages 165–174,
August 1993.
[11] L. G. Henyey and J. L. Greenstein. Diffuse radiation in the galaxy. Astrophysics
Journal, 93:70–83, 1941.
[12] A. Ishimaru. Wave Propagation and Scattering in Random Media, volume 1.
Academic Press, New York, 1978.

BSSRDF
Figure 11: A face rendered using the BRDF model (top) and the
BSSRDF model (bottom). We used our measured values for skin
(skin1) and the same lighting conditions in both images (the BRDF
image also includes global illumination). The face geometry has
been modeled by hand; the lip-bumpmap is handpainted, and the
bumpmap on the skin is based on a gray-scale macro photographof a piece of skin. Even with global illumination the BRDF gives a
hard appearance. Compare this to the faithful soft appearance of theskin in the BSSRDF simulation. In addition the BSSRDF captures
the internal color bleeding in the shadow region under the nose.

[13] G. Kortum. Reflectance Spectroscopy. Springer-Verlag, 1969.
[14] F. E. Nicodemus, J. C. Richmond, J. J. Hsia, I. W. Ginsberg, and T. Limperis.
Geometric considerations and nomenclature for reflectance. Monograph 161,
National Bureau of Standards (US), October 1977.
[15] M. Pharr and P. Hanrahan. Monte Carlo evaluation of non-linear scattering equations for subsurface reflection. In Computer Graphics Proceedings, Annual Conference Series, 2000, pages 75–84, July 2000.
[16] L. Reynolds, C. Johnson, and A. Ishimaru. Applied Optics, 15:2059, 1976.
[17] J. Stam. Multiple scattering as a diffusion process. In Eurographics Rendering
Workshop 1995. Eurographics, June 1995.

8