Diffusion-based Methods. The diffusion equation can be solved using finite element methods which can handle heterogeneous materials and complex boundaries, but which require building the finite elements and solving the diffusion equation numerically, which is usually impractical for complex scenes. (Arbree, Walter, Bala. Heterogeneous subsurface scattering using the finite element method.)

All of these methods use the classical diffusion theory approximation, which suffers from significant errors in near-source and high-absorption regimes.

  1. Diffusion dipole introduced by Jensen  (Jensen, Marschner,Levoy, Hanrahan. A practical model for subsurface light transport.)

  1. Multipole to accurately handle thin and multi-layer materials. (Donner, Jensen. Light diffusion in multilayered translucent materials)

  1. Quantized Diffuse.
  1. Modified diffusion-type equation from neutron transport which explicitly decouples singlyand multiply-scattered light.
  2. More accurate exitance calculation from optics.

D’Eon and Irving also promoted the use of an extended source term instead of approximating it as an impulse at a single depth like the dipole. Since no closed-form solution exists to the extended source integral, and available numerical approaches were expensive, they approximated the resulting diffusion profile as a sum of Gaussians. The quantization in QD has many practical drawbacks, including complicated and error-prone implementation, numerical instability, and high computational complexity due to the large number of Gaussians needed for accurate profiles.

Monte Carlo Methods. At generally higher cost, especially for highly scattering materials.

Brute-force Monte Carlo particle tracing (Donner, Lawrence, Ramamoorthi. An empirical bssrdf

model) to tabulate an empirical BSSRDF model which handles oblique illumination and a broader range of albedo values than possible with diffusion. Kulla and Fajardo (Importance sampling techniques for path tracing in participating media) recently proposed equiangular and decoupled importance sampling strategies to reduce noise when path tracing single-scattering.

Monte Carlo algorithms that simulate multiple scattering more efficiently typically rely on some form of light path caching:

Generalized photon mapping to participating media by caching photons in volumes. (Jensen and Christensen. Efficient simulation of light transport in scenes with participating media using photon maps.)

Beam radiance estimate. Showed significant quality and performance improvements by considering the entire length of camera rays (the beam radiance estimate) or light paths (photon beams) when performing density estimation (Later, Jarosz et al. The beam radiance estimate for volumetric photon mapping.  A comprehensive theory of volumetric radiance estimation using photon points and beams). This generalization relates to “track-length estimators” used in the neutron transport field for some time (Spanier. Two pairs of families of estimators for trans-

port problems. )

Turn photons into “virtual” light sources instead of using density estimation. The rendering pass then simulates an extra bounce of illumination by connecting these lights to shading points using visibility queries.  (Krivanek, Hasan, Arbree, … Optimizing realistic rendering with many-light methods.)

Virtual ray light. Novak et al. developed an improved many-light method for media which, instead of turning photons into virtual point sources, turns photon beams into extended “virtual ray light” (VRL) sources—providing faster convergence and reduced singularities for multiple scat-

tering (Novak, Virtual ray lights for rendering scenes with participating media.)

Hybrid Methods. 

A few hybrid methods have combined diffusion and Monte Carlo techniques.

Bounce Path Tracing approximation. Hybrid approach which performs a few bounces of path tracing from the eye, and approximates subsequent bounces with the classical dipole. (Li, Pellacini, Torrance.  A hybrid Monte Carlo method for accurate and efficient subsurface scattering)

Photon diffusion (PD) performs a photon tracing pass and then interprets the photons as diffusion sources. Like QD, this accounts for the extended incident beam of light, and further accounts for oblique incident illumination and internal blockers (approximately), but it forfeits accurate multilayered material support and uses the inferior classical diffusion model. (Donner and Jensen. Rendering translucent materials using photon diffusion.)

Photon Beam Diffusion computes diffusion from photon sources, but operate on continuous photon beams [JNSJ11] instead of discrete photon points. This is similar in principle to many-light methods, but instead of computing exact one-bounce illumination from the photon sources, these methods compute approximate multi-bounce illumination using the diffusion approximation from photon points (for PD) or from photon beams (for our method, PBD). Photon Beam Diffusion can compute the same extended source integral numerically in a straightforward and very efficient way—eliminating the need for quantization and the sum of Gaussians altogether. We show that with just 3–5 numerical samples we can obtain a profile of equal quality to QD using 20–66 Gaussians, providing more flexibility and allowing us to easily approximate more general participating media transport than possible with QD by combining with other Monte Carlo rendering methods. (Jarosz, Nowrouzezahraim Sadeghi, Jensen.  A comprehensive theory of volumetric radiance estimation using photon points and beams. )

Photon Beam Diffusion.  (Ralf Habel, Per H. Christensen, Wojciech Jarosz. Photon Beam Diffusion: A Hybrid Monte Carlo Method for Subsurface Scattering)

Light Transport in Scattering Media

The change in radiance L in direction ω is the sum of three terms: a decrease in radiance dictated by the extinction coefficient ( σt = σs + σa ), and an increase due to the source function Q and the in-scattering integral on the second line.

The BSSRDF and the Searchlight Problem

Outgoing radiance, Lo , at position and direction (xo , ωo ) as a convolution of the incident illumination, Li , and the BSSRDF, S , over all incident positions and directions

For efficiency, S is often decomposed into reduced radiance, single scattering, and multi-scattering terms, S = S(0) + S(1) + Sd

Consideration order:

  1. Beam of light is incident at the origin on a semi-infinite planar homogeneous medium.
  2. Photons refract through a Fresnel boundary and travel in the downward direction until they are scattered by the medium, and ultimately get absorbed or escape the material.
  3. The distribution of photons exiting the upper boundary forms a spatial reflectance profile R d (x) , which is radially symmetric (1D) for normally-incident light, Rd(x) = Rd( x )

The reflectance profile is now centered at the incident light position, xi

Ft is the Fresnel transmittance

4Cφ is a constant needed for normalization

To evaluate this expression efficiently, most methods have relied on diffusion theory to obtain an analytic approximation for Rd .

Improved Diffusion Model and Extended Source

The reflectance profile Rd can be obtained by considering an exponentially-decreasing source beam within the material and expressing its scattering and propagation using the diffusion


  - is the extended source function

  - is the radiant intensity on the surface due to a dipole along the beam at   xr(t) = tω

The radiant intensity depends on fluence and vector irradiance terms, defined as:

where we use the shorthand dr(t) = d(x, xr(t)) and analogously for dv , and where zr(t) = xr(t) · n is the depth of the real source.

In classical diffusion theory, the diffuse reflectance profile Rd (x) depends only on the surface flux (the gradient of the fluence along the surface normal) and not the fluence φ(x) itself due to Fick’s law!!!!

with factors Cφ = 0 and CE = 1 so only vector flux, and not fluence, contributes to the diffuse reflectance in the Classical Diffuse Theory!

Diffusion Equation -  Compared to the classical theory, new diffusion equation multiplies

the source term by the reduced albedo α . This multiplication is a result of the exclusion of single scattering, which provides more accurate solutions near the source and for low albedo.

Green’s Function Isotropic source - Grosjean proposed a different approximation for the fluence due to an isotropic point source (monopole) in an infinite medium:

This is the sum of the exact single scattering and approximate multi scattering using Grosjean’s modified diffusion coefficient.

The single-scattering term is not considered part of the diffusion theory and is therefore handled separately. That’s why:

Reflection Parameter - For semi-infinite configurations, the boundary condition can be handled in an approximate fashion with the “method of images” by placing a mirrored negative source outside the medium for every positive source inside the medium. The mirroring is performed above an extrapolated boundary such that the fluence is zero at a distance zb = 2AD above the surface.

This offset takes into account the index of refraction mismatch at the boundary through

the reflection parameter A .

 where η is the relative index of refraction at the surface and Ci is the i-th angular moment of the Fresnel function. Analytic solution to the angular moments are

defined. Approximations to the angular moments, including constant factors used in the calculations:

Diffusion Coefficient -

Radiant Exitance. Instead of using Fick’s Law, which relies solely on vector flux, d’Eon and Irving use a Robin boundary condition introduced (by Kienle and Patterson. Boundary conditions for diffusion of light.)  which uses a linear combination of the fluence and its derivative, the vector flux. The result is a more accurate radiant

exitance, defined as:


Radiant Exitance from Extended Source.

To obtain the diffuse reflectance profile on the surface we need to instead integrate equation  

over the beam:

D’Eon and Irving  noted that

 has no closed-form solution, and proposed to approximate it using a sum of Gaussians (Quantized Diffuse):

where x = x , (w R (i) w i ) are the weights, and v i are the variances of normalized 2D Gaussians G 2D . The equations necessary to obtain the weights and variances of the Gaussians are themselves fairly complex summations of integrals depending on further weights wφR (v, i), wER(v, i) .

In photon beam diffusion is showed that it is efficient to numerically approximate Equation:

using Monte Carlo integration with exponential and equiangular importance sampling (Kulla, Fajardo. Importance sampling techniques for path tracing in participating media.) (Novak. Virtual ray lights for rendering scenes with participating media.)