Physically based rendering: Monte Carlo light transport and variance reduction

Why this subject matters now

For thirty years, physically based rendering was a specialist craft practiced offline. A film frame took hours on a farm, and the mathematics of unbiased light transport was something a handful of studios and a few graphics labs cared about. Two things changed. First, real-time hardware ray tracing arrived (RTX in 2018, then every console and phone GPU), which turned path tracing from an offline luxury into a latency-bound problem that a driver has to solve in milliseconds, and made denoising and sample-efficient estimators a first-class engineering concern rather than a footnote. Second, and less obviously, the entire apparatus of Monte Carlo light transport turned out to be the correct language for a large slice of modern computer vision and generative modeling. Neural radiance fields solve a discretized volume rendering equation. Inverse and differentiable rendering propagate gradients through the exact transport integral to recover geometry, materials, and lighting from photographs. The variance-reduction vocabulary, importance sampling, control variates, Russian roulette, stratification, reappears verbatim in the training of diffusion models and in reinforcement learning. A practitioner today is expected to know why radiance and not irradiance is the quantity you propagate, why an unbiased estimator can still be useless if its variance is unbounded, why multiple importance sampling was the single most important idea in Veach's thesis, and how a renderer decides, per bounce, whether to keep tracing. Those are the questions this page answers, and they are the same questions whether the output is a film frame or the forward pass of a neural renderer.

The reference this page is written against is Pharr, Jakob, and Humphreys, Physically Based Rendering: From Theory to Implementation, whose fourth edition is free online and whose source is the pbrt renderer. Where a result originates in a paper, that paper is cited. The deepest single source for the variance-reduction half is Veach's 1997 dissertation, which introduced multiple importance sampling, bidirectional path tracing in its modern form, and Metropolis light transport.

Radiometry, the quantities you actually propagate

Light transport is bookkeeping of radiant power, and the first job is to define the quantities precisely, because the entire theory rests on one of them, radiance, being the right thing to compute. The base quantity is radiant flux (power) \( \Phi \), measured in watts, the energy per unit time crossing a surface or leaving a source. Every other radiometric quantity is a density of flux with respect to area, direction, or both.

Solid angle and the projected-area measure

Directions are points on the unit sphere, and the natural measure on directions is the solid angle. A surface patch of area \( dA \) at distance \( r \) from a point, whose normal makes angle \( \theta \) with the line of sight, subtends the solid angle

measured in steradians (sr). The numerator \( dA\cos\theta \) is the patch's area projected onto the plane perpendicular to the line of sight, and dividing by \( r^2 \) makes the measure scale so that the whole sphere has \( 4\pi \) sr and a hemisphere \( 2\pi \). In spherical coordinates \( (\theta,\phi) \) about a surface normal, the differential solid angle is \( d\omega = \sin\theta\, d\theta\, d\phi \), the Jacobian of the sphere. This is the form you integrate against. The projected form and the spherical form are the same measure written two ways, and moving between them, the substitution \( d\omega = \cos\theta\, dA / r^2 \), is exactly what converts the rendering equation between its hemispherical and its area forms below.

Irradiance and radiance

Irradiance \( E \) is flux per unit area arriving at a surface, \( E = d\Phi/dA \), in \( \mathrm{W\,m^{-2}} \). It is a scalar attached to a point and tells you nothing about the directional distribution of the arriving light. Radiance \( L \) is the finer quantity, flux per unit projected area per unit solid angle,

in \( \mathrm{W\,m^{-2}\,sr^{-1}} \). The projected area \( dA^{\perp} = \cos\theta\, dA \) in the denominator is what makes radiance a function of a beam rather than of a surface orientation. It measures power per unit area as seen from the direction \( \omega \). Irradiance is the directional integral of incident radiance weighted by the cosine that projects each beam onto the surface,

the cosine appearing because a beam arriving at grazing incidence spreads its power over a larger surface footprint. This cosine is not a modeling choice. It is the geometric fact that flux density on a surface depends on the angle of arrival, and it is the same cosine that will sit inside the rendering equation.

Invariance along a ray makes radiance the right quantity

The reason renderers propagate radiance and not irradiance or flux is that radiance is constant along a ray in a vacuum. This is what lets a path tracer treat "the radiance leaving point \( x \) toward \( x' \)" and "the radiance arriving at \( x' \) from \( x \)" as the same number, and it is what makes a camera pixel a well-defined integral of surface radiance. Consider two patches \( dA_1 \) at \( x_1 \) and \( dA_2 \) at \( x_2 \), separated by distance \( r \), exchanging power along the ray joining them, with normals at angles \( \theta_1, \theta_2 \) to that ray. The power leaving patch 1 toward patch 2 is, by the definition of radiance,

so \( d^2\Phi = L_1 \, \dfrac{\cos\theta_1\cos\theta_2\, dA_1\, dA_2}{r^2} \). This expression is symmetric in the two patches. Writing the same power as radiance leaving patch 2's solid angle toward patch 1 gives \( d^2\Phi = L_2 \, \dfrac{\cos\theta_1\cos\theta_2\, dA_1\, dA_2}{r^2} \). Conservation of energy in a vacuum (no absorption or scattering between the patches) forces the two to be equal, hence \( L_1 = L_2 \). Radiance does not attenuate with distance. The quantity that does fall off, irradiance, loses the factor \( \cos\theta_2\, dA_2/r^2 \) as the receiver recedes, which is the inverse-square law. Radiance is spared it because the same \( r^2 \) that shrinks the received solid angle is what the projected-area measure already divides out. A camera sensor integrates incident radiance over the pixel's area and over the solid angle of the aperture, weighted by a sensor importance function \( W_e \). Because radiance is constant along each ray, that integral can be evaluated at the first surface the ray hits, which is the whole reason ray tracing works. The symmetric factor \( \cos\theta_1\cos\theta_2 / r^2 \) that appeared here is the geometry term, and it will reappear the moment the rendering equation is written over surfaces.

The rendering equation

The rendering equation, introduced by Kajiya in 1986, is the steady-state energy balance for radiance at a surface point. It says the radiance leaving a point in a given direction is what the point emits plus what it reflects, and reflection is an integral over all incoming directions of incident radiance filtered by the surface's bidirectional reflectance distribution function (BRDF). The BRDF itself, its definition, reciprocity, and the common analytic models, is derived on the rendering foundations page. Here it is a given function \( f_r(x,\omega_i,\omega_o) \) with units \( \mathrm{sr^{-1}} \) that this page's job is to integrate against.

Hemispherical form

Fix a surface point \( x \) with outward normal \( \mathbf n \) and an outgoing direction \( \omega_o \). The outgoing radiance is emitted radiance plus reflected radiance,

where the integral runs over the hemisphere \( \mathcal H^2 \) about \( \mathbf n \), \( \theta_i \) is the angle between \( \omega_i \) and \( \mathbf n \), and \( \cos\theta_i \) is the same projected-solid-angle cosine from the irradiance definition. The BRDF is defined precisely so that this holds. It is \( f_r = dL_o / dE_i \), the ratio of the differential reflected radiance to the differential irradiance that caused it, which is why the integrand is \( f_r \cdot L_i \cos\theta_i \) rather than \( f_r \cdot L_i \). What makes this an equation rather than a formula is the coupling. The incident radiance is the outgoing radiance of whatever surface is visible along \( \omega_i \),

where \( r(x,\omega_i) \) is the ray-cast function returning the nearest surface point seen from \( x \) in direction \( \omega_i \). Substituting one into the other exposes the structure. \( L_o \) appears on both sides, so this is a Fredholm integral equation of the second kind, \( L = L_e + \mathcal T L \), with the light-transport operator \( \mathcal T \) doing one bounce of reflection. Its formal solution is the Neumann series \( L = \sum_{k=0}^{\infty} \mathcal T^k L_e = L_e + \mathcal T L_e + \mathcal T^2 L_e + \cdots \), whose \( k \)-th term is the light that has bounced exactly \( k \) times. Each application of \( \mathcal T \) adds one more nested hemispherical integral, so the \( k \)-bounce term is a \( 2k \)-dimensional integral, and this ballooning dimensionality is precisely what forces Monte Carlo rather than quadrature.

Three-point (area) form and the geometry term

The hemispherical form integrates over directions, which hides the fact that each direction corresponds to a point on some other surface. Converting the measure from solid angle to surface area, using \( d\omega_i = \cos\theta_y\, dA_y / \lVert x - y\rVert^2 \) from the radiometry section (now with \( \theta_y \) the angle at the far surface \( y \)), rewrites the reflected term as an integral over all surfaces \( \mathcal M \) in the scene. Writing radiance between two points as \( L(y\!\to\!x) \) for the radiance leaving \( y \) toward \( x \), the equation becomes

with the geometry term

where \( \theta_x \) is the angle at \( x \) between \( \mathbf n_x \) and the segment to \( y \), \( \theta_y \) the corresponding angle at \( y \), and \( V(x,y)\in\{0,1\} \) the visibility function that is \( 1 \) when \( x \) and \( y \) see each other and \( 0 \) otherwise. This is the same symmetric \( \cos\theta_x\cos\theta_y/r^2 \) factor that fell out of the radiance-invariance derivation. The area form simply names it. Three surface points appear, \( x' \) (where light goes), \( x \) (where it reflects), and \( y \) (where it came from), which is why this is called the three-point form. Its value for a renderer is that it makes the light source an explicit domain of integration. Sampling a point on a light and connecting to it, next-event estimation, is natural in the area form and awkward in the hemispherical form, because in the area form the light is a region of \( \mathcal M \) you can sample directly.

Path-integral formulation

Veach's path-integral formulation, from his 1997 dissertation, unrolls the Neumann series into a single integral over the space of all light-carrying paths, and it is the formulation that makes bidirectional methods and Metropolis light transport possible. A measurement (the value of one pixel) is written as an integral of a contribution function over path space,

where a path \( \bar x = x_0 x_1 \cdots x_k \) of length \( k \) is a sequence of \( k+1 \) surface vertices, \( x_0 \) on a light and \( x_k \) on the sensor, and \( d\mu \) is the product area measure \( dA(x_0)\cdots dA(x_k) \). The contribution of a single path factorizes into an emission term, a product of BRDF-and-geometry factors for the interior bounces, and a sensor-importance term,

Reading it left to right, light of radiance \( L_e \) leaves the source vertex \( x_0 \), the geometry term \( G(x_0,x_1) \) carries it to the first surface, each interior vertex reflects it with a BRDF and another geometry term, and the sensor response \( W_e^{(j)} \) reads off how much of the arriving radiance lands in pixel \( j \). The power of this form is conceptual and practical at once. Conceptually, the entire transport problem is now a single (very high-dimensional, variable-length) integral, so any Monte Carlo estimator for integrals applies directly. Practically, it exposes that a path can be built from either end, or from both and connected in the middle, because the contribution function is a symmetric product of local factors. This is exactly what bidirectional path tracing exploits, and what multiple importance sampling was invented to combine. The rest of this page is, in effect, the study of good probability densities on this path space \( \Omega \).

Monte Carlo integration

Because the transport integral is high-dimensional (a \( k \)-bounce path lives in a \( 2k \)-or-more-dimensional space) and the integrand has discontinuities from visibility, deterministic quadrature is hopeless. A tensor-product rule with \( m \) nodes per axis costs \( m^d \) points, and its error does not even converge for discontinuous integrands. Monte Carlo integration trades a deterministic error for a random one whose rate is independent of dimension, which is the entire reason it dominates rendering.

The estimator, unbiasedness, and variance

To estimate \( I = \int_{\mathcal D} f(x)\, dx \), draw independent samples \( X_1,\dots,X_N \sim p \) from a density \( p \) that is positive wherever \( f \) is nonzero, and form the estimator

Its expectation is the integral for any valid \( p \),

so the estimator is unbiased. The single condition is \( p(x) > 0 \) whenever \( f(x)\neq 0 \). A density that is zero on part of the support where \( f \) is nonzero makes the estimator biased, because that region contributes to \( I \) but can never be sampled. The variance of one sample is

and because the \( N \) samples are independent, the estimator's variance is that divided by \( N \),

The root-mean-square error falls as \( O(N^{-1/2}) \), and the crucial fact is that the rate has no \( d \) in it. Doubling the dimension does not change the exponent, only the constant \( \sigma \). This is why quadrature wins in one or two dimensions and loses badly by ten. The price is the slow rate, \( N^{-1/2} \) means halving the error costs four times the samples, so every trick in the rest of the page is about shrinking the constant \( \sigma \) rather than improving the exponent.

Efficiency and control variates

Variance alone is the wrong thing to minimize, because an estimator that halves variance but triples the cost per sample is a loss. The quantity that actually matters is efficiency, the reciprocal of the product of variance and cost,

where \( T \) is the time to draw and evaluate one sample. Two estimators are compared at equal time by comparing \( \Var\cdot T \), and a variance-reduction trick is worth adopting only when it lowers that product. This is why Russian roulette, which raises variance, can raise efficiency. It lowers \( T \) by more than it raises \( \Var \). A complementary variance-reduction tool that does not touch the sampling density is the control variate. If a function \( g \) with known integral \( \int g = G \) approximates the integrand \( f \), then for any constant \( c \)

is unbiased for \( I \), because the subtracted term has expectation \( c\,G \) which is added back. Its variance is \( \frac1N\Var[(f-cg)/p] \), minimized over \( c \) at \( c^\* = \Cov[f/p,\,g/p]/\Var[g/p] \), at which the variance is reduced by the factor \( (1-\rho^2) \) where \( \rho \) is the correlation between \( f/p \) and \( g/p \). A control variate that correlates strongly with the integrand, an analytic approximation to the direct lighting, say, removes most of the variance, and the residual is only the part of \( f \) that \( g \) failed to predict. Importance sampling reshapes where samples land, while control variates subtract a known baseline. The two compose, and both leave the estimator unbiased.

Importance sampling

The variance \( \sigma^2 = \int f^2/p - I^2 \) depends on the density \( p \), and the estimator is unbiased for any valid \( p \), so \( p \) is a free lever for reducing variance. Importance sampling is the choice of \( p \) to make the ratio \( f/p \) as flat as possible, concentrating samples where the integrand is large. There is a provably optimal choice, and it drives the variance to zero.

The optimal density

Minimize \( \int f^2/p \) (the only \( p \)-dependent part of the variance) subject to \( \int p = 1 \). Apply the Cauchy-Schwarz inequality to the pair \( |f|/\sqrt p \) and \( \sqrt p \),

So \( \int f^2/p \ge (\int|f|)^2 \) always, and equality holds when \( |f|/\sqrt p \propto \sqrt p \), that is \( p \propto |f| \). For a nonnegative integrand this gives the optimal density \( p^\*(x) = f(x)/I \), for which \( f/p^\* = I \) is constant and the variance is

Zero variance means a single sample gives the exact answer. This optimal density is unusable directly, because normalizing it requires \( I \), the very number being sought, but it is the target to approximate. It says the way to reduce variance is to sample from a density shaped like the integrand, and every practical importance-sampling scheme is an approximation to \( p\propto f \) using the factors of the integrand you can sample from analytically, such as the BRDF, the cosine, and the light's emission.

Cosine-weighted hemisphere sampling

The reflected-radiance integrand always carries a \( \cos\theta_i \) factor, so for a diffuse surface (constant BRDF \( f_r = \rho/\pi \)) the integrand is proportional to \( \cos\theta_i \), and the variance-optimal density is \( p(\omega) \propto \cos\theta \). Normalizing over the hemisphere,

so the cosine-weighted density is \( p(\omega) = \cos\theta/\pi \). To sample from it, reduce to the marginal in \( \theta \) by integrating out \( \phi \), which gives \( p(\theta) = \int_0^{2\pi} \frac{\cos\theta}{\pi}\sin\theta\, d\phi = 2\cos\theta\sin\theta = \sin 2\theta \). Its cumulative distribution is

Inverting by setting \( P(\theta) = \xi_1 \) with \( \xi_1\sim\mathrm{Unif}[0,1] \) gives \( \sin\theta = \sqrt{\xi_1} \), hence \( \cos\theta = \sqrt{1-\xi_1} \). The azimuth is uniform, \( \phi = 2\pi\xi_2 \). This is the analytic inversion. Malley's method reaches the same distribution geometrically by sampling a point uniformly on the unit disk and projecting it up to the hemisphere, because the uniform-disk radius \( r = \sqrt{\xi_1} \) equals \( \sin\theta \) automatically. The payoff is that for a diffuse surface with constant incident radiance, the estimator \( f_r L_i \cos\theta / p = (\rho/\pi) L_i \cos\theta / (\cos\theta/\pi) = \rho L_i \) is constant, so the estimator has zero variance, exactly the optimal-density result specialized to the cosine.

Sampling a BRDF lobe

For glossy surfaces the integrand is peaked around the specular direction, and the density should follow the BRDF lobe. Take the normalized (modified) Phong model, \( f_r(\omega_i,\omega_o) = \frac{n+2}{2\pi}\cos^n\alpha \), where \( \alpha \) is the angle between \( \omega_i \) and the mirror-reflection direction \( \omega_r \) of \( \omega_o \), and \( n \) is the shininess exponent. A density proportional to the lobe, \( p(\omega) \propto \cos^n\alpha \), normalizes over the hemisphere about \( \omega_r \) as

The marginal in \( \alpha \) is \( p(\alpha) = (n+1)\cos^n\alpha\sin\alpha \), with cumulative distribution \( P(\alpha) = 1 - \cos^{n+1}\alpha \). Inverting \( P(\alpha) = \xi_1 \) gives

a direction sampled about \( \omega_r \) that is then rotated into world space. The estimator weight is the integrand over the density,

so the peaked \( \cos^n\alpha \) cancels and only the mild ratio \( (n+2)/(n+1) \) and the geometric cosine \( \cos\theta_i \) remain. The variance that would have come from the sharp lobe is gone. The same recipe applies to the physically-based microfacet models. For the GGX (Trowbridge-Reitz) normal distribution \( D(\mathbf h) = \frac{\alpha^2}{\pi\,(\cos^2\theta_h(\alpha^2-1)+1)^2} \), one importance samples the half-vector \( \mathbf h \) from \( D(\mathbf h)\cos\theta_h \) via \( \cos\theta_h = \sqrt{\frac{1-\xi_1}{\xi_1(\alpha^2-1)+1}} \), then reflects \( \omega_o \) about \( \mathbf h \) to obtain \( \omega_i \), with a Jacobian \( \partial\omega_h/\partial\omega_i = 1/(4\, \omega_i\!\cdot\!\mathbf h) \) that must be folded into the density. The details are in the rendering foundations page and in pbrt. The structure is identical to Phong, sample the lobe, cancel the peak.

Multiple importance sampling

Importance sampling with a single density fails when the integrand is a product of factors that peak in different places, which is the normal case in rendering. The direct-lighting integrand is \( f_r\,L_i\cos\theta \), and the BRDF may be sharp while the light is small and elsewhere. A density shaped like the BRDF is efficient for a glossy surface under a large light but breaks down for a small bright light the BRDF does not point at. A density shaped like the light is the reverse. Multiple importance sampling (MIS), Veach's central contribution, combines several sampling strategies into one estimator that is nearly as good as the best strategy at every point, and provably not much worse than the best combination.

The estimator and the balance heuristic

Given \( n \) sampling techniques with densities \( p_1,\dots,p_n \), taking \( n_i \) samples from technique \( i \), the multi-sample MIS estimator is

with weighting functions \( w_i \). This is unbiased for any weights satisfying \( \sum_{i} w_i(x) = 1 \) wherever \( f(x)\neq 0 \) (and \( w_i(x)=0 \) where \( p_i(x)=0 \)), because then

The question is which weights minimize variance. Veach's balance heuristic sets

which weights each technique at a point by how likely that technique was to have generated the point. Substituting, the term for technique \( i \) becomes \( w_i f/p_i = f \cdot n_i p_i / (n_i \sum_k n_k p_k) \), so the effective density seen by the whole estimator is the mixture \( \bar p(x) = \frac{1}{N}\sum_k n_k p_k(x) \) with \( N=\sum_k n_k \). The balance heuristic is therefore exactly the estimator you would get by sampling from the mixture density, which is why it is close to optimal. It behaves like a single importance sampler whose density is the best available blend.

Veach's variance bound and the power heuristic

Veach proved that the balance heuristic is nearly optimal in a precise sense. No other unbiased combination of the same techniques can have much lower variance. Specifically, if \( \hat I_b \) is the balance-heuristic estimator and \( \hat I_w \) is the estimator with any other valid weights, then

The additive slack shrinks as the number of samples grows and vanishes relative to the signal, so the balance heuristic cannot be beaten by more than a controlled amount. In practice the residual variance comes from points where two techniques both have moderate density. The power heuristic sharpens the weights to push more weight onto the technique that is clearly best,

with \( \beta = 2 \) the standard choice. Setting \( \beta=1 \) recovers the balance heuristic and \( \beta\to\infty \) recovers a hard maximum (pick the single densest technique). The power heuristic with \( \beta=2 \) is the default in pbrt and Mitsuba because it noticeably reduces variance in the common case of two comparable techniques while remaining unbiased for exactly the same reason, the weights still sum to one.

Path tracing

Path tracing is the Monte Carlo estimator for the Neumann series of the rendering equation. To estimate the radiance toward the camera, trace a random walk of bounces through the scene, at each bounce sampling one outgoing direction, and accumulate the emitted radiance found along the way. The naive form, sampling one direction per bounce and only registering light when a path happens to hit an emitter, is unbiased but very noisy, because most random directions miss the lights entirely. The two ideas that make it practical are next-event estimation and Russian roulette.

The algorithm and next-event estimation

Consider building a path from the camera. At the first hit \( x_1 \), the outgoing radiance splits into an emitted term and a reflected term. The reflected term is a hemispherical integral estimated by sampling a direction \( \omega_i \) from a density (cosine or BRDF), recursing to get \( L_i(x_1,\omega_i) \), and weighting by \( f_r\cos\theta_i/p \). Pure path tracing stops there and relies on the recursion eventually striking a light. Next-event estimation (NEE), also called direct light sampling, adds at every bounce an explicit connection to a point sampled on a light source, evaluated in the three-point form so the light is an area to integrate over. At vertex \( x \) the direct term is

with \( y \) sampled from a density \( p_{\mathcal A} \) over the light's surface and \( G \) including the visibility test (a shadow ray). Because NEE samples the light directly, it captures small bright sources that random BRDF bounces would almost never hit. Because BRDF sampling captures sharp specular highlights that area-light sampling handles poorly, the two are combined with MIS, weighting each connection by the balance or power heuristic over \( p_{\mathcal A} \) and the BRDF density. This pairing, NEE plus BRDF sampling combined by MIS, is the workhorse of every production path tracer, and it is the canonical application of the MIS derivation above.

Russian roulette

The Neumann series is infinite, so a path could in principle bounce forever, and every renderer must terminate paths. Truncating at a fixed depth introduces bias (it drops the tail of the series). Russian roulette terminates paths stochastically while keeping the estimator exactly unbiased. Replace an estimator \( F \) (the remaining radiance to be gathered from continuing the path) by

for a survival probability \( q\in(0,1] \). Its expectation is unchanged,

so terminating paths this way costs no bias. The price is variance. The second moment is \( \E[F_{\mathrm{rr}}^2] = q\,(F/q)^2 = \E[F^2]/q \) for a deterministic \( F \), so the variance increases by the additive amount

which for a constant \( F=X \) is exactly \( X^2(1-q)/q \). A simulation with \( X=5 \), \( q=0.6 \), over \( 2\times 10^6 \) trials, returned mean \( 5.006 \) (confirming unbiasedness) and variance \( 16.66 \) against the predicted \( X^2(1-q)/q = 25\cdot 0.4/0.6 = 16.67 \). The tradeoff is the whole point. Setting \( q \) small kills paths early and saves work but inflates variance, while \( q=1 \) never terminates. The standard policy ties \( q \) to the path's accumulated throughput \( \beta \) (the product of BRDF-over-density weights so far), often \( q = \min(1,\,\max(\beta)) \), so that dim paths, which contribute little, are killed aggressively while bright paths survive. With a per-bounce survival probability \( q \), the expected number of additional bounces is geometric, \( \sum_{k\ge 0} q^k = 1/(1-q) \), so Russian roulette turns an unbounded random walk into one with a finite expected length while leaving the answer exactly right. Efficiency, \( 1/(\Var\cdot\text{cost}) \), is what is actually optimized. Roulette is worthwhile precisely when the work it saves outweighs the variance it adds.

  camera ray                                   sampled light point y
     |                                                    * L_e
     v                                                   /
   [x1] --BRDF sample--> [x2] --BRDF sample--> [x3] ----/ shadow ray (NEE)
     |                     |                     |     (G includes visibility V)
   emit?                 emit?      Russian roulette: continue with prob q,
     |                     |        weight surviving paths by 1/q  -> unbiased
   throughput beta *= f_r * cos(theta) / p    at each bounce
        

Bidirectional methods and density estimation

Path tracing builds paths only from the camera, which struggles with light that reaches the camera through a chain of hard-to-sample events, caustics (light focused through glass onto a diffuse surface), or a room lit only by light bouncing off a wall. The path-integral formulation shows why. A path's contribution is a symmetric product of local factors, so it can be built from the light side too, or from both. Bidirectional path tracing (BDPT), developed by Lafortune and Willems and independently by Veach and Guibas around 1993-1995, traces a subpath from the camera and a subpath from a light, then connects every prefix of one to every prefix of the other. A path of length \( k \) can be formed in \( k+2 \) different ways (how many vertices came from the camera versus the light), and each way is a distinct sampling technique for the same path. MIS is what combines all of them into one low-variance estimator. BDPT is strictly more robust than unidirectional path tracing because for any given light-carrying path at least one of its construction strategies is usually efficient, and MIS automatically down-weights the bad ones.

Photon mapping, from Jensen's work in the mid-1990s and his 2001 book, takes a different route. In a first pass it shoots photons from the lights and stores where they land in a spatial data structure (the photon map). In a second pass it renders from the camera and, at each visible point, estimates the incident illumination by a density estimate over the nearby photons, gathering the \( k \) nearest photons and dividing their total power by the area of the disk containing them. This is a kernel density estimate of the radiance, and it is biased but consistent. For a finite photon count the density estimate blurs illumination over a nonzero radius (bias), but as the number of photons grows and the gather radius shrinks, the estimate converges to the true radiance (consistency). This bias-versus-consistency stance is the opposite of path tracing's unbiased-but-noisy stance, and it is the defining tradeoff of density-estimation renderers. Progressive photon mapping (Hachisuka and colleagues, 2008) makes the bias vanish in the limit by shrinking the radius over successive passes, giving a consistent estimator whose error goes to zero. Photon mapping is strongest exactly where path tracing suffers, caustics and dense indirect light, because photons naturally accumulate where light concentrates.

Metropolis light transport

Metropolis light transport (MLT), from Veach and Guibas in 1997, applies the Metropolis-Hastings algorithm to path space. Instead of sampling paths independently, it runs a Markov chain over the space of light-carrying paths whose stationary distribution is proportional to each path's contribution to the image. New paths are proposed by locally mutating the current path (perturbing a vertex, adding a bounce) and accepted or rejected to preserve that stationary distribution. Because it explores path space by local moves, MLT is particularly effective on scenes where the important light-carrying paths form a thin, hard-to-find region, light squeezing through a slightly open door, since once the chain finds a bright path it stays in that neighborhood and samples it thoroughly. Its costs are the difficulty of designing good mutation strategies and non-uniform convergence (some image regions clear up faster than others). The primary-sample-space variant of Kelemen and colleagues (2002) reformulates the mutations in the unit hypercube of random numbers that drive a standard path tracer, which makes MLT much simpler to implement, at some cost in the locality that makes the path-space version so powerful.

Sampling patterns and quasi-Monte-Carlo

Everything above assumed independent uniform samples, but the samples that drive an estimator are points in a unit hypercube (one dimension per random decision), and how those points are placed changes the constant in front of \( N^{-1/2} \), sometimes the exponent itself. The goal is samples that are well spread, with no clumps and no gaps.

Stratified sampling

Stratified sampling partitions the domain into \( M \) equal cells and places one sample uniformly in each. To see why it never hurts, split the single-sample variance of plain Monte Carlo by the law of total variance over the cell index \( J \) (uniform on the \( M \) cells),

Plain Monte Carlo with \( M \) samples pays \( \tfrac1M \) of the whole thing. Stratified sampling, one sample per cell, is an independent sum across cells, so its estimator variance is \( \tfrac{1}{M^2}\sum_j \Var[f\mid J{=}j] = \tfrac1M\,\E_J[\Var[f\mid J]] \), exactly the within-cell term and nothing else. The between-cell variance, the part coming from the integrand's large-scale trend across cells, is removed for free, and since variances are nonnegative the stratified estimator is never worse. For a smooth integrand in one dimension the improvement can be large. In \( d \) dimensions a full grid needs \( M^d \) cells, so stratification is usually applied per-dimension or via Latin hypercube designs that stratify each axis independently. The catch is that stratification is only as good as the alignment between the cells and the integrand's structure, and it does not compose well across the many dimensions of a long path, which is what motivates low-discrepancy sequences.

Low-discrepancy sequences and star discrepancy

Quasi-Monte-Carlo replaces random points with deterministic sequences engineered to be more uniform than random. Uniformity is measured by the star discrepancy \( D^\*_N \), the worst-case gap between the fraction of points inside an axis-aligned box anchored at the origin and that box's volume,

the supremum over all such boxes \( B \). The Koksma-Hlawka inequality bounds the integration error by the product of the discrepancy and the integrand's total variation, \( |\hat I_N - I| \le V(f)\, D^\*_N \), so a lower-discrepancy point set integrates better for any function of bounded variation. Random points have \( D^\*_N = O(\sqrt{\log\log N / N}) \), giving the familiar \( N^{-1/2} \). Low-discrepancy sequences, the Halton sequence (radical-inverse in a different prime base per dimension) and the Sobol sequence (built from direction numbers over the base-2 digits), achieve \( D^\*_N = O((\log N)^d / N) \), which for fixed \( d \) is close to \( N^{-1} \), a full power faster than random sampling. In a renderer this shows up as noise that decreases roughly like \( N^{-1} \) rather than \( N^{-1/2} \) in the smooth parts of the image, which is why pbrt and every production renderer default to Sobol-based samplers. The catch is that the \( (\log N)^d \) constant grows with dimension, so the advantage erodes for very long paths, and the deterministic structure can produce visible aliasing patterns rather than pleasant noise.

Randomized QMC and blue noise

Owen scrambling randomizes a Sobol or Halton sequence by permuting its digits, which restores an unbiased estimator (the scrambled points are individually uniform) while preserving the low discrepancy. Owen proved scrambled nets can reach variance \( O(N^{-3}(\log N)^{d-1}) \) for smooth integrands, far below the random rate, and the scrambling turns structured aliasing back into noise that a denoiser handles well. Blue-noise sampling addresses a different objective. Rather than minimizing integration error at a single pixel, it distributes error across pixels so that the residual noise has its energy in the high frequencies, where the human visual system and modern denoisers tolerate it best. Techniques that assign correlated sample sequences to neighboring pixels so the per-pixel errors are spread as blue noise (Heitz and colleagues at Unity, 2019) produce images that look markedly cleaner at equal sample count even though the per-pixel variance is unchanged, which matters greatly for real-time rendering where the sample budget is one or two paths per pixel.

Microfacet, subsurface, and volumetric transport

The rendering equation as stated assumes light reflects at an opaque surface, but real materials transmit, scatter beneath the surface, and fill volumes. Each extension changes the integrand, not the Monte Carlo machinery, which is why the variance-reduction results above carry over unchanged.

Microfacet BSDFs and multiple scattering

Physically based surface models treat a rough surface as a distribution of tiny mirror facets, giving the microfacet BRDF \( f_r = \frac{D(\mathbf h)\, F(\omega_o,\mathbf h)\, G(\omega_i,\omega_o)}{4\cos\theta_i\cos\theta_o} \), where \( D \) is the normal distribution (GGX above), \( F \) the Fresnel reflectance, and \( G \) the shadowing-masking term that accounts for facets occluding one another. This single-scattering model conserves less energy than it should, because it ignores light that bounces between facets before leaving. At high roughness the missing multiple-scatter energy darkens the material visibly. Heitz and colleagues (2016) gave a stochastic multiple-scattering model that traces light through the microsurface and restores energy conservation, and Kulla and Conty (2017) an analytic energy-compensation term now standard in production. The extension to transmission, the BSDF, folds refraction into the same half-vector sampling with the Fresnel term deciding reflect-versus-transmit. The derivation on the rendering foundations page covers the single-bounce BRDF, and the transport machinery here integrates it.

Subsurface scattering and the BSSRDF

For translucent materials, skin, marble, milk, light enters at one point and leaves at another after scattering beneath the surface, so a BRDF (which assumes entry and exit at the same point) is wrong. The bidirectional scattering-surface reflectance distribution function (BSSRDF) \( S(x_i,\omega_i,x_o,\omega_o) \) relates outgoing radiance at \( x_o \) to incident flux at a different point \( x_i \), and the reflected radiance becomes a double integral over both the incident directions and the surface area. Jensen and colleagues (2001) introduced a practical dipole diffusion approximation that models the subsurface light distribution as a pair of point sources, reducing the expensive volumetric integral to an analytic surface term for optically dense media. Modern renderers use a photon-beam or path-traced subsurface for accuracy and the diffusion approximation for speed. The BSSRDF is the bridge between surface reflectance and full volumetric transport.

Participating media and the volume rendering equation

When light travels through fog, smoke, or tissue, it is absorbed and scattered along the way, not only at surfaces. The radiance along a ray obeys the volume rendering equation, which integrates emission and in-scattering along the ray weighted by transmittance,

where \( \sigma_a \) and \( \sigma_s \) are the absorption and scattering coefficients, \( \sigma_t = \sigma_a + \sigma_s \) the extinction, \( p \) the phase function (the volumetric analog of a BRDF), and the transmittance

is the fraction of light surviving from \( x \) to \( x_t \) without being absorbed or scattered out, the Beer-Lambert law. Monte Carlo estimators for this integral sample free-flight distances from the transmittance (which for homogeneous media inverts analytically, \( t = -\ln(1-\xi)/\sigma_t \)) and handle spatially varying media with delta tracking, a rejection scheme that adds fictitious "null" collisions to make the extinction constant and thus samplable. This exact equation is what neural radiance fields discretize. NeRF replaces the analytic media coefficients with a neural field and approximates the transmittance integral by a quadrature sum over samples \( t_i \) along the ray, \( C = \sum_i T_i\,(1 - e^{-\sigma_i \delta_i})\,\mathbf c_i \) with \( T_i = \exp(-\sum_{j<i}\sigma_j \delta_j) \), which is precisely the volume rendering equation with emission-only media and a learned \( (\sigma, \mathbf c) \). The connection is developed on the diffusion and large vision models page. The transport theory here is what that discretization approximates.

Differentiable rendering

Inverse rendering asks the reverse question. Given photographs, recover the scene parameters (geometry, materials, lights) that produced them, by gradient descent on a rendering loss. This requires differentiating the transport integral with respect to scene parameters, and the difficulty is that the integrand has discontinuities, at object silhouettes and shadow boundaries the visibility function \( V \) jumps, and those discontinuities move as the geometry moves. A naive pathwise derivative misses the contribution of the moving boundary, so the gradient is biased. Li and colleagues (2018) introduced the first correct differentiable path tracer (redner) by splitting the derivative into an interior term, the ordinary derivative of the integrand where it is smooth, and a boundary term, an integral over the silhouette edges that captures how the discontinuity sweeps across the domain, sampled by an edge-sampling procedure. Loubet and colleagues (2019) reformulated this with a reparameterization that moves the discontinuity out of the integrand so a standard automatic-differentiation pass gives the right gradient. Mitsuba 3 (Nimier-David and colleagues, 2019, built on the just-in-time Dr.Jit autodiff system) productionized differentiable and inverse rendering, with radiative backpropagation (Nimier-David and colleagues, 2020) and path replay backpropagation (Vicini and colleagues, 2021) reducing the memory of the reverse pass from linear in path length to constant by re-tracing rather than storing the path. This line of work is what connects classical light transport to the gradient-based optimization that drives modern vision and graphics learning, and it is why the same Monte Carlo variance-reduction results, now applied to gradient estimators, remain the active research frontier.

Implementation

The heart of this page is the set of Monte Carlo estimators, so the first code block is a runnable NumPy implementation of the four demonstrations whose numbers appear in the worked problems. These are the \( O(N^{-1/2}) \) convergence of a uniform estimator, the variance drop from importance sampling, the zero-variance cosine estimator, and a two-strategy MIS combination beating either strategy alone. Every number quoted above came from running exactly this code.

import numpy as np
rng = np.random.default_rng(20260724)

# ---- 1. uniform MC and O(1/sqrt N): I = int_0^1 e^x dx = e - 1 ----
true = np.e - 1.0
varf = (np.e**2 - 1) / 2 - true**2          # analytic single-sample variance = 0.24204
for N in [16, 256, 4096, 65536]:
    est = np.exp(rng.random((4000, N))).mean(axis=1)   # 4000 independent trials
    rmse = np.sqrt(np.mean((est - true) ** 2))
    print(N, "rmse", round(rmse, 6), "predicted", round(np.sqrt(varf / N), 6))

# ---- 2. importance sampling: I = int_0^1 e^(-8x) dx ----
a, b = 8.0, 6.0
true2 = (1 - np.exp(-a)) / a
Z = (1 - np.exp(-b)) / b                     # normalizer for p(x) ~ e^(-b x) on [0,1]
xu = rng.random((4000, 4096))                # uniform sampler
rmse_u = np.sqrt(np.mean((np.exp(-a * xu).mean(1) - true2) ** 2))
u = rng.random((4000, 4096))                 # inverse-CDF sample from truncated exponential
xis = -np.log(1 - u * (1 - np.exp(-b))) / b
p = np.exp(-b * xis) / Z
rmse_is = np.sqrt(np.mean(((np.exp(-a * xis) / p).mean(1) - true2) ** 2))
print("IS rmse ratio", round(rmse_u / rmse_is, 2), "(predicted ~6.81)")

# ---- 3. cosine-weighted hemisphere: E = int cos(theta) dw = pi, zero variance ----
xi = rng.random(100_000)
cos_theta = np.sqrt(1 - xi)                   # cosine-weighted: p(w) = cos/pi
est_cos = np.pi * np.ones_like(cos_theta)     # f/p = cos / (cos/pi) = pi, exact
print("cosine E", round(est_cos.mean(), 5), "var", est_cos.var())   # ~pi, ~0

# ---- 4. MIS: integrate product of two gaussians, balance heuristic ----
def g(x, m, s): return np.exp(-0.5 * ((x - m) / s) ** 2) / (s * np.sqrt(2 * np.pi))
mA, sA, mB, sB = 0.0, 0.6, 2.0, 0.4
f  = lambda x: g(x, mA, sA) * g(x, mB, sB)
pA = lambda x: g(x, mA, sA)
pB = lambda x: g(x, mB, sB)
true4 = np.exp(-0.5 * (mA - mB) ** 2 / (sA**2 + sB**2)) / np.sqrt(2 * np.pi * (sA**2 + sB**2))
n = 4096
xa = rng.normal(mA, sA, (3000, n)); xb = rng.normal(mB, sB, (3000, n))
wA = pA(xa) / (pA(xa) + pB(xa)); wB = pB(xb) / (pA(xb) + pB(xb))
est_mis = (wA * f(xa) / pA(xa)).sum(1) / n + (wB * f(xb) / pB(xb)).sum(1) / n
print("MIS rmse", round(np.sqrt(np.mean((est_mis - true4) ** 2)), 5), "true", round(true4, 6))

import jax, jax.numpy as jnp
from jax import random

# Vectorized, jit-compiled importance-sampling estimator; the same math as NumPy,
# but differentiable: grad flows through the estimator w.r.t. the pdf parameter b,
# which is the seed of differentiable rendering (grad of an integral estimate).
def is_estimate(key, a, b, N):
    Z = (1 - jnp.exp(-b)) / b
    u = random.uniform(key, (N,))
    x = -jnp.log(1 - u * (1 - jnp.exp(-b))) / b     # inverse CDF of e^(-b x) on [0,1]
    p = jnp.exp(-b * x) / Z
    return jnp.mean(jnp.exp(-a * x) / p)             # unbiased for any b > 0

key = random.PRNGKey(0)
a = 8.0
true2 = (1 - jnp.exp(-a)) / a
est = jax.jit(is_estimate, static_argnums=3)(key, a, 6.0, 1 << 20)
print("estimate", est, "true", true2)

# d(estimate)/db: the estimator value depends on the sampling density; its
# EXPECTATION does not (unbiased for all b), so E[d/db] = 0 -- a score-function
# identity. Nonzero for a finite sample, which is exactly the variance-of-gradient
# problem differentiable renderers must control.
g = jax.grad(lambda b: is_estimate(key, a, b, 1 << 18))(6.0)
print("d est / db (single sample)", g)

// Inner loop of a unidirectional path tracer with next-event estimation,
// MIS (power heuristic), and Russian roulette. Radiance for one camera ray.
// Vec3 has the usual operators; scene.intersect / sampleLight are elsewhere.
Vec3 radiance(Ray ray, const Scene &scene, RNG &rng) {
    Vec3 L(0.0), beta(1.0);         // accumulated radiance, path throughput
    bool specularBounce = true;     // count emission on the first / specular hits
    for (int depth = 0; ; ++depth) {
        Hit hit;
        if (!scene.intersect(ray, &hit)) break;         // ray escaped the scene

        // Emission: add it only when BRDF sampling brought us here (MIS pairs with NEE)
        if (hit.isEmitter && specularBounce)
            L += beta * hit.Le;

        // --- Next-event estimation: sample a light, connect with a shadow ray ---
        LightSample ls = scene.sampleLight(hit.p, rng);   // point y, Le, pdf over area
        if (ls.pdfArea > 0 && !scene.occluded(hit.p, ls.p)) {
            Vec3  wi   = normalize(ls.p - hit.p);
            float G    = fabs(dot(hit.n, wi) * dot(ls.n, -wi)) / distanceSq(hit.p, ls.p);
            Vec3  f    = hit.bsdf.eval(-ray.d, wi);        // f_r * cos folded in eval
            float pL   = ls.pdfArea;                       // light pdf (area measure)
            float pB   = hit.bsdf.pdf(-ray.d, wi) * G;     // BRDF pdf, area measure
            float w    = (pL * pL) / (pL * pL + pB * pB);  // power heuristic, beta = 2
            L += beta * w * f * G * ls.Le / pL;
        }

        // --- BRDF sampling: choose the next direction, update throughput ---
        BsdfSample bs = hit.bsdf.sample(-ray.d, rng);      // wi, f_r*cos, pdf
        if (bs.pdf <= 0.0f) break;
        beta *= bs.f / bs.pdf;                              // f_r * cos(theta) / p
        specularBounce = bs.isSpecular;
        ray = Ray(hit.p, bs.wi);

        // --- Russian roulette: unbiased stochastic termination past a few bounces ---
        if (depth > 3) {
            float q = std::min(0.95f, maxComponent(beta));  // survive w.p. q ~ brightness
            if (rng.uniform() >= q) break;                  // terminate
            beta /= q;                                       // keep the estimator unbiased
        }
    }
    return L;
}

The C++ tab is the loop that pbrt, Mitsuba, and every game-engine path tracer run per ray. Three details carry the theory. The specularBounce flag prevents double-counting emission. When NEE already sampled a light, BRDF sampling that then lands on the same light must not add the emission again, so the MIS weight on the emission term is folded in through this flag. The power-heuristic weight w combines the light and BRDF densities, both converted to the same area measure so they are comparable. And Russian roulette divides beta by the survival probability on the paths that live, which is the \( 1/q \) that keeps the estimator unbiased.

How it is done in practice

A production path tracer is the algorithm above wrapped in a great deal of engineering, and the gap between the derivation and a shipped renderer is where most of the work lives. Ray-scene intersection dominates the cost, so scenes are stored in a bounding-volume hierarchy and traversed with hand-tuned SIMD or hardware ray-tracing cores. Intel's Embree and NVIDIA's OptiX are the two libraries the field standardizes on, and on RTX hardware the BVH traversal and ray-triangle test run in dedicated silicon. The stochastic nature of Monte Carlo means the raw image is noisy at the sample counts a film or a game can afford (a film might spend thousands of paths per pixel, a real-time renderer one or two), so every modern pipeline ends in a denoiser, edge-aware filters or learned denoisers like NVIDIA's OptiX denoiser and Intel's Open Image Denoise, which turn a noisy few-sample estimate into a clean image by exploiting spatial and temporal coherence. The interaction between the sampler and the denoiser is why blue-noise sample distribution matters so much in practice. A denoiser removes high-frequency noise cleanly and low-frequency blotches poorly, so pushing the error into high frequencies is worth more than reducing its total magnitude.

Wavefront path tracing reorganizes the per-ray loop into per-stage kernels (generate rays, intersect, shade, extend) so that a GPU processes thousands of paths in lockstep without the divergence that stalls a naive megakernel. This is how real-time path tracers keep the SIMD lanes busy when different paths take different branches. Spectral rendering replaces RGB with sampled wavelengths to get dispersion and accurate color right, adding one more dimension to importance sample. And the material system, the layered, multi-lobe BSDFs of a production shader, is where a large fraction of the sampling code lives, because each lobe needs its own importance sampler and the lobes are combined with, once again, multiple importance sampling. The through-line is that the two ideas this page derives, unbiased Monte Carlo estimation and multiple importance sampling, are not academic. They are the load-bearing structure of every renderer in film and games.

Sample budgets are not spread uniformly across the image. Adaptive sampling estimates the variance per pixel or per tile as rendering proceeds and steers additional samples toward the regions that are still noisy, glossy reflections, caustics, high-contrast edges, while leaving flat, well-converged regions alone. Because the estimator variance falls as \( \sigma^2/N \), spending samples where \( \sigma \) is large maximizes the reduction in total image error per sample, which is the efficiency metric derived above applied at the granularity of a pixel. Firefly pixels, isolated bright spots from a rare high-throughput path that a low sample count failed to average out, are the visible symptom of unbounded per-sample variance, and production renderers clamp or specially handle them because a single \( f/p \) with a tiny \( p \) can dominate a pixel's mean. The recurring theme across all of this engineering is that the derivations of this page set the ceiling, an unbiased estimator with a known variance and a known cost, and the production work is the continual minimization of \( \Var\cdot T \) under that ceiling, pixel by pixel and bounce by bounce.

The current research frontier

Three directions are active. The first is neural and learned sampling, ML models that learn the importance-sampling density from the scene, so that a network proposes the next bounce direction better than any analytic BRDF sampler can. Neural importance sampling with normalizing flows (Müller and colleagues at Disney Research and ETH Zurich, 2019) learns a per-scene density that adapts to the integrand, and path guiding (Vorba and colleagues, 2014, and the practical online variant by Müller and colleagues, 2017) builds a spatial-directional data structure during rendering that steers samples toward the light, cutting variance on hard indirect lighting by large factors. The second is real-time global illumination. ReSTIR (reservoir spatiotemporal importance resampling, Bitterli and colleagues at NVIDIA and Dartmouth, 2020) reuses light samples across neighboring pixels and frames through a streaming resampling scheme, making many-light direct illumination tractable at one sample per pixel, and its path-tracing generalization (ReSTIR PT, 2022) extends the idea to full global illumination. The third is differentiable and inverse rendering, discussed above, where the frontier is unbiased and low-variance gradient estimators through discontinuities and volumes. Mitsuba 3 and Dr.Jit (Jakob and colleagues at EPFL) are the shared platform, and the open problem is estimating boundary gradients cheaply enough for practical inverse rendering of complex scenes. Across all three, the unifying observation is that the classical theory, unbiased estimation, importance sampling, MIS, is being combined with learning and with sample reuse, and the variance-reduction vocabulary of Veach's thesis is the language every one of these papers is written in.

Open source to read

• mmp/pbrt-v4 is the reference renderer accompanying the Pharr-Jakob-Humphreys book. Start with src/pbrt/cpu/integrators.cpp, where the path, bidirectional-path, and MLT integrators are implemented in the same clear style as the text. It is the single best way to see the derivations of this page turned into working code.
• mitsuba-renderer/mitsuba3 is a research renderer built for differentiable and inverse rendering on top of the Dr.Jit autodiff compiler. Read src/integrators/path.cpp for the forward path tracer and the prb integrators for path-replay backpropagation. This is the platform most differentiable-rendering papers build on.
• embree/embree is Intel's high-performance ray-triangle and BVH library. Open tutorials/pathtracer to see how a production intersection kernel is structured. It is the intersection layer under many renderers.
• appleseedhq/appleseed is a complete production-quality path tracer with a full material system and MIS-based direct lighting. It is a good middle ground between pbrt's pedagogy and a shipping renderer's complexity.
• tunabrain/tungsten is a compact, fast physically based renderer implementing path tracing, BDPT, progressive photon mapping, and primary-sample-space MLT. Its small codebase makes it the easiest place to read a real BDPT and a real MLT side by side.
• BachiLi/redner is the first differentiable path tracer with correct edge-sampling gradients. Read it alongside the 2018 paper to see how the boundary term of the differentiated rendering integral is actually sampled.

Common misconceptions

"An unbiased estimator is always better than a biased one." Bias and variance are both error. A biased-but-consistent estimator like photon mapping can produce a visually cleaner image at a fixed budget than an unbiased path tracer whose variance shows up as noise. The right metric is total error (mean squared error), and low-variance bias often wins on the caustics and dense-indirect scenes where path tracing struggles.

"More bounces always means a more correct image." The Neumann series converges, so the contribution of deep bounces decays. Past a handful of bounces the added energy is negligible for most scenes, and Russian roulette is exactly the tool that spends samples on deep bounces only in proportion to how much they contribute. Tracing to a huge fixed depth wastes work without visibly changing the image.

"Importance sampling reduces error by sampling more where the function is big." Close but backwards in emphasis. It reduces variance by making the ratio \( f/p \) flat, which means sampling in proportion to the integrand's magnitude including regions that are moderately large, not just the peak. A density concentrated only on the peak but missing a broad shoulder can have higher variance than uniform sampling, and a density that is zero where \( f \) is nonzero makes the estimator biased.

"MIS is just averaging several estimators." A plain average of a good and a bad estimator is dragged down by the bad one. MIS weights each technique per-sample by how likely it was to produce that sample, so the bad technique contributes almost nothing at points where the good one dominates. That per-point weighting, not the averaging, is what makes MIS provably close to the best combination.

"Radiance falls off with distance like everything else." Radiance is constant along a ray in a vacuum, as the two-patch derivation shows. It is irradiance and flux density that obey the inverse-square law. Confusing the two is the classic source of a renderer that is too dark or too bright by a factor of \( r^2 \).

"Quasi-Monte-Carlo is biased because the points are not random." A deterministic low-discrepancy sequence gives a deterministic (hence technically biased) estimate, but Owen-scrambled randomized QMC restores an unbiased estimator whose variance is far below plain Monte Carlo, so you keep the convergence benefit and the unbiasedness at once. Production renderers use the randomized form for exactly this reason.

Self-check

References

1. M. Pharr, W. Jakob, G. Humphreys. Physically Based Rendering: From Theory to Implementation, 4th ed. MIT Press, 2023. Free online at pbr-book.org.
2. E. Veach. Robust Monte Carlo Methods for Light Transport Simulation. PhD dissertation, Stanford University, 1997. graphics.stanford.edu/papers/veach_thesis.
3. J. T. Kajiya. The rendering equation. SIGGRAPH, 1986. doi:10.1145/15922.15902.
4. R. L. Cook, T. Porter, L. Carpenter. Distributed ray tracing. SIGGRAPH, 1984. doi:10.1145/800031.808590.
5. E. Veach, L. J. Guibas. Optimally combining sampling techniques for Monte Carlo rendering. SIGGRAPH, 1995. doi:10.1145/218380.218498.
6. E. P. Lafortune, Y. D. Willems. Bidirectional path tracing. Compugraphics, 1993.
7. E. Veach, L. J. Guibas. Metropolis light transport. SIGGRAPH, 1997. doi:10.1145/258734.258775.
8. H. W. Jensen. Realistic Image Synthesis Using Photon Mapping. A K Peters, 2001.
9. H. W. Jensen, S. R. Marschner, M. Levoy, P. Hanrahan. A practical model for subsurface light transport (the BSSRDF dipole). SIGGRAPH, 2001. doi:10.1145/383259.383319.
10. T. Hachisuka, S. Ogaki, H. W. Jensen. Progressive photon mapping. ACM Trans. Graphics (SIGGRAPH Asia), 2008. doi:10.1145/1409060.1409083.
11. C. Kelemen, L. Szirmay-Kalos, G. Antal, F. Csonka. A simple and robust mutation strategy for the Metropolis light transport algorithm. Eurographics, 2002. doi:10.1111/1467-8659.00703.
12. A. Keller. Quasi-Monte Carlo image synthesis in a nutshell. In Monte Carlo and Quasi-Monte Carlo Methods, Springer, 2013. doi:10.1007/978-3-642-41095-6_8.
13. A. B. Owen. Randomly permuted (t,m,s)-nets and (t,s)-sequences. In Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing, Springer, 1995. doi:10.1007/978-1-4612-2552-2_19.
14. E. Heitz, J. Hanika, E. d'Eon, C. Dachsbacher. Multiple-scattering microfacet BSDFs with the Smith model. ACM Trans. Graphics (SIGGRAPH), 2016. doi:10.1145/2897824.2925943.
15. B. Walter, S. R. Marschner, H. Li, K. E. Torrance. Microfacet models for refraction through rough surfaces (GGX). Eurographics Symposium on Rendering, 2007. doi:10.2312/EGWR/EGSR07/195-206.
16. T. Müller, B. McWilliams, F. Rousselle, M. Gross, J. Novák. Neural importance sampling. ACM Trans. Graphics, 2019. arXiv:1808.03856.
17. T. Müller, M. Gross, J. Novák. Practical path guiding for efficient light-transport simulation. Computer Graphics Forum (EGSR), 2017. doi:10.1111/cgf.13227.
18. B. Bitterli, C. Wyman, M. Pharr, P. Shirley, A. Lefohn, W. Jarosz. Spatiotemporal reservoir resampling for real-time ray tracing with dynamic direct lighting (ReSTIR). ACM Trans. Graphics (SIGGRAPH), 2020. doi:10.1145/3386569.3392481.
19. T.-M. Li, M. Aittala, F. Durand, J. Lehtinen. Differentiable Monte Carlo ray tracing through edge sampling (redner). ACM Trans. Graphics (SIGGRAPH Asia), 2018. doi:10.1145/3272127.3275109.
20. G. Loubet, N. Holzschuch, W. Jakob. Reparameterizing discontinuous integrands for differentiable rendering. ACM Trans. Graphics (SIGGRAPH Asia), 2019. doi:10.1145/3355089.3356510.
21. M. Nimier-David, D. Vicini, T. Zeltner, W. Jakob. Mitsuba 2: a retargetable forward and inverse renderer. ACM Trans. Graphics (SIGGRAPH Asia), 2019. doi:10.1145/3355089.3356498.
22. M. Nimier-David, S. Speierer, B. Ruiz, W. Jakob. Radiative backpropagation: an adjoint method for lightning-fast differentiable rendering. ACM Trans. Graphics (SIGGRAPH), 2020. doi:10.1145/3386569.3392406.
23. D. Vicini, S. Speierer, W. Jakob. Path replay backpropagation: differentiating light paths using constant memory and linear time. ACM Trans. Graphics (SIGGRAPH), 2021. doi:10.1145/3450626.3459804.
24. B. Mildenhall, P. P. Srinivasan, M. Tancik, J. T. Barron, R. Ramamoorthi, R. Ng. NeRF: representing scenes as neural radiance fields for view synthesis. ECCV, 2020. arXiv:2003.08934.