・・・なんと前回の記事から２年が経っている。何もしていないわけではないが記事にまとめる気力がなかったのだ。少しまとまった時間ができたので久しぶりに書こうと思う。
今回のソースはこちら。
1. 素のフォトンマッピングの課題
フォトンマッピング法による画像生成は従来の手法（パストレーシングなど）に比べかなり少ないレイ数でも間接光や集光模様を綺麗に表現できるとても素晴らしい手法だと思う。とはいえフォトンマッピング法も完璧ではない。ある点の色（輝度）を求めるため、有限のフォトンと有限の領域（面積）から放射輝度推定を行っているからである。完璧な画像を得るためには、無限個のフォトン、無限小の領域で推定すればいいのだが、現実的には不可能である。そのため、次のような課題がある（と思う）。
- フォトンデータをメモリに読み込む必要があり大量のフォトンを使った計算をするには大量のメモリを持つコンピュータが必要。
- 放射輝度推定に必要なフォトンを抽出するためkd木というデータ構造を使っていて比較的効率は良いようだが、フォトン数が増えると抽出にとても時間が掛かる。
- フォトン数を少なくすればノイズ（色の滲み）が多くなる。これは推定に使うフォトンが少ないので推定値と真の輝度との差が大きくなるため。
- ノイズを減らそうと思うと、より広い領域から多数のフォトンを集めれば良いが、余計なフォトンも推定に使われ、また輝度値が周りと平均化され全体的にぼやけた画像になる（影の縁が鋭くならない）。
このような課題に対応する手法はいろいろ提案されており、直系（？）の研究ではプログレッシブフォトンマッピング法がある。これはさらに改良されてStocastic Progressive Photon Maapingとなっているそうだ。以下は参考URL。
ほかにも、いろいろ研究されているようだ。今回はこの中から Probabilistic Approach という手法を選択した。理由は次の通り。
- 放射輝度推定に使う領域サイズ（参照半径というらしい）を小さくしながら通常のフォトンマッピングを実行するだけ＝フォトンマッピング処理をこのアルゴリズム用に修正する必要がない。
- 参照半径は繰り返し（イテレーション）毎に変えるが、i番目の繰り返しでの半径は数式で求められるため、前後のイテレーションとは独立して実行できる。そのため複数のコンピュータによる並列処理も容易。
- 1回のイテレーションは少ないフォトン数で実行するためメモリが少なくて済む（これはProgressive PM全般の特徴）
- 100回、1000回と多数のイテレーションを行うことで副次的な効果が簡単に実装できる。例えばアンチエリアシング、被写界深度、グロス面などだ。
以後、この手法(Progressive Photon Mapping: Probabilistic Approach)をPPM-PAと書くことにする。
2. PPM-PAの実装
-1. 輝度値出力と統合処理
上記の通り、このPPM-PAは従来のフォトンマッピングプログラムを（ほぼ）変更しなくて済むのがウリの一つだ。参照半径を変化させて生成した画像を最後に平均すれば良い。アルゴリズムは論文にも書いてあるが、簡単な擬似コードで書くと以下の通り。
r_1 ← 初期の半径
for i = 0 to N { # Nはイテレーション回数
参照半径r_iでフォトンマッピング法を実行し画像データを保管
i+1番目の参照半径(r_i+1)を求める
}
全画像データの各画素値の平均を計算し保管
今回はちょっとサボってこの統合処理はRubyで記述した（util/iterator.rb)。詳細はコードを見て貰えばわかると思う。ところで従来はPPM形式の画像ファイルを出力する。つまり各画素はRGB 256階調の整数値だ。そのため一定以上の輝度値以上は255にクリッピングされている。PPM-PAでは輝度値を集計して平均する必要がある。でないと何百回、何千回繰り返すと、ほとんどの画素は黒、所々に色が出るといった画像になるので、データを平均すると値が小さくなりすぎて真っ黒な画像が出来上がるのだ。最初そこに気づかず単純にPPMデータを足したため、ほぼ真っ黒の画像となってしまった。。。
ということで、フォトンマッピングプログラムには２点の小改造を実施した。
-
"progressive"フラグを追加する。
- 従来通りPPMの画像を出力したい場合と、PPM-PAのために一旦画像生成結果の輝度値を書き出す場合を分けるため、"progressive"フラグをスクリーン設定ファイルに追加。値は"yes"か"no"である。
-
progressiveフラグに応じて出力をPPMか輝度値で出力する。
- PPM-PAを使う場合は輝度値をそのまま記録するため、PPMフォーマットなら0-255の整数値を書き出すところ、浮動小数点数を文字列として書き出す。以下はその部分。アンチエリアシング処理は別途考えるため、"smooth"関数は呼び出さない。
if (progressive scr) == True
then -- 輝度値を出力
forM_ [0..(V.length image - 1)] $ \i -> do
putStrLn $ radianceToString (image V.! i)
else do
let pixels = V.map (radianceToRgb scr) image
forM_ [0..(V.length pixels - 1)] $ \i -> do
rgb <- smooth tracer scr pixels i
putStrLn $ rgbToString rgb
ここで出力される形式の画像フォーマットは無いが、あとで集計できるように一旦書き出すだけなので拘らない。なお、結果を輝度値で出力することになると、従来の「エセ」アンチエリアシングがうまく働かなくなる。隣り合う画素のRGB値の差でアンチエリアシング処理を行うかどうかを判断していたからだ。これについては後で触れよう。
-2. 評価
■ フォトン数を一定にしたとき
フォトンマッピング法の生成画像の品質はフォトン数による。そこで、総フォトン数を500万個に固定しイテレーションのフォトン数を変えて比較してみた。
| photon数 | iteration数 | 処理時間(s) | 処理時間(s)/iteration | メモリ量(MB)/iteration | 生成画像 |
|---|---|---|---|---|---|
| 10,000 | 500 | 3,632 | 7.27 | 65 |  |
| 50,000 | 100 | 1,034 | 10.34 | 120 |  |
| 100,000 | 50 | 835 | 16.70 | 169 |  |
| 200,000 | 25 | 875 | 34.99 | 359 |  |
| 500,000 | 10 | 1,661 | 166.05 | 360 |  |
そうフォトン数が同じでも処理時間と品質には結構ばらつきが出た。このシーンについては、フォトン数100,000個が最もバランスが取れているように見える。フォトン数が10,000個程度と少ない場合でもある程度のオーバーヘッドがあるようで、1回のイテレーションに掛かる時間は劇的に短いわけではないようだ。
品質面ではちょっと微妙だった。画像を見比べてもあまりわからないかもしれないが、イテレーション回数が多いと右下の集光模様の縁の鋭さがだいぶ違う。10,000個・500イテレーションが最も良い。ただ、色の滲みは一番ざらざらした感じである。一方、100,000個以上だと、滲みはほとんど違いがないので、縁の鋭さと処理時間を考えると200,000個以上にする意味はほぼないと思われる。
■ イテレーション回数を一定にしたとき
次にイテレーション回数を100に固定して比較した。
| photon数/iteration | 処理時間(s) | 処理時間(s)/iteration | 処理時間(ms)/photon | 生成画像 |
|---|---|---|---|---|
| 1,000 | 552 | 5.52 | 5.52 |  |
| 5,000 | 619 | 6.19 | 1.24 |  |
| 10,000 | 669 | 6.69 | 0.67 |  |
| 50,000 | 1,034 | 10.34 | 0.21 |  |
| 100,000 | 1,508 | 15.08 | 0.15 |  |
| 200,000 | 3,407 | 34.07 | 0.17 |  |
当然の結果だが、同一回数のイテレーションならイテレーション毎のフォトン数が多い方が品質が良いに決まっている。ただ、闇雲に増やすとメモリも時間も掛かりすぎる。ここでも100,000個以上になると品質的にはほぼ変わらないように見える。放射輝度推定に使うフォトンはkd木から検索する。アルゴリズムとしては計算量が O(n log n)になるらしいのでnが10,000以上になってくると処理時間はほぼnに比例するように思えるのだが、実際には指数関数的に処理時間が掛かるようだ。フォトン数100,000個まではオーバーヘッドのためかあまり時間差は無いが、それを超えると処理時間が急増する。ちなみに１回あたりフォトン数500万個で計算しようとしたら6時間経っても1イテレーションが完了しないので中断した。今回のアルゴリズムでは参照半径が決まっているため、フォトンが多いほどその領域に含まれるフォトンが多くなるわけで、一個のフォトンを検索する計算量xフォトン数増加量の時間が掛かることが一因ではと思う。
■ 並列処理
今回は統合処理をRubyで書いたが、並列処理（マルチスレッド）もRubyで行った。結論からいうとあまり効果が出なかった。RubyではマルチスレッドにしてもCPU処理は直列に実行されるらしい。
フォトン数100,000個、500イテレーションでマルチスレッド数を変えて時間を計測した。
| 並列度 | 1 | 2 | 4 | 10 |
|---|---|---|---|---|
| 処理時間(s) | 835 | 757 | 637 | 664 |
使用したPCは4コアなので、I/Oがある程度並列に実行されたのか並列度4までは若干時間短縮がはかれているが、それ以上だと逆にオーバーヘッドが大きくなるようだ。
-3. OpenEXR対応
先日図書館で"[digital] ライティング&レンダリング 第３版"という本を見つけたので読んでいたらいろいろ知らないことが書いてあったのだが、そこに High Dynamic Range Image (HDRI) の説明があった。従来はRGB計24bitで画像データを保管してきたが、画像生成時は輝度を計算している。上記の通りPPM-PAのためにわざわざ輝度値（浮動小数点数）で出力するように変えたのだから、せっかくならそのまま保管できた方が活用の幅が広がる。そこで、最終結果はPPMだけでなくOpenEXR形式でも保管できるようにした。util/averager2.rbがそれである。ファイル形式の詳細は以下のURLに詳しく書いてある。
特に最後にサンプルデータが説明されているのだが、バイト列とその意味が全部記載されており大いに参考になった。ひとまずこの形式のデータが作れれば良いので、Single-Part、非圧縮、チャネルはRGB３色のみ、各色は単精度浮動小数点数(32bit)とした。ファイルサイズが大きくなるのでは？とも思ったが、PPMでも各色が３桁の数字だと1ピクセルあたり12バイト、OpenEXRでも先の条件なら12バイトなのでそれ程の差はない。実際、256x256解像度のファイルだと約791KBである。ちなみにImageMagickでJPEGに変換してみると約17KBになった・・・。
なおOpenEXRで保存しておくと、あとで露出などを簡単に変えられてなかなか楽しいこともわかった。
3. PPM-PAの応用
PPM-PAの直接的なメリットは、レンダリング時のメモリ使用量を抑え、かつフォトンマッピング法に特有の色の滲みとぼやけを同時に緩和できることであると思う。しかし、独立したイテレーションを数百、数千回繰り返すアルゴリズムなので、副次的に次の効果も容易に実装できるのが強みである。
- アンチエリアシング
- 被写界深度によるボケ
- 光沢面
- エリア光源の半影（直接光に古典レイトレーシング法を用いた場合）
今回はアンチエリアシングを（簡単だが）実装したので書いておく。
■ anti-aliasing
最初に書いたように、PPM-PA法を実現するために従来のアンチエリアシング処理は諦めざるをえなくなったが、多数のイテレーションを行うことでより良い（本当の？）アンチエリアシングが可能になった。従来と今回の実装の違いは下図のとおり。
従来（左）は一度レンダリングした後の「追加処理」だったので追加するレイ(赤丸)を４つだけにし、レイの方向も偏らないようあえて固定した。PPM-PA（右）では100回、1000回とイテレーションを行うので、ランダムにレイを飛ばしても偏りは少ないだろう。
また、PPM-PAでのアンチエリアシングは「追加処理」ではないので計算コストが増えないのも大きなメリットである。イテレーション毎にレイの方向をランダムにずらしてあげれば良いのだ。ということで、実装も非常に簡単である。ここで r3,r4は画素の中心からのx,y座標のずれ具合（-0.5 〜 0.5）である。
(r3, r4) <- if prflag == True && aaflag == True
then do
r3' <- MT.randomIO :: IO Double
r4' <- MT.randomIO :: IO Double
return (r3' - 0.5, r4' - 0.5)
else
return (0.0, 0.0)
if文で"prflag == True"はPPM-PAを使うこと、"aaflag == True"はアンチエリアシングを行うことを
意味している。どちらかのフラグがFalseの場合はずらさないので0.0を返す。やることはたったこれだけである。
効果のほどは下の画像で比べてみればよくわかる。左がアンチエリアシングなし、右がありだ。右の方は球の上部や天井からの光の斜めになっている境界のギザギザ（ジャギー）が緩和されているのがわかる。計算コストが増えないので、もうこれからはアンチエリアシングありだけで良いだろう。(画像は解像度256x256、イテレーション500回、1回あたり100000フォトン)
| アンチエリアシングなし | アンチエリアシングあり |
|---|---|
 |
 |
4. まとめ
2年以上ぶりの記事となった。その間細かい変更をいくつかやっていて、記事にまとめておかないと備忘録にならないので書いてみたが、今回の記事はほとんど実装面での説明がなく面白味がないかもしれない。ただ、3で挙げたいくつかの重要な「ボケ」などの効果・レンダリング処理のための布石としては取り組んでよかったと思う。次はこれらの効果を追加していこう