(画面左下にMathJaxのタイプセットステータスが表示されます。)
環境によっては表示が崩れているかもしれません。その際はブラウザーを変更してみてください。
ボリュームレンダリング方程式 (Volume Rendering Equation) 1
通常のレンダリング方程式に加えて、空間中の媒質における散乱も考慮に加えた方程式。さらに写実的なCGをレンダリングすることができます。
※MRIやCTのような可視化の分野でもボリュームレンダリングという言葉がありますが、それとは異なる概念です。あちらは3次元中で着目したい部分を見やすく表示する事を言いますが、こちらは物理モデルに基づいて散乱現象等をレンダリングする事を言います。
関与媒質(participating media)と光学現象
これまでのページで紹介してきた話は、真空を仮定していて光が向きや強さを変えるのは物体表面だけでした。しかし，現実をリアルに再現するにあたって真空という仮定では再現できない現象が多くあります。例えば絵画などにおいて空気遠近法で表現される、遠くの物ほど色が霞んで見える現象がありますが、これは遠くから観測者に届く光が途中で大気中で散乱されてしまうことによって起こります。また、雲や煙のような明確に表面は持たないが明らかに光に作用する物質、炎のように空間が光を発する現象なども物体表面だけのモデルでは表現が難しい例として挙げられます。
レンダリングの世界では空間中で光と相互作用を起こす物質のことを関与媒質(participating media)と呼びます。関与媒質に入射した光は、吸収されたり散乱したりしてその向きや強さが変わります。また、炎のように自ら発光している媒質を考えるなら光エネルギーの増加も考えられます。以下ではそれぞれの現象による放射輝度の変化について見ていきます。
吸収(absorption)
光が関与媒質を微小距離 $ \d s $ 突き抜けるとき、光が媒質に吸収(absorption)されることによって放射輝度の減少が起こります。その放射輝度の変化 $ \d L $ は次の式で与えられます。
\begin{equation}
\d L(\vx, \vomega) = -\sigma_a(\vx) L(\vx, \vomega) \d s \label{eq:absorption}
\end{equation}
ここで、$ \sigma_a(\vx) $ は吸収係数(absorption coefficient, absorption cross section)で、正の値を持ちます。単位は $ [\mathrm{m}^{-1}] $ です。$ \vx $ の関数となっており、媒質の密度によって変化します。吸収は放射輝度が減る現象なので右辺にマイナスの記号がついています。また、減少する放射輝度は入射する放射輝度に比例します。
この微分方程式を変数分離を使って $ L $ について解くことで、$ \vomega $ にそってある距離 $ d $ 進んだときの放射輝度の減少量がわかります。
\begin{eqnarray*}
\frac{1}{L(\vx, \vomega)} \d L(\vx, \vomega) &=& -\sigma_a(\vx) \d s \\
\int_{L(\vx(0), \vomega)}^{L(\vx(d), \vomega)} \frac{1}{L(\vx, \vomega)} \d L(\vx, \vomega) &=& -\int_0^d \sigma_a(\vx(s)) \d s \\
\Brk{\ln{L(\vx(d), \vomega)} - \ln{L(\vx(0), \vomega)}} &=& -\int_0^d \sigma_a(\vx(s)) \d s \\
L(\vx(d), \vomega) &=& L(\vx(0), \vomega) \exp{\Prtc{-\int_0^d \sigma_a(\vx(s)) \d s}}
\end{eqnarray*}
この式より、最初に放射輝度 $ L(\vx(0), \vomega) $ が入射した場合、距離を進めるに従って吸収係数に応じて指数関数的に減衰していくことがわかります。吸収係数が定数 $ \sigma_a(\vx) = \sigma_a $ だった場合には本当に単純な指数関数的減衰になります。
\begin{equation*}
L(\vx(d), \vomega) = L(\vx(0), \vomega) \exp{(-\sigma_a d)}
\end{equation*}
放射輝度が減少する散乱(out-scattering)
光が関与媒質を微小距離 $ \d s $ 突き抜けるとき、光が媒質中の微粒子などで散乱されることによって放射輝度の減少が起こります。これをout-scatteringと呼び、放射輝度の変化 $ \d L $ は次の式で与えられます。
\begin{equation}
\d L(\vx, \vomega) = -\sigma_s(\vx) L(\vx, \vomega) \d s \label{eq:out_scattering}
\end{equation}
ここで、$ \sigma_s(\vx) $ は散乱係数(scattering coefficient, scattering cross section)で、正の値を持ちます。単位は $ [\mathrm{m}^{-1}] $ です。$ \vx $ の関数となっており、媒質の密度によって変化します。外へ向かって逸れる散乱も放射輝度が減る現象なので右辺にマイナスの記号がついています。また、同様に減少する放射輝度は入射する放射輝度に比例します。
散乱係数と吸収係数の和 $ \sigma_e(\vx) = \sigma_a(\vx) + \sigma_s(\vx) $ を消散係数(消滅係数, extinction coefficient, attenuation coefficient)と呼び、以下に示すような消散係数に対する散乱係数の割合 $ \Lambda $ を散乱アルベド(scattering albedo)と呼びます。
\begin{equation*}
\Lambda = \frac{\sigma_s(\vx)}{\sigma_e(\vx)}
\end{equation*}
吸収とout-scatteringを合わせた放射輝度の変化率は次の式になります。 \begin{equation*} \d L(\vx, \vomega) = -\sigma_e(\vx) L(\vx, \vomega) \d s \end{equation*} この微分方程式を吸収の場合と同様に解くことで次の式を得ます。 \begin{equation*} L(\vx(d), \vomega) = L(\vx(0), \vomega) \exp{\Prtc{-\int_0^d \sigma_e(\vx(s)) \d s}} \end{equation*} 結果として媒質内を進む放射輝度は、媒質の吸収係数と散乱係数の和である消散係数に応じて指数関数的に減衰します。ちなみに消散係数の逆数は、フォトンが一度散乱を起こして次の散乱・吸収まで進む距離の平均に等しく、その距離を平均自由行程(mean free path)と呼びます。
放射輝度が増加する散乱(in-scattering)
外へ向かう散乱があれば、当然入ってくる散乱もありin-scatteringと呼ばれます。放射輝度の変化 $ \d L $ は次の式で与えられます。 \begin{equation} \d L(\vx, \vomega) = \sigma_s(\vx) \int_{\mathcal{S}^2} f_p(\vx, \vomega', \vomega) L(\vx, \vomega') \d\vomega' \d s \label{eq:in_scattering} \end{equation} $ \sigma_s(\vx) $ はout-scatteringと同じ散乱係数です。$ f_p(\vx, \vomega', \vomega) $ は位相関数と呼ばれ、$ \vomega' $ の方向で入射した光が $ \vomega $ 方向にどの程度散乱するかを表します。詳細は後述します。$ \vx $ において全ての方向から入射する放射輝度に関して積分することで増加量が計算できます。この場合の放射輝度の変化量は元の放射輝度には依存しません。
発光(emission)
炎のように発光(emission)する媒質の場合は、放射輝度が増加します。放射輝度の変化 $ \d L $ は次の式で与えられます。 \begin{equation} \d L(\vx, \vomega) = L_e^V(\vx, \vomega) \d s \label{eq:emission} \end{equation} ここで、$ L_e^V(\vx, \vomega) $ はvolume emittance function(体積放射率関数？)で $ \vx $ から $ \vomega $ 方向に自発光として出射される単位長あたりの放射輝度を表します。単位は $ [\mathrm{W \cdot m^{-3} \cdot sr^{-1}}] $ です。発光の場合も変化量は元の放射輝度に依存しません。
位相関数(phase function)
物体表面の反射や透過を考えるときにはBSDFによって、入射方向と出射方向の組み合わせに対する光の伝達率が記述されていました。位相関数はその関与媒質版と言えるものですが、いくつかの点において違います。最も大きな違いは位相関数は正規化されているということです。式で表すと次のようになります。 \begin{equation*} \int_{\mathcal{S}^2} f_p(\vx, \vomega', \vomega) \d\vomega = 1 \end{equation*} BSDFでは積分値によってそれ自体が反射や吸収の比を表していましたが、関与媒質における散乱と吸収の量は散乱係数と吸収係数によって定められます(注意点として、散乱されなかったものが全て吸収されるというわけではなく、一部は何の反応も起こさず媒質を突き抜けることがあります)。位相関数の単位は $ \Brka{\mathrm{sr}^{-1}} $ です。
位相関数もBSDFと同じようにその散乱の仕方によって大きく3種類に分類されます。
あらゆる方向に対して同じような強さで散乱する場合、等方散乱と呼ばれます。完全な等方散乱の場合は全立体角に対して等しい散乱を起こすため、その値は $ 1 / 4\pi $ の定数になります。入射方向に優位的に散乱する場合は前方散乱、逆に入射方向とは逆方向が優位的になる場合は後方散乱と呼ばれます。なお、位相関数における入射方向はBSDFとは違い実際の光の向きで記述されることが多いようです。このウェブサイトでも明記しない限りそのように記述します。位相関数は大抵入射方向 $ \vomega' $ と出射方向 $ \vomega $ のなす角 $ \theta $ にのみ依存するようです。
基本形のボリュームレンダリング方程式
空間中のある点 $ \vx $ に $ -\vomega $ 方向から入射する放射輝度 $ L_i $ は次のように表されます。 \begin{eqnarray*} L_i(\vx, -\vomega) = T(\vx_S, \vx) L_o(\vx_S, \vomega) + \int_{0}^{\norma{\vx_S - \vx}} T(\vx', \vx) L_a^V(\vx', \vomega) \d s \end{eqnarray*} \begin{eqnarray*} \vx' := \vx - s \vomega \hspace{10mm} \vx_S := \raycast \Prta{\vx, -\vomega} \end{eqnarray*} ここで、$ \vx_S $ はレイキャスト関数 $ \raycast $ で表されているように、$ \vx $ から $ -\vomega $ 方向にレイを飛ばして最初に見つかる物体表面の点を表しています。$ L_o(\vx_S, \vomega) $ は、$ \vx_S $ を $ \vomega $ 方向に出射した放射輝度を表していて、$ \vx $ に到達するまでに散乱・吸収を起こして残った割合が $ T(\vx_S, \vx) $ で表されています。$ L_a^V(\vx', \vomega) $ は、$ \vx $ と $ \vx_S $ の間の点 $ \vx' $ で散乱や発光によって増加した単位距離あたりの放射輝度を表していて、同様に減衰後に残った割合が $ T(\vx', \vx) $ で表されています。そして区間で積分することで $ \vx_S \vx $ 間で増加した放射輝度の総量を表しています。 直感的には「物体表面を出た光が減衰したもの」と「途中で得た光が減衰したもの」の和となっています。増加する放射輝度を発光とin-scatteringによって表すことで基本形のボリュームレンダリング方程式(Volume Rendering Equation, Volume Light Transport Equation, 図6)を得ます。 \begin{eqnarray} L_i \Prta{\vx, -\vomega} = T \, \Prta{\vx_S, \vx} L_o \Prta{\vx_S, \vomega} + \int_0^{\norma{\vx_S - \vx}} T \, \Prta{\vx', \vx} \Brk{ L_e^V \Prta{\vx', \vomega} + \sigma_s \Prta{\vx'} \int_{\mathcal{S}^2} f_p \Prta{\vx', \vomega', \vomega} L \Prta{\vx', \vomega'} \d\vomega' } \d s \label{eq:VLTE_basic_form1} \end{eqnarray} ボリューメトリックな現象をレンダリングする場合はこのボリュームレンダリング方程式を解くことになります。この方程式の導出をこのページの最後のセクションで示しています。散乱・吸収による減衰度合いである透過率(transmittance) $ T $ は消散係数 $ \sigma_e $ を用いて次のように表されます。 \begin{eqnarray*} T \, \Prta{\Vec{a}, \Vec{b}} := \exp \Prtc{-\int_0^{\norma{\Vec{a} - \Vec{b}}} \sigma_e \Prta{\Vec{b} - t\vomega_{\Vec{a}\Vec{b}}} \d t} \end{eqnarray*} $ T $ の指数内の積分を光学的深度(optical depth)や光学的厚さ(optical thickness)と呼びます。ボリュームレンダリング方程式は境界条件として以下の物体表面のレンダリング方程式を持ちます。 \begin{eqnarray*} L_o \Prta{\vx_S, \vomega} = L_e^S \Prta{\vx_S, \vomega} + \int_{\mathcal{S}^2} f_s \Prta{\vx_S, \vomega_i, \vomega} L_i \Prta{\vx_S, \vomega_i} \absb{\vomega_i \! \cdot \vec{n}} \d\vomega_i \end{eqnarray*} 物体表面のレンダリング方程式は「レンダリング方程式(Rendering Equation)」で見た通りですね。物体表面のレンダリング方程式と表記を統一するために、式 \eqref{eq:VLTE_basic_form1} の放射輝度の表記を左辺は $ L_o $、右辺は $ L_i $ を使って書き換えると次のようになります。 \begin{eqnarray*} L_o \Prta{\vx, \vomega} = T \, \Prta{\vx_S, \vx} L_o \Prta{\vx_S, \vomega} + \int_0^{\norma{\vx_S - \vx}} T \, \Prta{\vx', \vx} \Brk{ L_e^V \Prta{\vx', \vomega} + \sigma_s \Prta{\vx'} \int_{\mathcal{S}^2} f_p \Prta{\vx', -\vomega', \vomega} L_i \Prta{\vx', \vomega'} \d\vomega' } \d s \end{eqnarray*} 媒質中では界面を考える必要がないため、放射輝度に関して入射と出射を区別する必要がありません。左辺の $ L_o $ に関しては方向の表記ルールが変わっただけです。そして右辺のin-scattering項に関しては、入射光の方向記述ルールを使用したため位相関数の入射方向は反転して考える必要があります。
ボリュームレンダリング方程式は物体表面のレンダリング方程式と同様に右辺に直接評価できない未知の関数 $ L_i $ を含んでいます。「ボリュームレンダリング方程式(Volume Rendering Equation) 2」では、この方程式をただの積分問題とする定式化について説明します。
＊ボリュームレンダリング方程式の導出
吸収 $ \eqref{eq:absorption} $、散乱 $ \eqref{eq:out_scattering} $、$ \eqref{eq:in_scattering} $、発光 $ \eqref{eq:emission} $ による放射輝度の変化率全てを考慮すると以下の式が成り立ちます。
\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*}
この式は微分積分方程式形式の放射伝達方程式(Radiative Transfer Equation, RTE)と呼ばれます。この式をよく見てみると、非同次一階線形微分方程式 $ y'(x) + f(x)y(x) = g(x) $ になっていることがわかります。この場で重要ではない文字を省略して、式を簡潔なものに変形すると以下のようになります。
\begin{equation}
L'(s) + \sigma_e(s)L(s) = L_a^V(s) \label{eq:simplified_RTE}
\end{equation}
$ \vx $ は直線上の距離 $ s $ の関数だったので、$ \vx $ の関数であった係数等は結局 $ s $ の関数と言い換えられます。また、in-scatteringと発光の項をまとめて $ L_a^V $ で表しています。この微分方程式を定数変化法を用いて解いてみます。
まず右辺 = 0、つまり同次形とみなして微分方程式を解きます。これはin-scatteringと発光が無く、吸収とout-scatteringだけの状態を意味します。
\begin{eqnarray}
\frac{\d L(s)}{\d s} &=& -\sigma_e(s) L(s) \nonumber \\
\int_{L(s_0)}^{L(s)} \frac{1}{L(s')} \d L(s') &=& -\int_{s_0}^{s} \sigma_e(s') \d s' \nonumber \\
L(s) &=& L(s_0) \exp{\Prtc{-\int_{s_0}^{s} \sigma_e(s') \d s'}} \nonumber \\
&=& L_0 \exp{\Prtc{-\int_{s_0}^{s} \sigma_e(s') \d s'}} \label{eq:RTE_homogeneous_general_solution}
\end{eqnarray}
同次形の一般解が計算できました。ここで簡単のため $ L_0 = L(s_0) $ と置き換えています。また、右辺の積分では区間の文字 $ s $ と積分内の $ s $ がかぶるとあまり良くないので $ s' $ を使っています。吸収のところで説明した式と同じ形式になっていますね。次に $ L_0 = L_0(s) $ と見なし、微分方程式 $ \eqref{eq:simplified_RTE} $ の左辺に代入します。
\begin{eqnarray*}
L'(s) + \sigma_e(s)L(s) &=& \frac{\d L_0(s)}{\d s} \exp{\Prtc{-\int_{s_0}^{s} \sigma_e(s') \d s'}} \\
& & - \sigma_e(s) L_0(s) \exp{\Prtc{-\int_{s_0}^{s} \sigma_e(s') \d s'}} + \sigma_e(s) L_0(s) \exp{\Prtc{-\int_{s_0}^{s} \sigma_e(s') \d s'}} \\
&=& \frac{\d L_0(s)}{\d s} \exp{\Prtc{-\int_{s_0}^{s} \sigma_e(s') \d s'}}
\end{eqnarray*}
そして元の微分方程式 $ \eqref{eq:simplified_RTE} $ の右辺との関係より $ L_0(s) $ に関する微分方程式がつくられます。
\begin{equation*}
\frac{\d L_0(s)}{\d s} \exp{\Prtc{-\int_{s_0}^{s} \sigma_e(s') \d s'}} = L_a^V(s)
\end{equation*}
この微分方程式を以下のようにして解きます。
\begin{eqnarray*}
\frac{\d L_0(s)}{\d s} &=& \exp{\Prtc{\int_{s_0}^{s} \sigma_e(s') \d s'}} L_a^V(s) \\
L_0(s) - L_0(s_0) &=& \int_{s_0}^{s} \exp{\Prtc{\int_{s_0}^{t} \sigma_e(s') \d s'}} L_a^V(t) \d t \\
\end{eqnarray*}
ここでも右辺の積分内では $ s $ の代わりに $ t $ という文字を使っています。これを式 $ \eqref{eq:RTE_homogeneous_general_solution} $ に代入することで以下の式を得ます。
\begin{eqnarray*}
L(s) &=& \exp{\Prtc{-\int_{s_0}^{s} \sigma_e(s') \d s'}} \Brk{L_0(s_0) + \int_{s_0}^{s} \exp{\Prtc{\int_{s_0}^{t} \sigma_e(s') \d s'}} L_a^V(t) \d t} \\
&=& \exp{\Prtc{-\int_{s_0}^{s} \sigma_e(s') \d s'}} L_0(s_0) + \int_{s_0}^{s} \exp{\Prtc{-\int_{t}^{s} \sigma_e(s') \d s'}} L_a^V(t) \d t
\end{eqnarray*}
これで導出はほとんど終わりです。この式の意味するところは、位置 $ s $ における放射輝度 $ L(s) $ は、初期位置 $ s_0 $ における放射輝度 $ L_0(s_0) $ が減衰したもの(第一項目)と、$ s_0 $ から $ s $ へ辿る過程で追加された放射輝度が減衰したもの(第二項目)の和である、ということです。最後に初期位置 $ s_0 $ を物体表面 $ \vx_S $ と考え、境界条件である物体表面のレンダリング方程式より $ L_0(s_0) = L_o(\vx_S, \vomega) $ と置き換え、省略した文字などを再度記述することでボリュームレンダリング方程式の基本形を得ます。
参考文献
- [Jensen2002] Henrik Wann Jensen(著), 苗村 健(訳) - "フォトンマッピング―実写に迫るコンピュータグラフィックス", 2002, オーム社
- [Marschner2009] Steve Marschner - "Volume light transport", 2009
- [Pauly1999] Mark Pauly - "Robust Monte Carlo Methods for Photorealistic Rendering of Volumetric Effects", 1999