6. 重点サンプリング付きマテリアル

Ray Tracing: The Rest of Your Life (v3.2.3): 6 Importance Sampling Materials / 3.6 重点サンプリング付きマテリアル

3.5 で Lambertian 表面の散乱 PDF $s(\omega )=\cos \theta /\pi$ を導入しました。ここではサンプリング PDF $p$ と散乱 PDF $s$ の関係を整理し，その違いが収束速度と最終的な積分値にどう影響するかを数値実験で確かめます。

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

• $p=s$（一致），$p\ne s$（不一致），$s$ 自体を変える（別マテリアル）の 3 通りの比較
• テスト関数 $f(\omega )={\cos }^2\theta$ を使った数値実験

scatter_pdf と pdf_value の枠組み

レンダリング方程式のモンテカルロ推定は

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

と書けます。ここで

• $s(\omega )$：マテリアルの散乱 PDF（scattering_pdf）
• $p(\omega )$：方向をサンプリングするためのサンプリング PDF（pdf_value）

$s$ が確定すると積分値（物理的な色）は確定します。$p$ の選び方は積分値を変えず，収束速度（分散）のみに影響します。

3 通りの組み合わせ

原典 6.1〜6.3 では以下の 3 通りを比較しています。

方法	散乱 PDF $s$	サンプリング PDF $p$	積分値	収束速度
A	$\cos \theta /\pi$	$\cos \theta /\pi$（完全一致）	$1/2$	高（分散 0 に近い）
B	$\cos \theta /\pi$	$1/(2\pi )$（不一致）	$1/2$	低
C	$1/(2\pi )$	$1/(2\pi )$	$1/3$	高（ただし別マテリアル）

方法 A と B は $s$ が同じ（Lambertian）なので積分値は等しく，$p$ の選択は収束速度だけに影響します。方法 C は $s$ 自体を変えているため，根本的に異なるマテリアル（一様半球拡散）になり，別の積分値を与えます。

原典ではこの 3 通りの違いをコーネルボックスのレンダリング画像で示しています。本章では数値実験で同じ結論（方法 A・B は $1/2$，方法 C は $1/3$）を確認するにとどめ，コーネルボックスのライブレンダリングは実装しません。コーネルボックスレンダラーは 3.9「光源の直接サンプリング」で初めて導入し，以降の章で段階的に改善していきます。

数値実験：${\cos }^2\theta$ の重み付き積分（方法 A・B・C の比較）

議論を具体化するため，テスト関数 $f(\omega )={\cos }^2\theta$ に対する重み付き積分

$$I_s={\int }_{\text{hemisphere}}f(\omega )\cdot s(\omega )d\omega$$

を評価します。

Lambertian の場合：

$$I_s^{\text{Lambert}}=\frac{1}{\pi }{\int }_{\text{hemisphere}}{\cos }^3\theta d\omega =\frac{1}{\pi }{\int }_0^{2\pi }{\int }_0^{\pi /2}{\cos }^3\theta \sin \theta d\theta d\varphi =2{\int }_0^1{u}^{3}du=\frac{1}{2}$$

一様半球の場合：

$$I_s^{\text{uniform}}=\frac{1}{2\pi }{\int }_{\text{hemisphere}}{\cos }^2\theta d\omega =\frac{1}{2\pi }{\int }_0^{2\pi }{\int }_0^{\pi /2}{\cos }^2\theta \sin \theta d\theta d\varphi ={\int }_0^1{u}^{2}du=\frac{1}{3}$$

Rust 実装

r306-importance-sampling-play クレートは 3 通りの推定器をそれぞれ独立した RNG で実行します。

sample_cosine は Malley's method の公式 $\cos \theta =\sqrt{1-r_2}$（$r_2$ は一様乱数）を使います。

r306-importance-sampling-play/src/lib.rs

// Sample a direction from the cosine-weighted hemisphere (z = cosθ).
fn sample_cosine(rng: &mut XorShift64) -> f64 {
    let r2 = rng.next_f64();
    rng.next_f64(); // consume r1 (we only need cosθ = sqrt(1-r2))
    (1.0 - r2).sqrt()
}

// Sample a direction from the uniform hemisphere (z = cosθ = 1 - r2).
fn sample_uniform(rng: &mut XorShift64) -> f64 {
    let r2 = rng.next_f64();
    rng.next_f64(); // consume r1
    1.0 - r2
}

両関数が 2 回 RNG を消費するのは，実際は 2D 方向（$\varphi$ と $\theta$）のサンプルである構造を反映しています。ここでは $\varphi$ 方向（r1）は不要なので消費のみ行います。

3 通りの推定器はそれぞれ $f(\omega )\cdot s(\omega )/p(\omega )$ を計算します。方法 A では $s/p=1$ となるため，寄与は $f(\omega )={\cos }^2\theta$ のみです。

r306-importance-sampling-play/src/lib.rs

// Method A: Lambertian scatter + cosine sampling (perfect match).
// contribution = f(d) · scatter_pdf(d) / sampling_pdf(d) = cos²θ
fn estimate_method_a(n: usize, rng: &mut XorShift64) -> f64 {
    let sum: f64 = (0..n)
        .map(|_| {
            let cos_theta = sample_cosine(rng);
            let s = scatter_pdf_lambertian(cos_theta);
            let p = pdf_cosine(cos_theta);
            f(cos_theta) * s / p
        })
        .sum();
    sum / n as f64
}

// Method B: Lambertian scatter + uniform sampling (imperfect match).
// contribution = f(d) · scatter_pdf_lambertian(d) / sampling_pdf_uniform(d)
//             = cos²θ · (cosθ/π) / (1/(2π)) = 2cos³θ
fn estimate_method_b(n: usize, rng: &mut XorShift64) -> f64 {
    let sum: f64 = (0..n)
        .map(|_| {
            let cos_theta = sample_uniform(rng);
            let s = scatter_pdf_lambertian(cos_theta);
            let p = pdf_uniform(cos_theta);
            f(cos_theta) * s / p
        })
        .sum();
    sum / n as f64
}

// Method C: Uniform hemisphere scatter + uniform sampling (different material).
// contribution = f(d) · scatter_pdf_uniform(d) / sampling_pdf_uniform(d) = cos²θ
// Expected value = 1/3 (different from Lambertian → different render result)
fn estimate_method_c(n: usize, rng: &mut XorShift64) -> f64 {
    let sum: f64 = (0..n)
        .map(|_| {
            let cos_theta = sample_uniform(rng);
            let s = scatter_pdf_uniform_hemi(cos_theta);
            let p = pdf_uniform(cos_theta);
            f(cos_theta) * s / p
        })
        .sum();
    sum / n as f64
}

generate_report は $N={10}^6$ の単発比較に続き，$N=100$ から $500000$ までの収束チェックポイントで方法 A と B を並べて出力します。

C++ と Rust の違い

C++ 版では ray_color() 関数内で

double scattering_pdf = rec.mat->scattering_pdf(r, rec, scattered);
double pdf_value = scattering_pdf; // or a different PDF
color color_from_scatter =
    (attenuation * scattering_pdf * ray_color(scattered, depth-1, world))
    / pdf_value;

として比較しています。Rust 版の estimate_method_a/b/c はこれと同じ $s/p$ の計算を数値デモとして実装したものです。

r306-importance-sampling-play/src/lib.rs

use std::fmt::Write;

const PI: f64 = std::f64::consts::PI;

// Exact values of the integrals:
//   ∫_hemi f(d) · scatter_pdf(d) dω  where  f(d) = cos²θ
// Lambertian scatter_pdf = cosθ/π → integral = 1/2
// Uniform   scatter_pdf = 1/(2π) → integral = 1/3
const LAMBERTIAN_INTEGRAL: f64 = 0.5;
const UNIFORM_INTEGRAL: f64 = 1.0 / 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)
    }
}

// Sample a direction from the cosine-weighted hemisphere (z = cosθ).
fn sample_cosine(rng: &mut XorShift64) -> f64 {
    let r2 = rng.next_f64();
    rng.next_f64(); // consume r1 (we only need cosθ = sqrt(1-r2))
    (1.0 - r2).sqrt()
}

// Sample a direction from the uniform hemisphere (z = cosθ = 1 - r2).
fn sample_uniform(rng: &mut XorShift64) -> f64 {
    let r2 = rng.next_f64();
    rng.next_f64(); // consume r1
    1.0 - r2
}

fn pdf_cosine(cos_theta: f64) -> f64 {
    cos_theta / PI
}

fn pdf_uniform(_cos_theta: f64) -> f64 {
    1.0 / (2.0 * PI)
}

fn scatter_pdf_lambertian(cos_theta: f64) -> f64 {
    cos_theta / PI
}

fn scatter_pdf_uniform_hemi(_cos_theta: f64) -> f64 {
    1.0 / (2.0 * PI)
}

// f(d) = cos²θ  (arbitrary smooth function for demonstration)
fn f(cos_theta: f64) -> f64 {
    cos_theta * cos_theta
}

// Method A: Lambertian scatter + cosine sampling (perfect match).
// contribution = f(d) · scatter_pdf(d) / sampling_pdf(d) = cos²θ
fn estimate_method_a(n: usize, rng: &mut XorShift64) -> f64 {
    let sum: f64 = (0..n)
        .map(|_| {
            let cos_theta = sample_cosine(rng);
            let s = scatter_pdf_lambertian(cos_theta);
            let p = pdf_cosine(cos_theta);
            f(cos_theta) * s / p
        })
        .sum();
    sum / n as f64
}

// Method B: Lambertian scatter + uniform sampling (imperfect match).
// contribution = f(d) · scatter_pdf_lambertian(d) / sampling_pdf_uniform(d)
//             = cos²θ · (cosθ/π) / (1/(2π)) = 2cos³θ
fn estimate_method_b(n: usize, rng: &mut XorShift64) -> f64 {
    let sum: f64 = (0..n)
        .map(|_| {
            let cos_theta = sample_uniform(rng);
            let s = scatter_pdf_lambertian(cos_theta);
            let p = pdf_uniform(cos_theta);
            f(cos_theta) * s / p
        })
        .sum();
    sum / n as f64
}

// Method C: Uniform hemisphere scatter + uniform sampling (different material).
// contribution = f(d) · scatter_pdf_uniform(d) / sampling_pdf_uniform(d) = cos²θ
// Expected value = 1/3 (different from Lambertian → different render result)
fn estimate_method_c(n: usize, rng: &mut XorShift64) -> f64 {
    let sum: f64 = (0..n)
        .map(|_| {
            let cos_theta = sample_uniform(rng);
            let s = scatter_pdf_uniform_hemi(cos_theta);
            let p = pdf_uniform(cos_theta);
            f(cos_theta) * s / p
        })
        .sum();
    sum / n as f64
}

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

    // Header
    writeln!(out, "Integral: ∫_hemi cos²θ · scatter_pdf(d) dω").unwrap();
    writeln!(out).unwrap();

    // Single-run comparison with N = 1_000_000
    let n = 1_000_000;
    let mut rng_a = XorShift64::new(0x1a2b_3c4d_0001);
    let mut rng_b = XorShift64::new(0x1a2b_3c4d_0002);
    let mut rng_c = XorShift64::new(0x1a2b_3c4d_0003);

    let result_a = estimate_method_a(n, &mut rng_a);
    let result_b = estimate_method_b(n, &mut rng_b);
    let result_c = estimate_method_c(n, &mut rng_c);

    writeln!(out, "=== Single Run (N = {n}) ===").unwrap();
    writeln!(out, "{:<10} {:<12} {:<12} {:<12}", "Method", "Estimate", "Exact", "Error").unwrap();
    writeln!(
        out,
        "{:<10} {:<12.8} {:<12.8} {:<+12.2e}",
        "A (cos/cos)",
        result_a,
        LAMBERTIAN_INTEGRAL,
        result_a - LAMBERTIAN_INTEGRAL
    )
    .unwrap();
    writeln!(
        out,
        "{:<10} {:<12.8} {:<12.8} {:<+12.2e}",
        "B (cos/uni)",
        result_b,
        LAMBERTIAN_INTEGRAL,
        result_b - LAMBERTIAN_INTEGRAL
    )
    .unwrap();
    writeln!(
        out,
        "{:<10} {:<12.8} {:<12.8} {:<+12.2e}",
        "C (uni/uni)",
        result_c,
        UNIFORM_INTEGRAL,
        result_c - UNIFORM_INTEGRAL
    )
    .unwrap();
    writeln!(out).unwrap();

    // Convergence table for A and B (same target, different PDFs)
    writeln!(out, "=== Convergence: Method A (cos/cos) vs B (cos/uni) ===").unwrap();
    writeln!(
        out,
        "{:<10} {:<14} {:<12} {:<14} {:<12}",
        "N", "A Estimate", "A Error", "B Estimate", "B Error"
    )
    .unwrap();

    let checkpoints = [100usize, 1_000, 10_000, 100_000, 500_000];
    let mut rng_a = XorShift64::new(0x5e6f_7a8b_0001);
    let mut rng_b = XorShift64::new(0x5e6f_7a8b_0002);

    for &n in &checkpoints {
        let a = estimate_method_a(n, &mut rng_a);
        let b = estimate_method_b(n, &mut rng_b);
        writeln!(
            out,
            "{:<10} {:<14.8} {:<12.2e} {:<14.8} {:<12.2e}",
            n,
            a,
            a - LAMBERTIAN_INTEGRAL,
            b,
            b - LAMBERTIAN_INTEGRAL
        )
        .unwrap();
    }
    writeln!(out).unwrap();

    // Method C comparison
    writeln!(out, "=== Note: Method C gives a different answer ===").unwrap();
    writeln!(
        out,
        "Lambertian scatter_pdf (A & B): exact = {}",
        LAMBERTIAN_INTEGRAL
    )
    .unwrap();
    writeln!(
        out,
        "Uniform scatter_pdf   (C):      exact = {:.8}",
        UNIFORM_INTEGRAL
    )
    .unwrap();

    out
}

ブラウザ上での実行結果

収束チェックポイントの表では，方法 A の A Error が方法 B の B Error より一貫して小さくなります。単発比較の Error 欄では方法 A と B の厳密値が同じ（1/2）であり，方法 C だけが 1/3 と異なることを確認できます。

まとめ

• レンダリング方程式の推定 $A\cdot s(\omega )/p(\omega )\cdot \text{color}(\omega )$ において，$s$ が積分値を決め，$p$ が収束速度を決める。
• $p=s$（方法 A）では分散が最小になり，$p\ne s$（方法 B）でも同じ積分値に収束する。
• $s$ 自体を変えると（方法 C）別マテリアルになり，積分値が変わる（$1/2\to 1/3$）。

次章では，この枠組みをレイトレーサ本体に接続します。