散乱理論とレンダリング方程式

December 22, 2020

これは レイトレ Advent Calendar 2020 の記事です。

はじめに #

散乱理論 [ 砂川 1977] では, 散乱現象を微分散乱断面積 \(\frac{d\bm{\sigma}}{d\bm{\omega}}\;[\text{m}^2\,\text{sr}^{-1}]\)

で表現します. この概念は, Rendering 方程式にある BRDF \(f_r(\bm{i},\bm{o})\;[\text{sr}^{-1}]\) とよく似ています.

微分散乱断面積は, レイリー散乱や, 剛体球に弾かれる粒子の散乱など, 様々な散乱現象を求めるのに使われます. これと BRDF との対応関係がわかれば, 散乱理論で培われた知見を, さまざまな BRDF の考察に応用することができるかもしれません.

この対応関係は, [ Arvo 1993] で示唆されているように, 既知であることが伺えます. しかし, このことを一貫して解説した文献を見つけることができなかったため, 自分なりにノートをまとめました.

散乱媒質と光 #

散乱媒質 #

空間を漂うホコリや, 酸素や窒素など大気を構成する分子, そして水滴などのエアロゾルを まとめて散乱媒質と呼ぶことにします. そしてこれらの粒子は, 空間中を位置 \(\bm{x}\), 運動量 \(\bm{p}\) で運動するものと考えます.

光 #

光とは, 電磁波の振動数 \(\nu\) が 可視光 の領域 400THz ～ 800THz 程度にあるものを指します. 電磁波にはその最小単位があり, それはエネルギー \(h\nu\) をもつ光子として知られています. 光子には, 空間を波動として伝搬し, 物質と粒子として相互作用する性質があります.

ここでは, 光子の波動的性質を無視します. つまり光子を, 空間を粒子として伝搬し, 物質と粒子として相互作用する物体とします.

この簡略化によって, 散乱媒質と光の両者は, 位置 \(\bm{x}\), 運動量 \(\bm{p}\) で空間を運動する 粒子として等しく扱うことができるようになります.

気体分子運動論 #

位置 \(\bm{x}\), 運動量 \(\bm{p}\) で空間を運動する粒子の集団があったとき, その集団の性質を考える体系は 気体分子運動論 として知られています. のちに, 散乱媒質による光の散乱を議論するため, ここで必要な知識をおさらいしておきます.

密度関数 #

時間 \(t\), 位置 \(\bm{x}(t)\), 運動量 \(\bm{p}(t)\) にある粒子の集団を, 密度関数 \(f(\bm{x},\bm{p},t)\) で表します¹.

そして, 密度関数を運動量体積要素 \(d\bm{p}=dp_xdp_ydp_z\) で積分すれば, 位置 \(\bm{x}\) の粒子数密度 \[ \frac{1}{h^3}\int f(\bm{x},\bm{p},t)d\bm{p} = n(\bm{x},t) \qquad[\text{m}^{-3}]\] が, 粒子数密度 \(n(\bm{x},t)\) を体積要素 \(d\bm{x}=dxdydz\) で積分すれば, 全粒子数 \[ \int n(\bm{x},t)d\bm{x} = N \qquad[\text{-}],\] と得られるとします. ここで, 全粒子数 \(N\) は粒子の運動に伴って変化しないとしました². また, \(d\bm{x}d\bm{p}\) の単位 \([\text{J}^3\,\text{s}^3]\) をキャンセルするため, プランク定数 \(h\,[\text{J}\,\text{s}]\) を使った係数 \(\frac{1}{h^3}\) を掛けました.

この密度関数を使い, 媒質の中を光がどう広がっていくかを考えます. 光が媒質の中を進むに従って衝突を繰り返し, その密度分布を変えていくことを表現するために, 密度関数の時間変化を \[ \frac{df}{dt} = \left(\frac{\partial f}{\partial t}\right)_\text{coll}\] と表します. 左辺は \(f(\bm{x}(t),\bm{p}(t),t)\) を \(t\) で全微分することで \[ \frac{df}{dt} = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial\bm{x}}\cdot\dot{\bm{x}} + \frac{\partial f}{\partial\bm{p}}\cdot\dot{\bm{p}}\] \[ \dot{\bm{x}}:=\frac{d\bm{x}}{dt},\; \dot{\bm{p}}:=\frac{d\bm{p}}{dt}\] と与えられます.

右辺は \(\left(\frac{\partial f}{\partial t}\right)_\text{coll}\) はこれから考えます. 以降, この項を衝突項と呼びます,

衝突率 #

粒子同士の衝突を考えるために, 単位時間当たりの衝突率 \(\eta\,[\text{s}^{-1}]\) を導入し, 衝突による密度関数の変化を \[ \left(\frac{\partial f}{\partial t}\right)_\text{coll} \sim f\,\eta\] と見積もることにします.

そして, 粒子へ空間的な広がりを持たせるために, 散乱断面積 \(\sigma\,[\text{m}^2]\) という概念を導入します.

これに空間中の粒子密度 \(n\,[\text{m}^{-3}]\) と, 粒子の相対速度 \(v\,[\text{m}\,\text{s}^{-1}]\) を掛ければ, 粒子同士の衝突率は \[ \eta = n \sigma v \qquad [\text{s}^{-1}]\] と求まります.

微分散乱断面積 #

次に, 位置 \(\bm{x}\) で \(\bm{p}\) と \(\bm{p}_*\) の粒子が衝突して, \(\bm{p}'\) と \(\bm{p}'_*\) になることを考えます.

このとき, 前節で導入した衝突率は方向依存性がないため, 球の形をした粒子同士の衝突しか表現することができません. そこで, 次のように微分散乱断面積 \(\left|\frac{d\bm{\sigma}}{d\bm{\omega}}\right|\) を導入します.

まず, 方向 \(\bm{i}\) と \(\bm{o}\) を \[ \bm{o}=\frac{\bm{p}}{|\bm{p}|},\quad \bm{i}=\frac{\bm{p}'}{|\bm{p}'|}\] とし, それらのなす角を \(\theta,\,\varphi\). 立体角を \(\bm{\omega}=\sin\theta d\theta d\varphi\) と書きます³.

そして, 方向 \(\bm{\omega}\) に垂直な微小断面積 \(d\bm{\omega}\) を, 微小立体角の関数であることを仮定して, \[ \sigma = \int d\bm{\sigma} = \int_{S^2}\left|\frac{d\bm{\sigma}}{d\bm{\omega}}\right|d\bm{\omega} \qquad[\text{m}^2\,\text{sr}^{-1}]\] と, 微分散乱断面積 \(\left|\frac{d\bm{\sigma}}{d\bm{\omega}}\right|\) を導入します.

衝突項 #

それでは, 衝突項の具体的な表現を与えます.

まず, 密度関数 \(f(\bm{x},\bm{p},t)\) が増えるときと減るときに分けて考え, 衝突項を獲得項と損失項の二つで表します. \[ \left(\frac{\partial f}{\partial t}\right)_\text{coll} = \left(\frac{\partial f}{\partial t}\right)_\text{gain} - \left(\frac{\partial f}{\partial t}\right)_\text{loss}\]

密度関数が減るときは, 位置 \(\bm{x}\) で \(\bm{p}\) と \(\bm{p}_*\) の粒子が衝突して, \(\bm{p}'\) と \(\bm{p}'_*\) になる場合です.

\(f(\bm{p})\) は, 運動量 \(\bm{p}_*\) を持つ粒子に衝突したぶんだけ減少します. その粒子数は \(f(\bm{p}_*)d\bm{p}_*\), 相対速度は \(|\bm{v}-\bm{v}_*|\), 散乱断面積は \(\left|\frac{d\bm{\sigma}}{d\bm{\omega}}\right|d\bm{\omega}\) となることから, 衝突率は \[ |\bm{v}-\bm{v}_*|\left|\frac{d\bm{\sigma}}{d\bm{\omega}}\right|f(\bm{p}_*)d\bm{p}_*d\bm{\omega} \qquad[\text{s}^{-1}]\] となります. 従って, 損失項は \[ \left(\frac{\partial f}{\partial t}\right)_\text{loss} = \iint|\bm{v}-\bm{v}_*|\left|\frac{d\bm{\sigma}}{d\bm{\omega}}\right|f(\bm{p})f(\bm{p}_*)d\bm{p}_*d\bm{\omega}\] と得られます.

密度関数が増えるときは, 減るときと逆に, \(\bm{p}'\) と \(\bm{p}'_*\) の粒子が衝突して \(\bm{p}\) と \(\bm{p}_*\) になる場合を考えます.

詳細は参考文献 [ Cercignani 1988; 藤田 1989; Kardar 2007] に譲って省略しますが, 衝突時のエネルギー保存と運動量保存を考慮することで, 損失項は \[ \left(\frac{\partial f}{\partial t}\right)_\text{gain} = \iint|\bm{v}-\bm{v}_*|\left|\frac{d\bm{\sigma}}{d\bm{\omega}}\right|f(\bm{p}')f(\bm{p}'_*)d\bm{p}_*d\bm{\omega}\] と得られます.

上記の損失項と獲得項を使い, 衝突項は \[ \left(\frac{\partial f}{\partial t}\right)_\text{loss} = \iint|\bm{v}-\bm{v}_*|\left|\frac{d\bm{\sigma}}{d\bm{\omega}}\right| \left[f(\bm{p}')f(\bm{p}'_*) - f(\bm{p})f(\bm{p}_*)\right]d\bm{p}_*d\bm{\omega}\] となります.

Boltzmann 方程式 #

以上の結果をまとめると \[\begin{aligned} &\left(\frac{\partial}{\partial t} + \dot{\bm{x}}\cdot\frac{\partial}{\partial\bm{x}} + \dot{\bm{p}}\cdot\frac{\partial}{\partial\bm{p}}\right)f(\bm{p})\\ &\qquad = \iint|\bm{v}-\bm{v}_*| \left|\frac{d\bm{\sigma}}{d\bm{\omega}}\right|\left[f(\bm{p}')f(\bm{p}'_*) - f(\bm{p})f(\bm{p}_*)\right]d\bm{p}_*d\bm{\omega} \end{aligned}\] が得られます. これは, Boltzmann 方程式として知られています. 詳しくは, [ Cercignani 1988; 藤田 1989; Kardar 2007] などを参考にしてください.

この節のまとめ #

この節では, 粒子の集団を 密度関数 \(f(\bm{x},\bm{p},t)\), 粒子同士の衝突を 微分散乱断面積 \(\left|\frac{d\bm{\sigma}}{d\bm{\omega}}\right|\) で書いたとき, その密度関数が空間を動く様子が Boltzmann 方程式 で与えられることを導きました.

Radiative Transfer #

光と媒質の Boltzmann 方程式 #

空間中に光と散乱媒質の2成分から成る粒子の集団があったとき, 光に関する Boltzmann 方程式を求めます.

光の密度関数を \(f_p(\bm{x},\bm{p},t)\), 媒質の密度関数を \(f_m(\bm{x},\bm{p},t)\) と書けば, Boltzmann 方程式は \[ \left( \frac{\partial }{\partial t} + \dot{\bm{x}}\cdot\frac{\partial}{\partial\bm{x}} + \dot{\bm{p}}\cdot\frac{\partial}{\partial\bm{p}} \right)f_p(\bm{p}) = \left(\frac{\partial f_p}{\partial t}\right)_{pm} + \left(\frac{\partial f_p}{\partial t}\right)_{pp}\] \[\begin{aligned} \left(\frac{\partial f_p}{\partial t}\right)_{pm} &= \int |\bm{v}-\bm{v}_*| \left|\frac{d\bm{\sigma}_{pm}}{d\bm{\omega}}\right|\left[f_p(\bm{p}')f_m(\bm{p}'_*) - f_p(\bm{p})f_m(\bm{p}_*)\right]d\bm{p}_*d\bm{\omega}\\ \left(\frac{\partial f_p}{\partial t}\right)_{pp} &= \int |\bm{v}-\bm{v}_*| \left|\frac{d\bm{\sigma}_{pp}}{d\bm{\omega}}\right|\left[f_p(\bm{p}')f_p(\bm{p}'_*) - f_p(\bm{p})f_p(\bm{p}_*)\right]d\bm{p}_*d\bm{\omega} \end{aligned}\] となります. \(\sigma_{pm}\) は光が媒質へ衝突するときの散乱断面積, \(\sigma_{pp}\) は光が光へ衝突するときの散乱断面積です. また, \(f(\bm{x},\bm{p},t)\) を \(f(\bm{p})\) と略記しています.

いま, 光に対して

• 外力を無視する \[ \dot{\bm{p}}=0\]
• 速度は \(c\,[\text{m}\,\text{s}^{-1}]\) で一定とする

と仮定します.

そして, 媒質に対して,

• 媒質の移動速度は, 光の速度 \(c\) に比べて無視できる \[ |\bm{v}-\bm{v}_*|\approx c\]
• 光との衝突による運動量変化は無視できる \[ f_m(\bm{p}'_*)\approx f_m(\bm{p}_*)\]
• 光同士の散乱断面積は, 媒質と光の散乱断面積に比べて無視できる \[ \sigma_{pm} \gg \sigma_{pp}\]

と仮定します.

これをさきほどの Boltzmann 方程式へ適用すると, \[ \left(\frac{\partial}{\partial t} + \dot{\bm{x}}\cdot\frac{\partial}{\partial\bm{x}} \right)f_p(\bm{p}) = c n_m \int \left|\frac{d\bm{\sigma}_{pm}}{d\bm{\omega}}\right|\left[f_p(\bm{p}') - f_p(\bm{p})\right]d\bm{\omega}\] となります. ここで, 媒質の粒子密度を \(n_m = \int f_m(\bm{p}_*)d\bm{p}_*\) と書きました.

いろいろと仮定を入れて方程式のかたちが変わってきたので, 先に進む前に, いったん記号を整理します.

放射輝度 #

微分散乱断面積 で導入した方向 \(\bm{i}\) と \(\bm{o}\) を使って, 運動量 \(\bm{p}\) と \(\bm{p}'\) を \[ \bm{p}=|\bm{p}|\bm{o},\quad \bm{p}'=|\bm{p}'|\bm{i}\] と書き直します. これに伴い, 密度関数 \(f_p(\bm{x},\bm{p},t)\) の引数は \(f_p(\bm{x},\bm{o},|\bm{p}|,t)\) となります.

以前に言及したように, いま光の粒子はエネルギー \(h\nu\) をもつ 光子であるとしています. そこで, Einstein - de Broglie の関係式 \(|\bm{p}|=\frac{h\nu}{c}\) を用いて, 密度関数 \(f_p(\bm{x},\bm{o},|\bm{p}|,t)\) の引数を \(f_p(\bm{x},\bm{o},\nu,t)\) と書き直します.

さらに, 放射輝度 \(L(\bm{x},\bm{\omega},\nu,t)\) とのあいだに \[ L(\bm{x},\bm{\omega},\nu,t) = \frac{2h\nu^3}{c^2}f_p(\bm{x},\bm{\omega},\nu,t) \quad[\text{W}\,\text{m}^{-2}\,\text{sr}^{-1}\,\text{Hz}^{-1}]\] という関係がある [ 藤田 1990] ことから, 以降は光の密度関数 \(f_p\) の代わりに放射輝度 \(L\) を用いることにします.

位相関数 #

微分散乱断面積を散乱断面積 \(\sigma\) で規格化したものを, 位相関数 \(p(\bm{i},\bm{o})\) と定義します. \[ p(\bm{i},\bm{o}) := \frac{1}{\sigma}\left|\frac{d\bm{\sigma}}{d\bm{\omega}_i}\right| \quad[\text{sr}^{-1}]\] 微分散乱断面積 の定義から, 位相関数は 1 に規格化されています. \[ \int p(\bm{i},\bm{o}) d\bm{\omega}_i = 1\]

散乱係数 #

粒子数密度 \(n\,[\text{m}^{-3}]\) と散乱断面積 \(\sigma\,[\text{m}^2]\) を使って, 散乱係数 \(\mu\) を \[ \mu:=n\sigma \qquad[\text{m}^{-1}]\] と定義しておきます.

光と媒質の Boltzmann 方程式 #

さいごに, 積分 \(\int d\bm{\omega}\) が, 方向 \(\bm{o}\) を軸として, 方向 \(\bm{i}\) をぐるっと回すことを表すことを強調するために, \(\bm{o}\) と \(\bm{i}\) のなす角を \(\theta_i,\,\varphi_i\), 微小立体角を \(d\bm{\omega}_i = \sin\theta_id\theta_id\varphi_i\) と書き直しておきます.

以上で導入した 放射輝度 と 位相関数 を使うと, さきほどの Boltzmann 方程式は, \[ \frac{\partial L(\bm{o})}{\partial t} + c\frac{\partial L(\bm{o})}{\partial\bm{x}}\cdot\bm{o} = c n_m \sigma_{pm} \int p(\bm{i},\bm{o}) L(\bm{i}) d\bm{\omega}_i - c n_m \sigma_{pm} L(\bm{o})\] となります. 両辺を整理し, 散乱係数 を使うと \[ \left(\frac{1}{c}\frac{\partial}{\partial t} + \bm{o}\cdot\frac{\partial}{\partial\bm{x}} + \mu_s\right)L(\bm{o}) = \mu_s\int p(\bm{i},\bm{o}) L(\bm{i}) d\bm{\omega}_i\] となります. ここで, 放射輝度 \(L(\bm{x},\bm{\omega}, \nu, t)\) を \(L(\bm{\omega})\) と略記しました.

これで, 媒質へ衝突しながら進む, 光の密度分布を表現する方程式が得られました.

Radiative Transfer #

局所熱平衡 #

Chandrasekhar は, 恒星大気内で局所的に熱平衡にある領域に対して, 媒質が吸収した光 \(\mu_a(\bm{x},\nu) L(\bm{x},\bm{o},\nu,t)\) は黒体輻射 \(B(\nu,T)\) で放出されるとしました [ Chandrasekhar 1960]. \[ \int \mu_a(\bm{x},\nu) L(\bm{x},\bm{o},\nu,t) d\nu = \int \mu_a(\bm{x},\nu) B(\nu, T) d\nu\] ここで, 媒質が光を吸収することを表現するために, 散乱係数 \(\mu_s\) に加えて吸収係数 \(\mu_a\) を導入しました.

Radiative Transfer 方程式 #

局所熱平衡 を考慮して, 媒質による光の吸収 \(\mu_aL\) と, 黒体輻射による放出 \(\mu_aB\) を先ほどの Boltzmann 方程式へ 導入すると, \[ \left(\frac{1}{c}\frac{\partial}{\partial t} + \bm{o}\cdot\frac{\partial}{\partial\bm{x}} + \mu_t\right)L(\bm{o}) = \mu_aB + \mu_s\int p(\bm{i},\bm{o}) L(\bm{i}) d\bm{\omega}_i\] となります. ここで全消散係数 \(\mu_t = \mu_a + \mu_s\) を導入しました.

この式は Radiative Transfer 方程式と呼ばれています [ Chandrasekhar 1960; Arvo 1993].

この節のまとめ #

この節では, 光と散乱媒質の Boltzmann 方程式 へ, 光が媒質に比べてとても小さく速いという仮定とともに, 局所熱平衡 を取り入れることで, Radiative Transfer 方程式 を導きました.

Rendering 方程式 #

Volume Rendering 方程式 #

Radiative Transfer方程式 に対して, 人間の感覚に合わせた近似や調整を行います.

まず, 人間にとって光は一瞬で伝わるものと考え, 放射輝度の時間依存性を落として \(L_\text{re}(\bm{x},\bm{\omega},\nu,t)\rightarrow L_\text{re}(\bm{x},\bm{\omega},\nu)\) とします. そして, 媒質が吸収した光の再放射は黒体輻射に限らないとして, 黒体輻射 \(B(\nu,T)\) を一般の再放射 \(L_\text{re}(\bm{x},\bm{\omega},\nu)\) で書き換えます.

これらを Radiative Transfer 方程式へ適用することで \[ \left(\bm{o}\cdot\frac{\partial}{\partial\bm{x}} + \mu_t\right)L(\bm{o}) = \mu_aL_\text{re}(\bm{o}) + \mu_s\int_{S^2} p(\bm{i},\bm{o}) L(\bm{i}) d\bm{\omega}_i\] が得られます. ここで, 放射輝度 \(L(\bm{x},\bm{\omega},\nu)\) を \(L(\bm{\omega})\) とするなど, 方向成分以外の引数を略して書きました. この記事では, この式を Volume Rendering 方程式と呼びます.

Rendering 方程式 #

Volume Rendering 方程式は, 位置 \(\bm{x}\) に関する微分方程式なので, これを解くには空間的な境界条件が必要になります.

そこで境界条件として, Volume Rendering 方程式から \(\bm{o}\cdot\frac{\partial L}{\partial\bm{x}}\) を落とすことで, 衝突率 100% の表面を用意します. \[ \mu_tL(\bm{o}) = \mu_aL_\text{re}(\bm{o}) + \mu_s\int_{H^+} f_r(\bm{i},\bm{o})(\bm{i}\cdot\bm{n}) L(\bm{i}) d\bm{\omega}_i\] なお, 表面の法線を \(\bm{n}\) とし, \(\bm{i}\cdot\bm{n}>0\) である半球面を \(H^+\) と書きました. また位相関数 \(p(\bm{i},\bm{o})\) の代わりに導入した \(f_r(\bm{i},\bm{o})(\bm{i}\cdot\bm{n})\) については, すぐ後で説明します.

ここで反射率 \(k_r:=\frac{\mu_s}{\mu_t}\), 吸収率 \(k_a:=\frac{\mu_a}{\mu_t}\) という量を導入すれば, \[ L(\bm{o}) = k_aL_\text{re}(\bm{o}) + k_r\int_{H^+} f_r(\bm{i},\bm{o})(\bm{i}\cdot\bm{n}) L(\bm{i}) d\bm{\omega}_i\] となります. なお, 散乱係数のあいだの関係 \(\mu_t=\mu_a+\mu_s\) から, 反射率と吸収率のあいだには \(k_r + k_s = 1\) という関係があります.

上式の \(k_aL_\text{re}(\bm{o})\) を自己発光 \(L_\epsilon\) に差し替えると, Kajiya の Rendering 方程式 [ Kajiya 1986] と同じものが得られることから, この記事では上式を Rendering 方程式と呼びます.

Bidirectional Reflectance Distribution Function (BRDF) #

法線 \(\bm{n}\) を持つ表面の微小面積を \(d\bm{\sigma}_n\) としたとき, それを方向 \(\bm{i}\) へ射影した面積は \((\bm{i}\cdot\bm{n})d\bm{\sigma}_n\) です. これを表面の上半球面 \(H_n^+\) で集めたものを, 散乱断面積 \(\sigma_r\) とします. \[ \sigma_r = \int_{H_n^+}(\bm{i}\cdot\bm{n})d\bm{\sigma}_n \qquad[\text{m}^2]\]

そして, 微分散乱断面積 で行った議論と同様にして, 単位立体角当たりの面積分布 \(\left|\frac{d\bm{\sigma}_n}{d\bm{\omega}_i}\right|\) を導入します. \[ \sigma_r = \int_{H_n^+}\left|\frac{d\bm{\sigma}_n}{d\bm{\omega}_i}\right|(\bm{i}\cdot\bm{n})d\bm{\omega}_i \qquad[\text{m}^2]\]

このとき, 面積分布 \(\left|\frac{d\bm{\sigma}_n}{d\bm{\omega}_i}\right|\) を散乱断面積 \(\sigma_r\) で規格化した量 \(f_r(\bm{i},\bm{o})\) を BRDF と呼びます. \[ f_r(\bm{i},\bm{o}):=\frac{1}{\sigma_r}\left|\frac{d\bm{\sigma}_r}{d\bm{\omega}_i}\right| \qquad[\text{sr}^{-1}]\]

なお, 方向 \(\bm{o}\) に対して法線 \(\bm{n}\) が傾いている場合, 散乱断面積は \((\bm{o}\cdot\bm{n})\sigma_r\) となります.

Microfacet モデル #

この節では, 粗い面を表現する BRDF である Microfacet BRDF を作ります. 作り方は [ Walter et al. 2007; Heitz 2014] と同様ですが, 微分散乱断面積に焦点を当てて話を進めます.

これまで考えていた, 法線 \(\bm{n}\) をもつ表面を macrosurface と呼び, その上の微細なデコボコを microsurface と呼びます そして, microsurface は完全鏡面で出来ているとします.

ハーフベクトル #

方向 \(\bm{i}\) と方向 \(\bm{o}\) を2分割する方向を, ハーフベクトル \(\bm{h}\) と呼びます. \[ \bm{h}:=\frac{\bm{i}+\bm{o}}{|\bm{i}+\bm{o}|}\] 完全鏡面の上では, その法線 \(\bm{m}\) とハーフベクトル \(\bm{h}\) は一致します.

そのため以後は, microsurface 上の法線を \(\bm{m}\) の代わりに \(\bm{h}\) と表記します.

法線分布関数 #

BRDF と同様にして, microsurface 上の微小面積 \(d\bm{\sigma}_h\) を考えます. これを方向 \(\bm{i}\) へ射影した面積 \((\bm{i}\cdot\bm{h})d\bm{\sigma}_h\) を microsurface 上の半球面 \(H_h^+\) で集めたものを, 散乱断面積 \(\sigma_r\) とします.

続いて, ハーフベクトルを使った微小立体角を \(d\bm{\omega}_h:=\sin\theta_hd\theta_hd\varphi_h\) と書き, 単位立体角当たりの面積分布 \(\left|\frac{d\bm{\sigma}_h}{d\bm{\omega}_h}\right|\) を \[ \sigma_r = \int_{H_h^+}\left|\frac{d\bm{\sigma}_h}{d\bm{\omega}_h}\right|(\bm{i}\cdot\bm{h})d\bm{\omega}_h \qquad[\text{m}^2]\] と導入します.

このとき, 面積分布 \(\left|\frac{d\bm{\sigma}_h}{d\bm{\omega}_h}\right|\) を散乱断面積 \(\sigma_r\) で規格化した量 \(D(\bm{h})\) を法線分布関数と呼びます. \[ D(\bm{h}):=\frac{1}{\sigma_r}\left|\frac{d\bm{\sigma}_h}{d\bm{\omega}_h}\right| \qquad[\text{sr}^{-1}]\]

ところで, BRDF では, 方向 \(\bm{o}\) に対して macrosurface が傾いている場合に, 散乱断面積が \((\bm{o}\cdot\bm{n})\sigma_r\) となるとしました.

Masking 関数 #

Microsurface の可視面積を議論するため, \(\bm{i}\) か \(\bm{o}\) のいずれかを表す方向を \(\bm{\omega}\) と書きます.

方向 \(\bm{\omega}\) へ射影した微小面積 \((\bm{i}\cdot\bm{h})d\bm{\sigma}_h\) を 集めるとき, Microsurfaceの表側がどれくらいあるかが記述された関数 \(G_1\) があったとして, \[\begin{aligned} (\bm{\omega}\cdot\bm{n})\sigma_r &= \int_{H_h^+}G_1(\bm{\omega},\bm{h}) (\bm{\omega}\cdot\bm{h})d\bm{\sigma}_h\\ &= \int_{H_h^+}G_1(\bm{\omega},\bm{h}) \left|\frac{d\bm{\sigma}_h}{d\bm{\omega}_h}\right|(\bm{\omega}\cdot\bm{h})d\bm{\omega}_h \end{aligned}\] と求まるとします. この関数 \(G_1\) は, Smith masking 関数として知られています. 詳細は [ Smith 1967; Walter et al. 2007; Heitz 2014] を参照してください⁴.

Shadowing-Masking 関数 #

隣接する microsurface のによる遮蔽は, 方向 \(\bm{i}\) と \(\bm{o}\) の両者に対して考える必要があります. ここでも詳細は [ Heitz 2014] に任せて, そのような遮蔽が考慮された関数 \(G_2\) があるとします.

この場合, 方向 \(\bm{i}\) から入射した光が, microsurface 上の 1 回反射で方向 \(\bm{o}\) へ向かわない領域が生じます.

そこで, 両者の遮蔽を考慮した場合の可視部分を \(\sigma_r^\text{single}\), それ以外を \(\sigma_r^\text{multi}\) と呼び分けておきます. \[ \sigma_r = \sigma_r^\text{single} + \sigma_r^\text{multi}\]

そうすると, \(\sigma_r^\text{single}\) は Shadowing-masking 関数 \(G_2\) を使って \[ (\bm{o}\cdot\bm{n})\sigma_r^\text{single} = \int_{H_h^+}G_2(\bm{i},\bm{o},\bm{h}) \left|\frac{d\bm{\sigma}_h}{d\bm{\omega}_h}\right|(\bm{i}\cdot\bm{h})d\bm{\omega}_h\] と書けます.

Fresnel の式 #

光が屈折率の異なる媒質同士の境界を通過するとき, 反射と屈折が起きます.

反射と屈折がどのように起こるかは Fresnel の式 によって記述されます. その詳細や 近似式 は参考文献 [ Born and Wolf 1999] [ Schlick 1994] を参照してください.

この記事では, 入射光が反射光と屈折光に分かれることを, 全断面積 \(\sigma_t\), 反射断面積 \(\sigma_r\), 吸収断面積 \(\sigma_a\) を使って \[ \sigma_t = \sigma_r + \sigma_a \qquad[\text{m}^2]\] と使って表すことにします.

そして, その反射成分の割合を表す Fresnel 反射率を \(F(\bm{i},\bm{h})\) と書きます.

Microfacet BRDF #

以上の結果を使って, Microfacet BRDF を導きます.

Shadowing-Masking 関数と Fresnel の式 で, 入射光が単一反射と屈折, そして多重反射に分かれることを \[ \sigma_t = \sigma_r^{\text{single}} + \sigma_r^{\text{multi}} + \sigma_a\] と表現しました. そして, 単一反射は \[ \sigma_r^{\text{single}} = \int_{H^+_h}\frac{F(\bm{i},\bm{h})G_2(\bm{i},\bm{o},\bm{h})}{(\bm{o}\cdot\bm{n})} \left|\frac{d\bm{\sigma}_h}{d\bm{\omega}_h}\right|(\bm{i}\cdot\bm{h})d\bm{\omega}_h\] と求まります.

\(d\bm{\omega}_h=\frac{1}{4(\bm{i}\cdot\bm{h})}d\bm{\omega}_i\) を使って, \(d\bm{\omega}_h\) を \(d\bm{\omega}_i\) へ変数変換すれば \[ \sigma_r^\text{single} = \int_{H^+_n}\frac{F(\bm{i},\bm{h})G_2(\bm{i},\bm{o},\bm{h})}{4(\bm{o}\cdot\bm{n})(\bm{i}\cdot\bm{n})} \left|\frac{d\bm{\sigma}_h}{d\bm{\omega}_h}\right|(\bm{i}\cdot\bm{n})d\bm{\omega}_i\] となり, これを BRDF の微分散乱断面積 \[ \sigma_r = \int_{H_n^+}\left|\frac{d\bm{\sigma}_n}{d\bm{\omega}_i}\right|(\bm{i}\cdot\bm{n})d\bm{\omega}_i\] と見比べると, 単一反射について \[ \left|\frac{d\bm{\sigma}_n}{d\bm{\omega}_i}\right| = \frac{F(\bm{i},\bm{h})G_2(\bm{i},\bm{o},\bm{h})}{4(\bm{o}\cdot\bm{n})(\bm{i}\cdot\bm{n})} \left|\frac{d\bm{\sigma}_h}{d\bm{\omega}_h}\right|\] という関係が得られます.

これにより, 単一反射を表す Microfacet モデルの BRDF は [ Walter et al. 2007; Heitz 2014] と同様に \[ f_r(\bm{i},\bm{o}) = \frac{F(\bm{i},\bm{h})G_2(\bm{i},\bm{o},\bm{h})D(\bm{h})}{4(\bm{o}\cdot\bm{n})(\bm{i}\cdot\bm{n})}\] と求まりました.

この節のまとめ #

[ Arvo 1993] と同様に, Radiative Transfer 方程式 で時間変化を無視できるとして Volume Rendering 方程式 を導き, その境界条件として Rendering 方程式 を考えました.

そして [ Walter et al. 2007; Heitz 2014] と同様にして Microfacet BRDF を導きました.

GGX 法線分布関数 #

前節で Microfacet BRDF を導きました. BRDF は微分散乱断面積で書いたので, 散乱理論の手法を応用できます.

この節では, その一例として 衝突パラメータ を用いた GGX 法線分布関数 [ Trowbridge and Reitz. 1975; Walter et al. 2007] の導出を行います.

衝突パラメータ #

はじめに, 衝突パラメータを用いた微分散乱断面積の表現方法をおさらいします [ Landau and Lifshitz 1976; 砂川 1977].

いま, microfacet が回転対称な形状をしているとして, その中心軸からの距離を衝突パラメータ \(b\) と呼びます.

このとき, 微小散乱断面積 \(d\bm{\sigma}_r\) と \(b\) の関係は. \[ d\bm{\sigma}_r = bdbd\varphi_h\] となります.

ここで衝突パラメータ \(b\) が \(\theta_h\) の関数であるとすれば, \[ d\bm{\sigma}_r = \left|\frac{db}{d\sigma_h}\right|d\theta_hd\varphi_h = \frac{1}{\sin\theta_h}\left|\frac{db}{d\theta_h}\right|d\bm{\omega}_h\] と書き換えることができます. 最後の等号で, \(d\bm{\omega}_h=\sin\theta_hd\theta_hd\varphi_h\) を使いました.

従って, 衝突パラメータ \(b\) を使って微分散乱断面積を \[ \left|\frac{d\bm{\sigma}_r}{d\bm{\omega}_h}\right| = \frac{1}{\sin\theta_h}\left|\frac{db}{d\theta_h}\right|\] と書くことができます.

衝突パラメータによる法線分布関数の表現 #

前節の 衝突パラメータ を使って, 法線分布関数 を表現します.

微小散乱断面積 \(d\bm{\sigma}_r\) と microfacet 上の微小面積 \(d\bm{\sigma}_h\) のあいだの関係, \[ d\bm{\sigma}_r = \cos\theta_h d\bm{\sigma}_h\] を, 前節で求めた微分散乱断面積へ用いると, 単位立体角当たりの面積分布が \[ \left|\frac{d\bm{\sigma}_h}{d\bm{\omega}_h}\right| = \frac{1}{\cos\theta_h}\frac{1}{\sin\theta_h}\left|\frac{db}{d\theta_h}\right|\] と求まります.

法線分布関数の定義は, \[ D(\bm{h}):=\frac{1}{\sigma_r}\left|\frac{d\bm{\sigma}_h}{d\bm{\omega}_h}\right| \qquad[\text{sr}^{-1}]\] としたので, 次は衝突パラメータ \(b\) で具体的な形状を表現することが必要です.

GGX 法線分布関数の導出 #

いま, microfacet の形状が \(z(b)\) なる関数で書けて, ハーフベクトルと \[ \frac{dz}{db} = \tan\theta_h\] という関係にあるとします[ Trowbridge and Reitz. 1975].

回転楕円体 #

\(z(b)\) が表す形状として, [ Trowbridge and Reitz. 1975] と同様に, 回転楕円体 を考えます.

長半径が \(1\), 短半径が \(\alpha\) の回転楕円体は \[ x^2 + y^2 + \frac{z^2}{\alpha^2} = 1\] です.

この式を \[ z = \alpha\sqrt{1-x^2-y^2}\] と変形し, 衝突パラメータ \(b\) を使って \[ b^2 = x^2 + y^2\] とすれば, \[ z(b) = \alpha\sqrt{1-b^2}\] と書き直せます.

GGXモデルの導出 #

\(z(b)\) と \(\theta_h\) の関係 \(\frac{dz}{db} = \tan\theta_h\) を 回転楕円体へ適用すれば, 衝突パラメータ \(b\) が \(\theta_h\) の関数として \[ b^2 = \frac{\tan^2\theta_h}{\alpha^2 + \tan\theta_h}\] と求まります.

これを使い, 単位立体角当たりの面積分布を求めると \[ \left|\frac{d\bm{\sigma}_h}{d\bm{\omega}_h}\right| = \frac{1}{\cos\theta_h}\frac{1}{\sin\theta_h}\left|\frac{db}{d\theta_h}\right| = \frac{\alpha^2}{\cos^4\theta_h(\alpha+\tan^2\theta_h)^2}\] となります.

従って散乱断面積は \[\begin{aligned} \sigma_r &= \int_{H_h^+}\left|\frac{d\bm{\sigma}_h}{d\bm{\omega}_h}\right|\cos\theta_hd\bm{\omega}_h\\ &= \int^{2\pi}_0d\varphi_h\int^{\frac{\pi}{2}}_0 \frac{\alpha^2}{\cos^4\theta_h(\alpha^2+\tan^2\theta_h)^2}\cos\theta_h\sin\theta_hd\theta_h\\ &= \pi \end{aligned}\] となり, 法線分布関数は \[ D(\bm{h}) = \frac{1}{\sigma_r}\left|\frac{d\bm{\sigma}_h}{d\bm{\omega}_h}\right| = \frac{\alpha^2}{\pi\cos^4\theta_h(\alpha+\tan^2\theta_h)^2}\] と求まります.

以上で GGX 法線分布関数 [ Trowbridge and Reitz. 1975; Walter et al. 2007] が求まりました.

この節のまとめ #

この節では, はじめに 法線分布関数 と 衝突パラメータ の関係を作りました.

そして, 回転楕円体 を衝突パラメータで表現することで, GGX 法線分布関数 を求めました.

まとめ #

この記事では, 主に [ Cercignani 1988] [ Arvo 1993] を参考にして Boltzmann 方程式 から Radiative Transfer 方程式 を導き, そこから Volume Rendering 方程式 と Rendering 方程式 を導きました.

次に [ Heitz 2014] を参考にして Microfacet BRDF を導きました. この際, 光と媒質の衝突部分を一貫して 微分散乱断面積 を用いて記述することで, 何らかの幾何形状から Microfacet BRDF を求められるようにしました.

最後に, Microfacet BRDF を求める例のひとつとして, 散乱理論 [ 砂川 1977] で知られている 衝突パラメータ を使って, 回転楕円体 から GGX 法線分布関数 を求める例を紹介しました.

参考文献 #

Lorentz, H. A. 1905. The motion of electrons in metallic bodies i. In KNAW, proceedings. Vol. 7, pp. 438-453. #

Chandrasekhar, S. 1960. Radiative transfer. Dover Publications Inc. #

Smith, B. G. 1967. Geometrical shadowing of a random rough surface. IEEE Trans. on Antennas and Propagation, pp. 668–671. #

Trowbridge, T. S. and Reitz, K. P. 1975. Average irregularity representation of a rough surface for ray reflection. J. Opt. Soc. Am., 65(5):531–536. #

Landau, L. D. and Lifshitz, E. M. 1976. Mechanics. Vol. 1 (3rd ed.). Butterworth-Heinemann. #

砂川, 重信. 1977. 散乱の量子論. 岩波書店. #

Cook, R. L. and Torrance, K. E. 1982. A reflectance model for computer graphics. ACM Trans. Graph., 1(1):7–24. #

Kajiya, J. T. 1986. The rendering equation. SIGGRAPH Comput. Graph., 20(4):143–150. #

Cercignani, C. 1988. The Boltzmann Equation and Its Applications. Springer-Verlaq. #

藤田, 重次. 原, 啓明, et al. 共訳. 1989. 統計熱物理学. 裳華房. #

藤田, 重次. 原, 啓明, et al. 共訳. 1990. 量子統計物理学. 裳華房. #

Arvo, J. 1993. Transfer equations in global illumination. SIGGRAPH’93 Course Notes. #

Schlick, C. 1994. An inexpensive BRDF model for physically‐based rendering. In Computer graphics forum, vol. 13, no. 3, pp. 233-246. #

Born, M and Wolf, E. 1999. Principles of Optics (7th ed.). Cambridge University Press. #

Kardar, M. 2007. Statistical Physics of Particles. Cambridge University Press. #

Walter, B., et al. 2007. Microfacet Models for Refraction through Rough Surfaces. In Proceedings of the 18th Eurographics Conference on Rendering Techniques, EGSR’07, pp. 195–206. #

Heitz, E. 2014. Understanding the masking-shadowing function in microfacet-based BRDFs. Journal of Computer Graphics Techniques (JCGT), 3(2):48–107 #