初めに

レイトレ(Raytracing)のカレンダー | Advent Calendar 2023 - Qiitaの23日の記事です。通常のMicrofacet BSDFは何度も反射する多重散乱光(Multiple Scatteirng)を無視しているため、実はエネルギー損失が起きていることが知られています。この記事はこうしたMultiple Scatteringを考慮した[Heitz et al. 2016]の手法と最近発表されたFast Multiple Scattering Approximationという手法の実装について解説していきます。また、最近のMultiple Scatteringの動向について簡単に話します。

かなり込み入ったMicrofacetの話をするのでUnderstanding the Masking-Shadowing Function in Microfacet-Based BRDFs[Heitz. 2014]やUshioさんのMicrofacet入門(1)、または自著の「GGXは何処から来たのか」を見ておくことをお勧めします。

https://jcgt.org/published/0003/02/03/

ushiostarfish.hatenablog.com

motimoti3d.booth.pm

初めに
Single, Multiple Scatteringとは
[Heitz, et al. 2016]のMultiple-Scattering
Fast Multiple Scattering Approximation
Multiple Scatteringの動向
- Position Free
- エネルギー保存 BSDF
参考文献

Single, Multiple Scatteringとは

粗さを持つ物体表面に光が入射してきた時、その光は表面の小さな凹凸(Microfacet)によって様々な方向へと反射させられます。その反射の仕方は複雑で、単純に一度反射しただけでたまたま他のMicrofacetに当たることなく出ていける光もあれば、何度も何度もMicrofacetに当たっては反射を繰り返してようやく出射するような光もあります。

それぞれ1度だけ反射して外に出ていくという反射過程はSingle-Scatteringと言い、対称的に何度も反射する過程はMultiple-Scatteringと呼びます。何回反射した光なのかというのも明確に表現するときはNth bounceと呼んだりします。

CGで粗さを持つ物体を表現するにはこのような反射を何とか統計的に処理して、BSDFという形で表現する必要があります。ですが、皆さんは既に粗さを持つSpecular反射のBSDFを知っています。いつもの例の式です。

$\displaystyle{ f(\omega_i,\omega_o) = \frac{D(\omega_m) G(\omega_i,\omega_o) F(\omega_i,\omega_h) }{4 |\omega_i \cdot n| |\omega_o \cdot n| } }$

この式の呼称はTorrance-Sparrow ModelだとかCook-Torrance ModelだとかGGXだのあまり明確に定められていませんが、ここではMicrofacet BSDFと表現することとします。

Microfacet BSDFは今ではPBRであればオフライン、リアルタイム問わず様々なところで見かけるデファクトスタンダードと化していますが、実はMicrofacet BSDFには大きな問題点が存在します。

それはMicrofacet BSDFはSingle-Scatteringのみの反射を表現しているということです。Multiple-Scatteringによる光を一切存在しないという仮定のもと導出されたものであり、本来はもっとあるような光 2nd,3rd,... bounceの光をバッサリ捨てるような式になっています。そのため、式の性質を明確に示すのであればSingle-Scattering Microfacet BSDFと呼んだ方が適切だと思います(以降この呼称を使います)。

Single-Scattering Microfacet BSDFはこの性質のため、エネルギー保存則を満たしません。反射率であるフレネルを1としてWhiteFurnaceTestにかけてみると、背景とは溶け込まずにRoughnessが高いほど暗くなる＝エネルギー損失が起きていることが分かります。Roughnessが高いことはMicrofacetの凹凸の激しさを意味し、Multiple-Scatteringが起きやすいという性質を持つのでこのように仮定の影響をもろに受けてしまっているという訳です。

Single Scattering Microfacet BSDFのWhiteFurnaceTest、左からRoughnessが0~1エネルギー保存を満たすなら完全に背景に溶け込む

このエネルギー損失は実用的にも問題とされていて、粗い金属やすりガラスの表現をしようとするとぱっと見ただけでも分かるぐらいに暗くなる形で現れます。特にガラスなんかは分かりやすいですね。

[Heitz, et al. 2016]より引用　Single-Scatteringだけだとガラスが確かに暗くなる

そのため、Multiple-Scatteringを考慮したMicrofacet BSDFはなんとしても欲しいところでありますが、やはり完全にMultiple-Scatteringを含めたBSDFをClosed-Formに表すのは難しく、現在に至るまで様々な手法が考えられています。

今回の記事では[Heitz, et al. 2016]の手法とSIGGRAPH ASIA 2023で発表されたFast Multiple Scattering Approximationという手法の実装について書いていこうと思います。ついでに近年の手法とか簡単にまとめてみようと思います。

[Heitz, et al. 2016]のMultiple-Scattering

[Heitz, et al. 2016] Multiple-Scattering Microfacet BSDFs with the Smith ModelはMicrofacetを一種のボリュームと見なしてボリュームレンダリングによってMultiple-Scatteringの散乱を実行するという手法です。この手法はアンバイアスであり、現在ではMultiple-Scatteringの基本理論として様々な手法の基礎、リファレンスとして扱われる重要な手法です。

初めて聞いた方だと「霧とか雲とかじゃなくて固体っぽいものにボリュームレンダリング？」と思うかもしれませんが、実はボリュームレンダリングは髪や毛などの複雑な形状を持つ固体の複雑な散乱の表現に使われたりと一般な多重散乱を正確に表現できる手法であります。うまいこと適切に係数さえ設定できれば今までClosed-Formにできなかった難しい散乱がボリュームレンダリングとして計算することが出来ます。

どうやってその係数とかを設定するのかというのはレイトレ合宿9のセミナー資料で書きましたので詳しい理論的な話はそちらの方をご参照ください。Microfacetの3次元版みたいなMicroflake理論から出発して、Microflake-Volume間の橋渡しをしてMicrofacetをMicroflakeの一つと見なして係数を設定して行っています。

speakerdeck.com

Microflake 理論

Microflake理論とはVolumeの散乱について「ボリュームの散乱はボリューム内の微細な粒子(Microflake)の散乱によって引き起こされる」という仮定のもとボリュームの性質を物理ベースに導出していく理論です[Jakob, et al. 2010]。ボリュームには消散係数 $\mu_t$ 、散乱係数 $\mu_s$ 、位相関数 $f_p$ の3つのパラメーター(関数)が必要であり、これをミクロなMicroflakeの散乱や配置から導いていきます。

まずはMicroflakeが具体的にどんなものであるかというのを知る必要がありますね。Microflakeはある方向を向いた平面の極小の粒子であり、いわばMicrofacetをそのまま3次元に配置したようなものです。Microflakeはそれぞれバラバラな法線を向いても良く、ある法線 $\omega_m$ 方向を向いているMicroflakeの面積比を法線分布関数 $D(\omega_m)$ で表現します(Microfacetと同様の定義です)。また、Volume内でのMicroflakeがどのように分布しているかを単位面積当たりのmicroflake数で定義される密度 $\rho$ で表すこととします。

Microflakeに当たった光はMicroflake単体の位相関数 $f_{micro}$ (BSDFと思って大丈夫です)に従い、様々な方向へと反射されます。法線 $\omega_m$ のMicroflakeが $\omega_i$ 方向から入射してきた光を $\omega_o$ 方向へと出射させる時の位相関数を $f_{micro}(\omega_m,\omega_i,\omega_o)$ と表記しましょう。ここで注意ですが入射方向ベクトルはVolumeに入って行く方向を指しています。BSDFとかで考える入射方向ベクトルとは逆なことに注意してください。

[Jakob, et al. 2010]や[Heitz, et al. 2015]によるとMicroflakeの密度 $\rho$ 、法線分布関数 $D(\omega_m)$ 、位相関数 $f_{micro}(\omega_m,\omega_i,\omega_o)$ とMicroflakeによって作られるボリュームの消散係数 $\mu_t$ 、散乱係数 $\mu_s$ 、位相関数 $f_p$ には次のような関係があることが求められています。

$\displaystyle{ \mu_t(\omega_i) = \rho \int_{\Omega} \langle - \omega_i,\omega_m \rangle D(\omega_m) d\omega_m }$

$\displaystyle{ \mu_s(\omega_i) = \rho \int_{\Omega} \langle - \omega_i,\omega_m \rangle D(\omega_m) d\omega_m }$

$\displaystyle{ f_p(\omega_i,\omega_o) = \frac{1}{\mu_s(-\omega_i)} \int_{\Omega} \rho f_{micro}(\omega_m,\omega_i,\omega_o) \langle \omega_i, \omega_m \rangle D(\omega_m) d\omega_m }$

式がややこしいのでここでMicroflakeの $-\omega_i$ 方向への投影面積 $\sigma(\omega_i)$

$\displaystyle{ \sigma(\omega_i) = \int_{\Omega} \langle -\omega_i,\omega_m \rangle D(\omega_m) }$

を導入すると次のように書くことが出来ます。

$\displaystyle{ \mu_t(\omega_i) = \rho \sigma(\omega_i) }$

$\displaystyle{ \mu_s(\omega_i) = \rho \sigma(\omega_i) }$

$\displaystyle{ f_p(\omega_i,\omega_o) = \frac{1}{\sigma(- \omega_i)} \int_{\Omega} f_{micro}(\omega_m,\omega_i,\omega_o) \langle \omega_i,\omega_m \rangle D(\omega_m) d\omega_m }$

Microvolume

Microfacet理論ではMicrofacetの法線と高さが非相関的である、要するにランダムに平面が配置されたような形状をしているとするSmithモデルが良く用いられており、シャドウマスキング関数などはこの仮定から導出されています。

先ほど話していたMicroflake理論ではMicroflakeは小さな平面でボリューム内にバラバラっと配置されているという形でモデルを作っていたことを考えると、SmithモデルのMicrofacetというのはMicroflakeの一種として見なすことがというアイディアが浮かび上がります。

[Heitz, et al. 2016]はこのアイディアを基にSmithモデルのMicrofacetをMicroflakeと見なし、Microvolumeというボリュームに対応付けることでMultiple-Scatteringをボリュームレンダリングの範疇へと落とし込むことに成功しています。

Microflake理論ではMicroflakeに密度 $\rho$ 、法線分布関数 $D(\omega_m)$ 、位相関数 $f_{micro}(\omega_m,\omega_i,\omega_o)$ というパラメーターを設定する必要がありましたが、Smith-Modelにおけるこれらのパラメーターはどうなっているか考えていきましょう。

その前に軽く散乱におけるレイの記号について決めておきます。Microvolumeの散乱ではMultiple-Scatteringと同じ動きをするため、ある高さ $h_r$ にあったレイが $\omega_r$ 方向に進んで高さ $h_{r+1}$ に衝突して反射、これを繰り返していって散乱を再現します。

法線分布関数

まずは法線分布関数ですが、そもそもMicrofacetとMicroflakeで使う法線分布関数は同じように定義されているので、当然ながらMicrofacetで良く用いられているBeckmann分布やGGX(Trowbridge Reitz分布)を使って問題ありません。ここでは等方的GGXを使うことにしましょう。

$\displaystyle{ D(-\omega_m) = \frac{1}{\pi \alpha^2(\frac{m_x^2}{\alpha^2} + \frac{m_y^2}{\alpha^2} + m_z^2)^2} }$

これで投影面積を求められますが、中々ややこしい話があります。投影面積は $D$ は $\omega_m \cdot n$ が0以下のの領域では0となるため、y軸に対して対称性を持ちません。なので $\sigma(\omega_r)$ と $\sigma(-\omega_r)$ で積分の値が異なってくるという話があります。ということなのでそれぞれ求めていきます。

$\sigma(-\omega_r)$ に関しては、これは[Heitz 2014]のAppendix Aに同様の式の導出があり、それに従い解いていくと次のように求めることが知られています。

$\displaystyle{ \sigma(\omega_r) = (1 + \Lambda(\omega_r)) \cos{\theta_r} }$

ここで使用されている $\Lambda$ 関数はSingle-Scattering Microfacet BSDFのシャドウマスキング関数で出てくるあの $\Lambda$ 関数ですが、ここではちょっと拡張されてsign関数が使われています。ここは実装で注意しないといけないポイントです。

$\displaystyle{ \Lambda(\omega) = \frac{-1 + sign(a) \sqrt{1 + \frac{1}{a^2}}}{2}, a = \frac{1}{\alpha \tan{\theta}}. }$

じゃあ $\sigma(\omega_r)$ はどうなのかという話ですが、[Heitz, et al. 2016]では法線分布の性質を用いて導出を行います。方向ベクトル $\omega_r$ に対して法線分布関数 $D(\omega_m)$ は

$\displaystyle{ \int_{\Omega} D(\omega_m) (\omega_r \cdot \omega_m) d\omega_m = \cos{\theta_r} }$

という定義を持ちます。この時内積 $(\omega_m \cdot n)$ をクランプ内積 $\langle \omega_m \cdot n \rangle$ に分解してみれば、投影面積に対する条件を導くことが出来ます。

$\begin{align} \int_{\Omega} D(\omega_m) (\omega_r \cdot \omega_m) d\omega_m &= \int_{\Omega} D(\omega_m) \langle \omega_r \cdot \omega_m \rangle d\omega_m + \int_{\Omega} D(\omega_m) \langle \omega_r \cdot \omega_m \rangle d\omega_m \\\ &= \sigma(\omega_r) + \sigma(-\omega_r) = \cos{\theta_r} \\\ \end{align}$

従って、 $\sigma(\omega_r)$ は次のように与えられます。

$\displaystyle{ \sigma(\omega_r) = \Lambda(\omega_r) \cos{\theta_r} }$

また、これより可視法線分布 $D_{\omega_r}(\omega_m)$ は投影面積を使って、次にように定義することが出来ます。

$\displaystyle{ D_{\omega_r}(\omega_m) = \frac{\langle \omega_r,\omega_m \rangle D(\omega_m) }{\sigma^{smith}(-\omega_i)} }$

密度

ある地点におけるMicrofacetの単位面積当たりの個数として密度 $\rho$ を定義する必要がありますが、Single-ScatteringとかではこうしたMicrofacetの分布を決めたりすることはないのであまり馴染みがないと思います。

[Heitz, et al. 2016]ではMicrofacetの高さの分布を表す[tex: P¹(h)]という分布を導入して密度 $\rho$ を定義します。単純に[tex: P¹(h)]をそのまま密度 $\rho(h)$ と定義しても良さそうなものですが、ややこしいことにレイの出発地点の高さ $h$ に対して条件付確率として以下の様に定義する必要があります。

$\displaystyle{ \rho^{smith}(h') = P^1(h | h' \le h) = \frac{\chi^{+}(h' \le h) P^1(h')}{C^1(h)} }$

なぜこんな定義にする必要があるのかというと、Microfacet上のレイというのは常にMicrofacetの上に存在しているためです。SmithモデルではMicrofacetの重なりを許さないので図のような雪庇状の凹凸などは存在せず、レイの真上方向には一切Microfacetは存在しません。そのため、レイより高いMicrofacetというのは考える必要がありません。

なんだかそれだと高さだけ見たらレイより高いとこにあるMicrofacetとの衝突が考えられないんじゃないかと思うかもしれませんが、あくまでこの密度はある時の瞬間的に考えるもので、実際レイが進んでいく様子は逐次的にレイを進めていき、Microfacetとの衝突はレイとMicrofacetの高さの一致によって判断することが出来ます。そのため、レイの今の高さより上のMicrofacetは感がる必要がなく、従って密度もレイより高い場所のものはないものとして扱っても良いという訳です。

レイの高さ $h$ を与えられた時、 $h$ にあるようなMicrofacetの密度は次のように表現することができ、これはレイがMicrofacetと当たる確率として扱うことが出来ます。

$\displaystyle{ \rho^{smith}(h) = \frac{P^1(h)}{C^1(h)} }$

こうして密度の定義はできましたが、 $P^1$ の分布って具体的にどうゆうの使ったらいいのかという疑問を持ったかもしれません。後述しますが、 $P^1$ は実は割と何でもよくちゃんと分布の条件(全積分したら1になる)さえ満たしていれば一様分布でもガウス分布でも問題なく、統計的にどの分布でも結果が同じになることが知られています。なので簡単のため適当な[-1,1]で[tex:P¹(h) = 1/2]のような一様の分布を考えてもらって大丈夫です。

Microfacetの位相関数

Microfacetの位相関数 $f_{micro}(\omega_m,\omega_i,\omega_o)$ ですがこれは単にMicrofacetでのBSDFを考えてあげれば問題ありません。Microfacet理論ではMicrofacet上の反射は理想反射、屈折が行われることとしているのでIdeal ReflectionやIdeal Refractionを位相関数としてあげれば問題ありません。今回は屈折までやると大変なのでIdeal Refractionのみの金属面だけを考えて話を進めていきます。この時のMicrofacetの位相関数 $f_{micro}(\omega_m,\omega_i,\omega_o)$ はIdeal Reflectionとして次のように書くことが出来ます。

$\displaystyle{ f_{micro}(\omega_m,\omega_i,\omega_o) = \frac{F(\omega_i,\omega_m) \delta (\omega_m - \omega_h)}{4 |\omega_h \cdot \omega_i|} }$

ここで使われている $\omega_h$ は $\omega_i,\omega_o$ のハーフベクトルです。

$\displaystyle{ \omega_h = \bf{normalize}(\omega_i + \omega_o) }$

普通のIdeal Reflectionとはちょっと式が異なりますが、Microvolumeの位相関数がちゃんとエネルギー保存できるように調整したものなのであまり気にしないでください。そうゆうものだと考えて頂ければ幸いです。

Microvolumeのパラメーター

こうしてMicrofacetの各パラメータの設定が出来ました。次はこれをMicrovolumeのパラメータに変換してみましょう。

$\omega_i$ 方向を向いたレイが高さ $h$ にある時、その位置における消散係数 $\mu_t(\omega_i,t)$ は以下の様に書くことが出来ます。

$\displaystyle{ \mu_t(\omega_i,h) = \sigma(\omega_i) \rho(h) = \Lambda(\omega_i) \cos{\theta_i} \frac{P^1(h)}{C^1(h)}. }$

また、Microflakeの位相関数の定義にMicrofacetの位相関数を代入してあげると次のように書くことが出来ます。

$\begin{align} f_p(\omega_r,\omega_o) &= \frac{1}{\sigma(-\omega_i)} \int_{\Omega} \frac{F(\omega_i,\omega_m) \delta (\omega_m - \omega_h)}{4 |\omega_h \cdot \omega_i| }\langle \omega_i,\omega_m \rangle D(\omega_m) d\omega_m. \\\ &= \frac{F(\omega_r,\omega_h) D_{\omega_i}(\omega_h) }{4 |\omega_h \cdot \omega_r|} \end{align}$

以上でMicrovolumeに必要なパラメーターを定義することが出来ました。こうしてMicrovolumeのボリュームレンダリングを行えば晴れてMultiple-Scatteringを実行することが出来るのですが、[Heitz, et al. 2016]の手法は普通のと少々毛色が異なる形でボリュームレンダリングを行います。

流れをざっと説明します。まず最初に外から入射してきたレイ $\omega_1$ をHeight Fieldの最大高さ $h_1$ に設定し、次の衝突点の位置 $h_2$ をサンプリングします。この際の高さ $h_2$ のサンプリングはFree-Path サンプリングと呼びます。その後、衝突点で反射をMicrovolumeの位相関数から次の方向 $\omega_2$ をサンプリングして、再び衝突点の判定を繰り返していきます。最終的にはFree-Path サンプリングで得られた高さ $h_r$ が $P^1(h)$ の最大の高さを超えた場合、そこから先はMicrofacetが存在しないわけですから脱出した光と見なして出射方向ベクトル $\omega_o$ として返します。以上の流れはRandom Walkと言い、Multiple-Scatteringによって散乱した光を統計的に追跡していきます。

そんなわけで、Free-Pathサンプリングと位相関数の方向サンプリングについて求めていきましょう。

Free-Path Sampling

Free-Path サンプリングはある高さ $h_r$ から放たれた $\omega_r$ 方向のレイが高さ $h_{h+1}$ で衝突する確率 $P^1_{h_r,\omega_r}(h_{r+1})$ に基づいて行分ければいけませんが、その確率というのはどのようなものなのか考えていきましょう。

遮蔽されない確率

まずは、手始めにレイ $\omega_r$ が高さ $h_r$ から出発して距離 $\tau$ まで進んでも衝突せずに済む確率 $G_1(\omega_r,h_r,\tau)$ を導出しましょう。 $\tau$ からちょっと進んで $\tau + d \tau$ 進んだ際、消散係数の意味を考えればレイが持つ光は消散係数 $\mu_t$ の割合だけ減っていきます。従って確率 $G_1(\omega_r,h_r,\tau)$ は次のような線形微分方程式を満たし、

$\displaystyle{ \frac{G_1(\omega_r,h_r,\tau)}{d\tau} = - \mu_t(\omega_r,h_r) G_1(\omega_r,h_r,\tau) }$

線形微分方程式の解から $[0,\tau$ ]までの積分として次のように確率 $G_1(\omega_r,h_r,\tau)$ を求めることが出来ます。ボリュームレンダリングにおける透過率の積分ですね。ただし $\mu_t$ がレイの距離 $\tau$ ではなく高さ $h$ に依存するので変数変換して $\tau$ の積分として作っています。

$\displaystyle{ G_1(\omega_r,h_r,\tau) = \exp(-\int_0^{\tau} \mu_t(\omega_r, h_r + \tau' \cot \theta_r) |\frac{\partial h}{\partial \tau}| d\tau) }$

この式は実はシャドウマスキング関数の導出でも使われおり、次のように解くことが出来ます。Ushioさんの記事に同義の関数 $S(\xi_0,\mu,\tau)$ の証明が詳しく載っておりますのでご参照ください。

Microfacet入門（１） - ushiostarfish’s diary

$\displaystyle{ G_1(\omega_r,h_r,\tau) = (\frac{C^1(h_r)}{C^1 (h_r + \tau \cot \theta_r)})^{\Lambda(\omega_r)} = (\frac{C^1(h_r)}{C^1 (h_{r+1})})^{\Lambda(\omega_r)} }$

ここでは $\tau$ 進んだ先の高さを $h_{r+1}$ としています。こうして「 $\tau$ だけ進んでも衝突しない確率」というのを求めることが出来ました。これを用いてSingle Scatteringのように $h_r$ から出発した光がMicrofacetに当たることなくMicrofacetから脱出できる確率 $G_1(\omega_r,h_r)$ というのは次のように定義することが出来ます。

$\displaystyle{ G_1(\omega_r,h_r) = G_1(\omega_r,h_r,\infty) = (\frac{C^1(h_r)}{C^1 (\infty)})^{\Lambda(\omega_r)} }$

CDFの定義から $\tau \rightarrow \infty$ の場合は1になるので次のように書くこともできます。

$\displaystyle{ G_1(\omega_r,h_r) = C^1(h_r)^{\Lambda(\omega_r)} }$

衝突確率

以上で衝突しない確率というのは分かりました。これより、余事象である「ある $\tau$ までの距離を進んだ時に遮蔽されるレイの確率」を表す $C^1_{h_r,\omega_r}(h_{r+1})$ は $G_1(\omega_r,h_r,\tau)$ を用いて

$\displaystyle{ C^1_{h_r,\omega_r}(h_{r+1}) = 1 - G_1(\omega_r,h_r,\tau) }$

というように表されます。ちょっとややこしいですが、 $C^1_{h_r,\omega_r}(h_{r+1})$ はあくまで $\tau$ に至るまでに減っていたレイの合計の割合を示すものであることから、 $C^1_{h_r,\omega_r}(h_{r+1})$ と $P^1_{h_r,\omega_r}(h_{r+1})$ はCDFとPDFの関係にあります。このため、 $P^1_{h_r,\omega_r}(h_{r+1})$ は $C^1_{h_r,\omega_r}(h_{r+1})$ の微分から次のように導出されます(Appendx C参照)。中々ややこしいPDFではありますが重点的サンプリングで打ち消されるので実装上ではあまり気にしなくても大丈夫です。

$\displaystyle{\begin{align} P^1\_{h_r,\omega_r}(h\_{r+1}) &= \left\{ \begin{array}{cc} 0 & h_{r+1} < h_r, \theta_r < \frac{\pi}{2}\\ 0 & h_{r+1} > h_r, \theta_r > \frac{\pi}{2}\\ |\Lambda(\omega_r)| P^1(h_{r+1}) \frac{C^1(h_r)}{C^1(h_{r+1})^{1+\Lambda(\omega_r)}} + G_1(\omega_r,h_r) \delta_{\infty}(h_{r+1}) & \end{array} \right. \end{align} }$

Free-Path サンプリングはこのPDFに従って衝突位置の高さ $h_{r+1}$ を求めていく必要があり、先ほど出たCDF $C^1_{h_r,\omega_r}(h_{r+1})$ を逆関数法に掛けることで導くことが出来ます。[0,1]のUniform乱数 $U$ を考えた時、

$\displaystyle{ u = 1 - G_1(\omega_r,h_r,\tau) }$

という式から $h_{r+1}$ についての式を求めると以下のようなサンプリング関数を得ることが出来ます。

$\displaystyle{ h_{r+1} = C^{-1} (\frac{C^1 (h_r) }{(1- u)^{1/\Lambda(\omega_r)} }) }$

これがFree-Pathサンプリングですが、１つ重要な点があります。CDFには上限値として $1 - G_1(\omega_r,h_r)$ があり、これを超える場合それはMicrofacetに当たらなかったことを表します。つまり Microfacetの脱出は $u > 1 - G_1(\omega_r,h_r)$ によって判定することができ、実際のサンプリングは次の疑似コードのように $u$ の値でまず脱出の判定をしてサンプリングを行います。脱出の条件が満たされる場合には位相関数のサンプリングに移らずにBSDFの値を返す処理を行います。

[Heitz. et al. 2016]より引用 Free-Path サンプリングの疑似コード

位相関数のサンプリング

Free-Pathサンプリングによって衝突が行われた場合、そこで光がMicrofacetに当たり散乱として位相関数のサンプリングを行います。肝心の位相関数ですが、金属の場合は非常にシンプルに可視法線分布関数とフレネルによって

$\displaystyle{ f_p(\omega_r,\omega_o) = \frac{F(\omega_r,\omega_h) D_{\omega_i}(\omega_h) }{4 |\omega_h \cdot \omega_r|} }$

と表されるわけでした。従って位相関数のサンプリングには可視法線分布関数 $D_{\omega_r}$ のサンプリングを採用すればいいわけです。ここで可視法線サンプリングのPDFというのが出射ベクトル $\omega_{r+1}$ を用いると

$\displaystyle{ pdf(\omega_r,\omega_o) =\frac{D_{\omega_i}(\omega_h)}{4 |\omega_h \cdot \omega_r|} }$

という形になるため、位相関数の評価は $w_r$ としては次のように殆どPDFと打ち消されてモンテカルロ法でのウェイトとしては単にフレネルの値を掛けていけば問題ないわけです。

$\displaystyle{ weight = \frac{f_p(\omega_r,\omega_o)}{pdf(\omega_r,\omega_o)} = F(\omega_r,\omega_h) }$

Random walk

以上でランダムウォークに必要な式が出揃いました。ランダムウォークの流れは先ほど述べましたが、今一度式を踏まえて確認して行きましょう。

最初に入射してきた方向 $\omega_1$ 、輝度 $L$ のレイは初期位置として $P^1(h)$ の範囲で最大の高さ $h_1$ から来たものとして扱います。このレイが次に衝突する位置 $h_2$ を先ほどのFree-Pathサンプリングによって求めます。一様乱数 $u$ を用意して

$\displaystyle{ h_{2} = C^{-1} (\frac{C^1 (h_r) }{(1- u)^{1/\Lambda(\omega_r)} }) }$

という形でサンプリングを行います。ただし $u$ が $u > 1 - G_1(\omega_1,h_1)$ を満たす場合は脱出の処理を行います。ここでは一旦脱出しなかったとして話を進めていきます。

衝突点のサンプリングが終わったら次は位相関数から次の反射方向を求めます。位相関数のサンプリングは可視法線サンプリングであり、一様乱数 $v_1,v_2$ を用意してMicrofacetの法線 $m_1$ を取得します。

$\displaystyle{ m_2 \leftarrow \bf{sample}(D_{\omega_1},v_1,v_2) }$

Microfacetの法線 $m_1$ から次の方向ベクトル $\omega_2$ は通常の法線サンプリングと同様に $m_1$ を法線とする $-\omega_1$ の反射ベクトルとして計算します。

$\displaystyle{ \omega_2 = \bf{reflect} (-\omega_1, m_1) }$

この仮定における光の減少率としてウェイト $w_1$ を考え、そのウェイト $w_1$ はこの反射におけるフレネルの値が入れられます。

$\displaystyle{ w_2 = F(\omega_1,m_1) }$

$\omega_2$ のサンプリング後は同様に次の衝突点の高さ $h_3$ をFree-Pathサンプリングで求めて生き、脱出条件を満たすまで一連の流れを繰り返します。r番目におけるウェイト $w_r$ は前の反射におけるウェイトに今回のフレネルをかけていく形で更新されます。

$\displaystyle{ w_r = F(\omega_r,m_r) w_{r-1} }$

こうしてr番目の反射でようやく脱出条件 $u > 1 - G_1(\omega_1,h_1)$ を満たした時、ランダムウォークを終了させ、最終的に出射ベクトル $\omega_o$ を $\omega_r$ として、その輝度 $L_r$ を

$\displaystyle{ L_r = w_r L }$

という形でパストレーシング側に返してあげればMultiple Scatteringの処理は完了します。

BSDFを一般化して処理を行っている場合はスループットの更新をBSDFの値とコサイン項、PDFで割る処理になっていると思いますが、Random walkだとボリュームレンダリングの式で行っているのでちゃんと $w_r$ になるように返す値を設定する必要があります

$\displaystyle{ throughput *= \frac{BSDF \cos{theta_0}}{pdf} \rightarrow w_r }$

Implementation

実装に関しては非常にわかりやすい資料がimplementation documentとして提供されているのでimplementationのコードをCUDAに書き換えた私の実装例を見せながらパストレーシングで必要最低限の実装部分だけ抜き出して話していこうと思います。自分で実装する方はHetizの論文では座標系にZ-upを使用しており、私の実装はY-upのものであることに注意してください（私はこれで結構時間を食いました)。

eheitzresearch.wordpress.com

まずMicrofacetの高さ分布ですが、最も簡単な例として一様な[-1,1]のHeightFieldを考えたHeight PDF[tex: P¹(h)]を用います。

 __device__ float P1(const float h)const {
        //Uniform
        const float value = (h >= -1.0f && h <= 1.0f) ? 0.5f : 0.0f;
        return value;
    };

これに従い、Height CDF[tex: C¹(h)]は以下の様に実装ができ、同様にinverse CDFも簡単に実装することが出来ます。

 __device__ float C1(const float h)const {

        //Uniform
        const float value = fmin(1.0f, fmax(0.0f, 0.5f * (h + 1.0f)));
        return value;

    };

    __device__ float invC1(const float U)const {
        const float h = fmax(-1.0f, fmin(1.0f, 2.0f * U - 1.0f));
        return h;
    }

Free Pathサンプリングの実装の方はこんな感じです。概ね式の値をそのままやっている感じですが、Multiple-Scatteringだと入射ベクトルが様々な方向を取りやすいのでNaNチェックを入れておくような実装になっています。脱出したときの判定は返し値の高さがFLT_MAX(フロートの最大値)であることを確認する形で行います。

 __device__ float G_1_Height(const float3& wi, const float h0) const {
        if (wi.y > 0.9999f)
            return 1.0f;
        if (wi.y <= 0.0f)
            return 0.0f;
        // height CDF
        const float C1_h0 = C1(h0);
        // Lambda
        const float Lambda = GGX_Lambda(wi);
        // value
        const float value = powf(C1_h0, Lambda);
        return value;
    }

 __device__ float sampleHeight(const float3& wr, const float hr, const float U) const {

        if (wr.y > 0.9999f)
            return FLT_MAX;
        if (wr.y < -0.9999f)
        {
            const float value = invC1(U * C1(hr));
            return value;
        }
        if (fabsf(wr.y) < 0.0001f)
            return hr;

        // probability of intersection
        const float G_1_ = G_1_Height(wr, hr);
        if (U > 1.0f - G_1_) // leave the microsurface
            return FLT_MAX;
        const float h = invC1(
            C1(hr) / powf((1.0f - U), 1.0f / GGX_Lambda(wr))
        );
        return h;
    }

特に注意すべき点は $\Lambda$ 関数の実装です。通常の $\Lambda$ 関数と違いsign関数が含まれた式であるため、これを忘れるとなんかもったりした感じの結果が出てしまいます。その割には分かりにくいのでなんかバグなのかみたいなことがあったらこれを疑ってください。

$\displaystyle{ \Lambda(\omega) = \frac{-1 + sign(a) \sqrt{1 + \frac{1}{a^2}}}{2}, a = \frac{1}{\alpha \tan{\theta}}. }$

 __device__ float GGX_Lambda(const float3& wi) const
    {
        if (wi.y > 0.9999f)
            return 0.0f;
        if (wi.y < -0.9999f)
            return -1.0f;

        // a
        const float theta_i = acosf(wi.y);
        const float a = 1.0f / (tanf(theta_i) * alpha);

        // value
        const float value = 0.5f * (-1.0f + sign(a) * sqrtf(1 + 1 / (a * a)));
        return value;
    }

また、接空間座標表示の $\Lambda$ 関数であれば次のように実装できます。

$\displaystyle{ \Lambda(\omega) = \frac{-1 + sign(y_{\omega}) \sqrt{1 + \frac{\alpha^2 x_{\omega}^2 + \alpha^2 z_{\omega}^2}{y_{\omega}^2 }}}{2}. }$

 __device__ float GGX_Lambda(const float3& v) const {
        if (v.y > 0.9999f)
        return 0.0f;
        if (v.y < -0.9999f)
            return -1.0f;
        float delta = 1.0f + (alpha * alpha * v.x * v.x + alpha * alpha * v.z * v.z) / (v.y * v.y);
        return (-1.0 + sign(v.y) * sqrtf(delta)) / 2.0f;
    }

位相関数のサンプリングはこんな感じです。元実装だと古いVisible Normal Samplingが使われているため、今では[Heitz 2018]の手法や[Dupuy, Benyoub. 2023]の手法の方が無難だと思います。ここでは[Dupuy, Benyoub. 2023]の手法を使っています。

 __device__ float3 sampleVisibleNormal(float2 uv, float3 wo) {
        float3 strech_wo = normalize(make_float3(wo.x * alpha, wo.y, wo.z * alpha));

        float phi = 2.0f * PI * uv.x;
        float z = fma((1.0f - uv.y), (1.0f + strech_wo.y), -strech_wo.y);
        float sinTheta = sqrtf(clamp(1.0f - z * z, 0.0f, 1.0f));
        float x = sinTheta * cos(phi);
        float y = sinTheta * sin(phi);

        float3 c = make_float3(x, z, y);

        float3 h = c + strech_wo;

        float3 wm = normalize(make_float3(h.x * alpha, h.y, h.z * alpha));

        return wm;
    }

以上でサンプリング周りが出来たので、ランダムウォークの処理は前節の流れのように書いてあげれば問題ありません。元実装だとwhile文のbreakがMicrofacetの脱出時のみしかないせいか中々処理が終わらないことがあったので、ここでは5回以上の反射は追わないようにbreak分を追加しています。

__device__ float3 sample(const float3& wi, float3& wo, int& scatteringOrder, CMJState& state)
{
    // init
    float3 wr = -wi;
    float hr = 1.0f + invC1(0.999f);
    // random walk
    scatteringOrder = 0;
    float3 weight = make_float3(1.0);
    while (true)
    {
        // next height
        float U = cmj_1d(state);
        hr = sampleHeight(wr, hr, U);
        // leave the microsurface?
        if (hr == FLT_MAX) {
            //printf("test");
            break;
        }
        else
            scatteringOrder++;

        if (scatteringOrder > 5) {
            wo = make_float3(0, 0, 1);
            return make_float3(0, 0, 0);
        }
        float3 weight_1;
        wr = samplePhaseFunction(-wr, state, weight_1);
        weight *= weight_1;

        if ((hr != hr) || (wr.z != wr.z))
            return make_float3(0, 0, 1);
    }
    wo = wr;
    return weight;
}

#define MAX_SCATTERING_ORDER 5
    __device__ float3 sampleBSDF(const float3& wo, float3& wi, CMJState& state, float& pdf) {
        int scatteringOrder;
        float3 bsdf = sample(wo, wi, scatteringOrder, state);
        if (wi.y < 0.0 || scatteringOrder > MAX_SCATTERING_ORDER) {
            return make_float3(0.0);
        }
        pdf = fabsf(wi.y);
        return bsdf;
    }
};

これで5000sppで実行した結果がこんな感じです。左からRoughness=0で0.1ずつ増やしていったマテリアルを示しています(alphaにはroughenss * roughnessを使っています)。Single-Scatteringの実装と比べると確かにRoughness 0.5以降のマテリアルは明確に明るくなっていることが分かります。

Multiple-Scatteringの実装例 F0 = 0.8 (ガンマ補正あり)

ちゃんとエネルギー保存が成立しているかについても見てみましょう。フレネルを1として設定し、WhiteFurnaceTestにかけてみるとちゃんと全てのラフネスでエネルギー保存が成り立っていることが分かります(というかウェイトが1になるので当然と言えば当然な気がしますが)。フレネルを付ければちゃんとエネルギー減少が見られるのでエネルギー増加とかもなくしっかりPBRに基づいたモデルであると確認できます。

[Heitz et al. 2016]でのWhiteFurnaceTest 完全に背景と溶け込んでいるのでエネルギー保存が保たれている

フレネルを付けた場合のWhiteFurnaceTest　背景より明るくなることがないのでエネルギーが増加することはないことが確認できる

Fast Multiple Scattering Approximation

この記事を書くちょっと前ぐらいにSIGGRAPH ASIA 2023が開催されまして、pepar fast forwardを見ていたのですが、Multiple Scattering周りで面白そうな手法が発表されているのを発見しました。

https://dl.acm.org/doi/fullHtml/10.1145/3610548.3618231

Fast Multiple Scattering Approximation [Rosales et al. 2023]という手法なのですが、これはMultiple-Scatteringをうまいこと近似的にClosed-Formな式として提供するとのことで、リアルタイムレンダリングでも可能！とのことでしたので気になって実装してみました。この手法はMultiple-Scatteringにおける2nd bounceの光をV-Cavity構造の仮定から求め、Single-Scatteringに加算する2nd bounceの調節項を求めています。

V-CavityというのはSmith-Modelとは異なり、Microfacetが常にVの字のように対称的な形状を持っているとする有名な仮定の1つです。V-Cavity自体は昔から考えられており、Smithモデルが一般的になる前はV-Cavityから求められたシャドウマスキング関数が使われていましたが、Smithモデルの方が正確な結果が出るのであまり使われなくなりました。

ですが、解析的に計算する上ではV-Cavityという形状は便利であり、V-Cavityが非対称的であったり向きを変えられるといった自由度を高くした仮定を用いてMultiple-ScatteringのClosed-Formな式を導出した手法が[Lee et al. 2018] [Xie and Hanrahan 2018]にあったみたいです。ですが、反射に特異点があり、Roughnessが高いのに何故かつるつるした反射を引き起こすなど結構弱点があったと述べられています。画像を見ると確かになかなかな弱点とわかりますね。

[Rosales, et al. 2023]より引用、一番左が[Lee et al. 2018]の結果一番右が[Heitz et al. 2016]の結果

[Rosales et al. 2023]では[Lee et al. 2018]の計算を改善して、また2nd bounceにのみ注目して調節項の導出を行っているみたいです。具体的な計算式についてはちゃんと追えてませんが、2nd bounceの式として次のような式を使えばいいとのことです。

$\displaystyle{ f_{2nd}(\omega_i,\omega_o) = \frac{D_I(\theta_s) G_I(\omega_i \cdot \omega_o) F(\omega_i,\omega_h)^2}{2 (\omega_c \cdot \omega_v)} }$

まずベクトル $\omega_c$ ですがこれはハーフベクトル $\omega_h$ と法線 $n$ のさらなるハーフベクトルです。

$\displaystyle{ \omega_c = normalize(\omega_h + n) }$

この時 $\omega_c$ と $n$ が成す角度を $\theta_c$ を定義しています。それで $\theta_s$ という角度なのですがこれはV-Cavityの半角を表しており、次のような値として求められます(元論文だと $\theta_m$ なのですが記事内のMicrofacetの法線と紛らわしいので $\theta_s$ と書いています)。

$\displaystyle{ \theta_s = \frac{(\pi - acos(\omega_i \cdot \omega_o))}{4} }$

ここで使用される $D_I$ は単にGGXなのですが、対象となる変数が角度 $\theta_s$ なので実装上ではこっちの近似式を使った方がやりやすいと思います。

$\displaystyle{ D_I(\theta_s) = \frac{\alpha^2}{\pi (\cos^2\theta_m (\alpha^2 - 1) +1)^2} }$

シャドウマスキング関数 $G_I$ はSmithモデルとは全く別物で次のように定義されています。

$\displaystyle{ G_I(\omega_i,\omega_o) = \frac{\sin(\theta_{ic} - \theta_s)}{\sin(\theta_{ic} + \theta_s)} }$

$\theta_{ic}$ はベクトル $\omega_i$ と $\omega_c$ のなす角です。

こうして定義された2nd bounceの式をSingle-Scattering Microfacet BSDFの調節項として次のように加算した式がMultiple-ScatteringのBSDFとなります。

$\displaystyle{ f(\omega_i,\omega_o) = \frac{D(\omega_m) G(\omega_i,\omega_o) F(\omega_i,\omega_h) }{4 |\omega_i \cdot n| |\omega_o \cdot n| } + \frac{D_I(\theta_s) G_I(\omega_i \cdot \omega_o) F(\omega_i,\omega_h)^2}{2 (\omega_c \cdot \omega_v)} }$

こうやってみるとかなりシンプルでほんとに？という感じではありますが、実際に実装してみて[Heitz et al. 2016]で2nd bounceまで取った結果と比べてみるとかなり近い結果が出ていることが分かりました。さすがに5th bounceと比べるとちょっと暗い感じはありますが、それでもSingle-Scatteringより遥かに明るいですしなかなかいい感じの結果が出ています

[Rosales et al. 2023]、[Heitz et al. 2016]とSingle Scatteringの結果の比較(ガンマ補正あり)

実装は以下のような感じでやっています。重点的サンプリングについては可視法線サンプリングを使用しています。

 __device__ float multipleG(const float3& wo, const float3& wi, const float3& wc) {
        float theta_c = acosf(dot(wo, wc));
        float theta_m = (PI - acosf(dot(wo, wi))) * 0.25;
        float OP = sinf(theta_c - theta_m) / sinf(theta_c + theta_m);
        return 1.0 - fmaxf(0.0, OP);
    }

    __device__ float GGX_D_Approximate(const float& mdot) const {
        float term1 = (mdot * mdot * (alpha * alpha - 1.0) + 1.0);
        return alpha * alpha / (PI * term1 * term1);
    }

__device__ float3 sampleBSDF(const float3& wo, float3& wi, float& pdf, CMJState& state) {
    float2 xi = cmj_2d(state);
    const float3 wm = sampleVisibleNormal(xi, wo);
    //const float3 wm = sampleD(xi);

    wi = reflect(-wo, wm);

    if (wi.y <= 0.0) {
        pdf = 1.0f;
        return { 0.0,0.0,0.0 };
    }

    float ggxD = GGX_D(wm);
    float ggxG2 = GGX_G2_HeightCorrelated(wi, wo);
    float3 ggxF = shlickFresnel(F0, wi, wm);

    float jacobian = 0.25f / absdot(wo, wm);
    //Walter PDF
    //pdf = ggxD * wm.y * jacobian;

    //Visible Normal PDF
    float ggxG1 = GGX_G1(wo);
    pdf = ggxD * ggxG1 * absdot(wo, wm) * jacobian / fabsf(wo.y);

    float3 bsdf = ggxD * ggxG2 * ggxF / (4.0 * wo.y * wi.y);
    float3 wc = normalize(make_float3(0, 1, 0) + wm);
    float Gi = multipleG(wo, wi, wc);
    float theta_m = (PI - acosf(dot(wo, wi))) * 0.25;
    float cos_theta_m = cos(theta_m);
    float Di = GGX_D_Approximate(cos_theta_m);
    bsdf += Di * Gi * ggxF * ggxF / (2.0 * dot(wc, wo));

    return bsdf;
}

Closed Formな式というだけで色々ありがたいですし、その結果も結構[Heitz et al. 2007]の結果と近いので安価ながら結構嬉しい手法であると感じます。弱点とかはまだよくわからないので要検証ですが、個人的には[Heitz et al. 2016]の代わりとして使っても良さそうないい手法だと期待しています。

Multiple Scatteringの動向

今回の記事ではMultiple Scatteringの基礎となった[Heitz et al. 2016]と気になっていた[Rosales et al. 2023]について話しましたが、Multiple Scatteringの手法はまだまだ他にもあり盛んに研究されています。ここでは知っている限りの手法を簡単にまとめて行こうと思います。

Position Free

[Heitz et al. 2016]の手法は収束が遅いという問題があったのですが、これをMicrovolumeのPosition Freeという性質を用いてより分散を減少させる手法として[Wang et al. 2022] Position-free Multiple-bounce Computations for Smith Microfacet BSDFsが提案されています。レイトレ合宿8では水鳥さんがこの手法の実装をしていました。

github.com

MicrovolumeではNDFと高さが非相関なSmithモデルを採用しているおかげか、どうやら高さのサンプリングは必要ではないという性質があるようです。これをPosition Freeと言い、高さ分布 $P^1$ が変わっても結果が変わらなかったのはここから来ているらしいです。

そんな感じで改善された手法なのですが、導出が大変かつどうやらバイアスがあるみたいで今年のSIGGRAPH ASIAではMultiple-bounce Smith Microfacet BRDFs using the Invariance Principle [Cui et al. 2023]という論文で改善策が提案されたそうです。

エネルギー保存 BSDF

今までMultiple Scatteringを直接考えるのではなく、エネルギー保存の性質からSingle-Scattering Microfacet BSDFを無理やりエネルギー保存させる調節項を加えることでMultiple Scatteringを実現するというアプローチを取っている手法もあります。

このアプローチの先駆けとなったのはA Microfacet Based Coupled Specular-Matte BRDF Model with Importance Sampling[Kelemen 2001]の手法でエネルギー損失を拡散項でうまいこと補填するというような手法です。この手法はUshioさんが詳しく記事にまとめてくれています。

ushiostarfish.hatenablog.com

他にも最近ではRevisiting Physically Based Shading at Imageworks [Kulla and Conty. 2016]ではエネルギー損失を事前にLUTにして置き、実際のレンダリング時に補填するというような手法が取られています。ここで問題になるのがフレネルの値で、Multiple-Scatteringでは何度も反射する光も存在するので反射率となるフレネルの値も実際は複雑な形になります。何らかの代表的なフレネル値をどうにか近似しようと色々考えられているみたいです。この辺はPractical multiple scattering compensation for microfacet modelsという資料にシンプルにまとまっています。

こうしたMultiple-Scatteringの補填はごり押し的な感じはありますが、シンプルかつ強力なのでFillamentを始めとしたリアルタイム向けのマテリアルではこちらの手法が使われているらしいです。

https://google.github.io/filament/Filament.md.html#materialsystem/specularbrdf/geometricshadowing(specularg)

参考文献

[Heitz 2014] Eric Heitz. 2014. Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs

[Heitz et al. 2016] Eric Heitz, Johannes Hanika, Eugene d’Eon and Carsten Dachsbacher. 2016. Multiple-Scattering Microfacet BSDFs with the Smith Model

[Wang et al. 2022] Beibei Wang, Wenhua Jin, Jiahui Fan, Jian Yang, Nicolas Holzschuch, Ling-Qi Yan. 2022. Position-free Multiple-bounce Computations for Smith Microfacet BSDFs

[Wang et al. 2022] Eric Heitz, Jonathan Dupuy, Cyril Crassin, Carsten Dachsbacher. 2015. The SGGX Microflake Distribution

[Jakob et al. 2010] Wenzel Jakob, Adam Arbree, Jonathan T. Moon, Kavita Bala, Steve Marschner. 2010. A radiative transfer framework for rendering materials with anisotropic structure

[Cui et al. 2023] Yuang Cui, Gaole Pan, Jian Yang, Lei Zhang, Ling-Qi Yan, Beibei Wang. 2023 Multiple-bounce Smith Microfacet BRDFs using the Invariance Principle

[Rosales et al. 2023] Enrique Rosales, Fatemeh Teimury, Joshua Horacsek, Aria Salari, Xuebin Qin, Adi Bar-Lev, Xiaoqiang Zhe, Ligang Liu. 2023. Fast-MSX: Fast Multiple Scattering Approximation

[Kulla and Conty. 2016] Revisiting Physically Based Shading at Imageworks

[Kelemen. 2001] Csaba Kelemen and László Szirmay-Kalos. 2001. A Microfacet Based Coupled Specular-Matte BRDF Model with Importance Sampling

[Lee et al. 2018] Joo Ho Lee, Adrian Jarabo, Daniel S. Jeon, Diego Gutierrez, and Min H. Kim. 2018. Practical Multiple Scattering for Rough Surfaces.

[Xie and Hanrahan. 2018] Feng Xie and Pat Hanrahan. 2018. Multiple Scattering from Distributions of Specular V-Grooves.

[Heitz. 2018] Eric Heitz. 2018. Sampling the GGX Distribution of Visible Normals

[Dupuy and Benyoub. 2023] Jonathan Dupuy, Anis Benyoub. 2023. Sampling Visible GGX Normals with Spherical Caps

ushiostarfish’s diary Microfacet入門（1） Microfacet入門（１） - ushiostarfish’s diary

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