7. ランダムな方向の生成

Ray Tracing: The Rest of Your Life (v3.2.3): 7 Generating Random Directions / 3.7 ランダムな方向の生成

3.2〜3.6 のモンテカルロ推定は「$z$ 軸周りの散乱角 $\theta$ だけ考える 1 次元の議論」でした。ここでは3 次元空間の方向ベクトルを確率的に生成する手法を体系化します。次章で学ぶ ONB（正規直交基底）と組み合わせると，任意の法線方向に対して散乱させることができます。

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

• 球面座標の PDF から逆変換法で方向を生成する一般的な手順
• 全球面一様・半球一様・cosine 重み付きの 3 種類の具体的な公式
• 同じ積分 ${\int }_{\text{hemi}}{\cos }^3\theta d\omega$ を 2 通りの PDF で推定して一致を確認

球面座標と立体角の PDF

単位球面上の面積素は

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

です。PDF $p(\text{direction})=f(\theta )$（$\varphi$ について対称）が与えられたとき，$\theta$ と $\varphi$ の周辺 PDF は

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

となります。逆変換法で乱数 $r_1,r_2\in [0,1)$ から方向をサンプリングします。

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

$r_2$ の式を $\cos \theta$ について解くと，各 PDF に対する具体的な式が得られます。

全球面一様サンプリング

$p(\text{direction})=1/(4\pi )$（球面積 $4\pi$ で正規化）のとき：

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

$\sin \theta =\sqrt{1-{\cos }^2\theta }=2\sqrt{r_2(1-r_2)}$ を使って Cartesian 座標に変換すると：

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

半球一様サンプリング

$p(\text{direction})=1/(2\pi )$（半球面積 $2\pi$ で正規化）のとき：

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

$z=1-r_2$ は $0$（水平）から $1$（真上）まで変化します。

このサンプリングで ${\int }_{\text{hemi}}{\cos }^3\theta d\omega$ を推定します。解析解は

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

です。モンテカルロ推定量は $f(\mathbf{d})/p(\mathbf{d})={\cos }^3\theta /(1/(2\pi ))=2\pi {\cos }^3\theta$ です。

Cosine 重み付きサンプリング

$p(\text{direction})=\cos \theta /\pi$（Lambertian 散乱 PDF）のとき：

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

Cartesian 座標に変換すると（$\sin \theta =\sqrt{r_2}$）：

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

同じ ${\int }_{\text{hemi}}{\cos }^3\theta d\omega$ に対する推定量は

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

となります。半球一様（分散大）と cosine 重み付き（Lambertian に最適化）の両方が同じ $\pi /2$ に収束することを確認します。

Rust 実装

r307-random-directions クレートは XorShift64 をローカルに実装し，3 種類の方向生成関数を提供します。

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 は 3 つのセクションを出力します。まず球面一様分布のベクトル例を 10 点表示し，続いて半球一様サンプリングと cosine サンプリングそれぞれで ${\int }_{\text{hemi}}{\cos }^3\theta d\omega$ を $N={10}^6$ サンプルで推定します。

C++ と Rust の違い

C++ 版では random_cosine_direction() は vec3 を返し，グローバルな乱数関数（random_double()）を使います。Rust 版はこの章では 3 タプル (f64, f64, f64) を返すシンプルな実装にしています。次章（3.8）で 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("=== 球面一様分布：単位ベクトルの例（先頭 10 点 / 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("... (残り 190 点省略)\n\n");

    let n = 1_000_000usize;

    // --- Section 2: MC estimate of ∫ cos³θ dA via uniform hemisphere --------
    out.push_str("=== 一様半球サンプリングで ∫_hemi cos³θ dA を推定 ===\n");
    out.push_str(&format!("  π/2 (解析解) = {:.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!(
        "  推定値 (N={}) = {:.12}\n\n",
        n,
        sum2 / n as f64
    ));

    // --- Section 3: MC estimate of ∫ cos³θ dA via cosine direction ----------
    out.push_str("=== cosine 方向サンプリングで ∫_hemi cos³θ dA を推定 ===\n");
    out.push_str(&format!("  π/2 (解析解) = {:.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!(
        "  推定値 (N={}) = {:.12}\n",
        n,
        sum3 / n as f64
    ));

    out
}

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

ブラウザ上での実行結果

2 つの推定値がどちらも解析解 $\pi /2\approx 1.5708$ に近い値を示すことを確認できます。cosine 方向サンプリングは $f/p=\pi {\cos }^2\theta$ と $\cos$ の因子が 1 つ消えるため，一様半球より分散が小さくなる傾向があります。

まとめ

• 球面座標の PDF $f(\theta )$ から逆変換法で $\varphi =2\pi r_1$，$\cos \theta$ を求める一般手順を確立した。
• 全球面一様（$\cos \theta =1-2r_2$），半球一様（$\cos \theta =1-r_2$），cosine 重み付き（$\cos \theta =\sqrt{1-r_2}$）の 3 公式を導出した。
• 半球一様と cosine 重み付きの両方で ${\int }_{\text{hemi}}{\cos }^3\theta d\omega =\pi /2$ が得られることを数値確認した。

次章では ONB（正規直交基底）を導入し，これらの方向を任意の法線座標系に変換する仕組みを作ります。