レイトレ合宿7 徒競走でやったこと

はじめに

今回のレイトレ合宿では、イベントの一環として BVH の構築/交差判定を高速化する徒競走というイベントが行われました。今いる会社的なあれで高速化に詳しくなりたいというのもありましたが、やってみるとこれがかなり楽しくて、 本題の CG 作品の提出を忘れて 取り組んでいました。そこそこ勉強込で時間をかけたので、何とか忘れないようにやったことを記したいと思います。

実際のソースコードこちら になります。
ソースコードの規約が統一されていなかったり、部分的におかしなコードがあったりしますが、昔書いたコードを必要最低限の改変で持ってきたりしてたので、許して下さい。

ルールは、BVH の構築 + 交差判定の時間が短い人の勝ちです。

とりあえず、効果が 0 でなさそうな、やったことの概要を記すと以下のようになります。
抜け漏れはあるかもしれません。あと、どうやらどこかがバグっているらしいのですが、特定できていません…

やったこと

3次元ベクトルを Expression Template 化

ソースコード

最初、SIMD 化もされていないコードだったので、BVH 構築/交差判定両方で一時オブジェクト作成を削減するためにやりました。最終的に、交差判定ではどのみち一時オブジェクトは作成されないので、最終的な貢献度は薄かったかもしれません。今見ると、きちんと作ろうとして途中で投げ出した後があって面白い。

ポリゴンの当たり判定を正規化した空間でやる(名前がわからない)

ソースコード

Polygon の当たり判定も、もっと速い方法絶対あるだろーと思って色々ググった結果、この記事 に行き当たりました。仕組みは割と簡単で、三角形を (0, 1, 0), (0, 0, 0), (1, 0, 0) の空間に写してから当たり判定するというものです。手元では、そこそこ効果があったように見えたので、採用しました。

三角形をなさない polygon の除去

ソースコード

本来は高速化の意図ではありませんでしたが、上記のポリゴン当たり判定方法を採用すると、obj の一部で nan を吐きまくったため、抑制のために入れました。誤差とはいえ、多少は影響あったかもしれません。

Top Down な SAH ベースの BVH の並列構築

ソースコード

この記事 を色々眺めた結果、とりあえず知ってるやつをやろうということでこうなりました。構築は Task Queue を使って並列化しました。詳細はアドベントカレンダーの記事 に書いたので、興味があれば見てもらえると嬉しいです。

std::stable_partition を使った O( n log n ) 構築(nはポリゴン数)

ソースコード

この内容は、 Bonsai の論文 で知りました。

TopDown SAH ベースの構築をする時、(僕の理解が間違っていなければ) polygon id を管理する配列を x / y / z 軸の3つ分用意し、 x / y / z 軸のそれぞれで sort した状態を保持します。 node を left / right で分けた時、この polygon id の配列を left / right の子ノードの分に分離しつつ、それぞれの子ノードの中では x / y / z 軸で sort されている状態を保ちたいです。
例えば x 軸で bounding box を 2分割したとすると、x 軸分の polygon id の分割は容易です(あるしきい値より左が left, 右が right なので)が、y / z 軸については left / right にいるべき polygon id がごちゃごちゃしているので、これを std::stable_partition を使って分離します。x 軸の polygon id list を見れば、「特定の polygon id が left / right にいるべき」というのは容易にわかるので、各 polygon が left / right のどちらにいるべきかフラグを持たせておき、 y / z 軸の polygon id を std::stable_partition で振り分けます。この時、stable 付きのものを呼ぶと、「 left どうし、right どうしの polygon ID の順番は元の配列のソート順を保つ」という性質があるので、y / z 軸で sort された状態も保つことが出来ます。

8分木の構築 + Node の当たり判定 8並列

ソースコード(2分木の展開)
ソースコード(8 並列 intersect)

定石がよくわからなかったので、2分木を作ってから展開して8分木にしました。
node の当たり判定については、reference 実装のものをそのまま SIMD に描き起こした感じなので、もうちょっと工夫が出来たかもしれません。

8 並列当たり判定で交差した node を効率的に抜き出す

ソースコード

この技術的課題は、この論文で知りました。

8 並列で node の当たり判定をした後、距離順に sort して 子ノードの traverse 順序を決めたいのですが、8要素のソートは重いので、交差した個数に応じて sort を分けるのがいいと記載がありました。これを実現するためには、交差した部分の index はどれかを高速に特定する必要があります。
ここで、「子ノードの bounding box との距離が INF でない bitmask を作成して、bitmask の位置にある index を 右詰めして抜き出す」という処理だと思うと、 pdep / pext 命令が使えることに気づきます。ざっくり書くと、

  1. 入力:1bit x 8 の bitmask当たり判定結果
  2. 8bit x 8 の bitmask を作成する
    • pdep で 8bit おきに bit を配置するように設定し、0xFF をかけて作る
    • avx512 でも大体似たような感じでいける
  3. 0x0706050403020100 は、64bit の数字に index を埋め込んだもので、ここから ↑の bitmask に被るところを pext で右詰めして保存
  4. 3.で作られた整数は実質 8bit integer なので、32bit interger に拡張して store
  5. 何個当たり判定したかは、入力の popcount を数える

みたいな感じになっています。この処理も参考にしたブログがあったのですが、メモに残っていませんでした。
なお、肝心のソート用関数を個数に応じて使い分けるのは、あんまり速くならなかったので、詰めた後は std::sort してます。

Polygon の当たり判定 8並列

ソースコード

これは本当に 1行ずつ SIMD にするだけ。

Polygon の(メモリ配置的な)並び替え

ソースコード

index の間接参照みたいなのを減らすためにやりました。バグらせるのが怖かったので愚直に配列作る + copy をしました。
こういうのが他にも複数あるので、多分構築は遅かったと思います。

(補足)やったけど、効果がなかったやつ(上2つはバグってた可能性大)

Bonsai

構築が明らかに遅かった+改善の余地があったので、とりあえず Bonsai しておけばいいだろうと思って実装したのですが、部分的に手抜きして作ったせいかめちゃくちゃ交差判定が遅くなり、時間切れで諦めました。今回唯一悔やまれる点があるとすると、これのきちんとした実装が出来なかったことだと思います。

insertion based optimization

論文

見た瞬間、これはやりてぇと思って実装しました。が、何か思うように速くならず…かつ、論文とか見ると結構構築が重そうで、終了判定タイミングも結構難しいなと思ったので、優先度をそこそこ落としてしまいました。

AVX512 化(16分木 + ポリゴン 16並列当たり判定)

ほとんど AVX2 の命令を移植するだけでしたが、AVX512 では bitmask が一部 __mmask16 のように圧縮されている形式だったので、その扱いを少し変える必要がありました。手元では全然速くならなかったので、結局諦めました。

SIMD の bitonic sort

ソートが遅いということで、子ノードの交差個数が多かったらこっちの方が速いかもなーと思って実装はしてみました。やっぱり交差するノード数の期待値の問題なのか、多少調整しても速くなったか微妙だったので、やめてしまいました。

intersectAny で leaf node を優先して当たり判定

シャドウレイを飛ばす場合などでは、 polygon に 1つでも当たったら終わりなので、距離よりも当たりやすさを優先すべきだよなぁと思っていました。というわけで、Leaf Node に当たったらそっちを優先的に調べることにしていたんですが、比較条件が複雑な割に改善ケースが少ないせいか、遅くなってしまいました。

開発過程でやったこと

visual studio のプロファイラを眺める

大体ボトルネックは全てこれで把握していました。途中 inline 展開された部分の把握が不可能だったので、もうちょっと何とかならなかったものか。

途中まで SIMD 化して残りを愚直に実装する

SIMD の処理を最初から最後まで書くと大抵 100% バグっていたので、「途中まで SIMD で書いて残りは愚直に書く」をして、愚直をどんどん減らすような方針で書いていきました。

intel intrinsics guide をざっくり眺める

https://software.intel.com/sites/landingpage/IntrinsicsGuide/

そもそも自明じゃないテーマで本格的に SIMD 高速化をするのが初めてだったので、とりあえず何ができるのかをざっくり把握することをやりました。1個ずつというのは到底出来ないので、カテゴリごとに似たような命令はすっ飛ばして、雰囲気を感じていました。

BVH の validator をたくさん書く

ソースコード

BVH の構築後に色々データをいじりまくっていたので、一度バグりだすと何が原因なのか特定するのは難しいです。この時、処理における事前条件/事後条件を書き出して、処理する度にチェックするコード debug build で入れてました。かなり役立ち度は高かったように思います。

開発過程でやらなかったこと

cache miss のプロファイルをとる

これは本当はやりたかったのですが、 VS 2019 でやる方法がわからず。わかる人がいたらぜひ教えて欲しいです。

BVH の visualizer とか評価尺度を見て考察する

割と定数倍高速化をしてたら終わってしまったので、BVH の構築は top-down SAH のものからほとんど進歩していませんでした。もうちょい構築の工夫を入れてたら、 SAH とか体積とか当たり判定の分布とかを真面目にとっても良かった気がします。

UT を書く

reference 実装がきっちり与えられていたので、あまり必要ではありませんでした。今回は validator があれば大体問題なかったような気がします。

最終的な速度

正確に計測するのは大変なので諦めましたが、本番環境(CUI 版) + 本番データセットの hairball で 3回計測したときの中央値を取ると、最初と最後でこれ位の差がありました。あと、embree みたいな現状最強格のやつとも比較したかったのですが、使い方がよくわからなかったので諦めました。まぁこれだけ速度差があるといっても、バグを修正しないと何とも言えないよなぁという気持ちはあります。

build(ms) intersect(ms)
reference 9,491 202,253
最終版 7,660 9,500
速度比(ref / 最終版) 1.24 21.29

終わりに

振り返ってみると、これは本当にレイトレ合宿の参加記なのかという気持ちになってきました。とはいえ、こうしてまとめてみると、意外と頑張ったなという気持ちになってきたので、とりあえずは良かったかなと思います。初めてガチで最適化を頑張ってみましたが、同じ言語でも、これだけ違うんだなぁみたいな感慨深さはあります。

来年こそは、綺麗な絵+CUDA を頑張りたい気持ちです。

役に立つかもしれない、Task Queueの作り方

はじめに

この記事は、レイトレ合宿7 アドベントカレンダーの記事です。

様々なタスクをマルチスレッドで処理するために Task Queue を使う場合があります。 OpenMP 等のように気軽に並列化できる処理もたくさんありますが、例えば BVH を構築する場合、 この記事 内に紹介されるような Task Queue が使われることも多いかなと思います。 BVH を何種類も作ったり、他の処理でも Task Queue は使いどころがある気がします(具体的にたくさん挙げられるとは言っていない)し、可読性向上の観点からも、モジュールとしては小さいですが、独立したモジュールを作りたい気持ちになります。

Task Queue を独立したモジュールとして作ろうと思って色々考えてみると、終了条件をどう表現するかが比較的難しいポイントと感じました。 例えば、今は Queue が空でも、処理中のスレッドが新しくタスクが放り込まれるかもしれないため、Queue が空だという条件は使えません。 また、処理が終了したら迅速にスレッドを終了する必要もあるため無限ループするわけにもいかず、終了条件を記述した関数を渡すのもあまりしっくり来ませんでした。

このような時、数値で表現される進捗を queue に登録する方法だと、比較的上手く表現できそうと感じたので、紹介してみることにします。 例えば、BVH の構築であれば、Leaf Node を作成したタイミングでその Leaf Node に含まれる primitive (Polygon とか) な object の個数を進捗として報告し、最終的に進捗が BVH 構築時に渡される object の個数に到達すれば終了して良い、といった感じになります。

実装

c++ での実装例はこちらになります。TaskQueue の実装と、BVH の構築とよく似た処理を行うサンプルコードを一緒に置いておきました。 この記事を書いた環境の都合上、std::optional の代わりに boost を使っています。

github.com

例えば std::queue との差分としては、

  • empty() や size() 相当の処理等はなく、進捗登録とタスク登録/取得しかできない
  • front() と pop() は同時にやる
  • 要素が取れる保証がないので、要素取得は optional で包む

といった感じになります。

おわりに

今回はこれ以上深追いしませんでしたが、もっといい設計はできるかもしれないので、いい方法等ありましたら教えて頂けると幸いです。また、この Queue に名前がついていたり、こんな使い方が出来そうみたいな例などあったりしたら、教えてもらえると嬉しいです。

レイトレ合宿6 参加記

要約

  • レイトレ合宿6 に参加した
  • めちゃくちゃ面白かった
  • アセットは適当に作っても美しくならない

はじめに

9/1(土)〜2(日)にかけて、レイトレ合宿6 に参加してきました。最終的に、僕はこんな絵を書いて、 16 / 19 位でした。

f:id:xyz600:20180903221808p:plain

スライド

speakerdeck.com

コード:https://bitbucket.org/xyz600600/xyz_renderer/src/v0.3/

コードは Rust で書かれていますが、初心者なのであんまり Rust っぽいコードとして参考になるかはわかりません… (一部 unsafe に const ptr の dereference とかやっている)(解決方法はわからず)

(当然)上位の方々に比べるとイケてる絵は出せなかったのですが、来年やりたいことを忘れないように、ブログにまとめることにしました。来年は僕もイケてる絵を出力してみせるぞ!!

本番前と本番とで、やったこと思ったことを中心にざっくり書いていこうと思います。
全部書く時間がなかったので、部分的に結構抜けてますが悪しからず。

本番前

5月

理屈を理解した上でレンダラを書いたことなかったので、今回の合宿の目標をあれこれ考えていた。最終的には

  1. オブジェクトをたくさん配置する
  2. 高速化を頑張る
  3. 理論を理解する
  4. アセットを自分で作ってみる

とかかなーと思っていた。この時は、適当にアセットを作ってもパストレだしきれいな絵が出るだろうくらいに思ってたので、アセットを考えるのは完全に後回しにしていた。そんなこと全くなかったんだけどね。

また、どうせなら新しい言語を習得したいと思って、Rust で書くことにした。そのため、5月は Rust の勉強をしつつ、Ray Tracing in One Weekend を読み進めていた。

6月

仕事が忙しい。レイトレ業はお休み。

7月

あれ、何もやってないけど残り2ヶ月しかないぞ…??
という訳で、焦り始めて割とここから頑張った。

色んな資料を見ながらレンダリング方程式の理解に挫折していたが、

www.slideshare.net

を見て疑問点が氷解し、念願のレンダリング方程式を会得した。@Shocker_0x15 さんの記事は割と最後までずっとお世話になっていた。

設計で悩んでいたこともあって、月末あたりでやっとコーネルボックスが動き始めたが、並列化をしようと思ったらライフタイム周りでたくさんコンパイラに怒られるようになり、厳しい気持ちになった。

あと、すっかり忘れてたけどアドベントカレンダーのネタを考えないと明らかにやばいので、月末あたりで書く内容をやっと決めた。結局、こんな感じのものを書いた。

xyz600.hatenablog.com

特に、Frustum Bounds とかめっちゃカッコイイと思ったので実装したかったけど、時間的に無理だった。残念。とはいえ、当日に数人から「アドベントカレンダー書いた人」という認知をされてよかったなぁと思ったので、次は実際に実装する技術を調べる感じで、積極的に書いていきたい。

途中でICFPC に出場していたり、その反動で3日位体調崩してたりしてた。

8月

BVH の実装とobj / mtl ファイルのローダが8月序盤で書き終わり、そろそろアセットを決めようと悩んでいたところ、交通渋滞のニュースを見て「車たくさん描画すれば、何かかっこよくね??」と思ってこれに決めた。そんなこと(ry

途中アメリカ出張が入り、時差ボケに苦しみながら書いたコードのあらゆる所にバグを埋め込んでいて、厳しい気持ちになった。

並列化の方法が全くわからなかった最中、プログラミング Rust が届いたことで解決した。rayon すごい。 またこの頃から「Rust PathTracing」でぐぐるようになり、Rust で パストレしてる先人の @gam0022さんの存在を知る。なぜもっと前からぐぐらなかったのか。

なんとなく道路っぽい雰囲気の何かを作り、車を並べてレンダリングしてみて、あまりのリアリティのなさに絶望した。思わず2度見して、3回くらい描画し直した。とはいえ、おそらく今年はこれ以上出来ることもないから、ノイズを出来る限り減らす方針で頑張ろうと決意した。

まずは強そうだし双方向パストレでも実装するか??と思って理屈を理解して頑張って実装していたけど、カメラパスで考慮した画素位置と別の画素位置にライトパスが飛び込むケース(伝われ)を完全に見落としていて、並列化の戦略と相まって実装困難になる。という訳で、諦めて合宿が終わってから実装することにした。

とはいえ、単純なパストレではノイズが全く取れず、 Next Event Estimation とか MIS を理解した上で実装したくて色々調べたり考えたりを繰り返していたのだけど、

rayspace.xyz computergraphics.stackexchange.com

あたりを見て割ときちんと理解できた。積分形式のままで考えることで、 MIS も割とちゃんと理解できた、はず…

よく考えたら、リアリティがない理由の1つは車も色も全部同じだからだよなぁと思って、せめて色だけでも変えることにした。適当に車体の色を表現するmtl ファイルの箇所を特定して、車体を6色の中からランダムに割り振ったり、車の位置をランダムに移動させたりしてごまかした。

このままフィニッシュできるかと思いきや、本番環境で動かない旨の通知を終了5時間前に受け取り、めっちゃ焦った。動的ライブラリが見つからないことが原因らしく、ここの一番最後の方法で --target=x86_64-unknown-linux-musl をつけてbuild したらOKだった。実行速度が15% 位落ちたけど、動いて本当に良かった…

こういうところは新しめの言語の方が強そうな気がしたけど、他のネイティブコードを吐く言語はどうなんだろうか。

本番

初日

CEDEC の話で盛り上がると過去の記事に書いてあったので、当日の資料 を見て意識を高める。なるほどわからん(ベイクが何だか知らなかったため)けど発表聞いたら絶対楽しかったんだろうなぁ。

到着するなり温泉に向かい、道中で色々話しながら親交を深めた。普段そんなにCGについて話す訳じゃないので、新鮮な話題ばっかりで楽しかった。途中、@gam0022 さんとお会い出来たので、Rust めっちゃ大変だったと伝えたら分かってもらえて嬉しかった。

特別イベント楽しみにしてたら、いきなり検定が始まるわ、三葉レイの応援メッセージが届くわで、かなり笑った。企画力高いなぁというのもそうだけど、これめっちゃ準備大変では… 運営の方々凄すぎる。 @ZinTwitt さんと @ishiyama_ さんのセミナーも普段聞けない話で面白かったんだけど、発表直前で自分のバイナリが果たして動いているのかが気になり始め、集中できなかった。ぐぬぬ

いよいよ発表会が始まったが、みんなめっちゃきれいな絵を次々に出しててめっちゃ焦った。自分でモデリングしている人もいるし、ただただ凄いなぁという感じ。僕のレンダラも無事動作したとわかってホッとしながら発表した。初めてにしては色々やってて凄いね〜 と@redqueenrenderさんに言ってもらえ、めっちゃ嬉しかった。発表会の後はご飯食べたり花火したりでワイワイしてた。

発表会後、色んな人に聞いてみると、Houdini や blender を使って1〜2日使ってモデル作ったぜって人がそれなりにいて、僕もすごく使ってみたくなった。自作モデルどう考えてもカッコイイんだよなぁ。僕も作れるようになりたくなった。

旅館の部屋で、@q_cinnamonさんにBCDというデノイザを詳しく教えてもらった。そもそもパストレの絵をデノイズする発想がなかったし、理論的にも面白いよと言って下さったので、落ち着いたら論文読むぞ。

絵に得点をつけるの凄い難しいよなぁと思いながら、寝る前までにみんなの絵に得点をつけた。みんなめっちゃきれいな絵なんだけど、どれか1つと言われたら @ushiostarfish さんの布が最高に好きだったなぁ。

2日目

いよいよ結果発表だったけど、流石に順当な順位だった。商品貰えなかったなーと思いきや、プレゼントが1個増えてじゃんけんで勝ったので、PS本体とレイトレ教材をもらった。翌日、きちんと教材をやりきって、無事にレイトレを完全に理解した。hardは無理ゲー。

船の出発まで時間があったので、海で飛び込みして遊んでた。天気が荒れたり晴れたりを繰り返してて、自分の荷物が危うく水没するところだった。次から島に行くときはゴミ袋を持っていくべきなのかも。

謎の力で@shuto_fmsさんと帰りの船の席が隣だったので、カッコイイ絵に関して色々話した。それまでは物理的に(ある程度)正しい絵が簡単に作れるってすげぇって思ってたのだけど、単純にきれいな絵が出せるのは楽しいよなぁという気持ちが強くなった。僕も「今回の作品は〇〇をこだわってみた。めっちゃカッコイイでしょ??」みたいなコメントが言えるようになりたい。

所感

  • 理論の勉強をきちんとしてよかった
    • 一部でいいから詳しい人の言ってることがわかると、めちゃくちゃ面白い
  • 今回の合宿に参加して、CGの絵がめっちゃ好きになった。
  • 来年こそはきれいな絵をレンダリング出来るようになりたい。

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次元)

f:id:xyz600:20180725200554j:plain

数字の 0〜6 は AABB につけたindexで、0〜2が内部ノード、3〜6が葉ノードです。葉ノードには Polygon が複数個含まれています。
Ray Packet 内の Ray は簡単のため a,b,c の3本とします。
BVH における AABB/Polygon の管理構造をグラフにすると下記のようになり、このグラフに対して Root のノード 0 から BFS を行います。

f:id:xyz600:20180725200236j:plain

BFS を行うための Queue を用意すると、Queue の初期状態はノード 0 が入っています。

f:id:xyz600:20180725201245j:plain

次に、 ノード0 の子であるノード1, 2と Packet 内の Ray a,b,c の交差判定を行います。
この時、ノード1, 2 は a,b,c の少なくとも1つと交差するので、 BFS で1, 2 の更に子ノードを辿る対象になります。 また、各ノードが a, b, c の順に交差判定する場合、aと交差した後で b, c とは交差判定を行いません。

f:id:xyz600:20180725201630j:plain

同様に、ノード1, 2の子であるノード3〜6と交差判定を行います。
この時、ノード6は Ray Packet 内のどの Ray とも交わらないため、Queueには入れられません。また、ノード3は Ray aの1つしか交わりませんが、Ray が1本しか交差しなくとも Queue に入れられます

f:id:xyz600:20180725202208j:plain

これで、交差判定させる葉ノードが列挙できたので、後は葉ノード内の 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 に番号をつけ、あるノードに当たった最も小さい/大きい番号を区間としてノード毎に覚えておきます。
子ノードを探索する時、この区間の外を探索する必要がないので、少し速くなります。

f:id:xyz600:20180725210628j:plain

交差しなかったRayを除外する

今度は、Ray番号のうち終端の番号だけ覚えておきます。また、先頭のRayから順次交差判定を行っていき、交わらなかった場合はその番号と終端の番号をswapし、交わらなかったRayの分だけ終端の番号を減らします。

この場合、探索するにつれて番号が入れ替わるため、まだ交わる可能性のある Ray 番号を全て子ノードに伝える必要があります。

f:id:xyz600:20180725205902j:plain

更なる効率化へ向けて

この手法がうまく行く前提として、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]に記載されていたので、こちらに書かれた内容を調べてみました。

ざっくりした処理の流れは、

  1. Frustum の構築
  2. Frustum の側面によって、Polygon が分断され得るかを判定
  3. Polygon が存在する平面上で Frustum と Polygon が交差しうるかを判定
  4. 枝刈り出来なければ、愚直に各 Ray と Polygon のペアを交差判定

という流れになります。

Frustum の構築

Frustum 構築の際、xyz 軸のうち 全ての Ray が同じ方向を向いている、最も Bounding Box の辺が長い軸 を選択します。以下では x 軸が選択されたとします。 Bounding Box の yz 平面に平行な 2 枚の長方形に対して各 Ray との交点を求め、交点を囲う xyz の軸と平行な2枚の長方形で Frustum を構築します。やや見にくいですが、図に表すと下記のようになります。

f:id:xyz600:20180812120411j:plain

全ての Ray が同じ方向を向く制約が必要な理由は、より強力な枝刈りが行えるためだと思っています。全ての Ray から見て、ある Polygon の後ろに隠れている Polygon の検出を行い、余分な交差判定を省ける場合があるようです(詳細は省略)。

Frustum の側面によって、Polygon が分断され得るかを判定

Frustum には x 軸方向に伸びる平面が4枚あり、ここでは側面と呼びます。Polygon の各頂点に対して、Frustum の側面から見て内外どちらにあるか判定します。 Polygon の全頂点が外側にあると判定される側面が1枚でもある場合は、Ray Packet と Polygon が交差することはないです。 側面の外側にあるか否かの判定は、①側面の法線、②側面上の点から Polygon の頂点までの方向ベクトル、の内積の正負を計算すれば良いです。

f:id:xyz600:20180812120536j:plain

Polygon が存在する平面上で Frustum と Polygon が交差しうるかを判定

2.のチェックにより Frustum と Polygon が交差する可能性がある場合、更なる交差可能性の判定を行います。ある Polygon が存在する平面上で、Frustum における 4 本の Corner Ray がどこに交わるかを計算します。ここで、Corner Ray とは、下図に示すような Frustum の側面にある辺のことです。

f:id:xyz600:20180812120604j:plain

Corner Ray と Polygon 平面との交点で作られる四角形の一部が Polygon の内部に含まれる場合、 Ray Packet 内の Ray が Polygon と接触する可能性があります。Polygon の外側にある場合は、Ray Packet が Polygon と交差することがないため、交差判定を省くことができます。

f:id:xyz600:20180812120628j:plain

おわりに

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 (お天気コン)に参加しました。

wn2017_1.contest.atcoder.jp

結果は4位でした。全体的にかなり強実装なのが辛かったですが、動くenc/decを初めて作ったし超楽しかった!

やったこと

大まかに、

  1. 画像を10x10のブロックに分割
  2. ブロック内の画素値を予測
  3. 予測誤差の相関除去
  4. 2.3.の情報+3で残った誤差を算術符号で圧縮

という流れでした。

1.画像を10x10のブロックに分割

予測誤差を抑えるゲームと思っていたので、付加情報を送ってでも誤差を下げようと考えてました。 予測方法や相関除去の付加情報が増えすぎないように、ブロック単位で取れる選択肢を制限しています。 大きさを2べきでなく10x10にした理由は、適当な大きさで画像をピッタリ埋め尽くせるためです。

付加情報を周囲の状況から推測する方法もあります(machyさんはそうしていたよう)が、 最初は素直なものを作ろうと思ってそうしました。 付加情報を送れば予測誤差はまぁ下がるはずなので、どっちがいいかはともかく、トレードオフにはなってるはずです。 きっと推測すべきなのだとは思います。

ブロックによる分割は、画像にとって自然な分割にならないから辛いことも多いのですが、目視で時間変化を見て大丈夫そうな雰囲気(適当)を感じたのと、自分にとって慣れてるということでブロック単位の管理を選択しました。

2.ブロック内の画素値を予測

ブロック単位で、以下の3つの選択肢から選びます。 どれを選択するかは、後述する相関除去後の差分が小さいやつが選択されます。

  1. 同じ画像内の1つ上の画素と1つ左の画素の平均値をピクセル単位で計算して予測値とする
  2. 同じ波長画像で1つ前の時間の画像でブロックマッチングしてコピー
  3. 同じ時刻の1つ短い波長の画像で、同じ位置のブロックを線形変換してコピー
    • 何も考えずにコピーできるように、圧縮が終わった画像は参照される画像にスケールを合わせて保存。
    • 1つ短い波長ブロック → 圧縮対象ブロックの線形変換をブロックごとにやった。必要な係数は周辺画素から動的に求めてた。
    • 同時刻に撮った画像同士のマッチングなので、動き量は0に限定してもよいだろうと判断した。動き量を送らなくて良くなるのもポイント高い。

ブロック単位で2bitの情報 + 2.での動き量を伝えます。 動き量は似てることが多いため、動き量どうしの差分を伝えます。

最終的に、最初の時刻に撮られた画像の圧縮率が他と比べて圧倒的に悪かった(2. の予測方法を使えないため)です。 なので、3. の部分をどうにか出来ないかなぁというので結構苦労してました。

3.予測誤差の相関除去

予測誤差を眺めてみると、誤差の乗り方も似ていることがよくあるので、予測誤差どうしもブロック内で差分をとります。 差分のとり方は

  1. 横方向に差分を取る
  2. 縦方向に差分を取る
  3. 斜め方向に差分を取る
  4. そもそも差分を取らない

を全通り試して、ブロック内で差分とった結果の絶対値和が一番小さいやつを選んでいます。 1〜3は、ジグザグする感じの訪問順序になります(伝われ)。 ブロック単位で2bit + 最終的な差分の情報を伝えます。

4. 2.3.の情報+3で残った誤差を算術符号で圧縮

まず前提として、ビット列は01の出現頻度の偏りが大きいほど圧縮しやすいです。直感的には、1000桁のビット列を覚える時、①全部0のデータと②01がランダムに混じるデータを比べると、前者が覚えやすい、位で十分と思います。

今回は、 付加情報として送付したビット情報を全部別々のビット列とみなして算術符号化にかけました(実際にビット列を分割したわけではなく、01の発生頻度の確率だけ分けた)。ビット列毎に01の出現頻度は、直前に受け取ったビット列の出現頻度を元に、動的に計算しています(適応的算術符号とかでググると色々出てくると思います)。直感的には、圧縮時と復元時に同じパラメータを使って同じ手順で処理を再現できるなら、特に正確な固定値の出現頻度値を使わなくてもOKという話です。

ビット列が生成された由来や状況に応じてビット列を別物として扱う(=コンテキスト分類する)ことで、01の偏りを強めることができる場合があります。が、今回は特に強いコンテキスト分類まで考えられていませんでした。

コンテスト中に考えていたけどやり残したこと

  1. ある波長画像を別の波長画像にマッピングする処理
    • 大局的に見て線形性がなかったので、参照画像を圧縮するときの予測誤差の大きさに応じて、ピクセル単位で画素位置をN分割して、各分割領域ごとに輝度値と輝度値の対応表をそのまま送ろうと思ってました。が、バグがとれず…
    • 画素位置を分割せずに、入力輝度値に対応する出力輝度値のヒストグラムを見ると、2山になってた。特に深く調べてなかったけど、海と陸で取れる波長の傾向が違うのかなーなんて妄想をしていた。画素位置を分割したらほぼ1山になってたから深くは追わなかったけど、実際のところは知らない。
  2. 予測方法1. で、平均を取らずに誤差最小の係数を線形回帰して計算
    • 予測方法2. 3. の方が圧倒的に効果があると思っていたので、結局手付かずになってしまった
  3. 一旦最適な予測方法が決まった後で、予測方法1. を使うブロックに限定して改めて線形回帰用の係数を計算し直す
    • 予測方法2. 3. のブロック内画素位置を回帰対象に含める必要がない
      1. と同じく、画像間の重ね合わせ部分が最重要だと思ってたので、手付かずになってしまった。
  4. 誤差最小ではなく、符号量最小の選択肢を選ぶ
    • 厳密な測定をしようと思うと、様々な状態に依存するので正直無理ですが、近似的には出来ます。最大サイズの画像は時間的にかなり厳しいですが、小さい画像の最適化ではやるチャンスがあったかもしれないです。
  5. 画像の拡大にBilinearを使っていたけど、もうちょいましな拡大方法を使っても良かったかも

余談

  • 圧縮率と関係無いですが、エンコーダとデコーダって同じ順番で付加情報なり誤差を読まなきゃいけないんだけど、どうにかして共通化する方法ないのだろうか(これのせいで、デバッグ時間的に微ロスを重ねていたので)

感想

  • 前職でチョットやってたので、いつか作りたいと思ってた。実際動くものが作れてからはテンションが5割増くらいだった。
  • いろんな処理を無限にバグらせたし、思ったことが一部できずに終わってしまったので、無限に実装力が欲しくなった。精進しなきゃ。
  • データの圧縮って個人的にすごい好きなんですが、「圧縮っていいよね!!」て言ってる人たちが少なからずいたようで、何か嬉しかった。
  • 今年からマラソンの参加記を書こうという気になったので、これからは終わってすぐ公開できるように準備しておきたいと思います。よろしくお願いします。