Hashed-Octree(ハッシュ化八分木)による空間分割を実装した

そろそろレイトレがしたい。具体的に言えばフォトンマッピングでGIを実装したい。

だが、フォトンマッピングを実装するには現状あまりにも足りないものが多すぎる。

その最たるものが、高速な衝突判定アルゴリズムだ。

シンプルに設計されたレイトレ専用のアプリケーションなら、レイの衝突判定は、シーン中の全てのオブジェクトに対して総当たり的にやってもいいだろう。メッシュも、プリミティブな形状に限定すれば、高速に処理できる。

しかし、現在私が作ってるのは曲がりなりにもリアルタイムレンダリングエンジンだ。シーンの大きさや配置されるオブジェクトの数は固定ではない。現実的な処理時間で動作させなくてはならない。また、メッシュも多くのポリゴンを有する複雑なモデルが存在するかもしれない。
そのためには、レイの衝突判定において、衝突しないことが自明なオブジェクトに対しては処理をスキップするような仕組みが必要である。

このような仕組みで最も直感的なものは、ユニフォームグリッドによる空間分割だろう。
シーン全体を一様なグリッドで分割し、所属するグリッド内のオブジェクトのみで衝突判定を行う。
ただし、この分割法には多くの問題がある。まず、グリッドの横断が想定されていない。空間の分割の粒度をどれくらいにするにせよ、位置、サイズによっては必ず複数のグリッドにまたがるオブジェクトが出現するはずである。基本的には一つのオブジェクトが複数のグリッドから参照される形になるだろう。また、メモリ効率も悪い。一般に、空間というものは疎なもので、基本的には何もない。そして、偏りがあるものである。例えば、空中にオブジェクトがあることはまれだが、あるテーブルに着目した時、そのテーブルの上には多くのオブジェクトがあることが期待される。このような分布の空間において、一様なグリッドで分割するのはあまりにも非効率である。

このようなケースにおける最適な空間分割法として、八分木(Octree)がある。
まず、ルートの空間を定義する。次に、X-Y-Zについてそれぞれ2分割する。分割後の空間を空間レベル1(L=1)として、同様の操作を任意の空間レベルまで行う。

ユニフォームグリッドとの大きな違いは、空間レベルという概念があること。「グリッドをまたぐ」というのは存在せず、仮にグリッドをまたぐ場合は、それは1ランク上の「親空間」と呼ばれる、より大きな分割単位の空間に所属する仕組みになる。また、グリッドの番号にはモートンオーダーを用いることで、ビット演算を駆使してあるオブジェクトの属するすべての空間レベルのグリッド番号を取得することも可能である。

この辺の細かい説明は、すでに必要十分な説明を記載しているサイトがいくつもあるので省略することにする。
その８ 4分木空間分割を最適化する！

ちなみに、アルゴリズム的には木構造であるが、モートンオーダーに対して空間レベルの等比数列の総和のオフセットを付加することで、任意の空間レベルの任意のグリッドのオーダーをハッシュ化することができる。
上記のサイトでもそれは行われているが、実装が配列になっているため、結局常に （8^L - 1 / 7）の要素数の配列の動的確保を要求していて効率が悪い。 std::unordred_mapのようなハッシュマップを使えば、疎な空間における効率的なメモリ確保が実現できるだろう。

実装にあたって、八分木分割後のビジュアライザーが必要だと感じたので、そちらもついでに作った。
空間座標から所属するモートンオーダーの算出方法は上記サイト様に記載のある通りだが、モートンオーダーから空間座標を逆算する計算が不明だったので、ここに記しておく。

AABB CalculateOctreeBoxAABBFromMortonNumber(uint32_t number) {
    int level = 0;
   // ハッシュ値から所属する最小空間のモートンオーダーに変換
    while (number >= std::pow(8, level)) {
        number -= std::pow(8, level);
        level++;
    }

    uint32_t s = 0;

    for (int i = level; i > 0; i--) {
        s = s | (number >> (3 * i - 2 - i) & (1 << i - 1));
    }
    uint32_t x = s;

    s = 0;
    for (int i = level; i > 0; i--) {
        s = s | (number >> (3 * i - 1 - i) & (1 << i - 1));
    }
    uint32_t y = s;

    s = 0;
    for (int i = level; i > 0; i--) {
        s = s | (number >> (3 * i - i) & (1 << i - 1));
    }
    uint32_t z = s;

    // _rootAABB.size: ルート空間のサイズ。空間レベルで割って所属する空間レベルの分割サイズを求める
    Size3D boxSize = _rootAABB.size() / (1 << level);
    // _rootAABB.bpos: ルート空間の開始座標
    Vector3D bpos = Vector3D(x * boxSize.w, y * boxSize.h, z * boxSize.d) + _rootAABB.bpos;
    // 所属するAABB
    return AABB(bpos, Vector3D(bpos.x + boxSize.w, bpos.y + boxSize.h, bpos.z + boxSize.d));
}

これでハッシュ値から、所属するAABBを求めることができる。
戦略的には3bitごとに区切られたX-Y-Zのビットに着目して、特定の次元について切り詰める。例えば、110111001というモートンオーダーをXについて切り詰めると、011だ。今回の場合は空間レベルは最大8に設定したので、24桁のビットがある。もう少し効率的な方法がある気もするが....

実際に計算して可視化してみた。
各オブジェクトが所属する最小八分木空間を計算した後は、それと同じサイズのボックスメッシュを生成し、ワイヤーフレーム描画してそれっぽく見せている。

同じオブジェクトでも、位置によってはグリッドをまたいでしまうため、より大きなAABB（より高いレベルの空間）に所属している様子が分かる。

これで空間内の衝突判定の準備は整った。次は、レイが空間を通過するときの衝突リストの作成だが...いまいち効率的な方法が思いついていない。