パストレ編ってなんだよ、って話ですが、昨日「パストレがしたい」という名前の記事を書きながら肝心のパストレまで到達できなかったので、今回の記事で補完しようと思います。
前回の記事です。
dora119.hateblo.jp
前回は直射光だけを考慮するDirectIntegratorを実装しました。これを拡張してパストレーサーであるPathIntegratorを実装しますが、とりあえず比較しやすいようにコーネルボックスを描画するようにしときました。
思ったよりぽくなって若干嬉しいです。
BRDFのサンプリング
パストレーサーを実装するときにBRDFのサンプリングをする必要があるので、先にこちらの話をしておこうと思います。
おさらいをすると、を分散を抑えつつ計算するには一般的にMultiple Importance Sampling(MIS)という手法によりBRDF側でサンプリングした結果とライト側でサンプリングした結果をいいかんじに合成するのが効果的だが、前回の実装だとライト側のサンプリングで十分なのでBRDFのサンプリングはしなくてもよいという話でした。なので今のところBRDF::f関数とLight::sampleLi関数しか実装されていません。
そこで今回はBRDF::sampleF関数を実装しましょう。BRDF::sampleF関数は反射角を受け取ったときに、入射角をサンプリングし、その角度とともに対応する確率密度を返す関数です。
単純に球面から一様分布でサンプリングしても良いですが、Importance Samplingにおいて重要なことは被積分関数に刺さるようなサンプリングをする必要がある、つまり確率密度関数が被積分関数の定数倍に近くなればなるほど良いということだったのでBRDF::sampleF関数もそのようなサンプリングを目指します。
今回はpbrtにおいてデフォルトで採用されているとなる確率分布を採用します。というのも、レンダリング方程式をよく見ると最後にの項があり、当然の値が1に近いほどモンテカルロ法において結果への寄与が大きいからです。
cosine-weighted sampling
前回の記事ではサンプリングまわりの話は全て回避してきましたが、今回はそうもいかないみたいです。
(は定数)となるようなサンプリングを行いたいわけですが、結論からいうと円から一様分布でサンプリングをし、それを半単球面に投射したものがその分布になるらしいです。
イメージ的には球面上の微小面がで円に投影されることを考えると筋が通っていそうな気もしますが、証明はpbrを参考にしてください。(2D Sampling with Multidimensional Transformations)
これによって、円から一様分布でサンプリングできればいいことになりました。間違えがちですが、これは極座標についてとを一様分布でサンプリングすればいいわけではありません。*1
確率密度関数が(は定数)となり、円全体で積分すると1になってほしいので、単位円の面積よりとなります。ヤコビアンを考えるととなるので、
が分かります。後はこれらの情報をもとにサンプリングを行いたいわけですが、ここではThe Inversion Method(逆変換法)というものによって]の一様分布から求める確率分布にしたがった確率変数をサンプリングしてみようと思います。
逆変換法のアイディアはとても単純で、確率密度関数ではなく累積分布関数、の逆関数を用います。
確率密度関数に対しでしたが、累積分布関数はのことを指します。この関数は当然単調増加し値域は]なので、累積分布関数の逆関数が求められれば]の一様乱数を用いてに基づく確率変数をサンプリングすることができます。
今回の場合に適用すると、をサンプリングするにはの累積分布関数がなのでこの逆関数を、をサンプリングするにはより累積分布関数が、よって逆関数を求めればよいことが分かります。(ここでは一様変数)
をサンプリングできれば後はそれを座標に変換し、半単位球面に投射するのでを計算するだけです。
virtual RGB sampleF(const Vector3f& _wo, Vector3f* _wi, float* _pdf, Sampler& _sampler) const { float r = std::sqrt(_sampler.generateCanonical()); float theta = 2 * Pi * _sampler.generateCanonical(); Vector<float, 2> p{r * std::cos(theta), r * std::sin(theta)}; float z = std::sqrt(std::max(0.0f, 1 - p.at(0) * p.at(0) - p.at(1) * p.at(1))); (*_wi)[0] = p.at(0); (*_wi)[1] = z; (*_wi)[2] = p.at(1); if (_wo.at(1) < 0.0) (*_wi)[1] *= -1; *_pdf = pdf(_wo, *_wi); return f(_wo, *_wi); } virtual float pdf(const Vector3f& _wo, const Vector3f& _wi) const { return (_wo.at(1) * _wi.at(1) > 0.0) ? std::abs(_wi.at(1)) / Pi : 0.0; }
_sampler.generateCanonical()関数は内部的にはstd::generate_canonicalを実行しているだけです。*2
パストレーシング(理論)
いよいよ本丸ですが、まずは理論的な話から入ろうと思います。
前回の手法では直射光のみを考慮していたため、コーネルボックスの画像を見ても球の影となる部分は真っ黒になっていることが分かると思います。これは本来ならば、光源でないあるオブジェクトの点も光源に照らされて光を反射しているため二次的な光源となりえますが、この間接光の効果を無視していたからです。よりリアルな画像生成のためには、間接光に照らされた地点は再び二次的な光源となることで他の地点を照らすという光の伝播をシミュレートする必要があります。
レンダリング方程式を用いて説明すると次のようになります。*3
上のは本来他の地点から放射されたRadianceであるため、地点からの方向にレイを飛ばしたときの衝突地点をと表すことにしたときと書けるはずです。*4さらにをと書くようにすれば、レンダリング方程式は
となります。この再帰的な式は一つを展開することで以下のようになります。
めちゃくちゃ読みづらくてアレなのでとしてまとめると
であると分かります。
同様に無限に展開しつづけると、
とすればレンダリング方程式はとなります。*5*6
このはから発射されたRadianceがどれくらいに届くかをあらわし、throughputと呼ばれるそうです。
この形に落とし込んだことでようやく各に関してはモンテカルロ法が適用できるようになりましたが、まだを求めるには無限和を計算する必要があります。そこで一般的にはRussian Rouletteという手法を用います。
Russian Roulette
基本的なアイディアはモンテカルロ法と変わりませんが、無限和を求めたいときに確率で計算を切り上げ、代わりに残りの和を0とすることを考えます。ただし、計算したときはで割るものとすると、全体の期待値は素直に総和を計算した場合と一致します。
上の総和をある地点で分割してとの足し算と考えてみると、F_1計算時の期待値はとなり、確実にを計算した場合と一致します。これが全ての項に関して言えるので総和の期待値も変わらないそうです。*7
パストレーシング(実装)
ようやく最後のパストレ実装になります。理論的には各をモンテカルロ法により求めれば良いわけですが、計算を省くためには毎回前回の結果を利用することを考えます。というのも、モンテカルロ法を行うときに確率変数をサンプリングしてを計算することになりますが、を増やしていくときに毎回を一から計算するのはあまりにしんどいので、前回の計算値(betaとします)に新しくサンプリングした確率変数を用いてを新しいbetaにすることでサボろうというわけです。*8
また、愚直に実装すると飛ばしたパスがライトに衝突しない限りが0となってしまうため非常に効率が悪いですが、各を計算するときにパスの最後の頂点だけライトのみからサンプリングすることで確実にライトの光を伝播させることができます、これをNext Event Estimationというそうです。
以下、クソ実装です。
RGB li(const Ray& _ray, const Scene& _scene) const { Sampler sampler; RGB l{0.0, 0.0, 0.0}; RGB beta{1.0, 1.0, 1.0}; Ray ray = _ray; for (int bounces = 0; bounces < 5; bounces++) { SurfaceInteraction si; float t; if (!_scene.intersect(ray, &si, &t)) break; auto light = _scene.getLight(); l += beta * estimateDirect(si, light, _scene); Vector3f l_wo = si.worldToLocal(-ray.getDir().normalize()), l_wi; float pdf; RGB f = si.brdf->sampleF(l_wo, &l_wi, &pdf, sampler); if (pdf == 0.0) break; beta *= f * std::abs(l_wi.dot(Vector3f{0.0, 1.0, 0.0})) / pdf; ray = Ray(si.point, si.localToWorld(l_wi)); if (beta.max() < 1.0 && bounces > 3) { float q = std::max(0.05f, 1.0f-beta.max()); if (sampler.generateCanonical() < q) break; beta /= 1.0 - q; } } return l; }
基本的な流れは
- ray変数に入っているレイを飛ばす
- オブジェクトと衝突しなかったら終了
- Next Event Estimationによりbeta(throughput項)と前回実装した直射光を計算するestimateDirectをかけあわせたものをlにたす
- パスの最後の頂点からBRDF::sampleFにより次に進む方向を決め、betaとrayを更新
- パスの長さが4以上となったらRussian Rouletteをはじめる*9
となっています。
beta(throughput)の値が小さいほど全体への寄与率は小さくなり、そのパスはこれ以上伸ばしてもあまり価値がないということを意味しているので、その値に応じてRussian Roulettteで死ぬ確率も高くなるような実装になっています。
間接光もシミュレートされるようになったので全体的に明るく、球の影も薄くなり壁の色を若干反射しています。
というか思ったより反射していて、ライトを明るくしすぎたからかもしれないです。
パストレーシングのデメリットとして画面が全体的にノイジーになりました。サンプリングを工夫することでノイズを減らすことができますが、今回は最短でパストレをしたかったのでこんくらいのノイズになっています。*10
そして、もう一つのデメリットですがパストレーシングは大量にシミュレートして平均をとる必要があるので、この画像でも1ピクセルにつき512回くらいシミュレートしていたと思います。
おまけ
今回の実装は前回と同じリポジトリにあります。なぜか単一ファイルに書いてしまいました。
github.com
レイトレーシング、およびCGは思ったより数学がでてきて結構重いですね。ただ、画像がうまく生成できたときはとても嬉しいです。
僕はひとまずパストレができたので満足しました。
今回省いた重要な事柄
- Stratified Sampling(Stratified Sampling)
*1:そうすると中心部の密度が高くなる
*2:実装したコードはy軸を高さにしてしまったため(*_wi[1]) = zという意味不明コードになった、ここでバグり散らかしていた
*3:は前回出てきませんでしたがオブジェクトが発しているRadianceです
*4:そもそもオブジェクトに衝突しなかった場合はは0とします
*5:簡単のためを反転させています
*6:ここらへん書きながら式変形したので間違っている可能性大です、もし間違えていたら教えていただけるとありがたいですm(_ _)m
*7:自信がない
*8:もちろんこの省略によってノイズが増えるなどの弊害がある、自分的にはに相関関係が生まれることで正しい結果とは微妙に異なる気がするけど分からず...