7. Generating Random Directions

Ray Tracing: The Rest of Your Life (v3.2.3): 7 Generating Random Directions / 3.7 Generating Random Directions

The Monte Carlo estimates in Sections 3.2 through 3.6 were one-dimensional discussions that considered only the scattering angle $\theta$ around the $z$ axis. I now systematize methods for generating random direction vectors in three-dimensional space. Combined with the ONB, or orthonormal basis, introduced in the next chapter, these methods allow scattering around any surface normal.

This chapter covers three topics:

• A general procedure for generating directions from a PDF in spherical coordinates using inverse transform sampling
• Specific formulas for uniform full-sphere, uniform hemisphere, and cosine-weighted sampling
• Estimating the same integral ${\int }_{\text{hemi}}{\cos }^3\theta d\omega$ with two PDFs and confirming agreement

Spherical Coordinates and Solid-Angle PDFs

The area element on the unit sphere is

$$dA=\sin \theta d\theta d\varphi$$

Given a PDF $p(\text{direction})=f(\theta )$ that is symmetric in $\varphi$, the marginal PDFs of $\theta$ and $\varphi$ are

$$a(\varphi )=\frac{1}{2\pi },b(\theta )=2\pi f(\theta )\sin \theta$$

A direction is sampled from random values $r_1,r_2\in [0,1)$ using inverse transform sampling.

$$r_1={\int }_0^{\varphi }\frac{1}{2\pi }dt=\frac{\varphi }{2\pi }\Rightarrow \varphi =2\pi r_1$$$$r_2={\int }_0^{\theta }2\pi f(t)\sin tdt$$

Solving the equation for $r_2$ in terms of $\cos \theta$ produces a specific formula for each PDF.

Uniform Full-Sphere Sampling

For $p(\text{direction})=1/(4\pi )$, normalized by the sphere's area $4\pi$:

$$r_2={\int }_0^{\theta }\frac{1}{2}\sin tdt=\frac{1-\cos \theta }{2}\Rightarrow \cos \theta =1-2r_2$$

Converting to Cartesian coordinates with $\sin \theta =\sqrt{1-{\cos }^2\theta }=2\sqrt{r_2(1-r_2)}$ gives

$$x=\cos (2\pi r_1)\cdot 2\sqrt{r_2(1-r_2)},y=\sin (2\pi r_1)\cdot 2\sqrt{r_2(1-r_2)},z=1-2r_2$$

Uniform Hemisphere Sampling

For $p(\text{direction})=1/(2\pi )$, normalized by the hemisphere's area $2\pi$:

$$r_2={\int }_0^{\theta }\sin tdt=1-\cos \theta \Rightarrow \cos \theta =1-r_2$$

The value $z=1-r_2$ ranges from $0$ at the horizon to $1$ directly overhead.

This distribution is used to estimate ${\int }_{\text{hemi}}{\cos }^3\theta d\omega$. The analytic solution is

$${\int }_0^{2\pi }{\int }_0^{\pi /2}{\cos }^3\theta \sin \theta d\theta d\varphi =2\pi {\int }_0^1{u}^{3}du=\frac{\pi }{2}\approx 1.5708$$

The Monte Carlo estimator is $f(\mathbf{d})/p(\mathbf{d})={\cos }^3\theta /(1/(2\pi ))=2\pi {\cos }^3\theta$.

Cosine-Weighted Sampling

For the Lambertian scattering PDF $p(\text{direction})=\cos \theta /\pi$:

$$r_2={\int }_0^{\theta }2\cos t\sin tdt=1-{\cos }^2\theta \Rightarrow \cos \theta =\sqrt{1-r_2}$$

Converting to Cartesian coordinates, with $\sin \theta =\sqrt{r_2}$, gives

$$z=\sqrt{1-r_2},x=\cos (2\pi r_1)\sqrt{r_2},y=\sin (2\pi r_1)\sqrt{r_2}$$

For the same integral ${\int }_{\text{hemi}}{\cos }^3\theta d\omega$, the estimator is

$$\frac{f(\mathbf{d})}{p(\mathbf{d})}=\frac{{\cos }^3\theta }{\cos \theta /\pi }=\pi {\cos }^2\theta$$

Both uniform-hemisphere sampling, with its higher variance, and cosine-weighted sampling, optimized for Lambertian scattering, converge to the same value $\pi /2$.

Rust Implementation

The r307-random-directions crate implements XorShift64 locally and provides three direction-generation functions.

r307-random-directions/src/lib.rs

/// Uniform distribution on the full unit sphere.
///
/// Inversion method: `φ = 2π·r1`, `cos(θ) = 1 - 2·r2`.
fn random_unit_sphere(rng: &mut XorShift64) -> (f64, f64, f64) {
    let r1 = rng.next_f64();
    let r2 = rng.next_f64();
    let z = 1.0 - 2.0 * r2;
    let r = (1.0 - z * z).sqrt(); // = 2·sqrt(r2·(1-r2))
    let phi = 2.0 * PI * r1;
    (phi.cos() * r, phi.sin() * r, z)
}

/// Uniform distribution on the hemisphere around +Z.
///
/// Inversion method: `φ = 2π·r1`, `cos(θ) = 1 - r2`.
fn random_uniform_hemi(rng: &mut XorShift64) -> (f64, f64, f64) {
    let r1 = rng.next_f64();
    let r2 = rng.next_f64();
    let z = 1.0 - r2;
    let r = (1.0 - z * z).sqrt();
    let phi = 2.0 * PI * r1;
    (phi.cos() * r, phi.sin() * r, z)
}

/// Cosine-weighted distribution on the hemisphere around +Z.
///
/// Inversion method: `φ = 2π·r1`, `cos(θ) = sqrt(1 - r2)`.
fn random_cosine_direction(rng: &mut XorShift64) -> (f64, f64, f64) {
    let r1 = rng.next_f64();
    let r2 = rng.next_f64();
    let z = (1.0 - r2).sqrt();
    let phi = 2.0 * PI * r1;
    let r = r2.sqrt(); // = sin(θ)
    (phi.cos() * r, phi.sin() * r, z)
}

generate_report produces three sections. It first displays ten example vectors from the uniform spherical distribution, then estimates ${\int }_{\text{hemi}}{\cos }^3\theta d\omega$ with $N={10}^6$ samples using both uniform-hemisphere and cosine-weighted sampling.

Differences Between C++ and Rust

In the C++ version, random_cosine_direction() returns a vec3 and uses the global random_double() function. For simplicity, the Rust version in this chapter returns a three-tuple, (f64, f64, f64). The next chapter uses the vector type when combining these directions with an ONB.

r307-random-directions/src/lib.rs

use std::f64::consts::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)
    }
}

/// Uniform distribution on the full unit sphere.
///
/// Inversion method: `φ = 2π·r1`, `cos(θ) = 1 - 2·r2`.
fn random_unit_sphere(rng: &mut XorShift64) -> (f64, f64, f64) {
    let r1 = rng.next_f64();
    let r2 = rng.next_f64();
    let z = 1.0 - 2.0 * r2;
    let r = (1.0 - z * z).sqrt(); // = 2·sqrt(r2·(1-r2))
    let phi = 2.0 * PI * r1;
    (phi.cos() * r, phi.sin() * r, z)
}

/// Uniform distribution on the hemisphere around +Z.
///
/// Inversion method: `φ = 2π·r1`, `cos(θ) = 1 - r2`.
fn random_uniform_hemi(rng: &mut XorShift64) -> (f64, f64, f64) {
    let r1 = rng.next_f64();
    let r2 = rng.next_f64();
    let z = 1.0 - r2;
    let r = (1.0 - z * z).sqrt();
    let phi = 2.0 * PI * r1;
    (phi.cos() * r, phi.sin() * r, z)
}

/// Cosine-weighted distribution on the hemisphere around +Z.
///
/// Inversion method: `φ = 2π·r1`, `cos(θ) = sqrt(1 - r2)`.
fn random_cosine_direction(rng: &mut XorShift64) -> (f64, f64, f64) {
    let r1 = rng.next_f64();
    let r2 = rng.next_f64();
    let z = (1.0 - r2).sqrt();
    let phi = 2.0 * PI * r1;
    let r = r2.sqrt(); // = sin(θ)
    (phi.cos() * r, phi.sin() * r, z)
}

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

    // --- Section 1: first 10 of 200 random unit vectors on the sphere -------
    out.push_str("=== Uniform Sphere: Example Unit Vectors (First 10 of 200) ===\n");
    out.push_str("         x           y           z\n");
    let mut rng1 = XorShift64::new(0x1a2b_3c4d_0001);
    for _ in 0..10 {
        let (x, y, z) = random_unit_sphere(&mut rng1);
        out.push_str(&format!("{:11.6} {:11.6} {:11.6}\n", x, y, z));
    }
    out.push_str("... (remaining 190 points omitted)\n\n");

    let n = 1_000_000usize;

    // --- Section 2: MC estimate of ∫ cos³θ dA via uniform hemisphere --------
    out.push_str("=== Estimate ∫_hemi cos³θ dA with Uniform Hemisphere Sampling ===\n");
    out.push_str(&format!("  π/2 (analytic solution) = {:.12}\n", PI / 2.0));
    let mut rng2 = XorShift64::new(0x1a2b_3c4d_0002);
    let sum2: f64 = (0..n)
        .map(|_| {
            let (_, _, z) = random_uniform_hemi(&mut rng2);
            // p(direction) = 1/(2π), f = cos³θ → f/p = 2π·cos³θ
            z * z * z / (1.0 / (2.0 * PI))
        })
        .sum();
    out.push_str(&format!(
        "  estimate (N={}) = {:.12}\n\n",
        n,
        sum2 / n as f64
    ));

    // --- Section 3: MC estimate of ∫ cos³θ dA via cosine direction ----------
    out.push_str("=== Estimate ∫_hemi cos³θ dA with Cosine-Weighted Sampling ===\n");
    out.push_str(&format!("  π/2 (analytic solution) = {:.12}\n", PI / 2.0));
    let mut rng3 = XorShift64::new(0x1a2b_3c4d_0003);
    let sum3: f64 = (0..n)
        .map(|_| {
            let (_, _, z) = random_cosine_direction(&mut rng3);
            // p(direction) = cos(θ)/π, f = cos³θ → f/p = π·cos²θ
            z * z * z / (z / PI)
        })
        .sum();
    out.push_str(&format!(
        "  estimate (N={}) = {:.12}\n",
        n,
        sum3 / n as f64
    ));

    out
}

pub fn report_random_directions() -> String {
    generate_report()
}

Results in the Browser

Both estimates are close to the analytic solution $\pi /2\approx 1.5708$. With cosine-weighted sampling, one cosine factor cancels to give $f/p=\pi {\cos }^2\theta$, so its variance tends to be lower than with uniform-hemisphere sampling.

Summary

• I established a general inverse-transform procedure for deriving $\varphi =2\pi r_1$ and $\cos \theta$ from a spherical-coordinate PDF $f(\theta )$.
• I derived formulas for uniform full-sphere sampling ($\cos \theta =1-2r_2$), uniform-hemisphere sampling ($\cos \theta =1-r_2$), and cosine-weighted sampling ($\cos \theta =\sqrt{1-r_2}$).
• I numerically confirmed that both uniform-hemisphere and cosine-weighted sampling produce ${\int }_{\text{hemi}}{\cos }^3\theta d\omega =\pi /2$.

The next chapter introduces an orthonormal basis (ONB) that transforms these directions into a coordinate system around any surface normal.