4. 球上の点を使ったモンテカルロ積分

Ray Tracing: The Rest of Your Life (v3.2.3): 4 MC Integration on the Sphere of Directions / 3.4 球上の点を使ったモンテカルロ積分

3.3 では 1 次元区間上の積分を推定しました。ここではレイトレーシングに直接つながる方向空間（単位球面 $S^2$）上の積分を扱います。

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

• 球面上の一様サンプリングと対応する PDF
• ${\cos }^2\theta$ の球面積分による実測確認

積分の厳密値

題材の積分は

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

です。ここで $\theta$ は $z$ 軸との角度，$d\omega$ は立体角要素です。球座標では $d\omega =\sin \theta d\theta d\varphi$ なので

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

$u=\cos \theta$，$du=-\sin \theta d\theta$ と置換すると

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

球面一様サンプリング

単位球面上の一様分布を生成するには，$d\omega =\sin \theta d\theta d\varphi$ が面積素であることを利用します。面積要素が $\theta$ に依存するため，$\theta$ を一様にとっても $\varphi$ を一様にとっても球面を一様に覆えません。

正しいアプローチは $z=\cos \theta$ を一様にとることです。$z\in [-1,1]$ の一様分布は $\varphi \in [0,2\pi )$ の一様分布と組み合わせることで，球面を一様に覆います。$u\sim U[0,1)$ から

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

として方向 $(r\cos \varphi ,r\sin \varphi ,z)$ を作れます。この方向の PDF は球面積の逆数

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

です。

Rust 実装

この章用に r304-sphere-directions-mc クレートを追加しました。XorShift64 は章内でローカルに実装します。

球面サンプリングと PDF 定数は次の通りです。

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

.max(0.0) は浮動小数点誤差で $1-z^2$ がわずかに負になるケースを防ぎます。

推定本体は ${\cos }^2\theta =z^2$ を PDF で割り，$N$ 個の平均を取ります。

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 は 2 段階の出力を構成します。まず $N={10}^6$ の単発推定を出力し，続いて $N=100$ から $500000$ までの各チェックポイントを独立した 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()
        ));
    }

各チェックポイントは独立したシードを持つため，3.2 の累積カウント方式とは異なり，小 $N$ の推定が大 $N$ の結果に影響しません。

C++ と Rust の違い

C++ 原典の random_unit_vector は構造的に同じです。

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++ は sqrt(1 - z*z) で .max(0.0) による保護を行っていません。Rust 実装では浮動小数点誤差への防御として .max(0.0) を追加しています。

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
}

ブラウザ上での実行結果

abs_error は厳密値 $4\pi /3$ との差の絶対値です。各チェックポイントは独立した RNG を使うため，単調には減らず前後しながらも全体として誤差は小さくなっていきます。

まとめ

• ${\int }_{S^2}{\cos }^2\theta d\omega =4\pi /3$ を球座標積分で導いた。
• $z=1-2u$，$\varphi =2\pi v$ による球面一様サンプリングを実装した。$z$ の一様性が面積要素 $\sin \theta d\theta$ の補正に相当する。
• モンテカルロ推定量 ${\hat{I}}_N=\frac{1}{N}\sum z_i^2/(1/4\pi )$ が $4\pi /3$ に収束することを実測した。

次章では，この球面 PDF の考え方をレイトレーサの散乱モデルに接続します。