Subsurface Scattering and the Diffusion Approximation

A skin surface, a marble slab, and a wax candle do not reflect light at the point where it arrives. Light enters at one place, scatters a few times over a distance comparable to a millimeter or a centimeter, and leaves somewhere else. A BRDF cannot express that, because a BRDF is a function of two directions at a single point. The bidirectional scattering-surface reflectance distribution function, or BSSRDF, is the generalization: a function of both points and both directions.

Writing the BSSRDF down is trivial. Evaluating it is not, and the machinery that makes it tractable is the diffusion approximation, which converts the transport problem into a Poisson-type equation with a closed-form solution in simple geometries.

What the extra degrees of freedom cost

The BSSRDF is a function of eight variables: two positions on the surface and two directions. To make it usable, the usual factorization separates the spatial profile from an ordinary BRDF evaluated at each endpoint:

$$S(x_i, \omega_i, x_o, \omega_o) = \frac{1}{\pi} F_t(\omega_i), R_d(|x_i - x_o|), F_t(\omega_o)$$

with F_t the Fresnel transmission coefficient at entry and exit, and R_d the diffuse spatial profile. Every term in that factorization is an approximation. Light leaving the medium in a direction far from the surface normal exits through a Fresnel boundary and is attenuated by T(θ) = 1 − F(θ), so the angular distribution of the exiting light is not Lambertian. It is Lambertian multiplied by a factor that grows toward grazing. The factor 1/π appears because the profile R_d is a reflectance distributed over area and the angular integral must normalize to the total.

The factorization assumes the entry and exit events are independent, which fails where the scattering distance is short relative to the surface curvature. A thin ear lobe or a small leaf has no room for the far-field profile to develop; the light sees the whole object.

The coefficients that describe a participating medium

Inside the medium, light is described by three coefficients per position and wavelength. The scattering coefficient σ_s is the probability per unit length of a scattering event; the absorption coefficient σ_a is the probability per unit length of absorption; their sum is the extinction coefficient σ_t. The units are inverse length, so a coefficient of 100 m⁻¹ means a mean free path of one centimeter.

The scattering is not isotropic in general. A single scattering event redirects light according to a phase function, and the Henyey–Greenstein phase function is the standard one-parameter choice, with g the mean cosine of the deflection angle. A g near one means strong forward scattering: light continues roughly in the direction it was going, as in a cloud or a dense suspension. A g near zero means isotropic scattering.

For multiple scattering, the direction of a single event matters much less than the amount of direction lost per unit length. This is what the reduced scattering coefficient expresses: σ’_s = σ_s (1 − g). Forward scattering does not randomize the direction, so it contributes only weakly to diffusion. The reduced extinction coefficient σ’_t = σ_a + σ’_s is the effective rate at which light forgets its original direction.

A transport mean free path of 1/σ’_t is the distance over which the light’s direction becomes uncorrelated with where it started. For media where that distance is much larger than the wavelength but much smaller than the object, diffusion is a good description; the geometric optics of the facets above the surface is then a separate, purely specular problem layered on top.

The diffusion equation and why it applies

In the diffusive regime, the radiance inside the medium is nearly isotropic: its angular variation is small compared with its mean. That single assumption supports an expansion of the radiance in spherical harmonics truncated after the first band, which leaves one scalar field, the fluence φ, and one vector field, the flux. Eliminating the flux gives

$$\nabla \cdot \left( D \nabla \phi \right) = \sigma_a \phi, \qquad D = \frac{1}{3\sigma’_t}$$

The diffusion coefficient D has units of length. In a homogeneous medium the equation reduces to ∇²φ = φ / (D σ_a), whose solutions are exponentials with a characteristic decay length δ = √(D/σ_a), one over the square root of the product of absorption and the reduced extinction. When absorption is small, light spreads far; when absorption is large, the fluence decays within a fraction of a millimeter and the whole medium behaves like a BRDF again.

The validity condition deserves emphasis: diffusion requires σ_a ≪ σ’_s, meaning the medium is scattering-dominated, and it requires the observation point to be several transport mean free paths from any boundary or source. Near a boundary, the fluence is not isotropic because light is escaping, and the approximation degrades. Rendering a thin translucent sheet with diffusion produces results that are visibly too dark or too broad, because the slab is thinner than the assumption’s range.

/* Henyey-Greenstein phase function, normalized so that the integral over
   the sphere is one. cos_t is the cosine between incoming and outgoing. */
static double hg_phase(double g, double cos_t) {
    double g2 = g * g;
    double d = 1.0 + g2 - 2.0 * g * cos_t;
    return (1.0 - g2) / (4.0 * M_PI * d * sqrt(d));
}

/* Single-scattering albedo: the fraction of an extinction event that scatters. */
static double albedo(double sigma_s, double sigma_a) {
    return sigma_s / (sigma_s + sigma_a);
}

The single-scattering albedo is the parameter that decides how many bounces light survives. At an albedo of one, no absorption occurs and the medium conserves energy exactly; a semi-infinite slab of such a medium has a reflectance that can be computed analytically in the diffusion limit and is high regardless of the details. At an albedo of 0.5, half the light is lost per event, and the effective penetration is governed by σ_a alone.

The dipole as a boundary solution

For a semi-infinite homogeneous medium illuminated by a narrow beam, the diffusion equation with a boundary condition has a closed-form solution: the dipole. The real light source is placed at a depth equal to one transport mean free path below the surface, and a negative image source is placed above the surface at the mirror position. The fluence at a surface point at distance r is the sum of the two Green’s functions

$$\phi(r) = \frac{1}{4\pi D}\left[\frac{e^{-\sigma_{tr} d_r}}{d_r} - \frac{e^{-\sigma_{tr} d_v}}{d_v}\right]$$

where σ_tr = √(σ_a/D) is the effective attenuation coefficient and d_r, d_v are the distances from the two sources. The negative image enforces an approximately absorbing boundary: light that returns to the surface escapes rather than being reflected back into the medium.

The result is a function that peaks at r = 0 and decays roughly exponentially, with a shape that has been fitted extensively to measured data on skin and other materials. Its virtue is that it costs a handful of exponentials and a square root per lookup, and its parameters can be estimated from a small number of measurements or from tabulated optical properties.

The profile is wrong in two identifiable ways. At small r, the diffusion solution rises too steeply, because the beam’s entry is treated as a point source and the actual entry has a finite extent; this produces a bright cusp. At large r, the exponential tail decays too fast compared with measured data on strongly forward-scattering media, because the diffusion approximation underestimates the contribution of long paths. Both errors are largest exactly where the visual signature is strongest, at the sharp highlight near the entry point.

Rendering with it

Because R_d is a spatial kernel rather than a directional one, evaluating it per shading point requires an integral over the surface. In a path tracer, that means a new sampling dimension: for each exit point, a sample from the profile with the correct area density. The profile’s normalization is over area, so the estimator divides by the density in area units, not in solid angle. Confusing the two produces a scattering contribution that scales with the surface’s projected area instead of its actual area.

The alternative used by many real-time systems is to precompute the spatial profile into a small texture. A sum-of-Gaussians fit to R_d is common, because a Gaussian blur of the subsurface albedo gives a reasonable approximation and runs in a few passes. The Gaussians are fitted to match the profile’s moments rather than its exact tail, so the approximation is good for soft materials and poor for highly absorptive ones where the profile is nearly a delta.

Physical plausibility checks apply here as they do for BRDFs. The total diffuse reflectance of a semi-infinite medium must be below one, and a directional albedo integration over the profile’s area must be consistent with it. A profile fitted purely for appearance can easily violate this, which then appears as a surface that seems brighter than the illumination around it.

The interaction with the specular layer is the part most often mishandled. A translucent material has a dielectric boundary on top: light reflects at the surface by Fresnel reflectance and only the transmitted share enters the medium. The exiting light is similarly attenuated. A model that adds the subsurface profile on top of a full-strength diffuse BRDF double-counts the energy that the surface already reflected, and the result is bright, waxy, and impossible to fix by tuning the profile.

Where the model gives way

Three regimes fall outside the diffusion model. Scattering that is weak relative to absorption produces single-scattering transport, where the phase function matters directly and the microfacet view of the surface is a better description than a volumetric one. Scattering that is comparable to the object size turns the object into a thin slab, where the boundary conditions dominate and a full transport solution is required. And scattering at sizes near the wavelength produces structural color, where the medium stops having a phase function at all.

For media with internal structure (layered skin, marble with veins, foliage), the homogeneous coefficients can still be assigned per region, and the diffusion solution applied piecewise, but the interfaces between regions need their own boundary conditions. The usual practice is to ignore internal interfaces and fit the parameters to the observed appearance, which works visually and destroys the connection between the parameters and the material’s actual composition. That trade is legitimate as long as it is understood: the diffusion model’s parameters are then descriptive, not physical.