このページではMathJaxによる数式表示を行っています。対応しているブラウザーであれば少し待てば数式が表示されます。
(画面左下にMathJaxのタイプセットステータスが表示されます。)
環境によっては表示が崩れているかもしれません。その際はブラウザーを変更してみてください。

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

このページはレイトレ Advent Calendar 2017のの記事として制作しました。
遅れましたが人口密度が低いので問題無いでしょう...

BRDFは物体表面における光の散乱を記述するモデルで、多くの質感を表現できますが、人間の肌を始めとする柔らかな印象を持つ質感はうまく表現することができません。このページではそれらの質感をBRDFよりも高いレベルかつ効率的に再現できるモデルであるBSSRDFについて紹介します。

表面下散乱 (Sub-Surface Scattering)とBSSRDF

表面下散乱無し
(a) 表面下散乱無し
表面下散乱
(b) 表面下散乱
図1. 表面下散乱の有無の違い
多くの材質では(b)のように(巨視的な)物体表面下において光の伝達が起こる。

BRDFでは光の入射位置と出射位置の距離が十分に近い(図1(a))という仮定を用います。この仮定は金属に対しては有効ですが、半透明(translucent)な材質では入射した光が出射するまでに表面下で長い距離を移動する(図1(b))ので無効です。また、ぱっと見て半透明に見えない材質でもBRDFを用いてレンダリングするとどこか硬い、いかにもCGといった質感になりがちです。物体に表面から入射した光が、内部で散乱して異なる点から出射する現象を表面下散乱 (Sub-Surface Scattering, SSS)と呼びます。表面下散乱を起こす材質例としては、大理石、布、紙、皮膚、牛乳、チーズ、パン、肉、果物、植物、魚、海水、雪など挙げればきりがないほどあります。表面下散乱を考慮した手法でないとこれらの質感をうまく表現することができません。

表面下散乱を正確に再現する方法として一つ考えられるのが、ボリュームレンダリングをモンテカルロ積分を使って真面目に解くことです。しかし人間の皮膚や牛乳などを筆頭に多くの表面下散乱物質は媒質の密度が高く、また散乱アルベドもかなり高いため、光が物質に入射してから出射するまでに非常に多くの散乱イベント(優に数百回〜数千回)が発生します。そのため、一般的にはこの手法は非常に重たくなります。表面下散乱物質は明確な物体表面を持ち、また身の回りの表面下散乱物体の媒質の内部に入って観察するような状況も無いため、実際に興味があるのはBRDFと同じく物体表面における光の挙動だけです。そこで、媒質内の複雑な散乱イベントをブラックボックスとして考え、物体表面位置と方向をパラメターとした反射モデルがBSSRDF (Bidirectional Scattering Surface Reflectance Distribution Function, 双方向散乱面反射率分布関数)です。CG界ではフルスペルが相当長い専門用語として有名?かもしれません(日本語も長いし、時には双方向深層散乱面反射率分布関数というさらに長い記述も見受けられる)。
このBSSRDFの内部の定義を簡単なモデルで記述し、本来の物体表面の光の分布を近似できれば、大量の散乱イベントを真面目に計算しなくとも、良い結果が得られるだろうという狙いでCGの世界に導入されました[Jensen2001]。

基礎理論

BSSRDF
図2. BSSRDF
BSSRDFは物体表面の入射位置、方向、出射位置、方向をパラメターとする8次元関数となる。表面下の個々の光の経路は考えず、2つの位置・方向間の関係を与える。

BSSRDFは物体表面のある位置 $ \vx_i $ に $ \vomega_i $ から入射した微小放射束 $ d\Phi(\vx_i, \vomega_i) $ と、それによって生じる別の位置 $ \vx_o $ から $ \vomega_o $ へと出射される微小放射輝度 $ dL_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{dL_o(\vx_o, \vomega_o)}{d\Phi_i(\vx_i, \vomega_i)} = \frac{dL_o(\vx_o, \vomega_o)}{L_i(\vx_i, \vomega_i) \abs{\vomega_i \cdot \vec{n}_i} d\vomega_i dA(\vx_i)} \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 dA(\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 $ に関する微小放射照度 $ dE_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 dA(\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*}

散乱光とBSSRDFの分解

単散乱
(a) 単散乱
多重散乱
(b) 多重散乱
図3. 散乱回数による分類

表面下散乱して物体表面から出射する光は散乱回数に応じて分類することができます。表面下で1回だけ散乱してから物体表面に戻って出ていくことを単散乱(single scattering) (図3(a))、2回の散乱なら二重散乱(double scattering) (図3(b))という具合にわけられます。単散乱して出てくる光は、物体表面の粗さにもよりますが高い指向性を持っており、2回以上散乱した光は指向性を失うため、2回以上の散乱を多重散乱(multiple scattering)として単散乱と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)の形で表すことができます。前述のdiffuse reflectanceや拡散方程式の導出などの詳細については「表面下散乱の拡散近似」で解説したいと思います。

\begin{equation} D \nabla^2 \phi(\vx) = \sigma_a \phi(\vx) - Q_0(\vx) \label{eq:diffusion_equation} \end{equation}

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

ダイポールモデル
図4. ダイポールモデル
ダイポールモデルでは、2つの点光源からの寄与の和でfluenceを表現する。

拡散方程式 $ \eqref{eq:diffusion_equation} $ は無限に拡がる媒質中に点光源があるケースでは解析的な解が求まります。しかし興味があるのは物体表面の光の分布なので、物体表面におけるエネルギーの収支を境界条件として拡散方程式を解きたいところです。任意の物体形状に対して解析解を出すのは困難なので、物体形状が半無限の板であるという仮定をまず使います。一応、この仮定の下では正確な解があるみたいなのですが、ベッセル関数の無限和を使っていたりと使い勝手が悪く、ダイポール(dipole)モデルと呼ばれる近似が開発されます。ダイポールモデルでは名前の通り、2つの点光源を境界条件を満たすように配置し、それらの寄与の和によってfluenceの分布が表現されます。式としては次の形になります。

\begin{equation} \phi(\vx) = \frac{\Phi}{4 \pi D} \Prt{\frac{\mathrm{e}^{-\sigma_{tr} d_r}}{d_r} - \frac{\mathrm{e}^{-\sigma_{tr} d_v}}{d_v}} \label{eq:fluence_dipole} \end{equation}

物体表面より下 $ z_r $ の距離に正の放射束を持つ光源(real source)が、物体表面より上 $ z_v = z_r + 4AD $ の位置に負の放射束を持つ光源(virtual source)が配置されます(図4)。 ここで、$ d_r = \norm{\vx - \vx_r} $ は正の光源と評価点 $ \vx $ 間の距離、$ d_v = \norm{\vx - \vx_v} $ は負の光源と評価点 $ \vx $ 間の距離です。$ \sigma_{tr} = \sqrt{3 \sigma_a \sigma_e'} $ はeffective transport coefficient (実効輸送係数?)と呼ばれます。正の光源の深さは $ z_r = 1 \;/\; \sigma_e' $ とすることが提案されています。これは表面下における平均自由行程になります。$ A $ は表面下の光が物体表面裏側で、フレネル反射によりいくらか反射して戻ってくるという境界条件を解くと求まる数値で次のように表されます。

\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*}

ここで、$ \eta $ は相対屈折率(ですが、ここでは外側に対する内側の屈折率として考えられていることに注意)です。次の多項式近似を使うことも可能です。

\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*}

ダイポールによるfluenceの式 $ \eqref{eq:fluence_dipole} $ をdiffuse reflectanceの式 $ \eqref{eq:diffuse_reflectance} $ に代入すると以下の式を得ます。(1本の光線として入射する放射束によってのみ生じるfluenceを考えているので、恐らく以下の式中の $ d\phi(\vx_o) \;/\; d\Phi_i(\vx_i) $ 部分に 式 $ \eqref{eq:fluence_dipole} $ を $ \Phi_i(\vx_i) $ で割った形が対応します。)

\begin{eqnarray*} R_d(\vx_i, \vx_o) &=& R(r) = -D \frac{d (\vec{n}_o \cdot \nabla \phi)(\vx_o)}{d\Phi_i(\vx_i)} \\ &=& \frac{\Lambda'}{4 \pi} \Brk{z_r (1 + \sigma_{tr} d_r) \frac{\mathrm{e}^{-\sigma_{tr} d_r}}{d_r^3} + z_v (1 + \sigma_{tr} d_v) \frac{\mathrm{e}^{-\sigma_{tr} d_v}}{d_v^3}} \end{eqnarray*}

ここで $ \Lambda' = \sigma_s' \;/\; \sigma_e' $ はreduced scattering albedoです。diffuse reflectanceは一般的には $ \vx_i $ と $ \vx_o $ それぞれの位置の関数 $ R(\vx, \vomega) $ でしたが、(無限平面の仮定を使っているため)結果としてはそれらの間の距離 $ r $ の関数 $ R(r) $ となります。この距離の関数をdiffusion profile (拡散プロファイル)やdiffusion kernel (拡散カーネル)と呼ぶことがあります。
diffuse reflectanceと物体表面のフレネル項を併せてBSSRDFの多重散乱項を次のように構成します。

\begin{equation} S_d(\vx_i, \vomega_i, \vx_o, \vomega_o) = \frac{1}{\pi} F_t(\eta, \vomega_i) R(\norm{\vx_i - \vx_o}) F_t(\eta, \vomega_o) \label{eq:BSSRDF_diffusion_term} \end{equation}

この式の解釈は次のようになります。
1. 物体表面 $ \vx_i $ に $ \vomega_i $ から入射した放射束がフレネル反射されて残った分が透過する。
2. diffuse reflectanceにより $ \vx_o $ の内側の放射照度へ変換される。
3. $ \pi $ で割ることにより $ \vx_o $ 内側に一様に入射する放射輝度となる。
4. フレネル反射して残った分が $ \vomega_o $ へと出射する放射輝度となる。

BSSRDFの単散乱項

BSSRDFの単散乱項
図5. BSSRDFの単散乱項

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

粗い表面への一般化

ここまでは、物体表面の材質は滑らか = 完全スペキュラーとしてBSSRDFを記述しましたが、任意のBSDFを使う場合の理論の拡張が[Donner2006]で述べられています。多重散乱項のフレネル透過を表す項の代わりにBSDFの積分で置き換えます。

\begin{eqnarray*} S_d(\vx_i, \vomega_i, \vx_o, \vomega_o) = \int_{\mathcal{H}_-^2} f_s(\vx_i, \vomega_i, \vomega'_i) \abs{\vomega'_i \cdot \vec{n}_i} d\vomega'_i \frac{1}{\pi} R(\norm{\vx_i - \vx_o}) \int_{\mathcal{H}_-^2} f_s(\vx_o, \vomega'_o, \vomega_o) \abs{\vomega'_o \cdot \vec{n}_o} d\vomega'_o \end{eqnarray*}

ここで、$ \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 $ も次に示すようにhemispherical-hemispherical reflectanceに置き換える必要があります。

\begin{equation*} A = \frac{1 + \rho_d}{1 - \rho_d} \end{equation*} \begin{equation*} \rho_d = \frac{1}{\pi} \int_{\mathcal{H}_-} \int_{\mathcal{H}_-} f_s(\vx, \vomega', \vomega) \abs{\vomega \cdot \vec{n}} \abs{\vomega' \cdot \vec{n}} d\vomega d\vomega' \end{equation*}

さいごに

このページでは表面下散乱物質の質感のモデル化に有用なBSSRDF、それを効率的に計算することを可能にするダイポール近似モデルについて紹介しました。CGの分野には、他の研究分野から流用された技術が多くありますが、BSSRDFで使われている考え方、特に拡散方程式やその解き方などは他の分野においても重要性が高く歴史的な厚み・重みを感じます。
ダイポール近似はBSSRDFとしては初期のモデルなので、物体形状や散乱パラメターに制約が強く、制約を破るとアーティファクトが出ます。制約を取り払うため、multipoleやdirectional dipole, photon beam diffusion, empiricalモデル等、数多くのモデルが開発されていますが、根底にある考え方は同様です。

参考文献

  1. [Donner2006] Craig Steven Donner - "Towards Realistic Image Synthesis of Scattering Materials", 2006
  2. [Habel2013] Ralf Habel, Per H. Christensen, Wojciech Jarosz - "Classical and Improved Diffusion Theory for Subsurface Scattering", 2013
  3. [Jensen2001] Henrik Wann Jensen, Stephen R. Marschner, Marc Levoy, Pat Hanrahan - "A Practical Model for Subsurface Light Transport", 2001
  4. [Pharr2016] Matt Pharr, Wenzel Jakob, Greg Humphreys - "PHYSICALLY BASED RENDERING: From Theory to Implementation Third Edition", 2016, Morgan Kaufmann

Today : 00000013 Total : 00065179
Copyright © 2017 shocker.