9. 直接的なライトのサンプリング

Ray Tracing: The Rest of Your Life (v3.2.3): 9 Sampling Lights Directly / 3.9 直接的なライトのサンプリング

前章まで，散乱方向は cosine 重み付き分布でサンプリングしてきました。しかしこれでは光源の方向に特別な重みを与えていないため，光源が小さいほど多くのレイがはずれて収束が遅くなります。光源の面積を直接使うサンプリング PDF を導入することでこの問題を解決します。

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

• 光源面積 PDF から方向 PDF への変換：$p(\text{direction})=d^2/(\cos \alpha \cdot A)$
• flip_face による片面発光：天井の隙間からのノイズを防ぐ
• cosine_pdf クラスの設計：ONB を使った PDF オブジェクト化（次章の混合 PDF への布石）

光源 PDF の導出

面積 $A$ の光源上の点 $q$ を一様にサンプリングするとき，面積の PDF は $p_q(q)=1/A$ です。これを方向の PDF $p(\text{direction})$ に変換します。

面積素 $dA$ と立体角素 $d\omega$ の関係は，光源面の法線と方向のなす角 $\alpha$ を使って：

$$d\omega =\frac{dA\cdot \cos \alpha }{d^2(p,q)}$$

サンプリングの確率が等しい（$p(\text{direction})d\omega =p_q(q)dA$）から：

$$p(\text{direction})\cdot \frac{dA\cdot \cos \alpha }{d^2}=\frac{dA}{A}$$$$\boxed{{\displaystyle p(\text{direction})=\frac{d^2(p,q)}{\cos \alpha \cdot A}}}$$

この $p$ を使う MC 推定器の 1 サンプル寄与は：

$$\frac{f(\text{direction})}{p(\text{direction})}=\frac{f(\text{direction})\cdot \cos \alpha \cdot A}{d^2}$$

片面発光（flip_face）

コーネルボックスの光源は天井に埋め込まれており，発光は下方向だけを向いています。DiffuseLight の emitted() が rec.front_face を参照して裏面は黒を返すようにし，光源面の法線を反転させる flip_face ラッパーで向きを揃えます。これにより天井の隙間からのライトアップを防げます。

cosine_pdf クラス

前章（3.8）の ONB を「PDF オブジェクト」として封じ込めたのが cosine_pdf クラスです。value(direction) で PDF 値を，generate() でサンプル方向を取得できます。pdf 抽象基底クラスは 2 つのメソッドを持つシンプルなインタフェースで，次章（3.10 混合 PDF）で hittable_pdf と組み合わせます。

Rust 実装

r309-light-sampling クレートは，コーネルボックスのフロア中央から光源の立体角 $\mathrm{\Omega }$ を数値推定し，一様半球サンプリングと光源 PDF サンプリングの収束速度を比較します。

立体角は面積積分の変換式から：

$$\mathrm{\Omega }={\int }_{\text{光源}}\frac{\cos \alpha }{d^2}dA$$

一様半球サンプリング： ランダム方向がたまたま光源に当たる頻度から推定します。

$${\hat{\mathrm{\Omega }}}_{\text{hemi}}=\frac{2\pi }{N}\sum _{i=1}^N\mathbf{1}[\text{レイが光源に当たる}]$$

光源 PDF サンプリング： 光源上の点を直接サンプリングして寄与を計算します。

$${\hat{\mathrm{\Omega }}}_{\text{light}}=\frac{1}{N}\sum _{i=1}^N\frac{\cos {\alpha }_i\cdot A}{d_i^2}$$

r309-light-sampling/src/lib.rs

/// Returns true if a ray from P in direction (dx, dy, dz) hits the light.
fn ray_hits_light(dx: f64, dy: f64, dz: f64) -> bool {
    if dy < 1e-10 { return false; }
    let t = (LIGHT_Y - P[1]) / dy;
    let x = P[0] + dx * t;
    let z = P[2] + dz * t;
    x >= LIGHT_X0 && x <= LIGHT_X1 && z >= LIGHT_Z0 && z <= LIGHT_Z1
}

光源 PDF の 1 サンプル寄与は，光源上の点 $q=(q_x,y_{\text{light}},q_z)$ をランダムに選んで計算します：

r309-light-sampling/src/lib.rs

let qx = LIGHT_X0 + rng_l.next_f64() * (LIGHT_X1 - LIGHT_X0);
let qz = LIGHT_Z0 + rng_l.next_f64() * (LIGHT_Z1 - LIGHT_Z0);
let tx = qx - P[0];
let ty = LIGHT_Y - P[1]; // always positive
let tz = qz - P[2];
let d2 = tx * tx + ty * ty + tz * tz;
let cos_alpha = ty / d2.sqrt();
if cos_alpha > 1e-10 {
    // contribution = cos(α) · A / d²  (= 1 / p_dir)
    sum_l += cos_alpha * LIGHT_AREA / d2;
}

generate_report は $N=10,100,\dots ,{10}^6$ の 6 チェックポイントで両推定値を並べて出力します。光源面積は半球面積の約 0.7% しかないため，一様半球では $N={10}^4$ 程度まで推定値がほぼゼロになる一方，光源 PDF では小さな $N$ から安定した収束が確認できます。

C++ と Rust の違い

C++ 版は ray_color() に光源サンプリングを直接ハードコードして動作を確認します：

auto on_light = point3(random_double(213,343), 554, random_double(227,332));
auto to_light = on_light - rec.p;
auto distance_squared = to_light.length_squared();
to_light = unit_vector(to_light);

if (dot(to_light, rec.normal) < 0)
    return emitted;

double light_area = (343-213)*(332-227);
auto light_cosine = fabs(to_light.y());
if (light_cosine < 0.000001)
    return emitted;

pdf = distance_squared / (light_cosine * light_area);
scattered = ray(rec.p, to_light, r.time());

flip_face と cosine_pdf は C++ 版で次のように定義されます：

// flip_face: front_face のときのみ発光
if (rec.front_face)
    return emit->value(u, v, p);
else
    return color(0,0,0);

// cosine_pdf: ONB を使った PDF オブジェクト
class cosine_pdf : public pdf {
  public:
    cosine_pdf(const vec3& w) { uvw.build_from_w(w); }
    double value(const vec3& direction) const override {
        auto cosine = dot(unit_vector(direction), uvw.w());
        return (cosine <= 0) ? 0 : cosine/pi;
    }
    vec3 generate() const override {
        return uvw.local(random_cosine_direction());
    }
    onb uvw;
};

Rust 版ではこの章は数値デモに留め，PDF トレイトの設計と flip_face の実装はコーネルボックスレンダラー実装時に行います。C++ の pdf 抽象基底クラスに対応する Rust のインタフェースは trait で表現します。

r309-light-sampling/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)
    }
}

// Rectangular light: x ∈ [213, 343], z ∈ [227, 332], y = 554
const LIGHT_X0: f64 = 213.0;
const LIGHT_X1: f64 = 343.0;
const LIGHT_Z0: f64 = 227.0;
const LIGHT_Z1: f64 = 332.0;
const LIGHT_Y: f64 = 554.0;
const LIGHT_AREA: f64 = (LIGHT_X1 - LIGHT_X0) * (LIGHT_Z1 - LIGHT_Z0); // = 13650.0

// Hit point: center of Cornell box floor
const P: [f64; 3] = [278.0, 0.0, 279.5];

fn random_uniform_hemi(rng: &mut XorShift64) -> (f64, f64, f64) {
    let r1 = rng.next_f64();
    let r2 = rng.next_f64();
    let y = 1.0 - r2; // cos(θ) = 1 - r2, so y ∈ (0, 1]
    let r = (1.0 - y * y).sqrt();
    let phi = 2.0 * PI * r1;
    (phi.cos() * r, y, phi.sin() * r)
}

/// Returns true if a ray from P in direction (dx, dy, dz) hits the light.
fn ray_hits_light(dx: f64, dy: f64, dz: f64) -> bool {
    if dy < 1e-10 {
        return false;
    }
    let t = (LIGHT_Y - P[1]) / dy;
    let x = P[0] + dx * t;
    let z = P[2] + dz * t;
    x >= LIGHT_X0 && x <= LIGHT_X1 && z >= LIGHT_Z0 && z <= LIGHT_Z1
}

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

    out.push_str("=== コーネルボックスの光源立体角を推定 ===\n");
    out.push_str(&format!("  ヒット点 p = ({}, {}, {})\n", P[0], P[1], P[2]));
    out.push_str(&format!(
        "  光源：x ∈ [{}, {}]，z ∈ [{}, {}]，y = {}\n",
        LIGHT_X0, LIGHT_X1, LIGHT_Z0, LIGHT_Z1, LIGHT_Y
    ));
    out.push_str(&format!("  光源面積 A = {:.1}\n\n", LIGHT_AREA));

    let checkpoints: [usize; 6] = [10, 100, 1_000, 10_000, 100_000, 1_000_000];
    let n_total = *checkpoints.last().unwrap();

    let mut hemi_est = Vec::with_capacity(checkpoints.len());
    let mut light_est = Vec::with_capacity(checkpoints.len());

    let mut rng_h = XorShift64::new(0x1a2b_3c4d_0001);
    let mut rng_l = XorShift64::new(0x5e6f_7a8b_0001);
    let mut count_h = 0usize;
    let mut sum_l = 0.0f64;
    let mut cp_idx = 0usize;

    for i in 0..n_total {
        let (dx, dy, dz) = random_uniform_hemi(&mut rng_h);
        if ray_hits_light(dx, dy, dz) {
            count_h += 1;
        }

        let qx = LIGHT_X0 + rng_l.next_f64() * (LIGHT_X1 - LIGHT_X0);
        let qz = LIGHT_Z0 + rng_l.next_f64() * (LIGHT_Z1 - LIGHT_Z0);
        let tx = qx - P[0];
        let ty = LIGHT_Y - P[1]; // always positive
        let tz = qz - P[2];
        let d2 = tx * tx + ty * ty + tz * tz;
        let d = d2.sqrt();
        let cos_alpha = ty / d; // angle at light face
        if cos_alpha > 1e-10 {
            // contribution = cos(α) · A / d²  (= 1 / p_dir)
            sum_l += cos_alpha * LIGHT_AREA / d2;
        }

        let n = i + 1;
        if cp_idx < checkpoints.len() && n == checkpoints[cp_idx] {
            hemi_est.push(2.0 * PI * count_h as f64 / n as f64);
            light_est.push(sum_l / n as f64);
            cp_idx += 1;
        }
    }

    out.push_str("=== 収束比較：立体角 Ω の推定値 (sr) ===\n");
    out.push_str(&format!(
        "{:>10}  {:>22}  {:>22}\n",
        "N", "一様半球サンプリング", "光源 PDF サンプリング"
    ));
    out.push_str(&format!("{:->10}  {:->22}  {:->22}\n", "", "", ""));
    for (i, &n) in checkpoints.iter().enumerate() {
        out.push_str(&format!(
            "{:>10}  {:>22.8}  {:>22.8}\n",
            n, hemi_est[i], light_est[i]
        ));
    }

    out
}

#[cfg(test)]
mod tests {
    use super::*;

    /// Light PDF estimator converges to the correct solid angle.
    #[test]
    fn light_pdf_converges() {
        let n = 500_000;
        let mut rng = XorShift64::new(0xdead_beef_0001);
        let mut sum = 0.0f64;
        for _ in 0..n {
            let qx = LIGHT_X0 + rng.next_f64() * (LIGHT_X1 - LIGHT_X0);
            let qz = LIGHT_Z0 + rng.next_f64() * (LIGHT_Z1 - LIGHT_Z0);
            let tx = qx - P[0];
            let ty = LIGHT_Y - P[1];
            let tz = qz - P[2];
            let d2 = tx * tx + ty * ty + tz * tz;
            let d = d2.sqrt();
            let cos_alpha = ty / d;
            if cos_alpha > 1e-10 {
                sum += cos_alpha * LIGHT_AREA / d2;
            }
        }
        let omega = sum / n as f64;
        // Solid angle should be roughly 0.04–0.06 sr (small light far away)
        assert!(
            omega > 0.03 && omega < 0.07,
            "light PDF solid angle estimate {omega:.6} out of expected range"
        );
    }

    /// Hemisphere estimator converges to same solid angle (needs many samples).
    #[test]
    fn hemi_converges_same_as_light_pdf() {
        let n = 1_000_000;
        let mut rng_l = XorShift64::new(0xdead_beef_0002);
        let mut sum_l = 0.0f64;
        for _ in 0..n {
            let qx = LIGHT_X0 + rng_l.next_f64() * (LIGHT_X1 - LIGHT_X0);
            let qz = LIGHT_Z0 + rng_l.next_f64() * (LIGHT_Z1 - LIGHT_Z0);
            let tx = qx - P[0];
            let ty = LIGHT_Y - P[1];
            let tz = qz - P[2];
            let d2 = tx * tx + ty * ty + tz * tz;
            let d = d2.sqrt();
            let cos_alpha = ty / d;
            if cos_alpha > 1e-10 {
                sum_l += cos_alpha * LIGHT_AREA / d2;
            }
        }
        let omega_light = sum_l / n as f64;

        let mut rng_h = XorShift64::new(0xdead_beef_0003);
        let count_h: usize = (0..n)
            .filter(|_| {
                let (dx, dy, dz) = random_uniform_hemi(&mut rng_h);
                ray_hits_light(dx, dy, dz)
            })
            .count();
        let omega_hemi = 2.0 * PI * count_h as f64 / n as f64;

        // Both should agree within 5%
        let rel_err = (omega_hemi - omega_light).abs() / omega_light;
        assert!(
            rel_err < 0.05,
            "hemi={omega_hemi:.6} vs light={omega_light:.6}, rel_err={rel_err:.4}"
        );
    }
}

ブラウザ上での実行結果

光源面積は上半球面積 $2\pi$ の約 0.7% しかないため，一様半球サンプリングでは $N={10}^4$ 程度まで推定値がほぼゼロになります。光源 PDF サンプリングは $N=10$ の段階からすでに安定した値を示し，同じ精度を達成するのに必要なサンプル数が大幅に少ないことが確認できます。

コーネルボックスのレンダリング

r309-cornell-box クレートは，FlipFace を適用した片面発光でコーネルボックスをレンダリングします。DiffuseLight::emitted が rec.front_face を参照するため，天井の隙間からのライトアップが発生しません。

まとめ

• 光源の面積 PDF $p_q=1/A$ を方向 PDF $p(\text{direction})=d^2/(\cos \alpha \cdot A)$ に変換する公式を導出した。
• 光源が小さいほど一様半球サンプリングは非効率になり，光源直接サンプリングで分散を劇的に削減できる。
• flip_face で片面発光を実現し，cosine_pdf クラスで ONB ベースのサンプリングを PDF オブジェクトとして封じ込める設計を確認した。

次章では cosine_pdf と光源 PDF を混合する 混合 PDF を導入し，コーネルボックスのレンダリングに統合します。