双方向散乱面反射率分布関数 (BSSRDF)

BSSRDFは物体表面のある位置 $ \vx_i $ に $ \vomega_i $ から入射した微小放射束 $ \d\Phi(\vx_i, \vomega_i) $ と、それによって生じる別の位置 $ \vx_o $ から $ \vomega_o $ へと出射される微小放射輝度 $ \d L_o(\vx_o, \vomega_o) $ の関係を与えます(図2)。 一般的に文字 $ S $ で表され、数式で書くと次のような8次元関数となります。単位は $ \Brk{\mathrm{m^{-2} \cdot {sr}^{-1}}} $ です。 \begin{equation*} S(\vx_i, \vomega_i, \vx_o, \vomega_o) = \frac{\d L_o(\vx_o, \vomega_o)}{\d\Phi_i(\vx_i, \vomega_i)} = \frac{\d L_o(\vx_o, \vomega_o)}{L_i(\vx_i, \vomega_i) \abs{\vomega_i \cdot \vec{n}_i} \d\vomega_i \d A(\vx_i)} \hspace{5mm} \Brk{\mathrm{m^{-2} \cdot {sr}^{-1}}} \end{equation*} ここで考えている微小放射束は立体角(方向)と面積(位置)に関する二重の微小量なのでイメージ的には、分子は $ \d^2 L_o(\vx_o, \vomega_o) $、分母は $ \d\Phi_i^2(\vx_i, \vomega_i) $ と書いたほうが良いかもしれません(二乗じゃないですよ)。

BSSRDFを考えている物体表面からの反射光は次のように計算されます。 \begin{equation*} L_o(\vx_o, \vomega_o) = \int_{\mathcal{M}} \int_{\mathcal{H}_+^2} S(\vx_i, \vomega_i, \vx_o, \vomega_o) L_i(\vx_i, \vomega_i) \abs{\vomega_i \cdot \vec{n}_i} \d\vomega_i \d A(\vx_i) \end{equation*} ここで、$ \vec{n}_i $ は入射点 $ \vx_i $ における法線、$ \mathcal{M} $ と $ \mathcal{H}_+^2 $ はそれぞれ物体表面と半球の積分範囲(外側)を表しています。
BRDFによる反射光の計算式を比較のために書いてみます。 \begin{equation*} L_o(\vx, \vomega_o) = \int_{\mathcal{H}_+^2} f_r(\vx, \vomega_i, \vomega_o) L_i(\vx, \vomega_i) \abs{\vomega_i \cdot \vec{n}} \d\vomega_i \end{equation*} BRDFの場合、併せて考える物理量はある入射方向 $ \vomega_i $ に関する微小放射照度 $ \d E_i(\vx, \vomega_i) = L_i(\vx, \vomega_i) \abs{\vomega_i \cdot \vec{n}} \d\vomega_i $ で、入射方向に関して積分することで反射する放射輝度が求まりました。BSSRDFの場合は、併せて考える物理量はある入射方向 $ \vomega_i $、さらにある入射位置 $ \vx_i $ に関する微小放射束 $ \d\Phi_i(\vx_i, \vomega_i) = L_i(\vx_i, \vomega_i) \abs{\vomega_i \cdot \vec{n}_i} \d\vomega_i \d A(\vx_i) $ なので、入射方向と入射位置に関して積分することで反射する放射輝度が求まります。
またBSSRDFはBRDFと同じく相反性を満たします(実装が従っているかは別...)。つまり次のように入射と出射のペアを入れ替えても同じ値となります。 \begin{equation*} S(\vx_i, \vomega_i, \vx_o, \vomega_o) = S(\vx_o, \vomega_o, \vx_i, \vomega_i) \end{equation*}

表面下散乱して物体表面から出射する光は散乱回数に応じて分類することができます。表面下で1回だけ散乱してから物体表面に戻って出ていくことを単散乱(single scattering) (図3(a))、2回の散乱なら二重散乱(double scattering)という具合にわけられます。単散乱して出てくる光は、物体表面の粗さにもよりますが高い指向性を持っており、2回以上散乱した光は指向性を失うため、2回以上の散乱を多重散乱(multiple scattering) (図3(b))として単散乱と2つにわけて扱い、BSSRDFのモデルとしては次のように単散乱と多重散乱を分けて扱うことがよく行われます。

\begin{equation*} S(\vx_i, \vomega_i, \vx_o, \vomega_o) = S^{(1)}(\vx_i, \vomega_i, \vx_o, \vomega_o) + S_d(\vx_i, \vomega_i, \vx_o, \vomega_o) \end{equation*}

散乱性の非常に高い物質内では、例え入射光の分布や位相関数形状が指向性を持っていた(anisotropic)としても、何度も散乱しているうちに、光の方向に関する分布は一様になっていきます。そのため、媒質を等方的な位相関数と適切に修正された散乱パラメターを使って光の分布を大きく変えることなくモデル化し直すことができます。これをprinciple of similarityと呼びます。修正された散乱パラメターは次のように定義されます。 \begin{equation*} \sigma_e' = \sigma_s' + \sigma_a \hspace{5mm} \sigma_s' = \sigma_s (1 - g) \end{equation*} $ \sigma_e' $, $ \sigma_s' $ はそれぞれreduced extinction/scattering coefficientと呼ばれ、通常の係数に位相関数の方向性を表すパラメター $ g $ で補正をかける形で表されます。$ g $ は位相関数の散乱方向の平均コサインであり、次のように定義されます。 \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 $ の場合は完全等方散乱です。直感的な解釈としては、前方散乱が強くなるほど($ g \RAR 1 $)、散乱しない($ \sigma_s' \RAR 0 $)場合と区別がつかなくなることを表しています。逆に、後方散乱が強くなれば($ g \RAR -1 $)、光が同じ方向に続けて進みづらくなるので係数が大きくなるようにはたらきます。このような解釈から"reduced"の訳を考えるならば「換算」とかになるのかな、というのが個人的な意見です(力学で出てくるreduced massの和訳も換算質量)。

光の方向に関する分布が一様な状況では、放射輝度分布を少数の球面調和関数を用いて近似することができます。球面調和関数近似と位相関数の等方性の仮定より、次に示すBSSRDFのdiffuse reflectanceの式が導かれます。 \begin{equation} R_d(\vx_i, \vx_o) = -D \frac{\d (\vec{n}_o \cdot \nabla \phi)(\vx_o)}{\d\Phi_i(\vx_i)} \hspace{5mm} \Brk{\mathrm{m^{-2}}} \label{eq:diffuse_reflectance} \end{equation} ここで、$ D = 1 \;/\; 3 \sigma_e' $ は拡散定数(diffusion constant)、$ \phi(\vx) $ はradiant fluence (scalar irradiance)と呼ばれ、光の位置に関する分布を表しています。$ R_d $ は位置 $ \vx_i $ に入射した(界面通過後の)放射束と、それによって位置 $ \vx_o $ の外側に法線方向 $ \vec{n}_o $ に生じる放射発散度の関係を表します。入射光によって生じる(微小)fluenceの式さえ用意できれば出射点の放射発散度がわかるので出射する放射輝度まで関連づけることができます。それではfluenceの式はどのようにして決めるのでしょうか。
fluenceの式は、前述の式と同様に球面調和関数近似と位相関数の等方性の仮定などを用いることで、以下に示す古典的な定常状態の拡散方程式(diffusion equation)の形で表すことができます。 \begin{equation} D \nabla^2 \phi(\vx) = \sigma_a \phi(\vx) - Q_0(\vx) \label{eq:diffusion_equation} \end{equation} また、物体表面位置 $ \vx_s $ におけるエネルギーの収支(簡単に言えば、物体表面内側でフレネル反射によりある程度光が媒質内部に戻ってくる)を考えると次に示す境界条件を得ます。 \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*} ここで、$ F_{dr, -} $ は物体表面「内側」のdiffuse Fresnel reflectanceで次のように定義されます。(元論文[Jensen2001]では恐らく $ 1 \;/\; \pi $ が抜けているようです。$ \pi $ で割ることで後述の多項式近似とも一致する。) \begin{equation*} F_{dr, -} = \frac{1}{\pi} \int_{\mathcal{H}_-} F_r(\eta, \vomega') \abs{\vomega' \cdot \vec{n}} \d\vomega' \end{equation*} 積分範囲 $ \mathcal{H}_- $ は媒質内側の半球、$ \eta $ は相対屈折率(ですが、ここでは外側に対する内側の屈折率として考えられていることに注意)、$ F_r(\eta, \vomega') $ はフレネル反射率です。次の多項式近似を使うことも可能です。 \begin{equation*} F_{dr, -} \approx \begin{cases} \displaystyle -0.4399 + \frac{0.7099}{\eta} - \frac{0.3319}{\eta^2} + \frac{0.0636}{\eta^3} & \eta < 1 \\ \displaystyle -\frac{1.4399}{\eta^2} + \frac{0.7099}{\eta} + 0.6681 + 0.0636 \eta & \eta > 1 \\ \end{cases} \end{equation*} diffuse reflectance (式 $ \eqref{eq:diffuse_reflectance} $) や拡散方程式 $ \eqref{eq:diffusion_equation} $、境界条件 $ \eqref{eq:boundary_condition} $ の導出などの詳細については「表面下散乱の拡散近似」で解説したいと思います。

ダイポールモデル (dipole model)

diffuse reflectanceと物体表面のフレネル項を併せてBSSRDFの多重散乱項を次のように構成します。 \begin{equation*} S_d(\vx_i, \vomega_i, \vx_o, \vomega_o) = \frac{1}{\pi} F_t(\eta, \vomega_i) R_d(\norm{\vx_i - \vx_o}) F_t(\eta, \vomega_o) \end{equation*} ここで、$ F_t(\eta, \vomega) = 1 - F_r(\eta, \vomega) $ はフレネル透過率です。 この式の解釈は次のようになります。
1. 物体表面 $ \vx_i $ に $ \vomega_i $ から入射した放射束がフレネル反射されて残った分が透過する。
2. diffuse reflectanceにより $ \vx_o $ の外側の放射発散度へ変換される。
3. $ \pi $ で割ることにより $ \vx_o $ において外側に一様に出射する放射輝度となる。
4. フレネル項をかけることで出射する放射輝度がフレネル反射(透過)を考慮した分布となる。

正規化係数の導入

単散乱項は特に近似計算は行わず明示的に散乱イベントを考えます。単散乱の光の経路を表面下で構築するためには、入射するレイ、出射するレイそれぞれが透過・屈折したレイが交わる必要があります。そのため物体表面が滑らか、完全スペキュラーとして考えている場合は、任意に選んだ入射・出射レイの組み合わせに対して寄与を持つことは無く、単散乱項は暗にデルタ関数を含むことになります。

ここまでは、物体表面の材質は滑らか = 完全スペキュラーとしてBSSRDFを記述しましたが、任意のBSDFを使う場合の理論の拡張が[Donner2006]で述べられています。多重散乱項のフレネル透過を表す項の代わりにBSDFの積分 = directional-hemispherical transmittance $ \rho_t(\vomega) $で置き換えます。また、正規化係数においても該当の項を同様にhemispherical-hemispherical transmittance $ \rho_{dt} $ で置き換えます。 \begin{eqnarray*} S_d(\vx_i, \vomega_i, \vx_o, \vomega_o) &=& \rho_t(\vomega_i) \; R_d(\norm{\vx_i - \vx_o}) \; C \; \rho_t(\vomega_o) \\ &=& \int_{\mathcal{H}_-^2} f_s(\vx_i, \vomega_i, \vomega'_i) \absb{\vomega'_i \cdot \vec{n}_i} \d\vomega'_i \; R_d(\norm{\vx_i - \vx_o}) \; C \int_{\mathcal{H}_-^2} f_s(\vx_o, \vomega_o, \vomega'_o) \absb{\vomega'_o \cdot \vec{n}_o} \d\vomega'_o \end{eqnarray*} \begin{equation*} C = \frac{1}{\pi \rho_{dt}} = \frac{1}{\int_{\mathcal{H}_+^2} \int_{\mathcal{H}_-^2} f_s(\vx_o, \vomega_o, \vomega'_o) \absb{\vomega'_o \cdot \vec{n}_o} \absb{\vomega_o \cdot \vec{n}_o} \d\vomega'_o \d\vomega_o} \end{equation*} ここで、$ \vomega'_i $, $ \vomega'_o $ は位置 $ \vx_i $, $ \vx_o $ における内側の方向を表しています。$ \vec{n}_i $, $ \vec{n}_o $ も同様にそれぞれの位置における法線を表しています。BSDFを積分することで、放射束としての透過率を表しています。完全スペキュラーBSDFを代入すると元の式 $ \eqref{eq:BSSRDF_diffusion_term} $ と等しくなります(デルタ関数の扱いに注意)。
また、diffuse Fresnel reflectanceから計算されていた数値 $ A $ も次に示すように内側のbihemispherical reflectance $ \rho_{dr, -} $ に置き換える必要があります。 \begin{equation*} A = \frac{1 + \rho_{dr, -}}{1 - \rho_{dr, -}} \end{equation*} \begin{equation*} \rho_{dr, -} = \frac{1}{\pi} \int_{\mathcal{H}_-} \int_{\mathcal{H}_-} f_s(\vx, \vomega', \vomega) \absb{\vomega \cdot \vec{n}} \absb{\vomega' \cdot \vec{n}} \d\vomega \d\vomega' \end{equation*}