4. Monte Carlo Integration on the Sphere of Directions

Ray Tracing: The Rest of Your Life (v3.2.3): 4 MC Integration on the Sphere of Directions / 3.4 Monte Carlo Integration on the Sphere of Directions

Section 3.3 estimated an integral over a one-dimensional interval. I now integrate over the space of directions—the unit sphere $S^2$—which connects directly to ray tracing.

This chapter covers two topics:

• Uniform sampling on a sphere and its corresponding PDF
• Empirical verification using the spherical integral of ${\cos }^2\theta$

The Exact Value of the Integral

The integral of interest is

$$I={\int }_{S^2}{\cos }^2\theta d\omega$$

where $\theta$ is the angle from the $z$ axis and $d\omega$ is an element of solid angle. In spherical coordinates, $d\omega =\sin \theta d\theta d\varphi$, so

$$I={\int }_0^{2\pi }{\int }_0^{\pi }{\cos }^2\theta \sin \theta d\theta d\varphi =2\pi {\int }_0^{\pi }{\cos }^2\theta \sin \theta d\theta$$

Substituting $u=\cos \theta$ and $du=-\sin \theta d\theta$ gives

$$I=2\pi {\int }_{-1}^1{u}^{2}du=2\pi \cdot \frac{2}{3}=\frac{4\pi }{3}$$

Uniform Sampling on the Sphere

Generating a uniform distribution over the unit sphere uses the fact that $d\omega =\sin \theta d\theta d\varphi$ is its area element. Because this element depends on $\theta$, choosing both $\theta$ and $\varphi$ uniformly does not cover the sphere uniformly.

The correct approach is to choose $z=\cos \theta$ uniformly. A uniform distribution over $z\in [-1,1]$, combined with a uniform distribution over $\varphi \in [0,2\pi )$, covers the sphere uniformly. Given $u,v\sim U[0,1)$, set

$$z=1-2u,\varphi =2\pi v,r=\sqrt{1-z^2}$$

to construct the direction $(r\cos \varphi ,r\sin \varphi ,z)$. The PDF for this direction is the reciprocal of the sphere's surface area:

$$p(\omega )=\frac{1}{4\pi }$$

Rust Implementation

I add the r304-sphere-directions-mc crate for this chapter, with another local implementation of XorShift64.

The spherical sampling function and PDF constant are as follows.

r304-sphere-directions-mc/src/lib.rs

fn sample_unit_vector(rng: &mut XorShift64) -> (f64, f64, f64) {
    let z = 1.0 - 2.0 * 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 sphere_pdf() -> f64 {
    1.0 / (4.0 * PI)
}

The .max(0.0) guard prevents floating-point error from making $1-z^2$ slightly negative.

The estimator divides ${\cos }^2\theta =z^2$ by the PDF and averages the result over $N$ samples.

r304-sphere-directions-mc/src/lib.rs

fn estimate_integral(samples: usize, rng: &mut XorShift64) -> f64 {
    let pdf = sphere_pdf();
    let mut sum = 0.0;
    for _ in 0..samples {
        let (_, _, z) = sample_unit_vector(rng);
        let cosine_squared = z * z;
        sum += cosine_squared / pdf;
    }
    sum / samples as f64
}

generate_report constructs its output in two stages. It first reports a single estimate with $N={10}^6$, then estimates each checkpoint from $N=100$ through $500000$ using an independent RNG.

r304-sphere-directions-mc/src/lib.rs

    let checkpoints = [100usize, 1_000, 10_000, 100_000, 500_000];
    for (index, n) in checkpoints.iter().copied().enumerate() {
        let mut rng = XorShift64::new(0x3040_1000u64 + index as u64);
        let est = estimate_integral(n, &mut rng);
        output.push_str(&format!(
            "{},{:.12},{:.12}\n",
            n,
            est,
            (est - EXACT_INTEGRAL).abs()
        ));
    }

Because each checkpoint has an independent seed, estimates at small $N$ do not affect the results at larger $N$, unlike the cumulative-counting approach used in Section 3.2.

Differences Between C++ and Rust

The original C++ random_unit_vector has the same structure.

inline vec3 random_unit_vector() {
    auto z = random_double(-1, 1);
    auto a = random_double(0, 2*pi);
    auto r = sqrt(1 - z*z);
    return vec3(r*cos(a), r*sin(a), z);
}

C++ calculates sqrt(1 - z*z) without the .max(0.0) guard. The Rust implementation adds .max(0.0) as protection against floating-point error.

r304-sphere-directions-mc/src/lib.rs

const PI: f64 = std::f64::consts::PI;
const EXACT_INTEGRAL: f64 = 4.0 * PI / 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)
    }
}

fn sample_unit_vector(rng: &mut XorShift64) -> (f64, f64, f64) {
    let z = 1.0 - 2.0 * 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 sphere_pdf() -> f64 {
    1.0 / (4.0 * PI)
}

fn estimate_integral(samples: usize, rng: &mut XorShift64) -> f64 {
    let pdf = sphere_pdf();
    let mut sum = 0.0;
    for _ in 0..samples {
        let (_, _, z) = sample_unit_vector(rng);
        let cosine_squared = z * z;
        sum += cosine_squared / pdf;
    }
    sum / samples as f64
}

pub fn generate_report() -> String {
    let mut output = String::new();
    output.push_str("MC Integration on the Sphere of Directions\n");
    output.push_str("========================================\n\n");
    output.push_str("Target integral: integral_{S^2} cos^2(theta) d_omega\n");
    output.push_str("Sampling PDF: p(direction) = 1/(4*pi)\n");
    output.push_str(&format!("Exact value: {:.12}\n\n", EXACT_INTEGRAL));

    let mut single_rng = XorShift64::new(0x3040_0001u64);
    let single_estimate = estimate_integral(1_000_000, &mut single_rng);
    output.push_str("Single run (N=1,000,000):\n");
    output.push_str("samples,estimate,abs_error\n");
    output.push_str(&format!(
        "1000000,{:.12},{:.12}\n\n",
        single_estimate,
        (single_estimate - EXACT_INTEGRAL).abs()
    ));

    output.push_str("Convergence checkpoints:\n");
    output.push_str("samples,estimate,abs_error\n");

    let checkpoints = [100usize, 1_000, 10_000, 100_000, 500_000];
    for (index, n) in checkpoints.iter().copied().enumerate() {
        let mut rng = XorShift64::new(0x3040_1000u64 + index as u64);
        let est = estimate_integral(n, &mut rng);
        output.push_str(&format!(
            "{},{:.12},{:.12}\n",
            n,
            est,
            (est - EXACT_INTEGRAL).abs()
        ));
    }

    output
}

Results in the Browser

abs_error is the absolute difference from the exact value $4\pi /3$. Because each checkpoint uses an independent RNG, the error does not decrease monotonically, but it becomes smaller overall as it fluctuates.

Summary

• I derived ${\int }_{S^2}{\cos }^2\theta d\omega =4\pi /3$ by integration in spherical coordinates.
• I implemented uniform spherical sampling using $z=1-2u$ and $\varphi =2\pi v$. Choosing $z$ uniformly accounts for the $\sin \theta d\theta$ area element.
• I empirically confirmed that the Monte Carlo estimator ${\hat{I}}_N=\frac{1}{N}\sum z_i^2/(1/4\pi )$ converges to $4\pi /3$.

The next chapter connects the idea of a spherical PDF to the ray tracer's scattering model.