マイクロファセット法線分布関数を作ってみる

これは レイトレ Advent Calendar 2018 - Qiita の記事として作成しました .

はじめに

マイクロファセットモデルには, 法線分布関数というものがあります[Heitz2014].
"GGX"と唱えれば, いい感じにラフな見た目を作り出してくれるアレです.f:id:kenha:20181209173133j:plain

{ \begin{align} D_{\mathrm GGX}({\bf m}) = \frac{\alpha^2}{\pi\cos^4\theta_m(\alpha^2+\tan^2\theta_m)^2}\qquad[1/{\mathrm sr}] \tag{1} \end{align} }
この関数は, 入射光がどのくらいの範囲に散乱されるかを表す, 単位立体角当たりの(規格化された)面積分布を表しています.
マイクロファセット法線とハーフベクトルが等しいという仮定から, これは, 物体表面の微細構造による散乱の統計的な結果を, 鏡面で出来た3次元形状の面積分布に置き換えていると見てもいいのではないでしょうか.
この見かたが正しければ, 適当な形状の鏡を持ってきて, 単位立体角当たりの面積分布を作ることで, それに対応する法線分布関数が求まるはずです.
この機会にちょっと試してみます.

マイクロファセット分布関数を作ってみる

話を簡単にするために, 3次元形状は回転体に限るとしておきます.
法線 {\bf m} をもつ微小面 d\sigma_m [{\mathrm m^2}] に対して, 方向 {\bf i} から光のビームを打ち込むことを考えます.



ビームの広がりを d\sigma\,[{\mathrm m^2}], 放射照度を E_i\,[{\mathrm W}/{\mathrm m}^2] とすれば, 微小面が受ける放射束は

{ \begin{align*} d\Phi = E_i d\sigma = E_i \cos\theta_m d\sigma_m \qquad[{\mathrm W}] \end{align*} }


となります.
次に, 光の当たった微小面 {d\sigma_m} の位置を決めます.光は, 角度 {\theta_m} に応じて回転軸から距離 {b} だけずれた場所に当たるとし, さらに形状が {b} を使って {z(b)} と書けるとします. その上で, 形状の傾きはその位置の傾斜角に等しい

{ \begin{align} \frac{dz(b)}{db} = \tan\theta_m \tag{2} \end{align} }

という仮定を行えば, {z(b)} と {\theta_m} から {b} を求めることが可能になります. これにより, 微小面 {d\sigma_m} と距離 {b} のあいだに

{ \begin{align} d\sigma = \cos\theta_m d\sigma_m = bdbd\phi_m \tag{3} \end{align} }

という関係を作ることができます.
また, マイクロファセットの法線 {\bf m} がハーフベクトル {{\bf h}=\frac{{\bf i}+{\bf o}}{|{\bf i}+{\bf o}|}} と一致することから, 方向 {\bf o} に向かって反射した光の広がる微小立体角を {d\omega_m=\sin\theta_md\theta_md\phi_m} と表せるため,
{ \begin{align} bdbd\phi_m = b\left|\frac{db}{d\theta_m}\right|d\theta_m d\phi_m = \frac{b}{\sin\theta_m}\left|\frac{db}{d\theta_m}\right|d\omega_m \tag{4} \end{align} }
と書き直すことができます. ここで, {\frac{db}{d\theta_m}} は負の値を取ることがあるため, 絶対値を取るようにしました*1. 式(3)と(4)から, 単位立体角当たりの面積は
{ \begin{align} \frac{d\sigma_m}{d\omega_m} = \frac{1}{\cos\theta_m}\frac{b}{\sin\theta_m}\left|\frac{db}{d\theta_m}\right|\qquad[{\mathrm m^2}/{\mathrm sr}] \tag{5} \end{align} }
と得られます. この {\frac{d\sigma_m}{d\omega_m}} を法線分布関数 {D({\bf m})} に対応するものとみなして, その規格化について考えます. 法線分布関数の規格化条件は

{ \begin{align*} \int_\Omega D({\bf m})\cos\theta_m d\omega_m=1 \end{align*} }
でした[Heitz2014]. なお, {\Omega} は積分範囲の 単位半球面です. このことを考慮して,

{ \begin{align*} \sigma^{\mathrm(total)}=\int_\Omega\frac{d\sigma_m}{d\omega_m}\cos\theta_md\omega_m\quad[{\mathrm m^2}] \end{align*} }
なる全断面積 {\sigma^{\mathrm (total)}} を導入することにします.
この名前を"全断面積"としたのは, {d\sigma=\cos\theta_md\sigma_m} より, この量が入射光からマイクロファセットを眺めたときの断面積を表しているためです. これを規格化係数に用いることで, 法線分布関数を

{ \begin{align} D({\bf m}) = \frac{1}{\sigma^{\mathrm (total)}}\frac{d\sigma_m}{d\omega_m}\quad[1/{\mathrm sr}] \tag{6} \end{align} }
と定義することができます.

回転楕円体のマイクロファセット法線分布関数

以上の議論が正しければ, [TR1975] と同様にして, 回転楕円体からGGX分布が得られるはずです.
はじめに, 長半径を{1}, 短半径を{\alpha} とする回転楕円体


{ \begin{align*} x^2 + y^2 + \frac{z^2}{\alpha^2} = 1 \end{align*} }
を考えます. {x} と {y} が対称なことから, {b^2 = x^2 + y^2} としてこの形状を

{ \begin{align*} z(b) = \alpha\sqrt{1-b^2} \end{align*} }
と {b} の関数として書き直せば, 式(2) より

{ \begin{align*} \frac{dz(b)}{db} = -\frac{\alpha b}{\sqrt{1-b^2}} = \tan\theta_m \end{align*} }
となり, これを変形すれば {b} が

{ \begin{align*} b^2=\frac{\tan^2\theta_m}{\alpha^2 + \tan^2\theta_m} \end{align*} }
と 求まります. これを 式(5) に対して用いることで,

{ \begin{align*} \frac{d\sigma_m}{d\omega_m} &= \frac{1}{\cos\theta_m}\frac{b}{\sin\theta_m}\left|\frac{db}{d\theta_m}\right|\\ &= \frac{1}{\cos\theta_m\sin\theta_m}\frac{\tan\theta_m}{\sqrt{\alpha^2+\tan^2\theta_m}}\frac{\alpha^2}{\cos^2\theta_m(\alpha^2+\tan^2\theta_m)^{\frac{3}{2}}}\\ &= \frac{\alpha^2}{\cos^4\theta_m(\alpha^2+\tan^2\theta_m)^2} \end{align*} }
と得られます. ここで

{ \begin{align*} \frac{db}{d\theta_m} = \frac{\alpha^2}{\cos^2\theta_m(\alpha^2+\tan^2\theta_m)^{\frac{3}{2}}} \end{align*} }
を使いました. さらに, 全断面積も

{ \begin{align*} \sigma^{\mathrm(total)} &= \int_\Omega\frac{d\sigma_m}{d\omega_m}\cos\theta_m d\omega_m\\ &= \int^{2\pi}_0d\phi_m\int^{\frac{\pi}{2}}_0\frac{\alpha^2}{\cos^4\theta_m(\alpha^2+\tan^2\theta_m)^2}\cos\theta_m\sin\theta_md\theta_m\\ &= \pi \end{align*} }
と求まります. 以上の結果を 式(6) に代入すれば, 回転楕円体のマイクロファセット法線分布関数は

{ \begin{align*} D({\bf m}) = \frac{1}{\sigma^{\mathrm(total)}}\frac{d\sigma_m}{d\omega_m} = \frac{\alpha^2}{\pi\cos^4\theta_m(\alpha^2+\tan^2\theta_m)^2} \end{align*} }
と得られます. 確かに, これはGGX分布関数(1)になっています.

さいごに

この記事では, 鏡面を持つ回転体形状の, 単位立体角当たりの面積分布を考察することで, 法線分布関数を定義してみました. その一例として回転楕円体を考えたとき, GGX分布関数が得られることがわかりました.

微分散乱断面積との関連

この記事で行った光反射の考察はラザフォード散乱(https://en.wikipedia.org/wiki/Rutherford_scattering )を参考にしました. そして法線分布関数 {D({\bf m})} と, 微分散乱断面積{\frac{d\sigma}{d\omega_m}} は,


{ \begin{align*} D({\bf m})\cos\theta_m = \frac{1}{\sigma^{\mathrm (total)}}\frac{d\sigma}{d\omega_m} \end{align*} }
という関係にあるとしました.

参考文献

[Heitz2014] Eric Heitz. Understanding the masking-shadowing function in microfacet-based brdfs. Journal of Computer Graphics Techniques (JCGT), 3(2):48–107, June 2014.
[TR1975] T. S. Trowbridge and K. P. Reitz. Average irregularity representation of a rough surface for ray reflection. J. Opt. Soc. Am., 65(5):531–536, May 1975.

*1:いま, 面積を求める計算を行っているからです.