エネルギー保存Microfacet BRDF - ushiostarfish’s diary

この記事はレイトレ合宿6アドベントカレンダー第3週目の記事です。

概要

前回 Microfacet入門（１） - ushiostarfish’s diary では、多重拡散を考慮しないナイーヴなマイクロファセットモデルではエネルギーが完全には保存せず、一部が失われてしまう問題について触れました。そこで今回はその問題への対処として、ややアドホックではありますが、とてもシンプルで軽量な解決方法を提案している、 Csaba Kelemen and László Szirmay-Kalos, "A Microfacet Based Coupled Specular-Matte BRDF Model with Importance Sampling" (以降 [Kelemen 2001] のアイディアに触れてみようと思います。

[Kelemen 2001] のアプローチ

[Kelemen 2001] のアプローチは一言でいうと、単一拡散モデルで失われたエネルギーを、ランバート反射のような均質な拡散反射で補填する、というものです。 [Kelemen 2001] では鏡面反射成分をよくある単一拡散マイクロファセットモデルの以下の形式を用いて記述します。

$\displaystyle { f_{ r,spec }\left( \omega _{ i },\omega _{ o } \right) =\frac { F\left( \lambda ,\omega _{ h }\cdot \omega _{ i } \right) G_{ 2 }\left( \omega _{ h },\omega _{ i },\omega _{ o } \right) D\left( \omega _{ h } \right) }{ 4\left| \omega _{ i },\omega _{ g } \right| \left| \omega _{ o },\omega _{ g } \right| } }$

※論文中では $D\left( \omega _{ h } \right)$ の代わりに $P_{ H }\left( \omega _{ h } \right)$ が使われておりますが、同じです。

Fはフレネル項であり、波長依存性を明示的にλを引数に取ることで主張しつつ、Gはシャドウイング・マスキング関数で、論文中はV-Cavityが用いられています。Dは法線分布関数で、論文中はBeckmann分布が用いられています。
自分も今回はそれにならってV-Cavity, Beckmannを使っていくことにします。（この記事の中ではその詳細には立ち入らず、この部分は後日別記事にするかもしれません）ただし、論文中のコアアイディアとしては、各関数に何を用いるのかはさほど重要なポイントではありません。

当然ながら、これはエネルギーを完全には保存しません。

$\displaystyle { \int _{ \Omega _{ o } }^{ }{ f_{ r,spec }\left( \omega _{ i },\omega _{ o } \right) \cos { \theta _{ o } } d\omega _{ o } } \le 1\forall \omega _{ i } }$

そこで、失われたエネルギーを補うマット成分の項を追加して、

$\displaystyle { f_{ r }\left( \omega _{ i },\omega _{ o } \right) =f_{ r,spec }\left( \omega _{ i },\omega _{ o } \right) +f_{ r,matte }\left( \omega _{ i },\omega _{ o } \right) \\ \int _{ \Omega _{ o } }^{ }{ f_{ r }\left( \omega _{ i },\omega _{ o } \right) \cos { \theta _{ o } } d\omega _{ o } } =1\forall \omega _{ i } }$

これならエネルギーが保存されます。これは強引に解釈するとすれば、多重拡散した結果指向性を失い、マット成分として放出された、というふうに言えなくもないかもしれません。

スペキュラ反射能

実際のところ、エネルギーを補うためには、そもそもどの程度エネルギーが失われているのかを測る必要があります。そこで全体のエネルギー反射率を、 $\theta '$ から入射するときのスペキュラ反射能として、以下のように定義します。

$\displaystyle { a_{ spec }\left( \theta ' \right) =a_{ spec }\left( \vec { L } \right) =\int _{ \Omega _{ o } }^{ }{ f_{ r,spec }\left( \omega _{ i },\omega _{ o } \right) \cos { \theta _{ o } } d\omega _{ o } } }$

スペキュラ反射能は、単位放射照度あたりの、出射の放射輝度を全部かき集めたもの(=放射発散度)、と解釈できます。BRDFの単位は1/srですので、積分した値は無次元数であり、エネルギーの保存率と考えても良さそうです。また、きちんと相反するBRDFを考えているなら、反射能についても出射と入射どちらを考えても同じ関数になります。指向性アルベドとも呼ばれるそうです。

$\displaystyle { a_{ spec }\left( \theta ' \right) =a_{ spec }\left( \vec { L } \right) =a_{ spec }\left( \vec { V } \right) =a_{ spec }\left( \theta \right) }$

マット項

理想的には、多重拡散項がClosed-formに求まるのがベストですが、現状は難しく、Eric Heitz, Johannes Hanika, Eugene d’Eon and Carsten Dachsbacher, "Multiple-Scattering Microfacet BSDFs with the Smith Model" ですら、確率的な評価にならざるを得ず、さらに理論はヘビーでとても難しいです。[Kelemen 2001] では、一昔前の論文ではありますが、もっとシンプルに以下のようにマット成分を記述することにしてます。

$\displaystyle { f_{ r,matte }\left( \omega _{ i },\omega _{ o } \right) =k\left( \lambda \right) \cdot s\cdot r\left( \theta ' \right) \cdot r\left( \theta \right) }$

kは適当な吸収係数であり、rはそれぞれ入射・出射に依存する関数を個別に定義し、sは正規化項として後ほど明らかにします。 kにはλが引数になっているのは、ここに適当な色をつけたりなどが想定されています。

ここでこのマット項の反射能を考えてみると、スペキュラ反射能と同様に、

$\displaystyle { a_{ matte }\left( \theta ' \right) =a_{ matte }\left( \vec { L } \right) =\int _{ \Omega _{ o } }^{ }{ f_{ r,matte }\left( \omega _{ i },\omega _{ o } \right) \cos { \theta _{ o } } d\omega _{ o } } }$

となります。ここで積分を重積分に展開して、 $\theta '$ に依存しない項を外に出すと、

$\displaystyle { a_{ matte }\left( \theta ' \right) =\int _{ \Omega _{ o } }^{ }{ f_{ r,matte }\left( \omega _{ i },\omega _{ o } \right) \cos { \theta _{ o } } d\omega _{ o } } \\ =\int _{ 0 }^{ 2\pi }{ \int _{ 0 }^{ \frac { \pi }{ 2 } }{ k\left( \lambda \right) \cdot s\cdot r\left( \theta ' \right) \cdot r\left( \theta \right) \cos { \theta } \sin { \theta } d\theta } d\phi } \\ =k\left( \lambda \right) \cdot s\cdot 2\pi \cdot r\left( \theta ' \right) \int _{ 0 }^{ \frac { \pi }{ 2 } }{ r\left( \theta \right) \cos { \theta } \sin { \theta } d\theta } }$

この式をよく見てみると、 $a_{ matte }\left( \theta ' \right)$ は、 $r\left( \theta ' \right)$ に比例しており、他の項は定数であることがわかります。ここで $a_{ matte }\left( \theta ' \right)$ は別の観点、スペキュラ反射能のロスがマット反射能だという仮定からみると、

$\displaystyle { a_{ matte }\left( \theta ' \right) =1-a_{ spec }\left( \theta ' \right) }$

もまた同時に言え、 $a_{ matte }\left( \theta ' \right)$ は、 $1-a_{ spec }\left( \theta ' \right)$ に比例しており、さらにこれは入射と出射の立場を入れ替えても同じことが言えるわけなのでもはや、

$\displaystyle { r\left( \theta ' \right) =1-a_{ spec }\left( \theta ' \right) \\ r\left( \theta \right) =1-a_{ spec }\left( \theta \right) }$

としてしまっても差し支えないのではないでしょうか。すると、

$\displaystyle { f_{ r,matte }\left( \omega _{ i },\omega _{ o } \right) =k\left( \lambda \right) \cdot s\cdot r\left( \theta ' \right) \cdot r\left( \theta \right) \\ f_{ r,matte }\left( \omega _{ i },\omega _{ o } \right) =k\left( \lambda \right) \cdot s\cdot \left( 1-a_{ spec }\left( \theta ' \right) \right) \cdot \left( 1-a_{ spec }\left( \theta \right) \right) }$

となりました。ここで正規化項が何になるか考えるために、

$\displaystyle { a_{ matte }\left( \theta ' \right) =\int _{ \Omega _{ o } }^{ }{ f_{ r,matte }\left( \omega _{ i },\omega _{ o } \right) \cos { \theta _{ o } } d\omega _{ o } } }$

から初めて、

$\displaystyle { 1-a_{ spec }\left( \theta ' \right) =\int _{ \Omega _{ o } }^{ }{ f_{ r,matte }\left( \omega _{ i },\omega _{ o } \right) \cos { \theta _{ o } } d\omega _{ o } } \\ 1-a_{ spec }\left( \theta ' \right) =\int _{ \Omega _{ o } }^{ }{ s\cdot \left( 1-a_{ spec }\left( \theta ' \right) \right) \left( 1-a_{ spec }\left( \theta \right) \right) \cos { \theta _{ o } } d\omega _{ o } } \\ \left( 1-a_{ spec }\left( \theta ' \right) \right) =s\cdot \left( 1-a_{ spec }\left( \theta ' \right) \right) \int _{ \Omega _{ o } }^{ }{ \left( 1-a_{ spec }\left( \theta \right) \right) \cos { \theta _{ o } } d\omega _{ o } } \\ 1=s\int _{ \Omega _{ o } }^{ }{ \left( 1-a_{ spec }\left( \theta \right) \right) \cos { \theta _{ o } } d\omega _{ o } } \\ s=\frac { 1 }{ \int _{ \Omega _{ o } }^{ }{ \left( 1-a_{ spec }\left( \theta \right) \right) \cos { \theta _{ o } } d\omega _{ o } } } \\ s=\frac { 1 }{ \int _{ \Omega _{ o } }^{ }{ \cos { \theta _{ o } } d\omega _{ o } } -\int _{ \Omega _{ o } }^{ }{ a_{ spec }\left( \theta \right) \cos { \theta _{ o } } d\omega _{ o } } } \\ s=\frac { 1 }{ \pi -\int _{ \Omega _{ o } }^{ }{ a_{ spec }\left( \theta \right) \cos { \theta _{ o } } d\omega _{ o } } } \\ s=\frac { 1 }{ \pi -\pi \frac { 1 }{ \pi } \int _{ \Omega _{ o } }^{ }{ a_{ spec }\left( \theta \right) \cos { \theta _{ o } } d\omega _{ o } } } }$

ここで、

$\displaystyle { a_{ spec }^{ avg }=\frac { 1 }{ \pi } \int _{ \Omega }^{ }{ a_{ spec }\left( \theta \right) \cos { \theta } d\omega } }$

と平均のアルベドを定義すると、

$\displaystyle { s=\frac { 1 }{ \pi -\pi a_{ spec }^{ avg } } \\ s=\frac { 1 }{ \pi \left( 1-a_{ spec }^{ avg } \right) } }$

とスッキリと記述できてしまいました。これをsの代わりに使うと、

$\displaystyle { f_{ r,matte }\left( \omega _{ i },\omega _{ o } \right) =k\left( \lambda \right) \cdot \frac { \left( 1-a_{ spec }\left( \theta ' \right) \right) \cdot \left( 1-a_{ spec }\left( \theta \right) \right) }{ \pi \left( 1-a_{ spec }^{ avg } \right) } }$

とマット項が定まりました。なんともシンプルな結果であり、かつ相反性を満たし、エネルギーを完全に保存します
ただ残念なことは、 $a_{ spec }^{ avg }$ も $a_{ spec }\left( \theta \right)$ も解析的に求めることができず、事前計算でテーブルに格納しなければならない点です。等方的なら、前者は２次元テーブル（α, θ）、後者は１次元テーブル（α）で済みます。

フレネル項

論文中では、 $a_{ spec }\left( \theta \right)$ を求めるとき、スペキュラ項に誘導体フレネルを含んでいます。この実装では、以下のような見た目になります。

f:id:ushiostarfish:20180624201701p:plain

※Figure 7, [Kelemen 2001]

多層マテリアルのような見た目です。これはある意味当たり前で、誘導体の中に潜り込んだ分のエネルギーを拡散項として放出しているからです。

これはこれで綺麗だと思う反面、物理的な意味に乏しいという欠点はあります。ただこれはかならずしも誘導体フレネルでやらねければならないわけではなく、導体マテリアルのエネルギー補填として使うことを考え、フレネル項を省いたスペキュラ項で $a_{ spec }\left( \theta \right)$ を求めておいて使うことが可能です。こちらも物理的意味には乏しいものの、Eric Heitz, Johannes Hanika, Eugene d’Eon and Carsten Dachsbacher, "Multiple-Scattering Microfacet BSDFs with the Smith Model" をGround Truthとした近似手法が一部で試みがあるようで、Self Shadow, "A Multi-Faceted Exploration (Part 1)" や PataBlog, "Multiple-Scattering BRDF" で確認することができます。以降は導体マテリアルに使うことをメインよりにすすめていこうと思います。

$a_{ spec }^{ avg }$ と $a_{ spec }\left( \theta \right)$ のベイク

[Kelemen 2001] の手法は、アイディアと式はとてもシンプルな反面、テーブルへの格納とサンプリングを余儀なくされるため、実装側はやや面倒です。ここでは計算結果を Cereal を使ってデータをまるっとシリアライズすることにします。また、テーブルを引くとときに、補間処理が必要なので、とりあえずwikipediaからコピペしてきたbicubic補間を使うことにしました。汎用的に使えるようにテンプレートにしつつ、以下のようにしました。

これを使って、まずは $a_{ spec }\left( \theta \right)$ です。

保存は以下のようにします

$a_{ spec }\left( \theta \right)$ は、なんだったかというと、

でした。したがって、V-Cavity, Beckmann を使ってベイクするコードは、

のようにしました。導体と誘導体の２種類をここではベイクする想定です。積分方法は、 $a_{ spec }\left( \theta \right)$ は単純にモンテカルロ法でやりました。合成シンプソンとかでもいいのかな？とも思ったのですが、αが小さいとき、誤差が大きくなる問題もあったのと、効率的な重点サンプリングが存在しているというのが大きいです。ベイクした結果を可視化するために、exrで保存しておくと、デバッグに便利でした。一部 openframeworks のコードが混在してしまっていますが、やっていることは伝わるのではないでしょうか

ベイクするサイズはどのくらいが良いのか？というのを決めなくてはなりませんが、実際 $a_{ spec }^{ avg }$ も $a_{ spec }\left( \theta \right)$ もかなり緩やかでなめらかな形をしているため、それほど大きな解像度は必要なさそうです。imageworksのスライドでは、32x32を使っているようですが、自分は今回は大きめの64x64にしました。

実際どんなかんじになるのかというと、導体のフレネルのケースでは例えば以下のようになります（画像は拡大しています）。もちろん使ったDやG（もし誘導体ならFも）にも影響されるでしょう

f:id:ushiostarfish:20180701164033p:plain

一番右下の暗いところは、0.7くらいの数値になっていました。粗さαが高くて、入射（出射）角度が0に近い部分は多重拡散のエネルギーロスが激しいことがわかります。

次に $a_{ spec }\left( \theta \right)$ ですが、これは、

$\displaystyle { a_{ spec }^{ avg }=\frac { 1 }{ \pi } \int _{ \Omega }^{ }{ a_{ spec }\left( \theta \right) \cos { \theta } d\omega } }$

でした。しかし実際のところ、

$\displaystyle { a_{ spec }^{ avg }=\frac { 1 }{ \pi } \int _{ \Omega }^{ }{ a_{ spec }\left( \theta \right) \cos { \theta } d\omega } \\ =\frac { 1 }{ \pi } \int _{ 0 }^{ 2\pi }{ \int _{ 0 }^{ \frac { \pi }{ 2 } }{ a_{ spec }\left( \theta \right) \cos { \theta } \sin { \theta } d\theta } d\phi } \\ =\frac { 1 }{ \pi } 2\pi \int _{ 0 }^{ \frac { \pi }{ 2 } }{ a_{ spec }\left( \theta \right) \cos { \theta } \sin { \theta } d\theta } \\ =2\int _{ 0 }^{ \frac { \pi }{ 2 } }{ a_{ spec }\left( \theta \right) \cos { \theta } \sin { \theta } d\theta } }$

と一次元の積分に落とすことができるため、こちらを使うのが良いでしょう。こちらは以下のようにコードにしました。

ベイクのほうは、合成シンプソンを使って、

ベイクは、

のようにしました。これは例えば導体の場合、

f:id:ushiostarfish:20180701154436p:plain

のような結果になりました。

[Kelemen 2001] の近似について

[Kelemen 2001] の2.1.1. Simplification of the geometric term ～あたりで、G項の一部を簡略化することで、BRDFやその重点サンプリングの式を軽量化・テーブル化できることを提案しています。しかしながら、その結果は視覚誤差もあり、後により良い重点サンプリングが提案されているため、この部分の有用度は低いと思われます。導出も別段すごく面白いというわけでもなかったので、ここではこれは使わないことにしました。

実装

あとはストレートに、

をベーシックなマイクロファセットのBRDFに足せば良いです。実装は例えば、

glm::dvec3 brdf_diff = kLambda
    * (1.0 - CoupledBRDFConductor::specularAlbedo().sample(alpha, cos_term_wo))
    * (1.0 - CoupledBRDFConductor::specularAlbedo().sample(alpha, cos_term_wi))
    / (glm::pi<double>() * (1.0 - specularAlbedo));

のようになるでしょう。

重点サンプリング

次に拡散項の重点サンプリングについてです。同時に考えると難しいので、やはりここは $f_{ r,spec }$ と $f_{ r,matte }$ を分離して考えます。出射放射輝度のそれぞれの割合は、

$\displaystyle { 1=\int _{ \Omega _{ i } }^{ }{ \left\{ f_{ r,spec }\left( \omega _{ i },\omega _{ o } \right) +f_{ r,matte }\left( \omega _{ i },\omega _{ o } \right) \right\} \cos { \theta _{ i } } d\omega _{ i } } \\ =\int _{ \Omega _{ i } }^{ }{ f_{ r,spec }\left( \omega _{ i },\omega _{ o } \right) \cos { \theta _{ i } } d\omega _{ i } } +\int _{ \Omega _{ i } }^{ }{ f_{ r,matte }\left( \omega _{ i },\omega _{ o } \right) \cos { \theta _{ i } } d\omega _{ i } } \\ =a_{ spec }\left( \theta \right) +a_{ matte }\left( \theta \right) \\ =a_{ spec }\left( \theta \right) +\left( 1-a_{ spec }\left( \theta \right) \right) }$

となるので、このエネルギー比率をそれぞれをサンプルする確率にすれば良さそうです。ただこれはあくまで吸収係数 $k \left( \lambda \right)$ や、導体想定で $a_{ spec }\left( \theta \right)$ のベイク時にフレネル項が含まれていない場合は実際の $f_{ r,spec }$ には追加でフレネル項が入るわけです。なのでここはちょっと場合分けが必要のように思います。

ベイクするときに誘導体フレネルを含んだ場合

この場合、吸収係数 $k \left( \lambda \right)$ は任意の色をつけるために使う用途が想定されます。そして、したがって例えば以下のような比率でサンプリングを組み合わせる方法があります。

$\displaystyle { a_{ spec }\left( \theta \right) :\widetilde { k } \left( 1-a_{ spec }\left( \theta \right) \right) }$

ここで $\widetilde { k }$ は例えば波長それぞれの平均が使えるでしょう。ただそれは最適かというとそうでもなく、平均の代わりにXYZ色空間へ変換して輝度が集約されたY成分を使ったりする実装もありなのではないでしょうか。まあ $\widetilde { k }$ の実装方法は一旦おいておいて、実際の確率は、

$\displaystyle { P_{ spec }=\frac { a_{ spec }\left( \theta \right) }{ a_{ spec }\left( \theta \right) +\widetilde { k } \left( 1-a_{ spec }\left( \theta \right) \right) } \\ P_{ matte }=\frac { \widetilde { k } \left( 1-a_{ spec }\left( \theta \right) \right) }{ a_{ spec }\left( \theta \right) +\widetilde { k } \left( 1-a_{ spec }\left( \theta \right) \right) } =1-P_{ spec } }$

となります。

導体想定の、ベイクするときにフレネル項を含まなかった場合

$f_{ r,spec }$ にはベイク時には考慮していないフレネル項を追加してレンダリングします。しかし後述しますが、吸収係数 $k \left( \lambda \right)$ を多重拡散時の色を近似する項として使用する予定です。なので結局どちらの項にも同一ではないもののフレネル項を使います。なのでここでは簡単のため吸収係数 $k \left( \lambda \right)$ を気にしないことにします。よって確率は、

$\displaystyle { P_{ spec }=a_{ spec }\left( \theta \right) \\ P_{ matte }=1-a_{ spec }\left( \theta \right) }$

を使うことにします。

組み合わせる場合の確率密度

組み合わせる、とはどういうことかというと、 $P_{ spec }$ の確率で $PDF_{ spec }$ (確率密度) のサンプリングが起こり、 $P_{ matte}$ の確率で $PDF_{ matte}$ (確率密度)のサンプリング起こるようにします。したがってその場合の確率密度は、

$\displaystyle { P_{ spec }\times PDF_{ spec }+P_{ matte }\times PDF_{ matte } }$

とすれば良いことがわかります。

$f_{ r,matte }$ の重点サンプリング

$f_{ r,matte }$ は、以下の式でした。

ここで $\omega_i$ をサンプリングすることを考えた場合、 $\omega_i$ に依存しない項を定数項として全部まとめ上げると、

$\displaystyle { f_{ r,matte }\left( \omega _{ i },\omega _{ o } \right) =T\left( 1-a_{ spec }\left( \theta ' \right) \right) =T\left( 1-a_{ spec }\left( \theta \right) \right) }$

となり、LTEとしては、

$\displaystyle { L_{ o }=\int _{ \Omega _{ i } }^{ }{ L_{ i }T\left( 1-a_{ spec }\left( \theta _{ i } \right) \right) \cos { \theta _{ i } } d\omega _{ i } } }$

とこのようになるので、BRDFのサンプリングとしては、

$\displaystyle { \left( 1-a_{ spec }\left( \theta _{ i } \right) \right) \cos { \theta _{ i } } }$

に比例するように $\omega_i$ がサンプリングされるのが良いと思われます。ここでもし完璧な確率密度関数 $p_{ exact }\left( \omega \right)$ が手に入ったとしたら、以下のような形になるでしょう。

$\displaystyle { p_{ exact }\left( \omega \right) =T\left( 1-a_{ spec }\left( \theta \right) \right) \cos { \theta } \\ 1=\int _{ \omega }^{ }{ T\left( 1-a_{ spec }\left( \theta \right) \right) \cos { \theta } d\omega } }$

T は適切に求められた正規化定数とします。これを $\theta$ の積分に持っていくと、

$\displaystyle { 1=\int _{ 0 }^{ 2\pi }{ \int _{ 0 }^{ \frac { \pi }{ 2 } }{ T\left( 1-a_{ spec }\left( \theta \right) \right) \cos { \theta } \sin { \theta } d\theta } d\phi } \\ 1=\int _{ 0 }^{ 2\pi }{ d\phi \int _{ 0 }^{ \frac { \pi }{ 2 } }{ T\left( 1-a_{ spec }\left( \theta \right) \right) \cos { \theta } \sin { \theta } d\theta } } \\ 1=2\pi \int _{ 0 }^{ \frac { \pi }{ 2 } }{ T\left( 1-a_{ spec }\left( \theta \right) \right) \cos { \theta } \sin { \theta } d\theta } \\ 1=\int _{ 0 }^{ \frac { \pi }{ 2 } }{ 2\pi T\left( 1-a_{ spec }\left( \theta \right) \right) \cos { \theta } \sin { \theta } d\theta } }$

なので実際には、

$\displaystyle { \left( 1-a_{ spec }\left( \theta \right) \right) \cos { \theta } \sin { \theta } }$

に比例するように $\theta$ をサンプリングするのが良いことがわかります。

[Kelemen 2001] ではこれを、解析的ではないものの事前計算テーブルを使った逆関数法が紹介されています。しかしどうにもこの論文では $\sin \theta$ を無視しているようなのと、個人的にはサンプリングについてはunbiasedなものを使いたいと考えたのでこの部分はちょっと変えることにしました。

まずは単純に $\theta_{ i }$ の範囲を均等にn分割して、

f:id:ushiostarfish:20180704001859p:plain

そこの代表点を $\left( 1-a_{ spec }\left( \theta \right) \right) \cos { \theta } \sin { \theta }$ で評価します。あとはその値に比例して範囲が選ばれるようにします。次にその範囲内均一に $\theta$ をサンプリングします。あとは $\phi$ もまた均一にサンプリングすれば、半球に欲しかった性質を持つサンプルを生成できます。

しかし残念ながらこれだけではだめで、我々は実際の確率密度関数が必要です。なのでこの場合どんな確率密度になるのか考えてみます。

まずn分割した区間個別は、

$\displaystyle { b-a=\frac { \frac { \pi }{ 2 } }{ n } =\frac { \pi }{ 2n } }$

であり、それぞれの範囲の中では均一にサンプリングすることにしたため、個別の範囲の確率密度関数 $p_{ x_{ i } }$ と表記すると、

$\displaystyle { p_{ x_{ i } }=\frac { 2n }{ \pi } }$

であると言えます。そして実際のところある特定の一箇所が選ばれる確率密度は、その範囲がまず選ばれていなければいけません。なので $\theta$ の確率密度関数 $p\left( \theta \right)$ は、離散的な確率との積になるはずです。

$\displaystyle { p\left( \theta \right) =p_{ x_{ i } }p_{ discrete }\left( x_{ i } \right) \\ p\left( \theta \right) =\frac { 2n }{ \pi } p_{ discrete }\left( x_{ i } \right) }$

例えば n=1、つまり完全に均一にサンプリングする場合、

$\displaystyle { p\left( \theta \right) =\frac { 2 }{ \pi } }$

と確かに辻褄が合っています。また他に簡単にテストするとして、積分して1になるか確かめてみます。

$\displaystyle { X=\int _{ 0 }^{ \frac { \pi }{ 2 } }{ \frac { 2n }{ \pi } p_{ discrete }\left( \theta \right) d\theta } \\ \frac { \pi }{ 2n } X=\int _{ 0 }^{ \frac { \pi }{ 2 } }{ p_{ discrete }\left( \theta \right) d\theta } }$

ここで $p_{ discrete }\left( \theta \right)$ はきれいな短冊形の関数なので、

$\displaystyle { \frac { \pi }{ 2n } X=\sum _{ i=1 }^{ n }{ \frac { b-a }{ n } p_{ discrete }\left( x_{ i } \right) } \\ \frac { \pi }{ 2n } X=\sum _{ i=1 }^{ n }{ \frac { \frac { \pi }{ 2 } -0 }{ n } p_{ discrete }\left( x_{ i } \right) } \\ \frac { \pi }{ 2n } X=\sum _{ i=1 }^{ n }{ \frac { \pi }{ 2n } p_{ discrete }\left( x_{ i } \right) } \\ X=\sum _{ i=1 }^{ n }{ p_{ discrete }\left( x_{ i } \right) } =1 }$

と確かに1になりました。一方 $\phi$ の確率密度関数はどうなるかというと、こちらも均一にサンプリングするため、

$\displaystyle { p\left( \phi \right) =\frac { 1 }{ 2\pi } }$

になります。これを２次元のサンプリングで捉えると、

$\displaystyle { 1=\int _{ 0 }^{ 2\pi }{ \int _{ 0 }^{ \frac { \pi }{ 2 } }{ p\left( \theta \right) p\left( \phi \right) d\theta } d\phi } \\ 1=\int _{ 0 }^{ 2\pi }{ \int _{ 0 }^{ \frac { \pi }{ 2 } }{ \frac { p\left( \theta \right) p\left( \phi \right) }{ \sin { \theta } } \sin { \theta } d\theta } d\phi } \\ 1=\int _{ \omega }^{ }{ \frac { p\left( \theta \right) p\left( \phi \right) }{ \sin { \theta } } d\omega } }$

したがって、立体角での確率密度関数は、

$\displaystyle { p\left( \theta \right) =\frac { 2n }{ \pi } p_{ discrete }\left( \theta \right) \\ p\left( \phi \right) =\frac { 1 }{ 2\pi } \\ p\left( \omega \right) =\frac { p\left( \theta \right) p\left( \phi \right) }{ \sin { \theta } } \\ p\left( \omega \right) =\frac { n }{ { \pi }^{ 2 }\sin { \theta } } p_{ discrete }\left( \theta \right) }$

となります。

実装例

まずは任意の個数の、数値に比例する離散的なサンプルを生成できるクラスを以下のようにしました。

そして実際にサンプルを生成するところは以下のようにしました。テーブルサイズは32x128ですが、これにあまり根拠はありません。もうちょっと小さくても良いかもしれません。

これらを使って、あとはサンプルと確率密度は、以下のようにしました。

これで重点サンプリングの実装は完了です。この構築アプローチは積分を解く必要がなく、テーブルサイズも小さいため、わざわざベイクする必要はないと判断し、レンダラー起動時に計算することにしました。

White Furnace Test

完全にエネルギーを保存するBRDFでは、最初に触れたとおり、

$\displaystyle { \int _{ \Omega _{ o } }^{ }{ f_{ r }\left( \omega _{ i },\omega _{ o } \right) \cos { \theta _{ o } } d\omega _{ o } } =1\quad \forall \omega _{ i } }$

となるべきです。これはiとoを入れ替えても

$\displaystyle { \int _{ \Omega _{ i } }^{ }{ f_{ r }\left( \omega _{ i },\omega _{ o } \right) \cos { \theta _{ i } } d\omega _{ i } } =1\quad \forall \omega _{ o } }$

これは、LTEにおいて、 $L_{ e }=1$ かつ $L_{ i }=1$ のケースであり、

$\displaystyle { L_{ o }=L_{ e }+\int _{ \Omega _{ i } }^{ }{ L_{ i }f_{ r }\left( \omega _{ i },\omega _{ o } \right) \cos { \theta _{ i } } d\omega _{ i } } \\ =\int _{ \Omega _{ i } }^{ }{ f_{ r }\left( \omega _{ i },\omega _{ o } \right) \cos { \theta _{ i } } d\omega _{ i } } }$

これは、"周囲からの放射輝度がすべて1の場合、出射放射輝度は常に1であるべきだ" という解釈ができます。これをテストする手法を "White Furnace Test" と呼びます。これは単にコードベースの単体テストとしても有効ですが、例えばテストしたいマテリアルを適用した球体を、出射放射輝度が全方位で1、反射0の理想拡散光源で囲んで、レンダリングする。という方法があります。（反射0でなくても、直接照明だけに計算を限定しても良いでしょう）球体自身はいろんな出射角度を持つものの、やはりその出射放射輝度は1であるはずで、球体の後ろ側の背景もやっぱり1になるべき、つまり真っ白な画像が出てくるべきだ、というわけです。

実際、例えば通常の単一拡散マイクロファセット導体マテリアルのフレネルがないケースにおいて、α=1としてこのテストを実施すると、

f:id:ushiostarfish:20180701165833p:plain

のようになり、中心付近の数値が0.7くらいになります。これは $a_{ spec }\left( \theta \right)$ のベイクと結果が一致しています。本質的に同じ方程式を解いているからですね。

一方で完全にエネルギー保存するようなケース、例えば、 Eric Heitz, Johannes Hanika, Eugene d’Eon and Carsten Dachsbacher, "Multiple-Scattering Microfacet BSDFs with the Smith Model" の場合を試してみます。これは、tutorial implementation (documentation, code and unit tests)にソースコードがあるので、ちょっとこれを借りることにします。

MicrosurfaceConductorクラスのevalを呼ぶと f * cosθoが、確率的な評価ではありますが、BRDFを評価することができます。重点サンプリングも提供してくれていますが、あとで使いたいフレネルを含ませるやり方がわからなかったため、ここではevalを使いつつ、重点サンプリングは通常の単一拡散のものを使いつつ、uniformなサンプリングを織り交ぜることにしました。

glm::dvec3 sample(PeseudoRandom *random, const glm::dvec3 &wo) const override {
    glm::dvec3 wi;
    double singleScattering = 0.8;
    if (random->uniform() < singleScattering) {
        wi = VCavityBeckmannVisibleNormalSampler::sample(random, alpha, wo, Ng);
    }
    else {
        wi = UniformHemisphereSampler::sample(random, Ng);
    }
    return wi;
}
double pdf(const glm::dvec3 &wo, const glm::dvec3 &sampled_wi) const override {
    double singleScattering = 0.8;
    double pdf_omega =
        0.8 * VCavityBeckmannVisibleNormalSampler::pdf(sampled_wi, alpha, wo, Ng)
        +
        (1.0 - singleScattering) * UniformHemisphereSampler::pdf(sampled_wi, Ng);
    return pdf_omega;
}

これを使ってWhite Furnace Testを実行すると、

f:id:ushiostarfish:20180701171303p:plain

多少の分散こそありますが、見事に成功しています。では多重拡散の代わりに、新しく作った、エネルギー補填項を持つCouple BRDFはどうでしょうか？

f:id:ushiostarfish:20180701171748p:plain

素晴らしい、確かにWhite Furnace Testに成功しています。

"Multiple-Scattering Microfacet BSDFs with the Smith Model" との比較

もうちょっと複雑なシーンではCouple BRDFはどのような見た目になるのかテストしてみましょう。また、同時に"Multiple-Scattering Microfacet BSDFs with the Smith Model"はclosed-formではないものの、多重拡散の"Ground Truth" と言っても差し支えないと思いますので、それとの比較を積極的にやっていくことにします。

シーンはこんなかんじです。モデルは、McGuire Computer Graphics Archive のMori Knobをお借りしました。

f:id:ushiostarfish:20180701172718p:plain

rouphnessについては、 $\alpha = \left( rouphness \right) ^2$ としました。

ここでまずはフレネル項を除いたケースを考えてみます。

まずは通常の単一拡散の場合です。 f:id:ushiostarfish:20180708003229p:plain

確かにラフネス値が高いとエネルギーが失われた結果、不必要に暗くなっているのがわかります。

一方でCouple BRDFはどうなるかというと、 f:id:ushiostarfish:20180708003236p:plain

"Multiple-Scattering Microfacet BSDFs with the Smith Model" の場合はどうなるかというと、 f:id:ushiostarfish:20180708003252p:plain

物理的な意味には乏しいにもかかわらず、かなり見た目が近いです！では今度はフレネル項を追加することを考えます。

"Multiple-Scattering Microfacet BSDFs with the Smith Model" の場合、Supplemental Material: Implementationを読むと、以下のように追加すれば良さそうです。

float MicrosurfaceConductor::evalPhaseFunction(const vec3& wi, const vec3& wo) const
{
    // half vector 
    const vec3 wh = normalize(wi+wo);
    if(wh.z < 0.0f)
        return 0.0f;
    
#if USE_FRESNEL
    // Fresnel
    const float F = fresnel_unpolarized(n, k, glm::dot(wi, wh));

    // value
    const float value = 0.25f * F * m_microsurfaceslope->D_wi(wi, wh) / dot(wi, wh);
#else 
    // value
    const float value = 0.25f * m_microsurfaceslope->D_wi(wi, wh) / dot(wi, wh);
#endif
    return value;
}

ここでフレネル項は、Ole Gulbrandsen Framestore, "Artist Friendly Metallic Fresnel" でも使われている、複素屈折率 $\eta = n+ik$ をパラメータにもつフレネル反射の近似式（元はPBRT v2で、v3では見当たらない・・・？）である、

$\displaystyle { r_{ \bot }\left( n,k,\theta \right) =\frac { n^{ 2 }+k^{ 2 }-2\eta \cos \theta +\cos ^{ 2 }{ \theta } }{ n^{ 2 }+k^{ 2 }+2\eta \cos \theta +\cos ^{ 2 }{ \theta } } \\ r_{ \parallel }\left( n,k,\theta \right) =\frac { \left( n^{ 2 }+k^{ 2 } \right) \cos ^{ 2 }{ \theta } -2n\cos \theta +1 }{ \left( n^{ 2 }+k^{ 2 } \right) \cos ^{ 2 }{ \theta } +2n\cos \theta +1 } \\ F\left( n,k,\theta \right) =\frac { r_{ \bot }\left( n,k,\theta \right) +r_{ \parallel }\left( n,k,\theta \right) }{ 2 } }$

を使うことにしました。今回はナイーヴなRGBレンダラーを想定して、代表的な波長を

color	wavelength
R	650nm
G	550nm
B	450nm

としました。あまり強い根拠はなく、Complex Fresnel shader - OSL Shaders - Chaos Group Help のマネをしただけです。

例えば金の場合は、

color	n	k
R	0.15557	3.6024
G	0.42415	2.4721
B	1.3821	1.9155

※ Refractive index of Au (Gold) - Johnson より

を使って、これを元に "Multiple-Scattering Microfacet BSDFs with the Smith Model" でレンダリングすると以下のようになりました。

f:id:ushiostarfish:20180708003448p:plain

綺麗です。Coupled でも試してみます。しかしこの場合は一つ問題があり、

の

$k\left( \lambda \right)$ をどうするのか？という点です。普通に考えると多重拡散項もフレネル項の影響を受けるはずです。しかし一旦忘れて $k\left( \lambda \right)$ を1でレンダリングしてみます。すると、

f:id:ushiostarfish:20180708003509p:plain

うーん、近いことは近いです。しかしαが1に近くなるにつれて白くなりすぎています。もうちょっとなんとかならないでしょうか？

多重拡散項の近似

ここの多重拡散項をより近づけるための手法として、今回は siggraph 2017のimageworksの資料および、A Multi-Faceted Exploration (Part 2) - Self Shadow を参考にしてみることにします。

まずある一回の衝突では、全体のエネルギーを1としたとき、平均的に見て $a_{ spec }^{ avg }$ が単一拡散、 $1-a_{ spec }^{ avg }$ が多重拡散として衝突を続けると考えられます。ただ多重拡散項はまた次の衝突で、 $a_{ spec }^{ avg }$ と $1-a_{ spec }^{ avg }$ の割合で分断されます。このプロセスが再帰的に繰り返されると考えられます。そして衝突した分だけフレネル項の影響を受けるでしょう。ただフレネル項も多重拡散分は角度が不明であるため、cos重みづけ平均をここでのフレネル項とします。

$\displaystyle { F_{ avg }\left( \eta ,k \right) =\int _{ \Omega }^{ }{ \frac { 1 }{ \pi } F\left( \eta ,k,\cos { \theta } \right) \cos { \theta } d\omega } }$

実際これは、

$\displaystyle { F_{ avg }\left( \eta ,k \right) =\int _{ \Omega }^{ }{ \frac { 1 }{ \pi } F\left( \eta ,k,\cos { \theta } \right) \cos { \theta } d\omega } \\ =\int _{ 0 }^{ 2\pi }{ \int _{ 0 }^{ \frac { \pi }{ 2 } }{ \frac { 1 }{ \pi } F\left( \eta ,k,\cos { \theta } \right) \cos { \theta } \sin { \theta } d\theta } d\phi } \\ =2\pi \int _{ 0 }^{ \frac { \pi }{ 2 } }{ \frac { 1 }{ \pi } F\left( \eta ,k,\cos { \theta } \right) \cos { \theta } \sin { \theta } d\theta } \\ =\int _{ 0 }^{ \frac { \pi }{ 2 } }{ 2F\left( \eta ,k,\cos { \theta } \right) \cos { \theta } \sin { \theta } d\theta } }$

と単純化でき、さらに

$\displaystyle { \mu =\cos { \theta } \\ \frac { d\mu }{ d\theta } =-\sin { \theta } \\ -d\mu =\sin { \theta } d\theta }$

と変数変換すると、

$\displaystyle { =-\int _{ \cos { 0 } }^{ \cos { \frac { \pi }{ 2 } } }{ 2F\left( \eta ,k,d \right) \mu d\mu } \\ =-\int _{ 1 }^{ 0 }{ 2F\left( \eta ,k,d \right) \mu d\mu } \\ =2\int _{ 0 }^{ 1 }{ F\left( \eta ,k,d \right) \mu d\mu } }$

とさらに簡単になり、なめらかでゆるやかな関数のため、合成シンプソンなどで少ない評価点でも近似できます。

以上をふまえ、 $E=a_{ spec }^{ avg }$ 、フレネルを $F$ として図にしてみます。

f:id:ushiostarfish:20180701180926p:plain

これを式にすると、

$\displaystyle { FE+F\left( 1-E \right) +{ \left\{ F\left( 1-E \right) \right\} }^{ 2 }+{ \left\{ F\left( 1-E \right) \right\} }^{ 3 }+...\\ =FE\sum _{ k=0 }^{ \infty }{ { \left\{ F\left( 1-E \right) \right\} }^{ k } } }$

となり、これは等比数列の和として解くことができるので、

$\displaystyle { a=FE\\ r=F\left( 1-E \right) \\ a+ar+a{ r }^{ 2 }+...a{ r }^{ n-1 }={ \frac { a\left( 1-r^{ n } \right) }{ 1-r } }\\ =\frac { FE\left( 1-{ { \left\{ F\left( 1-E \right) \right\} }^{ \infty } } \right) }{ 1-F\left( 1-E \right) } =\frac { FE }{ 1-F\left( 1-E \right) } }$

ただ、実際今考えているのは多重拡散項であり、最初の $FE$ は不要です。したがって、

$\displaystyle { \frac { FE }{ 1-F\left( 1-E \right) } -FE=\frac { FE }{ 1-F\left( 1-E \right) } -\frac { 1-F\left( 1-E \right) }{ 1-F\left( 1-E \right) } FE\\ =\frac { FE }{ 1-F\left( 1-E \right) } -\frac { FE-EFF\left( 1-E \right) }{ 1-F\left( 1-E \right) } =\frac { E{ F }^{ 2 }\left( 1-E \right) }{ 1-F\left( 1-E \right) } }$

ところが実際知りたかったのは、余剰分のエネルギーがフレネルの影響がどの程度あるかということであります。それは、 $\frac { E_{ fresnel } }{ E_{ no-fresnel } }$ と考えられます。なのでつまりは、

$\displaystyle { \frac { E_{ fresnel } }{ E_{ no-fresnel } } =\frac { \frac { E{ F }^{ 2 }\left( 1-E \right) }{ 1-F\left( 1-E \right) } }{ \left( 1-E \right) } =\frac { E{ F }^{ 2 } }{ 1-F\left( 1-E \right) } }$

が、欲しかった多重拡散用のフレネル項であると言えます。これをコードにすると、

inline double fresnel_avg(double n, double k) {
    return 2.0 * rt::composite_simpson<double>([&](double cosTheta) {
        return fresnel_unpolarized(n, k, cosTheta) * cosTheta;
    }, 10, 0.0, 1.0);
}

glm::dvec3 F(
    fresnel_avg(eta.r, k.r),
    fresnel_avg(eta.g, k.g),
    fresnel_avg(eta.b, k.b)
);
glm::dvec3 E = glm::dvec3(specularAlbedo);
glm::dvec3 ONE = glm::dvec3(1.0);
kLambda = E * F * F / (ONE - F * (ONE - E));

これをレンダリングすると、

f:id:ushiostarfish:20180708003629p:plain

お、わりと近いのではないでしょうか？しかしながら、ほんの僅かにα=1のとき、色が強く出すぎているような気がします。感覚的ですが、自分はいささか唐突にも感じます。どうにもF項を過剰評価しているのではないかという気がします。いっそのこと、 $F_{ avg }$ をかけるだけ、にしてみてはどうか？というアイディアがあったのでそれも試してみることにしました。

glm::dvec3 F(
    fresnel_avg(eta.r, k.r),
    fresnel_avg(eta.g, k.g),
    fresnel_avg(eta.b, k.b)
);
kLambda = F;

これらを比較してみると、 f:id:ushiostarfish:20180708003709p:plain

$F_{ avg }$ が近すぎる・・・！

別な材質、銅を試してみます。

color	n	k
R	0.23780	3.6264
G	1.0066	2.5823
B	1.2404	2.3929

※ Refractive index of Cu (Copper) - Johnson より

f:id:ushiostarfish:20180708003722p:plain

やはり $F_{ avg }$ が近すぎる・・・！自分とんでもないポカをやっているのだろうか・・・？気づいた方いらっしゃいましたらご一報いただけると嬉しいです。。

自分のレンダラーには、高ラフネス値時の色の変化がやはり唐突に感じるのと、手元のテストの結果を優先し、 $F_{ avg }$ を使うことにしました。

まとめ

[Kelemen 2001] のアプローチは大変シンプルなアイディアで、理論も軽量でありながら、相反性、エネルギー保存の問題を解消でき、さらにはエネルギー保存に一歩届かないマテリアルのエネルギー補填として、とても汎用的に使えるテクニックであるため、古い論文ながらも有用性は高いのではないでしょうか。リアルタイムレンダリングでも使えるのではないかと思います。一方でテーブルへの格納するコードは若干面倒臭く、コードも増えてしまう欠点もありますが、このあたりは解析的な関数でのフィッティングなどがあれば、いくらか良いのかもしれません。また、ここではBRDFのみにフォーカスしておりBSDFへの拡張は、おざなりでした。こちらもチャレンジの余地があります。Revisiting Physically Based Shading at Imageworks, (Christopher Kulla and Alejandro Conty) SIGGRAPH 2017 Course: Physically Based Shading in Theory and Practice - Self Shadow の[supplemental] では、Alejandro Conty Estevez, Christopher Kulla, "Production Friendly Microfacet Sheen BRDF" にて velvet-like material の中でもこの手法が使われており、こちらもチャレンジしてみるのも良さそうです。