ボリュームレンダリング方程式 (Volume Rendering Equation) 2

このページはレイトレ Advent Calendar 2019の記事として制作しました。当初はトラッキングアルゴリズムのPDFについての記事を書こうかと思ってましたが、そういえばこのページもずっと書かずに放置していたことに気がついたのでこの機会に埋めておきます。

このページではボリュームレンダリング方程式の純粋な積分問題への変形や、それに伴う被積分関数について紹介します。

ボリュームレンダリング方程式のおさらい

ボリュームレンダリング方程式は、微分積分方程式である放射伝達方程式(RTE)を変形して得られる純粋な積分方程式で、次の形で表されます。各項の説明や導出に関しては「ボリュームレンダリング方程式1」を参照してください。 \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*} \begin{equation*} \vx_S = \raycast(\vx, -\vomega) \end{equation*} ここで $ \raycast(\vx, \vomega) $ はレイキャスト関数で、点 $ \vx $ から $ \vomega $ 方向にレイを飛ばして見つかる最も近い物体表面上の点を返します。ボリュームレンダリング方程式の境界条件として以下に示す物体表面のレンダリング方程式が設定され、第一項目の $ L_o(\vx_S, \vomega) $ に対応します。 \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', \vomega} L_i \Prta{\vx_S, \vomega'} \absb{\vomega' \! \cdot \vec{n}} \d\vomega' \end{eqnarray*} これをボリュームレンダリング方程式に代入してみます。また媒質中の場合、ある点に入射する光と出射する光は方向が同じであれば区別する必要が無いので、左辺の $ L_o(\vx, \vomega) $ を $ -\vomega $ 方向「から」入射する光 $ L_i(\vx, -\vomega) $ として表します。以下の式を得ます。 \begin{eqnarray} L_i \Prta{\vx, -\vomega} &=& T \, \Prta{\vx_S, \vx} \Brk{ L_e^S \Prta{\vx_S, \vomega} + \int_{\mathcal{S}^2} f_s \Prta{\vx_S, \vomega', \vomega} L_i \Prta{\vx_S, \vomega'} \absb{\vomega' \! \cdot \vec{n}} \d\vomega' } + \nonumber \\ && \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 \label{eq:substituted_VLTE} \end{eqnarray} 右辺の一項目と二項目がかなり似た表現になっていることに気がつくと思います。物体表面か媒質中かの違いはありますが、いずれも発光に関する項と入射光の散乱を積分する項を持っています。そしてそれら二つの項に減衰がかかっていることも共通しています。この式を見れば、空間中のある点への入射光は、別の点における入射光によって再帰的に定義されていることがわかります。そのためボリュームレンダリング方程式は「入射光の散乱」と「減衰」という2つのステップの繰り返しと捉えることができます。

記号の定義

オペレーターの定義とノイマン級数展開

図1. 散乱オペレーター $ \mathbf{K} $

物体表面上・媒質中の入射光の場に作用して散乱による出射光の場に変換する。

図2. 輸送オペレーター $ \mathbf{G} $

$ \vomega $ 方向から来る(減衰した)出射光(発光・散乱を区別しない)の場を積分して入射光の場に変換する。

ボリュームレンダリング方程式も物体表面のレンダリング方程式と同様に、オペレーターを使った表記が可能です。前節で述べた2つのステップに対応するオペレーターを定義します。1つ目の散乱オペレーター $ \mathbf{K} $ (図1)を以下のように定義します。 \begin{equation} \newcommand\opK{\mathbf{K}} \newcommand\opG{\mathbf{G}} \newcommand\entspace{\Omega} \newcommand\volspace{\mathcal{V}} < \opK g > (\vx, \vomega) := \left\{ \begin{array}{ll} \displaystyle \int_{\mathcal{S}^2} f_s(\vx, \vomega', \vomega) g(\vx, \vomega') \abs{\vec{n} \cdot \vomega'} \d\vomega' & \vx \in \mathcal{M} \\ \displaystyle \sigma_s(\vx) \int_{\mathcal{S}^2} f_p(\vx, -\vomega', \vomega) g(\vx, \vomega') \d\vomega' & \vx \in \volspace \end{array} \right. \label{eq:scattering_operator} \end{equation} このオペレーターは入射光 $ g $ に作用して、$ \vx \in \mathcal{M} $、つまり点 $ \vx $ が物体表面上の場合は表面で散乱(反射・透過)して出射する放射輝度に、媒質中の場合は媒質から散乱して出射する単位長あたりの放射輝度に変換します。
次に2つ目の減衰に対応する輸送オペレーター $ \mathbf{G} $ (図2)を以下のように定義します。 \begin{equation} < \opG g > (\vx, \vomega) := \int_0^{\norm{\vx_S - \vx}} T(\vx', \vx) g(\vx', -\vomega) \d s + T(\vx_S, \vx) g(\vx_S, -\vomega) \label{eq:transport_operator} \end{equation} このオペレーターは媒質中・表面上両方の出射光 $ g $ に作用して、点 $ \vx $ に入射する放射輝度に変換します。入射光の方向 $ \vomega $ を基準としているため、出射光 $ g $ の方向は逆転して考える必要があります。2つのオペレーターの定義中の $ g $ に関して、方向の記述ルールに疑問を感じるかもしれません。これらのオペレーターは関数自体を引数として受け取るため $ g $ はただの変数であり、二者の定義中の $ g $ は別物です。

これらのオペレーターを使って式 \eqref{eq:substituted_VLTE} を以下のように簡潔に表すことができます。 \begin{equation} L_i = \opG (L_e + \opK L_i) \label{eq:Li_operator_notation} \end{equation} ここで $ L_e $ は媒質中の単位長あたりの発光輝度、もしくは物体表面上の発光輝度を表しています。 \begin{equation*} L_e(\vx, \vomega) := \left\{ \begin{array}{ll} \displaystyle L_e^S(\vx, \vomega) & \vx \in \mathcal{M} \\ \displaystyle L_e^V(\vx, \vomega) & \vx \in \volspace \end{array} \right. \end{equation*} \eqref{eq:Li_operator_notation} の左辺を右辺に再起的に代入することで次のノイマン級数を得ます。 \begin{equation*} L_i = \opG L_e + \opG \opK \opG L_e + (\opG \opK)^2 \opG L_e + \cdots = \sum_{n=0}^\infty (\opG\opK)^n \opG L_e \end{equation*} この式は、ある点に届く放射輝度 $ L_i $ は光源から出た光 $ L_e $ が0, 1, 2, ...回散乱・減衰したものの和で表されることを意味しています。一回も散乱せずに届く場合でも減衰は必ず発生するため $ \opG $ がひとつ追加でかかっています。

変数変換

物体表面のレンダリング方程式を扱う際、次に示す立体角測度と面積測度の変数変換 \begin{equation} \int_{\mathcal{S}^2} f(\raycast(\vx, \vomega)) \d\vomega = \int_{\mathcal{M}} f(\vx') \frac{\abs{\vec{n}' \cdot \vomega_{\vx \vx'}}}{\norm{\vx' - \vx}^2} V(\vx' \LRAR \vx) \d A(\vx') \label{eq:change_of_variables_surface} \end{equation} を用いました。この変換によって、任意の物体表面上の点 $ \vx' $ を「点 $ \vx $ から $ \vomega $ 方向に見える点」と表現する代わりに「物体表面上の点 $ \vx' $」と直接的に表現した式を構築することが可能になります。

ボリュームレンダリング方程式における、媒質中の位置 $ \vx' $ に関しても「点 $ \vx $ から $ \vomega $ 方向に距離 $ r $ にある点」という表現と「媒質中の点 $ \vx' $」という直接的な表現を用いることが可能で、それらの間には次に示す変数変換が成り立ちます。 \begin{eqnarray} \int_{\mathcal{S}^2} \int_0^{\norm{\raycast(\vx, \vomega) - \vx}} f(\vx + r\vomega) \d r \d\vomega = \int_\volspace \frac{1}{\norm{\vx' - \vx}^2} V(\vx' \LRAR \vx) f(\vx') \d V(\vx') \label{eq:change_of_variables_volume} \end{eqnarray}

経路積分による定式化

光センサーの測定値は物体表面上の光輸送を考えたときと同様に、次に示すmeasurement equationによって表されます。 \begin{equation*} I_j = \int_{\mathcal{M}} \int_{\mathcal{S}^2} W_e^j (\vx, \vomega) L_i (\vx, \vomega) \abs{\vec{n} \cdot \vomega} \d\vomega \d A(\vx) \end{equation*} 「レンダリング方程式2」に示したmeasurement equationに対して、ボリュームレンダリングの場合は減衰を考える必要があるため、センサーから離れた場所の出射光ではなく、センサーへの入射光を用いて表現する必要があります。また、センサーに入射する光の散乱点を統一的に扱うため、内側の積分範囲は物体表面上ではなく一旦入射光の方向として記述しています。 式 \eqref{eq:Li_operator_notation} をこのmeasurement equationに代入します。 \begin{equation*} I_j = \int_{\mathcal{M}} \int_{\mathcal{S}^2} W_e^j (\vx, \vomega) \Prt{\opG (L_e + \opK L_i)} (\vx, \vomega) \abs{\vec{n} \cdot \vomega} \d\vomega \d A(\vx) \end{equation*} 輸送オペレーター $ \opG $ を展開すると次の式を得ます。 \begin{eqnarray*} I_j &=& \int_{\mathcal{M}} \int_{\mathcal{S}^2} W_e^j (\vx, \vomega) \Bigl[ \int_0^{\norm{\raycast(\vx, \vomega) - \vx}} T(\vx', \vx) \Brk{L_e + \opK L_i}(\vx', -\vomega) \d s \\ && \hspace{30mm} + \, T(\raycast(\vx, \vomega), \vx) \Brk{L_e + \opK L_i}(\raycast(\vx, \vomega), -\vomega) \Bigr] \abs{\vec{n} \cdot \vomega} \d\vomega \d A(\vx) \end{eqnarray*} 式 \eqref{eq:change_of_variables_surface}, \eqref{eq:change_of_variables_volume} をそれぞれ適用すると、二つの積分を媒質中と物体表面上の位置に関する積分に変換することができます。 \begin{eqnarray*} I_j &=& \int_{\mathcal{M}} \left[ \int_\volspace W_e^j (\vx' \RAR \vx) T(\vx', \vx) \frac{1}{\norm{\vx' - \vx}^2} V(\vx' \LRAR \vx) \Brk{L_e + \opK L_i}(\vx' \RAR \vx) \d V(\vx') \right. \\ && \left. \hspace{6mm} + \int_{\mathcal{M}} W_e^j (\vx' \RAR \vx) T(\vx', \vx) \frac{\abs{\vec{n}' \cdot \vomega_{\vx \vx'}}}{\norm{\vx' - \vx}^2} V(\vx' \LRAR \vx) \Brk{L_e + \opK L_i}(\vx' \RAR \vx) \d A(\vx') \right] \abs{\vec{n} \cdot \vomega} \d A(\vx) \end{eqnarray*} 式の簡潔な記述のために、媒質中と物体表面上両方に一般化し可視関数 $ V $ や減衰も含めた新たな幾何項を定義します。 \begin{equation*} G(\Vec{a} \LRAR \Vec{b}) := T(\Vec{a}, \Vec{b}) V(\Vec{a} \LRAR \Vec{b}) \frac{D_\Vec{a}(\vomega_{\Vec{a}\Vec{b}}) D_\Vec{b}(\vomega_{\Vec{b}\Vec{a}})}{\norm{\Vec{b} - \Vec{a}}^2} \end{equation*} ここで、$ D $ はcos項を一般化したもので、物体表面上でのみcos項としてはたらき、媒質中では $ 1 $ になります。 \begin{equation*} D_\Vec{a}(\vomega) := \left\{ \begin{array}{ll} \displaystyle \abs{\vec{n}_\Vec{a} \cdot \vomega} & \Vec{a} \in \mathcal{M} \\\ \displaystyle 1 & \Vec{a} \in \volspace \end{array} \right. \end{equation*} 一般化・拡張された幾何項を用いると次の式を得ます。 \begin{eqnarray*} I_j &=& \int_{\mathcal{M}} \left[ \int_\volspace W_e^j (\vx' \RAR \vx) G(\vx' \LRAR \vx) \Brk{L_e + \opK L_i}(\vx' \RAR \vx) \d V(\vx') \right. \\ && \left. \hspace{6mm} + \int_{\mathcal{M}} W_e^j (\vx' \RAR \vx) G(\vx' \LRAR \vx) \Brk{L_e + \opK L_i}(\vx' \RAR \vx) \d A(\vx') \right] \d A(\vx) \end{eqnarray*} $ \opK L_i = \sum_{n = 0}^\infty \opK (\opG \opK)^n \opG L_e = \sum_{n = 1}^\infty (\opK \opG)^n L_e $ であることを利用すると次のように経路長ごとの総和に変形できます。 \begin{eqnarray} I_j &=& \sum_{n = 0}^\infty \int_{\mathcal{M}} \left[ \int_\volspace W_e^j (\vx' \RAR \vx) G(\vx' \LRAR \vx) \Brk{(\opK \opG)^n L_e}(\vx' \RAR \vx) \d V(\vx') \right. \nonumber \\ && \left. \hspace{14mm} + \int_{\mathcal{M}} W_e^j (\vx' \RAR \vx) G(\vx' \LRAR \vx) \Brk{(\opK \opG)^n L_e}(\vx' \RAR \vx) \d A(\vx') \right] \d A(\vx) \label{eq:measurement_equation_with_KG} \end{eqnarray}

図3. 結合されたオペレーター $ \opK \opG $ ($ \vx \in \mathcal{M} $ の場合)

出射光の場を(減衰を考慮して)収集・積分して散乱光の場に変換する。

ここで、散乱オペレーターと輸送オペレーターの連結 $ \opK \opG $ をまとめてみます(図3)。散乱点が物体表面上 $ \vx \in \mathcal{M} $ の場合は次の式で表されます。 \begin{eqnarray*} < \opK \opG g > (\vx, \vomega) &=& \int_{\mathcal{S}^2} f_s(\vx, \vomega', \vomega) \left[ \int_0^{\norm{\raycast(\vx, \vomega') - \vx}} T(\vx', \vx) g(\vx', -\vomega') \d s \right. \\ && \hspace{50mm} \left. + \, T(\raycast(\vx, \vomega'), \vx) g(\raycast(\vx, \vomega'), -\vomega') \vphantom{\int} \right] \abs{\vec{n} \cdot \vomega'} \d\vomega' \end{eqnarray*} 変数変換を行って媒質中と物体表面上の位置に関する積分に変換します。 \begin{eqnarray*} < \opK \opG g > (\vx \RAR \vx'') &=& \int_\volspace f_s(\vx' \RAR \vx \RAR \vx'') G(\vx' \LRAR \vx) g(\vx' \RAR \vx) \d V(\vx') \\ && \hspace{50mm} + \int_\mathcal{M} f_p(\vx' \RAR \vx \RAR \vx'') G(\vx' \LRAR \vx) g(\vx' \RAR \vx) \d A(\vx') \end{eqnarray*} $ \vx'' $ は散乱後に光が進む方向 $ \vomega $ 上にある任意の点ですが、位置だけの式にする上で必要だっただけです。同様に、散乱点が媒質中 $ \vx \in \volspace $ の場合も変数変換を経て次の式を得ます。 \begin{eqnarray*} < \opK \opG g > (\vx \RAR \vx'') &=& \sigma_s(\vx) \left[ \int_\volspace f_s(\vx' \RAR \vx \RAR \vx'') G(\vx' \LRAR \vx) g(\vx' \RAR \vx) \d V(\vx') \right. \\ && \left. \hspace{50mm} + \int_\mathcal{M} f_p(\vx' \RAR \vx \RAR \vx'') G(\vx' \LRAR \vx) g(\vx' \RAR \vx) \d A(\vx') \right] \end{eqnarray*} 物体表面上の散乱と媒質中の散乱を次に示す一般化された散乱の分布関数 $ \bar{f} $ にまとめます。 \begin{equation*} \bar{f}(\vx' \RAR \vx \RAR \vx'') := \left\{ \begin{array}{ll} \displaystyle f_s(\vx' \RAR \vx \RAR \vx'') & \vx \in \mathcal{M} \\ \displaystyle \sigma_s(\vx) \, f_p(\vx' \RAR \vx \RAR \vx'') & \vx \in \volspace \end{array} \right. \end{equation*} これを用いると上記の結合オペレーターは次のように、全ての領域 $ \entspace $ においてひとつの定義で対応できます。 \begin{eqnarray*} < \opK \opG g > (\vx \RAR \vx'') &=& \int_\volspace \bar{f}(\vx' \RAR \vx \RAR \vx'') G(\vx' \LRAR \vx) g(\vx' \RAR \vx) \d V(\vx') \\ && \hspace{50mm} + \int_\mathcal{M} \bar{f}(\vx' \RAR \vx \RAR \vx'') G(\vx' \LRAR \vx) g(\vx' \RAR \vx) \d A(\vx') \end{eqnarray*} これを式 \eqref{eq:measurement_equation_with_KG} に代入することで経路長に対応した次元の多重積分(の総和)を得ることができます。光源の点と散乱点が物体表面上の位置か媒質中の位置か、に対応して散乱回数 $ n $、つまり経路長 $ k = n + 1 $ の場合は $ 2^k $ の項数の多重積分になります。

ここで、経路長 $ k $ の各散乱が物体表面上だったか媒質中だったかに応じて個別の経路空間 $ \mathcal{P}_k^{\Vec{c}} $ を次のように定義します。 \begin{equation*} \newcommand\bigtimes{\mathop{\vcenter{\huge\times}}} \mathcal{P}_k^{\Vec{c}} := \bigtimes_{i = 0}^k \left\{ \begin{array}{ll} \displaystyle \mathcal{M} & \mbox{if} \,\, \Vec{c}_i = 0 \\ \displaystyle \mathcal{V} & \mbox{if} \,\, \Vec{c}_i = 1 \end{array} \right. \end{equation*} ここで、$ \Vec{c} \in \Brc{0, 1}^{k + 1} $ は経路長 $ k $ における頂点 $ \vx_i $ ($ i \in [0, k] $) が物体表面上か媒質中かを $ \Vec{c}_i $ によって表すベクトルです。$ \bigtimes $ は直積集合、デカルト積を表しており、各集合から取り出した元の組み合わせを元として持つ集合を作ります。経路空間 $ \mathcal{P}_k^{\Vec{c}} $ の表記を用いるとボリュームレンダリング方程式の経路積分定式化は次の形で表すことができるようになります。 \begin{equation*} I_j = \sum_{k = 1}^\infty \sum_{\Vec{c} \in \Brc{0, 1}^{k + 1}} \int_{\mathcal{P}_k^{\Vec{c}}} f_j(\bar{\vx}_k) \d\mu_k^{\Vec{c}}(\bar{\vx}_k) \end{equation*} 経路 $ \bar{\vx}_k $ に対応する微小測度 $ \d\mu_k^{\Vec{c}}(\bar{\vx}_k) $ は次の形で表されます。 \begin{equation*} \d\mu_k^{\Vec{c}}(\bar{\vx}_k) := \prod_{i = 0}^k \left\{ \begin{array}{ll} \displaystyle \d A(\vx_i) & \mbox{if} \,\, \Vec{c}_i = 0 \\ \displaystyle \d V(\vx_i) & \mbox{if} \,\, \Vec{c}_i = 1 \end{array} \right. \end{equation*} そして被積分関数であるmeasurement contribution functionは次のようになります。 \begin{eqnarray*} f_j(\bar{\vx}_k) = W_e^j(\vx_1 \RAR \vx_0) G(\vx_1 \LRAR \vx_0) \prod_{i = 1}^{k - 1} \Brc{\bar{f}(\vx_{i + 1} \RAR \vx_i \RAR \vx_{i - 1}) G(\vx_{i + 1} \LRAR \vx_i)} L_e(\vx_k \RAR \vx_{k - 1}) \end{eqnarray*} それぞれの項が一般化されていたり、拡張されていたりしますが、物体表面間の光輸送におけるmeasurement contribution functionと同様の形になったことがわかります。さいごに、あらゆる経路長、あらゆる散乱の組み合わせの経路空間 $ \mathcal{P} $ を次のように定義します。 \begin{eqnarray*} \mathcal{P} := \bigcup_{k = 1}^\infty \bigcup_{\Vec{c} \in \Brc{0, 1}^{k + 1}} \mathcal{P}_k^{\Vec{c}} \end{eqnarray*} これによってボリュームレンダリング方程式も(無事？)純粋な積分問題に帰着させることができます。 \begin{eqnarray*} I_j = \int_{\mathcal{P}} f_j(\bar{\vx}) \d\mu(\bar{\vx}) \end{eqnarray*}

参考文献

• [Jakob2013] Wenzel Jakob - "LIGHT TRANSPORT ON PATH-SPACE MANIFOLDS", 2013
• [Pauly1999] Mark Pauly - "Robust Monte Carlo Methods for Photorealistic Rendering of Volumetric Effects", 1999