ラベル 3DCG の投稿を表示しています。 すべての投稿を表示
ラベル 3DCG の投稿を表示しています。 すべての投稿を表示

2014年3月27日木曜日

フォトンマップの実装

イントロダクション

いままで本ブログでは、3DCGのレンダリング手法であるレイトレーシングとパストレーシングの記事をいくつか書いてきました。レイトレーシングはフォンシェーディングなどのレンダリング手法からすれば、反射・屈折などが扱え、格段に写実的になり高品質な画像を合成できるようになっています。ただし、レイトレーシングにも弱点がありました。フォンシェーディングはリアルタイムレンダリングに使える手法ですが、レイトレーシングでは、リアルタイムに描画するのは難しいでしょう。また、間接光や集光模様などの影響が扱えません。全体的に暗めの画像が生成され、影も違和感があります。写真のような画質とはとても言えません。そんなに時間をかけずに、そこそこの品質の画像を生成するのにはよいでしょう。写真のような画質の画像を生成するには、すべての光の軌道を計算する必要があります。現実的にその計算を行うことは不可能ですが、それに近いことを行うのがパストレーシングです。パストレーシングは、全ての光の軌道を計算する代わりに、その中のいくつかの光の軌道をランダムに選び、全体像を推定する統計的な手法です。選ぶ光の軌道の数が少なければ、正確に推定することができずに、生成された画像にムラができてしまいます。なので、十分な数だけ光の軌道を選ぶ必要があります。こうして生成された画像は、写真とは見分けがつかないほどリアルなものになります。ただ、十分な数の光の軌道を追跡するのは、かなり時間がかかるのが欠点と言えるでしょう。本記事では、それらの欠点を補うために、フォトンマップを使った手法を取り上げます。


図1. レイトレーシングとパストレーシングの比較

パストレーシング

フォトンマップに行く前に、パストレーシングの問題点を振り返っておきましょう。パストレーシングでは、まず視点(カメラ)からスクリーン上のあるピクセルに向かって光線を放ちます。その光線はシーン中の視点にもっとも近いポリゴン上のある1点と衝突します。その表面の材質によって光線のその後の軌道が変化します。鏡面や屈折面だった場合は、正確にその後の軌道が計算されます(フレネルの法則)。拡散面(ざらざらした面)だった場合は、光は乱反射します。乱反射するということは、その後の光の軌道は、様々な方向へ光は分散していくということです。


図2. スクリーンのあるピクセルに放った光線


図3. 鏡面反射と乱反射

ちなみに、分散レイトレーシングでは、拡散面に衝突した時は、実際に複数(例えば、100本だったり1000本だったり)の光線を追跡します。ただ、この方法だと拡散面に衝突する度に、100本や1000本の光線が生まれ、結果的に膨大な光線(100のn乗や1000のn乗、nは拡散面との衝突回数)を追跡する羽目になります。

一方、パストレーシングでは、光線が拡散面に衝突した時、考えうる反射方向の中から、1本だけランダムに軌道を選びます。したがって、視点から放たれた光線は、最初から最後まで1本です。拡散面では、光線の反射方向はランダムに選ばれるので、状況によって光線の反射方向は毎回違います。つまり、スクリーン上のあるピクセルの色は、視点から放たれた光線の経路によって変わってくるので、もしかしたら、ある時は赤、またある時は青といったように状況によって(乱数生成環境によって)変わってきます。

この状況を現実世界に置き換えて考えてみると、光源から放たれた1つの光(光子)が、いろんな経路を通って目に届くことを意味しています。パストレーシングは(レイトレーシングも)、光源から来た目に届くはずの光を、逆に目から光源に向かって追跡しているのです。別々の経路を通った光は、様々な材質の表面を通ってくるので、やはり、それぞれ色(波長)が違うものになるはずです。ただ、現実世界では、光源から放たれた光線は、たったの1本だけが目に届くわけではなく、複数の光線(それこそ無限近いの光線)がいろんな経路を通って目に届くはずです。目に届く光は、無数の光線の合計の色として認識されます。

パストレーシングも現実世界の例に習い、光線の方向は逆方向ではありますが、視点から1つのピクセルに向かって複数の光線を放ちます。ある光線は、いくつかの反射・屈折の後、光源まで届くでしょう。また、ある光線は、十分な回数の反射・屈折を経た後でも、光源に到達することはないでしょう。実際に、光源から目に届く複数の光の経路というのは、視点から放った光線の内、光源に到達した光線の経路と同じになるはずです。つまり、パストレーシングは、十分な数の光線を使用すれば、ほぼ完全に現実世界をシミュレートすることができます。

しかし、十分な数の光線を使用しなかった場合は、合成された画像には、ノイズが生じます。このノイズを説明するために、サイコロを振ったときに出る目を例にあげましょう。例えば、サイコロを5回振ったとき、出た目が順に、3、2、2、1、2だったとしましょう。この場合、出た目の平均は、(3 + 2 + 2 + 1 + 2)/5 = 2 になります。もし十分な回数だけサイコロを振っていれば、おそらく平均値は3.5になるはずです。ある程度の回数まで試行すれば、平均値が3.5からあまりズレなくなるでしょう。ですが、試行回数が少ないと平均値が3.5から大きくずれることが多いでしょう。 パストレーシングにおいても、同様のことが言えます。視点からあるピクセルに放出する光線の数が少ないと、実際の色から大きく外れた色になる確率が高くなります。光線の数は十分に大きくとらないといけないのです。


図4. ノイズの生じたパストレーシング

パストレーシングの最大の難点は、この部分にあると言えるでしょう。十分な数の光線を使用すると、描画に膨大な時間がかかります。しかも、いくつかの場合は、ほとんどの光線が光源に到達することがありません。光源に到達しない光線というのは、無駄に時間を消費してしまっており、パストレーシングによるレンダリングの時間がかかる原因の1つと言えるでしょう。もっと効率的に、計算する必要がありそうです。

この問題を解決しようと試みるのが、本記事の主題であるフォトンマップを使用したレンダリング手法です。この手法を使えば、パストレーシングに比べ画像生成は格段に早くなります。ただし、近似計算的要素が含まれているため、完全な現実世界の再現にはなっていない部分もあると思いますが、その影響は小さくすることができます。

フォトンマップ

パストレーシングの問題点を解決するには、拡散面で反射した光線が光源に到達するまでの経路探索を、より効率よくすることです。そこで、拡散面上のある点Pにおける光の入出力に注目して、人の目に届く光の量を考えてみましょう。

点Pを経由して目に届く光の量は、光源から点Pに届く光の量に依存しているはずです。ここで、点Pに関する光の量を2つに分類します。1つめは、光源から点Pに届く光の量。2つめは、点Pから目に届く光の量。まず、2番目に分類した光の量について考えましょう。つまり、点Pに届いた光の量のうち、どの程度が目に届いているのかを考えます。完全な拡散面の場合は、点Pに入射した光の方向に依らずに、どの方向にも一定の割合で反射するので、点Pからは一部の光だけが目に届きます。また、目に届く光の量は、点Pから目までの距離に依存します。例えば、点Pに入射した光子の個数が100個だったとすると、反射して出て行く光子の数も100個です。そして、この100個の光子は、この後どこにいくのか? 簡単の為に、この光子は1秒間の間に1m進むものと仮定すると(本当は秒速30万キロ)、1秒後には、100個の光子は、点Pから1mだけ離れた場所にいるはずです。ある点から定距離だけ離れた点の集合が作る曲面は、球面です。なので、100個の光子は1秒後には、半径1mの半球面上(半球なのは、光子は裏側へは反射しないから)に一様に分布しているはずです。この時の球面の面積は、球面積の公式を使って、4 x π x 1 x 1 = 4πの半分の2πです。光子の数の密度は、100/2πとなります。点Pの光の密度の1/2πになります。光の明るさは、目に入ってくる光子の量に比例するので、点Pのすぐ近くにいる人と点Pから1mだけ離れている点にいる人が感じる光の明るさの比は、2π:1となります。では、rメートルだけ離れた人が感じる明るさは、どの程度でしょうか? 点Pから射出された100個の光は、rメートル進んだ時には、半径rメートルの半球面上に一様に分布しているので、rメートルだけ離れた場所の光の密度は、100/4πr^2になります。これが光の明るさです。つまり、光の明るさは、距離の二乗に反比例して弱くなっていきます。余談ですが、重力や電磁気力(そもそも光は電磁場の振動)などの力にも同じことが当てはまります。これらの力も同じような理由(この場合はフォトンの密度ではなくて、球面状の場の密度)で距離の二乗に反比例します。

次に、光源から複数回の反射・屈折を経て点Pに入射してくる光の量を考えてみましょう。パストレーシングでは、この量を点Pから光源に向かって追跡して求めました。前述したように、点Pからどの方向へ光線を飛ばせば、複数回の反射・屈折を経て、光源に衝突するのか分からないので、ランダムに光線を飛ばします。どこにも反射せずに、直接光源に届く経路を見つけるならまだしも、複数回反射・屈折して光源に到達する経路を見つけるのは、ランダムな方向に光線を射出していたのでは、偶然に衝突するまで待たなければなりません。あまり効率が良くないのです。しかも、現実に近い輝度を推定するには、十分な数の光線が光源に衝突しなければなりません。その十分な数の光線を得るには、膨大な時間がかかります。この点が、パストレーシングの難点でした。

パストレーシングのような方法とは別に、点Pに入射してくる光の量を推定する良い方法はないでしょうか。この問いに対する一つの解決案は、とても単純です。点Pから光源に至る経路を探索するのではなく、現実世界と同じように、光源から光線(光子=フォトン)を追跡すればよいのです。光源からフォトンを追跡すれば、全てのフォトンが輝度をもっているので、シーンのどこかで衝突すれば、それは意味のあるものになります。パストレーシングでは、放った光線が光源に衝突するのを待つのに比べたら、だんぜん効率がよさそうです。また、シーン全体にどのようにフォトンが分布しているかを知ることができるので、一度フォトンの分布を調べておけば、点Pの輝度を推定するだけでなく、他の点の輝度も推定できることがパストレーシングの場合と違います。パストレーシングでは、物体上の各点から、ランダムに光線を放って、光源までの経路を推定しなければなりませんでした。いろいろな場所で同じことを繰り返さなければならないのです。逆に光源から追跡する方法だと、一度だけシーン全体のフォトンの分布を計算すれば、輝度を推定したい他の場所でも使い回すことができます。この点は、大幅な計算量の削減に貢献するでしょう。このようなシーン全体のフォトンの分布はフォトンマップと呼ばれています。では、なぜ今までこの安直な方法(誰でもまず最初に考える)を試してこなかったのでしょうか?それは、光源からフォトンをランダムに放っても、点Pに届く確率がとても小さくて、非常に効率が悪いからだと思います。点Pには、ほとんど面積がないので、フォトンが点Pとほとんど衝突しません。これでは、点Pに入射する光量を推定することが困難です。ただ、これが点Pの周辺まで含めれば話は変わります。点Pの周辺も点Pとだいたい似たような輝度だろうという大胆な仮定を取り入れるならば、点Pの周辺(点Pを中心とする球)に入射する複数のフォトンから点Pの輝度を計算することができます。ただし、点Pの周辺は大体同じような輝度だという仮定は、常に成り立つわけではありません。壁と床の境目などの場所では急激に輝度が変化するので、そのような場所にはこの方法は、うまく働きません。


図5. 光源からフォトンを放つ・点Pの周辺のフォトンを収集

とりあえず、ここまでの方法をまとめましょう。
1. 光源からフォトンを放ち、シーン全体のフォトンの分布フォトンマップを構築する。
2. 視点からスクリーン上のあるピクセルに向かって放った光線と拡散面との交点をPとする。
3. 点P周辺のフォトンを収集し輝度を推定する。
4. 点Pの輝度から目に届く輝度を計算する。(2πr^2で割る)

フォトンマップの構築

物体上のある点(とその近傍)に入射する光の量を求めるには、フォトンマップを構築しなければなりません。フォトンマップは、シーンの中の拡散面に衝突する度にフォトンの情報を保持しています。保持する情報としては、衝突した位置、フォトンの色などです。フォトンの発生源はシーンにある全ての光源です。光源の種類によって違いますが、基本的に物理的に可能な方向へランダムに発射します。たとえば、平面光源の場合は、光源上からランダムに1つの点を選び、その点を中心に、天頂角が180度以内の方向に(半球を形作るように)ランダムにフォトンを発射します。一度拡散面に衝突したフォトンは、フォトンマップに格納され、そこで追跡をやめるのではなく、その後もフォトンを追跡していきます。

フォトンマップのデータ構造にも、注意を払う必要があります。というのも、ある点に入射する光の量を計算するときに、その点の近傍にあるフォトンを、フォトンマップの中から収集しなければならないからです。ただの配列にフォトン情報を詰めてもいいですが、ある点の近傍にあるフォトンを見つけるには、配列にある全てのフォトンとの距離を計算する羽目になります。探索が容易なデータ構造を選ぶのよいでしょう。ちなみに、今回、僕は八分木を使いました。


図7. 光源からシーンに放たれたフォトンの分布

レンダリング

フォトンマップを構築したら、レイトレーシングでフォトンを収集していきます。視点からスクリーン上のあるピクセルに向かって放たれた光線は、拡散面に衝突するまで反射屈折を繰り返します。拡散面に到達したら、その点の周辺のフォトンをフォトンマップから収集して輝度を計算します。点の周辺が意味するところは、さまざまあると思います。基本的には、その点を中心とした球の内側に入るフォトンを輝度の計算に使います。例えば、ある点中心にして、フォトンが100個収まるような球の使うという方法があります。100個が収まるには、ある点では、大きな球が必要かもしれませんし、別の点では、小さな球で済むかもしれません。いずれにせよ、各点で球の半径が変わります。そして、100個のフォトンが半球面を通過して、点Pの周辺に到達したことになります。基本的に、その点の入射輝度は、単位半球面を通過するフォトンの合計です。なので輝度を計算するには、球の大きさが単位半球だったら、フォトンがいくつ入射するかを推定しなければなりません。それは100個入射したときの球の表面積と単位球の表面積の比から求めることができます。例えば、100個のフォトンが収まる球の半径が10だったとしたら、表面積は4 x π x 10 x 10 = 400πで、単位球の表面積は、4πです。したがって、単位球の表面積は、半径が10の球の表面積の1/100になります。おそらく、入射するフォトンの個数も1/100になり、同時に輝度も1/100になります。なので、その点の入射輝度は、100個のフォトンの合計輝度の1/100です。これが各点における入射輝度の計算方法です。

別の入射輝度の計算方法として、各点で球の特定の数のフォトンの個数を収集できる大きさの球を使用するのではなくて、半径を固定して、その球の中に入ってくるフォトンで輝度を計算することもできます。どの点でも、同じ大きさの球を使用しているので、球に入って来たフォトンの合計の輝度の比がそのまま、各点の輝度の比と等しくなります。生成された画像に多少ムラができやすかもしれませんが、こちらの方法の方が簡単です。

先ほども述べましたが、壁と床の境なのどが輝度計算に使用する球の中に含まれてしまうと、正確に輝度を計算することができません。床のある点の輝度を計算しているのに、壁に到達したフォトンも輝度計算に入れちゃうからです。この現象を避けるには、ある点の近傍の範囲を決めるのを球を使わずに、円盤上のものを使うのがよいでしょう。輝度を計算したい点があるポリゴンの法線方向に球を潰したような立体です。こうすることで、法線と直交する平面(ポリゴン)に到達したフォトンを収集できる割合が高くなり、関係ない物体に届いたフォトンの収集を少なくすることができます。

フォトンが楕円体の中にあるか外にあるかの判定は、楕円体の中心とフォトンまでの距離と楕円体の中心からフォトンがある方向の楕円の表面までの距離を比較することで判定します。中心から表面までの距離は、
a * b / sqrt( b * cos(θ) * b * cos(θ) + a * sin(θ) * a * sin(θ) )
で求めることができます。aは長軸の長さ、bは短軸の長さ、θは法線と中心からフォトンまでの線分のなす角です。


図8. 円盤状の球体でフォトンを収集する


図9. 楕円を使わなかった場合、面と面との境付近で余分なフォトンが含まれている

このようにして収集したフォトンから、入射輝度(色)計算します。1つのフォトンが運ぶ光の量は同じなので、入射輝度は、単純にフォトンを足し合わせることで計算できます。入射輝度が計算できたら、その点から視点に届く光の量を計算します。この点は拡散面上の点なので、届いた光は一様に、様々な方向に反射して生きます。その一部が目に届くのです。この場合の目に届く光の強さは、先ほど述べたように1/2πr^2です。

あとは、スクリーン上の全てのピクセルに大して同じことを繰り返すことで、画像を生成することができます。

レンダリング結果

フォトンマップを使った方法では、光源から放射するフォトンの量によって、生成される画像の精度が変化します。その変化の様子の例を以下に示します。


100個フォトンを使用し、半径が0.5の場合。


1000個フォトンを使用し、半径が0.3の場合。


10000個フォトンを使用し、半径が0.1の場合。


100000個フォトンを使用し、半径が0.3の場合。


100000個フォトンを使用し、半径が0.1の場合。


500000個フォトンを使用し、半径が0.1の場合。


1000000個フォトンを使用し、半径が0.2の場合。

2014年3月17日月曜日

モンテカルロ積分を使ったレイトレーシング

本記事では、表面積分を統計的に行って光源からの寄与を計算するレイトレーシングを行います。積分を統計的処理で計算する手法をモンテカルロ積分というらしいです(本記事がそれに該当するかは不明。ただ、似たようなことをやっていると思います)。解析的に(数値計算や近似を使わずに式変形だけで正確に)積分計算をするのが難しい積分値を計算する時の手法です。ランダムに値を生成することから、カジノで有名なモンテカルロの名前を与えられたと聞いています。前回の記事では、スクリーン上のあるピクセルの色を計算するとき、光源からの寄与は、ある場所(輝度を計算したい場所)から光源を見たときの視野角に比例すると仮定して計算しました。今回は、もっと正確に計算するために、ある場所を中心とする単位半球面上に投影される光源の面積を統計的に求めます。つまり、光源がありそうな方向一帯に光線をいくつもランダムに放ち、光線が光源と衝突した回数と衝突しなかった回数との比から面積を求めます。もちろん、放つ光線の数が多いほど正確に面積を求めることができます。

面積の統計的推定

輝度を計算したい場所(視点から放たれた光線と視点に最も近いポリゴンとの交点)から見て、光源がありそうな方向の決定には、前回計算した視野角を使用します。ここで言う視野角は、図のように、視点から光源の外接円(または外接球)が収まる角度です。つまり、ある場所と光源の中心とを結ぶ線分と単位半球面との交点を中心として、半径が視野角の半分となるような単位半球面上に描かれる円の内側のみに光線を放ちます。円の内側のみ放つというのは、この円の面積は解析的に求めることができるので、光源がないとわかっているところに光線を放っても意味ないからです。今回は100本の光線を放ち面積を推定しました。100本のうち、ある光線は、光源にぶつかり、ある光線は、光線にぶつからないでしょう。そして、何本が光源に衝突するかがわかると、衝突した回数を100(放った光線の本数)で割れば図の単位半球面上の円の中に占める光源の割合が分かります。球面上の円の面積は、解析的に求めることができるので、その面積と先ほど求めた割合をかければ、求める光源の面積が計算できます。もちろん、光線の本数が多ければ多いほど正確に面積が求めることができるわけです。100本とか200本とか1000本とかで割合を計算してみて、結果にばらつきがないくらいの本数だったら十分だと思います。単位球面上の円の面積は、sinθdθdφを0≤θ≤'視野角/2', 0≤φ≤πの範囲で積分した値になります。
ここまでの過程で、ある場所の光源からの輝度が計算できます。輝度は、先ほど求めた光源の面積に比例するはずです。あとは、視点からスクリーン上の全てのピクセルに向かって光線を放ち、視点に最も近いポリゴン(ポリゴンでなくともよいが)との交点について、それぞれ同じ計算を施せば、1つの画像が生成されることになります。当然、前回に比べて計算量が増えるので、描画処理は遅くなります。

結果比較

前回の結果と比較してみます。前回は、単位半球面上に投影される光源の面積は視野角に比例すると仮定して輝度を計算しました。以下の最初の図(左側の図)が、今回行った統計的に面積を求めた時の結果です。次の図(右側の図)が前回の結果です。今回の結果の方が壁面の影のエッジが丸っこくなっている感じがします。おそらく、今回の結果の方が正しい結果なはずなんだけど、果たして実際にそうなんだろうか。現実は、間接光の影響があるので、直接光だけの影響を受けた物体の写真あるはずもないので確かめられないのだが。うーむ。

2014年3月1日土曜日

レイトレーシングの拡散面の輝度計算

3次元グラフィックの描画方法の一つであるレイトレーシング。光の反射・屈折などが正確に計算できる手法で結構キレイに画像生成ができます。ただし、間接光が計算できないので、暗めの画像になってしまいます。本ブログでも、Scalaで実装したレイトレーシングで描画した画像を載せていましたが、実は今までレイトレーシングよって生成される画像の各ピクセルの色の計算を結構適当に計算していました。今回は、この色の計算について考えてみたいと思います。

レイトレーシングの復習

簡単にレイトレーシングの仕組みについて復習してみましょう。現実の世界で、光源から発せられた光が目に届くまでには、いろんな物体に何回も衝突し、その一部が目に届きます。写真のような画像を作りたければ、光源が発せられた全ての光の道筋を計算すればよいことになります。でも、光源からの光の道筋を全て計算しても、最終的に目に届く光というのは、ごく一部です。ほとんど無駄な計算になります。コンピュータの使えるリソースは、限られているのに無駄な計算を何日もかけて計算する気にはなりません。もっと効率的に計算できる方法はないか。ちょっとぐらい写真画質じゃなくてもいいよ。と思った人が、たぶんレイトレーシングという効率的な方法を考えだしたんでしょう。

光源から光を追跡したら、膨大な計算量を必要とする上に、目に届くのは一部。じゃあ、目に届いた光だけを逆に追跡すればいいんじゃね。そう考えたんですね。目に届いた光というのは、スクリーン上のある点と目のある点を結ぶ半直線です。つまり、光源からの光を追跡するんじゃなくて、目がある位置からスクリーン上のある点を通過する半直線(光線=レイ)を追跡(トレース)しようってわけです。これがレイトレーシング(光線追跡)です。光源からの光は無数にありますが、目からの光線は、目がある場所からスクリーン上の各ピクセルとを結ぶ半直線のみなので、ピクセル数の光線を追跡するだけで済みます。

これで、目に届く光だけを計算することができます。目から出た光線は、シーンの中にあるいろんな物体にぶつかって、材質によって反射したり屈折したりして、ある光線は、光源にぶつかり、ある光線は、どこにもぶつからずに闇に消えていきます。光源にぶつかった光線は、その光線が通過したスクリーン上の点を明るくするでしょう。どこにもぶつからなかった光線は、その光線が通過したスクリーン上の点を明るくすることはありません。このようにして、各ピクセルの色が決められていきます。なので、物体表面の反射や屈折を計算して、どの光線が光源とぶつかるかを計算するのが重要になってきます。光線の反射や屈折の計算は、鏡や透明なガラスなどの材質だったら簡単です。では、ざらざらな表面の場合は、どうでしょう?ざらざらした表面は、光が乱反射します。ちなみに、このような性質の表面を拡散面と呼ぶようです。拡散面では、光線は、いろんな方向に反射します。鏡みたいに一方向に反射しません。ん?ちょっと待てよ。色んな方向に反射する?ということは、無数に分散する光線を追跡していかなきゃならんのですね。これじゃ、光源から発せられた光線を追跡するのと変わらんやん!だめやん。振り出しに戻ったやん。と思ってしまうかもしれませんが、レイトレーシングでは、拡散面にぶつかったら、そこで光線の追跡終えます。これがレイトレーシングが正確にピクセルの色を計算できないとこです。では、どうするのかというと、拡散面にぶつかった光線のピクセルの色は、拡散面の位置と光源の位置から近似計算します。間接光の計算はしません。間接光を計算するには、拡散面で分散した光線を追跡しなければならないからです。ちなみに、この追跡を統計的に行うのがパストレーシングです。ということで、レイトレーシングは、拡散面は近似計算で鏡やガラスなどの材質は、正確に計算するレンダリング手法です。

拡散面に届く光の量と目に届く光の量

先程も、拡散面では近似計算して、ピクセルの色を計算すると言いました。ここでは、具体的にその計算方法を考えてみます。以下のピクセル値の計算は、僕が勝手に近似計算しているだけなので、よりよい方法があるかもしれませんし、間違っているかもしれません。あしからず。

拡散面から来る光の量を計算を簡単にするために、光の量の種類を2つに分類します。1つは、光源から拡散面に届く光の量で、2つめは、拡散面から目に届く光の量です。

拡散面に届く光の量の近似計算

光源から拡散面上のある点に入射してくる光の量を正確に計算するには、その点を中心にした単位半球面を考えます。半球面は、その点での法線方向側にあります。つまり、拡散面の表に半球面があります。その単位半球面を横切る光線の数が、この点に入射してくる光の量です。これは、入射光線ベクトルを半球面に沿って積分することと同じことです。半球面に入射してくる光線ベクトルは、光源から直接入射してくるのもあれば、何回も反射してきて、入射してくるのもあります。前者の光線は直接光、後者の光線は間接光です。先程も言及したとおり、この入射ベクトルを全て計算するのは、大変なのでレイトレーシングでは、この積分計算を正確に計算せずに近似計算します。
まず、間接光は計算しません。間接光を計算するには、何回も反射した光線を考慮しなければならないので、非常に計算量を要するからです。直接光だけを考えましょう。半球面上を横切る直接光の数は、半球面上に投影される光源の面積に比例するはずです。半球面に投影される光源の面積は、図にあるように、魚眼レンズで見た時の光源が形作る図形の面積のようになると考えるとイメージしやすいでしょう。それでは、この面積を計算しましょう、と言いたいところですが、この積分計算もややこしそうなので、思い切って端折って、この直接校の計算も近似することにします。
では、どのようにして近似しましょう。たぶん、半球面上に投影される光源の面積は、図のように、拡散面上の点から光源を見た時の角度にだいたい比例するんだと思います。きっと。

以下にその角度を計算するコードを示します。
  val v0 = p -> center // 拡散面から光源の中心までのベクトル
  val d = radius // 光源を囲む円の半径
  val nt = normal * (v0 * normal) // v0の光源の法線方向成分
  val vr = v0 + (nt -> v0).normalize * d // 拡散面から光源の端までのベクトル
  val dt = vr angle v0 // 拡散面から光源を見た時の角度の半分

目に届く光の量の近似計算

次に、拡散面から目に届く光の量を考えます。拡散面に入射した光の量の一部が目に届きます。拡散面なので、この面に入射した光が全て目に届くはずがありません。そのごく一部の光が目にとどきます。正確には目に届く光は、ちょうど目のある位置の方向に偶然反射した光だけです。したがって、拡散面で反射する光が、どの方向にも一様に反射すると仮定すると、目に届く光の量は、拡散面に入射した光の量を半球面の面積で割った量だけになるはずです。ここで半球面の半径は、拡散面上の点から目の位置までの距離です。

それに対して、鏡で反射して目に入ってくる光は、鏡のその点に入射した光が全て同じ方向(目がある方向)に反射した光です。なので、この場合は、半球面の面積で割る必要がありません。

描画結果

レイトレーシングで描画した結果です。この画像をみる限り、光の量の計算は、よい近似になってるのではないでしょうか。

2014年2月20日木曜日

フォトンマップのための八分木の実装


フォトンマップとは?

今までの記事では、パストレーシングを利用してレンダリングしてきましたが、一番の問題は、なんといってもレンダリングに時間がかかりすぎることです。フォトンマップを使ったレンダリングは、この問題を解決してくれます。フォトンマップとは、光源から放射されたフォトン(レイ、光線)が拡散面に衝突した位置(と入射方向や反射方向)とその衝突面の材質を格納するデータ構造です。このフォトンマップを使えば、拡散面における輝度推定が容易になります。つまり、フォトンマップを使ったレンダリング過程は、大きくわけると2段階にわけられます。第1段階は、フォトンマップの構築です。第2段階は、レイトレーシングまたはパストレーシングを使って、拡散面の輝度推定値を構築されたフォトンマップから計算してレンダリングします。

フォトンマップのデータ構造

拡散面に衝突したフォトンの情報が格納されてれば、どんなデータ構造を使おうが問題ありません。ただ、輝度推定の過程で、格納された全てのフォトンの中から目的のフォトンを効率よく見つけ出す必要があるので、この目的に合ったデータ構造を使えば、より短い時間でレンダリングできます。フォトンマップによる輝度推定では、拡散面上の対象の点に近い位置にあるフォトンを集めて輝度を計算します。このような計算に対して、フォトンを効率よく見つけるデータ構造をいろいろ考えだされているようですが、ここでは、パストレーシングに比べてかなり早くなっていれば、他のデータ構造に比べてちょっと遅いくらいは問題ないので、比較的単純な構造である八分木を実装することにします。二分木の3次元版みたいなものだから、イメージしやすいし。

八分木の仕様

八分木は空間を再帰的に8分割して、その子空間にデータを格納する木構造です。この木構造をフォトンマップとして使えるように、次のような仕様にしました。
  • 空間として直方体を使う
  • 直方体の内側に位置するフォトンは、その直方体に格納される
  • 1つの直方体はフォトンに対して容量を持っていて、フォトンを格納できる数が決まっている
  • 直方体のフォトンが容量に達したら、直方体を8分割し、該当する子の直方体にフォトンを格納する

八分木の実装

八分木クラスの主要なメソッドを概観します。

コンストラクタ

コンストラクタでは、直方体を定義します。
  abstract class Octree[B, C <: Octree[B, C]](
    val xmin: Float, val xmax: Float,
    val ymin: Float, val ymax: Float,
    val zmin: Float, val zmax: Float,
    val capacity: Int,
    val parent: Option[C]
  )
型パラメータBは、格納するデータの型です。型パラメータCは、サブクラスの型です。子の直方体を生成する際には、C型のインスタンスを生成するようにします。これは、サブクラスで定義したメンバーを使用できるようにするためです。

子を生成するメソッド

直方体を8分割して子の直方体を生成するときに使用されるメソッドです。コンストラクタと同じ引数を取ります。サブクラスで実装します。
    protected def create(
      xmin: Float, xmax: Float,
      ymin: Float, ymax: Float,
      zmin: Float, zmax: Float,
      capacity: Int,
      parent: Option[C]
    ): C
戻り値の型は、上述したようにサブクラスの型です。

フォトンの追加

フォトンを追加するメソッドです。既に最大容量に達していれば、子を作成し、子に格納する。
10 
11 
    def add( a: Vector3, b: B )(implicit tag: ClassTag[C]): Unit = if( ndata < capacity ) {
      _data = (a, b) :: _data
      ndata += 1
      onUpdate( b )
    } else {
      createChildren // 既に子を持っている場合は子を作成しない
      // このOctreeが持っているデータを子に移す
      _data.foreach { case (p, q) => findChild( p ).foreach( _.add( p, q ) ) }
      findChild( a ).foreach( _.add( a, b ) )
      _data = Nil
    }

フォトンの探索

指定した球に含まれる、または接している最下層の直方体を見つけ、クロージャに渡す。
10 
11 
    def filterChildren( b: BoundingSphere, f: (Int, C) => Unit ): Unit = boundingSphere.contains( b ) match {
      // 離れている場合
      case 3 => ()
      // bがこのOctreeを含んでいる場合
      case 2 => f( 2, this )
      // このOctreeがbを含んでいるか接している場合
      case n => _children match {
        case None => f( n, this )
        case Some( cs ) => cs.foreach { _.filterChildren( b, f ) }
      }
    }

2014年2月13日木曜日

魚眼レンズの実装

大域照明で3次元モデルをレンダリングする場合、モデル中のある点における輝度を計算するために、図1のようにその点を中心とする単位半球面を通過する光の量を求める必要があります。つまり、その点からモデル全体を見渡した時に、その点から見て光源はどの方向にありどのくらい入射してくるのか、また光源じゃない部分からもその点に間接的に光は届くので、その部分からの光の量はどのくらいなのかということを計算しなければいけません。前回のパストレーシングは、ある点に光源から直接届く光と光源から1度別の点で反射して届く光のみを考慮して輝度を計算しました。したがって、完全な大域照明と言えないでしょう。

図1 ある点に入射する光。球面に沿って面積分して光の総和を求める

モデル中のある点に入射してくる光の量というのは、言い換えれば、その点を中心とする単位半球面に沿った光の量の面積分のことだと言えるでしょう。パストレーシングは、その面積分を行っているわけですね。つまり、直接光に関して言えば、半球面上に投影されるモデル中の光源の面積が直接光の量に比例します。であるからには、ある点からモデル全体を見渡した時、その半球面上に描かれるモデルを視覚化してみたくなりました。そこで、魚眼レンズのようなものを実装し、魚眼レンズに移った像を平面に描いてみました。球面を平面で表現する方法は、地球を地図で表現する方法と同じようにすれば良さそうですね。学生のころ地理の授業でいろいろならった記憶があります。メルカトル図法とか。ここでは、ランベルト正積方位図法を使って2次元に映し出してみることにしました。ランベルト正積方位図法っていうのは、面積比が実際と比べてだいたい同じくらいになる図法です。正確に面積比が同じなるわけではないようです。ですが、簡単に球面上の面積を視覚化するには、ちょうどいい方法だと思い、この図法を使うことにしました。

平面(スクリーン)上のピクセル値から球面座標を求める

球面上の点をスクリーンに投影するには、球面上の点とスクリーン上のピクセル値を結び付けなければなりません。図2に示すように、ランベルト正積方位図法では、球面上の点P'がスクリーン上に投影する点P''は、天頂Aと点P'が作る線分AP'を半径とする円とスクリーン2と交わる点のうち近い方の点とします。なので、スクリーン上のピクセル値を決めれば、その点に対応する球面上の点を決めることができる。


図2 ランベルト方位図法

球面上の点は、極座標で表現するのが簡単なので、天頂角AOP'を求めるのがよいでしょう。この天頂角AOP'は、余弦定理で計算できます。そんな定理ありましたね。懐かしいですね。ほとんど忘れてるのでWikipediaで調べてみました。
cos(x) = (b^2 + c^2 - a^2) / (2bc)
三角形の三辺の長さが分かれば、角度が分かるよっていう定理ですね。
ここで三角形AOP'を見ると、辺AP'は分かってます。でも、辺AOと辺OP'は分かってないですね。ただ、辺AOと辺OP'は円の半径OBなので同じなので、求めるべきは円の半径OBということになります。

ところで、点Bは、魚眼レンズの一番端の部分です。この点をスクリーン上で考えると、やはりスクリーンの一番端になります。スクリーンが正方形だと仮定して、その幅をwとすると線分AB'の長さは、w/2となります。線分AB'と線分ABの長さは等しいので、線分ABの長さもw/2になります。三角形AOBは、直角二等辺三角形なので、各辺の長さの比AO:OB:ABは、1:1:sqrt(2)となります。つまり、辺AOの長さは、辺ABの長さを1/sqrt(2)倍したw/(2 x sqrt(2))です。先ほど述べたように、線分AOと線分P'Oの長さは等しいので、線分P'Oの長さもw/(2 x sqrt(2))です。

これで、三角形AOP'の全ての辺の長さが分かりました。余弦定理に当てはめると天頂角AOP'を求めることができます。次のソースコードは、Scalaでスクリーン上のピクセル値から半球上の点の座標を求めるコードです。

10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
    /** スクリーンの幅 */
    private lazy val s = screen.ymax - screen.ymin

    /** スクリーンの半径 */
    private lazy val sr = s/( 2f * sqrt( 2f ) )

    def position( x: Int, y: Int ) = {
      // ピクセル値からシステム座標系に変換
      val sx = screen.screenX( x )
      val sy = screen.screenY( y )
      // OBの長さの2乗
      val ss = s * s
      // AP''の長さ
      val rr = sx * sx + sy * sy
      val r = sqrt( rr )
      // 天頂角AOP'
      val theta = acos( 1f - 4 * rr / ss )
      // OP'のxy成分
      val xy = sr * sin( theta )
      // 点P'の座標
      if( x == screen.width/&& y == screen.width/) Vector3( 0, 0, sr * cos( theta ) )
      else Vector3( sx * xy / r, sy * xy / r, sr * cos( theta ) )
    }

半球面の点P'に投影されるモデル中の点を求める

スクリーン上のピクセル値から球面上の点P'の座標は求めることができました。では、この点P'にはモデル中のどの点が投影されるのでしょうか? この計算はさほど難しくありません。なぜなら、ベクトルOP'とシーンの中の全てのポリゴンとの交点のうち点P'と最も距離の近いものが、求める点になります。これは、レイトレーシングの時に登場する計算と同じです。

魚眼レンズを使ってレンダリング

前回のシーンを実装した魚眼レンズを使ってレイトレーシングでレンダリングしました。比較の為に普通のカメラを使ってレンダリングした画像を載せておきます。左が普通のカメラ。右が魚眼レンズ。



拡散面から見たモデル全体

下左図の赤丸の中の白い点から魚眼レンズでシーン全体を見渡してみました。下右図は、レイトレーシングのレンダリング結果です。


以下はその時のコード。
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
package examples {

  import net.ruffy.marble.graphics2d.{ Color }
  import net.ruffy.marble.graphics3d.{ Screen, DefaultCamera, LambertCamera, RayTracer }
  import net.ruffy.marble.math.{ Math, Vector3 }
  import Math._

  object LambertMap {

    def main( args: Array[String] ) {

      val renderer = new RayTracer {
        val width = 300
        val height = 300
        val screen = Screen( width, height )
        val root = cornelBox
        val camera = {
          val cam = DefaultCamera( screen, 17f, 88f.toRadians, 90f.toRadians, 5 )
          root.toCamera( cam ).nearestIntersection( cam.createPhoton(width/8, height/2) ) match {
            case Some( i ) =>
              val wp = cam.toWorld( i.point )
              val wq = cam.toWorld( i.point + i.normal )
              LambertCamera( screen, wp, wp -> wq, 0f )
            case _ => LambertCamera( screen, 5f, 88f.toRadians, 90f.toRadians )
          }
        }
        val numberOfTimesToTrace = 5
      }

      renderer.render
      renderer.screen.write( "png""images/raytrace_lambert.png" )

    }

  }

}

2014年2月3日月曜日

パストレーシングの収束2

前回は、数字でパストレーシングの収束を考えました。今回は実際にレンダリングして画像のノイズの減り具合を目で比べてみました。光源から拡散面に1度反射してカメラに到達したレイ(直接光)と光源から拡散面に2度反射してカメラに到達したレイ(環境光)を使って輝度を計算しました。環境光として2度拡散面に反射してカメラに到達したレイしか考慮してないためか、全体的に暗めな印象になってます。各ピクセルの輝度計算に使用したレイは、30000 paths/pixelくらいまでです。かなり時間かかりました。下の画像のように、10000 paths/pixelくらいでノイズが目立たなくなって来ています。

2014年1月17日金曜日

パストレーシングはいつ収束するのか?

前回はパストレーシングで光源からの物体に届く光の量を推定しました。各ピクセルから500本の光線を飛ばしましたが、ノイズが残っていて滑らかに明るさが変化していませんでした。そこで今回は、各ピクセルにどのくらいの光線を飛ばせば、ノイズが減り、滑らかに明るさがへ変化するかを調査してみました。

どのように調査するのか?

画像全体で調査するとかなり時間がかかるので、とりあえず、前回の図の一部分(下の図の赤い線に囲まれた部分)だけの明るさに注目して、その部分の隣り合うピクセル同士が滑らかに変化しているかを調査することにします。

ヒストグラムにして比べてみる


この赤い部分の各ピクセルに届く直接光の本数を左のピクセルからカウントした結果をヒストグラムにしました。ヒストグラムは4つあります。これら4つの違いは各ピクセルから発射した光線の本数の違いです。それぞれ発射した本数は、順に1000本/ピクセル、10000本/ピクセル、30000本/ピクセル、50000本/ピクセルです。各ヒストグラムの横軸は、画像の場所を表しています。縦軸は、そのピクセルに到達した光源からの直接光の本数です。例えば、最初のヒストグラム(1000 photons/pixel)の横軸が60のとこの縦軸の値は60くらいを示しています。これは、赤い部分の左から60ピクセル目には、直接光は60くらい届いたよってことを表しています。横軸の0~10の間は、縦軸が0になっていますが、これは赤い部分の左端の真っ黒な部分を表しています。

では、各ヒストグラムを比べてみます。最初のヒストグラム(1000 photons/pixel)は、ガタガタしていて決して滑らかであるとは言えません。これは各ピクセルの輝度が滑らかに変化していないこと示しています。つまり、画像にノイズが多いってことですね。1000 photons/pixel程度では、きれいな画像は生成できないのでしょう。次のヒストグラム(10000 photons/pixel)は、射出した光線数が10倍になっているだけあって、結構滑らかになっています。最後のヒストグラム(50000 photons/pixel)になれば、かなり滑らかになっていることがわかります。これくらい光線を飛ばせば、きっときれいな画像を生成できるんじゃないでしょうか(かなり時間はかかりますが)。

いつ収束するのか?

では、各ピクセルから飛ばす光線の数によって、それぞれ生成された画像がどの程度滑らかになっていくかを数値化してみます。こうすることで、どの程度の光線数で充分かがある程度推測しやすくなります。

滑らか度を数値化すると言っても、どういった方法で数値するのが正しいのか分からないのですが、ここでは、隣のピクセル同士の直接光から光線数の差の割合dyの平均値と、dyの値自体にばらつき(標準偏差)によって数値化することにしました。ヒストグラムの横軸をx、縦軸をyxとすると、dy、平均値μ、標準偏差σは、次のようになります。
dyx = |yx - yx-1|/yx
μ = Σxdyx/n
σ = √1/(n-1)Σx(m - dyx)2
平均値μ、標準偏差σは、1つのヒストグラムにつきそれぞれ1つづつ計算されます。上のヒストグラムは4つしか載せてないですが、実際には各ピクセルから飛ばす光線の数は、1~50000本まで1000本づつ増やしていきました。なのでヒストグラムは50個生成されています。この50個のヒストグラムに対して、平均値μ、標準偏差σを計算し、それをグラフにしたのが下の図です。

横軸は各ピクセルから飛ばす光線の数、縦軸は平均値μ、誤差棒は標準偏差σです。光線の数が多くなるにつれて収束していってるのがわかります。30000本あたりからあまり変化がなくなってきているので、これ以上光線を増やしても、時間がかかる割に生成された画像は、それほどキレイにならないんだと思います。

まとめ

この記事では、パストレーシングに充分な各ピクセル数から射出する光線数は、
30000本
とします。
でも、これじゃ画像生成するのに時間かかりすぎだよなぁ。そろそろフォトンマップに移った方がいいのかも。

関連記事
Scalaでパストレーシング実装中

2014年1月8日水曜日

[3DCG] Scalaでパストレーシング実装中

最近は3DCGのプログラミングをしています。Scalaを使ってます。まだまだ実装途中ですが、ここまでの結果をログとして記録しておこうと思います。

今はパストレーシングを実装中。光源から光の届き具合を調べる為に、直接光が当たる部分と間接的に光が当たる部分の画像に分けて出力しました。ピクセル毎に500パスを使用して描画しました。以下は、その結果の画像です。ちなみに、まだ色は考えてないので、白黒の画像です。輝度のみです。

直接光による寄与

各ピクセルから500本の光線を飛ばし、拡散面で1度だけ反射して光源に到達した光線を描いた画像。手前の球に光源が映っているところがおかしい。ノイズがある。光線の本数が足りないのかなぁ。


間接光による寄与

直接光の画像と同じように、各ピクセルから500本の光線を飛ばしました。違いは、2度拡散面で反射し光源に到達した光線を描いたところが違います。全体的に光が分散してる感じなので悪くないのでは!?


直接光+間接光

上の2つの画像を足しあわせた画像。うーん、いまいち。何かが足りない感じ。参考にしている本によれば、100paht/pixelくらいで結構キレイな画像生成しているだけどなぁ。


あとは、これをもうちょっと改造して色をつけてって感じです。最終的にはフォトンマップまで実装しようと思ってます。