PathTracingでのCausticsレンダリング手法
初めに
この記事はレイトレ合宿8アドベントカレンダー7/17の2個目の記事です。
パストレーシングにおいて効率的にCausticsをレンダリングすることができるManifold Next Event Estimation(MNEE)という手法とそれの発展手法である Specular Manifold Sampling(SMS)を先行手法のManifold Exploration(Manifold Walk)と共に大まかに解説する記事となっています。
- 初めに
- CausticsとPathTracing
- PathTracingでCausticsを出す手法
- Manifold Exploration
- Manifold Next Event Estimation
- Specular Manifold Sampling
- 実装の資料
- おわりに
- 参考文献、引用論文
CausticsとPathTracing
Caustics(集光模様)とは複雑な曲面による反射、屈折によってできる光の模様のことです。Causticsは日常生活でも簡単に見られ、水の入ったグラスやペットボトルに光を当ててみると非常に複雑で美しい光の模様(Caustics)を見ることができます。ほかにもきれいな海の底にはよくCauticsが現れていて、絵的に映えるエフェクトの一つです。
実際に水が入ったグラスやペットボトルに光を当ててみるとCauticsが出てきます。
Causticsが生まれる流れを見てみましょう。ガラスを通した光というのは屈折によりその経路を曲げられてしまいます。ガラス面が複雑だと例え平行光を当てたとしても、屈折によってさまざまな方向へと光が曲げられてしまうわけです。そうすると地面には光が集まりやすい場所と逆に光が全然集まらない場所というように光の偏りが生まれ、結果としてCauticsとして我々は光の模様を見ることになります。すごい簡単なイメージとしては虫眼鏡で光を集めるようなものです。
Cauticsは水やガラスなど屈折する物体が支配的なシーンだと絵的な部分で非常に重要なエフェクトになります。なので是非ともレイトレでこのCauticsを出したいところなのですが、実は問題があります。凹凸が激しい理想反射面が存在するシーンを用意し、純PathTracingを1000spp程度で行うと以下のようにファイアフライが強く現れます
これはBRDFからサンプリングする方向がCausticsの経路と全くの無関係であるためです。レイがスペキュラー面を通った上で”偶然”光源に当たることでしか寄与を得ることができません。しかも、間にスペキュラー面が入ることで寄与が高い経路が指向性を持つことになるため、このようなひどいノイズが生まれるわけです。明示的に光源とサンプリング点をつなげる手法としてNEEがありますが、通常のNEEだとまっすぐにパスをつなげることになるのでこのような屈折反射を考慮したパスというのは生成できません。このように単純なPathTracingでは実用的にはCauticsを出ことができません。
Cauticsは視点からの探索では難しいのですが、一方で光源側から探索する場合には簡単です。そのため、光源側からもパスを作るBDPT(BiDirectional Path Tracing)や光源からフォトンを飛ばすフォトンマッピングはこのCauticsを綺麗に出すことが可能な手法として知られています。
#raytracing #cpp #3dcg #pbr #rendering
— yumcyawiz🍵 (@yumcyawiz) 2021年12月19日
自作レンダラーにVolumetric progressive photon mappingを実装して水入りコーネルボックスをレンダリングしてみた。Causticsと水の質感がめちゃくちゃ綺麗に出てくれて楽しい pic.twitter.com/3BumMWhUkS
これらの手法は静止画においては非常に強力な手法ですが、アニメーションの分野に入ってくると話が変わります。この手法をアニメーションのレンダリングを行うとメモリや目立つノイズなどのアニメーション特有の問題点が現れ、これら手法をそのままアニメーションレンダリングに使うことは難しいとされ[Hanika et al. 2015]、多くのアニメーションを前提とする商用レンダラーではあまり採用されていません(多分)。対照的に PathTracingはアニメーションレンダリングにおいて安定したレンダリングを行うことができ、多くの商用レンダラーではPathTracingを採用しています。こうした理由から動的なCausticsレンダリングを行うには、新しくPathTracingの枠組みの中で効率的にCausticsを出す手法を見つけ出す必要があります。
PathTracingでCausticsを出す手法
話はちょっと変わって最近Blenderがバージョンアップして3.2になり、目玉機能としてShadowCausticsと呼ばれる機能が追加されたとのことでした。これは何だろうとみてみたところなんとCausticsを出すことができるものでした。色々とShadowCausticsに関する動画はいっぱいありましたが、見た感じちゃんとしたCausticsを出すことができるっぽいです。
実際にやってみると以下のようになりました。同じシーンで大体同じ秒数でShadowCausticsなしとありで比べています。
こうしてみるとちゃんと実用的な形でCausticsを出していることがわかります。BlenderのCyclesはPathTracerのはずなので遂にPathTracerでCausticsが出せるようになったんだと感じました。
こんなことがあり、どうやってこんなん実装してるんだろうということで調べたところ[Hanika et al. 2015]"Manifold Next Event Estimation"(MNEE)という手法を使っているということがわかり、レイトレ合宿もあるので折角なので実装目当てで周辺論文含め色々調査してみました。今現在もまだいまいち実装まで漕ぎつけていませんが、大まかな内容をまとめるため記事にしたという次第です。
結構元論文[Hanika et al. 2015]は先行論文を読んでいることが前提のような書き方であり、実装も(なぜかBlenderも含め)見当たらず中々探るのは難しかったのですが、MNEEの改良を行った[Zeltner et al. 2020]"Specular Manifold Sampling"(SMS)という論文ではMNEEを含め先行研究をまとめており、なおかつ実装などの資料が豊富であったためこちらの論文を主に探りつつほかの論文を調べていきました。
そういう感じでMNEE及びSMSの周辺論文をまとめると
- [Jakob Marschner et al. 2012] Manifold exploration
- [Jakob 2013] Light Transport On Path-Space Manifolds
- [Kaplanyan et al. 2014] The natural-constraint representation of the paath space for efficient light transport simulation
- [Hanika et al. 2015b] Improved Half Vector Space Light Transport
- [Hanika et al. 2015a] Manifold Next Event Estimation
- [Zeltner et al. 2020] Specular Manifold Sampling
というような流れで作られていると思われます。特にManifold explorationについてはMNEE及びSMSの根幹的な手法であるため、この記事ではManifold explorationを含めMNEEとSMSについて簡単にまとめて行こうと思います。
Manifold Exploration
Manifold Explorationは[Jakob Marschner et al. 2012]にて発表された手法です。この論文では純モンテカルロ法ではなく、MCMC(マルコフ連鎖モンテカルロ法)で扱う形になっていますが、ここで提唱されたManifold Exploration自体は後々の論文で基礎的手法として使われています。ちなみにJakobさんの博論[Jakob 2013]はこれをより詳細に述べているので調べる際はこちらの方が良いかと思います。
この手法はSpecular Manifoldと呼ばれるものを定義し、ガラスなどを通るような現在のパスを少し動かし、物理的に成立する新たなパスを発見する方法です。
Specular Manifold
図のようなLDSDEパスを作ったものとして、スペキュラー面は理想的な反射、屈折を行うもののみ(glossyはなし)を考えます。この際、スペキュラー面Sではから来た光はに当たり、反射することでに向かいます。この反射は当然皆さんが知る通り、スペキュラー面の法線との入射角、反射角が等しくなるように設定されています。
当然と言えば当然なのですが、これは(理想反射の)スペキュラー面を通るパスに必ず課せられる条件となります。これを逸脱することは物理的にあり得ない光の軌道となってしまいます(フェルマーの原理に違反する)。
この制約から反射先であるというのはから決定されることから、このパスというのはの部分を取ったより低い次元で考えることができます。従って、スペキュラーが入っているパス空間はスペキュラーの制限が入ったより低い空間、Specular Manifoldで考えることができるというわけです。
なので、(あくまでスペキュラー面を通る)上記の例のパス空間は以下のようなSpecular Manifoldに落とし込むことができます。次元はからへとなっているはずです。(各点は何かしらの表面つまり多様体上にあるため、座標の自由度は2になる)
[Jakob Marschner et al. 2012]ではこのSpecular Manifoldを一般的に定義し、そこから近傍のパスの探索へと発展させています。
一般の場合を考え、長さnのパスがあるとして、各頂点をとしておきます。各スペキュラーS上の点には先ほどの反射、屈折の法則から前の点と次の点の間には位置座標に制約がかけられます。この制約をある関数として表すこととして、の制約を
と関数を使って表現します。この制約関数はが反射、屈折の条件を満たしていれば0,満たしていない時0以外を返します。こうしておくことで、スペキュラーを通るパスというのが制約関数の解として表現できるというわけです。
それはよいのですが、実際このような制約関数はどうやって作るのかというとハーフベクトルと法線の関係を用います。[Jakob Marschner et al. 2012]では以下のように制約関数を定義しています。
は単にのハーフベクトルを返す関数で二番目の式のように作られています。nはそれぞれ境界面の屈折率で屈折におけるハーフベクトルも同様に扱えるようになっています。
要としてはですが、これは点の接空間のtangentとbinormalの基底ベクトルで作られた接平面への投影行列です。具体的な表現としてはのtangent基底を,binormalの基底を、ハーフベクトルをとすれば
となっています。なので制約関数は二次元ベクトルとなっています。図的なイメージとしては以下のようにtangent,binormalで作るtangent planeにハーフベクトルをが落とした影のベクトルがというわけです。
接空間は法線が一つの直行基底ベクトルなわけですから、tangent plane上に投影すると0ベクトルになることがわかります。従ってが以上の投影で0ベクトルになるということは法線と平行であることを示し、ちゃんと制約関数は(1.1)を満たしていることになります。
ここで長さのパスの内、p個の頂点はスペキュラー面Sに存在するとします。この時、スペキュラー面にある各点はおのおの制約関数の条件を満たしている必要があり、これをパス全体の制約としてパスの制約関数を考えることができます。ということで一般的なSpecular Manifoldは以下のように定義することができます。
長々とSpecular Manifoldの定義をやってきましたが、このような定義にしたのはSpecular Manifoldのパラメトライズを行うためです。
パス制約関数は実質的に各頂点の陰関数として表現されており、その陰関数が作る多様体がSpecular Manifoldとなっています。また、ジオメトリ上の点の集まりでしかないためパス制約関数は微分が可能です。 このことからパス制約関数に対して陰関数定理を適用でき、ある特定のパス近傍でパラメトライズが保証されることになります。これは何が言いたいかというとある頂点を使ってほかの頂点をの関数としてとして表現できることが保証されていることを示しています(多分)。
今あるパスの頂点をちょっと動かした新しいパスが欲しいなというときにパス制約関数をちゃんと満たすパスの他の頂点を数値的に計算することが可能である、つまり近傍のパスの探索を行うことができるわけです。
SpecularManifold上ではDiffuseの点が決定すると、Specular上の点はそれらの点から自動的に決まります。このため、このパラメトライズは個のDiffuse点からSpecular点を求める関数とみなすことができるわけです(ただし現在のパス近傍限定)。
の微分は多様体の接空間を与えており、これはの微分から得ることができます。そのため、の微分を計算していく必要があります。
具体的な例
もう少し具体的な状況として図のようなパスを考えます。 このパスはと各頂点がつなげられており、その内からはスペキュラー上に存在します。このようなパスに対してパス制約関数はそれぞれの各点の制約関数を集めた行列として考えることができます。
このに対して微分を考えます。の中身は各頂点の微分が入っており、その微分は以下のような5×7行列として書き下すことができます。
各で依存しない頂点の微分を0にすれば、
というように書くことができます。この微分はSpecular Manifoldにおいて様々な情報を持ち、2 ~ 6列というのは特にスペキュラー面に関する情報を持ちます。ということでDiffuse点の微分である1,7列目部分を、Specular点の微分であるで行列を分けることとします。
これらを用いると、SpecularManifoldのtangent spaceを以下のように求めることができます。(ここの理由はちょっとわかんないです)。
これは、の行列で各列はDiffuse点に関する各点の微分の値を示しています(の微分とも取れる?)
この結果は、が少し動いたとき、はどれだけ動くのかということを(1.8)によって求めることができることを示しています。すなわちちゃんとSpecular Manifoldのパラメトライズができていることになります。ちなみにについて逆行列を求めていますが、は式(1.7)のように三重対角行列行列のため、実装にそこまで計算の負担にはなりません。
Manifold Exploration
以上によってSpecular Manifoldのパラメトライズなどをお話していきましたが、結局のところ何がしたいのかというと、こうしたパラメトライズをすることでパスの探索に役立てることができるというわけです。
例えばですが、終点のDiffuse点ではなく、特定の点にパスの終点を移動させたい場合を考えます。この場合、通常のパストレだとに狙って当たるようにのサンプリングを変えることはできないので、この問題は難しいものとなります。
しかし、今まで考えてきたManifold Specularとの値から我々はに関するパラメトライズができるようになっています。をに移動させた時、Specular上の点はどのように移動すればいいかは(1.8)より、
として求めることができるわけです( を挟んでいるのは、あくまでこのパラメトライズが"tangent空間でのパラメトライズ"のためで))。これにて新しくできた方向にからサンプリングするとに到着するパスが出来上がるわけです。
このようにSpecular Manifoldを使うことで新たなパスを効率的に求めることができ、この手法をManifold Explorationと呼びます。論文内ではアルゴリズムそのものはManifold Walkと呼んでたりします。
それでは実際にManifold Walkの疑似コードを見ながら、どのような処理でこれができるのかを見ていきます。
の頂点があるとして、終点であるはDiffuse点で他の点はすべてスペキュラー上の点として考えています。このManifold Walkは終点をという新しい点へと移動させるということを考えています。
はイテレーションであり、はManifold Walkにおけるステップの度合を管理するものでこれが大きければステップにおける移動量が大きくなります(初期値は1)。
- Line2~3
2よりManifold Walkの処理が開始します。このWhile文は更新していったの位置がターゲットに非常に近くなった時に成功とみなし、回行ってもこれを満たさなかった場合は失敗という形でループを構成しています。
まず、Line3にて複雑な式で新しい点を生成していますが、一つ一つほどくとの意味が分かります。まず、マイナス部分と合わせての部分というのはこれは式(1.8)と(1.9)から""を意味しています。つまりは何らかの頂点を移動させた結果なわけであります。
今回、Ideal Specularのみを考慮しているため、Specular上の点は1つでも決まればそこから明示的に他が決定するため、についてのみの移動さえ考えればいいわけです。ということで先ほどの移動量からに関するもののみを取り出すという行列をかけて、の移動量を求めているわけです。この移動量はタンジェント空間におけるものなのでワールド空間上への移動量にするためを挟んで最終的に必要なを得ているわけです。最後に移動量にβをかけていますが、単に移動量を調節するために入れているだけです。
以上のことから、が意味することは"点がへと移動すると考えられる新しいの位置"であることがわかると思います。
- Line4
ただし、の計算はあくまでManifold Specular上で行われているため、実際のGeometry上の点ではありません。なので正しい値にするため、から方向にレイを飛ばし、新たなパスをLine4で生成しています。
この新しいパスの終端はおおよその場合で完璧ではないもののに近づきます(Manifoldの計算のため)が、必ずしもというわけではなく、移動量が多すぎてを通り過ぎてしまう場合があります。この場合分けをLine5,Line8で行っています。
- Line5 ~ 9
がにちゃんと近づいている場合は現在のパスの値をへと更新して、同様の探索を再び行っていきます。
一方でがから離れてしまうということはこれは移動量が多すぎて、がを飛び越してしまう状況です。なので移動量の度合の値を小さくし、再び同じように探索を行わせることで、別のパスを調べていくことになります。
以上が大まかなアルゴリズムの説明で、こうすることでSpecular Manifold上で新しいパスをある程度指定して探索することができます。[Jakob 2013]ではあくまでこの探索法をMCMCで使用していましたが、後年ではPathTraceにこの方法が応用されMNEE,SMSへと使われていくこととなります(多分)。また、この手法はの解を数値的に求めているニュートン法に近く、後の論文ではこの手法をニュートン法として述べているものもありました。
Cの微分の計算
Manifold walkを行うのに一番重要なのはの微分の計算で、これをどうやって計算するかという問題がまだ残っています。例として、最も簡単なSpecular chainを考えます。の頂点があり、はDiffuse点、はSpecular点とみなします。この際のは点に関するとなります。
この時のの微分は、以下のように示されます。
について、式(1.2)より書き下すと各微分について依存する変数で以下のようになります。
この時、それぞれの点のTangent Spaceの基底をそれぞれU,V(UV座標とかではなく任意の基底です)と考えます。それぞれの頂点には基底U,Vの位置微分、法線微分を得られるとします。それぞれの意味はU、V方向にx_iを少しだけ動かした時の移動量と法線の変わり具合を示しており、実際のパストレで簡単に得られるとされています。
次に各頂点をU,V空間上の座標(u,v)で表すこととします。上に存在するの座標というのはを用いて、以下のように表すことができます。また、その場所の法線も同様に表すことができます。の時は頂点の位置となります。
以上の論議は頂点をTangent Space上の座標(u,v)として表現しようということで、Tangent Space上でちょっとだけ動いたら(=微分)がどうなるかを表すのに便利なのでこのようにしたわけです。ちなみになどが得られることを前提としていましたが、実際の実装でも例が当たったジオメトリから求められるのでこの辺は問題ありません。
では、の微分部分についてみていきましょう。については特に微分には無関係なので置いておいて、ハーフベクトル部分の計算をしていきます。先ほどのパラメーター化ではとなりますので、ハーフベクトルの微分というのは(u,v)各座標での微分として分けることができます。
ハーフベクトルの計算は実際には
となりますので、ベクトルをに関する微分公式
より求めることができます。はの位置が求められるため、この式を求めることができ、と併せると2×2の行列として表現することができます。のtangetnSpaceの基底を とすると、
と書き下すことができます。この式はでも同様であり、この式によっての両端は計算することができます。
一方でに関する微分はの微分に加え更に、の微分も入ってきます。まず、の微分の前にに関する関数としてについて考えていきます。これは単純にがそれぞれ微小に動いた時の法線から新たにTangent Spaceを作る関数として考えればよいわけです。[Jakob 2013]では直交座標の作り方としてグラム・シュミットの正規直交化法を採用しており、は
となります。これについて、について微分すると式(1.12)より
というように求めることができます。これによってについても求めることができるというわけです。
以上のことから、頂点のを得ることができればの計算をすることができます。ではこの値とはどうやって得るのかというと3角ジオメトリでは基底方向のベクトルがそのままとなり、頂点の位置をそれぞれで微小にずらした頂点の法線からの値を得ることが可能です。(この辺は[Zeltner et al. 2020]の実装などを見た方が良いと思います)
[Jakob 2013]では付録でのCプログラムが載せられていまして、以下のような形になります。
struct Vertex { Point p; // Vertex position Vector dpdu, dpdv; // Tangent vectors Normal n; // Normal vector Vector dndu, dndv; // Normal derivatives float eta; // Relative index of refraction Matrix2x2 A, B, C; // Matrix blocks of ∇C in the row // associated with the current vertex }; void computeDerivatives(Vertex *v) { /* Compute relevant directions and a few useful projections */ Vector wi = v[-1].p - v[0].p; Vector wo = v[1].p - v[0].p; float ili = 1 / wi.length(); float ilo = 1 / wo.length(); wi *= ili; wo *= ilo; Vector H = wi + v[0].eta * wo; float ilh = 1 / H.length(); H *= ilh; float dot_H_n = dot(v[0].n, H), dot_H_dndu = dot(v[0].dndu, H), dot_H_dndv = dot(v[0].dndv, H), dot_u_n = dot(v[0].dpdu, v[0].n), dot_v_n = dot(v[0].dpdv, v[0].n); /* Local shading tangent frame */ Vector s = v[0].dpdu - dot_u_n * v[0].n; Vector t = v[0].dpdv - dot_v_n * v[0].n; ilo *= v[0].eta * ilh; ili *= ilh; /* Derivatives of C with respect to x_{i-1} */ Vector dH_du = (v[-1].dpdu - wi * dot(wi, v[-1].dpdu)) * ili, dH_dv = (v[-1].dpdv - wi * dot(wi, v[-1].dpdv)) * ili; dH_du -= H * dot(dH_du, H); dH_dv -= H * dot(dH_dv, H); v[0].A = Matrix2x2( dot(dH_du, s), dot(dH_dv, s), dot(dH_du, t), dot(dH_dv, t)); /* Derivatives of C with respect to x_i */ dH_du = -v[0].dpdu * (ili + ilo) + wi * (dot(wi, v[0].dpdu) * ili) + wo * (dot(wo, v[0].dpdu) * ilo); dH_dv = -v[0].dpdv * (ili + ilo) + wi * (dot(wi, v[0].dpdv) * ili) + wo * (dot(wo, v[0].dpdv) * ilo); dH_du -= H * dot(dH_du, H); dH_dv -= H * dot(dH_dv, H); v[0].B = Matrix2x2( dot(dH_du, s) - dot(v[0].dpdu, v[0].dndu) * dot_H_n - dot_u_n * dot_H_dndu, dot(dH_dv, s) - dot(v[0].dpdu, v[0].dndv) * dot_H_n - dot_u_n * dot_H_dndv, dot(dH_du, t) - dot(v[0].dpdv, v[0].dndu) * dot_H_n - dot_v_n * dot_H_dndu, dot(dH_dv, t) - dot(v[0].dpdv, v[0].dndv) * dot_H_n - dot_v_n * dot_H_dndv); /* Derivatives of C with respect to x_{i+1} */ dH_du = (v[1].dpdu - wo * dot(wo, v[1].dpdu)) * ilo; dH_dv = (v[1].dpdv - wo * dot(wo, v[1].dpdv)) * ilo; dH_du -= H * dot(dH_du, H); dH_dv -= H * dot(dH_dv, H); v[0].C = Matrix2x2( dot(dH_du, s), dot(dH_dv, s), dot(dH_du, t), dot(dH_dv, t)); }
まず、頂点のVertexの情報ですが普段使う「位置座標p」、「法線n」、「IOR eta]の他に各微分情報「位置の微分dpdu,dpdv」「法線の微分dndu,dndv」を持っているとしています。これらの情報はジオメトリから計算することができます。なお、ここではの結果を入れるためのA,B,Cという行列を用意しています。
struct Vertex { Point p; // Vertex position Vector dpdu, dpdv; // Tangent vectors Normal n; // Normal vector Vector dndu, dndv; // Normal derivatives float eta; // Relative index of refraction Matrix2x2 A, B, C; // Matrix blocks of ∇C in the row // associated with the current vertex };
の計算はcomputeDerivativesの方で行っています。vertexにはSpecular点をv[0]、その前後の頂点をv[-1],v[1]として持っており、これを入力として受け取っています。
最初の部分では入射、出射、ハーフベクトルなどの計算をしており、wi,wo,Hとして保存しています。また、各長さはili,ilo,ilhとして保存しており、これを用いて各ベクトルは正規化されています。その後、後々に使う計算用に各ベクトルの内積を計算しています。
/* Compute relevant directions and a few useful projections */ Vector wi = v[-1].p - v[0].p; Vector wo = v[1].p - v[0].p; float ili = 1 / wi.length(); float ilo = 1 / wo.length(); wi *= ili; wo *= ilo; Vector H = wi + v[0].eta * wo; float ilh = 1 / H.length(); H *= ilh; float dot_H_n = dot(v[0].n, H), dot_H_dndu = dot(v[0].dndu, H), dot_H_dndv = dot(v[0].dndv, H), dot_u_n = dot(v[0].dpdu, v[0].n), dot_v_n = dot(v[0].dpdv, v[0].n);
次に行うのがのTangent Spaceを計算しています。この時の正規直交基底の生成にはグラム・シュミットの正規直交化法を使用しており、それぞれs,tとして持っています。そのあとのilo,iliの計算はハーフベクトルの計算などに使用するために屈折率で補正を行っています。
/* Local shading tangent frame */ Vector s = v[0].dpdu - dot_u_n * v[0].n; Vector t = v[0].dpdv - dot_v_n * v[0].n; ilo *= v[0].eta * ilh; ili *= ilh;
次に行うことはの微分部分でハーフベクトルの微分を計算しています。ここではまず最初にdH_du, dH_dvで式(1.15)の第一項部分を計算しています。その後、次の式で式(1.15)の第二項目について計算している形です。最終的に得られる結果は行列Aへと式(1.16)の計算結果を入れているわけです。
/* Derivatives of C with respect to x_{i-1} */ Vector dH_du = (v[-1].dpdu - wi * dot(wi, v[-1].dpdu)) * ili, dH_dv = (v[-1].dpdv - wi * dot(wi, v[-1].dpdv)) * ili; dH_du -= H * dot(dH_du, H); dH_dv -= H * dot(dH_dv, H); v[0].A = Matrix2x2( dot(dH_du, s), dot(dH_dv, s), dot(dH_du, t), dot(dH_dv, t));
そのあとに続くものはx_2の微分、つまりSpecular点の微分部分です。これは先ほどと同様ハーフベクトルの微分を計算していますが、最後の行列部分でさらにTangent Spaceの微分式(1.18)を計算している形で行列を求めているというわけです。
/* Derivatives of C with respect to x_i */ dH_du = -v[0].dpdu * (ili + ilo) + wi * (dot(wi, v[0].dpdu) * ili) + wo * (dot(wo, v[0].dpdu) * ilo); dH_dv = -v[0].dpdv * (ili + ilo) + wi * (dot(wi, v[0].dpdv) * ili) + wo * (dot(wo, v[0].dpdv) * ilo); dH_du -= H * dot(dH_du, H); dH_dv -= H * dot(dH_dv, H); v[0].B = Matrix2x2( dot(dH_du, s) - dot(v[0].dpdu, v[0].dndu) * dot_H_n - dot_u_n * dot_H_dndu, dot(dH_dv, s) - dot(v[0].dpdu, v[0].dndv) * dot_H_n - dot_u_n * dot_H_dndv, dot(dH_du, t) - dot(v[0].dpdv, v[0].dndu) * dot_H_n - dot_v_n * dot_H_dndu, dot(dH_dv, t) - dot(v[0].dpdv, v[0].dndv) * dot_H_n - dot_v_n * dot_H_dndv);
最後は終点の計算でこれはと同様に行っています。
/* Derivatives of C with respect to x_{i+1} */ dH_du = (v[1].dpdu - wo * dot(wo, v[1].dpdu)) * ilo; dH_dv = (v[1].dpdv - wo * dot(wo, v[1].dpdv)) * ilo; dH_du -= H * dot(dH_du, H); dH_dv -= H * dot(dH_dv, H); v[0].C = Matrix2x2( dot(dH_du, s), dot(dH_dv, s), dot(dH_du, t), dot(dH_dv, t));
以上での計算をプログラムとして計算することができ、Manifold Walkを実装することができるというわけです。なお、多分ですがTangent Spaceの基底は任意であるため、重心座標系やUV座標系をTangent Spaceの基底として使用することも可能と思われます(いくつか実装を見たところそれぞれの実装で使いやすい物を使用している節があります)。
Manifold Next Event Estimation
Manifold ExplorationはSpecular面が存在するシーンで新たなパスを構築するのに便利でしたが、[Jakob Marschner et al. 2012]ではこの手法はMCMC上で使うことが目的でした。そこで、Manifold ExplorationをPathtracingへと応用したのがManifold Next Event Estimation(MNEE)[Hanika et al. 2015]になります。
元論文では本当は[Jakob Marschner et al. 2012]のSpecular Manifoldをさらに一般化した[Kaplanyan et al. 2014],[Hanika et al. 2015b]のハーフベクトル条件を使った探索で論議しているのですが、多分Ideal reflection,refractionなら多分Manifold Explorationとして考えてもいいと思われます(この辺はちょっとまだ理解不足です)。
以下のような床の上に雫があるようなシーンを考え、その内部に入っていったパスについて考えます。このシーンでは点光源があり、雫の中にが今考えてるパスの点があるとします。
通常のNEEでは点から点を直線的に繋げたようなパスを考えますが、今回の状況では中間に屈折鏡面が存在しているため、図上で点線で表されたパスというのはフェルマーの原理に反し、あり得ない経路となってしまいます。そのため、通常であればこのような状況ではNEEは全く使うことができません。
[Hanika et al. 2015]はこれに対して、NEEで作った不可能なパス(Seed Path)からManifold Explorationによる探索で新たにフェルマーの原理を満たすようなパス(Admissible Path)を生成することで屈折境界面に対応したNEEを行うManifold NEEを提案しました。
MNEEでは、まずとを直接NEEのように繋ぎ、その中間となるSpecular点を得ます。このような直線的なパスはSeed Pathと言い、図上で点線で表される経路です。この時のSeed Pathをと表現します。
SpecularManifoldの方で論議した制約関数にこのSeed Pathを入れるとその値というのはSpecular制限を逸脱しているため、
となります。一方でSpecular上ではどこかしらにはちゃんとフェルマーの原理を満たしたうえでとをつなぐことができるようなパスというのは存在するはずです。このパスはAdmissible Pathと言い、と表記します。この時、は当然以下のようにの解として与えられます。
これは逆を返せば、におけるの解を求めることができれば、をつなぐAdmissible Pathを生成することが可能になると言えるわけです。
ではどうやって解くのかというと、[Hanika et al. 2015]ではManifold Walkを使うことで解となるAdmissible Pathを探索します。今回、Manifold explorationの時とは状況が違い、を固定したうえで、Specular上の点を動かすような形でManifold Walkを行います。
SeedPathのを初期値としてManifold Walkを行い、よりの値が小さくなるようなの移動量を求め、を徐々に動かしていきが0になる位置を求めることでAdmissible Pathを求めます。
イメージとしてはに関するの多次元ニュートン法で解を徐々に求めていくという形です。
Manifold WalkでAdmissible Pathさえ見つかれば、それによる寄与というのは簡単に計算でき(Ideal Specularにおいては)、アルゴリズムとしては以下のようになります。
はスペキュラー上の点であり、これらに対しManifold Walkをすることで新たに解としてを取得します。この経路がAmissible Pathであり、このパスの寄与をそのまま計算するという形でMNEEのアルゴリズムが働きます。なお、Admissible Pathの確率密度はSeedPathの確率密度に等しいですので寄与の計算の際はというSeedPathの確率密度を使用しています(詳しい話は後述)。
残念ながらMNEEの直接的な実装はまだ見つけていないため、MNEEにおけるManifold Walkの具体的な実装法はよくわかっていませんが、[Zeltner et al. 2020]のManifold Walkと思われる実装を見ると恐らくManifold Explorationの時のManifold Walkとそこまで変わらないものと思われます。
以上がMNEEの簡単な説明になります。MNEEは特にシンプルな形状のスペキュラー面では非常にきれいにサンプリングすることができ、以下の図のような中にガラス内部に物体があるようなシーンでは非常に効果的です。また、ガラスによるCausticsもうまく出すことができ、Specular面が強く出るシーンでは非常に強力な手法となっています。
MNEEの元論文[Hanika et al. 2015]ではglossyな物体に対するMNEEなどを考えており、非常に一般的な手法として発表しています。しかも、最近では先ほど述べたようにBlender3.2でShadowCausticsという形で実装されており、その実用性は非常に高いようです。
MNEEの問題点
MNEEはPathTraceでCauticsを出す数少ない有効な手法でした。しかしながら、1つのSeed PathからAdmissible Pathは1つしか作れない、つまりSeed PathとAdmissible Pathは1対1の関係にあるという問題があります。
これはなぜかというと、Admissible Pathの確率密度の計算のためです。今1つのSeed Pathを用意したものとしましょう。本当であればとを通るようなパス、つまりAdmissible Pathは雫の下部分を通っているものもあれば上部分から来たりするパスなどいくつもの経路を考えることができます。
Seed PathはこれらのAdmissible Pathからどれか一つを探索するわけなのですが、あるAdmissibleパスを見つけた際にこののパスの確率密度はSeed Pathに依存し、なおかつ他のAdmissible Pathを取る可能性もあったため
と表されます。については実際にやってみないと計算することが不可能であるため、この積分というのは計算することができません。
MNEEではこの問題に対して、あるSeed Pathに対してただ一つのAddmissible Pathのみを取るという条件を課すことで解決しています。このようにすることではのみ1となり、それ以外は0ということになり、
というようにSeed Pathの確率をそのまま使うことができます。実装上ではManifold Explorationの初期値に乱数を入れなければ固定のAddmissible Pathに収束するのでSeedPathをそのまま与えて探索しています。
これの何が問題かというと、裏を返せばこれは複数個あるAddmissible Pathを1つ以外破棄しているということでもあります。そのため、本来ならいくつも経路があるのにそれらが一切取られないことによってCausticsにノイズが走ってしまう場合があります。これは経路の多さによって生じるため、より複雑な形状を持つような高周波Causticsでは特にMNEEは対応できない場合もあります。これがMNEEの大きな問題点として挙げられています[Zeltner et al. 2020]
Specular Manifold Sampling
MNEEによってパストレでもCausticsを効率的に扱えるようになりましたが、以上のように高周波に限ってはまだまだ難しい問題でありました。しかしながら、比較的最近の2020年の論文[Zeltner et al. 2020] Specular Manifold Sampling for Rendering High-Frequency Caustics and Glintsにて高周波Causticsを効率的に扱えるSpecular Manifold Sampling(SMS)という手法が提案されました。SMSという手法はMNEEの確率密度問題を解決し、1つのSeedPathから複数のAddmissible Pathを扱えるようにした手法となっています。[Zeltner et al. 2020]ではunbiasな方法とbiasな方法の2つを述べています。
unbias SMS
MNEEの問題とは先ほども述べたように、他の解がとられる確率というのが求められないため、以下の確率密度の積分を計算できないということでした。
MNEEではこれを解決するためニュートン法における初期値を固定し、SeedPathとAdmissiblePathに1対1の対応を与えることで制限していました。
一方でSMSは確率密度の逆数をうまい具合に計算することでどうにか以上の積分を数値計算する手法を見つけ、1つのSeedPathに対し複数のAdmissiblePathを生成することが可能しています。
まずはAdmissiblePathを複数サンプルする方法について述べていきます。
ニュートン法によって得られる解というのは与えた初期値に完全に依存しています。例えば、あるという初期値を与えるとその値から出る解というのは必ず1つの解だけになります。もし、他の解が得たいとなると初期値をではなく、別の値で与えて別途ニュートン法を行う他ありません。ですが、必ずしもから以外の解が出るわけではなく、同じ解に収束する可能性もあり、その予測は非常に困難です。
このような性質をよく表したのがニュートンフラクタルです。ニュートンフラクタルはある複素関数の解について、複素平面上の点をニュートン法の初期値として与え、どの解に収束したかで色分けしたものです。ニュートンフラクタルは一般に複雑な形状をしており、初期値と解の複雑な関係性が見て取れると思います。
手法の方に話を戻しまして、MNEEはニュートン法の初期値をDiffuse点と光源点を直接つないだ時の間の点のみとしており、SeedPathから得られる解を1つだけに絞っていたわけです。一方でSMSでは初期値として与えるに対してランダム性を与え、ニュートン法を行うことで他の解についても探索できるようにしています。
のサンプリング方法としてはDiffuse面のBSDFサンプリングの乱数を変えるなど色々考えられますがとにかく2つの一様乱数によって生成されるものと考えます([0,1]じゃないのは後の確率密度の条件を満たすため)。この際、の各点はサンプル空間の点として表現することができます。
こうして作ったをニュートン法に掛けるとそれぞれ別々の解に収束して行きます。解に収束する乱数点をとして、これの集合をと書きます。ニュートンフラクタルのように収束した解で色を付けると以下の図(d)のようになります。それぞれの色の領域がとなります。
こうすることで複数のAdmissiblePathをサンプルできるようにはなりましたが、問題はどうやってその確率密度を求めるかです。は一様乱数であり、上のどの点も等しい確率で選ばれるため、ある解tex:kが選ばれる確率というのはの面積に等しいことになります
はの収束する解がであった時に1になりそれ以外は0となる関数です。この式も積分ではありますが、不偏推定量として以下のように求めることが可能です。
これはつまり、いくつかサンプリングして採用する解に収束した時に1、しなかったら0として総和を取り、平均すればに近づくということです。一応、これで計算は不可能ではなくなったわけです。
しかし、これには大きな問題があります。Pathtraceにおいては確率密度はそのままではなく、逆数で用いられます。ここで期待値の性質として逆数の期待値は元の期待値を逆数にしたものと一致しません。
すなわち
であり、(3.3)を求めたところで実際には使うことができません。ではどうするかというと、[Zeltner et al. 2020]は幾何級数展開を使うことで逆数の期待値を求める手法を取りました。
積分の逆数の不偏推定を求める問題は別のとこでも取り組まれており、[Booth 2007]から研究が始まりました。最近で[Qin et al. 2015]にてフォトンマッピングのファイナルギャザリングで使用する逆数の不偏推定量を求める方法を新たに開発していました。[Zeltner et al. 2020]ではこのアイディアを元に逆数の不偏推定量ついて求める手法を取っています。
このアイディア自体はとても単純でについて、以下のように考えます。
ここでははであり、以外の解が選ばれる確率を意味します。この時である([0,1)の範囲なので最大でも1以下になる)ことから、この値をの幾何級数の総和としてみなすことができます。
これより、逆数の不偏推定量については以下のように求められます。この際、級数部分の各は回目にサンプリングした時にに収束しない確率の推定量という形でみなして表し、それぞれ独立なため各推定量は分けることができます。
[Zeltner et al. 2020]では、この式の推定量について、回目のManifold walkで以外の解に収束した時の部分の値を1、解に収束したときはの部分の値を0として与えることで計算することができると述べています(理由についてはいまいち理解できてません)。こうすることでUnbiasな確率密度の逆数を求めることができるというわけです。
この式の良いところとして、無限級数として表現されていますが、どれか一つでもの解に収束するとそれ以降の項はすべて0になるため、計算を打ち切ることが可能というところです。
以上の話をまとめ、Unbias SMSのアルゴリズムは以下のようになります。
まず初めにランダムにサンプルしたスペキュラー上の点を初期値としてmanifold walkを行い、解を得られたとします(Line1~2)。その次にこの解がとられる確率の逆数を求めていきます。Line3 ~ 9は式3.7の計算であり、別でサンプルしたスペキュラー上の点でmanifold walkを行い新たな解を取得します。この解がであった場合(判定は距離)はこの時点で計算を終了、Line10へと向かいパスの計算を行います。そうでなかった場合はの値に+1して、同様の処理を繰り返していきます。
以上がUnbias SMSのアルゴリズムです。こうすることでMNEEでは扱えなかった高周波Causticsを効率的に扱うことができるようになりました。[Zeltner et al. 2020]でMNEEとUnbias SMSの比較として以下のような図が出されていましたが、 その差はかなり違うことが見て取れ、かなり強力な手法のように見えます。
特に顕著なのが2番目の金属板のCusticsです。MNEEではほとんどCuasticsが見られませんがUnbiasedのものを見るとかなりCausticsが出ていることがわかります。これは特にこのCausticsが高周波(複数のAdmissiblePathが存在する)であるためで、大きな差があることが読み取れます。
以上のようにSMSは高周波Causticsに対してかなり強い手法であるとわかります。
Bias SMS
Unbias SMSは確率密度を求めるために複数の解を求めているのですが、これらは確率密度の計算以外には使われず破棄されています。Manifold walkはそこそこ重い処理であるため、せっかく解を求めたのに全くその寄与を得ないのは中々にもったいない気がします。
せっかく得た解の寄与をどうにか計算しようと思うとMISが思いつきますが、他の解の確率密度を求めなくてはならず、複数の解を探索する必要がまた出てきます。結局のところMISでまとめて寄与を得ることは難しいです。
[Zeltner et al. 2020]はこの問題に対して、Biasではありますがより実用的なSMSとしてBias SMSという手法を提案しています。
M回サンプリングを行い、各解を得たとします。この際、ある解が得られた回数を[tex: nl]と表記し、各解のパスの寄与を平均することで寄与を得ることとします。この時の式は
となります。この式の問題は各解の確率密度を求めるのが難しいということです。ここでを以下のように近似することとします。
単純にこれは確率の定義から考えられることで、これを(3.8)に代入すると
となります。
式(3.10)はであれば同値となりますが、中々にアバウトな近似となっていてある程度バイアスが生じることになります。[Zeltner et al. 2020]によると、このバイアスはエネルギーロスとして現れ、レンダリング結果がある程度暗くなってしまうことが述べられています。しかし、先ほどの図ではunbias、biasの結果がともに載っていますがそこまでバイアスが目立つことはなく、実用上あまり気にならない範囲ではあります。むしろ目立つノイズが消されるため、絵的にはより良い結果であるともいえます。
Bias SMSのアルゴリズムは以下のように書くことができます。
アルゴリズムはとても簡単で、サンプル数Mを指定しその回数分Manifold Walkを行います。この際出てきた解を集め、それぞれの寄与を単純に加算したものを返すというだけです。
Bias SMSによる最初のプールのレンダリング結果を見ると非常にきれいなCausticsが出ていることがわかります。
SMSの説明は大まかにこんな感じになっています。ここでは省略しますが、元論文では他にもより効率的にManifold Walkができるようにパラメーターを変更する手法やGlint(凹凸が非常に激しい鏡面)のレンダリングなど様々な応用についても述ベています。
実装の資料
MNEEの実装は今のところBlenderのリリースノートから飛んだ先のDeveloperのフォーラムに載っているものがありました(まだちょっと全部を見切れていないので確証はありませんが多分MNEEの実装であると思います)。
https://developer.blender.org/rB1fb0247 ちなみに、2021年12月ぐらいにBlenderのShadow Causticsの実装に関するスレッドがあり、ここでMNEEについてやSMSの話題が出ていました(SMSを知ったのはここからでした)。この辺を探ってもいいかもしれないです。
一方でSMSの実装については[Zeltner 2020]でGithub上でMitubaレンダラーの実装として公開されています。一部はMNEEと同様の処理を行うため、これよりMNEEの実装を見てもよいかと思います。
[Zeltner et al. 2020]は他にも先行研究の手法などをまとめた動画を載せており、以上の記事で紹介した手法についてわかりやすくまとめられています。もし、興味があるのでしたら[Zeltner et al. 2020]の論文を一読してから他の論文をたどってみるとわかりやすいと思います。
Specular Manifoldについては元論文の場所の実装リンクが途切れているので今のとこ実装を見つけていません。CausticsをレンダリングできるLuxCoreRendererがMLTで実装されているっぽいことに最近気づいたので、もしかしたら内部組み込まれているのかなと思ったりしています。
おわりに
長い記事でしたがご拝読ありがとうございます。
話を前々から聞いていて参加したかったレイトレ合宿が遂に開催されるということで、何かしら独自性のあるレンダラー作りたいなーとなんかしら探っていたところ先ほどのBlenderのShadow Causticsが発表されました。発表を見てみるとCausticsがちゃんと出てて、これ実装すれば結構良さそうと思い立ちMNEEを実装しようと色々情報を調べていきました。
ですが調べてみるとMNEEは日本語情報が全くないし、理論がめちゃくちゃ難しく、しかも実装が見当たらないという壁に当たり、四苦八苦どうにか先行研究やBlenderのフォーラム漁ってなんとか調べた結果がこの記事で紹介した論文達です。
それらを調べて大体こんな感じなのかなぁっとわかってきた(気がするので)全部まとめたろとやった結果がこの記事ですが、正直言って解説した内容がちゃんとあってるかが不安です。ただこれらに関する資料をまとめた記事とかはないと思うので、もし記事を見て手法に興味を持った人がこれを足掛かりにしてくれれば幸いです。
ここではCausticsに関わる部分だけを抜き出して解説していきましたが、それぞれの論文はCaustics以外にも様々な応用を考えてあります。Manifold ExplorationではGlossyやVolumeの物体にも対応できるように考えられており、MNEEではGlossyの物体にも一般化して理論を組み立てています。興味があれば是非見てみてください(誰かやって解説してください)。
ちなみにSSSのNEEに応用した論文もあるそうで(中身はちゃんと見てないです)、Manifold Explorationを使う手法って意外と多そうな雰囲気がある気がします。
レイトレ合宿までに実装が間に合う気がしない
参考文献、引用論文
[Booth 2007]Unbiased Monte Carlo Estimation of the Reciprocal of an Integral
[Jakob 2013] Light Transport On Path-Space Manifolds
[Hanika et al. 2015a] Manifold Next Event Estimation
[Hanika et al. 2015b] Improved Half Vector Space Light Transport
[Qin et al. 2015]Unbiased photon gathering for light transport simulation
[Zeltner et al. 2020] Specular Manifold Sampling for Rendering High-Frequency Caustics and Glints]