5. Light Scattering

Ray Tracing: The Rest of Your Life (v3.2.3): 5 Light Scattering / 3.5 Light Scattering

Sections 3.3 and 3.4 developed the framework for Monte Carlo estimation and integration over a sphere. I now introduce a probabilistic model for light scattering at a surface, establishing the ideas that will reduce noise in the ray tracer.

This chapter covers three topics:

• Describing a scattering model with albedo and a scattering PDF
• The cosine density of a Lambertian surface and importance sampling based on it
• Numerical verification by estimating ${\int }_{\text{hemisphere}}\cos \theta d\omega =\pi$

Albedo and the Scattering PDF

Let $A$, the albedo, be the probability that light incident on a surface scatters. When scattering occurs, the directional distribution can be described by a PDF $s(\mathbf{d})$, which I call the scattering PDF.

The radiance of the reflected light is given by the following integral:

$$\text{Color}={\int }_{\text{hemisphere}}A\cdot s(\mathbf{d})\cdot \text{color}(\mathbf{d})d\omega$$

With a sampling PDF $p(\mathbf{d})$, a Monte Carlo estimator approximates it as

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

Setting $p=s$ makes $s/p=1$, reducing the expression to $\text{Color}\approx A\cdot \text{color}(\mathbf{d})$ with zero variance. This is importance sampling applied to scattering.

Lambertian Scattering and Cosine Density

The scattering PDF of a Lambertian, or diffusely reflecting, surface is

$$s(\mathbf{d})=\frac{\cos \theta }{\pi }$$

where $\theta$ is the angle from the normal. This is the cosine density. Its normalization follows from

$${\int }_{\text{hemisphere}}\frac{\cos \theta }{\pi }d\omega =\frac{1}{\pi }{\int }_0^{2\pi }{\int }_0^{\pi /2}\cos \theta \sin \theta d\theta d\varphi =\frac{2\pi }{\pi }\cdot \frac{1}{2}=1$$

The Exact Value of the Integral

The radiance integral for a Lambertian surface with albedo $A=1$ under uniform incident light is

$$I={\int }_{\text{hemisphere}}\cos \theta d\omega$$

Expanding it in spherical coordinates gives

$$I={\int }_0^{2\pi }{\int }_0^{\pi /2}\cos \theta \sin \theta d\theta d\varphi =2\pi {\left[\frac{{\sin }^2\theta }{2}\right]}_0^{\pi /2}=\pi$$

Two Methods for Sampling a Hemisphere

Uniform Hemisphere Sampling

Because the hemisphere has area $2\pi$, the uniform PDF is

$$p_{\text{uniform}}(\omega )=\frac{1}{2\pi }$$

Setting $z=u\in [0,1)$ and $\varphi =2\pi v$ covers the upper hemisphere uniformly, with $z=\cos \theta \ge 0$.

Cosine-Density Sampling (Malley's Method)

The marginal distribution of $p(\cos \theta )$ is $p(c)=2c$, where $c=\cos \theta$. Its CDF is $P(c)=c^2$, so inverse transform sampling gives

$$c=\cos \theta =\sqrt{u},u\sim U[0,1)$$

Each sample in the estimator is then

$$\frac{\cos \theta }{p(\cos \theta /\pi )}=\frac{\cos \theta }{\cos \theta /\pi }=\pi$$

which is constant, giving zero variance.

Rust Implementation

The r305-scattering-pdf crate implements XorShift64 locally and compares the two estimators.

r305-scattering-pdf/src/lib.rs

fn sample_uniform_hemisphere(rng: &mut XorShift64) -> (f64, f64, f64) {
    let z = rng.next_f64();
    let phi = 2.0 * PI * rng.next_f64();
    let r = (1.0 - z * z).max(0.0).sqrt();
    (r * phi.cos(), r * phi.sin(), z)
}

fn sample_cosine_hemisphere(rng: &mut XorShift64) -> (f64, f64, f64) {
    let cos_theta = rng.next_f64().sqrt();
    let sin_theta = (1.0 - cos_theta * cos_theta).sqrt();
    let phi = 2.0 * PI * rng.next_f64();
    (sin_theta * phi.cos(), sin_theta * phi.sin(), cos_theta)
}

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

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

The estimation function accumulates $\cos \theta /p(\omega )$ over $N$ samples and returns the mean. The if cos_theta > 0.0 condition guards against floating-point error making $z$ slightly negative.

r305-scattering-pdf/src/lib.rs

fn estimate_lambert_integral(samples: usize, use_cosine: bool, rng: &mut XorShift64) -> f64 {
    let mut sum = 0.0;
    for _ in 0..samples {
        let (_, _, z) = if use_cosine {
            sample_cosine_hemisphere(rng)
        } else {
            sample_uniform_hemisphere(rng)
        };
        let cos_theta = z;
        if cos_theta > 0.0 {
            let pdf = if use_cosine {
                cosine_pdf(cos_theta)
            } else {
                uniform_pdf(cos_theta)
            };
            sum += cos_theta / pdf;
        }
    }
    sum / samples as f64
}

After a single estimate with $N={10}^6$, generate_report displays the uniform and cosine errors side by side at checkpoints from $N=100$ through $500000$.

Differences Between C++ and Rust

Cosine-hemisphere sampling in the original C++ is known as Malley's method and can be implemented by projecting a uniform sample on the unit disk onto the hemisphere.

vec3 random_cosine_direction() {
    auto r1 = random_double();
    auto r2 = random_double();
    auto phi = 2*pi*r1;
    auto x = cos(phi)*sqrt(r2);
    auto y = sin(phi)*sqrt(r2);
    auto z = sqrt(1-r2);
    return vec3(x, y, z);
}

Here $r=\sqrt{r_2}$ produces a uniform point on the disk, while $z=\sqrt{1-r_2}=\sqrt{1-r^2}=\cos \theta$. The Rust implementation uses the equivalent inverse transform, cos_theta = rng.next_f64().sqrt().

r305-scattering-pdf/src/lib.rs

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

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)
    }
}

fn sample_uniform_hemisphere(rng: &mut XorShift64) -> (f64, f64, f64) {
    let z = rng.next_f64();
    let phi = 2.0 * PI * rng.next_f64();
    let r = (1.0 - z * z).max(0.0).sqrt();
    (r * phi.cos(), r * phi.sin(), z)
}

fn sample_cosine_hemisphere(rng: &mut XorShift64) -> (f64, f64, f64) {
    let cos_theta = rng.next_f64().sqrt();
    let sin_theta = (1.0 - cos_theta * cos_theta).sqrt();
    let phi = 2.0 * PI * rng.next_f64();
    (sin_theta * phi.cos(), sin_theta * phi.sin(), cos_theta)
}

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

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

fn estimate_lambert_integral(samples: usize, use_cosine: bool, rng: &mut XorShift64) -> f64 {
    let mut sum = 0.0;
    for _ in 0..samples {
        let (_, _, z) = if use_cosine {
            sample_cosine_hemisphere(rng)
        } else {
            sample_uniform_hemisphere(rng)
        };
        let cos_theta = z;
        if cos_theta > 0.0 {
            let pdf = if use_cosine {
                cosine_pdf(cos_theta)
            } else {
                uniform_pdf(cos_theta)
            };
            sum += cos_theta / pdf;
        }
    }
    sum / samples as f64
}

pub fn generate_report() -> String {
    let mut output = String::new();
    output.push_str("Lambertian Scattering PDF Importance Sampling\n");
    output.push_str("==========================================\n\n");
    output.push_str("Target integral: integral_{hemisphere} cos(theta) d_omega\n");
    output.push_str("(Lambertian surface radiance, albedo=1)\n");
    output.push_str(&format!("Exact value: {:.12}\n\n", EXACT_INTEGRAL));

    let mut uniform_rng = XorShift64::new(0x3050_0001u64);
    let mut cosine_rng = XorShift64::new(0x3050_0002u64);

    let uniform_est = estimate_lambert_integral(1_000_000, false, &mut uniform_rng);
    let cosine_est = estimate_lambert_integral(1_000_000, true, &mut cosine_rng);

    output.push_str("Single run (N=1,000,000):\n");
    output.push_str("method,estimate,abs_error\n");
    output.push_str(&format!(
        "uniform,{:.12},{:.12}\ncosine_pdf,{:.12},{:.12}\n\n",
        uniform_est,
        (uniform_est - EXACT_INTEGRAL).abs(),
        cosine_est,
        (cosine_est - EXACT_INTEGRAL).abs()
    ));

    output.push_str("Convergence checkpoints:\n");
    output.push_str("samples,uniform_error,cosine_error\n");

    let checkpoints = [100usize, 1_000, 10_000, 100_000, 500_000];
    for (index, n) in checkpoints.iter().copied().enumerate() {
        let seed_base = 0x3050_2000u64 + index as u64;
        let mut u_rng = XorShift64::new(seed_base ^ 0xaaaa);
        let mut c_rng = XorShift64::new(seed_base ^ 0xbbbb);

        let u_est = estimate_lambert_integral(n, false, &mut u_rng);
        let c_est = estimate_lambert_integral(n, true, &mut c_rng);

        output.push_str(&format!(
            "{},{:.12},{:.12}\n",
            n,
            (u_est - EXACT_INTEGRAL).abs(),
            (c_est - EXACT_INTEGRAL).abs()
        ));
    }

    output
}

Results in the Browser

Because every sample theoretically equals $\pi$, the cosine_error column should be nearly zero. Floating-point rounding leaves a small error in practice, but it is overwhelmingly smaller than uniform_error.

Summary

• I introduced a scattering model based on albedo $A$ and a scattering PDF $s(\mathbf{d})$. The fundamental quantity in its Monte Carlo estimator is $A\cdot s(\mathbf{d})/p(\mathbf{d})$.
• The scattering PDF of a Lambertian surface is the cosine density $s(\mathbf{d})=\cos \theta /\pi$, and the exact value of ${\int }_{\text{hemisphere}}\cos \theta d\omega$ is $\pi$.
• With $p=s$, cosine-PDF sampling makes the estimator constant and its variance zero, realizing ideal importance sampling.

The following chapters extend this idea to mixtures of multiple PDFs—light-source sampling and surface sampling—to build a ray tracer with less noise.