3. 包含立体階層

Ray Tracing: The Next Week (v3.2.3): 3 Bounding Volume Hierarchies / 2.3 包含立体階層

レイ・オブジェクトのヒット判定はレイトレーサ最大のボトルネックであり，素朴な実装ではオブジェクト数 $N$ に対して $O (N)$ かかります。同じシーンに対して数百万本のレイを発射するという反復構造を活用すれば，二分探索の発想で $O (\log N)$ に落とせます。包含立体階層（Bounding Volume Hierarchy, BVH）はこれを実現するデータ構造で，本章のレイトレーサで最も込み入った部分でもあります。

新規の道具立ては次の 4 つです。

Aabb：軸平行境界ボックス（Axis-Aligned Bounding Box）と，それとレイの衝突判定（Andrew Kensler 最適化版）
surrounding_box：2 つの AABB を含む最小の AABB
Hittable トレイトに bounding_box メソッドを既存の r1XX クレートを壊さないよう加算的に追加
BvhNode：自身が Hittable でもある BVH ノード

基本的なアイデア

オブジェクトの集合の包含立体（bounding volume）とは，集合内の全オブジェクトを完全に包む立体のことです。たとえば 10 個の物体の包含球を計算しておけば，包含球と交わらないレイは 10 個の物体のどれとも交わりません。逆にレイが包含球と交わったときは，10 個のうちのどれか，あるいは複数と交わる可能性があります。コードの形は次のようになります。

text

レイが包含立体と交わる
    → 内側の物体それぞれと交わるか調べる
そうでない
    → 衝突なしを返す

ここで重要なのは，画面や空間ではなくオブジェクト集合を分割している点です。各オブジェクトはちょうど 1 つの包含立体に属しますが，包含立体同士は重なって構いません。

包含立体の階層

劣線形^[sublinear] $O (\log N)$ など，線形 $O (N)$ より小さいオーダーのこと。対義語は超線形（ $O (N^{2})$ など）。性を達成するには，包含立体を階層化します。たとえばオブジェクト集合を赤・青の 2 グループに分け，それぞれに矩形包含立体を被せると，両者を含む親の包含立体が定まります。

text

親と交わる
    hit_blue  = 青のグループ内の物体と交わる
    hit_red   = 赤のグループ内の物体と交わる
    どちらかでもヒットなら衝突情報（より近い方）を返す
そうでない
    衝突なしを返す

赤と青は親の中で重なっていてもよく，順序関係もありません。木構造としてみたとき，左右の子に意味的な区別はなく単に「親の中にいる」だけです。各レイは木をルートからたどり，親の AABB と交わらないサブツリーを丸ごと枝刈りします。葉に達すれば最終的に通常のヒット判定を実行します。

軸平行包含直方体（AABB）

包含立体には種々の選択肢がありますが，多くのモデルで実用的に優れているのが軸平行包含直方体（axis-aligned bounding box, AABB）です。AABB は $n$ 次元の AABB を $n$ 個の軸平行な区間（スラブ）の共通部分として捉えるのが標準的な見方です。2 次元なら $x$ の区間と $y$ の区間の共通部分が長方形を構成します。

レイ $P (t) = A + t b$ がスラブ $x \in [x_{0}, x_{1}]$ と交わるパラメータ $t$ は

t_{x 0} = \frac{x_{0} - A_{x}}{b_{x}}, t_{x 1} = \frac{x_{1} - A_{x}}{b_{x}}

で求まります。3 次元の AABB と交わる条件は， $x$ / $y$ / $z$ それぞれで得た $t$ の区間の共通部分が空でないこと，すなわち

max (t_{x 0}, t_{y 0}, t_{z 0}, t_{min}) < min (t_{x 1}, t_{y 1}, t_{z 1}, t_{max})

です。区間の重なりさえ判定できればよく，衝突点や法線は不要なので非常に高速に書けます。

注意点が 2 つあります。

レイが軸方向に逆向き（ $b_{x} < 0$ など）に進むと， $t_{x 0} > t_{x 1}$ という反転区間が出てくる。
$b_{x} = 0$ で 0 除算となる。IEEE 754 では結果が $\pm \infty$ になり，区間外のレイは両端とも同符号の無限大になる。NaN も出ることがある。

これらを fmin/fmax での吸収＋反転時の swap で扱うのが定石です。

AABB の衝突判定

Pixar の Andrew Kensler が提案した最適化版を採用します。除算を 1 回にまとめ，分岐は反転検出の swap だけです。

common/src/aabb.rs

rust

pub struct Aabb {
    pub minimum: Point3,
    pub maximum: Point3,
}

impl Aabb {
    pub fn new(a: Point3, b: Point3) -> Self {
        Aabb { minimum: a, maximum: b }
    }
    pub fn min(&self) -> Point3 { self.minimum }
    pub fn max(&self) -> Point3 { self.maximum }

    pub fn hit(&self, r: &Ray, mut t_min: f64, mut t_max: f64) -> bool {
        for a in 0..3 {
            let inv_d = 1.0 / r.direction()[a];
            let mut t0 = (self.minimum[a] - r.origin()[a]) * inv_d;
            let mut t1 = (self.maximum[a] - r.origin()[a]) * inv_d;
            if inv_d < 0.0 {
                std::mem::swap(&mut t0, &mut t1);
            }
            if t0 > t_min { t_min = t0; }
            if t1 < t_max { t_max = t1; }
            if t_max <= t_min {
                return false;
            }
        }
        true
    }
}

3 軸ループで区間 $[t_{min}, t_{max}]$ を逐次更新し，空区間になった時点で false を返します。

surrounding_box は 2 つの AABB を包む最小の AABB を返す自由関数です。BVH 構築時に親ノードの AABB を計算するのに使います。

common/src/aabb.rs

rust

pub fn surrounding_box(box0: Aabb, box1: Aabb) -> Aabb {
    let small = Point3::new(
        box0.min().x().min(box1.min().x()),
        box0.min().y().min(box1.min().y()),
        box0.min().z().min(box1.min().z()),
    );
    let big = Point3::new(
        box0.max().x().max(box1.max().x()),
        box0.max().y().max(box1.max().y()),
        box0.max().z().max(box1.max().z()),
    );
    Aabb::new(small, big)
}

`Hittable` トレイトを拡張する

任意の Hittable から AABB を取り出せるよう，トレイトに bounding_box メソッドを追加します。動くオブジェクトは時刻区間 $[t_{0}, t_{1}]$ にわたって全時刻の自身を含むAABB を返さねばならないので，引数として time0 / time1 を取ります。

common/src/hittable.rs

rust

pub trait Hittable {
    fn hit(&self, r: &Ray, t_min: f64, t_max: f64) -> Option<HitRecord>;

    /// AABB が定義できないオブジェクト（無限平面など）では `None` を返す。
    /// デフォルト実装は `None` を返すので，BVH に投入しないオブジェクト
    /// （第 1 編で実装した `Sphere::new` 経由の球など）は何もしなくてよい。
    fn bounding_box(&self, _time0: f64, _time1: f64) -> Option<Aabb> {
        None
    }
}

C++ と Rust の違い

C++ 原典では bounding_box を純粋仮想関数として宣言し，全派生クラスに実装を強制します。これに対して本稿では Option<Aabb> を返すデフォルト実装つきメソッドとしました。理由は次の 2 つです。

第 1 編の r1XX クレートで定義した Hittable 実装（Sphere::new 経由の素朴な球など）がそのままビルドできるように，加算的拡張に留めたい。
「AABB が存在しない」をブール値＋出力引数ではなく Option で表現するのが Rust の慣用。

無限平面のように真に AABB を持たないオブジェクトもあるため，トレイト全体としては Option<Aabb> を返すのが自然な型でもあります。

Sphere，MovingSphere，HittableList の 3 つはこのデフォルトを上書きして実装します。Sphere は中心 $\pm (r, r, r)$ の AABB です。

common/src/sphere.rs

rust

impl Hittable for Sphere {
    fn hit(...) -> Option<HitRecord> { /* 既存 */ }

    fn bounding_box(&self, _time0: f64, _time1: f64) -> Option<Aabb> {
        let r = self.radius.abs();
        let radius_vec = Vec3::new(r, r, r);
        Some(Aabb::new(self.center - radius_vec, self.center + radius_vec))
    }
}

MovingSphere は時刻 time0 の球の AABB と時刻 time1 の球の AABB を surrounding_box で包みます。

common/src/moving_sphere.rs

rust

fn bounding_box(&self, time0: f64, time1: f64) -> Option<Aabb> {
    let radius_vec = Vec3::new(self.radius, self.radius, self.radius);
    let box0 = Aabb::new(
        self.center(time0) - radius_vec,
        self.center(time0) + radius_vec,
    );
    let box1 = Aabb::new(
        self.center(time1) - radius_vec,
        self.center(time1) + radius_vec,
    );
    Some(surrounding_box(box0, box1))
}

HittableList は子の AABB を順に surrounding_box で畳み込んだものを返します。空リストや子が None を返す場合は None を伝播します。

common/src/hittable_list.rs

rust

fn bounding_box(&self, time0: f64, time1: f64) -> Option<Aabb> {
    if self.objects.is_empty() {
        return None;
    }
    let mut output: Option<Aabb> = None;
    for object in &self.objects {
        let temp = object.bounding_box(time0, time1)?;
        output = Some(match output {
            None => temp,
            Some(prev) => surrounding_box(prev, temp),
        });
    }
    output
}

BVH ノード

BVH も Hittable の一種です。それ自身は他の Hittable のコンテナに過ぎませんが，「このレイと交わるか？」に答えられます。原典に倣い，木とノードを別クラスに分けずノード 1 つだけで表します（ルートも単なるノード）。

common/src/bvh.rs

rust

pub struct BvhNode {
    left: Arc<dyn Hittable>,
    right: Arc<dyn Hittable>,
    bbox: Aabb,
}

子ポインタの型は具体型でなく Arc<dyn Hittable> です。BvhNode 自身も Sphere も MovingSphere も入ります。Arc を使うのは BVH 構築時に同じ葉オブジェクトを左右両方の子に複製する場合（要素数 1 のサブツリー）があるためです。

hit は素直で，親の AABB に外れたら子に降りないだけです。

common/src/bvh.rs

rust

impl Hittable for BvhNode {
    fn hit(&self, r: &Ray, t_min: f64, t_max: f64) -> Option<HitRecord> {
        if !self.bbox.hit(r, t_min, t_max) {
            return None;
        }
        let hit_left = self.left.hit(r, t_min, t_max);
        // 左でヒットした t を右の上限に使うことで，
        // 右側の探索範囲を更に絞れる。
        let new_t_max = match &hit_left {
            Some(rec) => rec.t,
            None => t_max,
        };
        let hit_right = self.right.hit(r, t_min, new_t_max);
        hit_right.or(hit_left)
    }

    fn bounding_box(&self, _time0: f64, _time1: f64) -> Option<Aabb> {
        Some(self.bbox)
    }
}

hit_right.or(hit_left) は，右でヒットしたらそれ（より近い側）を返し，さもなければ左の結果を返す慣用句です。new_t_max を縮めることで，右側探索の早期打ち切りも自動で機能します。

BVH の分割

効率化用データ構造でいちばん複雑なのは構築です。BVH 構築は BvhNode::new 内で行います。子の集合をうまく分けられれば hit は正しく動くので，正しさだけなら適当な分割で構いません。当然，子の AABB が親より小さいほうが速いものの，それは性能の話です。

ここでは原典に従い，シンプルさと性能の妥協点として次のアルゴリズムを採用します。

ランダムに 1 軸を選ぶ
その軸に沿ってオブジェクトをソートする
中央で 2 分割し，左右の子それぞれに対して再帰

要素数が 1 のときは左右の子に同じオブジェクトを入れて再帰を打ち切ります（Arc::clone のおかげで重複保持はゼロコストに近い）。要素数 2 のときは比較器で 1 回の swap だけ，3 以上で sort_by＋分割します。

common/src/bvh.rs

rust

fn build(
    objects: &mut Vec<Arc<dyn Hittable>>,
    start: usize, end: usize,
    time0: f64, time1: f64,
) -> Self {
    let axis = random_int(0, 2) as usize;
    let span = end - start;

    let (left, right): (Arc<dyn Hittable>, Arc<dyn Hittable>) = if span == 1 {
        (objects[start].clone(), objects[start].clone())
    } else if span == 2 {
        if box_compare(objects[start].as_ref(),
                       objects[start + 1].as_ref(), axis) == Ordering::Less {
            (objects[start].clone(), objects[start + 1].clone())
        } else {
            (objects[start + 1].clone(), objects[start].clone())
        }
    } else {
        objects[start..end]
            .sort_by(|a, b| box_compare(a.as_ref(), b.as_ref(), axis));
        let mid = start + span / 2;
        let l = Arc::new(Self::build(objects, start, mid, time0, time1));
        let r = Arc::new(Self::build(objects, mid, end, time0, time1));
        (l as Arc<dyn Hittable>, r as Arc<dyn Hittable>)
    };

    let box_left = left.bounding_box(time0, time1)
        .expect("No bounding box in BvhNode constructor");
    let box_right = right.bounding_box(time0, time1)
        .expect("No bounding box in BvhNode constructor");
    let bbox = surrounding_box(box_left, box_right);

    BvhNode { left, right, bbox }
}

公開 API は BvhNode::new と，HittableList を直接受けて消費する BvhNode::from_list の 2 つです。from_list は Box<dyn Hittable> を Arc<dyn Hittable> に変換する糖衣です（標準ライブラリの impl<T: ?Sized> From<Box<T>> for Arc<T> を使用）。

common/src/bvh.rs

rust

pub fn from_list(list: HittableList, time0: f64, time1: f64) -> Self {
    let objects: Vec<Arc<dyn Hittable>> =
        list.objects.into_iter().map(Arc::<dyn Hittable>::from).collect();
    Self::new(objects, time0, time1)
}

random_int(min, max) は [min, max] 上の一様乱数整数を返すヘルパで，common::utils に追加しました。random_double_range(min as f64, (max + 1) as f64).floor() as i32 で実装しています。

AABB の比較関数

ソートに使う比較器です。指定軸の AABB の最小角を比較するだけで，原典どおり「軸引数つきの汎用関数」を 1 つ書けば足ります（C++ 版のように軸ごとに 3 つの単項関数に分けるのは std::sort の関数ポインタ受け渡し都合で，Rust ではクロージャを使えるため不要）。

common/src/bvh.rs

rust

fn box_compare(a: &dyn Hittable, b: &dyn Hittable, axis: usize) -> Ordering {
    let box_a = a.bounding_box(0.0, 0.0)
        .expect("No bounding box in box_compare");
    let box_b = b.bounding_box(0.0, 0.0)
        .expect("No bounding box in box_compare");
    box_a.min()[axis]
        .partial_cmp(&box_b.min()[axis])
        .unwrap_or(Ordering::Equal)
}

f64 には Ord が実装されないため partial_cmp で取り，NaN 由来の None は Equal に潰します（健全な AABB ならば発生しません）。

common/src/bvh.rs

rust

use std::cmp::Ordering;
use std::sync::Arc;

use crate::aabb::{Aabb, surrounding_box};
use crate::hittable::{HitRecord, Hittable};
use crate::hittable_list::HittableList;
use crate::ray::Ray;
use crate::utils::random_int;

pub struct BvhNode {
    left: Arc<dyn Hittable>,
    right: Arc<dyn Hittable>,
    bbox: Aabb,
}

impl BvhNode {
    pub fn new(mut objects: Vec<Arc<dyn Hittable>>, time0: f64, time1: f64) -> Self {
        let len = objects.len();
        Self::build(&mut objects, 0, len, time0, time1)
    }

    pub fn from_list(list: HittableList, time0: f64, time1: f64) -> Self {
        let objects: Vec<Arc<dyn Hittable>> =
            list.objects.into_iter().map(Arc::<dyn Hittable>::from).collect();
        Self::new(objects, time0, time1)
    }

    fn build(
        objects: &mut Vec<Arc<dyn Hittable>>,
        start: usize, end: usize,
        time0: f64, time1: f64,
    ) -> Self {
        let axis = random_int(0, 2) as usize;
        let span = end - start;

        let (left, right): (Arc<dyn Hittable>, Arc<dyn Hittable>) = if span == 1 {
            (objects[start].clone(), objects[start].clone())
        } else if span == 2 {
            if box_compare(objects[start].as_ref(),
                           objects[start + 1].as_ref(), axis) == Ordering::Less {
                (objects[start].clone(), objects[start + 1].clone())
            } else {
                (objects[start + 1].clone(), objects[start].clone())
            }
        } else {
            objects[start..end]
                .sort_by(|a, b| box_compare(a.as_ref(), b.as_ref(), axis));
            let mid = start + span / 2;
            let l = Arc::new(Self::build(objects, start, mid, time0, time1))
                as Arc<dyn Hittable>;
            let r = Arc::new(Self::build(objects, mid, end, time0, time1))
                as Arc<dyn Hittable>;
            (l, r)
        };

        let box_left = left.bounding_box(time0, time1)
            .expect("No bounding box in BvhNode constructor");
        let box_right = right.bounding_box(time0, time1)
            .expect("No bounding box in BvhNode constructor");
        let bbox = surrounding_box(box_left, box_right);

        BvhNode { left, right, bbox }
    }
}

impl Hittable for BvhNode {
    fn hit(&self, r: &Ray, t_min: f64, t_max: f64) -> Option<HitRecord> {
        if !self.bbox.hit(r, t_min, t_max) {
            return None;
        }
        let hit_left = self.left.hit(r, t_min, t_max);
        let new_t_max = match &hit_left {
            Some(rec) => rec.t,
            None => t_max,
        };
        let hit_right = self.right.hit(r, t_min, new_t_max);
        hit_right.or(hit_left)
    }

    fn bounding_box(&self, _time0: f64, _time1: f64) -> Option<Aabb> {
        Some(self.bbox)
    }
}

fn box_compare(a: &dyn Hittable, b: &dyn Hittable, axis: usize) -> Ordering {
    let box_a = a.bounding_box(0.0, 0.0)
        .expect("No bounding box in box_compare");
    let box_b = b.bounding_box(0.0, 0.0)
        .expect("No bounding box in box_compare");
    box_a.min()[axis]
        .partial_cmp(&box_b.min()[axis])
        .unwrap_or(Ordering::Equal)
}

デモ：r203-bvh

r203-bvh クレートでは 2.2 と同じ「跳ねる球」シーンを，最後に BvhNode::from_list で BVH に包んでから渡します。

r203-bvh/src/lib.rs

rust

let world = BvhNode::from_list(random_scene(), 0.0, 1.0);

シーンの中身もカメラパラメータも 2.2 と同じです。違うのはヒット判定の計算量だけで，画面サイズ 300×169 を維持したままサンプル数を 30 → 100 に上げています。素朴な線形走査では 100 サンプルは現実的でない時間がかかりますが，BVH なら（おおよそ $\log_{2} (N)$ 倍の高速化が得られるため）十分許容できる時間で描画できます。

ページ冒頭の WASM デモで確認できます。出力画像は 2.2 とほぼ同等の傾向ですが，サンプル数が増えたぶんノイズが目に見えて減っています。

まとめ

Aabb と Andrew Kensler 最適化版の hit，および surrounding_box を common/ に追加した。
Hittable トレイトに bounding_box メソッドをデフォルト実装つきで加えた。これにより第 1 編の r1XX クレートはそのままビルドできる。
Sphere・MovingSphere・HittableList で bounding_box を実装した。
BvhNode を実装し，乱数軸選択＋ソート＋中央分割で再帰的に構築した。
r203-bvh で実用効果（同じシーンを高サンプル数で描画）を確認した。

これで第 2 編の最も難しい仕掛けは終わりです。次章ではテクスチャの導入に進み，BVH の上にいよいよ「絵作り」の自由度を載せていきます。

$O (\log N)$ など，線形 $O (N)$ より小さいオーダーのこと。対義語は超線形（ $O (N^{2})$ など）。 ↩︎

3. 包含立体階層 ​

基本的なアイデア ​

包含立体の階層 ​

軸平行包含直方体（AABB） ​

AABB の衝突判定 ​

Hittable トレイトを拡張する ​

C++ と Rust の違い ​

BVH ノード ​

BVH の分割 ​

AABB の比較関数 ​

デモ：r203-bvh ​

まとめ ​

3. 包含立体階層

基本的なアイデア

包含立体の階層

軸平行包含直方体（AABB）

AABB の衝突判定

`Hittable` トレイトを拡張する

C++ と Rust の違い

BVH ノード

BVH の分割

AABB の比較関数

デモ：r203-bvh

まとめ