このページは レイトレ Advent Calendar 2024 の記事として制作しました
水や肌など、絶縁体の媒質を伝播する光の近似的な振る舞いは Volume Rendering 方程式によって記述されます
[Arvo 1993][Cercignani 1988]。
媒質の中を光が多数回散乱するとき、その散乱方向はほぼ等方的になります。
このとき、Volume Rendering 方程式に対する近似として拡散近似が用いられます
[Jensen et al. 2001]
。
拡散近似のもとでは、Volume Rendering 方程式は単純な拡散方程式に帰着し、その解は畳み込み積分で表されます。
ゲームなどリアルタイムレンダリングで用いられる subsurface scattering では、
diffusion profile を用いたガウシアンの畳み込み積分がしばしば実装されます
[d’Eon and Luebke. 2007]
[Jimenez et al. 2015]。
この記事では Volume Rendering 方程式から出発し、 subsurface scattering に用いられるガウシアンの畳み込み積分に至る道筋を眺めます。
静止した媒質に光が入射する現象を考えます。光を単純な点粒子として扱う近似のもとでは、3次元空間上の位置 $\bm{r}$ と方向 $\bm{\omega}$ に対応する放射輝度 $L(\bm{x},\bm{\omega})$ の変化は、次式で記述されます [Arvo 1993][Cercignani 1988]。
$$ \bm{\omega}\cdot\nabla L(\bm{r},\bm{\omega}) = \mu_s\int_{\mathcal{S}^2}p(\bm{\omega}\cdot\bm{\omega}')L(\bm{r},\bm{\omega}')d\bm{\omega}' -\mu_tL(\bm{r},\bm{\omega}) + S(\bm{r},\bm{\omega}) $$
この式は Volume Rendering 方程式、あるいは Radiative Transfer 方程式と呼ばれています。
ここでは、$\mu_a$ と $\mu_s$ は一定、かつ位相関数 $p$ は軸対称、 つまり $\bm{\omega}$ と $\bm{\omega}'$ のなす角度のみに依存すると仮定します。 位相関数 $p$ は 1 に規格化されているとし
$$ \int_{\mathcal{S}^2}p(\bm{\omega}\cdot\bm{\omega}')d\bm{\omega}'=1 $$
異方性因子 $g$ を$$ g=\int_{\mathcal{S}^2}(\bm{\omega}\cdot\bm{\omega}')\,p(\bm{\omega}\cdot\bm{\omega}')d\bm{\omega} $$
によって定義します。 このような位相関数の例として、Mie散乱の近似として有名な Henyey-Greenstein 関数があります。
輝度 $L(\bm{r},\bm{\omega})$ を球面調和関数で展開すると、以下のように書けます
$$ L(\bm{r},\bm{\omega}) = \frac{1}{4\pi}\mathcal{M}_0[L(\bm{r},\bm{\omega})] +\frac{3}{4\pi}\bm{\omega}\cdot\mathcal{M}_1[L(\bm{r},\bm{\omega})] +\cdots $$
ここで、$\mathcal{M}_l$ は「$l$ 次のモーメント」と呼ばれる量を表します [Pharr et al. 2016]
$$ \mathcal{M}_l\left[f(\bm{\omega})\right] := \int_{\mathcal{S}^2}\underbrace{\bm{\omega}\otimes\bm{\omega}\otimes\cdots\otimes\bm{\omega}}_{l\,\text{factors}}\,f(\bm{\omega})d\bm{\omega} $$
輝度の0次、1次、2次のモーメントは、それぞれ照度の次元をもち、 以下のように定義されます [Jensen et al. 2001]
$$ \bm{E}_1(\bm{r}):=\mathcal{M}_1\left[L(\bm{r},\bm{\omega})\right] = \int_{\mathcal{S}^2}\bm{\omega}L(\bm{r},\bm{\omega})d\bm{\omega} $$
$$ \bm{E}_2(\bm{r}):=\mathcal{M}_2\left[L(\bm{r},\bm{\omega})\right] = \int_{\mathcal{S}^2}(\bm{\omega}\otimes\bm{\omega})L(\bm{r},\bm{\omega})d\bm{\omega} $$
同様に、光源項 $S$ のモーメントを
$$ Q(\bm{r}) :=\mathcal{M}_0[S(\bm{r},\bm{\omega})] =\int_{\mathcal{S}^2}S(\bm{r},\bm{\omega})d\bm{\omega} ,\quad Q_1(\bm{r}) :=\mathcal{M}_1[S(\bm{r},\bm{\omega})] =\int_{\mathcal{S}^2}\bm{\omega}S(\bm{r},\bm{\omega})d\bm{\omega} $$
などと定義します。
Volume Rendering 方程式の0次モーメントから、以下の式が得られます
$$ \int_{\mathcal{S}^2} \bm{\omega}\cdot\nabla L(\bm{r},\bm{\omega})d\bm{\omega} = \int_{\mathcal{S}^2}\left[ -\mu_tL(\bm{r},\bm{\omega}) +\mu_s\int_{\mathcal{S}^2}p(\bm{\omega}\cdot\bm{\omega}')L(\bm{r},\bm{\omega}')d\bm{\omega}' +S(\bm{r},\bm{\omega}) \right]d\bm{\omega} $$
左辺は $\int_{\mathcal{S}^2}\bm{\omega}\cdot\nabla Ld\bm{\omega}=\nabla\cdot\bm{E}_1$ となり、右辺第2項は $\mu_s\int_{\mathcal{S}^2}p(\bm{\omega}\cdot\bm{\omega}')L(\bm{r},\bm{\omega}')d\bm{\omega}' = \mu_sE$、 光源項の球面積分は $Q$ であるため、
$$ \nabla\cdot\bm{E}_1(\bm{r}) = -\mu_tE(\bm{r}) + \mu_sE(\bm{r}) + Q_0(\bm{r}) = -\mu_a(\bm{r})E(\bm{r}) + Q(\bm{r}) $$
これが0次モーメント方程式です。
Volume Rendering 方程式の1次モーメントから、以下の式が得られます
$$ \int_{\mathcal{S}^2} \bm{\omega}\left[\bm{\omega}\cdot\nabla L(\bm{r},\bm{\omega})\right]d\bm{\omega} = \int_{\mathcal{S}^2}\bm{\omega}\left[ -\mu_tL(\bm{r},\bm{\omega}) +\mu_s\int_{\mathcal{S}^2}p(\bm{\omega}\cdot\bm{\omega}')L(\bm{r},\bm{\omega}')d\bm{\omega}' +S(\bm{r},\bm{\omega}) \right]d\bm{\omega} $$
左辺は $\int_{\mathcal{S}^2}\bm{\omega}\left[\bm{\omega}\cdot\nabla L\right]d\bm{\omega}=\nabla\cdot\bm{E}_2$。 右辺第1項は $\int_{\mathcal{S}^2}\bm{\omega}\left(\mu_tL\right)=\mu_t\bm{E}_1$、 第2項は $\mu_s\int_{\mathcal{S}^2}L(\bm{r},\bm{\omega}')\left[\int_{\mathcal{S}^2}\bm{\omega}p(\bm{\omega}\cdot\bm{\omega}')d\bm{\omega}\right]d\bm{\omega}' = \mu_s g \bm{E}_1$、 第3項は $\bm{Q}_1$ です。よって
$$ \nabla\cdot\bm{E}_2(\bm{r}) = -\mu_t\bm{E}_1+\mu_sg\bm{E}_1(\bm{r})+\bm{Q}_1(\bm{r}). $$
$\mu'_s=\mu_s(1-g)$ とすれば、
$$ \nabla\cdot\bm{E}_2(\bm{r}) = -\left(\mu_a+\mu'_s\right)\bm{E}_1(\bm{r}) + \bm{Q}_1(\bm{r}) $$
これが1次モーメント方程式です。 以上の計算の詳細は、例えば [Jakob et al. 2010] の Expanded Technical Report に記載されています。
Volume Rendering方程式から得られる 0次と1次モーメント方程式は
$$ \left\{ \, \begin{aligned} & \nabla\cdot\bm{E}_1(\bm{r})= -\mu_a(\bm{r})E(\bm{r}) + Q(\bm{r}) \\ & \nabla\cdot\bm{E}_2(\bm{r})= -\left(\mu_a+\mu'_s\right)\bm{E}_1(\bm{r}) + \bm{Q}_1(\bm{r}) \end{aligned} \right. $$
です。 ただし、0次には $\bm{E}_1$、1次には $\bm{E}_2$ と高次モーメントが入って来るため、 この方程式の階層は閉じません。ここで拡散近似を導入し、高次のモーメントを打ち切ることを行います。
輝度の球面調和関数展開を1次で打ち切る近似を、拡散近似 [Jensen et al. 2001] あるいは$P1$近似と呼びます
$$ L(\bm{r},\bm{\omega}) = \frac{1}{4\pi}\mathcal{M}_0[L(\bm{r},\bm{\omega})] +\frac{3}{4\pi}\bm{\omega}\cdot\mathcal{M}_1[L(\bm{r},\bm{\omega})] $$
$\bm{E}_2$ の定義式に代入すると
$$ \bm{E}_2(\bm{r}) = \int_{\mathcal{S}^2}\left(\bm{\omega}\otimes\bm{\omega}\right) \left[\frac{1}{4\pi}E(\bm{r}) + \frac{3}{4\pi}\left(\bm{E}_1(\bm{r})\cdot\bm{\omega}\right)\right]d\bm{\omega} $$
$\left(\bm{\omega}\otimes\bm{\omega}\right)$ の球面積分は$\frac{4\pi}{3}I$ ($I$は単位テンソル)、 第2項の積分は対称性からゼロになるため、
$$ \bm{E}_2(\bm{r})=\frac{E(\bm{r})}{3}I. $$
これを1次モーメント方程式 $\nabla\cdot\bm{E}_2= -\left(\mu_a+\mu'_s\right)\bm{E}_1 + \bm{Q}_1$ に代入すると、
$$ \frac{1}{3}\nabla E(\bm{r}) = -\left(\mu_a+\mu'_s\right)\bm{E}_1(\bm{r}) + \bm{Q}_1(\bm{r}) $$
等方的な光源($\bm{Q}_1=0$)を仮定すれば、
$$ \bm{E}_1(\bm{r})=-DE(\bm{r}), \quad D:=\frac{1}{3(\sigma_a+\sigma'_s)} $$
これは $D$ を拡散係数とする、Fickの法則の形になります。
0次モーメント方程式
$$ \nabla\cdot\bm{E}_1(\bm{r}) = -\mu_a E(\bm{r}) + Q(\bm{r}) $$
に対して、拡散近似で得られた $\bm{E}_1=-DE$ を代入すると
$$ \nabla\cdot\left[-D\nabla E(\bm{r})\right] = -\mu_aE(\bm{r}) + Q(\bm{r}), $$
すなわち
$$ D\nabla^2E(\bm{r}) = \mu_aE(\bm{r}) - Q(\bm{r}) $$
これが拡散方程式です。
拡散方程式は、別名 Screened Poisson 方程式としても知られ、以下のように書くことが多いです
$$ (\nabla^2-\mu^2_\text{eff})E(\bm{r})=-\frac{Q(\bm{r})}{D} $$
ただし $\mu^2_\text{eff}:=\frac{\mu_a}{D}=3\mu_a(\mu_a+\mu'_s)$ です。
この方程式は Green 関数 $\mathcal{G}$ を用いて解が与えられることが知られています。 例えば $Q$ を点光源とすれば
$$ E(\bm{r}) = \frac{1}{D}\int Q(\bm{r}')\mathcal{G}(|\bm{r}-\bm{r}'|)d\bm{r}'. $$
Screened Poisson 方程式の Green 関数は Yukawa カーネルと呼ばれ、以下のように与えられます
$$ \mathcal{G}(|\bm{r}-\bm{r}'|) =\frac{1}{4\pi}\frac{e^{-\mu_\text{eff}|\bm{r}-\bm{r}'|}}{|\bm{r}-\bm{r}'|} $$
3次元空間上の2次元表面を $\mathcal{A}$ とし、そこでの座標を $\bm{s}$ と書きます。 Subsurface scattering の手法の一つとして、表面上の照度の初期分布 $E_\text{in}$ と diffusion profile $R_d$ の畳み込み積分によって、表面上の出射照度 $E_\text{out}$ を求めるものがあります。 典型的 [Jimenez et al. 2015] には
$$ E_\text{out}(\bm{s}) = \int_{\mathcal{A}}E_\text{in}(\bm{s}')R_d(|\bm{s}-\bm{s}'|)d\bm{s}'. $$
ここで用いられる $R_d$ は、3次元のグリーン関数 $\mathcal{G}$ を、深さ方向に対して積分し2次元化したものとみなせます [d’Eon and Luebke. 2007]。
拡散近似による照度
$$ E(\bm{r})=\int Q(\bm{r}')\mathcal{G}(|\bm{r}-\bm{r}'|)d\bm{r}' $$
を考え、$\bm{r}=(x,y,z)$ のうち表面 $z=0$ の値に着目します。 [Farrell et al. 1992] と同様に、 媒質内部の光源項 $Q(\bm{r}')$ を表面成分 $E_\text{in}(\bm{s}')$ と深さ方向の減衰成分 $d(z')$ で構成します
$$ Q(\bm{r}') = E_\text{in}(\bm{s}')d(z'). $$
そして $\mathcal{G}$ を深さ方向に対して積分すると、表面上で有効な応答関数
$$ R_d(|\bm{s}-\bm{s}'|) :=\int d(z')\mathcal{G}\left(\sqrt{(x-x')^2+(y-y')^2+(z')^2}\right)dz' $$
が得られます。 この記事では、これを2Dカーネルとしての diffusion profile として定義します。 よって最終的に
$$ E_\text{out}(\bm{s}) =\int E_\text{in}(\bm{s}')R_d(|\bm{s}-\bm{s}'|)d\bm{s}' $$
となります。
一部のリアルタイムレンダリングにおいて、 Yukawaカーネルより取り扱いやすいガウシアンカーネルが用いられることがあります。 例えば
$$ G(\bm{s};\,\nu_i) :=\frac{1}{2\pi\nu_i}\exp\left(-\frac{|\bm{s}|}{2\nu_i}\right) $$
というガウシアンで diffusion profile $R_d$ を近似する手法 (sum of Gaussians)が使われます [d’Eon and Luebke. 2007]。 ガウシアンカーネルは、座標軸ごとに積分を分離できるなど、計算を効率化しやすい性質をもちます [Jimenez et al. 2015] たとえば複数のガウシアンを重ね合わせることで、
$$ R_d(\bm{s}') \sim \sum_iw_iG(\bm{s};\,\nu_i) $$
という表現を用いることができます。このとき、
$$ E_\text{out}(\bm{s}) = \int_{\mathcal{A}}E_\text{in}(\bm{s}')R_d(|\bm{s}-\bm{s}'|)d\bm{s}' \sim \sum_kE_\text{in}(\bm{s}_k)R_d(|\bm{s}-\bm{s}_k|) $$
となります。
以上のように、Volume Rendering 方程式から出発して、 subsurface scattering で使われるガウシアンの畳み込み積分に至ることができました。
この記事では、Radiative Transfer 方程式から 0 次と 1 次モーメント方程式を導き、 拡散近似を通して 拡散方程式 が得られる流れを書きました。 そして拡散方程式のグリーン関数を2次元表面に射影したものを、 subsurface scattering の diffusion profile $R_d$ として定義し、 それがガウシアンを用いた sum-of-Gaussians 表現に繋がることを示しました。 このような、拡散近似で得た畳み込み積分は実装が簡単で物理的な妥当性もあり、 subsurface scattering の第一歩として広く応用されている手法のひとつと言えます。
2024.12.23. ishiyama