BFS / Frustum Bounds による交差判定高速化の紹介
要約
- BVH のような AABB を用いた Polygon 管理用のデータ構造と、複数の Ray をまとめた Ray Packet の交差判定を高速化する方法を調査した
- [Overbeck 08] Ray Packet 内の Ray が 1つでも AABB と交差するか否かを基準に、 幅優先探索(BFS) により BVH のノードを辿る
- [Reshetov 07] 葉ノード内の Polygon と Ray Packetの交差判定の前に、葉ノードの AABB 内の Ray のBounding Box相当のものを角錐台(Frustum)を使って表現し、枝刈りする
はじめに
この記事は、レイトレ合宿6 アドベントカレンダー 7週目の記事です。
この記事のきっかけは「パストレってどうやって高速化してるんだろう」と気になって色々調べたことでした。調査にあたっていくつかの論文の Previous works を見ると「幅優先探索(BFS)のように交差判定していく」という記載をいくつか見かけることに気づきました。BFS は馴染み深いのに概要が全く想像できず、しかも多くの手法のベースになっているように見えたので、調べてみました。
今回は、僕が理解できた中で一番基礎的に見えた方法を紹介します。ベースは[Overbeck 08]ですが、[Overbeck 08] における Frustum Bounds による枝刈りの詳細がわからなかったので、別途 [Reshetov 07]を参考にしました。各文献に記載された内容を網羅出来ているかと言われると怪しいですが、根本的なアイディアは伝えられていると信じます。
なお、レイトレ等は仕事や研究等で触れたことはなく、間違っている可能性は十分あるので、
間違い等含まれていましたらご指摘等頂けると幸いです。
BFS による BVH の探索
Ray と Polygon の交差判定は、 BVH(Bounding Volume Hierarchy) 等の Polygon 管理用のデータ構造を用いてもとても重いため、高速化したいです。ここでは、AABB(Axis Aligned Bounding Box) を用いた BVH に対して、複数(256個位)のRayが入ったRay Packet と交差判定する設定を考えます。
基本的な考えは BVH のノードを BFS の順で辿るのですが、Ray Packet 内の Ray が1つでも子ノードのAABBに交わる場合、子ノードはBFSの探索対象になる というのが重要です。
具体例(2次元)
数字の 0〜6 は AABB につけたindexで、0〜2が内部ノード、3〜6が葉ノードです。葉ノードには Polygon が複数個含まれています。
Ray Packet 内の Ray は簡単のため a,b,c の3本とします。
BVH における AABB/Polygon の管理構造をグラフにすると下記のようになり、このグラフに対して Root のノード 0 から BFS を行います。
BFS を行うための Queue を用意すると、Queue の初期状態はノード 0 が入っています。
次に、 ノード0 の子であるノード1, 2と Packet 内の Ray a,b,c の交差判定を行います。
この時、ノード1, 2 は a,b,c の少なくとも1つと交差するので、 BFS で1, 2 の更に子ノードを辿る対象になります。
また、各ノードが a, b, c の順に交差判定する場合、aと交差した後で b, c とは交差判定を行いません。
同様に、ノード1, 2の子であるノード3〜6と交差判定を行います。
この時、ノード6は Ray Packet 内のどの Ray とも交わらないため、Queueには入れられません。また、ノード3は Ray aの1つしか交わりませんが、Ray が1本しか交差しなくとも Queue に入れられます
これで、交差判定させる葉ノードが列挙できたので、後は葉ノード内の Polygon と愚直に交差判定していきます。
交差する Ray の絞り込み
先の例を少し考えると、Ray のどれか1つでも AABB と交差すると探索を全く省略できないので、結局あんまり効率化出来てない気がします。なので、もう少し効率化を考えたいです。
このアルゴリズムでは 1個でも AABB と交わる Ray が見つかれば後続の Ray の交差判定が打ち切れるため、 AABB と交差しないことがわかっている Ray を除外できれば、AABB と Ray Packet の交差判定を効率化できます。また、BVH のノード毎に交差しうる Ray を覚えておけば、葉ノード到達時に Polygon と交差判定する Ray の本数も減らせます。
重要なのは、親の AABB で交わらないRayは、子の AABB でも交わらないという点です(BVH なので当然ではありますが)。Rayを除外する方針はいくつか考えられます。
交差判定対象の Ray Packet の区間を管理する
Ray Packet 内の Ray に番号をつけ、あるノードに当たった最も小さい/大きい番号を区間としてノード毎に覚えておきます。
子ノードを探索する時、この区間の外を探索する必要がないので、少し速くなります。
交差しなかったRayを除外する
今度は、Ray番号のうち終端の番号だけ覚えておきます。また、先頭のRayから順次交差判定を行っていき、交わらなかった場合はその番号と終端の番号をswapし、交わらなかったRayの分だけ終端の番号を減らします。
この場合、探索するにつれて番号が入れ替わるため、まだ交わる可能性のある Ray 番号を全て子ノードに伝える必要があります。
更なる効率化へ向けて
この手法がうまく行く前提として、Ray Packet 内の Ray が同じような方向である等、Ray のコヒーレンシが十分に抽出出来ていることがあると思います。 [Overbeck 08]では、コヒーレンシの抽出まで焦点が当てられていないようだったので、コヒーレンシをどうやって抽出するかについては別の文献を参照する必要があるようです。
Frustum Bounds による枝刈り
Ray Packet がある葉ノード(直接ポリゴンを管理しているノードで、前の例だとノード3~6) に到達した場合を考えます。素直な実装であれば、Ray Packet 内の交差しうる Ray と、BFS で辿れた葉ノード内の Polygon の全パターンを総当たりで交差判定する必要があり、これを高速化したいです。
高速化の方法として、葉ノードの AABB 内部において Rayが通りうる領域を算出し、それとPolygon を事前に交差判定させる ことを考えます。Ray の Bounding Boxに相当するものになります。通常であれば 1つのPolygon に対して Ray Packet 内の Ray と1つずつ交差判定の必要がある所を、Frustum と交差判定するだけで全て省略できる可能性があります。
[Overbeck 08] でも紹介されているのですが、より詳しい説明が[Reshetov 07]に記載されていたので、こちらに書かれた内容を調べてみました。
ざっくりした処理の流れは、
- Frustum の構築
- Frustum の側面によって、Polygon が分断され得るかを判定
- Polygon が存在する平面上で Frustum と Polygon が交差しうるかを判定
- 枝刈り出来なければ、愚直に各 Ray と Polygon のペアを交差判定
という流れになります。
Frustum の構築
Frustum 構築の際、xyz 軸のうち 全ての Ray が同じ方向を向いている、最も Bounding Box の辺が長い軸 を選択します。以下では x 軸が選択されたとします。 Bounding Box の yz 平面に平行な 2 枚の長方形に対して各 Ray との交点を求め、交点を囲う xyz の軸と平行な2枚の長方形で Frustum を構築します。やや見にくいですが、図に表すと下記のようになります。
全ての Ray が同じ方向を向く制約が必要な理由は、より強力な枝刈りが行えるためだと思っています。全ての Ray から見て、ある Polygon の後ろに隠れている Polygon の検出を行い、余分な交差判定を省ける場合があるようです(詳細は省略)。
Frustum の側面によって、Polygon が分断され得るかを判定
Frustum には x 軸方向に伸びる平面が4枚あり、ここでは側面と呼びます。Polygon の各頂点に対して、Frustum の側面から見て内外どちらにあるか判定します。 Polygon の全頂点が外側にあると判定される側面が1枚でもある場合は、Ray Packet と Polygon が交差することはないです。 側面の外側にあるか否かの判定は、①側面の法線、②側面上の点から Polygon の頂点までの方向ベクトル、の内積の正負を計算すれば良いです。
Polygon が存在する平面上で Frustum と Polygon が交差しうるかを判定
2.のチェックにより Frustum と Polygon が交差する可能性がある場合、更なる交差可能性の判定を行います。ある Polygon が存在する平面上で、Frustum における 4 本の Corner Ray がどこに交わるかを計算します。ここで、Corner Ray とは、下図に示すような Frustum の側面にある辺のことです。
Corner Ray と Polygon 平面との交点で作られる四角形の一部が Polygon の内部に含まれる場合、 Ray Packet 内の Ray が Polygon と接触する可能性があります。Polygon の外側にある場合は、Ray Packet が Polygon と交差することがないため、交差判定を省くことができます。
おわりに
Frustum Bounds による枝刈りはそれ自体で完結している雰囲気がありそう(適当)ですが、[Overbeck 08] の方がどう発展していくかは個人的に結構気になるところです。これを基に色んな文献の高速化手法がわかるようになるといいなーと期待していますが、流石に厳しそうな気もしています。
また、本当は性能評価などもやってみたかったのですが、時間が思うように取れずにアドベントカレンダーでは断念してしまいました。
せっかく調べたので実装やりたいですが、当日までに組み込む余裕があるだろうか…
より詳細や正確な情報を知りたいと思った方は、直接各文献に当たることをおすすめします。ベンチマーク結果とかも載せようと思いましたが、一言で言うのが難しかったので諦めました。
肝心の合宿当日の話ですが、当日に向けてレンダリングしたいシーンも最近やっと決めたので、何とか面白い絵になることを祈ります。残り日数も少ないですが、間に合うといいなー
参考文献
[Overbeck 08] R. Overbeck, R. Ramamoorthi, W. R. Mark, "Large Ray Packets for Real-time Whitted Ray Tracing", 2008
[Reshetov 07] A. Reshetov, "Faster Ray Packets-Triangle Intersection through Vertex Culling", 2007
Weathernews Programming Contest 参加記
だいぶ遅れてしまいましたが、自分の記録用に書きました。
先日、 Weathernews Programming Contest (お天気コン)に参加しました。
結果は4位でした。全体的にかなり強実装なのが辛かったですが、動くenc/decを初めて作ったし超楽しかった!
やったこと
大まかに、
- 画像を10x10のブロックに分割
- ブロック内の画素値を予測
- 予測誤差の相関除去
- 2.3.の情報+3で残った誤差を算術符号で圧縮
という流れでした。
1.画像を10x10のブロックに分割
予測誤差を抑えるゲームと思っていたので、付加情報を送ってでも誤差を下げようと考えてました。 予測方法や相関除去の付加情報が増えすぎないように、ブロック単位で取れる選択肢を制限しています。 大きさを2べきでなく10x10にした理由は、適当な大きさで画像をピッタリ埋め尽くせるためです。
付加情報を周囲の状況から推測する方法もあります(machyさんはそうしていたよう)が、 最初は素直なものを作ろうと思ってそうしました。 付加情報を送れば予測誤差はまぁ下がるはずなので、どっちがいいかはともかく、トレードオフにはなってるはずです。 きっと推測すべきなのだとは思います。
ブロックによる分割は、画像にとって自然な分割にならないから辛いことも多いのですが、目視で時間変化を見て大丈夫そうな雰囲気(適当)を感じたのと、自分にとって慣れてるということでブロック単位の管理を選択しました。
2.ブロック内の画素値を予測
ブロック単位で、以下の3つの選択肢から選びます。 どれを選択するかは、後述する相関除去後の差分が小さいやつが選択されます。
- 同じ画像内の1つ上の画素と1つ左の画素の平均値をピクセル単位で計算して予測値とする
- 同じ波長画像で1つ前の時間の画像でブロックマッチングしてコピー
- 同じ時刻の1つ短い波長の画像で、同じ位置のブロックを線形変換してコピー
- 何も考えずにコピーできるように、圧縮が終わった画像は参照される画像にスケールを合わせて保存。
- 1つ短い波長ブロック → 圧縮対象ブロックの線形変換をブロックごとにやった。必要な係数は周辺画素から動的に求めてた。
- 同時刻に撮った画像同士のマッチングなので、動き量は0に限定してもよいだろうと判断した。動き量を送らなくて良くなるのもポイント高い。
ブロック単位で2bitの情報 + 2.での動き量を伝えます。 動き量は似てることが多いため、動き量どうしの差分を伝えます。
最終的に、最初の時刻に撮られた画像の圧縮率が他と比べて圧倒的に悪かった(2. の予測方法を使えないため)です。 なので、3. の部分をどうにか出来ないかなぁというので結構苦労してました。
3.予測誤差の相関除去
予測誤差を眺めてみると、誤差の乗り方も似ていることがよくあるので、予測誤差どうしもブロック内で差分をとります。 差分のとり方は
- 横方向に差分を取る
- 縦方向に差分を取る
- 斜め方向に差分を取る
- そもそも差分を取らない
を全通り試して、ブロック内で差分とった結果の絶対値和が一番小さいやつを選んでいます。 1〜3は、ジグザグする感じの訪問順序になります(伝われ)。 ブロック単位で2bit + 最終的な差分の情報を伝えます。
4. 2.3.の情報+3で残った誤差を算術符号で圧縮
まず前提として、ビット列は01の出現頻度の偏りが大きいほど圧縮しやすいです。直感的には、1000桁のビット列を覚える時、①全部0のデータと②01がランダムに混じるデータを比べると、前者が覚えやすい、位で十分と思います。
今回は、 付加情報として送付したビット情報を全部別々のビット列とみなして算術符号化にかけました(実際にビット列を分割したわけではなく、01の発生頻度の確率だけ分けた)。ビット列毎に01の出現頻度は、直前に受け取ったビット列の出現頻度を元に、動的に計算しています(適応的算術符号とかでググると色々出てくると思います)。直感的には、圧縮時と復元時に同じパラメータを使って同じ手順で処理を再現できるなら、特に正確な固定値の出現頻度値を使わなくてもOKという話です。
ビット列が生成された由来や状況に応じてビット列を別物として扱う(=コンテキスト分類する)ことで、01の偏りを強めることができる場合があります。が、今回は特に強いコンテキスト分類まで考えられていませんでした。
コンテスト中に考えていたけどやり残したこと
- ある波長画像を別の波長画像にマッピングする処理
- 予測方法1. で、平均を取らずに誤差最小の係数を線形回帰して計算
- 予測方法2. 3. の方が圧倒的に効果があると思っていたので、結局手付かずになってしまった
- 一旦最適な予測方法が決まった後で、予測方法1. を使うブロックに限定して改めて線形回帰用の係数を計算し直す
- 予測方法2. 3. のブロック内画素位置を回帰対象に含める必要がない
- と同じく、画像間の重ね合わせ部分が最重要だと思ってたので、手付かずになってしまった。
- 誤差最小ではなく、符号量最小の選択肢を選ぶ
- 厳密な測定をしようと思うと、様々な状態に依存するので正直無理ですが、近似的には出来ます。最大サイズの画像は時間的にかなり厳しいですが、小さい画像の最適化ではやるチャンスがあったかもしれないです。
- 画像の拡大にBilinearを使っていたけど、もうちょいましな拡大方法を使っても良かったかも
余談
- 圧縮率と関係無いですが、エンコーダとデコーダって同じ順番で付加情報なり誤差を読まなきゃいけないんだけど、どうにかして共通化する方法ないのだろうか(これのせいで、デバッグ時間的に微ロスを重ねていたので)
感想
- 前職でチョットやってたので、いつか作りたいと思ってた。実際動くものが作れてからはテンションが5割増くらいだった。
- いろんな処理を無限にバグらせたし、思ったことが一部できずに終わってしまったので、無限に実装力が欲しくなった。精進しなきゃ。
- データの圧縮って個人的にすごい好きなんですが、「圧縮っていいよね!!」て言ってる人たちが少なからずいたようで、何か嬉しかった。
- 今年からマラソンの参加記を書こうという気になったので、これからは終わってすぐ公開できるように準備しておきたいと思います。よろしくお願いします。