6. Importance-Sampled Materials

Ray Tracing: The Rest of Your Life (v3.2.3): 6 Importance Sampling Materials / 3.6 Importance-Sampled Materials

Section 3.5 introduced the scattering PDF $s(\omega )=\cos \theta /\pi$ for a Lambertian surface. I now clarify the relationship between the sampling PDF $p$ and scattering PDF $s$, then use a numerical experiment to examine how their differences affect convergence and the final integral value.

This chapter covers two topics:

• Comparing three cases: $p=s$ (matched), $p\ne s$ (mismatched), and changing $s$ itself (a different material)
• A numerical experiment using the test function $f(\omega )={\cos }^2\theta$

The scattering_pdf and pdf_value Framework

A Monte Carlo estimate of the rendering equation can be written as

$$\text{Color}\approx A\cdot \frac{s(\omega )}{p(\omega )}\cdot \text{color}(\omega )$$

where

• $s(\omega )$ is the material's scattering PDF (scattering_pdf).
• $p(\omega )$ is the sampling PDF used to generate a direction (pdf_value).

Once $s$ is fixed, the integral value—the physical color—is fixed. The choice of $p$ does not change that value; it affects only the convergence rate, or variance.

Three Combinations

Sections 6.1 through 6.3 of the original compare the following three cases.

Method	Scattering PDF $s$	Sampling PDF $p$	Integral	Convergence
A	$\cos \theta /\pi$	$\cos \theta /\pi$ (perfect match)	$1/2$	Fast (variance near zero)
B	$\cos \theta /\pi$	$1/(2\pi )$ (mismatch)	$1/2$	Slow
C	$1/(2\pi )$	$1/(2\pi )$	$1/3$	Fast (but a different material)

Methods A and B have the same Lambertian $s$, so their integral values are equal and the choice of $p$ affects only convergence. Method C changes $s$ itself, creating a fundamentally different material with uniform-hemisphere scattering and therefore a different integral value.

The original illustrates the difference among these three cases with rendered Cornell-box images. This chapter confines itself to confirming the same conclusion numerically—Methods A and B yield $1/2$, while Method C yields $1/3$—and does not implement a live Cornell-box renderer. I first introduce that renderer in Section 3.9, “Directly Sampling Light,” and improve it incrementally in subsequent chapters.

Numerical Experiment: A Weighted Integral of ${\cos }^2\theta$

To make the discussion concrete, I evaluate the weighted integral

$$I_s={\int }_{\text{hemisphere}}f(\omega )\cdot s(\omega )d\omega$$

for the test function $f(\omega )={\cos }^2\theta$.

For Lambertian scattering:

$$I_s^{\text{Lambert}}=\frac{1}{\pi }{\int }_{\text{hemisphere}}{\cos }^3\theta d\omega =\frac{1}{\pi }{\int }_0^{2\pi }{\int }_0^{\pi /2}{\cos }^3\theta \sin \theta d\theta d\varphi =2{\int }_0^1{u}^{3}du=\frac{1}{2}$$

For uniform-hemisphere scattering:

$$I_s^{\text{uniform}}=\frac{1}{2\pi }{\int }_{\text{hemisphere}}{\cos }^2\theta d\omega =\frac{1}{2\pi }{\int }_0^{2\pi }{\int }_0^{\pi /2}{\cos }^2\theta \sin \theta d\theta d\varphi ={\int }_0^1{u}^{2}du=\frac{1}{3}$$

Rust Implementation

The r306-importance-sampling-play crate runs the three estimators with independent RNGs.

sample_cosine uses the formula $\cos \theta =\sqrt{1-r_2}$ from Malley's method, where $r_2$ is uniformly random.

r306-importance-sampling-play/src/lib.rs

// Sample a direction from the cosine-weighted hemisphere (z = cosθ).
fn sample_cosine(rng: &mut XorShift64) -> f64 {
    let r2 = rng.next_f64();
    rng.next_f64(); // consume r1 (only cosθ = sqrt(1-r2) is needed)
    (1.0 - r2).sqrt()
}

// Sample a direction from the uniform hemisphere (z = cosθ = 1 - r2).
fn sample_uniform(rng: &mut XorShift64) -> f64 {
    let r2 = rng.next_f64();
    rng.next_f64(); // consume r1
    1.0 - r2
}

Both functions consume the RNG twice because a direction is actually a two-dimensional sample in $\varphi$ and $\theta$. This demonstration does not need the $\varphi$ component, r1, but still consumes it.

Each of the three estimators calculates $f(\omega )\cdot s(\omega )/p(\omega )$. In Method A, $s/p=1$, so the contribution is simply $f(\omega )={\cos }^2\theta$.

r306-importance-sampling-play/src/lib.rs

// Method A: Lambertian scatter + cosine sampling (perfect match).
// contribution = f(d) · scatter_pdf(d) / sampling_pdf(d) = cos²θ
fn estimate_method_a(n: usize, rng: &mut XorShift64) -> f64 {
    let sum: f64 = (0..n)
        .map(|_| {
            let cos_theta = sample_cosine(rng);
            let s = scatter_pdf_lambertian(cos_theta);
            let p = pdf_cosine(cos_theta);
            f(cos_theta) * s / p
        })
        .sum();
    sum / n as f64
}

// Method B: Lambertian scatter + uniform sampling (imperfect match).
// contribution = f(d) · scatter_pdf_lambertian(d) / sampling_pdf_uniform(d)
//             = cos²θ · (cosθ/π) / (1/(2π)) = 2cos³θ
fn estimate_method_b(n: usize, rng: &mut XorShift64) -> f64 {
    let sum: f64 = (0..n)
        .map(|_| {
            let cos_theta = sample_uniform(rng);
            let s = scatter_pdf_lambertian(cos_theta);
            let p = pdf_uniform(cos_theta);
            f(cos_theta) * s / p
        })
        .sum();
    sum / n as f64
}

// Method C: Uniform hemisphere scatter + uniform sampling (different material).
// contribution = f(d) · scatter_pdf_uniform(d) / sampling_pdf_uniform(d) = cos²θ
// Expected value = 1/3 (different from Lambertian → different render result)
fn estimate_method_c(n: usize, rng: &mut XorShift64) -> f64 {
    let sum: f64 = (0..n)
        .map(|_| {
            let cos_theta = sample_uniform(rng);
            let s = scatter_pdf_uniform_hemi(cos_theta);
            let p = pdf_uniform(cos_theta);
            f(cos_theta) * s / p
        })
        .sum();
    sum / n as f64
}

After a single comparison with $N={10}^6$, generate_report displays Methods A and B side by side at convergence checkpoints from $N=100$ through $500000$.

Differences Between C++ and Rust

The C++ version compares the PDFs inside ray_color() as follows:

double scattering_pdf = rec.mat->scattering_pdf(r, rec, scattered);
double pdf_value = scattering_pdf; // or a different PDF
color color_from_scatter =
    (attenuation * scattering_pdf * ray_color(scattered, depth-1, world))
    / pdf_value;

The Rust estimate_method_a/b/c functions implement the same $s/p$ calculation as a numerical demonstration.

r306-importance-sampling-play/src/lib.rs

use std::fmt::Write;

const PI: f64 = std::f64::consts::PI;

// Exact values of the integrals:
//   ∫_hemi f(d) · scatter_pdf(d) dω  where  f(d) = cos²θ
// Lambertian scatter_pdf = cosθ/π → integral = 1/2
// Uniform   scatter_pdf = 1/(2π) → integral = 1/3
const LAMBERTIAN_INTEGRAL: f64 = 0.5;
const UNIFORM_INTEGRAL: f64 = 1.0 / 3.0;

struct XorShift64 {
    state: u64,
}

impl XorShift64 {
    fn new(seed: u64) -> Self {
        let state = if seed == 0 { 0x9e3779b97f4a7c15 } else { seed };
        Self { state }
    }

    fn next_u64(&mut self) -> u64 {
        let mut x = self.state;
        x ^= x << 13;
        x ^= x >> 7;
        x ^= x << 17;
        self.state = x;
        x
    }

    fn next_f64(&mut self) -> f64 {
        (self.next_u64() >> 11) as f64 * (1.0 / (1u64 << 53) as f64)
    }
}

// Sample a direction from the cosine-weighted hemisphere (z = cosθ).
fn sample_cosine(rng: &mut XorShift64) -> f64 {
    let r2 = rng.next_f64();
    rng.next_f64(); // consume r1 (only cosθ = sqrt(1-r2) is needed)
    (1.0 - r2).sqrt()
}

// Sample a direction from the uniform hemisphere (z = cosθ = 1 - r2).
fn sample_uniform(rng: &mut XorShift64) -> f64 {
    let r2 = rng.next_f64();
    rng.next_f64(); // consume r1
    1.0 - r2
}

fn pdf_cosine(cos_theta: f64) -> f64 {
    cos_theta / PI
}

fn pdf_uniform(_cos_theta: f64) -> f64 {
    1.0 / (2.0 * PI)
}

fn scatter_pdf_lambertian(cos_theta: f64) -> f64 {
    cos_theta / PI
}

fn scatter_pdf_uniform_hemi(_cos_theta: f64) -> f64 {
    1.0 / (2.0 * PI)
}

// f(d) = cos²θ  (arbitrary smooth function for demonstration)
fn f(cos_theta: f64) -> f64 {
    cos_theta * cos_theta
}

// Method A: Lambertian scatter + cosine sampling (perfect match).
// contribution = f(d) · scatter_pdf(d) / sampling_pdf(d) = cos²θ
fn estimate_method_a(n: usize, rng: &mut XorShift64) -> f64 {
    let sum: f64 = (0..n)
        .map(|_| {
            let cos_theta = sample_cosine(rng);
            let s = scatter_pdf_lambertian(cos_theta);
            let p = pdf_cosine(cos_theta);
            f(cos_theta) * s / p
        })
        .sum();
    sum / n as f64
}

// Method B: Lambertian scatter + uniform sampling (imperfect match).
// contribution = f(d) · scatter_pdf_lambertian(d) / sampling_pdf_uniform(d)
//             = cos²θ · (cosθ/π) / (1/(2π)) = 2cos³θ
fn estimate_method_b(n: usize, rng: &mut XorShift64) -> f64 {
    let sum: f64 = (0..n)
        .map(|_| {
            let cos_theta = sample_uniform(rng);
            let s = scatter_pdf_lambertian(cos_theta);
            let p = pdf_uniform(cos_theta);
            f(cos_theta) * s / p
        })
        .sum();
    sum / n as f64
}

// Method C: Uniform hemisphere scatter + uniform sampling (different material).
// contribution = f(d) · scatter_pdf_uniform(d) / sampling_pdf_uniform(d) = cos²θ
// Expected value = 1/3 (different from Lambertian → different render result)
fn estimate_method_c(n: usize, rng: &mut XorShift64) -> f64 {
    let sum: f64 = (0..n)
        .map(|_| {
            let cos_theta = sample_uniform(rng);
            let s = scatter_pdf_uniform_hemi(cos_theta);
            let p = pdf_uniform(cos_theta);
            f(cos_theta) * s / p
        })
        .sum();
    sum / n as f64
}

pub fn generate_report() -> String {
    let mut out = String::new();

    // Header
    writeln!(out, "Integral: ∫_hemi cos²θ · scatter_pdf(d) dω").unwrap();
    writeln!(out).unwrap();

    // Single-run comparison with N = 1_000_000
    let n = 1_000_000;
    let mut rng_a = XorShift64::new(0x1a2b_3c4d_0001);
    let mut rng_b = XorShift64::new(0x1a2b_3c4d_0002);
    let mut rng_c = XorShift64::new(0x1a2b_3c4d_0003);

    let result_a = estimate_method_a(n, &mut rng_a);
    let result_b = estimate_method_b(n, &mut rng_b);
    let result_c = estimate_method_c(n, &mut rng_c);

    writeln!(out, "=== Single Run (N = {n}) ===").unwrap();
    writeln!(out, "{:<10} {:<12} {:<12} {:<12}", "Method", "Estimate", "Exact", "Error").unwrap();
    writeln!(
        out,
        "{:<10} {:<12.8} {:<12.8} {:<+12.2e}",
        "A (cos/cos)",
        result_a,
        LAMBERTIAN_INTEGRAL,
        result_a - LAMBERTIAN_INTEGRAL
    )
    .unwrap();
    writeln!(
        out,
        "{:<10} {:<12.8} {:<12.8} {:<+12.2e}",
        "B (cos/uni)",
        result_b,
        LAMBERTIAN_INTEGRAL,
        result_b - LAMBERTIAN_INTEGRAL
    )
    .unwrap();
    writeln!(
        out,
        "{:<10} {:<12.8} {:<12.8} {:<+12.2e}",
        "C (uni/uni)",
        result_c,
        UNIFORM_INTEGRAL,
        result_c - UNIFORM_INTEGRAL
    )
    .unwrap();
    writeln!(out).unwrap();

    // Convergence table for A and B (same target, different PDFs)
    writeln!(out, "=== Convergence: Method A (cos/cos) vs B (cos/uni) ===").unwrap();
    writeln!(
        out,
        "{:<10} {:<14} {:<12} {:<14} {:<12}",
        "N", "A Estimate", "A Error", "B Estimate", "B Error"
    )
    .unwrap();

    let checkpoints = [100usize, 1_000, 10_000, 100_000, 500_000];
    let mut rng_a = XorShift64::new(0x5e6f_7a8b_0001);
    let mut rng_b = XorShift64::new(0x5e6f_7a8b_0002);

    for &n in &checkpoints {
        let a = estimate_method_a(n, &mut rng_a);
        let b = estimate_method_b(n, &mut rng_b);
        writeln!(
            out,
            "{:<10} {:<14.8} {:<12.2e} {:<14.8} {:<12.2e}",
            n,
            a,
            a - LAMBERTIAN_INTEGRAL,
            b,
            b - LAMBERTIAN_INTEGRAL
        )
        .unwrap();
    }
    writeln!(out).unwrap();

    // Method C comparison
    writeln!(out, "=== Note: Method C gives a different answer ===").unwrap();
    writeln!(
        out,
        "Lambertian scatter_pdf (A & B): exact = {}",
        LAMBERTIAN_INTEGRAL
    )
    .unwrap();
    writeln!(
        out,
        "Uniform scatter_pdf   (C):      exact = {:.8}",
        UNIFORM_INTEGRAL
    )
    .unwrap();

    out
}

Results in the Browser

In the convergence table, Method A's A Error is consistently smaller than Method B's B Error. The Error column in the single-run comparison confirms that Methods A and B have the same exact value, 1/2, while only Method C differs at 1/3.

Summary

• In the rendering-equation estimate $A\cdot s(\omega )/p(\omega )\cdot \text{color}(\omega )$, $s$ determines the integral value and $p$ determines the convergence rate.
• With $p=s$ in Method A, variance is minimized. Method B, where $p\ne s$, still converges to the same integral value.
• Changing $s$ itself, as in Method C, creates a different material and changes the integral value from $1/2$ to $1/3$.

The next chapter connects this framework to the ray tracer itself.