Microfacetモデルの重点サンプリング

こんにちはtatsyです。

今回の記事は レイトレ合宿5 アドベントカレンダー の第8回の記事になります。

今回の記事では、Microfacetモデルに用いられる法線分布関数の代表例であるBeckmann分布とGGX分布について、重点サンプリングに用いる公式を導出し、それをsmallpt ベースのレンダラーに組み込むところまでをご紹介しようと思っています。今回の記事の内容は主に以下の2つの論文の内容をまとめたものになりますので、ご興味のある方は原著の方もご覧ください。

  • Walter et al. 2007, “Microfacet Models for Refraction through Rough Surfaces”, EGSR 2007.
  • Heitz et al. 2014, “Importance Sampling Microfacet-Based BSDFs with the Distribution of Visible Normals”, EGSR 2014.

Microfacetモデルの基本的な紹介については、昨年のレイトレ合宿アドベントカレンダーでPheemaさんが書かれた 「 脱・完全鏡面反射~GGXについて調べてみた~ 」にも紹介がされていますので、こちらも合わせてお読みいただければと思っております。
 

Microfacetモデルとは?


Microfacetモデルは物体表面にある細かな凹凸によって光の反射方向の分布が定まるというBSDF (BRDF + BTDF) の計算モデルの一つです。Microfacetモデルでは物体の表面は完全な鏡面反射あるいは鏡面透過であると仮定をして、いわゆる光沢面の効果は物体表面に現れる細かな凹凸に起因するものであるという風に考えます。

このとき物体表面のある位置に存在する細かな凹凸の法線をMicrofacet法線と呼び\omega_mで表します。また、その法線の分布を表す関数を法線分布関数 D(\omega_m)で表します。詳細は割愛しますが、この法線分布関数を用いると物体表面での反射および透過は次の式で表せます [Walter et al.2007]。

反射の場合 (BRDF)

(1)    \begin{equation*} f_r(\omega_i, \omega_o, n) = \frac{F(\omega_i, \omega_m) G(\omega_i, \omega_o, \omega_m) D(\omega_m)}{4 | \omega_o \cdot n | | \omega_i \cdot n |} \end{equation*}

透過の場合 (BTDF)

(2)    \begin{equation*} f_t(\omega_i, \omega_o, n) = \frac{|\omega_i \cdot \omega_m||\omega_o \cdot \omega_m|}{|\omega_i \cdot n||\omega_o \cdot n|} \frac{\eta_o^2 (1 - F(\omega_i, \omega_m)) G(\omega_i, \omega_o, \omega_m) D(\omega_m)}{(\eta_i (\omega_o \cdot \omega_m) + \eta_o (\omega_i \cdot \omega_m))^2} \end{equation*}

 

上記の式においてF(\omega_i, \omega_m)はフレネル反射の影響を表す項、G(\omega_i, \omega_m, \omega_m)はGeometry attenuation function (GAF)と呼ばれる、物体表面の凹凸による遮蔽の影響を考慮した項になります。GAFには多くの場合、Smith-Bourliear GAFあるいはSmithのMasking-shadowing関数と呼ばれるものを使うのが一般的なようです。こちらのついては後ほどもう少し掘り下げて説明いたします。

Microfacetモデルを使う方向サンプリングの場合、基本的にサンプリングする対象はMicrofacet法線です。従って実際に積分をするときには\omega_oから\omega_mへの変数変換に関する項が入ってきます。実はここの導出はあまり直観的ではないのですが[Walter et al. 2007]にある説明によると以下のようにヤコビアンが求まります([Walter et al. 2007]の図6および図7を参照)。

反射の場合

(3)    \begin{equation*} \left\| \frac{\partial \omega_m}{\partial \omega_o} \right\| = \frac{1}{4 | \omega_o \cdot \omega_m |} \end{equation*}

透過の場合

(4)    \begin{equation*} \left\| \frac{\partial \omega_m}{\partial \omega_o} \right\| = \frac{\eta_o^2 | \omega_o \cdot \omega_m |}{(\eta_i (\omega_i \cdot \omega_m) + \eta_o (\omega_o \cdot \omega_m))^2} \end{equation*}

 

これらのヤコビアンは以下で寄与に対する重みの計算に使いますので、頭の片隅に置いておいていただけると助かります。

 

SmithのMasking-shadowing関数


SmithのMasking-shadowing関数はより簡素な二つの関数の積として、以下のように近似できることが知られています。

(5)    \begin{equation*} G(\omega_i, \omega_o, \omega_m) = G_1(\omega_i, \omega_m) G_1(\omega_o, \omega_m) \end{equation*}

この式に現れるG_1は、とある方向\omegaから見える法線の影響を考慮したものになっていて、
式で表現すると以下のようになります。

(6)    \begin{equation*} G_1(\omega, \omega_m) = \frac{\chi^{+}(\omega \cdot \omega_m) (\omega \cdot n)}{\int_{\Omega} \omega \cdot \omega_m' D(\omega_m') d \omega_m } \end{equation*}

この式を実際に計算していくと、(解析的に解けるかどうかはさておき) 以下の形に変形できます。

(7)    \begin{equation*} G_1(\omega, \omega_m) = \frac{1}{1 + \Lambda(\omega)} \end{equation*}

Beckmann分布やGGX分布を扱う場合には、式中の\Lambda(\omega)が式(6)を変形していくことで解析的に得られます [Heitz et al. 2014a]。しかしながら、この導出にはスロープ表現を用いる必要がありますので、こちらについては後ほどご説明することにします。
 

基本的な方向サンプリングの流れ


記事の冒頭で申し上げた通り、今回の記事では[Walter et al. 2007]と[Heitz et al. 2014b]の方法の2つの方法をご紹介いたします。いずれの方法の場合にも方向サンプリングのやり方はとても似ていて、それぞれ以下のような流れになります。

Walter et al. 2007の場合

  1. 法線分布D(\omega_m) (\omega_i \cdot \omega_o)からMicrofacet法線\omega_mをサンプリングする
  2. 得られたMicrofacet法線を基に完全鏡面での反射および透過を計算する
  3. 方向とは別に最終的な寄与に乗ずる重みを計算する
     

Heitz et al. 2014の場合

  1. 可視法線分布D_{\omega_i}(\omega_m)からMicrofacet法線\omega_mをサンプリングする
  2. 得られたMicrofacet法線を基に完全鏡面での反射および透過を計算する
  3. 方向とは別に最終的な寄与に乗ずる重みを計算する
     

法線分布関数とサンプリング公式の導出


Microfacetモデルにより定まる光沢面 (Glossy surface) をBSDFにより表現する場合、
レンダリングの分野ではBeckmann分布とGGX分布 (Trowbridge-Reitz分布) という分布がよく使われます。
もし光沢面が方向に依らない等方的な法線分布を持っている場合にはBeckmann分布とGGX分布に対する
方向サンプリングの公式は解析的に求まります。これが、この二つの分布がよく用いられる理由です。
一方で、非等方的な法線分布を持っている場合には完全に解析的にはサンプリング公式が決まらないため、
数値的な処理を加えて方向をサンプリングすることになります。

Beckmann分布とGGX分布はそれぞれ等方の場合、非等方の場合で以下のように定義されます。
二つの分布は良く見るとBeckmann分布もGGX分布も似た形をしていることに気づきます。

Beckmann分布 (等方)

(8)    \begin{equation*} D(\omega_m) = \frac{\chi^+(\omega_m \cdot \omega_g)}{\pi \alpha^2 \cos^4 \theta} \exp \left( -\frac{\tan^2 \theta_n}{\alpha^2} \right) \end{equation*}

 

Beckmann分布 (非等方)

(9)    \begin{equation*} D(\omega_m) = \frac{\chi^+(\omega_m \cdot \omega_g)}{\pi \alpha_x \alpha_y \cos^4 \theta} \exp \left( -\tan^2 \theta_n (\frac{\cos^2 \phi_n}{\alpha_x^2} + \frac{\sin^2 \phi_n}{\alpha_y^2}) \right) \end{equation*}

 

GGX分布 (等方)

(10)    \begin{equation*} D(\omega_m) = \frac{\chi^+ (\omega_m \cdot \omega_g)}{\pi \alpha^2 \cos^4 \theta_n \left( 1 + \frac{\tan^2 \theta_n}{2} \right)^2} \end{equation*}

 

GGX分布 (非等方)

(11)    \begin{equation*} D(\omega_m) = \frac{\chi^+ (\omega_m \cdot \omega_g)}{\pi \alpha_x \alpha_y \cos^4 \theta_n \left( 1 + \tan^2 \theta_n (\frac{\cos^2 \phi_n}{\alpha_x^2} + \frac{\sin^2 \phi_n}{\alpha_y^2})  \right)^2} \end{equation*}

 

これらの分布の定義はご覧の通りかなり複雑です。実際、これら分布をそのままの形でサンプルすることは難しく、[Walter et al. 2007]の論文では、この分布に光の入射方向 \omega_i とMicrofacet法線 \omega_mの内積を付け加えることで、解析的にサンプリング公式を導出しています。等方的な分布については[Walter et al. 2007]でもサンプリング公式が示されていますし、導出方法もそれほど複雑ではありませんので、この記事では非等方の分布の方でサンプリングの公式を導出してみます。

Beckmann分布の場合

確率分布に従う変数をサンプリング公式は、確率分布の累積密度分布を求めた後、逆関数法と呼ばれる方法で導出するのが一般的です。まずはBeckmann分布について \omega_m = (\theta_m, \phi_m) に関する累積密度分布を求めてみます。

今回の分布は二つの角度を取る分布になっていますので、欲しい変数にだけ注目するために、それ以外の変数を定義域全体で積分して削除します。最初は\phi_mのサンプリング公式を導出するために \theta_m を削除します。Beckmann分布の式(9)で\theta_mに関係ない変数を適当な文字で置くと、Beckmann分布は以下のように書き直せます。

(12)    \begin{equation*} D(\omega_m) = \frac{1}{A \cos^4 \theta_n} \exp \left( -B \tan \theta_n \right) \end{equation*}

[Walter et al. 2007]に従い、 \omega_i \cdot \omega_mを掛け算し、これを半球上で積分します。実際に今積分する変数は \theta_m だけなので動く範囲は [0, \pi / 2] になります。以下の積分では立体角\omega_mを角度(\theta_m, \phi_m)に変数変換した影響で\sin \theta_mが入ってくることにご注意ください。

(13)    \begin{eqnarray*} \int_0^{\frac{\pi}{2}} \frac{1}{A \cos^4 \theta_n} \exp \left( -B \tan^2 \theta_n \right) \cos \theta_n \sin \theta_n d \theta_n &=& \int_0^{\frac{\pi}{2}} \frac{\sin \theta_n}{A \cos^3 \theta_n} \exp \left(-B \tan^2 \theta_n \right) d \theta_n \\ &=& \int_0^{\frac{\pi}{2}} \left( -1 / (2 A B) \right) \left( \exp \left( -B \tan^2 \theta_n \right) \right)' d \theta_n \\ &=& \left. -\frac{1}{2AB} \exp \left(-B \tan^2 \theta_n \right) \right\vert_{\theta_n = 0}^{\theta_n=\frac{\pi}{2}} \\ &=& \frac{1}{2AB} \\ &=& \frac{1}{2 \pi \alpha_x \alpha_y} \left( \frac{\cos^2 \phi_n}{\alpha_x^2} + \frac{\sin^2 \phi_n}{\alpha_y^2} \right)^{-1} \end{eqnarray*}

これで \phi_n に関する周辺分布が求まりました。ここからさらに\phi_mに関する累積密度分布を求めるため、式を整理して簡略化します。

(14)    \begin{equation*} D(\phi_n) = \frac{1}{\pi} \frac{\alpha_y}{\alpha_x} \frac{1}{1 + ((\frac{\alpha_y}{\alpha_x})^2 - 1) \cos^2 \phi_n} := \frac{P}{1 + Q cos^2 \phi_n} \end{equation*}

この表現を用いると累積密度分布は以下のように計算できます。かなり複雑な積分ですのでWalframAlphaを使って不定積分を求めていますが、もし手計算したいのであればt = \tan \phi_m / \sqrt{Q + 1}と置いて変数変換すると正しく計算できます。

(15)    \begin{eqnarray*} C(\phi_n) &=& \int_0^{\phi_n} \frac{P}{1 + Q cos^2 \phi_n'} d \phi_n' \\ &=& \frac{P}{\sqrt{Q + 1}} \arctan \left( \frac{\tan \phi_n}{\sqrt{Q+1}} \right) \\ &=& \frac{1}{2 \pi} \arctan \left( \frac{\alpha_x}{\alpha_y} \tan \phi_n \right) \end{eqnarray*}

この累積密度分布を使うと逆関数法から以下のサンプリングの公式が得られます。なお\mathcal{U}[0, 1][0, 1]の範囲における一様乱数を表します。

(16)    \begin{equation*} \phi_m = \arctan \left( \frac{\alpha_y}{\alpha_x} \tan 2 \pi u_1 \right) \qquad u_1 \sim \mathcal{U}[0, 1] \end{equation*}

ここで \phi_m を定数だと思って D(\theta_m) の累積密度分布を求めます。式(13)の最後に出てくる数はD(\theta_n)の規格化定数とみなせますので、以下のように累積密度分布が求まります。

(17)    \begin{eqnarray*} C(\theta_m) &=& \frac{ \int_{0}^{\theta_m} D(\omega_m) d \theta_m' }{ \int_{0}^{\frac{\pi}{2}} D(\omega_m) d \theta_m } \\ &=& 1 - \exp \left(-\tan^2 \theta_m \left( \frac{\cos^2 \phi_m}{\alpha_x^2} + \frac{\sin^2 \phi_m}{\alpha_y^2} \right) \right) \end{eqnarray*}

この累積密度分布から、逆関数法により\theta_mのサンプリング公式が得られます。

(18)    \begin{equation*} \theta_m = \arctan \left(- \log(1 - u_2) \left( \frac{\cos^2 \phi_m}{\alpha_x^2} + \frac{\sin^2 \phi_m}{\alpha_y^2} \right)^{-1} \right) \qquad u_2 \sim \mathcal{U}[0, 1] \end{equation*}

以上、2つのサンプリング公式からMicrofacet法線が求まりますので、ここから方向を決めることが可能です。
 

GGX分布の場合

上記のBeckmann分布に対する計算の流れは、GGX分布の場合にも概ね同じです。まずはGGX分布を扱いやすいように以下のように書き直します。

(19)    \begin{equation*} D(\omega_m) = 1 - \frac{1}{A \cos^4 \theta_m (1 + B \tan^2 \theta_m)^2} \end{equation*}

この表現を使い、\theta_mに関して定義域全体を積分します。これにより\phi_mに関する周辺分布が以下のように求まります。

(20)    \begin{eqnarray*} \int_0^{\frac{\pi}{2}} \frac{1}{A \cos^4 \theta_n (1 + B \tan^2 \theta_n)} \cos \theta_n \sin \theta_n d \theta_n &=& \int_0^{\frac{\pi}{2}} \frac{\sin \theta_n}{A \cos^3 \theta_n (1 + B \tan^2 \theta_n)^2} d \theta_n \\ &=& \int_0^{\frac{\pi}{2}} (-1 / (2AB)) \left(\frac{1}{1 + B \tan^2 \theta_n} \right)' d \theta_n \\ &=& \frac{1}{AB} \\ &=& \frac{1}{2 \pi} \arctan \left( \frac{\alpha_x}{\alpha_y} \tan \phi_n \right) \end{eqnarray*}

この結果は式(16)と同じですので、Beckmann分布同様に \phi_m のサンプリング公式は以下のようになります。

(21)    \begin{equation*} \phi_m = \arctan \left( \frac{\alpha_y}{\alpha_x} \tan 2 \pi u_1 \right) \qquad u_1 \sim \mathcal{U}[0, 1] \end{equation*}

続いて D(\theta_m) の累積密度関数は、

(22)    \begin{equation*} C(\theta_m) = \frac{1}{1 + \tan^2 \theta_m \left( \frac{\cos^2 \phi_m}{\alpha_x^2} + \frac{\sin^2 \phi_m}{\alpha_y^2} \right)} \end{equation*}

となるので \theta_m のサンプリング公式は逆関数法により以下のように求まります。

(23)    \begin{equation*} \theta_n = \arctan \sqrt{ \left( \frac{u_2}{1 - u_2} \left( \frac{\cos^2 \phi_n}{\alpha_x^2} + \frac{\sin^2 \phi_n}{\alpha_y^2} \right) \right)^{-1} } \qquad u_2 \sim \mathcal{U}[0, 1] \end{equation*}

以上から、GGX分布についてもMicrofacet法線をサンプリングする公式が得られました。
 

重みの計算

Beckmann分布、GGX分布のいずれを使う場合も得られる寄与に式(1), (2)の表現に応じた重みを乗じる必要があります。上記のサンプリングの計算、および出射方向を決める計算において法線分布関数の影響であるD(\omega_m)およびフレネル反射の影響であるF(\omega_i, \omega_m) は考慮されています。また積分をする対象が\omega_oから\omega_hに変更されていることを考慮し、式(3), (4)のヤコビアンの影響を考慮します。サンプリング時に cos \theta_m = (\omega_m \cdot n)を掛けてサンプリングした影響を打ち消すよう、 1 / |\omega_m \cdot n |を乗じます。最後に実装の都合上、BSDFはコサイン項をあらかじめ考慮に入れたサンプルをするので、あらかじめ |\omega_o \cdot n| を掛け算しておきます。

これらを全て考慮すると、サンプルに対する重みは反射、透過の場合ともに以下のようになります。

(24)    \begin{equation*} \text{weight}(\omega_o) = \frac{|\omega_i \cdot \omega_m| G(\omega_i, \omega_o, \omega_m)}{|\omega_i \cdot n ||\omega_m \cdot n|} \end{equation*}

あとはMasking-shadowing関数G(\omega_i, \omega_o, \omega_m)が計算できれば重みが決まるのですが、これを解析的に求めるにはスロープ表現の説明が必要になりますので、この詳細については後ほどご説明いたします。
 

スロープ表現を用いた可視法線分布のサンプリング


スロープ表現を用いた可視法線のサンプリングは2014年にHeitzらによって発表されました[Heitz et al. 2014b]。このサンプリング手法は、最新版のPBRTにも採用されています。先ほどご紹介したWalterらの研究に基づくサンプリング手法は光がどの方向から入射してくるかに関係なく同じようにMicrofacet法線をサンプリングしていました。ですがHeitzらの手法はスロープ表現を使い、光の入射方向から見える法線、すなわち可視法線だけをサンプリングします。

この可視法線だけをサンプリングすることは、物体表面の細かな凹凸による遮蔽の影響をサンプリング時に考慮しているものとみなせます。Walterらの手法では式(24)に示したように重みの分母に入射方向と法線の内積およびMicrofacet法線と面法線の内積が表れます。これらの内積は限りなく0に近い値を取る可能性があり、重みがとても大きな値になりえます。そうするとレンダリング結果に輝点が目立つようになるため、この項をサンプリング時に考慮できれば、より輝点の少ない画像が得られるだろう、というのが可視法線サンプリングの狙いです [Heitz et al. 2014b]。

可視法線だけをサンプリングするために、Microfacet法線の方向と光の入射方向の内積により重みづけた可視法線分布関数D_{\omega_i}(\omega_m)を考えます。

(25)    \begin{equation*} D_{\omega_i}(\omega_m) = \frac{ \chi^{+} (\omega_i \cdot \omega_m) (\omega_i \cdot \omega_m) D(\omega_m)}{\int_\Omega \langle \omega_i, \omega_m \rangle D(\omega_m) d \omega_m} \end{equation*}

上記の分布をスロープ空間と呼ばれる空間でサンプリングするのですが、ここでいうスロープというのは物体表面の細かな凹凸による表面の傾きのことです。物体表面の法線を\omega_m = (x_m, y_m, z_m)とするとき、スロープ\tilde{m}は二次元の変数で、\tilde{m} = (x_{\tilde{m}}, y_{\tilde{m}}) = (-x_m / z_m, -y_m / z_m) と表します。

このスロープの分布は法線分布の関数を変数変換することで、

(26)    \begin{eqnarray*} P^{22}(x_{\tilde{m}}, y_{\tilde{m}}) &=& \left\| \frac{\partial \omega_m}{\partial \tilde{m}} \right\| D(\omega_m) (\omega_m \cdot n) \\ &=& (\omega_m \cdot n)^4 D(\omega_m) = \cos^4 \theta_m D(\omega_m) \end{eqnarray*}

と書くことができます。入射方向に基づいた法線分布関数の場合にはMicrofacet法線と面法線の関係が式(25)の分子ですでに考慮されているので、D_{\omega_i}(\omega_m)は以下のように求められます。

(27)    \begin{equation*} P^{22}_{\omega_i}(x_{\tilde{m}}, y_{\tilde{m}}) = \cos^3 \theta_m D_{\omega_i}(\omega_m) \end{equation*}

と表せます。ここで式(26)を用いて、上記のスロープ空間における分布を展開すると、

(28)    \begin{eqnarray*} P^{22}_{\omega_i}(x_{\tilde{m}}, y_{\tilde{m}}) &=& \frac{1}{c} \langle \omega_i, \omega_m \rangle D(\omega_m) \\ &=& \frac{1}{c} \chi^{+}(x_i x_{\tilde{m}} + y_i y_{\tilde{m}} + z_i) (x_i x_{\tilde{m}} + y_i y_{\tilde{m}} + z_i) P^{22} (x_{\tilde{m}}, y_{\tilde{m}})  \end{eqnarray*}

と書けます。ただし、上記の式においてcは以下のような規格化定数を表します。

(29)    \begin{equation*} c = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \chi^{+}(x_i x_{\tilde{m}} + y_i y_{\tilde{m}} + z_i) (x_i x_{\tilde{m}} + y_i y_{\tilde{m}} + z_i) P^{22} (x_{\tilde{m}}, y_{\tilde{m}}) d x_{\tilde{m}} d y_{\tilde{m}} \end{equation*}

このスロープ分布関数を上手くサンプリングすることができれば、得られたスロープから一意に法線が決まりそうです。実際にはスロープの分布は表面のラフネスに依存していますので、以後はパラメータ空間の各軸方向におけるラフネス\alphax_, \alpha_yを明示してP^{22} (x_{\tilde{m}}, y_{\tilde{m}}, \alpha_x, \alpha_y)のように書きなおします。

一方で、ラフネスを変数とした場合には、上記の関数はサンプリングが難しいので、スロープ分布のスケール非依存性 (scale invariance) を利用してサンプリングがし易い形に持っていきます。スケール非依存性というのは以下のような性質を指します。

(30)    \begin{equation*} P^{22} (x_{\tilde{m}}, y_{\tilde{m}}, \alpha_x, \alpha_y) = \frac{1}{\lambda_x \lambda_y} P^{22} \left( \frac{x_{\tilde{m}}}{\lambda_x}, \frac{y_{\tilde{m}}}{\lambda_y}, \frac{\alpha_x}{\lambda_x}, \frac{\alpha_y}{\lambda_y} \right)  \end{equation*}

この式で\lambda_x, \lambda_yはスケーリングのための変数で自由に設定が可能です。ということは\lambda_x = \alpha_x, \lambda_y = \alpha_yと設定することで、後半の2つの引数を1にすることができます。つまり、ラフネスが1の場合にのみ上手くスロープ分布をサンプリングできれば良いことになります。

この性質は、可視法線だけを考えた場合にもそのまま受け継がれます。実際には入射方向をスケーリングする操作を

(31)    \begin{equation*} \mathcal{S}^{\lambda_x, \lambda_y} = \text{normalize} \left(\lambda_x x_i, \lambda_y y_i, z_i \right) \end{equation*}

として、可視法線に対するスロープの分布は以下のように書き直せます。

(32)    \begin{equation*} P_{\omega_i}^{22} (x_{\tilde{m}}, y_{\tilde{m}}, \alpha_x, \alpha_y) = \frac{1}{\lambda_x \lambda_y} P_{\mathcal{S}^{\lambda_x, \lambda_y} (\omega_i)}^{22} \left( \frac{x_{\tilde{m}}}{\lambda_x}, \frac{y_{\tilde{m}}}{\lambda_y}, \frac{\alpha_x}{\lambda_x}, \frac{\alpha_y}{\lambda_y} \right)  \end{equation*}

このスケーリング操作のことを元論文である[Hetiz et al. 2014]ではStretchと呼んでいます。

Stretchingにより、より簡易な形にもっていくことができましたが、サンプリングのためには、さらなる工夫が必要です。上記の可視法線に対するスロープの分布は入射方向に依存する形になっていますので、その依存関係を少しだけ緩めます。実際には、サンプリングするMicrofacet法線を面法線と入射方向が張る局所座標系の中でサンプリングして、それを法線を軸として回転させることでワールド座標系における法線のサンプリングを行います。

具体的には、入射方向が \omega_i = (x_i, 0, z_i)のようになっていると考えてサンプリングし、これを\phi_iだけ法線を軸として回転させることで実際のMicrofacet法線をサンプリングします。このような操作が可能になるのは、ラフネスがどの方向においても1になるようにStretchの操作をしているためです。この法線を回転する操作を[Heitz et al. 2014b]の論文ではrotateと呼んでいます。

これらをまとめるとスロープ表現を用いた可視法線のサンプリングは以下のような流れになります ([Heitz et al. 2014]のAlgorithm 4と対応)。

Step 1. 入射方向のスケーリング
 \omega_i' = \mathcal{S}^{\alpha_x, \alpha_y}(\omega_i)

Step 2. スロープのサンプリング
 (x_{\tilde{m}}', y_{\tilde{m}}') \sim \frac{1}{\alpha_x \alpha_y} P_{\omega_i'}^{22} (x_{\tilde{m}}', y_{\tilde{m}}', 1, 1)

Step 3. スロープ方向の回転
 \begin{pmatrix} x_{\tilde{m}}'' \ y_{\tilde{m}}'' \end{pmatrix} = \begin{pmatrix} \cos \phi_i & -\sin \phi_i \ \sin \phi_i \cos \phi_i \end{pmatrix} \begin{pmatrix} x_{\tilde{m}}' \ y_{\tilde{m}}' \end{pmatrix}

Step 4. スロープ方向のスケーリングを打ち消す
 x_{\tilde{m}} = \alpha_x x_{\tilde{m}}'', y_{\tilde{m}} = \alpha_y y_{\tilde{m}}''

Step 5. スロープから法線を計算
 \omega_m = \text{normalize} (-x_{\tilde{m}}, -y_{\tilde{m}}, 1)
 

Step 2.のスロープのサンプリングについてはx_{\tilde{m}}'y_{\tilde{m}}'のそれぞれで、以下の要領でサンプリング公式を得ます。

  • P^{22}_{\omega_i}x_{\tilde{m}} に関する周辺分布P^{2-}_{\omega_i}(x_{\tilde{m}}, 1, 1) および、その累積密度分布 C^{2-}_{\omega_i}(x_{\tilde{m}}, 1, 1) を求めて x_{\tilde{m}} のサンプリング公式を得る
  • x_{\tilde{m}}を定数として扱うことでy_{\tilde{m}}に関する条件付き密度分布関数P^{2|2}_{\omega_i}(y_{\tilde{m}}, 1, 1)および、その累積密度分布関数C^{2|2}_{\omega_i}(y_{\tilde{m}}, 1, 1)を求めてy_{\tilde{m}}のサンプリング公式を得る

以下ではBeckmann分布とGGX分布のそれぞれで (x_{\tilde{m}}', y_{\tilde{m}}') のサンプリング公式を計算します。
 

Beckmann分布の場合

Beckmann分布において\alpha_x = 1, \alpha_y = 1と置くと

(33)    \begin{equation*} D(\omega_m) = \frac{1}{\pi \cos^4 \theta_m} \exp(-\tan^2 \theta_m) \end{equation*}

と簡略化できます。これに式(27)の変数変換を施すとスロープ分布関数が以下のように得られます。

(34)    \begin{equation*} P^{22}(x_{\tilde{m}}, y_{\tilde{m}}) = \frac{1}{\pi} \exp(-x_{\tilde{m}}^2 - y_{\tilde{m}}^2) \end{equation*}

この式に対してy_{\tilde{m}}を定義域全体で積分して周辺化するとx_{\tilde{m}}だけの分布が以下のように得られます。

(35)    \begin{eqnarray*} P^{2-}(x_{\tilde{m}}, 1, 1) &=& \int_{-\infty}^{\infty} P^{22}(x_{\tilde{m}}, y_{\tilde{m}}, 1, 1) d y_{\tilde{m}} \\ &=& \frac{1}{\pi} \exp(-x_{\tilde{m}}^2) \int_{-\infty}^{\infty} \exp(-y_{\tilde{m}}^2) d y_{\tilde{m}} \\  &=& \frac{1}{\sqrt{\pi}} \exp(-x_{\tilde{m}}^2)  \end{eqnarray*}

この式とy_i = 0という仮定を用いるとP^{2-}(x_{\tilde{m}}, 1, 1)は以下のように表されます。

(36)    \begin{eqnarray*} P^{2-}_{\omega_i} (x_{\tilde{m}}, 1, 1) &=& \frac{1}{c} \int_{-\infty}^{\infty} \chi^{+}(x_i x_{\tilde{m}} + z_i)(x_i x_{\tilde{m}} + z_i) P^{22} (x_{\tilde{m}}, y_{\tilde{m}}, 1, 1) d y_{\tilde{m}} \\ &=& \frac{1}{c} \chi^{+} (x_i x_{\tilde{m}} + z_i) (x_i x_{\tilde{m}} + z_i) \int_{-\infty}^{\infty} P^{22}(x_{\tilde{m}}, y_{\tilde{m}}, 1, 1) d y_{\tilde{m}} \\ &=& \frac{1}{c} \chi^{+} (x_i x_{\tilde{m}} + z_i) (x_i x_{\tilde{m}} + z_i) P^{2-}(x_{\tilde{m}}, 1, 1) \end{eqnarray*}

上記の式は式(6)に示した関数G_1を用いて以下のように書き換えられます。

(37)    \begin{equation*} P^{2-}(x_{\tilde{m}}, 1, 1) = \frac{G_1(\omega_i)}{cos \theta_i} \chi^{+} (x_i x_{\tilde{m}} + z_i) (x_i x_{\tilde{m}} + z_i) P^{2-} (x_{\tilde{m}}) \end{equation*}

ここから更にx_{\tilde{m}}に関する累積密度分布を求めます。

(38)    \begin{eqnarray*} C^{2-}_{\omega_i}(x_{\tilde{m}}, 1, 1) &=& \frac{G_1(\omega_i)}{\cos \theta_i} \int_{\infty}^{x_{\tilde{m}}} \chi^{+} (x_i x_{\tilde{m}}' + z_i) (x_i x_{\tilde{m}}' + z_i) P^{2-} (x_{\tilde{m}}) d x_{\tilde{m}}' \\ &=& \frac{G_1(\omega_i)}{\sqrt{\pi} \cos \theta_i} \left[ x_i \int_{-\infty}^{x_{\tilde{m}}} x_{\tilde{m}}' \exp (-x_{\tilde{m}}'^2) d x_{\tilde{m}}' + z_i \int_{-\infty}^{\infty} \exp(-x_{\tilde{m}}'^2) d x_{\tilde{m}}' \right] \\ &=& \frac{G_1(\omega_i)}{\sqrt{\pi} \cos \theta_i} \left[ (-\frac{1}{2} ) (-\sin \theta_i) \int_{-\infty}^{x_{\tilde{m}}} (-2 x_{\tilde{m}}) \exp(-x_{\tilde{m}}'^2) d x_{\tilde{m}}' + \cos \theta_i \int_{-\infty}^{x_{\tilde{m}}} \exp(-x_{\tilde{m}}) d x_{\tilde{m}}' \right] \\ &=& \frac{G_1(\omega_i)}{\sqrt{\pi} \cos \theta_i} \left[ \left. \frac{\sin \theta_i}{2} \exp(-x_{\tilde{m}}'^2) \right\vert_{x_{\tilde{m}}' = -\infty}^{x_{\tilde{m}}' = x_{\tilde{m}}} + \left. \cos \theta_i \frac{\sqrt{\pi}}{2} \text{erf} (x_{\tilde{m}}') \right\vert_{x_{\tilde{m}}' = -\infty}^{x_{\tilde{m}}' = x_{\tilde{m}}}  \right] \\ &=& \frac{G_1(\omega_i) \tan \theta_i}{\sqrt{\pi}} \exp(-x_{\tilde{m}}^2) + G_1(\omega_i) \left( \frac{1}{2} + \frac{1}{2} \text{erf}(x_{\tilde{m}}) \right)  \end{eqnarray*}

本来ならば、この累積密度関数から逆関数法によりx_{\tilde{m}} のサンプリング公式を得たいのですが、上記の関数は逆関数が解析的にも止まらないので、実装の際には二分法でC^{2-}_{\omega_i}(w_{tilde{m}}, 1, 1) = u_1, u_1 \sim \mathcal{U}[0, 1] を解いて、数値的に逆関数の値を評価します。詳細は後ほど説明するプログラムをご参照ください。

これで x_{\tilde{m}}がサンプリングできましたので続いては y_{\tilde{m}} のサンプリングについて考えます。このとき、入射方向 \omega_i において y_i = 0 という仮定をおいていますので、y_{\tilde{m}} のサンプリングはP^{22}の条件付き分布を考えれば良いことになります。

(39)    \begin{eqnarray*} P_{\omega_i}^{2|2}(y_{\tilde{m}} | x_{\tilde{m}}, 1, 1) &=& P_{\omega_i}^{2|2}(y_{\tilde{m}} | x_{\tilde{m}}, 1, 1) \\ &=& \frac{P^{22}(x_{\tilde{m}}, y_{\tilde{m}}, 1, 1)}{\int_{-\infty}^{\infty} P^{22}(x_{\tilde{m}}, y_{\tilde{m}}, 1, 1) d y_{\tilde{m}}} \\ &=& \frac{1}{\pi} \exp(-y_{\tilde{m}}^2) \cdot \left(\sqrt{\pi}}\right)^{-1} \\ &=& \frac{1}{\sqrt{\pi}} \exp(-y_{\tilde{m}}^2) \end{eqnarray*}

また、この分布の累積密度分布は

(40)    \begin{equation*} C_{\omega_i}^{2|2}(y_{\tilde{m}}) = \frac{1}{\sqrt{\pi}} \int_{-\infty}^{y_{\tilde{m}}} \exp(-y_{\tilde{m}}'^2) d y_{\tilde{m}}' = \frac{1}{2} + \frac{1}{2} \text{erf} (y_{\tilde{m}}) \end{equation*}

のようになります。従って、y_{\tilde{m}} のサンプリング公式は以下のようになります。

(41)    \begin{equation*} y_{\tilde{m}} = \text{erf}^{-1}(2 u_2 - 1) \qquad u_2 \sim \mathcal{U}[0, 1] \end{equation*}

エラー関数 \text{erf}(y_{\tilde{m}}) の逆関数は解析的に求まりませんので、実装においては二分探索で逆関数を評価することとします。ちなみにPBRTの実装ではエラー関数の逆関数が多項式の分数により近似 (Rational approximation) されていますので、今回もその方法を使います。

以上でBeckmann分布からスロープをサンプリングすることができるようになりました。
 

GGX分布の場合

続いてはGGX分布です。GGX分布において\alpha_x = 1, \alpha_y = 1とおくと以下のように簡略化できます。

(42)    \begin{equation*} D(\omega_m) = \frac{1}{\pi} \frac{1}{(1 + x_{\tilde{m}}^2 + y_{\tilde{m}}^2)^2} \end{equation*}

この形を使いD_{\omega_i}x_{\tilde{m}}に関する周辺分布を求めます。この際、Beckmann分布のときの計算と同様にy_i = 0という仮定を利用すると、GGX分布に対するP^{2-}(x_{\tilde{m}}, 1, 1)は以下のように計算できます。

(43)    \begin{eqnarray*} P^{2-} (x_{\tilde{m}}) &=& \int_{-\infty}^{\infty} P^{22}(x_{\tilde{m}}, y_{\tilde{m}}, 1, 1) d y_{\tilde{m}} \\ &=& \frac{1}{\pi} \int_{-\infty}^{\infty} \frac{1}{(1 + x_{\tilde{m}}^2 + y_{\tilde{m}}^2)^2} d y_{\tilde{m}} \\ &=& \frac{2}{\pi} \int_{0}^{\infty} \frac{1}{(1 + x_{\tilde{m}}^2 + y_{\tilde{m}}^2)^2} d y_{\tilde{m}} \\ &=& \frac{2}{\pi} \left. \frac{1}{2} \left( \frac{y_{\tilde{m}}}{(1 + x_{\tilde{m}}^2)(1 + x_{\tilde{m}}^2 + y_{\tilde{m}}^2)} + \frac{\arctan \left( \frac{y_{\tilde{m}}}{\sqrt{1 + x_{\tilde{m}}^2}} \right)}{(1 + x_{\tilde{m}}^2)^{3/2}} \right) \right\vert_{y_{\tilde{m}} = 0}^{y_{\tilde{m}} = \infty} \\ &=& \frac{1}{2 (1 + x_{\tilde{m}}^2)^{3/2}} \end{eqnarray*}

式途中の定積分はWalframAlphaを使って計算していますが、もし手計算をしたい場合にはt = y_{\tilde{m}} / \sqrt{1 + x_{\tilde{m}}^2}とおくと正しく計算できます。Beckmann分布の計算の時に示した式(37)は分布の内容に依らない式ですので、こちらをGGX分布に対しても適用するとP^{2-}(x_{\tilde{m}}, 1, 1)は以下のように表されます。

(44)    \begin{eqnarray*} P^{2-}_{\omega_i}(\omega_m) = \frac{G_1(\omega_i)}{\cos \theta_i} \chi^{+} (x_i x_{\tilde{m}} + z_i) (x_i x_{\tilde{m}} + z_i) \frac{1}{2 (1 + x_{\tilde{m}}^2)^{3/2}} \end{eqnarray*}

ここからさらに累積密度分布を求めると以下のようになります。

(45)    \begin{eqnarray*} C_{\omega_i}(x_{\tilde{m}}) &=& \int_{-\infty}^{x_{\tilde{m}}} \frac{G_1(\omega_i)}{\cos \theta_i} \chi^{+} (x_i x_{\tilde{m}}' + z_i) (x_i x_{\tilde{m}}' + z_i) \frac{1}{2 (1 + x_{\tilde{m}}'^2)^{3/2}} d x_{\tilde{m}}' \\ &=& \frac{G_1(\omega_i)}{\cos \theta_i} \left[ x_i \int_{-\infty}^{x_{\tilde{m}}} \frac{x_{\tilde{m}}'}{2 (1 + x_{\tilde{m}}'^2)^{3/2}} d x_{\tilde{m}}' + z_i \int_{-\infty}^{x_{\tilde{m}}} \frac{1}{2 (1 + x_{\tilde{m}}'^2)^{3/2}} d x_{\tilde{m}}' \right] \\ &=& \frac{G_1(\omega_i)}{\cos \theta_i} \left[ (-\sin \theta_i) \left. \left(- \frac{1}{2  \sqrt{1 + x_{\tilde{m}}'^2}}} \right) \right\vert_{x_{\tilde{m}}' = -\infty}^{x_{\tilde{m}}' = x_{\tilde{m}}} + \left. \cos \theta_i \frac{x_{\tilde{m}}'}{2 \sqrt{1 + x_{\tilde{m}}'^2}} \right\vert_{x_{\tilde{m}}' = \infty}^{x_{\tilde{m}}' = x_{\tilde{m}}} \right] \\ &=& \frac{G_1(\omega_i)}{2 \cos \theta_i} \left[ \frac{\sin \theta_i}{\sqrt{1 + x_{\tilde{m}}}} + \frac{x_{\tilde{m}} \cos \theta_i}{\sqrt{1 + x_{\tilde{m}}}} + \cos \theta_i  \right] \end{eqnarray*}

あとは逆関数法より、サンプリング公式を得ます。今回は式が少し複雑ですので、以下のように式を変換していきます。

(46)    \begin{eqnarray*} u_1 &=& \frac{G_1(\omega_i)}{2 \cos \theta_i} \left[ \frac{\sin \theta_i}{\sqrt{1 + x_{\tilde{m}}^2}} + \frac{x_{\tilde{m}} \cos \theta_i}{\sqrt{1 + x_{\tilde{m}}^2}} + \cos \theta_i  \right] \qquad u_1 \sim \mathcal{U}[0, 1] \\ \frac{2 u_1}{G_1} - 1 &=& \frac{\tan \theta_i + x_{\tilde{m}}}{\sqrt{1 + x_{\tilde{m}}^2}} \\ \end{eqnarray*}

これを仮に

(47)    \begin{equation*} A = \frac{B + x_{\tilde{m}}}{\sqrt{1 + x_{\tilde{m}}^2}} \end{equation*}

とおいてx_{\tilde{m}}について解くと、

(48)    \begin{eqnarray*} x_{\tilde{m}} &=& \begin{cases} x_{\tilde{m}1} \qquad A < 0 \quad\text{or}\quad x_{\tilde{m}2} > \cot \theta_i \\ x_{\tilde{m}2} \qquad \text{otherwise} \end{cases} \\ x_{\tilde{m}1} &=& \frac{B}{A^2 - 1} - \sqrt{\frac{B^2}{(A^2 - 1)^2} - \frac{A^2 - B^2}{A^2 - 1}} \\ x_{\tilde{m}2} &=& \frac{B}{A^2 - 1} + \sqrt{\frac{B^2}{(A^2 - 1)^2} - \frac{A^2 - B^2}{A^2 - 1}} \end{eqnarray*}

となります。これでGGX分布についてもx_{\tilde{m}}がサンプリングできます。

続いてBeckmann分布の時と同様にy_{\tilde{m}}についてもサンプリング公式を出したいのですが、GGX分布の場合には条件付きの累積分布関数C^{2|2}(y_{\tilde{m}} | x_{\tilde{m}})は解析的に求まりません。そこで、変数変換を行い計算負荷の少ない形に変形します。

(49)    \begin{eqnarray*} C_{\omega_i}^{2|2} &=& \frac{ \int_{-\infty}^{y_{\tilde{m}}} P^{22}(x_{\tilde{m}}, y_{\tilde{m}}, 1, 1) d y_{\tilde{m}} }{ \int_{-\infty}^{\infty} P^{22}(x_{\tilde{m}}, y_{\tilde{m}}, 1, 1) d y_{\tilde{m}} } \\ &=& \frac{ \int_{-\infty}^{y_{\tilde{m}}} \frac{1}{(1 + x_{\tilde{m}}^2 + y_{\tilde{m}}^2)^2} d y_{\tilde{m}} }{ \int_{-\infty}^{\infty} \frac{1}{(1 + x_{\tilde{m}}^2 + y_{\tilde{m}}^2)^2} d y_{\tilde{m}} } \\ &=& \frac{ \int_{-\infty}^{y_{\tilde{m}}} \frac{1}{(a + y_{\tilde{m}}^2)^2} d y_{\tilde{m}} }{ \int_{-\infty}^{\infty} \frac{1}{(a + y_{\tilde{m}}^2)^2} d y_{\tilde{m}} } \\ &=& \frac{\displaystyle \int_{-\infty}^{y_{\tilde{m}}} \frac{1}{\left(1 + \frac{y_{\tilde{m}}^2}{a} \right)^2} d y_{\tilde{m}} }{\displaystyle \int_{-\infty}^{\infty} \frac{1}{\left(1 + \frac{y_{\tilde{m}}^2}{a} \right)^2} d y_{\tilde{m}} } \end{eqnarray*}

式の途中でa = 1 + x_{\tilde{m}}と置き換えていることにご注意ください。ここでさらにz' = y_{\tilde{m}} / \sqrt{a}と変数変換すると、

(50)    \begin{eqnarray*} C^{22}(y_{\tilde{m}} | x_{\tilde{m}}) &=& C_z(z) \\ &=& \frac{\int_{-\infty}^{z} \frac{1}{(1 + z'^2)^2} dz'}{\int_{-\infty}^{\infty} \frac{1}{(1 + z'^2)^2} dz'} \end{eqnarray*}

当然ながらこの式中のC_z(z)も解析的には求まらないので逆関数を求めることはできないのですが、これを多項式の分数として近似すると以下のような式が得られます [Heitz et al. 2014b]。

(51)    \begin{equation*} C_z^{-1}(\eta) = \frac{0.46341 \eta - 0.73369 \eta^2 + 0.27385 \eta^3}{0.597999 - \eta + 0.309420 \eta^2 + 0.0934073 \eta^3} \end{equation*}

これを利用すれば、y_{\tilde{m}} のサンプリング公式が以下のように得られます。

(52)    \begin{equation*} y_{\tilde{m}} = C_z^{-1}(u_2) \sqrt{1 + x_{\tilde{m}}^2} \qquad u_2 \sim \mathcal{U}[0, 1] \end{equation*}

Shockerさん より[Heitz 2016a]のテクニカルレポートにて逆関数のフィッティングを使わない方法が紹介されているとの助言をいただきました。

重みの計算

上記の可視法線のサンプリングでは[Walter et al. 2007]とは違い  (\omega_m \cdot n) を掛け算しておらず、また入射方向に依存する成分である  G_1(\omega_i, \omega_m) |\omega_i \cdot \omega_m| / |\omega_i \cdot n| が考慮されています。以上から、式(24)の重みは以下のように書き直せます。

(53)    \begin{equation*} \text{weight}(\omega_o) = G_1(\omega_i, \omega_o) \end{equation*}

最後にG_1をBeckmann分布、GGX分布の両方について求めておきます。式(7)に現れる\Lambda(\omega_o)はスロープ表現を用いると積分の形で以下のように示されます (導出については [Heitz et al. 2014a]のAppendix Aを参照)。

(54)    \begin{equation*} \Lambda(\omega_o) = \frac{1}{\cot \theta_o} \int_{\cot \theta_o}^{\infty} (x_{\tilde{m}} - \cot \theta_o) P^{2-} (x_{\tilde{m}}) d x_{\tilde{m}} \end{equation*}

この式に式(35)および式(43)を代入すると\Lambda (\omega_o)はBeckmann分布、GGX分布のそれぞれに対して以下のように求まります。

Beckmann分布の場合

(55)    \begin{equation*} \Lambda(\omega_o) = \frac{\text{erf}(a) - 1}{2} + \frac{1}{2 a \sqrt{\pi}} \exp(-a^2)} \end{equation*}

GGX分布の場合

(56)    \begin{equation*} \Lambda(\omega_o) = \frac{-1 + \sqrt{1 + \frac{1}{a^2}}}{2} \end{equation*}

いずれの場合も式中のaは以下のような式で表されるものです。

(57)    \begin{eqnarray*} a &=& \frac{1}{\alpha \tan \theta_o} \\ \alpha &=& \sqrt{\alpha_x^2 \cos^2 \phi_o + \alpha_y^2 \sin \phi_o} \end{eqnarray*}

これでようやく全ての式の説明が終わりました。それでは、これを使ってsmallpt を改良していきます。
 

Microfacetモデルのためのプログラム


ここからご紹介するプログラムは私のGitHub (https://github.com/tatsy/MicrofacetDistribution) から全体を確認できます。記事の中では断片だけを示しての説明のとなって今いますので、GitHubの方も合わせてご参照いただければ幸いです。

まずはMicrofacetモデルを表現するためのインターフェースを示します。

コンストラクタの引数に渡されている alphax および alphay はパラメータ座標系の各軸方向に対するラフネスを表し、sampleVis は光の入射方向から見えている法線だけをサンプリングに使うかどうかを表します。

こちらのインターフェースはMicrofacet法線 \omega_m をサンプルするための sampleWm と を \Lambda(\omega) を計算するのに使用する lambda の2つの関数だけが抽象化されており、このクラスを継承してBeckmann分布およびGGX分布における対応するメソッドを実装する仕組みにしています。重みの計算については可視法線をサンプルするかどうかだけに依存していますので、式(24)および式(53)を使用して定義をしています。

それでは、Beckmann分布とGGX分布のそれぞれについて実装の内容を見ていきます。

Beckmann分布のサンプリング

入射方向に依らない法線サンプリング

まずは光の入射方向に関係なくサンプリングする場合( sampleVisble = false の場合)をプログラムします。

このプログラムでは alphax == alphay が成り立つ場合とそうでない場合とで処理を分けています。
alphax == alphay のケースでは[Walter et al. 2007]に与えられているサンプリング公式をそのまま利用します。
今回の記事ではこちらのサンプリング公式は導出していませんが、非等方の場合の導出を参考にしていただければ簡単に導出できると思います。

続いては非等方 (alphax != alphay) の場合です。こちらのプログラムでは式(16)および式(18)を使用しています。
最初に式(16)を使って\phiをサンプルし、式(18)を使って \tan^2 \theta_m をサンプルしています。

プログラムの内容は以下の通りです。

 

可視法線分布によるサンプリング

可視法線をサンプリングする場合には記事の中で5つのステップに分けて説明した手順をまずプログラムします。そのソースコードが以下のようなものです。

 

続いてはスロープ分布から(x_{\tilde{m}}, y_{\tilde{m}})をサンプリングする処理です。この処理では式(38)の結果を使い、二分法によりx_{\tilde{m}}をサンプリングする処理、および式(41)を使いy_{\tilde{m}}をサンプリングする処理が含まれます。

なおx_{\tilde{m}}(-\infty, \infty)の範囲を動く変数であり、二分法と相性が悪いため、今回のプログラムでは\text{erf}(x_{\tilde{m}})を二分法により求めています。

 

以上で方向のサンプリングは可能ですので、最後に重み計算に使う\Lambda(\omega)をプログラムします。こちらには式(55)を使用します。

 

以上でBeckmann分布のプログラムは終了です。

GGX分布のサンプリング

入射方向に依らないサンプリング

こちらも最初は光の入射方向に関係なくサンプリングする場合( sampleVisble = false の場合)です。

Beckmann分布の時と同じく alphax == alphay が成り立つ場合とそうでない場合とで処理を分けています。
alphax == alphay のケースでは[Walter et al. 2007]に与えられているサンプリング公式をそのまま利用します。
GGX分布についても、非等方の場合の導出を参考にしていただければ簡単に導出できると思います。

続いては非等方 (alphax != alphay) の場合です。こちらのプログラムでは式(21)および式(23)を使用しています。
最初に式(21)を使って\phiをサンプルし、式(23)を使って \tan^2 \theta_m をサンプルしています。

プログラムの内容は以下の通りです。

 

可視法線分布によるサンプリング

可視法線分布をサンプルする場合にはGGX分布でも5つのステップで示したサンプリング処理を使います。この部分に関してはBeckmann分布のものと全く同じです。

 

続いてはスロープ分布のサンプリングです。ここでは式(48)および式(52)を使用します。

 

最後に\Lambda(\omega)です。こちらは式(56)を使用します。

 

これでGGXについても全てのプログラムが紹介し終わりました。
 

レンダリング結果


それでは作成したプログラムを実行してレンダリング結果を確認してみます。今回は1画素につき2000サンプルを取って画像を作成しました。画像を見るとラフネスが大きくなると[Walter et al. 2007]の方は相当輝点が出てしまっています。これに対してHeitzの方法はほとんど輝点がないので、概ね期待通りの結果と言っていいと思います。またBeckmann分布とGGX分布を見比べると同じラフネス値を使ってもBeckmann分布の方がやや拡散の効果が強くなるように見えます。非等方の場合に球の輪郭が明るくなってしまっているのは、もしかするとバグかもしれません。

\alpha_x = 0.1, \alpha_y=0.1 の場合 (クリックで拡大できます)

Beckmann分布(Walter+07) GGX分布(Walter+07) Beckmann分布(Heitz+14) GGX分布(Heitz+14)

 

\alpha_x = 0.5, \alpha_y=0.5 の場合 (クリックで拡大できます)

Beckmann分布(Walter+07) GGX分布(Walter+07) Beckmann分布(Heitz+14) GGX分布(Heitz+14)

 

\alpha_x = 0.5, \alpha_y=0.1 の場合 (クリックで拡大できます)

Beckmann分布(Walter+07) GGX分布(Walter+07) Beckmann分布(Heitz+14) GGX分布(Heitz+14)

 

計算時間

計算時間は可視法線分布かつBeckmann分布の場合を除いては概ね同じくらいでした。可視法線分布を使用してBeckmann分布をサンプリングするときには二分法の処理が入るので、その影響が強そうです。ただし、この二分法の部分については初期解のある程度良いものにしておくなどの改善も考えられるので、あくまで参考程度に考えていただければ幸いです。

PBRT v3では[Jakob 2014]のテクニカルレポートで紹介された高速な二分計算による逆関数の評価が使われているようです。

Beckmann分布(Walter+07) GGX分布(Walter+07) Beckmann分布(Heitz+14) GGX分布(Heitz+14)
475 秒 453 秒 700 秒 493 秒

※ Core i7 3.6 GHzのプロセッサで7スレッド並列にして 960×720 の画像をレンダリングした場合

まとめ


今回の記事では[Walter et al. 2007]および[Heitz et al. 2014b]で紹介されているMicrofacet法線のサンプリング公式の導出、およびプログラムへの組み込みをご紹介させていただきました。今年のレイトレ合宿では光沢面のレンダリングがたくさん見られるといいなと思っています(笑)。

本当はここからの発展として、多重散乱を考慮したSmithモデルの研究 [Heitz et al. 2016b] やBeckmann分布とGGX分布を一般化したStudentのt分布を使う方法 [Ribardiere et al. 2016] などを試したかったのですが、今回は記事の分量が多くなりすぎました。もし気が向いたらこの2つについても実装して公開しようと思っていますので、機会があればご覧いただければ嬉しいです。

それでは、長くなりましたが最後までお読みいただきありがとうございました (全部読んでくれる方はどのくらいいるのだろうか…?)。
 

参考文献


  • Walter et al. 2007, “Microfacet Models for Refraction through Rough Surfaces”, Eurographics Symposium on Rendering.
  • Heitz et al. 2014a, “Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs”, Journal of Computer Graphics Tools.
  • Heitz et al. 2014b, “Importance Sampling Microfacet-Based BSDFs with the Distribution of Visible Normals”, Eurographics Symposium on Rendering.
  • Heitz 2016a, “A Simpler and Exact Sampling Routine for the GGX Distribuion of Visible Normals”, Tech. Report.
  • Jakob 2014, “An improved visible normal sampling routine for the beckmann distribution”, Tech. Report.
  • Heitz et al. 2016b, “Multiple-Scattering Microfacet BSDFs with the Smith Model”, ACM Trans. Graph.
  • Ribardiere et al. 2016, “STD: Student’s t-Distribution of Slopes for Microfacet Based BSDFs”, Eurographics.

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です