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法線と呼び、その法線の分布を$D(\omega_m)$で表します。法線分布関数を用いると入射方向に対する物体表面での反射および透過の方向分布は次の式で表せます [Walter et al.2007]。

反射の場合 (BRDF)

$$ \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)

$$ \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} $$

この式で、$\omega_i$, $\omega_o$, $\omega_m$, $n$はそれぞれ、光の入射方向、光の出射方向、Microfacet法線、物体法線を表します。また$F(\omega_i, \omega_o)$はフレネル項の影響を$G(\omega_i, \omega_m, \omega_m)$はGeometry attenuation function (GAF)と呼ばれる、物体表面の凹凸による遮蔽の影響を表します。GAFについては、Smith-Bourliear GAFあるいはSmithのMasking-shadowing関数と呼ばれるものを使うのが一般的なようです。こちらのついては後ほどもう少し掘り下げて説明いたします。

レンダリング処理においては、物体法線方向の上半球面から入ってくる光を、上記の式で表されるBSDFの影響を考慮して積分します。この計算は以下のレンダリング方程式を用いて表されます。なお、以下の表示は[Walter et al. 2007]や[Heitz et al. 2014b]に基づいていて、$\omega_i$がカメラ側、$\omega_o$がライト側になっていることに注意してください。

$$ L(\omega_i) = \int_\Omega L(\omega_o) f(\omega_i, \omega_o) | \omega_o \cdot n | d \omega_o $$

Microfacetモデルを使う方向サンプリングの場合、離散的な積分を計算するためにサンプリングする対象はMicrofacet法線です。上記の積分は入射方向$\omega_i$での積分ですので、これをMicrofacet法線での積分に変数変換します。

この変数変換は反射の場合と透過の場合で異なってきます。[Walter et al. 2007]にある説明によると以下のようにヤコビアンが求まります([Walter et al. 2007]の図6および図7を参照)。

反射の場合

$$ \begin{equation} \left\| \frac{\partial \omega_m}{\partial \omega_o} \right\| = \frac{1}{4 | \omega_o \cdot \omega_m |} = \frac{1}{4 | \omega_i \cdot \omega_m |} \end{equation} $$

透過の場合

$$ \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} = \frac{\eta_i^2 | \omega_i \cdot \omega_m |}{(\eta_i (\omega_i \cdot \omega_m) + \eta_o (\omega_o \cdot \omega_m))^2} \end{equation} $$

BSDFの式ならびに上記の変数変換の式をレンダリング方程式に代入すると、反射の場合を例として、以下のように変形できます。

$$ \begin{aligned} L(\omega_i) &= \int_\Omega L(\omega_o) \frac{F(\omega_i, \omega_m) D(\omega_m) G(\omega_i, \omega_o, \omega_m)}{4 |\omega_i \cdot n| |\omega_o \cdot n|} | \omega_o \cdot n | \left\| \frac{\partial \omega_m}{\partial \omega_o} \right\|^{-1} d \omega_m \\ &= \int_\Omega L(\omega_o) \frac{F(\omega_i, \omega_m) D(\omega_m) G(\omega_i, \omega_o, \omega_m)}{4 |\omega_i \cdot n|} 4 | \omega_i \cdot \omega_m | d \omega_m \end{aligned} $$

この式を離散的に計算する際、$\omega_m$は法線分布関数$D(\omega_m)$に従ってサンプリングされるのですが、$D(\omega_m)$はその立体角の物体表面への射影が確率密度となるような分布関数なので、以下の式が成り立ちます。

$$ \int_\Omega D(\omega_m) | \omega_m \cdot n | d \omega_m = 1 $$

従って、$\omega_m$がサンプルされる確率$D(\omega_m) | \omega_m \cdot n |$を用いて、先ほどの積分式を離散化します。

$$ \begin{aligned} L(\omega_i) &= \sum_{\omega_m \sim D(\omega_m) | \omega_m \cdot n |} L(\omega_o) \frac{F(\omega_i, \omega_m) D(\omega_m) G(\omega_i, \omega_o, \omega_m)}{4 |\omega_i \cdot n|} 4 | \omega_i \cdot \omega_m | \frac{1}{D(\omega_m) | \omega_m \cdot n |} \\ &= \sum_{\omega_m \sim D(\omega_m) | \omega_m \cdot n |} L(\omega_o) F(\omega_i, \omega_m) \frac{ G(\omega_i, \omega_o, \omega_m) | \omega_i \cdot \omega_m |}{|\omega_i \cdot n| | \omega_m \cdot n |} \end{aligned} $$

以上の通り、実際のレンダリングでは、Microfacet法線を$D(\omega_m) |\omega_m \cdot n|$に従ってサンプリングし、入射光 $L(\omega_i)$にフレネル項と重みの項を掛け算して足し上げて行けば良いことが分かります。なお、この式は反射の場合の式ですが、透過の場合でも全く同じ式が得られます。

SmithのMasking-shadowing関数

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

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

この式は、式中の$\omega_i$ないし、$\omega_o$の方向から物体の表面を見た時に、物体表面の凹凸によって隠れてしまう面の影響を考慮しています。この時、$G_1$は法線分布関数$D(\omega_m)$と内積の符号に対する対する指示関数$\chi^{+}$を用いて以下のように書けます。

$$ \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’ } \label{eq:smith-g1} \end{equation} $$

この式を実際に計算していくと、

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

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

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

記事の冒頭で申し上げた通り、今回の記事では[Walter et al. 2007]と[Heitz et al. 2014b]の方法の2つの方法をご紹介いたします。いずれの方法の場合にも方向サンプリングのやり方はとても似ていて、Microfacet法線を通常の法線分布関数からサンプリングするか、光の入射方向$\omega_i$から見えているMicrofacet面上の法線の分布を表す可視法線分布関数からサンプリングするか、という違いがあるのみです。

Walter et al. 2007の場合

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

Heitz et al. 2014の場合

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

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

Microfacetモデルにより定まる光沢面 (Glossy surface) をBSDFにより表現する場合、 レンダリングの分野ではBeckmann分布とGGX分布 (Trowbridge-Reitz分布) という分布がよく使われます。

もし光沢面が光路方向に依らない等方的な法線分布を持っている場合にはBeckmann分布とGGX分布に対する方向サンプリングの公式は解析的に求まります。これが、この二つの分布がよく用いられる理由です。Beckmann分布とGGX分布はそれぞれ等方の場合、非等方の場合で以下のように定義されます。二つの分布は良く見るとBeckmann分布もGGX分布も似た形をしていることに気づきます。

Beckmann分布 (等方)

$$ \begin{equation} D(\omega_m) = \frac{\chi^+(\omega_m \cdot n)}{\pi \alpha^2 \cos^4 \theta} \exp \left( -\frac{\tan^2 \theta_m}{\alpha^2} \right) \end{equation} $$

Beckmann分布 (非等方)

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

GGX分布 (等方)

$$ \begin{equation} D(\omega_m) = \frac{\chi^+ (\omega_m \cdot n)}{\pi \alpha^2 \cos^4 \theta_m \left( 1 + \frac{\tan^2 \theta_m}{\alpha^2} \right)^2} \end{equation} $$

GGX分布 (非等方)

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

これらの分布の定義はご覧の通りかなり複雑ですが、実際に計算してみると解析的なサンプリング公式を得ることが可能です。等方的な分布については[Walter et al. 2007]でもサンプリング公式が示されていますし、導出方法もそれほど複雑ではありませんので、この記事では非等方の分布の方でサンプリングの公式を導出してみます。

Beckmann分布の場合

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

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

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

ただし、

$$ \frac{1}{A} = \frac{\chi^{+}(\omega_m \cdot n)}{\pi \alpha_x \alpha_y}, \quad B = \frac{\cos^2 \phi_n}{\alpha_x^2} + \frac{\sin^2 \phi_n}{\alpha_y^2} $$

となります。前述の通り、法線分布関数は$D(\omega_m) | \omega_m \cdot n |$を半球面上で積分すると、その値が1となるような確率密度分布になっているので、この値を積分していきます。

まずは、$\phi_m$に関する積分の公式を得るために、$\theta_m \in [0, \pi / 2]$を周辺化して変数から削除します。この積分は$\omega_m$に関する積分なので、これを$\theta_m$と$\phi_m$に関する積分に変換すると、ヤコビアンとして$\sin \theta_m$が現れます。また、局所座標系として$n = (0, 0, 1)$ならびに$\omega_m = (\cos \phi_m \sin \theta_m, \sin \phi_m \sin \theta_m, \cos \theta_m)$であり、それらの内積が$\cos \theta_m$になることに注意すると、$\theta_m$に関する積分は以下のようになります。

$$ \begin{eqnarray} \int_0^{\frac{\pi}{2}} \frac{1}{A \cos^4 \theta_m} \exp \left( -B \tan^2 \theta_m \right) \cos \theta_m \sin \theta_m d \theta_m &=& \int_0^{\frac{\pi}{2}} \frac{\sin \theta_m}{A \cos^3 \theta_m} \exp \left(-B \tan^2 \theta_m \right) d \theta_m \nonumber \\ &=& \int_0^{\frac{\pi}{2}} \left( -1 / (2 A B) \right) \left( \exp \left( -B \tan^2 \theta_m \right) \right)’ d \theta_m \nonumber \\ &=& \left. -\frac{1}{2AB} \exp \left(-B \tan^2 \theta_m \right) \right\vert_{\theta_m = 0}^{\theta_m=\frac{\pi}{2}} \nonumber \\ &=& \frac{1}{2AB} \nonumber \\ &=& \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} \label{eq:beckmann-theta-marginal} \end{eqnarray} $$

これで$\phi_m$に関する確率密度分布が得られましたが、逆関数法を用いたサンプリングでは累積密度分布が必要なので、さらに$\phi_m$に関する累積密度を積分により求めていきます。まずは$D(\phi_m)$を以下のように簡略化します。

$$ \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}$と置いて変数変換すると正しく計算できます。

$$ \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]$の範囲における一様乱数を表します。

$$ \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] \label{eq:beckmann-sample-phi} \end{equation} $$

次は$\theta_m$のサンプリング公式を求めていきますが、ここでは$\phi_m$がすでに決まっている前提で、以下の式で累積密度分布を求めます。

$$ \begin{equation} 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 } \end{equation} $$

これは式\eqref{eq:beckmann-theta-marginal}の途中式において$\theta_m = \frac{\pi}{2}$を代入するところをうまく利用すると以下のように計算できます。

$$ \begin{equation} C(\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{equation} $$

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

$$ \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] \label{eq:beckmann-sample-theta} \end{equation} $$

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

GGX分布の場合

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

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

この表現を使いると、$\phi_m$に関する周辺分布が以下のように求まります。

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

この結果は式\eqref{eq:beckmann-theta-marginal}と同じですので、Beckmann分布同様に$\phi_m$ のサンプリング公式は以下のようになります。

$$ \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] \label{eq:ggx-sample-phi} \end{equation} $$

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

$$ \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$ のサンプリング公式は逆関数法により以下のように求まります。

$$ \begin{equation} \theta_m = \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] \label{eq:ggx-sample-theta} \end{equation} $$

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

重みの計算

レンダリングの計算において$\omega_m$をサンプルした場合、そこから求めた$\omega_i$の方向に次のレイを飛ばします。この時、サンプルにかかる重みは、上記のレンダリング方程式による議論から以下のようになります (フレネル項は別途考慮します)。

$$ \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|} \label{eq:sample-weight} \end{equation} $$

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

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

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

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

可視法線だけをサンプリングするために、Microfacet法線の方向と光の入射方向の内積により重みづけた可視法線分布関数$D_{\omega_i}(\omega_m)$を考えます。これは$D(\omega_m)$に対して、Microfacet法線が$\omega_i$の方向から見えているかを表す指示関数$\chi^{+}(\omega_i \cdot \omega_m)$と、どのくらい見えているかの程度を表す$\omega_i \cdot \omega_m$を乗じた上で正規化すると、次のように定義できます。

$$ \begin{equation} D_{\omega_i}(\omega_m) = \frac{ \chi^{+} (\omega_i \cdot \omega_m) (\omega_i \cdot \omega_m) D(\omega_m)}{\int_\Omega (\omega_i \cdot \omega_m) D(\omega_m) d \omega_m} \label{eq:visible-normal-distrib} \end{equation} $$

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

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

$$ \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*} $$

と書くことができます。一方で、可視性を考慮したスロープの分布は、式\eqref{eq:visible-normal-distrib}と同様にして、次のように定義されます。

$$ \begin{equation} P^{22}_{\omega_i}(x_{\tilde{m}}, y_{\tilde{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{equation} $$

ただし、上記の式において規格化定数$c$は以下のようになります。

$$ \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} $$

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

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

$$ \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の場合にのみ上手くスロープ分布をサンプリングできれば良いことになります。

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

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

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

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

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

Beckmann分布の場合

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

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

と簡略化できます。この式をスロープ空間へ変数変換すると

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

となります。この式に対して$x_{\tilde{m}}$だけの分布を得るために密度分布を周辺化すると、以下の式が得られます。

$$ \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}} \nonumber \\ &=& \frac{1}{\pi} \exp(-x_{\tilde{m}}^2) \int_{-\infty}^{\infty} \exp(-y_{\tilde{m}}^2) d y_{\tilde{m}} \nonumber \\ &=& \frac{1}{\sqrt{\pi}} \exp(-x_{\tilde{m}}^2) \label{eq:beckmann-slope-x-marginal} \end{eqnarray} $$

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

$$ \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*} $$

上記の式は式\eqref{eq:smith-g1}に示した$G_1$を用いて以下のように書き換えられます。

$$ \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}}) \label{eq:slope-p2-minus} \end{equation} $$

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

$$ \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}}’ \nonumber \\ &=& \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] \nonumber \\ &=& \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] \nonumber \\ &=& \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] \nonumber \\ &=& \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) \label{eq:beckmann-c2-minus} \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}$の条件付き分布を考えれば良いことになります。

$$ \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*} $$

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

$$ \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}}$ のサンプリング公式は以下のようになります。

$$ \begin{equation} y_{\tilde{m}} = \text{erf}^{-1}(2 u_2 - 1) \qquad u_2 \sim \mathcal{U}[0, 1] \label{eq:beckmann-sample-slope-y} \end{equation} $$

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

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

GGX分布の場合

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

$$ \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)$は以下のように計算できます。

$$ \begin{eqnarray} P^{2-} (x_{\tilde{m}}) &=& \int_{-\infty}^{\infty} P^{22}(x_{\tilde{m}}, y_{\tilde{m}}, 1, 1) d y_{\tilde{m}} \nonumber \\ &=& \frac{1}{\pi} \int_{-\infty}^{\infty} \frac{1}{(1 + x_{\tilde{m}}^2 + y_{\tilde{m}}^2)^2} d y_{\tilde{m}} \nonumber \\ &=& \frac{2}{\pi} \int_{0}^{\infty} \frac{1}{(1 + x_{\tilde{m}}^2 + y_{\tilde{m}}^2)^2} d y_{\tilde{m}} \nonumber \\ &=& \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} \nonumber \\ &=& \frac{1}{2 (1 + x_{\tilde{m}}^2)^{3/2}} \label{eq:ggx-slope-x-marginal} \end{eqnarray} $$

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

$$ \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*} $$

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

$$ \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*} $$

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

$$ \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*} $$

これを仮に

$$ \begin{equation} A = \frac{B + x_{\tilde{m}}}{\sqrt{1 + x_{\tilde{m}}^2}} \end{equation} $$

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

$$ \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} \label{eq:ggx-sample-slope-x} \\ 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}} \nonumber \\ 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}} \nonumber \end{eqnarray} $$

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

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

$$ \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}$と変数変換すると、

$$ \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]。

$$ \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}}$ のサンプリング公式が以下のように得られます。

$$ \begin{equation} y_{\tilde{m}} = C_z^{-1}(u_2) \sqrt{1 + x_{\tilde{m}}^2} \qquad u_2 \sim \mathcal{U}[0, 1] \label{eq:ggx-sample-slope-y} \end{equation} $$

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

重みの計算

SmithのShadowing-masking関数に現れる$G_1$を用いると、式\eqref{eq:visible-normal-distrib}から、可視法線分布関数は次のように書けます。

$$ D_{\omega_i}(\omega_m) = \frac{\chi^{+}(\omega_i \cdot \omega_m) (\omega_i \cdot \omega_m) D(\omega_m)}{\int_\Omega (\omega_i \cdot \omega_m) D(\omega_m) d \omega_m} = G_1(\omega_i, \omega_m) \frac{(\omega_i \cdot \omega_m)}{(\omega_i \cdot n)} D(\omega_m) $$

この確率密度を用いてレンダリング方程式を式変形すると、次のようになります。

$$ \begin{aligned} L(\omega_i) &= \sum_{\omega_m \sim D_{\omega_i}(\omega_m)} L(\omega_o) \frac{F(\omega_i, \omega_m) D(\omega_m) G(\omega_i, \omega_o, \omega_m)}{4 |\omega_i \cdot n|} 4 | \omega_i \cdot \omega_m | \left( G_1(\omega_i, \omega_m) \frac{(\omega_i \cdot \omega_m)}{(\omega_i \cdot n)} D(\omega_m) \right)^{-1} \\ &= \sum_{\omega_m \sim D_{\omega_i}(\omega_m)} L(\omega_o) F(\omega_i, \omega_m) G_1(\omega_o, \omega_m) \end{aligned} $$

したがって、可視法線分布関数を用いる場合のサンプルにかかる重みは次のように書けます。

$$ \begin{equation} \text{weight}(\omega_o) = G_1(\omega_o, \omega_m) \label{eq:sample-weight-visible} \end{equation} $$

最後に$\Lambda(\omega_o)$はスロープ表現を用いると積分の形で以下のように示されます (導出については [Heitz et al. 2014a]のAppendix Aを参照)。

$$ \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} $$

この式に式\eqref{eq:beckmann-slope-x-marginal}および式\eqref{eq:ggx-slope-x-marginal}を代入すると$\Lambda (\omega_o)$はBeckmann分布、GGX分布のそれぞれに対して以下のように求まります。

Beckmann分布の場合

$$ \begin{equation} \Lambda(\omega_o) = \frac{\text{erf}(a) - 1}{2} + \frac{1}{2 a \sqrt{\pi}} \exp(-a^2) \label{eq:beckmann-lambda} \end{equation} $$

GGX分布の場合

$$ \begin{equation} \Lambda(\omega_o) = \frac{-1 + \sqrt{1 + \frac{1}{a^2}}}{2} \label{eq:ggx-lambda} \end{equation} $$

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

$$ \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モデルを表現するためのインターフェースを示します。

// ----------------------------------------------------------------------------
// The interface for microfacet distribution
// ----------------------------------------------------------------------------
class MicrofacetDistribution {
public:
    MicrofacetDistribution(double alphax, double alphay, bool sampleVisible)
        : alphax_{alphax}
        , alphay_{alphay}
        , sampleVisible_{sampleVisible} {
    }

    virtual ~MicrofacetDistribution() {
    }

    // Sample microfacet normal
    virtual Vec sampleWm(const Vec &wo) const = 0;

    // The lambda function, which appears in the slope-space sampling
    virtual double lambda(const Vec &wo) const = 0;

    // Smith's masking function
    double G1(const Vec &wo) const {
        return 1.0 / (1.0 + lambda(wo));
    }

    // Smith's masking-shadowing function
    double G(const Vec &wo, const Vec &wi) const {
        return G1(wo) * G1(wi);
    }

    // Weighting function
    double weight(const Vec &wo, const Vec &wi, const Vec &wm) const {
        if (!sampleVisible_) {
            return std::abs(wi.dot(wm)) * G(wo, wi) /
                   std::max(1.0e-8, std::abs(cosTheta(wi) * cosTheta(wm)));
        } else {
            return G1(wo);
        }
    }

    // Private parameters
    double alphax_, alphay_;
    bool sampleVisible_;
};

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

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

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

Beckmann分布のサンプリング

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

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

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

続いては非等方 (alphax != alphay) の場合です。こちらのプログラムでは式\eqref{eq:beckmann-sample-phi}および式\eqref{eq:beckmann-sample-theta}を使用しています。最初に式\eqref{eq:beckmann-sample-phi}を使って$\phi_m$をサンプルし、その後で式\eqref{eq:beckmann-sample-theta}を使って$\tan^2 \theta_m$ をサンプルしています。

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

if (!sampleVisible_) {
    // Sample half-vector without taking normal visibility into consideration
    double tan2Theta, phi;
    if (alphax_ == alphay_) {
        // Isotropic case
        // Following smapling formula can be found in [Walter et al. 2007]
        // "Microfacet Models for Refraction through Rough Surfaces"
        double alpha = alphax_;
        double logSample = std::log(1.0 - U1);
        tan2Theta = -alpha * alpha * logSample;
        phi = 2.0 * Pi * U2;
    } else {
        // Anisotropic case
        // Following sampling strategy is analytically derived from 
        // P.15 of [Heitz et al. 2014]
        // "Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs"
        double logSample = std::log(1.0 - U1);
        phi = std::atan(alphay_ / alphax_ * std::tan(2.0 * Pi * U2 + 0.5 * Pi));
        if (U2 > 0.5) {
            phi += Pi;
        }

        double sinPhi = std::sin(phi);
        double cosPhi = std::cos(phi);
        double alphax2 = alphax_ * alphax_;
        double alphay2 = alphay_ * alphay_;
        tan2Theta = -logSample / (cosPhi * cosPhi / alphax2 + sinPhi * sinPhi / alphay2);
    }

    // Compute normal direction from angle information sampled above
    double cosTheta = 1.0 / std::sqrt(1.0 + tan2Theta);
    double sinTheta = std::sqrt(std::max(0.0, 1.0 - cosTheta * cosTheta));
    Vec wm(sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta);

    if (wm.z < 0.0) {
        wm = -wm;
    }

    return wm;
}  

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

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

// Sample microfacet normal through slope distribution
static Vec sampleBeckmann(const Vec &wi, double alphax, double alphay) {
    // 1. stretch wi
    Vec wiStretched = Vec(alphax * wi.x, alphay * wi.y, wi.z).normalize();

    // 2. simulate P22_{wi}(x_slope, y_slope, 1, 1)
    double slopex, slopey;
    sampleBeckman_P22_11(cosTheta(wiStretched), &slopex, &slopey);

    // 3. rotate
    double tmp = cosPhi(wiStretched) * slopex - sinPhi(wiStretched) * slopey;
    slopey = sinPhi(wiStretched) * slopex + cosPhi(wiStretched) * slopey;
    slopex = tmp;

    // 4. unstretch
    slopex = alphax * slopex;
    slopey = alphay * slopey;

    // 5. compute normal
    return Vec(-slopex, -slopey, 1.0).normalize();    
}

続いてはスロープ分布から$(x_{\tilde{m}}, y_{\tilde{m}})$をサンプリングする処理です。この処理では式\eqref{eq:beckmann-c2-minus}の結果を使い、二分法により$x_{\tilde{m}}$をサンプルする処理、ならびに式\eqref{eq:beckmann-sample-slope-y}を用いて$y_{\tilde{m}}$をサンプリングする処理が含まれます。

なお$x_{\tilde{m}}$は$(-\infty, \infty)$の範囲を動く変数であり、二分法と相性が悪いので、その代わりに$\text{erf}(x_{\tilde{m}})$を二分法により求めています。

// Sample slope distribution with alphax and alphay equal to one
static void sampleBeckman_P22_11(double cosThetaI, double *slopex, double *slopey) {
    using ::erf;

    const double U1 = random.next();
    const double U2 = random.next();

    // The special case where the ray comes from normal direction
    // The following sampling is equivarent to the sampling of
    // microfacet normals (not slopes) on isotropic rough surface
    if (cosThetaI > 0.9999) {
        double r = std::sqrt(-std::log(1.0f - U1));
        double sinPhi = std::sin(2 * Pi * U2);
        double cosPhi = std::cos(2 * Pi * U2);
        *slopex = r * cosPhi;
        *slopey = r * sinPhi;
        return;
    }

    double sinThetaI = std::sqrt(std::max((double)0, (double)1 - cosThetaI * cosThetaI));
    double tanThetaI = sinThetaI / cosThetaI;
    double cotThetaI = 1 / tanThetaI;
    const double erfCotThetaI = erf(cotThetaI);

    // Initialize lower and higher bounds for bisection method
    double lower = -1.0;
    double higher = 1.0;

    // Normalization factor for the CDF
    // This is equivarent to (G1 / 2)^{-1}
    static const double SQRT_PI_INV = 1.f / std::sqrt(Pi);
    double normalization = 1.0 / (1.0 + erfCotThetaI + SQRT_PI_INV * tanThetaI * std::exp(-cotThetaI * cotThetaI));

    // The following bisection method acquires "erf(x_m)"
    int it = 0;
    double samplex = std::max(U1, 1e-6);
    while (std::abs(higher - lower) > 1.0e-6) {
        // Bisection
        double mid = 0.5 * (lower + higher);

        // Evaluation for current estimation
        double x_m = erfinv(mid);
        double value = normalization * (1 + mid + SQRT_PI_INV * tanThetaI * std::exp(-x_m * x_m)) - samplex;

        /* Update bisection intervals */
        if (value > 0.0) {
            higher = mid;
        } else {
            lower = mid;
        }
    }

    // Compute slopex from erf(x_m) given by above bisection method
    *slopex = erfinv(lower);

    // Sample y_m
    *slopey = erfinv(2.0f * std::max(U2, 1e-6) - 1.0);
}

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

double lambda(const Vec &wo) const override {
    using ::erf;

    double cosThetaO = cosTheta(wo);
    double sinThetaO = sinTheta(wo);
    double absTanThetaO = std::abs(tanTheta(wo)); 
    if (std::isinf(absTanThetaO)) return 0.;

    double alpha = std::sqrt(cosThetaO * cosThetaO * alphax_ * alphax_ +
                             sinThetaO * sinThetaO * alphay_ * alphay_);
    double a = 1.0 / (alpha * absTanThetaO);
    return (erf(a) - 1.0) / 2.0 + std::exp(-a * a) / (2.0 * a * std::sqrt(Pi) + 1.0e-6);
}

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

GGX分布のサンプリング

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

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

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

続いては非等方 (alphax != alphay) の場合です。こちらのプログラムでは式\eqref{eq:ggx-sample-phi}および式\eqref{eq:ggx-sample-theta}を使用しています。最初に式\eqref{eq:ggx-sample-phi}を使って$\phi_m$をサンプルし、その後、式\eqref{eq:ggx-sample-theta}を使って$\tan^2 \theta_m$ をサンプルしています。

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

if (!sampleVisible_) {
    double tan2Theta, phi;
    if (alphax_ == alphay_) {
        // Isotropic case
        double alpha = alphax_;
        tan2Theta = alpha * alpha * U1 / (1.0 - U1);
        phi = 2.0 * Pi * U2;                
    } else {
        phi = std::atan(alphay_ / alphax_ * std::tan(2.0 * Pi * U2 + 0.5 * Pi));
        if (U2 > 0.5) {
            phi += Pi;
        }

        double sinPhi = std::sin(phi);
        double cosPhi = std::cos(phi);
        double alphax2 = alphax_ * alphax_;
        double alphay2 = alphay_ * alphay_;
        double alpha2 = 1.0 / (cosPhi * cosPhi / alphax2 + sinPhi * sinPhi / alphay2);
        tan2Theta = U1 / (1.0 - U1) * alpha2;
    }

    // Compute normal direction from angle information sampled above
    double cosTheta = 1.0 / std::sqrt(1.0 + tan2Theta);
    double sinTheta = std::sqrt(std::max(0.0, 1.0 - cosTheta * cosTheta));
    Vec wm = Vec(sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta);

    if (wm.z < 0.0) {
        wm = -wm;
    }

    return wm;
}

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

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

static Vec sampleGGX(const Vec &wi, double alphax, double alphay) {
    // 1. stretch wi
    Vec wiStretched = Vec(alphax * wi.x, alphay * wi.y, wi.z).normalize();

    // 2. simulate P22_{wi}(x_slope, y_slope, 1, 1)
    double slopex, slopey;
    sampleGGX_P22_11(cosTheta(wiStretched), &slopex, &slopey);

    // 3. rotate
    double tmp = cosPhi(wiStretched) * slopex - sinPhi(wiStretched) * slopey;
    slopey = sinPhi(wiStretched) * slopex + cosPhi(wiStretched) * slopey;
    slopex = tmp;

    // 4. unstretch
    slopex = alphax * slopex;
    slopey = alphay * slopey;

    // 5. compute normal
    return Vec(-slopex, -slopey, 1.0).normalize();       
}

続いてはスロープ分布のサンプリングです。ここでは式\eqref{eq:ggx-sample-slope-x}および式\eqref{eq:ggx-sample-slope-y}を使用します。

static void sampleGGX_P22_11(double cosThetaI, double *slopex, double *slopey) {
    const double U1 = random.next();
    double U2 = random.next();

    // The special case where the ray comes from normal direction
    // The following sampling is equivarent to the sampling of
    // microfacet normals (not slopes) on isotropic rough surface
    if (cosThetaI > 0.9999) {
        double r = std::sqrt(U1 / (1.0 - U1));
        double sinPhi = std::sin(2 * Pi * U2);
        double cosPhi = std::cos(2 * Pi * U2);
        *slopex = r * cosPhi;
        *slopey = r * sinPhi;
        return;
    }

    double sinThetaI = std::sqrt(std::max(0.0, 1.0 - cosThetaI * cosThetaI));
    double tanThetaI = sinThetaI / cosThetaI;
    double a = 1.0 / tanThetaI;
    double G1 = 2.0 / (1.0 + std::sqrt(1.0 + 1.0 / (a * a)));

    // Sample slopex
    double A = 2.0 * U1 / G1 - 1.0;
    double B = tanThetaI;
    double tmp = std::min(1.0 / (A * A - 1.0), 1.0e12);

    double D = std::sqrt(B * B * tmp * tmp - (A * A - B * B) * tmp);
    double slopex1 = B * tmp - D;
    double slopex2 = B * tmp + D;
    *slopex = (A < 0.0 || slopex2 > 1.0 / tanThetaI) ? slopex1 : slopex2;

    // Sample slopey
    double S;
    if (U2 > 0.5) {
        S = 1.0;
        U2 = 2.0 * (U2 - 0.5);
    } else {
        S = -1.0;
        U2 = 2.0 * (0.5 - U2);
    }

    double z = (U2 * (U2 * (U2 * 0.27385f - 0.73369f) + 0.46341f)) /
               (U2 * (U2 * (U2 * 0.093073f + 0.309420f) - 1.000000f) + 0.597999f);
    *slopey = S * z * std::sqrt(1.0 + (*slopex) * (*slopex));
}

最後に$\Lambda(\omega)$です。こちらは式\eqref{eq:ggx-lambda}を使用します。

double lambda(const Vec &wo) const override {
    double absTanThetaO = std::abs(tanTheta(wo)); 
    if (std::isinf(absTanThetaO)) return 0.;
    
    double alpha =
        std::sqrt(cos2Phi(wo) * alphax_ * alphax_ + sin2Phi(wo) * alphay_ * alphay_);
    double alpha2Tan2Theta = (alpha * absTanThetaO) * (alpha * absTanThetaO);
    return (-1.0 + std::sqrt(1.f + alpha2Tan2Theta)) / 2;
}

これで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]のテクニカルレポートで紹介された高速な二分計算による逆関数の評価が使われているようです。

※ 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.