8. 正規直交基底

Ray Tracing: The Rest of Your Life (v3.2.3): 8 Orthonormal Bases / 3.8 正規直交基底

前章では $z$ 軸周りのランダム方向生成を学びました。実際のレンダリングでは散乱方向を任意の法線ベクトル $\mathbf{n}$ を基準に生成する必要があります。そのための道具が正規直交基底（ONB: Orthonormal Basis）です。

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

• ONB の定義と構築：法線 $\mathbf{n}$ から互いに直交する単位ベクトル $\mathbf{u},\mathbf{v},\mathbf{w}$ を求める手順
• ローカル座標からワールド座標への変換：$x\mathbf{u}+y\mathbf{v}+z\mathbf{w}$
• ONB 変換を介した cosine 方向サンプリングの検証（${\int }_{\text{hemi}}{\cos }^3\theta d\omega \approx \pi /2$）

相対座標系

原点 ${\mathbf{O}}^{\prime }$ と正規直交基底ベクトル $\mathbf{u},\mathbf{v},\mathbf{w}$ による座標系では，ある位置は

$${\mathbf{O}}^{\prime }+u\mathbf{u}+v\mathbf{v}+w\mathbf{w}$$

で表現されます。ONB の三本のベクトルは互いに直交する単位ベクトルです：

$$\mathbf{u}\cdot \mathbf{v}=\mathbf{u}\cdot \mathbf{w}=\mathbf{v}\cdot \mathbf{w}=0,|\mathbf{u}|=|\mathbf{v}|=|\mathbf{w}|=1$$

法線から ONB を構築する

$\mathbf{w}=\text{normalize}(\mathbf{n})$ を固定したとき，接平面ベクトル $\mathbf{v},\mathbf{u}$ は次の外積で得られます：

$$\mathbf{v}=\text{normalize}(\mathbf{w}\times \mathbf{a}),\mathbf{u}=\mathbf{w}\times \mathbf{v}$$

$\mathbf{a}$ は $\mathbf{w}$ と平行でなければ何でもよく，数値安定性のため

$$\mathbf{a}=\{\begin{array}{ll}(0,1,0)& |w_x|>0.9\\ (1,0,0)& \text{otherwise}\end{array}$$

とします。$|w_x|>0.9$ のとき $(1,0,0)$ を使うと $\mathbf{w}$ とほぼ平行になって外積がゼロベクトルに近くなるため，$y$ 軸に切り替えます。

ローカル座標からワールド座標への変換

前章の random_cosine_direction() は $z$ 軸周りのローカル座標 $(x,y,z)$ を返します。これをワールド座標に変換するには：

$$\text{方向}=x\mathbf{u}+y\mathbf{v}+z\mathbf{w}$$

local(a, b, c) はこの線形結合を，local_vec(d) は配列 $\mathbf{d}=(d_0,d_1,d_2)$ をそのまま渡す薄いラッパーです。

Rust 実装

r308-onb クレートは Onb 構造体と random_cosine_direction 関数を実装します。

r308-onb/src/lib.rs

/// Orthonormal basis built from a single surface normal.
///
/// `w` is aligned to the normal; `u` and `v` span the tangent plane.
struct Onb {
    u: [f64; 3],
    v: [f64; 3],
    w: [f64; 3],
}

impl Onb {
    /// Build an ONB whose `w` axis is aligned to `n`.
    fn build_from_w(n: [f64; 3]) -> Self {
        let w = normalize(n);
        // Pick an arbitrary vector not parallel to w.
        let a = if w[0].abs() > 0.9 {
            [0.0, 1.0, 0.0]
        } else {
            [1.0, 0.0, 0.0]
        };
        let v = normalize(cross(w, a));
        let u = cross(w, v);
        Self { u, v, w }
    }

    /// Transform local coordinates `(a, b, c)` to world space.
    ///
    /// Returns `a·u + b·v + c·w`.
    fn local(&self, a: f64, b: f64, c: f64) -> [f64; 3] {
        [
            a * self.u[0] + b * self.v[0] + c * self.w[0],
            a * self.u[1] + b * self.v[1] + c * self.w[1],
            a * self.u[2] + b * self.v[2] + c * self.w[2],
        ]
    }

    /// Transform a local direction vector `d` to world space.
    fn local_vec(&self, d: [f64; 3]) -> [f64; 3] {
        self.local(d[0], d[1], d[2])
    }
}

Lambertian 散乱では ONB を次のように使います。cos_theta は $\mathbf{w}$ との内積で得られるため，明示的な角度計算は不要です。

r308-onb/src/lib.rs

let local_dir = random_cosine_direction(&mut rng);
let world_dir = onb.local_vec(local_dir);
// cos(θ) relative to normal w = dot(world_dir, w)
let cos_theta = dot(world_dir, onb.w);
let pdf = cos_theta / PI;

generate_report は 2 つのセクションを出力します。まず normalize([1, 2, 3]) を法線として ONB を構築し，$\mathbf{u},\mathbf{v},\mathbf{w}$ の各ベクトルと直交性・単位長さの確認値を表示します。続いて $N={10}^6$ サンプルで ${\int }_{\text{hemi}}{\cos }^3\theta d\omega \approx \pi /2$ を推定します。前章（3.7）との違いは，random_cosine_direction の結果を onb.local_vec でワールド座標に変換してから $\cos \theta ={\mathbf{d}}_{\text{world}}\cdot \mathbf{w}$ を取る点です。ONB 変換は内積を保存するため，変換後でも $\cos \theta$ は正しく計算されます。

C++ と Rust の違い

C++ 版の onb クラスは axis[3] 配列を vec3 で保持し，operator[] で添字アクセスします。Rust 版では u/v/w フィールドを明示的に持つ構造体にしました。

C++ 版は build_from_w が void メソッドとして axis を書き換えます。Rust 版は build_from_w が Self を返すコンストラクタスタイルにしています（Rust では可変参照より値を返す方が自然なため）。

r308-onb/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)
    }
}

fn dot(a: [f64; 3], b: [f64; 3]) -> f64 {
    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
}

fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
    [
        a[1] * b[2] - a[2] * b[1],
        a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0],
    ]
}

fn normalize(v: [f64; 3]) -> [f64; 3] {
    let len = dot(v, v).sqrt();
    [v[0] / len, v[1] / len, v[2] / len]
}

/// Orthonormal basis built from a single surface normal.
///
/// `w` is aligned to the normal; `u` and `v` span the tangent plane.
struct Onb {
    u: [f64; 3],
    v: [f64; 3],
    w: [f64; 3],
}

impl Onb {
    /// Build an ONB whose `w` axis is aligned to `n`.
    fn build_from_w(n: [f64; 3]) -> Self {
        let w = normalize(n);
        // Pick an arbitrary vector not parallel to w.
        let a = if w[0].abs() > 0.9 {
            [0.0, 1.0, 0.0]
        } else {
            [1.0, 0.0, 0.0]
        };
        let v = normalize(cross(w, a));
        let u = cross(w, v);
        Self { u, v, w }
    }

    /// Transform local coordinates `(a, b, c)` to world space.
    ///
    /// Returns `a·u + b·v + c·w`.
    fn local(&self, a: f64, b: f64, c: f64) -> [f64; 3] {
        [
            a * self.u[0] + b * self.v[0] + c * self.w[0],
            a * self.u[1] + b * self.v[1] + c * self.w[1],
            a * self.u[2] + b * self.v[2] + c * self.w[2],
        ]
    }

    /// Transform a local direction vector `d` to world space.
    fn local_vec(&self, d: [f64; 3]) -> [f64; 3] {
        self.local(d[0], d[1], d[2])
    }
}

fn random_cosine_direction(rng: &mut XorShift64) -> [f64; 3] {
    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();
    [phi.cos() * r, phi.sin() * r, z]
}

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

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

    // --- Section 1: ONB construction and orthogonality check ----------------
    let n = normalize([1.0, 2.0, 3.0]);
    let onb = Onb::build_from_w(n);

    out.push_str("=== ONB の構築：n = normalize(1, 2, 3) ===\n");
    out.push_str(&format!(
        "  w (法線方向): [{:.6}, {:.6}, {:.6}]\n",
        onb.w[0], onb.w[1], onb.w[2]
    ));
    out.push_str(&format!(
        "  v (接線):     [{:.6}, {:.6}, {:.6}]\n",
        onb.v[0], onb.v[1], onb.v[2]
    ));
    out.push_str(&format!(
        "  u (従法線):   [{:.6}, {:.6}, {:.6}]\n",
        onb.u[0], onb.u[1], onb.u[2]
    ));
    out.push_str("\n  直交性の確認（≈ 0 であれば OK）：\n");
    out.push_str(&format!("  |u·v| = {:.2e}\n", dot(onb.u, onb.v).abs()));
    out.push_str(&format!("  |u·w| = {:.2e}\n", dot(onb.u, onb.w).abs()));
    out.push_str(&format!("  |v·w| = {:.2e}\n", dot(onb.v, onb.w).abs()));
    out.push_str("\n  単位長さの確認（≈ 1 であれば OK）：\n");
    out.push_str(&format!("  |u| = {:.12}\n", dot(onb.u, onb.u).sqrt()));
    out.push_str(&format!("  |v| = {:.12}\n", dot(onb.v, onb.v).sqrt()));
    out.push_str(&format!("  |w| = {:.12}\n\n", dot(onb.w, onb.w).sqrt()));

    let n2 = 1_000_000usize;

    // --- Section 2: MC estimate via ONB-transformed cosine direction ---------
    out.push_str("=== ONB 変換後の cosine 方向で ∫_hemi cos³θ dA を推定 ===\n");
    out.push_str("    （θ は法線 w に対する角度）\n");
    out.push_str(&format!("  π/2 (解析解) = {:.12}\n", PI / 2.0));

    let onb2 = Onb::build_from_w(n);
    let mut rng2 = XorShift64::new(0x1a2b_3c4d_0001);
    let sum2: f64 = (0..n2)
        .map(|_| {
            let local_dir = random_cosine_direction(&mut rng2);
            let world_dir = onb2.local_vec(local_dir);
            // cos(θ) relative to normal w = dot(world_dir, w)
            let cos_theta = dot(world_dir, onb2.w);
            // p = cos(θ)/π (cosine PDF), f = cos³θ → f/p = π·cos²θ
            cos_theta * cos_theta * cos_theta / (cos_theta / PI)
        })
        .sum();
    out.push_str(&format!(
        "  推定値 (N={}) = {:.12}\n",
        n2,
        sum2 / n2 as f64
    ));

    out
}

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

    /// The three ONB vectors are mutually orthogonal (dot products ≈ 0).
    #[test]
    fn onb_vectors_are_orthogonal() {
        let normals: &[[f64; 3]] = &[
            normalize([1.0, 0.0, 0.0]),
            normalize([0.0, 1.0, 0.0]),
            normalize([0.0, 0.0, 1.0]),
            normalize([1.0, 2.0, 3.0]),
            normalize([-1.0, 0.5, 0.2]),
        ];
        for &n in normals {
            let onb = Onb::build_from_w(n);
            assert!(dot(onb.u, onb.v).abs() < 1e-12);
            assert!(dot(onb.u, onb.w).abs() < 1e-12);
            assert!(dot(onb.v, onb.w).abs() < 1e-12);
        }
    }

    /// The three ONB vectors are unit length.
    #[test]
    fn onb_vectors_are_unit_length() {
        let onb = Onb::build_from_w(normalize([1.0, 2.0, 3.0]));
        assert!((dot(onb.u, onb.u).sqrt() - 1.0).abs() < 1e-12);
        assert!((dot(onb.v, onb.v).sqrt() - 1.0).abs() < 1e-12);
        assert!((dot(onb.w, onb.w).sqrt() - 1.0).abs() < 1e-12);
    }

    /// ONB-transformed cosine sampling estimates ∫_hemi cos³θ dA ≈ π/2.
    #[test]
    fn onb_cosine_sampling_estimates_pi_over_2() {
        let onb = Onb::build_from_w(normalize([1.0, 2.0, 3.0]));
        let n = 500_000;
        let mut rng = XorShift64::new(0xdead_beef_0001);
        let sum: f64 = (0..n)
            .map(|_| {
                let local_dir = random_cosine_direction(&mut rng);
                let world_dir = onb.local_vec(local_dir);
                let cos_theta = dot(world_dir, onb.w);
                cos_theta * cos_theta * cos_theta / (cos_theta / PI)
            })
            .sum();
        let est = sum / n as f64;
        let exact = PI / 2.0;
        assert!(
            (est - exact).abs() < 0.005,
            "ONB cosine estimate {est:.6} too far from π/2={exact:.6}"
        );
    }
}

ブラウザ上での実行結果

直交性チェックは内積の絶対値が浮動小数点誤差（≈ ${10}^{-16}$）であることを確認します。推定値は解析解 $\pi /2\approx 1.5708$ に近い値を示します。

まとめ

• 法線 $\mathbf{n}$ から ONB を構築する手順を確立した：$\mathbf{w}=\text{normalize}(\mathbf{n})$，$\mathbf{v}=\text{normalize}(\mathbf{w}\times \mathbf{a})$，$\mathbf{u}=\mathbf{w}\times \mathbf{v}$。
• local(a, b, c) で $z$ 軸周りのローカル方向をワールド座標に変換できる。
• ONB 変換後の cosine 方向サンプリングでも ${\int }_{\text{hemi}}{\cos }^3\theta d\omega =\pi /2$ が得られることを数値確認した。

次章では PDF クラスを導入し，複数の PDF を混合してサンプリングする仕組みを作ります。