QBVHを実装した - 明日ではないから

空間構造の一種であるQBVH(Quad-tree Bounding Volume Hierarchy)を実装した。
元の論文はhttp://www.uni-ulm.de/fileadmin/website_uni_ulm/iui.inst.100/institut/Papers/QBVH.pdf
QBVHはレイトレーサの交差判定を高速化する。
以前syoyoさんがイイヨとおっしゃっていたので実装してみた。

この論文で語られるQBVHの要点は以下。
１．４分木構造
 木構造においてひとつのノードに４つの子ノードを格納する。
２．ノードのAABB(バウンド)を親が持つ
一つ一つのノードがそのAABBを持つのではなくその親が４ついっぺんにもつ
３．SIMDによる判定
４つの子ノードのAABBの交差判定をSIMD演算を利用していっぺんに行う
４．面情報を格納するノードは作らない
面情報を格納するノードは実際には作らない。その代わり、ノードの子が枝ノードなのか面なのかをインデクッスの一部をフラグとし格納する。

■実装

普通BVHといえば２分木構造を用いるがQBVHでは４分木構造を用いる。これはSIMD(インテルCPUではSSE)ではfloat*4の演算を一度にすることができるため、４分木であることが最適であるからだ。
木構造のあるノードのAABBはそのノード自身が持つのではなく、その親のノードが持つ。ゆえにひとつのノードはその子の４つのAABBを持つ。
枝ノードはSSEでの実装では以下のような構造体で記述できる。

struct SIMDBVHNode{
    __m128 bboxes[2][3];//4 float min-max xyz
    size_t children[4]; //4 children
    int axis_top;       //top axis
    int axis_left;      //left axis
    int axis_right;     //right axis
    int reserved;       //padding 
};

SIMDによる判定はAABBに対して行う。
今回はその先の面の交差には用いない。（論文では実装はしてるようなのだ・・・けど・・・the original triangle layout is left untouched.となってるがはてはて・・・）
SSEでは以下のような関数となる。

inline int test_AABB(
    const __m128 bboxes[2][3],  //4boxes : min-max[2] of xyz[3] of boxes[4]
    const __m128 org[3],        //ray origin
    const __m128 idir[3],       //ray inveresed direction
    const int sign[3],          //ray xyz direction -> +:0,-:1
    __m128 tmin, __m128 tmax    //ray range tmin-tmax 
)
{
    // x coordinate
    tmin = _mm_max_ps(
        tmin,
        _mm_mul_ps(_mm_sub_ps(bboxes[sign[0]][0],org[0]), idir[0])
    );
    tmax = _mm_min_ps(
        tmax,
        _mm_mul_ps(_mm_sub_ps(bboxes[1 - sign[0]][0], org[0]), idir[0])
    );

    // y coordinate
    tmin = _mm_max_ps(
        tmin,
        _mm_mul_ps(_mm_sub_ps(bboxes[sign[1]][1],org[1]), idir[1])
    );
    tmax = _mm_min_ps(
        tmax,
        _mm_mul_ps(_mm_sub_ps(bboxes[1 - sign[1]][1], org[1]), idir[1])
    );

    // z coordinate
    tmin = _mm_max_ps(
        tmin,
        _mm_mul_ps(_mm_sub_ps(bboxes[sign[2]][2],org[2]), idir[2])
    );
    tmax = _mm_min_ps(
        tmax,
        _mm_mul_ps(_mm_sub_ps(bboxes[1 - sign[2]][2], org[2]), idir[2])
    );

    return _mm_movemask_ps(_mm_cmpge_ps(tmax, tmin));//tmin<tmaxとなれば交差
}

これで各枝ノードでの判定は終わり、面との交差判定となる。面と枝ノードはインデックスの一部をフラグとする。
元の論文では以下のようにレイアウトしている。

今回の実装では最上位ビットのみフラグとし、面の数はインデックスに含めないこととした。その代わり、面の最後にNULLを置くレイアウトとした。

以上のように実装した。

■性能
QVBHの性能を以下に示す。交差判定をSSEを使うもの(SIMD)とそうでないもの(SISD)をそれぞれ実装した。またQBVHではfloat*4なのでdouble*2のBBVHも実装してみた。
比較対象として、通常のBVH（ノード自身がAABBを持つ）,BIH,KD-Tree,Octreeを用意した。
構築時間

交差判定時間(1000000rays)

消費メモリー

グラフも示す。

QBVHを含め、BVH系構築時間は長めになっている。これは今回の実装が構築するとき空間で振り分けたのではなく、sortで並べることによるものだ。
QBVHでは他の構造と比べ消費メモリーも少なく、肝心の交差判定も非常に速いことがわかる。自分の環境(Core2Duo T8100 2.1GHzh)で約300000rays/sec程度処理できている。
SIMDとSISDではQBVHで約2倍、BBVHで約1.5倍の高速化が実現できていることがわかる。
QBVHではfloat精度なのだが、今回のQBVHではAABBの交差判定のみSIMD化を行った。そのため、最終的な精度は面での交差判定が行うため、こちらがdouble精度の場合、AABBには適切なマージン（最低でもstd::numeric_limits::epsilon()*2）分膨らませてあげれば精度については問題なくなるはずだ。（マージンより狭い範囲に面が集中するようなメッシュだと交差判定が遅くなることはあるがまあ到底考えなくていいだろう）

■結論
QBVHのSIMD化はメッシュの空間構造のみ手を入れればいいため、既存のシステムにマージしやすい。
交差判定スピードは非常に速いため、SSEのようなSIMD命令が使える環境では実装してみる価値はあるはずだ。

■追記
以前実装したときそんなに速くないと吹聴しましたが、参考にした実装(lux render-qbvhaccel.cpp)が間違っていました。
具体的にはAABBの交差判定後、交差するノードをスタックに積む順番が逆になっていました。

■ソース等
以下のURLからソース等の一部をダウンロードできるようにしておきました。参考になるかどうかはわかりませんが。
http://kaze-renderer.googlecode.com/files/acc_test20090925.zip
また全ソースは以下のURLにホスティングしてあります。
http://code.google.com/p/kaze-renderer/