Radiative Transfer方程式の拡散近似

これは レイトレ Advent Calendar 2019(https://qiita.com/advent-calendar/2019/raytracing)の記事です.

はじめに

BSSRDFモデル や subsurface scattering の理論的背景には, Radiative Transfer方程式の拡散近似があります[JMLH01]. ただし[JMLH01]ではこのことを軽く触れているのみで, 詳細は[Ish78]へ委ねられています. 一方で, Radiation Hydrodynamicsの本[Cas04]を見ると, この話題はモーメント定式化[Tho81]を用いてわかりやすく整理されています. そこで, この記事では [Tho81, Cas04]を参考にしてRadiative Transfer方程式のモーメント定式化を行い, [JMLH01]にある拡散方程式を求めます.

Radiative Transfer方程式

位置 {\mathbf x}, 時刻 t で方向 {\mathbf \omega} へ向かう光線の放射輝度を, L({\mathbf x},\,{\mathbf \omega},t) [\text{W}\cdot\text{m}^{-2}\cdot\text{sr}^{-1}] と書きます. このとき, 媒質内を進行する光線の放射輝度の変化は, Radiative Transfer方程式
{ \begin{align} \frac{1}{c}\frac{\partial L({\mathbf x},{\mathbf \omega},t)}{\partial t} + {\mathbf \omega}\cdot\nabla L({\mathbf x},{\mathbf \omega},t) = -\mu_e({\mathbf x}) L({\mathbf x},{\mathbf \omega},t) +\mu_s({\mathbf x}) L_s({\mathbf x},{\mathbf \omega},t) +S({\mathbf x},{\mathbf \omega},t), \label{eq:radiative_transfer}\tag{1} \end{align} }
{ \begin{align*} \qquad\qquad\qquad\qquad\qquad L_s({\mathbf x},{\mathbf \omega},t) := \int_{S^2} d{\mathbf \omega}'\,p({\mathbf \omega},\,{\mathbf \omega}')L({\mathbf x},{\mathbf \omega}',t) \end{align*} }

によって与えられます. ここで \int_{S^2}d{\mathbf \omega}=\int^{2\pi}_{\varphi=0}d\varphi\int^{\pi}_{\theta=0}d\theta\sin\theta は単位球面 S^2 上の積分を表します. そして c [\text{m}\cdot\text{s}^{-1}] は光速, \mu_s [\text{m}^{-1}], \mu_a [\text{m}^{-1}], \mu_e = \mu_s + \mu_a は媒質の散乱係数, 吸収係数, 消散係数です. p({\mathbf \omega},{\mathbf \omega}') は散乱の位相関数で,
{ \begin{align} \int_{S^2}d{\mathbf \omega}\,p({\mathbf \omega}\cdot{\mathbf \omega}')=1 \label{eq:phase_function_normalization_condition}\tag{2} \end{align} }
と規格化されているとします. また散乱は軸対称 p({\mathbf \omega},{\mathbf \omega}')=p({\mathbf \omega}\cdot{\mathbf \omega}')の場合に限るとして, その異方性因子を
{ \begin{align} g := \int_{S^2}d{\mathbf \omega}({\mathbf \omega}\cdot{\mathbf \omega}')p({\mathbf \omega}\cdot{\mathbf \omega}') \label{eq:asymmetry_factor}\tag{3} \end{align} }
と定義しておきます[JMLH01].
なお, [JMLH01]にある Radiative Transfer方程式は放射輝度の時間変化を無視しているのですが, この記事では残して導出を進めます. ほかの分野と比較したくなった時に便利だからです.
Radiative Transfer方程式(\ref{eq:radiative_transfer})は, 放射輝度 L に関する積分微分方程式となっていることに加えて, 位置 {\mathbf x}\in\mathbb{R}^3, 時間 t\in\mathbb{R}, 方向 {\mathbf \omega}\in S^2 の6つの独立変数を持っています. 従って簡単には解けません. そこで次節では, Radiative Transfer方程式を方向 {\mathbf \omega} に依存しない方程式群へ展開する, モーメント定式化を導入します.

モーメント定式化

モーメント展開

単位球面上に値をとる関数 f({\mathbf \omega}) はモーメント展開によって,

{ \begin{align} f({\mathbf \omega}) = \frac{1}{4\pi}\int d{\mathbf \omega}f({\mathbf \omega}) + \frac{3}{4\pi}{\mathbf \omega}\cdot\int d{\mathbf \omega}{\mathbf \omega}f({\mathbf \omega}) + \cdots \end{align} }
と展開できます[Tho81]. この式の第1項, 第2項, \cdots は球面調和関数展開の 0次, 1次, \cdots に対応しています. 展開の各項に現れる量 \int_{S^2} d{\mathbf \omega}\,f({\mathbf \omega}), \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}f({\mathbf \omega}) を一般的に表すために, 関数f({\mathbf \omega})の l次のモーメントを
{ \begin{align} \mu_l\left[f({\mathbf \omega})\right] := \int_{S^2}d{\mathbf \omega}\,\underbrace{{\mathbf \omega}\otimes{\mathbf \omega}\otimes\cdots\otimes{\mathbf \omega}}_{l\text{ factors}}f({\mathbf \omega}) \label{eq:define_moment}\tag{4} \end{align} }
と導入します[JAM+10]. \otimes はテンソル積です.
これを使うと, 放射輝度 L({\mathbf x},{\mathbf \omega},t) は
{ \begin{align} L({\mathbf x},{\mathbf \omega},t) = \frac{1}{4\pi}E({\mathbf x},t) + \frac{3}{4\pi}{\mathbf \omega}\cdot{\mathbf E}_1({\mathbf x},t) + \cdots \label{eq:moment_expansion}\tag{5} \end{align} }
と展開されます. ここで L({\mathbf x},{\mathbf \omega},t) の 0次のモーメントを scalar irradiance E({\mathbf x},t), 1次のモーメントを vector irradiance {\mathbf E}_1({\mathbf x},t) と書きました[JMLH01]:
{ \begin{align} E({\mathbf x},t) &:= \mu_0\left[L({\mathbf x},{\mathbf \omega},t)\right] = \int_{S^2}d{\mathbf \omega}\,L({\mathbf x},{\mathbf \omega},t)\qquad[\text{W}\cdot\text{m}^{-2}] \label{eq:scalar_irradiance}\tag{6}\\ {\mathbf E}_1({\mathbf x},t) &:= \mu_1\left[L({\mathbf x},{\mathbf \omega},t)\right] = \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\,L({\mathbf x},{\mathbf \omega},t)\qquad[\text{W}\cdot\text{m}^{-2}] \label{eq:vector_irradiance}\tag{7} \end{align} }
後のために, この記事では, 2次のモーメントも tensor irradiance {\mathbf E}_2({\mathbf x},t)
{ \begin{align} {\mathbf E}_2({\mathbf x},t) &:= \mu_2\left[L({\mathbf x},{\mathbf \omega},t)\right] = \int_{S^2}d\omega\,{\mathbf \omega}\otimes{\mathbf \omega}\,L({\mathbf x},{\mathbf \omega},t)\qquad[\text{W}\cdot\text{m}^{-2}] \label{eq:tensor_irradiance}\tag{8} \end{align} }
と書いておきます. また, 自己放射 S({\mathbf x},{\mathbf \omega},t) の0次と1次のモーメントも
{ \begin{align} Q_0({\mathbf x},t) &:= \mu_0\left[S({\mathbf x},{\mathbf \omega},t)\right] = \int_{S^2}d{\mathbf \omega}\,S({\mathbf x},{\mathbf \omega},t) \label{eq:scalar_emission}\tag{9}\\ {\mathbf Q}_1({\mathbf x},t) &:= \mu_1\left[S({\mathbf x},{\mathbf \omega},t)\right] = \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\,S({\mathbf x},{\mathbf \omega},t) \label{eq:vector_emission}\tag{10} \end{align} }
と書いておきます.

Radiative Transfer方程式のモーメント定式化

それでは, (\ref{eq:define_moment})を使って Radiative Transfer方程式の0次と1次のモーメント式を求めていきます.

0次のモーメント式

Radiative Transfer方程式(\ref{eq:radiative_transfer})の0次のモーメントをとると,

{ \begin{align*} \frac{1}{c}\frac{\partial}{\partial t}\underbrace{\int_{S^2}d{\mathbf \omega}\,L}_{ = E\text{ (\ref{eq:scalar_irradiance})} } + \underbrace{\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\cdot\nabla L}_{ = \nabla\cdot{\mathbf E}_1 \text{ (\ref{eq:nablaL_0th_moment})} } = -\mu_t({\mathbf x})\underbrace{\int_{S^2}d{\mathbf \omega}\,L}_{=E\text{ (\ref{eq:scalar_irradiance})}} + \mu_s({\mathbf x})\underbrace{\int_{S^2}d{\mathbf \omega}\,L_s}_{=E\text{ (\ref{eq:pL_0th_moment})}} + \underbrace{\int_{S^2}d{\mathbf \omega}\,S}_{=Q_0\text{ (\ref{eq:scalar_emission})}}. \end{align*} }
そして下かっこで示した式の整理を行うと,
{ \begin{align} \frac{1}{c}\frac{\partial E({\mathbf x},t)}{\partial t} + \nabla\cdot{\mathbf E}_1({\mathbf x},t) = -\mu_a({\mathbf x})E({\mathbf x},t) + Q_0({\mathbf x},t) \label{eq:rte_0th_moment}\tag{11} \end{align} }
となります.

1次のモーメント式

同様に, Radiative Transfer方程式(\ref{eq:radiative_transfer})の1次のモーメントをとると,

{ \begin{align*} \frac{1}{c}\frac{\partial}{\partial t} \underbrace{\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}L}_{ = {\mathbf E}_1\text{ (\ref{eq:vector_irradiance})} } + \underbrace{\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\left({\mathbf \omega}\cdot\nabla L\right)}_{ = \nabla\cdot{\mathbf E}_2 \text{ (\ref{eq:nablaL_1st_moment})} } = -\mu_t({\mathbf x}) \underbrace{\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}L}_{ = {\mathbf E}_1\text{ (\ref{eq:vector_irradiance})} } + \mu_s({\mathbf x})\underbrace{\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}L_s}_{=g{\mathbf E}_1\text{ (\ref{eq:pL_1st_moment})}} + \underbrace{\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}S}_{=Q_1\text{ (\ref{eq:vector_emission})}} \end{align*} }
そして下かっこで示した式の整理を行うと,
{ \begin{align} \frac{1}{c}\frac{\partial{\mathbf E}_1({\mathbf x},t)}{\partial t} + \nabla\cdot{\mathbf E}_2({\mathbf x},t) = -(\mu_a({\mathbf x}) + \mu'_s({\mathbf x})){\mathbf E}_1({\mathbf x},t) + {\mathbf Q}_1({\mathbf x},t) \label{eq:rte_1st_moment}\tag{12} \end{align} }
となります.ここで \mu'_s:=\mu_s(1-g) としました.

Radiative Transfer方程式の拡散近似

前節で得た 0次のモーメント式には1次の量 {\mathbf E}_1 が, 1次のモーメント式には2次の量 {\mathbf E}_2 が表れます. このままではモーメント式の階層が無限に続いてしまうので, 放射輝度のモーメント展開(\ref{eq:moment_expansion})を1次で打ち切る近似[JMLH01]を行います
{ \begin{align} L({\mathbf x},{\mathbf \omega},t) = \frac{1}{4\pi}E({\mathbf x},t) + \frac{3}{4\pi}{\mathbf \omega}\cdot{\mathbf E}_1({\mathbf x},t) \label{eq:moment_expansion_1st_approximation}\tag{13} \end{align} }
この近似は, 拡散近似と呼ばれています[Cas04]. これを, 1次のモーメント式(\ref{eq:rte_1st_moment})の第2項へ適用すると, 下記が得られます.

{ \begin{align*} \nabla\cdot{\mathbf E}_2({\mathbf x},t) = \frac{1}{3}\nabla E({\mathbf x},t) \end{align*} }
詳細
{ \begin{align*} \nabla\cdot{\mathbf E}_2({\mathbf x},t) &= \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\left[{\mathbf \omega}\cdot\nabla L({\mathbf x},{\mathbf \omega},t)\right]\\ &= \frac{1}{4\pi}\underbrace{\int_{S^2}d{\mathbf \omega}{\mathbf \omega}\left({\mathbf \omega}\cdot\nabla E({\mathbf x},t)\right)}_{ = \frac{4\pi}{3}\nabla E({\mathbf x},t)\text{ (\ref{eq:linear_function_1st_moment})} } +\frac{3}{4\pi}\int_{S^2}d{\mathbf \omega}{\mathbf \omega}\left[ {\mathbf \omega}\cdot\nabla\left({\mathbf \omega}\cdot{\mathbf E}_1({\mathbf x},t)\right) \right]\\ &= \frac{1}{3}\nabla E({\mathbf x},t) +\frac{3}{4\pi}\int_{S^2}d{\mathbf \omega}{\mathbf \omega} [{\mathbf \omega}^T\underbrace{\left(\nabla{\mathbf \omega}\right)}_{=0}{\mathbf E}_1({\mathbf x},t)] +\frac{3}{4\pi} \underbrace{ \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\left[{\mathbf \omega}^T\left(\nabla{\mathbf E}_1({\mathbf x},t)\right){\mathbf \omega}\right] }_{=0\text{ (\ref{eq:quadratic_form_1st_moment})}}\\ &= \frac{1}{3}\nabla E({\mathbf x},t) \end{align*} }
\Box

さらに, 平均自由行程 1/(\mu_a + \mu'_s) [\text{m}] のあいだで {\mathbf E}_1 の時間変化が無視できるくらい小さい \left|\frac{1}{c}\frac{\partial{\mathbf E}_1}{\partial t}\right| \ll \left|(\mu_a + \mu'_s){\mathbf E}_1\right| とすれば, 1次のモーメント式 (\ref{eq:rte_1st_moment}) は
{ \begin{align*} \frac{1}{3}\nabla E({\mathbf x},t) = -(\mu_a({\mathbf x}) + \mu'_s({\mathbf x})){\mathbf E}_1({\mathbf x},t) + {\mathbf Q}_1({\mathbf x},t) \end{align*} }
となります. これを変形して
{ \begin{align*} {\mathbf E}_1({\mathbf x},t) = 3D({\mathbf x}){\mathbf Q}_1({\mathbf x},t) - D({\mathbf x})\nabla E({\mathbf x},t) ,\quad D({\mathbf x}):=\frac{1}{3(\mu_a({\mathbf x})+\mu'_s({\mathbf x}))} \end{align*} }
と書き, 0次のモーメント式 (\ref{eq:rte_0th_moment}) へ代入すれば
{ \begin{align*} \frac{1}{c}\frac{\partial E({\mathbf x},t)}{\partial t} - \nabla\cdot\left[D({\mathbf x})\nabla E({\mathbf x},t)\right] = -\mu_a({\mathbf x})E({\mathbf x},t) + Q_0({\mathbf x},t) - 3\nabla\cdot\left[D({\mathbf x}){\mathbf Q}_1({\mathbf x},t)\right] \end{align*} }
が得られます. 同様にして E の時間変化も無視でき, 媒質が一様である \mu_a({\mathbf x})\rightarrow\mu_a, \mu_s({\mathbf x})\rightarrow\mu_s と仮定すれば,
{ \begin{align} D\nabla^2 E({\mathbf x},t) = \mu_aE({\mathbf x},t) - Q_0 + 3D\nabla\cdot{\mathbf Q}_1({\mathbf x},t) \label{eq:isotropic_diffusion_equation}\tag{14} \end{align} }
となります. これで, 文献[Ish78, JMLH01, JAM+10]の拡散方程式が求まりました.

まとめ

この記事では, 放射輝度の低次近似(\ref{eq:moment_expansion_1st_approximation})を取ることで, Radiative Transfer方程式(\ref{eq:radiative_transfer})から拡散方程式(\ref{eq:isotropic_diffusion_equation})が得られることを確認しました. ここで書いた内容が, Radiative Transfer方程式の拡散近似を用いた文献[JMLH01, JAM+10]などを読む際の助けになれば幸いです.

付録

この記事内で使った式を挙げます.

内積 {\mathbf a}\cdot{\mathbf \omega} と 2次形式 {\mathbf \omega}^T{\mathbf A}{\mathbf \omega} の 0次と1次のモーメント

[JAM+10] から, これらの関係を引用します.
{ \begin{align} \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\cdot{\mathbf a}=0 \label{eq:linear_function_0th_moment}\tag{A.1} \end{align} }
{ \begin{align} \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}({\mathbf \omega}\cdot{\mathbf a})=\frac{4\pi}{3}{\mathbf a} \label{eq:linear_function_1st_moment}\tag{A.2} \end{align} }
{ \begin{align} \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}^T{\mathbf A}{\mathbf \omega} = \frac{4\pi}{3}\text{Tr}({\mathbf A}) \label{eq:quadratic_form_0th_moment}\tag{A.3} \end{align} }
{ \begin{align} \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\left({\mathbf \omega}^T{\mathbf A}{\mathbf \omega}\right) = 0 \label{eq:quadratic_form_1st_moment}\tag{A.4} \end{align} }

位相関数 p({\mathbf \omega}\cdot{\mathbf \omega}') の1次のモーメント

{ \begin{align} \int_{S^2}d{\mathbf \omega}{\mathbf \omega}p({\mathbf \omega}\cdot{\mathbf \omega}') = {\mathbf \omega}'g \label{eq:phase_function_1st_moment}\tag{A.5} \end{align} }
p({\mathbf \omega}\cdot{\mathbf \omega}') の1次のモーメントをとると,

{ \begin{align*} \mu_1\left[p({\mathbf \omega},{\mathbf \omega}')\right] &= \int_{S^2}d{\mathbf \omega}{\mathbf \omega}p({\mathbf \omega}\cdot{\mathbf \omega}')\\ &= \int_{S^2}d{\mathbf \omega}\left[ {\mathbf \omega}'({\mathbf \omega}\cdot{\mathbf \omega}') + {\mathbf \omega}'\times({\mathbf \omega}\times{\mathbf \omega}') \right]p({\mathbf \omega}\cdot{\mathbf \omega}')\\ &= {\mathbf \omega}'\underbrace{ \int_{S^2}d{\mathbf \omega}({\mathbf \omega}\cdot{\mathbf \omega}')p({\mathbf \omega}\cdot{\mathbf \omega}') }_{=g \text{ (\ref{eq:asymmetry_factor})}} + {\mathbf \omega}'\times\left[ \int_{S^2}d{\mathbf \omega}{\mathbf \omega}p({\mathbf \omega}\cdot{\mathbf \omega}') \right]\times{\mathbf \omega}' \end{align*} }
初めの等号でベクトル3重積 {\mathbf \omega}={\mathbf \omega}'({\mathbf \omega}\cdot{\mathbf \omega}') + {\mathbf \omega}'\times({\mathbf \omega}\times{\mathbf \omega}') を使いました. また,
{ \begin{align*} {\mathbf \omega}'\times\int_{S^2}d{\mathbf \omega}{\mathbf \omega}p({\mathbf \omega}\cdot{\mathbf \omega}') \underbrace{\sim}_{\text{(\ref{eq:linear_function_1st_moment})}} {\mathbf \omega}'\times{\mathbf \omega}' = 0 \end{align*} }
を使うと,
{ \begin{align*} \int_{S^2}d{\mathbf \omega}{\mathbf \omega}p({\mathbf \omega}\cdot{\mathbf \omega}') = {\mathbf \omega}'g \end{align*} }
となります. \Box

放射輝度 \nabla L({\mathbf x},{\mathbf \omega}) の0次と1次のモーメント

{ \begin{align} \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\cdot\nabla L({\mathbf x},{\mathbf \omega}) = \nabla\cdot\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}L \label{eq:nablaL_0th_moment}\tag{A.6} \end{align} }
{ \begin{align} \int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\left[{\mathbf \omega}\cdot\nabla L({\mathbf x},{\mathbf \omega})\right] = \nabla\cdot\int_{S^2}d{\mathbf \omega}\,{\mathbf \omega}\otimes{\mathbf \omega}L({\mathbf x},{\mathbf \omega}) \label{eq:nablaL_1st_moment}\tag{A.7} \end{align} }
これらは, \nabla が {\mathbf \omega} に作用しないことから得られます.

L_s({\mathbf x},{\mathbf \omega}) = \int_{S^2}d{\mathbf \omega}'p({\mathbf \omega}\cdot{\mathbf \omega}')L({\mathbf x},{\mathbf \omega}') の0次と1次のモーメント

{ \begin{align} \int_{S^2}d{\mathbf \omega}L_s({\mathbf x},{\mathbf \omega}) = E({\mathbf x}) \label{eq:pL_0th_moment}\tag{A.8} \end{align} }
{ \begin{align} \int_{S^2}d{\mathbf \omega}{\mathbf \omega}L_s({\mathbf x},{\mathbf \omega}) = g{\mathbf E}_1({\mathbf x}) \label{eq:pL_1st_moment}\tag{A.9} \end{align} }
(\ref{eq:pL_0th_moment}): L_s の0次のモーメントを取ると

{ \begin{align*} \mu_1\left[L_s({\mathbf x},{\mathbf \omega})\right] &= \int_{S^2}d{\mathbf \omega}'L({\mathbf x},{\mathbf \omega}') \underbrace{ \int_{S^2}d{\mathbf \omega}p({\mathbf \omega}\cdot{\mathbf \omega}') }_{=1 \text (\ref{eq:phase_function_normalization_condition})}\\ &= E({\mathbf x}) \end{align*} }
となります. \Box
(\ref{eq:pL_1st_moment}): L_s の1次のモーメントを取ると
{ \begin{align*} \mu_1\left[L_s({\mathbf x},{\mathbf \omega})\right] &= \int_{S^2}d{\mathbf \omega}'L({\mathbf x},{\mathbf \omega}') \underbrace{\int_{S^2}d{\mathbf \omega}{\mathbf \omega}p({\mathbf \omega}\cdot{\mathbf \omega}')}_{ =g{\mathbf \omega}'\text{ (\ref{eq:phase_function_1st_moment})} }\\ &= g\int_{S^2}d{\mathbf \omega}'{\mathbf \omega}'L({\mathbf x},{\mathbf \omega}')\\ &= g{\mathbf E}_1({\mathbf x}) \end{align*} }
となります. \Box

参考文献

[Cas04] John I. Castor. Radiation Hydrodynamics. Cambridge University Press (2004).
[Ish78] Akira Ishimaru. Wave Propagation and Scattering in Random Media, Chapter 9 - diffusion approximation. Academic Press (1978).
[JAM+10] Wenzel Jakob, Adam Arbree, Jonathan T. Moon, Kavita Bala, and Steve Marschner. A radiative transfer framework for rendering materials with anisotropic structure. In SIGGRAPH 2010, ACM, New York, NY, USA, 53:1–53:13 (2010).
[JMLH01] Henrik Wann Jensen, Stephen R. Marschner, Marc Levoy, and Pat Han-rahan. A practical model for subsurface light transport. SIGGRAPH ’01, pages 511–518 (2001).
[Tho81]Kip S. Thorne. Relativistic radiative transfer: moment formalisms. Monthly Notices of the Royal Astronomical Society, 194(2):439–473, 02 (1981).

光子の集団としての光線について

これは レイトレ合宿7(https://sites.google.com/site/raytracingcamp7/ )のアドベントカレンダーの記事として作成しました .

はじめに

パストレーシングにおいて, 光は輝度 L を運ぶ光線として表され, その光散乱をRendering方程式
{ \begin{align} L({\mathbf x},{\mathbf\omega}_o) = L_e({\mathbf x},{\mathbf \omega}_o) +  \int f_r({\mathbf x},{\mathbf \omega}_i,{\mathbf \omega}_o) L_i({\mathbf x},{\mathbf \omega}_i)({\mathbf n}\cdot{\mathbf \omega}_i)d{\mathbf \omega}_i, \tag{1} \end{align} }
ならびに Volume Rendering 方程式に基づいて計算することで, 写実的な映像を作り出すことができます. これらの基礎となるものは Radiative Transfer方程式であり[Arvo1996], 更にそれはBoltzmann方程式の一種であることが知られています[Cercignani1988]. Boltzmann方程式とは希薄気体の運動方程式です. その運動は, 相空間分布関数\rho({\mathbf x},{\mathbf p},t)として表現した粒子の集団の時間変化として表現されるのですが, それならその分布関数と光線の概念は何らかの関係で結びついているはずです.
この記事では, Boltzmann方程式とRendering方程式の関係を考えるための準備として, 光子の集団を表す分布関数 \rho_\text{photon}({\mathbf x},{\mathbf p},t) と, 輝度 L を運ぶ光線との間の対応関係について考えます.

相空間分布関数

いろいろな運動量を持って動き回る点粒子の集団を考えます. このとき, 位置 {\mathbf x}, 時間 tにおいて, 運動量が {\mathbf p} である粒子数の密度を, 相空間分布関数 \rho({\mathbf x},{\mathbf p},t) として表現します. この関数は, 運動量で積分すれば粒子数密度 n({\mathbf x})
{ \begin{align} \frac{1}{h^3}\int \rho({\mathbf x},{\mathbf p},t)d{\mathbf p} = n({\mathbf x}) \quad[\text{m}^{-3}] \label{eq:number_density}\tag{2} \end{align} }
を, 相空間体積要素で積分すれば, 粒子数 N
{ \begin{align} \frac{1}{h^3}\iint \rho({\mathbf x},{\mathbf p},t)d{\mathbf x}d{\mathbf p} = N \quad [\text{-}]\tag{3} \end{align} }
を与えるように規格化されているとします.
ここで, d{\mathbf x}d{\mathbf p}=dxdydzdp_xdp_ydp_z の単位は [\text{J}^3\cdot\text{s}^3] なので, 分布関数を無次元にするために, プランク定数 h [\text{J}\cdot\text{s}] を使って \frac{1}{h^3} を掛けました.

光線中の光子の数

光線を一本の光のビームと考えて, それを構成する光子の分布関数が \rho_\text{photon} で与えられたとします. 光は振動の方向に応じて2つの偏光状態を取ることができるので, これを2つの異なる状態と数えて, 分布関数は
{ \begin{align} \frac{2}{h^3}\iint \rho_\text{photon}({\mathbf x},{\mathbf p},t)d{\mathbf x}d{\mathbf p} = N \quad[\text{-}]\tag{4} \end{align} }
と規格化されているとします.
(\ref{eq:number_density})より, この式は(光子の個数密度)\times(空間体積要素)という意味を持っているので, 光子のエネルギー \epsilon({\mathbf p}) から, 領域 d{\mathbf x}d{\mathbf p} 内の放射エネルギーdQ[\text{J}]を,
{ \begin{align} dQ=\frac{2}{h^3}\epsilon({\mathbf p}) \rho_\text{photon}({\mathbf x},{\mathbf p},t)d{\mathbf x}d{\mathbf p} \quad[\text{J}]\tag{5} \end{align} }
と計算できるはずです.
そして, この量を次のように書き換えます.

f:id:kenha:20190707190236p:plain
図: 立体角 d\Omega 内を進み, 面積dA, 法線 {\mathbf n} の微小面を通過する光子
光のビームが, 微小立体角 d\Omega 内の方向 {\mathbf \omega} へ光速 cで進み, 面積 dA, 法線 {\mathbf n} の微小面を通過するとします(図). このとき, 相空間体積素片 d{\mathbf x} は, dt の間にビームが掃き進んだ体積として, d{\mathbf x}=cdt({\mathbf n}\cdot{\mathbf\omega})dA [\text{m}^{3}]と書けます. さらに運動量体積要素 d{\mathbf p} も, 極座標表示を使って |{\mathbf p}|^2d|{\mathbf p}|d\Omega [\text{J}^{3}\cdot\text{m}^{-3}\cdot\text{s}^3\cdot\text{sr}]とします. そうすると, 放射エネルギーは
{ \begin{align} dQ=\frac{2c}{h^3}\epsilon({\mathbf p}) \rho_\text{photon}({\mathbf x},{\mathbf p},t) ({\mathbf n}\cdot{\mathbf \omega})|{\mathbf p}|^2d|{\mathbf p}|d\Omega dAdt \quad[\text{J}]\tag{6} \end{align} }
となります.
de Broglieの関係, および特殊相対性理論におけるエネルギーと運動量の関係によれば, 振動数が \nu である光子のエネルギー \epsilon({\mathbf p}) には,
{ \begin{align} \epsilon({\mathbf p}) = c|{\mathbf p}| = h\nu \quad[\text{J}]\tag{7} \end{align} }
という関係があります. 従って \epsilon({\mathbf p})=h\nu, そして |{\mathbf p}|^2d|{\mathbf p}|d\Omega = \frac{h^3\nu^2}{c^3}d\nu d\Omega
と書けば, 放射エネルギーは
{ \begin{align} dQ = L_\nu({\mathbf x},{\mathbf\omega},t)({\mathbf n}\cdot{\mathbf \omega}) d\nu d\Omega dAdt \quad[\text{J}] \label{eq:radiant_energy}\tag{8} \end{align} }
と得られます. ここで導入した分光放射輝度 L_\nu
{ \begin{align} L_\nu({\mathbf x},{\mathbf \omega},t) := \frac{2h\nu^3}{c^2}\rho_\text{photon}({\mathbf x},{\mathbf p},t)\quad [\text{W}\cdot\text{m}^{-2}\cdot\text{sr}^{-1}\cdot\text{Hz}^{-1}] \label{eq:spectral_radiance}\tag{9} \end{align} }
によって, 光のビームが運ぶ輝度と, 光子の分布関数の関係がわかりました.

黒体輻射

(\ref{eq:spectral_radiance}) により, 光子の分布関数を与えれば, その分光放射輝度を求めることができます. 温度 T [\text{K}]で熱平衡状態にある光子の分布関数, Planck分布
{ \begin{align} \rho_\text{photon}({\mathbf x},{\mathbf p},t) = \frac{1}{e^{\frac{h\nu}{kT}}-1} \quad[\text{-}]\tag{10} \end{align} }
を(\ref{eq:spectral_radiance})へ代入すれば, 分光放射輝度はただちに
{ \begin{align} L_\nu({\mathbf x},{\mathbf \omega},t) = \frac{2h\nu^3}{c^2}\frac{1}{e^{\frac{h\nu}{kT}}-1} \quad [\text{W}\cdot\text{m}^{-2}\cdot\text{sr}^{-1}\cdot\text{Hz}^{-1}] \label{eq:blackbody_radiation}\tag{11} \end{align} }
と黒体輻射の式として求まります. ここで, k はボルツマン定数です.

さいごに

この記事では, パストレーシングにおける光線を光子の集団によるビームと考えて, 放射エネルギーと分光放射輝度を導入(\ref{eq:radiant_energy})し, 光子の分布関数との対応関係(\ref{eq:spectral_radiance})を求めました. この内容は, 主に[藤田1990]をもとにしています. 興味のある方は, そちらも併せて参考にしてください.

参考文献

[Arvo1996] James Arvo. Transfer equations in global illumination. Global Illumination, SIGGRAPH’93 Course Notes, 01 1996.
[Cercignani1988] Carlo Cercignani. The Boltzmann Equation and Its Applications. Springer-Verlaq, 1988.
[藤田1990] 藤田重次. 量子統計物理学 第4版. 裳華房, 1990. 原啓明, 小幡常啓, 岡村好庸, 鈴木彰共訳.

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

これは レイトレ 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:いま, 面積を求める計算を行っているからです.