Haskellでいってみよう

日曜プログラマにも満たないレベルでもHaskellで何かソフトウェアを作りたい!

レイトレーシング(13): Progressive Photon Mapping (ただし本家ではない)

・・・なんと前回の記事から2年が経っている。何もしていないわけではないが記事にま とめる気力がなかったのだ。少しまとまった時間ができたので久しぶりに書こうと思う。

今回のソースはこちら。

1. 素のフォトンマッピングの課題

フォトンマッピング法による画像生成は従来の手法(パストレーシングなど)に比べかな り少ないレイ数でも間接光や集光模様を綺麗に表現できるとても素晴らしい手法だと思う。 とはいえフォトンマッピング法も完璧ではない。ある点の色(輝度)を求めるため、有限 のフォトンと有限の領域(面積)から放射輝度推定を行っているからである。完璧な画像 を得るためには、無限個のフォトン、無限小の領域で推定すればいいのだが、現実的には 不可能である。そのため、次のような課題がある(と思う)。

  • フォトンデータをメモリに読み込む必要があり大量のフォトンを使った計算をするには 大量のメモリを持つコンピュータが必要。
  • 放射輝度推定に必要なフォトンを抽出するためkd木というデータ構造を使っていて比較 的効率は良いようだが、フォトン数が増えると抽出にとても時間が掛かる。
  • フォトン数を少なくすればノイズ(色の滲み)が多くなる。これは推定に使うフォトン が少ないので推定値と真の輝度との差が大きくなるため。
  • ノイズを減らそうと思うと、より広い領域から多数のフォトンを集めれば良いが、余計 なフォトンも推定に使われ、また輝度値が周りと平均化され全体的にぼやけた画像にな る(影の縁が鋭くならない)。

このような課題に対応する手法はいろいろ提案されており、直系(?)の研究ではプログ レッシブフォトンマッピング法がある。これはさらに改良されてStocastic Progressive Photon Maapingとなっているそうだ。以下は参考URL。

プログレッシブフォトンマッピング1

プログレッシブフォトンマッピング2

Adaptive Photon Mapping

Probabilistic Approach

Probabilistic Approach(日本語訳)

ほかにも、いろいろ研究されているようだ。今回はこの中から Probabilistic Approach という手法を選択した。理由は次の通り。

以後、この手法(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データを足したため、ほぼ真っ黒の画像となってしまった。。。

ということで、フォトンマッピングプログラムには2点の小改造を実施した。

  • "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万個に固定しイテレーションフォトン数を変えて比較してみた。

フォトン イテレーション 処理時間(s) 処理時間(s)/イテレーション メモリ量(MB)/イテレーション 生成画像
10,000 500 3,632 7.27 65 f:id:eijian:20200822174619j:plain
50,000 100 1,034 10.34 120 f:id:eijian:20200822174626j:plain
100,000 50 835 16.70 169 f:id:eijian:20200822174845j:plain
200,000 25 875 34.99 359 f:id:eijian:20200822174907j:plain
500,000 10 1,661 166.05 360 f:id:eijian:20200822174920j:plain

そうフォトン数が同じでも処理時間と品質には結構ばらつきが出た。このシーンについて は、フォトン数100,000個が最もバランスが取れているように見える。フォトン数が 10,000個程度と少ない場合でもある程度のオーバーヘッドがあるようで、1回のイテレー ションに掛かる時間は劇的に短いわけではないようだ。

品質面ではちょっと微妙だった。画像を見比べてもあまりわからないかもしれないが、イ テレーション回数が多いと右下の集光模様の縁の鋭さがだいぶ違う。10,000個・500イテ レーションが最も良い。ただ、色の滲みは一番ざらざらした感じである。一方、100,000 個以上だと、滲みはほとんど違いがないので、縁の鋭さと処理時間を考えると200,000個 以上にする意味はほぼないと思われる。

イテレーション回数を一定にしたとき

次にイテレーション回数を100に固定して比較した。

フォトン数/イテレーション 処理時間(s) 処理時間(s)/イテレーション 処理時間(ms)/フォトン 生成画像
1,000 552 5.52 5.52 f:id:eijian:20200822175054j:plain
5,000 619 6.19 1.24 f:id:eijian:20200822175107j:plain
10,000 669 6.69 0.67 f:id:eijian:20200822175124j:plain
50,000 1,034 10.34 0.21 f:id:eijian:20200822175140j:plain
100,000 1,508 15.08 0.15 f:id:eijian:20200822175154j:plain
200,000 3,407 34.07 0.17 f:id:eijian:20200822175206j:plain

当然の結果だが、同一回数のイテレーションならイテレーション毎のフォトン数が多い方 が品質が良いに決まっている。ただ、闇雲に増やすとメモリも時間も掛かりすぎる。ここ でも100,000個以上になると品質的にはほぼ変わらないように見える。放射輝度推定に使 うフォトンはkd木から検索する。アルゴリズムとしては計算量が O(n log n)になるらし いのでnが10,000以上になってくると処理時間はほぼnに比例するように思えるのだが、実 際には指数関数的に処理時間が掛かるようだ。フォトン数100,000個まではオーバーヘッ ドのためかあまり時間差は無いが、それを超えると処理時間が急増する。ちなみに1回あ たりフォトン数500万個で計算しようとしたら6時間経っても1イテレーションが完了しな いので中断した。今回のアルゴリズムでは参照半径が決まっているため、フォトンが多い ほどその領域に含まれるフォトンが多くなるわけで、一個のフォトンを検索する計算量x フォトン数増加量の時間が掛かることが一因ではと思う。

■ 並列処理

今回は統合処理をRubyで書いたが、並列処理(マルチスレッド)もRubyで行った。結論か らいうとあまり効果が出なかった。RubyではマルチスレッドにしてもCPU処理は直列に実 行されるらしい。

RubyでCPUコアをフル活用してみた

フォトン数100,000個、500イテレーションでマルチスレッド数を変えて時間を計測した。

並列度 1 2 4 10
処理時間(s) 835 757 637 664

使用したPCは4コアなので、I/Oがある程度並列に実行されたのか並列度4までは若干時間 短縮がはかれているが、それ以上だと逆にオーバーヘッドが大きくなるようだ。

-3. OpenEXR対応

先日図書館で"[digital] ライティング&レンダリング 第3版"という本を見つけたので読 んでいたらいろいろ知らないことが書いてあったのだが、そこに High Dynamic Range Image (HDRI) の説明が あった。従来はRGB計24bitで画像データを保管してきたが、画像生成時は輝度を計算している。 上記の通りPPM-PAのためにわざわざ輝度値(浮動小数点数)で出力するように変えたのだから、 せっかくならそのまま保管できた方が活用の幅が広がる。そこで、最終結果はPPMだけでなくOpenEXR 形式でも保管できるようにした。util/averager2.rbがそれである。ファイル形式の詳細は以下の URLに詳しく書いてある。

OpenEXR File Layout

特に最後にサンプルデータが説明されているのだが、バイト列とその意味が全部記載されており大いに参考になった。 ひとまずこの形式のデータが作れれば良いので、Single-Part、非圧縮、チャネルはRGB3色のみ、 各色は単精度浮動小数点数(32bit)とした。ファイルサイズが大きくなるのでは?とも思ったが、 PPMでも各色が3桁の数字だと1ピクセルあたり12バイト、OpenEXRでも先の条件なら12バイトなので それ程の差はない。実際、256x256解像度のファイルだと約791KBである。ちなみにImageMagickJPEGに変換してみると約17KBになった・・・。

なおOpenEXRで保存しておくと、あとで露出などを簡単に変えられてなかなか楽しいこともわかった。

3. PPM-PAの応用

PPM-PAの直接的なメリットは、レンダリング時のメモリ使用量を抑え、かつフォトンマッピング法に 特有の色の滲みとぼやけを同時に緩和できることであると思う。しかし、独立したイテレーションを数百、 数千回繰り返すアルゴリズムなので、副次的に次の効果も容易に実装できるのが強みである。

今回はアンチエリアシングを(簡単だが)実装したので書いておく。

■ anti-aliasing

最初に書いたように、PPM-PA法を実現するために従来のアンチエリアシング処理は諦めざるをえなく なったが、多数のイテレーションを行うことでより良い(本当の?)アンチエリアシングが可能になった。 従来と今回の実装の違いは下図のとおり。

f:id:eijian:20200822175351p:plain

従来(左)は一度レンダリングした後の「追加処理」だったので 追加するレイ(赤丸)を4つだけにし、レイの方向も偏らないようあえて固定した。 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フォトン)

アンチエリアシング無し アンチエリアシング有り
f:id:eijian:20200822175440j:plain f:id:eijian:20200822175453j:plain

4. まとめ

2年以上ぶりの記事となった。その間細かい変更をいくつかやっていて、記事にまとめておかないと 備忘録にならないので書いてみたが、今回の記事はほとんど実装面での説明がなく面白味がないかも しれない。ただ、3で挙げたいくつかの重要な「ボケ」などの効果・レンダリング処理のための布石としては 取り組んでよかったと思う。次はこれらの効果を追加していこう。