5. 光の散乱

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

3.3 〜 3.4 でモンテカルロ推定の枠組みと球面積分を学びました。ここでは光が表面で散乱する際の確率モデルを導入し，レイトレーサのノイズ低減につながる考え方を作ります。

扱うのは 3 つのトピックです。

• アルベードと散乱 PDF による散乱モデルの記述
• Lambertian 表面の cosine density と，それを使った importance sampling
• ${\int }_{\text{hemisphere}}\cos \theta d\omega =\pi$ の推定による数値確認

アルベードと散乱 PDF

表面に入射した光が散乱する確率を $A$（アルベード）とします。散乱した場合，方向分布は PDF $s(\mathbf{d})$ で記述できます。これを散乱 PDFと呼びます。

反射光の放射輝度は次の積分で与えられます。

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

モンテカルロ推定では，サンプリング PDF $p(\mathbf{d})$ を使って

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

と近似します。$p=s$ に設定すると $s/p=1$ となり，$\text{Color}\approx A\cdot \text{color}(\mathbf{d})$ と分散がゼロになります。これが散乱における importance sampling です。

Lambertian 散乱と Cosine density

Lambertian（拡散反射）表面の散乱 PDF は

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

です（$\theta$ は法線からの角度）。これを cosine density と呼びます。PDF が正規化されていることは

$${\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$$

で確認できます。

積分の厳密値

Lambertian 表面の放射輝度積分（アルベード $A=1$，一様入射光）は

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

です。球座標で展開すると

$$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$$

半球サンプリングの 2 手法

一様半球サンプリング

半球面積は $2\pi$ なので，一様 PDF は

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

です。$z=u\in [0,1)$，$\varphi =2\pi v$ とすると上半球を一様に覆えます（$z=\cos \theta \ge 0$）。

Cosine density サンプリング（マリーの方法）

$p(\cos \theta )$ の周辺分布を求めると $p(c)=2c$（$c=\cos \theta$）です。CDF は $P(c)=c^2$ なので逆関数法により

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

でサンプルできます。このとき推定量の各サンプルは

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

と定数になり，分散はゼロです。

Rust 実装

r305-scattering-pdf クレートは XorShift64 をローカルに実装し，2 手法の推定を比較します。

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
}

推定関数は $\cos \theta /p(\omega )$ を $N$ 回加算して平均を取ります。if cos_theta > 0.0 は浮動小数点誤差で $z$ がわずかに負になるケースへの防御です。

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
}

generate_report は $N={10}^6$ の単発推定に続き，$N=100$ から $500000$ までのチェックポイントで uniform と cosine の誤差を並べて出力します。

C++ と Rust の違い

C++ 原典の cosine hemisphere サンプリングは Malley's method と呼ばれ，単位円板上の一様サンプルを半球へ投影する方法で実装されることがあります。

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

これは $r=\sqrt{r_2}$（円板上の一様点）として $z=\sqrt{1-r_2}=\sqrt{1-r^2}=\cos \theta$ になります。Rust 実装も同様の逆関数法を使っています（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
}

ブラウザ上での実行結果

cosine_error の列は，理論上は各サンプルで $\pi$ に等しいため，誤差がほぼゼロになるはずです。実際は浮動小数点丸め誤差が残りますが，uniform_error と比較すると圧倒的に小さくなります。

まとめ

• アルベード $A$ と散乱 PDF $s(\mathbf{d})$ による散乱モデルを導入した。モンテカルロ推定では $A\cdot s(\mathbf{d})/p(\mathbf{d})$ が基本単位となる。
• Lambertian 表面の散乱 PDF は $s(\mathbf{d})=\cos \theta /\pi$（cosine density）であり，${\int }_{\text{hemisphere}}\cos \theta d\omega =\pi$ が厳密値。
• $p=s$（cosine PDF サンプリング）のとき推定量は定数になり分散がゼロとなる理想的な importance sampling を実現できる。

次章以降では，この考え方を複数 PDF の混合（光源サンプリングと表面サンプリング）へ拡張して，ノイズの少ないレイトレーサを実現します。