Shading and Monte Carlo light transport

A ray tracer that can find the closest hit (previous skill) still has to answer the question that makes images: how much light leaves this point towards the camera? The answer is an integral over every direction light can arrive from, and almost every production renderer — offline path tracers, real-time GI, and the "ground truth" renders that radiance fields and Gaussian splats are trained on — evaluates it the same way: by Monte Carlo sampling.

This lesson builds the kernel of a path tracer: the units of light, the rendering equation, BRDFs (Lambert, Blinn–Phong, GGX), Monte Carlo estimators and importance sampling, a local frame around a normal, Russian roulette, and the linear-vs-sRGB distinction that decides whether the final pixels are right.

1. Radiometry in five definitions

Radiance is the one to internalise: it is constant along a ray in empty space, it is what a camera measures, and it is what NeRF's MLP outputs. Irradiance gathers radiance from the hemisphere, weighting each direction by cosθ (Lambert's cosine law — light arriving at a slant spreads over more area):

E(x) = ∫_Ω L(x, ω) cosθ dω, with dω = sinθ dθ dφ.

For a uniform sky of radiance 1, E = ∫₀^{2π}∫₀^{π/2} cosθ sinθ dθ dφ = π. That π will keep showing up.

2. The rendering equation

Kajiya (1986) wrote down what every renderer approximates. Outgoing radiance is emitted radiance plus reflected radiance:

Lₒ(x, ωₒ) = Lₑ(x, ωₒ) + ∫_Ω f(x, ωᵢ, ωₒ) · Lᵢ(x, ωᵢ) · cosθᵢ dωᵢ

and Lᵢ is itself the Lₒ of whatever point the ray from x in direction ωᵢ hits — so the equation is recursive. A path tracer unrolls the recursion one random direction per bounce:

3. BRDFs

A BRDF must be non-negative, reciprocal (f(ωᵢ, ωₒ) = f(ωₒ, ωᵢ)) and energy conserving: the directional albedo ∫ f cosθᵢ dωᵢ ≤ 1 for every ωₒ — a surface can't reflect more than arrives.

Lambert. A perfectly matte surface reflects equally in all directions: f = albedo / π. The π is not decoration: ∫ (albedo/π) cosθ dω = albedo. Leave it out and a white wall under a sky of radiance 1 glows with radiance π.

Blinn–Phong (1977) is the classic real-time highlight: f = (n + 8)/(8π) · (n·h)ⁿ with the half vector h = normalize(l + v). The (n + 8)/(8π) factor is the approximate normalisation that keeps the highlight's energy constant as the exponent sharpens it; without it, shiny materials go dark.

Microfacet models treat the surface as tiny perfect mirrors whose normals follow a distribution D. Only facets with normal h reflect l into v; the Cook–Torrance form is

f(l, v) = D(h) · G(l, v) · F(v·h) / (4 (n·l)(n·v))

• D — GGX / Trowbridge–Reitz (Walter et al. 2007): D(h) = α² / (π ((n·h)²(α² − 1) + 1)²). It is a distribution: projected onto the macro surface it integrates to one, ∫ D(h) (n·h) dω_h = 1. Its long tail is why GGX highlights have the soft glow real materials show.
• G — Smith: G = G₁(l) G₁(v), G₁(v) = 2(n·v) / (n·v + √(α² + (1 − α²)(n·v)²)). 1 at normal incidence, falling to 0 at grazing angles, where facets hide each other.
• F — Schlick (1994): F = F₀ + (1 − F₀)(1 − cosθ)⁵. F₀ ≈ 0.04 for dielectrics, the tinted base colour for metals. Everything reflects like a mirror at grazing angles.

The 4 (n·l)(n·v) converts from the half-vector domain to the light-direction domain (dω_h = dω_l / (4 (v·h))) and from projected to true facet area.

#include <cmath>
#include <numbers>
#include <print>

const double pi = std::numbers::pi;

double ggx_d(double nh, double a) {
    double d = nh * nh * (a * a - 1) + 1;
    return nh > 0 ? a * a / (pi * d * d) : 0.0;
}

int main() {
    // Check the normalisation  2π ∫₀¹ D(μ) μ dμ = 1  numerically (μ = n·h).
    for (double alpha : {0.05, 0.2, 0.5, 1.0}) {
        const int steps = 200000;
        double sum = 0;
        for (int i = 0; i < steps; ++i) {
            double mu = (i + 0.5) / steps;
            sum += ggx_d(mu, alpha) * mu;
        }
        std::println("alpha {:4}: peak D = {:9.2f}   integral = {:.4f}",
                     alpha, ggx_d(1.0, alpha), 2 * pi * sum / steps);
    }
}

Small α means a tall narrow peak (a near-mirror); the area stays one. That is the property the hidden tests check, and the first thing to verify when you port a BRDF from a paper.

4. Monte Carlo integration

To estimate I = ∫ g(x) dx, draw xᵢ from any density p that is non-zero where g is, and average:

Î = (1/N) Σ g(xᵢ) / p(xᵢ), E[Î] = ∫ (g/p) p dx = I.

The estimator is unbiased, its variance is Var[g/p] / N, so the error (standard deviation) falls as 1/√N — in any dimension. Halving the noise costs 4× the samples. The only other lever is the variance per sample: make p look like g (importance sampling). If p ∝ g exactly, every sample is I and the variance is zero.

#include <cmath>
#include <cstdint>
#include <numbers>
#include <print>

struct Rng {  // xorshift32: deterministic everywhere
    std::uint32_t s = 2463534242u;
    float next() {
        s ^= s << 13; s ^= s >> 17; s ^= s << 5;
        return static_cast<float>(s >> 8) * 0x1p-24f;
    }
};

int main() {
    const double pi = std::numbers::pi;
    // I = ∫₀^π sin x dx = 2. RMS error over 64 independent runs per N.
    Rng rng;
    for (int n : {16, 256, 4096, 65536}) {
        double se_uni = 0, se_imp = 0;
        for (int run = 0; run < 64; ++run) {
            double uni = 0, imp = 0;
            for (int i = 0; i < n; ++i) {
                double x = pi * rng.next();               // p(x) = 1/π
                uni += std::sin(x) / (1 / pi);
                double y = std::acos(1 - 2 * rng.next()); // p(y) = sin(y)/2  (inverse CDF)
                imp += std::sin(y) / (std::sin(y) / 2);   // exactly 2 every time
            }
            se_uni += (uni / n - 2) * (uni / n - 2);
            se_imp += (imp / n - 2) * (imp / n - 2);
        }
        std::println("N = {:6}: uniform RMS error {:.5f}   importance RMS error {:.5f}",
                     n, std::sqrt(se_uni / 64), std::sqrt(se_imp / 64));
    }
}

The uniform column's error shrinks by about 4× for every 16× more samples (√16 = 4); the importance-sampled column is exact from the first sample because p ∝ sin (up to float rounding in acos). (The acos(1 − 2u) is inverse-CDF sampling: if P(y) is the cumulative distribution, P⁻¹(u) with uniform u is distributed by p.)

5. Sampling the hemisphere

Work in a local frame where the normal is +z, so cosθ = z.

Uniform: p(ω) = 1/(2π). Since the area of a sphere band is proportional to its height (Archimedes), z = u₁ uniform and φ = 2π u₂ give uniform directions: (√(1 − z²) cos φ, √(1 − z²) sin φ, z).

Cosine-weighted: p(ω) = cosθ / π — it cancels the cosine in the integrand. Malley's method: pick a uniform point on the unit disk (r = √u₁, φ = 2π u₂; the √ because area grows as r²) and lift it onto the hemisphere, z = √(1 − r²). The disk's density 1/π per unit area, divided by the projection factor dA = cosθ dω, is exactly cosθ/π per steradian.

With cosine samples, a Lambertian sample is (albedo/π) · L · cosθ / (cosθ/π) = albedo · L. Under a uniform sky that is a constant — the estimator has zero variance. This is the white furnace test: put a white (albedo 1) object in a uniform white environment; with no absorption it must disappear (every pixel exactly 1). Any missing π, cosθ or mismatched pdf shows up immediately:

#include <cmath>
#include <cstdint>
#include <numbers>
#include <print>

const float pi = std::numbers::pi_v<float>;

struct Rng {
    std::uint32_t s = 0x9E3779B9u;
    float next() {
        s ^= s << 13; s ^= s >> 17; s ^= s << 5;
        return static_cast<float>(s >> 8) * 0x1p-24f;
    }
};

int main() {
    Rng rng;
    const int n = 100000;
    double good = 0, wrong_pdf = 0, no_cos = 0, no_pi = 0, uniform = 0;
    for (int i = 0; i < n; ++i) {
        float u1 = rng.next(), u2 = rng.next();
        float z = std::sqrt(1 - u1);                   // Malley: cos θ
        float pdf = z / pi;
        float f = 1.0f / pi, L = 1.0f;
        good      += f * L * z / pdf;
        wrong_pdf += f * L * z / (1 / (2 * pi));       // uniform pdf, cosine samples
        no_cos    += f * L / pdf;
        no_pi     += 1.0f * L * z / pdf;
        float zu = u2;                                  // a uniform sample instead
        uniform   += f * L * zu / (1 / (2 * pi));
    }
    std::println("correct {:.4f}  wrong pdf {:.4f}  no cos {:.4f}  no 1/pi {:.4f}  uniform {:.4f}",
                 good / n, wrong_pdf / n, no_cos / n, no_pi / n, uniform / n);
}

Expected: 1, 4/3, 2, π, and ≈ 1 (uniform sampling is unbiased, just noisy).

6. A frame around the normal

Samples are generated around +z and must be rotated into world space around the surface normal n: find unit t, b with (t, b, n) orthonormal, then world = x·t + y·b + z·n. The classic recipe "cross n with whichever axis is least parallel" branches; Frisvad (2012) gave a branch-free formula with a = 1/(1 + n.z) — which is infinite at n = −z and inaccurate near it. Duff et al. (2017) fixed it with one sign:

#include <cmath>
#include <print>

struct V { float x, y, z; };
float dot(V a, V b) { return a.x * b.x + a.y * b.y + a.z * b.z; }

void frisvad(V n, V& t, V& b) {
    float a = 1.0f / (1.0f + n.z);
    float c = -n.x * n.y * a;
    t = {1.0f - n.x * n.x * a, c, -n.x};
    b = {c, 1.0f - n.y * n.y * a, -n.y};
}

void duff(V n, V& t, V& b) {
    float sign = std::copysign(1.0f, n.z);
    float a = -1.0f / (sign + n.z);
    float c = n.x * n.y * a;
    t = {1.0f + sign * n.x * n.x * a, sign * c, -sign * n.x};
    b = {c, sign + n.y * n.y * a, -n.y};
}

int main() {
    for (float eps : {1e-1f, 1e-3f, 1e-4f}) {
        float len = std::sqrt(1.25f * eps * eps + 1);
        V n{eps / len, 0.5f * eps / len, -1 / len};   // almost straight down
        V t, b;
        frisvad(n, t, b);
        float e1 = std::fabs(dot(t, n)) + std::fabs(dot(t, b)) + std::fabs(dot(t, t) - 1);
        duff(n, t, b);
        float e2 = std::fabs(dot(t, n)) + std::fabs(dot(t, b)) + std::fabs(dot(t, t) - 1);
        std::println("eps {:g}: Frisvad error {:.2e}   Duff error {:.2e}", eps, e1, e2);
    }
}

Duff's sign + n.z is never smaller than 1 in magnitude, so there is no cancellation and no special case. (PBRT-v4's CoordinateSystem and most production renderers use it.) Normals near −z are not exotic: every floor seen from below, every inward-facing sphere and every back face hits them.

7. Russian roulette and the path tracer

A path could bounce forever; truncating at a fixed depth loses energy (bias). Russian roulette (Arvo & Kirk 1990) terminates paths at random: continue with probability q, and divide the survivor's throughput by q. The expected contribution is q · (C/q) + (1 − q) · 0 = C — unbiased, and dim paths (small throughput, so small q) rarely survive. Here it is in a complete, tiny path tracer: a camera inside a sphere whose walls emit Lₑ = 1 and reflect with albedo a. Light bounces forever, so the answer is the geometric series Σ aᵏ = 1/(1 − a) — a closed-form check for every piece above.

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <numbers>
#include <print>

const float pi = std::numbers::pi_v<float>;
struct V { float x, y, z; };
V operator*(float s, V v) { return {s * v.x, s * v.y, s * v.z}; }
V operator+(V a, V b) { return {a.x + b.x, a.y + b.y, a.z + b.z}; }
float dot(V a, V b) { return a.x * b.x + a.y * b.y + a.z * b.z; }

struct Rng {
    std::uint32_t s = 0x9E3779B9u;
    float next() { s ^= s << 13; s ^= s >> 17; s ^= s << 5; return float(s >> 8) * 0x1p-24f; }
};

// Leaving the unit sphere from inside: o + t d with |o| <= 1, d unit.
float exit_t(V o, V d) {
    float b = dot(o, d), c = dot(o, o) - 1;
    return -b + std::sqrt(std::max(0.0f, b * b - c));
}

int main() {
    const float albedo = 0.5f, Le = 1.0f;
    Rng rng;
    double total = 0;
    const int paths = 20000;
    for (int p = 0; p < paths; ++p) {
        V o{0, 0, 0}, d{0, 0, 1};
        float beta = 1, L = 0;
        for (int bounce = 0;; ++bounce) {
            V x = o + exit_t(o, d) * d;                 // hit the wall
            V n = -1.0f * x;                            // inward normal
            L += beta * Le;                             // emission
            // frame around n (Duff et al.), cosine-weighted direction
            float sign = std::copysign(1.0f, n.z), a = -1 / (sign + n.z), b = n.x * n.y * a;
            V t{1 + sign * n.x * n.x * a, sign * b, -sign * n.x}, bt{b, sign + n.y * n.y * a, -n.y};
            float u1 = rng.next(), u2 = rng.next(), r = std::sqrt(u1), phi = 2 * pi * u2;
            float z = std::sqrt(1 - u1);
            d = (r * std::cos(phi)) * t + (r * std::sin(phi)) * bt + z * n;
            beta *= (albedo / pi) * z / (z / pi);       // f cosθ / pdf = albedo
            if (bounce >= 3) {                          // Russian roulette
                float q = std::min(0.95f, beta);
                if (rng.next() >= q) break;
                beta /= q;
            }
            o = 0.9999f * x;                            // step off the surface
        }
        total += L;
    }
    std::println("estimate {:.3f}, exact {:.3f}", total / paths, Le / (1 - albedo));
}

Note o = 0.9999f * x: the next ray starts a hair inside the wall, the self-intersection ("shadow acne") guard from the ray-intersection skill.

8. Linear light, sRGB pixels

Radiance adds: two lights give the sum of their contributions, a pixel covered half by white and half by black receives half the light. All of that holds for linear values. Displays and 8-bit images instead store sRGB-encoded values, a curve close to x^(1/2.2) with a linear toe:

encode(x) = 12.92 x for x ≤ 0.0031308, else 1.055 x^(1/2.4) − 0.055.

So the pipeline is: decode sRGB inputs (textures, colour pickers) to linear → light, accumulate samples, blend, filter in linear → tone-map (compress HDR radiance into [0, 1], e.g. Reinhard x/(1 + x) or ACES) → encode once → quantise to 8 bits. Encoding before averaging is the classic bug: an anti-aliased edge between black and white comes out as 128 (which the display shows as 21 % of the light) instead of 188.

#include <cmath>
#include <print>

float encode(float x) { return x <= 0.0031308f ? 12.92f * x : 1.055f * std::pow(x, 1 / 2.4f) - 0.055f; }
float decode(float y) { return y <= 0.04045f ? y / 12.92f : std::pow((y + 0.055f) / 1.055f, 2.4f); }

int main() {
    float samples[] = {0, 0, 0, 1};            // 25 % coverage of a white object
    float lin = 0, enc = 0;
    for (float s : samples) { lin += s / 4; enc += encode(s) / 4; }
    std::println("average then encode: {}", std::lround(encode(lin) * 255));
    std::println("encode then average: {}  (the display shows only {:.0f}% of the light)",
                 std::lround(enc * 255), 100 * decode(enc));
}

Training data for neural rendering raises the same question: NeRF and 3D Gaussian Splatting fit colours in whatever space the photos are in (usually sRGB-encoded LDR); RawNeRF (Mildenhall et al. 2022) trains on linear raw images instead to get correct HDR and blur. Know which space your loss sees.

In Python (NumPy)

u1, u2 = rng.random(N), rng.random(N)
r, phi = np.sqrt(u1), 2 * np.pi * u2
d = np.stack([r * np.cos(phi), r * np.sin(phi), np.sqrt(1 - u1)], -1)   # cosine samples
E = np.mean(sky(d @ frame) * d[:, 2] / (d[:, 2] / np.pi))               # irradiance
srgb = np.where(x <= 0.0031308, 12.92 * x, 1.055 * x ** (1 / 2.4) - 0.055)

In an interview

• "Why albedo/π?" — ∫ cosθ dω = π; it makes a white Lambertian surface reflect exactly what it receives.
• "How do you know your path tracer is right?" — white furnace test (objects vanish), analytic cases like the emissive sphere above, convergence at 1/√N, reciprocity and energy checks on each BRDF.
• "Variance reduction?" — importance sampling (BRDF, light, both with MIS — Veach 1997), stratification / low-discrepancy sequences, Russian roulette, next-event estimation, and finally a denoiser.

Recap

• Radiance is what rays carry; irradiance integrates radiance × cosθ.
• The rendering equation: Lₒ = Lₑ + ∫ f · Lᵢ · cosθ dω; path tracing unrolls it.
• Lambert = albedo/π; GGX = D·G·F/(4 (n·l)(n·v)); BRDFs are reciprocal and conserve energy.
• Monte Carlo: average g/p, unbiased iff p is the density you sampled; error ∝ 1/√N.
• Cosine sampling (Malley) has pdf cosθ/π and makes Lambert samples constant.
• Build frames with Duff et al.'s branchless ONB; terminate with Russian roulette and divide by q.
• Light and average in linear space; tone-map and sRGB-encode once at the end.