このページはレイトレ Advent Calendar 2017の記事として制作しました。
今回はひたすら数式です。現時点では数式の導出過程の式変形などの検証が主で、「なぜそうしたか」といった論理の展開については深く追っていません。
BSSRDFのダイポール(dipole)モデルなどはその理論の根底に、拡散方程式による表面下の光の分布の表現・近似を使っています。ダイポールモデル等を単に実装するだけであれば示された理論をブラックボックス的に扱えば良く、あまり深く追う必要は無いかもしれませんが、発展的な論文などを読む上では自分の中に確たる理解を持っていれば安心して読み進められます。また個人的な経験として、BSSRDF関連の資料では結構数式が誤まって記述されていたり、表記のルール・前提が資料ごとに異なっていたりすることが散見されたため、正しい実装をする上では自分で理論をしっかり追ってみるのが重要と感じています。このページでは、表面下の光の分布を記述する拡散方程式の導出を行います。また、BSSRDFのコアとなるdiffuse reflectanceの式の導出も行います。
「ボリュームレンダリング方程式」でも解説したように、関与媒質中で光が微小距離 $ ds $ 進んだ際の放射輝度 $ L $ の変化は次の放射伝達方程式(RTE)で表されます。 \begin{equation*} \frac{\d L(\vx, \vomega)}{\d s} = -\sigma_e(\vx) L(\vx, \vomega) + \sigma_s(\vx)\int_{\mathcal{S}^2} f_p(\vx, \vomega', \vomega) L(\vx, \vomega') \d\vomega' + L_e^V(\vx, \vomega) \end{equation*} これを勾配 $ \nabla $ を使うことで次のようにも表すことができます。 \begin{equation} (\vomega \cdot \nabla) L(\vx, \vomega) = -\sigma_e(\vx) L(\vx, \vomega) + \sigma_s(\vx)\int_{\mathcal{S}^2} f_p(\vx, \vomega', \vomega) L(\vx, \vomega') \d\vomega' + L_e^V(\vx, \vomega) \label{eq:RTE} \end{equation} 若干難しそうな記述になりましたが、勾配をとることで3次元の各方向ごとの変化率を表すベクトルが得られるので、それと方向ベクトル $ \vomega $ の内積をとればある方向 $ \vomega $ への変化率がわかる、という話です。
各パラメターに関しておさらいしておくと、$ \sigma_e(\vx) = \sigma_a(\vx) + \sigma_s(\vx) $ は消散係数で媒質の密度を表します。$ \sigma_a(\vx) $ は吸収係数で吸収する媒質の密度、$ \sigma_s(\vx) $ は散乱係数で散乱する媒質の密度を表します。$ f_p(\vx, \vomega', \vomega) $ は位相関数で散乱時の方向分布を表します。位相関数の平均コサインが次の形で定義されます。 \begin{equation*} g = \int_{\mathcal{S}^2} (\vomega \cdot \vomega') \, f_p(\vomega', \vomega) \d\vomega' \end{equation*} $ g $ は $ [-1, 1] $ の範囲の値をとり、後方散乱の場合 $ -1 $ に、前方散乱の場合 $ 1 $ に近づきます。$ g = 0 $ の場合は完全等方散乱です。
「双方向散乱面反射率分布関数 (BSSRDF)」でも述べたように、表面下散乱物質中では光の散乱回数に応じて光を分類することが意味を持ちます。媒質中の光 $ L_(\vx, \vomega) $ を次のように2つの光に分解します。 \begin{equation*} L(\vx, \vomega) = L_{ri}(\vx, \vomega) + L_d(\vx, \vomega) \end{equation*} $ L_{ri}(\vx) $ は媒質に入射してから一回も散乱を起こしていない光を表しており、reduced incident radiance (reduced intensity, 減衰入射放射輝度?)と呼ばれます。一方 $ L_d(\vx, \vomega) $ は一回以上散乱した光でdiffuse radiance (拡散放射輝度?)と呼びます。一様な媒質の場合の場合、reduced incident radianceは次の式で表されます。 \begin{equation*} L_{ri}(\vx, \vomega) = \exp(-\sigma_e s) L(\raycast(\vx, -\vomega), \vomega) \end{equation*} ここで、$ \raycast(\vx, -\vomega) $ はレイキャスト関数で、$ \vx $ から $ -\vomega $ 方向にレイを飛ばして最初に見つかる物体表面の位置を得る関数です。つまり $ L(\raycast(\vx, -\vomega), \vomega) $ は物体に入射した直後の $ \vomega $ 方向への放射輝度を表しています。$ s $ は $ \vx $ と $ \raycast(\vx, -\vomega) $ の間の距離を表しています。つまり、媒質に入射した光は $ \vx $ に至るまで消散係数 $ \sigma_e $ に応じて指数的に減衰します。
diffuse radianceに関して 式 $ \eqref{eq:RTE} $ と同様の式を立てることができます。 \begin{equation} (\vomega \cdot \nabla) L_d(\vx, \vomega) = -\sigma_e(\vx) L_d(\vx, \vomega) + \sigma_s(\vx)\int_{\mathcal{S}^2} f_p(\vx, \vomega', \vomega) L_d(\vx, \vomega') \d\vomega' + Q(\vx, \vomega) \label{eq:diffuse_RTE} \end{equation} ここで $ Q(\vx, \vomega) $ はreduced incident radianceが初めて散乱してdiffuse radianceの一部となることを表す項で、次のように表されます。 \begin{equation*} Q(\vx, \vomega) = \sigma_s(\vx)\int_{\mathcal{S}^2} f_p(\vx, \vomega', \vomega) L_{ri}(\vx, \vomega') \d\vomega' \end{equation*}
式 $ \eqref{eq:diffuse_RTE} $ の両辺を全方位に関して積分します。まず左辺について考えます。 \begin{eqnarray*} \int_{\mathcal{S}^2} (\vomega \cdot \nabla) L_d(\vx, \vomega) \d\vomega &=& \int_{\mathcal{S}^2} \vomega \cdot \begin{pmatrix} \frac{\partial}{\partial x} L_d(\vx, \vomega) \\ \frac{\partial}{\partial y} L_d(\vx, \vomega) \\ \frac{\partial}{\partial z} L_d(\vx, \vomega) \\ \end{pmatrix} \d\vomega \\ &=& \int_{\mathcal{S}^2} \Prt{ \omega_x \frac{\partial}{\partial x} L_d(\vx, \vomega) + \omega_y \frac{\partial}{\partial y} L_d(\vx, \vomega) + \omega_z \frac{\partial}{\partial z} L_d(\vx, \vomega) } \d\vomega \\ &=& \int_{\mathcal{S}^2} \nabla \cdot \begin{pmatrix} \omega_x L_d(\vx, \vomega) \\ \omega_y L_d(\vx, \vomega) \\ \omega_z L_d(\vx, \vomega) \\ \end{pmatrix} \d\vomega \\ &=& \nabla \cdot \int_{\mathcal{S}^2} L_d(\vx, \vomega) \vomega \d\vomega \\ \end{eqnarray*} 次に右辺について考えます。 \begin{eqnarray*} && -\sigma_e(\vx) \int_{\mathcal{S}^2} L_d(\vx, \vomega) \d\vomega + \sigma_s(\vx) \int_{\mathcal{S}^2} \int_{\mathcal{S}^2} f_p(\vx, \vomega', \vomega) L_d(\vx, \vomega') \d\vomega' \d\vomega + \int_{\mathcal{S}^2} Q(\vx, \vomega) \d\vomega \\ &=& -\sigma_e(\vx) \int_{\mathcal{S}^2} L_d(\vx, \vomega) \d\vomega + \sigma_s(\vx) \int_{\mathcal{S}^2} L_d(\vx, \vomega') \d\vomega' + \int_{\mathcal{S}^2} Q(\vx, \vomega) \d\vomega \\ &=& -\sigma_a(\vx) \int_{\mathcal{S}^2} L_d(\vx, \vomega) \d\vomega + \int_{\mathcal{S}^2} Q(\vx, \vomega) \d\vomega \\ \end{eqnarray*} 1行目から2行目への2項目の変形では重積分の積分範囲が変数に依存していないことと、位相関数が正規化されていることを利用しています。以上により次の関係を得ます。 \begin{eqnarray*} \nabla \cdot \int_{\mathcal{S}^2} L_d(\vx, \vomega) \vomega \d\vomega &=& -\sigma_a(\vx) \int_{\mathcal{S}^2} L_d(\vx, \vomega) \d\vomega + \int_{\mathcal{S}^2} Q(\vx, \vomega) \d\vomega \end{eqnarray*} この式を次のように表します。 \begin{eqnarray} \nabla \cdot \vec{E}(\vx) &=& -\sigma_a(\vx) \phi(\vx) + Q_0(\vx) \label{eq:divergence_of_vector_irradiance} \end{eqnarray} ここで $ \phi(\vx) $ と $ \vec{E}(\vx) $ はそれぞれradiant fluence (scalar irradiance)とvector irradiance (vector radiant flux)と呼ばれ、次の式のように定義されます。 \begin{eqnarray*} \phi(\vx) = \int_{\mathcal{S}^2} L_d(\vx, \vomega) \d\vomega, \hspace{5mm} \vec{E}(\vx) = \int_{\mathcal{S}^2} L_d(\vx, \vomega) \vomega \d\vomega \end{eqnarray*} また、$ Q_0(\vx) $ は光源項 $ Q(vx, \vomega) $ の0次の項で次のように表されます。また、後で使うことになる光源項の1次の項 $ \vec{Q}_1(\vx) $ も併記します。 \begin{eqnarray*} Q_0(\vx) = \int_{\mathcal{S}^2} Q(\vx, \vomega) \d\vomega, \hspace{5mm} \vec{Q}_1(\vx) = \int_{\mathcal{S}^2} Q(\vx, \vomega) \vomega \d\vomega \end{eqnarray*} $ Q_0(\vx) $ は $ Q(\vx, \vomega) $ のうち方向によらず常に存在する成分を表しています。
散乱性の非常に高い物質内では、例え入射光の分布や位相関数形状が指向性を持っていた(anisotropic)としても、何度も散乱しているうちに、光の方向に関する分布は一様になっていきます。diffuse radianceは一様な分布に近づいていると考えられるため、実球面調和関数展開し、1次の球面調和関数までで近似します。球面調和関数展開についてはここでは省略しますが(そのうち書くかもしれない...)、フーリエ級数展開と似たような概念です。フーリエ級数展開の球面版のようなものです。 \begin{eqnarray*} L_d(\vx, \vomega) \approx \frac{1}{4\pi} \phi(\vx) + \frac{3}{4\pi} \vec{E}(\vx) \cdot \vomega \end{eqnarray*} これをRTE 式 $ \eqref{eq:diffuse_RTE} $ に代入し、両辺に $ \vomega $ をかけて全方位で積分します。 \begin{eqnarray*} \int_{\mathcal{S}^2} \vomega (\vomega \cdot \nabla) \Prtc{\frac{1}{4\pi} \phi(\vx) + \frac{3}{4\pi} \vec{E}(\vx) \cdot \vomega} \d\vomega &=& -\int_{\mathcal{S}^2} \sigma_e(\vx) \vomega \Prtc{\frac{1}{4\pi} \phi(\vx) + \frac{3}{4\pi} \vec{E}(\vx) \cdot \vomega} \d\vomega \\ && + \int_{\mathcal{S}^2} \sigma_s(\vx) \vomega \int_{\mathcal{S}^2} f_p(\vx, \vomega', \vomega) \Prtc{\frac{1}{4\pi} \phi(\vx) + \frac{3}{4\pi} \vec{E}(\vx) \cdot \vomega} \d\vomega' \d\vomega \\ && + \int_{\mathcal{S}^2} \vomega Q(\vx, \vomega) \d\vomega \end{eqnarray*} 両辺を整理すると次の式を得ます。 \begin{eqnarray*} \underbrace{ \frac{1}{4\pi} \int_{\mathcal{S}^2} \vomega (\vomega \cdot \nabla) \phi(\vx) \d\vomega }_{(a)} + \underbrace{ \frac{3}{4\pi} \int_{\mathcal{S}^2} \vomega (\vomega \cdot \nabla) \Prtb{\vec{E}(\vx) \cdot \vomega} \d\vomega }_{(b)} &=& -\underbrace{ \frac{\sigma_a(\vx)}{4\pi} \int_{\mathcal{S}^2} \vomega \phi(\vx) \d\vomega }_{(c)} - \underbrace{ \frac{3 \sigma_e(\vx)}{4\pi} \int_{\mathcal{S}^2} \vomega \Prtb{\vec{E}(\vx) \cdot \vomega} \d\vomega }_{(d)} \\ && + \underbrace{ \frac{3 \sigma_s(\vx)}{4\pi} \int_{\mathcal{S}^2} \vomega \int_{\mathcal{S}^2} f_p(\vx, \vomega', \vomega) \Prtb{\vec{E}(\vx) \cdot \vomega'} \d\vomega' \d\vomega }_{(e)} \\ && + \underbrace{ \int_{\mathcal{S}^2} \vomega Q(\vx, \vomega) \d\vomega }_{(f)} \end{eqnarray*} この式の各項 $ (a) \sim (f) $ をさらに整理していきます。それぞれの項を解く際には方向成分 $ \vomega $ に関連した積分を何度も解くことになりますが、それにあたっては次に示す性質が使えます。
方向 $ \vomega $ は球面座標を用いて次のように表される。 \begin{eqnarray*} \vomega = (\omega_x, \omega_y, \omega_z) = (\cos\phi \sin\theta, \sin\phi \sin\theta, \cos\theta) \end{eqnarray*} 方向の各成分の積の組み合わせに関する積分は単純に次のように書くことができる。 \begin{eqnarray*} \int_{\mathcal{S}^2} \omega_x^a \omega_y^b \omega_z^c \d\vomega = \int_{0}^{2\pi} \int_{0}^{\pi} \omega_x^a \omega_y^b \omega_z^c \sin\theta \d\theta \d\phi \end{eqnarray*} ある $ a, b, c $ の組み合わせの積分の解き方は、至って普通の置換積分や部分積分になるので省略するが、0次から3次までの答えは次のようになる。
上記の性質を用いて $ (a) \sim (f) $ を整理すると $ (a) \sim (d) $ は特筆すべきことも無く次の形となります。 \begin{eqnarray*} (a) &\cdots& \frac{1}{3} \nabla \phi(\vx) \\ (b) &\cdots& 0 \\ (c) &\cdots& 0 \\ (d) &\cdots& -\sigma_e(\vx) \vec{E}(\vx) \end{eqnarray*} $ (e) $ を整理するにあたって、位相関数はHenyey-Greensteinを仮定します。Henyey-Greenstein位相関数は次の形で表されます。 \begin{eqnarray*} f_p(\vx, \vomega', \vomega) = \frac{1 - g^2}{4\pi(1 + g^2 - 2g \vomega \cdot \vomega')^{1.5}} \end{eqnarray*} これを $ (e) $ に代入すると次のようになります。 \begin{eqnarray*} \frac{3 \sigma_s(\vx)}{4\pi} \int_{\mathcal{S}^2} \vomega \underbrace{\int_{\mathcal{S}^2} \frac{1 - g^2}{4\pi(1 + g^2 - 2g \vomega \cdot \vomega')^{1.5}} \Prtb{\vec{E}(\vx) \cdot \vomega'} \d\vomega'}_{(e')} \d\vomega \end{eqnarray*} まずは内側の積分 $ (e') $ を整理することを考えます。が、$ \vomega $ は外側の積分変数であり内側の積分を解く段階においては具体的な値が無いという問題があります。しかし $ \vomega \cdot \vomega' = \cos\theta $ ($ \theta $ は2つの方向ベクトルがなす角)という事に着目し、$ \vomega $ を $ \theta = 0 $ の極方向と捉えることによって内側の積分を解くことができます。幸い積分範囲に関しては全方位であるため、球面座標の方向に依存せず特に変換する必要はありません。もう一つの重要な注意点として $ \vec{E}(\vx) $ の座標系があります。今、内側の積分を解くにあたって球面座標の極方向を $ \vomega $ としているため $ \vec{E}(\vx) $ もそれに対応した座標系へと変換する必要があります。$ \vomega $ を極方向とした座標系は、$ \vomega_u, \vomega_v, \vomega $ の3つの正規直交基底ベクトルを使って表すことができます。また、$ \vec{E}(\vx) $ をその座標系で表したものを $ \vec{E}_{\vomega}(\vx) $ とすると2つは次の関係を持ちます。 \begin{eqnarray*} \vec{E}_{\vomega}(\vx) = \Prt{E_{\vomega x}(\vx),\; E_{\vomega y}(\vx),\; E_{\vomega z}(\vx)} = \Prt{\vec{E}(\vx) \cdot \vomega_u,\; \vec{E}(\vx) \cdot \vomega_v,\; \vec{E}(\vx) \cdot \vomega} \end{eqnarray*} したがって内側の積分 $ (e') $ は次のように表せます。 \begin{eqnarray*} \int_{0}^{2\pi} \int_{0}^{\pi} \frac{1 - g^2}{4\pi(1 + g^2 - 2g \cos\theta)^{1.5}} \Prtb{E_{\vomega x}(\vx) \cos\phi \sin\theta + E_{\vomega y}(\vx) \sin\phi \sin\theta + E_{\vomega z}(\vx) \cos\theta} \sin\theta \d\theta \d\phi \end{eqnarray*} これを解くと \begin{eqnarray*} (e') \cdots g \vec{E}(\vx) \cdot \vomega \end{eqnarray*} が得られ、$ (e) $ の外側の積分に代入して解けば、最終的に $ (e) $ の結果として次を得ます。 \begin{eqnarray*} (e) \cdots g \sigma_s(\vx) \vec{E}(\vx) \end{eqnarray*} $ (f) $ は今のところ $ Q $ の具体的な関数を考えていないのでこれ以上整理はできませんが、前述の光源項の1次の項 $ Q_1 $ として表せます。以上、各項の整理により次の式を得ます。 \begin{equation} \nabla \phi(\vx) = - 3 \sigma_e'(\vx) \vec{E}(\vx) + 3 \vec{Q}_1(\vx) \label{eq:gradient_of_fluence} \end{equation} ここで $ \sigma_e'(\vx) $ はreduced extinction coefficientであり、次のように定義されます。 \begin{eqnarray*} \sigma_e'(\vx) = \sigma_s'(\vx) + \sigma_a(\vx) \\ \sigma_s'(\vx) = (1 - g)\sigma_s(\vx) \end{eqnarray*} $ \sigma_s'(\vx) $ はreduced scattering coefficient と呼ばれ、位相関数の平均コサイン $ g $ を使って元の散乱係数 $ \sigma_s(\vx) $ をスケールすることで得られます。媒質を等方的な位相関数とこれらの修正された散乱パラメターを使って光の分布を大きく変えることなくモデル化し直すことができます。
式 $ \eqref{eq:divergence_of_vector_irradiance} $ に式 $ \eqref{eq:gradient_of_fluence} $ を代入することで次の式を得ます。 \begin{equation*} \frac{1}{3 \sigma_e'(\vx)} \nabla^2 \phi(\vx) = \sigma_a(\vx) \phi(\vx) - Q_0(\vx) \label{eq:diffusion_equation} + \frac{1}{\sigma_e'(\vx)} \nabla \cdot \vec{Q}_1(\vx) \end{equation*} 今位相関数は等方的であると考えられるため、$ \vec{Q}_1(\vx) = 0 $ となり右辺第三項は消え、散乱パラメター $ \sigma_e'(\vx) $, $ \sigma_a(\vx) $ が定数とみなせばdipoleモデルなどで使用する、次に示す拡散方程式が得られます。 \begin{equation*} D \nabla^2 \phi(\vx) = \sigma_a \phi(\vx) - Q_0(\vx) \end{equation*} ここで、$ D = 1 \;/\; 3\sigma_e' $ は拡散定数です。
BSSRDFでは物体表面における光の分布に興味を持っているため、それを考慮に入れた境界条件の下、拡散方程式を解く必要があります。物体表面の位置 $ \vx_s $ では、両側の屈折率が異なる界面の場合、反射によりいくらか光のエネルギーが媒質内に戻ってきます。式で表すと次のようになります。
\begin{equation*}
\int_{\mathcal{H}_-} L(\vx_s, \vomega) (\vomega \cdot \vec{n}_-) \d\vomega = F_{dr, -} \int_{\mathcal{H}_+} L(\vx_s, \vomega) (\vomega \cdot \vec{n}_+) \d\vomega
\end{equation*}
ここで、積分範囲 $ \mathcal{H}_+ $, $ \mathcal{H}_- $ はそれぞれ媒質外側、下側の半球を表しています。 $ \vec{n}_+ $, $ \vec{n}_- $ は位置 $ \vx_s $ における法線を表しており、それぞれ外側向き、内側向きに対応します。$ F_{dr, -} $ は物体表面「内側」のdiffuse Fresnel reflectanceで次のように定義されます。
\begin{equation*}
F_{dr, -} = \frac{1}{\pi} \int_{\mathcal{H}_-} F_r(\eta, \vomega') \abs{\vomega' \cdot \vec{n}} \d\vomega'
\end{equation*}
左辺の積分部は位置 $ \vx_s $ から媒質内側へと向かう放射束、右辺の積分部は逆に、位置 $ \vx_s $ において媒質外側へと向かう放射束を表しています。(どちらの放射束もあくまで媒質内部の放射束を表しているため、界面を通過した外側の放射束ではありません。)簡単に言えば、媒質から出ようとする放射束のうちdiffuse Fresnel reflectance分は反射して返ってくるということです。 本来であれば、フレネル反射率は方向に依存するため、diffuse Fresnel reflectanceではなく、積分内部でフレネル反射率をかけるべきではありますが、簡単化のためにこうなっているのでしょう。
ここでも放射輝度を1次までの球面調和関数で近似すれば次の式を得ます。
\begin{equation*}
\phi(\vx_s) + 2D (\vec{n} \cdot \nabla) \phi(\vx_s) = F_{dr, -} \Brk{\phi(\vx_s) - 2D (\vec{n} \cdot \nabla) \phi(\vx_s)}
\end{equation*}
この式を整理すれば、以下に示す境界条件を得ます。境界における解と解の微分を組み合わせた式の値を指定されたロビン境界条件(第三種境界条件)と呼ばれるものになります。
\begin{equation*}
\phi(\vx_s) + 2AD (\vec{n} \cdot \nabla) \phi(\vx_s) = 0 \label{eq:boundary_condition}
\end{equation*}
\begin{equation*}
A = \frac{1 + F_{dr, -}}{1 - F_{dr, -}}
\end{equation*}
この式は物体表面から $ 2AD $ の位置においてfluenceがゼロになることを表しています。
位相関数が等方的な場合 $ \vec{Q}_1(\vx) $ が消えるため、式 $ \eqref{eq:gradient_of_fluence} $ は次のように変形できます。 \begin{equation*} \vec{E}(\vx) = - D \nabla \phi(\vx) \end{equation*} この式にはマイナスの符号がついていることから、エネルギー密度(fluence)の高い領域から低い領域へとエネルギーが流れることを示唆しています。この式と任意のベクトルの内積をとることで位置 $ \vx $ における放射発散度を知ることができます。BSSRDFでは物体表面の法線ベクトル $ \vec{n} $ と内積がとられるでしょうから次のようになります。 \begin{equation*} E(\vx) = - D (\vec{n} \cdot \nabla) \phi(\vx) \end{equation*} ここでは $ \phi(\vx) $ は"完全な"fluenceですが、例えばある点 $ \vx_i $ に入射した放射束 $ \Phi_i(\vx_i) $ によって生じたfluenceによる放射発散度を知るためには $ \phi(\vx) $ を $ \Phi_i(\vx_i) $ で微分すればよいでしょう。つまり次の式のようになります。 \begin{equation*} R(\vx_i, \vx) = - D \frac{\d (\vec{n} \cdot \nabla) \phi(\vx)}{\d\Phi_i(\vx_i)} \end{equation*} これをdiffuse reflectanceと呼びます。