8分木空間でのレイトレースのコリジョンリスト収集処理を見直して10倍速くした

　フォトンマッピングの実装で、フォトンを1万～10万くらいバラまく処理が死ぬほど重かったので、8分木空間におけるレイトレースの処理を見直した。

　8分木はもとより衝突判定を高速にするための空間分割アルゴリズムだが、レイトレースをする際に、衝突判定をすべき空間の収集に最適化の余地があったのでいろいろコードを弄っていたら、10倍速くなった。

　まず簡潔に8分木の分割のロジックを説明する。が、簡略化のために2次元で考える。つまり、4分木として考える。
　例えば、分割レベル3の空間では、ルート空間、インデックス0～3の空間レベル1、インデックス0～15の空間レベル2の3つのレベルが作られる。この時、あるオブジェクトPが空間レベル2のインデックス3に属するとする。Pの衝突判定を行う場合、その空間と、そして親空間である、空間レベル1のインデックス0、ルート空間、この3つの空間に属するオブジェクトと比較をすればよい。これがN分木空間における衝突判定の高速化の原理である。ちなみに、空間レベルのグリッドをまたぐ場合は属する空間が親空間に移る。この辺の説明は以前の記事（Hashed-Octree(ハッシュ化八分木)による空間分割を実装した - My life accelerated）を参照されたい。
　

　さて、図を示すと上記のようなインデックスになる。特殊な配番に注意すること。これはシフト演算による処理の高速化のための、モートンオーダーを使った配番である。衝突処理の仕組みについては、前述の通りだが、ここで、レイトレースを行う場合の衝突処理を考える。原理的には、レイのポジションに対し、レイの方向ベクトルを加算し、そのたびにレイのポジションからN分木空間のインデックスを算出、衝突判定処理をすればよいが、それは効率的ではない。図を挙げてみよう。

　
　

　右上のオブジェクトはレイに衝突しないのは自明である。右下のオブジェクトは、グリッドをまたいでいるため、空間レベル2ではなく、空間レベル1のインデックス3に属している。私の従来の実装では、レイの衝突判定は、空間レベル2ベースでグリッドを移動する。インデックスのビジット順序は、空間レベル2の「0,2,3,9,12,14,15」、空間レベル1の「0,2,3」、そしてルート空間である。
　
　このような衝突リスト収集処理をエンジンから抜粋する。

std::set<uint32_t> GetColliderMortonList(SpaceOctree::OctreeFactoryBase* factory, Ray ray) {
    auto size = factory->GetMinBoxSize(); // 最大の空間レベルの分割サイズ
    auto rayForward = Vector3D(ray.dir.x * size.w, ray.dir.y * size.h, ray.dir.z * size.d); // レイが1ステップに進む距離
    auto rootAABB = factory->GetRootAABB(); // ルート空間
    
    _Vector3D<int16_t> grid = factory->CalculateGridCoordinate(ray.pos); // レイの初期位置から空間のグリッド座標を算出
    _Vector3D<int8_t> gridForward = _Vector3D<int8_t>( // レイ方向ベクトルの符号から1ステップにおけるグリッドの移動データを算出
        ray.dir.x >= 0.0f ? 1 : -1,
        ray.dir.y >= 0.0f ? 1 : -1,
        ray.dir.z >= 0.0f ? 1 : -1
        );

    Vector3D pos = Vector3D(grid.x * size.w, grid.y * size.h, grid.z * size.h) + rootAABB.bpos; // 初期位置
    _Vector3D<int16_t> nextGrid = grid;
    std::set<uint32_t> colliderList; // 衝突リスト（リストの中身は空間ハッシュ）

    while (rootAABB.Contains(pos)) {
　　 // グリッドから空間ハッシュ算出
        uint32_t number = SpaceOctree::Get3DMortonOrder(grid);

        // 空間ハッシュを、ルート空間まで遡って、衝突リストに格納していく（存在する場合のみ）
        for (int i = 0; i <= factory->GetSplitLevel(); i++) {
            uint32_t idx = static_cast<uint32_t>((number >> i * 3) + PrecomputedConstants::PowNumbers<8, 8>::Get(factory->GetSplitLevel() - i) / 7);
            if (factory->BoxExists(idx)) {
                colliderList.insert(idx);
            }
        }

        // 次のグリッド
        nextGrid = grid + gridForward;
        // 次の座標
        Vector3D nextpos = Vector3D(nextGrid.x * size.w, nextGrid.y * size.h, nextGrid.z * size.h) + rootAABB.bpos;

        // レイベクトルから、X方向、Y方向、Z方向のグリッドに到達する時のレイベクトルの係数を算出 
        float ax = ray.dir.x != 0.0f ? std::abs((nextpos.x - pos.x) / rayForward.x) : FLT_MAX;
        float ay = ray.dir.y != 0.0f ? std::abs((nextpos.y - pos.y) / rayForward.y) : FLT_MAX;
        float az = ray.dir.z != 0.0f ? std::abs((nextpos.z - pos.z) / rayForward.z) : FLT_MAX;

        // 最短で到達するグリッドの探索
        if (ax < ay && ax < az) {
            pos += rayForward * ax;
            grid.x += gridForward.x;
        }
        else if (ay < ax && ay < az) {
            pos += rayForward * ay;
            grid.y += gridForward.y;
        }
        else if (az < ax && az < ay) {
            pos += rayForward * az;
            grid.z += gridForward.z;
        }
        else {
            pos += rayForward;
            grid += gridForward;
        }
    }

    return colliderList;
}

　
　これはうまく動く。ただし、レイの1ステップにおける移動距離が、必ず最大空間レベル（例えば、上の4分木の例だとレベル2）の分割サイズにしかならない。これは、最適ではない。

　（少なくとも私の実装では）オブジェクトの存在しない空間は、ハッシュリストに登録されない。（あるオブジェクトを登録するとき、その親空間、更にその親空間...と遡って登録はする）。
　つまり、上記の例だと、ハッシュリストに登録されている空間は、空間レベル2の「5」、空間レベル1の「1,3」、そしてルート空間だけである。ならば、空間レベル1のインデックス0、左上の空間はまとめて無視できることがハッシュリストの構造から推測ができて、レイの距離は一気に9まで進めていいことがわかる。これを具体的にロジックで考えるならば、あるレイの点が属するすべての空間レベルのうち、「実際にハッシュリストに登録されている空間のレベル+1」の「分割サイズ」の距離だけレイを進めることができる。
　上記の例の場合、レイの始点がレベル2インデックス0空間とすると、レベル2のインデックス0及びレベル1のインデックス0はハッシュリストに存在せず、ルート空間（空間レベル0）のみが存在する。そのため、レイは空間レベル1の分割サイズだけレイを進めることができる。

　上記の疎空間におけるレイのステップ距離の最適化を施した後のコードが以下。

std::set<uint32_t> GetColliderMortonList(SpaceOctree::OctreeFactoryBase* factory, Ray ray) {
    auto min_size = factory->GetMinBoxSize();
    auto rootAABB = factory->GetRootAABB();
    
    _Vector3D<int16_t> grid = factory->CalculateGridCoordinate(ray.pos);
    _Vector3D<int16_t> gridForward = _Vector3D<int16_t>(
        ray.dir.x >= 0.0f ? 1 : -1,
        ray.dir.y >= 0.0f ? 1 : -1,
        ray.dir.z >= 0.0f ? 1 : -1
        );

    Vector3D pos = Vector3D(grid.x * min_size.w, grid.y * min_size.h, grid.z * min_size.h) + rootAABB.bpos;
    Vector3D next_pos = pos;
    std::set<uint32_t> colliderList;

    while (true) {
        uint32_t number = SpaceOctree::Get3DMortonOrder(grid);

        int exists_max_split_level = 0;
        for (int i = 0; i <= factory->GetSplitLevel(); i++) {
            int split_level = factory->GetSplitLevel() - i;
            uint32_t idx = static_cast<uint32_t>((number >> i * 3) + PrecomputedConstants::PowNumbers<8, 8>::Get(split_level) / 7);
            if (factory->BoxExists(idx)) {
                colliderList.insert(idx);
                // 存在していた空間レベルを保存
                exists_max_split_level = std::max(exists_max_split_level, split_level);
            }
        }

        // 探索すべき空間レベルの決定（+1する）
        exists_max_split_level = std::min(exists_max_split_level + 1, factory->GetSplitLevel());
        // 探索空間レベル基準の座標系で次のグリッド座標を決定する 
        auto next_grid = gridForward + factory->CalculateGridCoordinate(pos, exists_max_split_level);
        // 探索空間における分割サイズを計算する
        auto size = rootAABB.size() / static_cast<float>(1 << exists_max_split_level);
        // 次のグリッドサイズから、探索空間における座標を算出する
        next_pos = Vector3D(next_grid.x * size.w, next_grid.y * size.h, next_grid.z * size.h) + rootAABB.bpos;

        // 次のグリッドの座標がシーンから出ていたら終了
        if (!rootAABB.Contains(next_pos)) {
            break;
        }

        float ax = ray.dir.x != 0.0f ? std::abs((next_pos.x - pos.x) / ray.dir.x) : FLT_MAX;
        float ay = ray.dir.y != 0.0f ? std::abs((next_pos.y - pos.y) / ray.dir.y) : FLT_MAX;
        float az = ray.dir.z != 0.0f ? std::abs((next_pos.z - pos.z) / ray.dir.z) : FLT_MAX;

        if (ax < ay && ax < az) {
            pos += ray.dir * ax;
            grid.x = next_grid.x;
        }
        else if (ay < ax && ay < az) {
            pos += ray.dir * ay;
            grid.y = next_grid.y;
        }
        else if (az < ax && az < ay) {
            pos += ray.dir * az;
            grid.z = next_grid.z;
        }
        else {
            pos += Vector3D(ray.dir.x * ax, ray.dir.y * ay, ray.dir.z * az);
            grid = next_grid;
        }
    }

    return colliderList;
}

　これは小さなロジックの変更だが、非常に疎な空間（特に、空に向けて放たれるようなレイ）を進むレイのステップを著しく減らすことができる。これにより、私の開発中のレンダリングエンジンのベースのシーンであるコーネルボックスにおける、フォトン散布処理は10倍程度高速化された。

　しかし依然として現在はフォトンが1万個、バウンス制限3回という制約下でも、フォトン散布処理に10sec、ライトマップベイク処理に40secかかっている。特に、シーンの複雑度が増す場合、ライトマップのベイク処理はより増加することが想定される。高い品質のフォトンマッピングのために、まだまだチューニングの余地はあるだろう。（事前計算処理とはいえ、毎度数分待たされるとげんなりする）

　フォトンマッピングの実装は大体終わっているが、ノイズというかにじみがとれない。フォトンの量を十分に増やせばいずれなくなるのか、それとも何かミスがあるのか、判断がつかないでいる。全然わからない。私は雰囲気でGI実装をしている。