Sampling BRDF Lobes and Choosing Between Them

Monte Carlo evaluation of the reflectance integral replaces a sum over infinitely many directions with an average over finitely many samples. Which directions are chosen determines how many samples are needed. A uniform choice over the hemisphere is unbiased and nearly useless: for a sharp specular lobe, almost every sample lands where the BRDF is negligible, and the estimator’s variance is dominated by the rare samples that hit the lobe.

Importance sampling fixes this by drawing directions in proportion to the integrand. The ideal density is proportional to the BRDF times the cosine, because that product is the integrand; when the density matches it exactly, every sample contributes the same value and the variance is zero. Real densities approximate it, and the quality of a renderer’s sampling is largely the quality of that approximation.

The estimator and the weight it produces

With a density p(ω) over outgoing directions, the reflected radiance estimate for one sample is

$$\hat{L}_o = \frac{f_r(\omega_i, \omega), L_i(\omega), (\omega \cdot n)}{p(\omega)}$$

The value f_r cos / p is the sample weight. For a density proportional to f_r cos, the weight is a constant and every sample carries the same information. For a density that only resembles the integrand, weights vary and the variance of the estimator is the variance of the weights.

The practical test of a sampling routine is to evaluate the weight distribution, not to look at the rendered image. A routine with occasional weights in the hundreds produces fireflies, and no amount of denoising fixes the underlying problem: those samples are correct estimates of a rare bright event, and suppressing them biases the result.

For diffuse surfaces the natural density is cosine-weighted, p(ω) = cos θ / π. It is generated cheaply from two uniforms by the concentric-disk or polar mapping. Above a Lambertian BRDF the weight becomes albedo / π, which is constant, so a diffuse surface requires no more than one sample per light for exactness modulo the light’s own solid angle.

Why the half-vector is the right variable

For a specular lobe the integrand is concentrated near the mirror direction, and the width of the concentration is set by the microfacet distribution. The standard technique samples the half-vector rather than the outgoing direction: draw h from a density proportional to D(h)(h·n), then compute ω_o by reflecting ω_i about h. This is the VNDF-family idea, and the reason it works is that the distribution of visible facet normals, not the distribution of all facet normals, is what the BRDF actually weights.

Sampling all normals is biased against facets that face away from the incident direction; those facets exist on the surface but cannot be seen by the incoming light, so including them in the density under-weights the ones that matter. The correction is to multiply by the visibility factor and renormalize.

For GGX the resulting density has a closed form and an invertible cumulative distribution, which is why it is nearly universal. With α the roughness and cos θ_i the incidence cosine, the visible-normal distribution in the slope parameterization has a cumulative function along the stretched azimuth that can be inverted with square roots, producing a sample in a few arithmetic operations and no rejection loop.

/* Sample GGX visible normals. u1, u2 in [0,1); V is the view direction in a
   local frame where the macrosurface normal is +Z. Result: half-vector in Z. */
static void sample_ggx_vndf(double alpha, double u1, double u2,
                            const double V[3], double H[3]) {
    /* stretch the view vector into the hemisphere of the slope space */
    double Vh[3] = { alpha * V[0], alpha * V[1], V[2] };
    double len = sqrt(Vh[0]*Vh[0] + Vh[1]*Vh[1] + Vh[2]*Vh[2]);
    Vh[0] /= len; Vh[1] /= len; Vh[2] /= len;

    double lensq = Vh[0]*Vh[0] + Vh[1]*Vh[1];
    double T1[3] = { lensq > 0.0 ? -Vh[1]/sqrt(lensq) : 1.0,
                     lensq > 0.0 ?  Vh[0]/sqrt(lensq) : 0.0, 0.0 };
    double T2[3] = { Vh[1]*T1[2] - Vh[2]*T1[1],
                     Vh[2]*T1[0] - Vh[0]*T1[2],
                     Vh[0]*T1[1] - Vh[1]*T1[0] };

    double r = sqrt(u1), phi = 2.0 * 3.14159265358979 * u2;
    double t1 = r * cos(phi), t2 = r * sin(phi);
    double s = 0.5 * (1.0 + Vh[2]);
    t2 = (1.0 - s) * sqrt(1.0 - t1*t1) + s * t2;

    double Nh[3];
    for (int i = 0; i < 3; i++)
        Nh[i] = t1*T1[i] + t2*T2[i] + sqrt(fmax(0.0, 1.0 - t1*t1 - t2*t2))*Vh[i];

    /* unstretch back to slope space, then renormalize the scaled vector */
    double sx = alpha * Nh[0], sy = alpha * Nh[1], sz = fmax(0.0, Nh[2]);
    double nlen = sqrt(sx*sx + sy*sy + sz*sz);
    H[0] = sx / nlen;
    H[1] = sy / nlen;
    H[2] = sz / nlen;
}

The stretch by α, the projection onto the view hemisphere, and the unstretch back to slope space are the three structural steps. Omitting the stretch produces samples from the full normal distribution, which is the biased variant mentioned above; omitting the hemisphere projection produces samples that sometimes reflect into the surface.

The density in outgoing-direction space is not the same as the density in half-vector space. The Jacobian is 1/(4 |ω_i · h|), so the pdf of the outgoing direction is p_h(h) / (4 |ω_i · h|). A shader that returns the half-vector density as the direction density is off by that factor, and the error is largest at grazing angles, exactly where the specular lobe is most visible.

Selecting among lobes

A layered material has several components: a diffuse base, a rough specular coat, maybe a sheen term. Each has its own sampling routine and its own density. The estimator must use one density overall, so the components are chosen by a discrete probability and the returned pdf is the mixture:

$$p(\omega) = \sum_k w_k, p_k(\omega), \qquad \sum_k w_k = 1$$

with w_k the probability of choosing component k. Any set of weights gives an unbiased estimator. The choice that minimizes variance is roughly proportional to each component’s total energy, the integral of its contribution, because a component contributing little light should get few samples.

The trap is using fixed weights across all materials. A metal with no diffuse component should assign the diffuse lobe a probability of zero; otherwise a share of the samples lands where the integrand is zero and the estimator averages over contributions that never arrive, raising variance for nothing. A common compromise computes the weights from the material’s albedo values at shading time, which costs a few operations and adapts to the material.

A second trap is the density bookkeeping. When the mixture is used, the denominator must be the mixture density, not the selected component’s density. Using the component density alone produces an estimator that is correct in expectation only if the selection probability is folded into the weight, which is equivalent, but only if done consistently. Half-applying the correction is a frequent source of a lobe that renders at the wrong brightness only where a certain component dominates.

When the sampled density is not the BRDF’s

The Fresnel factor makes the true integrand narrower or wider than the distribution alone, depending on the angle. At grazing incidence, Fresnel reflectance rises steeply with the angle to the half-vector, so the integrand favors half-vectors closer to normal incidence than D does. Sampling D alone over-weights the grazing half-vectors, and the weights become noisy.

The magnitude of this effect depends on the material. For a dielectric with F₀ near 0.04 it is negligible away from the silhouette. For a conductor with a strongly colored F it changes the effective lobe shape enough to matter for low sample counts. Sampling a distribution that includes the Fresnel factor has no closed form in general, so practical renderers either accept the mismatch or apply a roughness adjustment that approximates the narrowing.

The masking-shadowing term presents the mirror-image problem: at grazing angles it suppresses the lobe, so the true integrand is again not proportional to D. Since G depends on the outgoing direction as well as the incoming one, it cannot be folded into the half-vector density at all; the sample is drawn first, then G is evaluated. This is one reason why very rough surfaces at grazing angles remain noisy even with VNDF sampling.

Measuring whether sampling works

Two diagnostics distinguish a sampling bug from a variance problem. The first is a white furnace test: put the surface in a uniform environment of radiance one and check that the mean estimate converges to the directional albedo, and that the weights for individual samples cluster around that value. A mean that converges but with a long weight tail means the density has the right mass in the wrong places; a mean that converges to the wrong value means the pdf is misnormalized.

The second is to evaluate the same integral with two independent sampling strategies and compare. When the two estimates disagree beyond their stated confidence intervals, one of them has a normalization error. This is cheap to automate and catches Jacobian mistakes that are invisible in rendered images, where the wrong pdf often manifests as a slight brightness offset that the eye attributes to material tuning.