トラッキングアルゴリズムの積分定式化

このページではまず、null-collisionを考慮したボリュームレンダリング方程式を紹介して、デルタトラッキングにモンテカルロ積分としての簡単な解釈を与えます。次に、デルタトラッキングを始めとする自由行程サンプリングアルゴリズム(が扱う問題)の積分形式の定式化について紹介します。単純な積分問題として定式化することで、特定の領域に捉われないクリアな視点が得られ、新たなアルゴリズムを開発する上での強力な礎となります。

null-collisionを考慮したボリュームレンダリング方程式

デルタトラッキングでは偽の粒子との衝突(= null-collision)が起こった場合には再度自由行程を延長するサンプリングを行うため、ひとつの自由行程のサンプリングのために複数かつ不定個の乱数を使用しました。純粋な物理モデルである放射伝達方程式(RTE)もしくはボリュームレンダリング方程式には当然偽の粒子に関する記述は無く、反射方向のサンプリング等、積分(や確率的な和)と一対一対応した乱数の使用と比べると少し扱いにくさがあります。
そこで、RTEにnull-collisionを説明する項を追加してみます。まずは通常のRTEを以下に示します。 \begin{equation*} \frac{\d L_i(\vx, \vomega)}{\d s} = -\sigma_e(\vx) L_i(\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*} この式は、位置 $ \vx $ における $ \vomega $ 方向への放射輝度 $ L_i $ が微小距離 $ \d s $ 進んだ時に変化する量を表しています。右辺第一項目は吸収とout-scatteringによる放射輝度の減少、第二項目、三項目がそれぞれin-scatteringと発光による放射輝度の増加を表しています(詳しくは「ボリュームレンダリング方程式1」参照)。null-collisionを表す項を追加すると次のようになります。 \begin{equation} \begin{split} \frac{\d L_i(\vx, \vomega)}{\d s} &= -\sigma_e(\vx) L_i(\vx, \vomega) - \sigma_n(\vx) L_i(\vx, \vomega) + \hphantom{A} \\ & \quad \sigma_s(\vx) \int_{\mathcal{S}^2} f_p(\vx, \vomega', \vomega) L(\vx, \vomega') \d\vomega' + \sigma_n(\vx)\int_{\mathcal{S}^2} \delta(\vomega' - \vomega) L(\vx, \vomega') \d\vomega' + \\ & \quad L_e^V(\vx, \vomega) \end{split} \label{eq:RTE_null_collision} \end{equation} ここで $ \sigma_n(\vx) $ は(デルタトラッキングでは一様化のために空間を埋めた)偽の粒子の係数を表しており、増えた項は $ \sigma_n(\vx) $ を含む二項目と四項目です。偽の粒子と光が衝突することによる放射輝度の減少(out-scattering)が第二項目で表されています。一方で、偽の粒子によるin-scatteringで増加する放射輝度が第四項目で表されています。これら二つの項はお互いにキャンセルしあうため、次の関係を持ちます。 \begin{equation*} -\sigma_n(\vx) L_i(\vx, \vomega) + \sigma_n(\vx)\int_{\mathcal{S}^2} \delta(\vomega' - \vomega) L(\vx, \vomega') \d\vomega' = 0 \end{equation*} 偽の粒子の位相関数は式中のデルタ関数で表されており、out-scattering = in-scatteringとなるため、これらの項をRTEに追加しても何も影響を与えません。

短く表記するために $ L_s(\vx, \vomega) = \int_{\mathcal{S}^2} f_p(\vx, \vomega', \vomega) L(\vx, \vomega') \d\vomega' $ を導入し、 式 $ \eqref{eq:RTE_null_collision} $ の両辺を積分、そしてデルタ関数を含む積分を整理するとnull-collisionを考慮したボリュームレンダリング方程式が得られます。 \begin{equation} 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(\vx') L_s(\vx', \vomega) + \sigma_n(\vx') L(\vx', \vomega) } \d s \label{eq:VLTE_null_collision} \end{equation} ここで $ \vx' = \vx - s\vomega $ です。偽の粒子によるin-scatteringで増加する放射輝度を表す項が増えています。また、この場合の透過率 $ T $ は消散係数 $ \sigma_e(\vx) $と偽の粒子の係数 $ \sigma_n(\vx) $ によって次のように表されます。 \begin{eqnarray*} T \, \Prta{\Vec{a}, \Vec{b}} := \exp \Prtc{ -\int_0^{\norma{\Vec{a} - \Vec{b}}} \Prtb{\sigma_e \Prta{\Vec{b} - t\vomega_{\Vec{a}\Vec{b}}} + \sigma_n \Prta{\Vec{b} - t\vomega_{\Vec{a}\Vec{b}}}} \d t } \end{eqnarray*} 式 $ \eqref{eq:VLTE_null_collision} $ の下では、デルタトラッキングにおける自由行程延長のサンプリングは、偽の粒子によるin-scattering項を確率的に選択して評価する処理に相当します。また、この定式化の下では偽の粒子の係数 $ \sigma_n(\vx) $ が正の数である必要はなく、負の数だとしても成り立ちます。これによって自由行程サンプリングアルゴリズムの自由度が向上します。ただし、もはや物理的に妥当な解釈を与えることは難しくなります。

トラッキングアルゴリズムの積分定式化

図1. パスのウェイト

経路 $ \bar{x} $ はnull-collisionやscatteringを経て最終的に関与媒質の発光へと到達している。この経路のウェイトは各散乱イベントにおけるウェイトを積算することで計算される。

自由行程をサンプリングする際のPDFを次のように定義します。 \begin{equation} p_d(s) = \bar{\sigma}(\vx') \exp\Prt{-\int_0^s \bar{\sigma}(\vx'') \d s'} \label{eq:PDF_free_path_sampling} \end{equation} ここで、$ \vx'' = \vx - s'\vomega $ です。$ \bar{\sigma}(\vx) = \sigma_e(\vx) + \sigma_n(\vx) $ は、自由行程を解析的にサンプリングできるように偽の粒子で空間を埋めることで作られます。そのため自由行程サンプリング係数(free-path-sampling coefficient)と[Kutz2017]らは呼んでいます。デルタトラッキングでは $ \bar{\sigma}(\vx) $ が一様となるようにしていますが、解析的にサンプリングできれば一様である必要はありません。
上記PDF $ \eqref{eq:PDF_free_path_sampling} $ はトランスミッタンスを用いて次のように変形できます。 \begin{eqnarray*} p_d(s) &=& \bar{\sigma}(\vx') \exp\Prt{-\int_0^s \bar{\sigma}(\vx'') \d s'} \\ &=& \bar{\sigma}(\vx') T(\vx', \vx) \\ T(\vx', \vx) &=& \frac{p_d(s)}{\bar{\sigma}(\vx')} \end{eqnarray*} 式 $ \eqref{eq:VLTE_null_collision} $ にこれを代入することで次を得ます。 \begin{equation*} L_i \Prta{\vx, \vomega} = \frac{p_d(\norma{\vx_S - \vx})}{\bar{\sigma}(\vx_S)} L_o \Prta{\vx_S, \vomega} + \int_0^{\norma{\vx_S - \vx}} p_d(s) \Brk{ \frac{1}{\bar{\sigma}(\vx')} L_e^V \Prta{\vx', \vomega} + \frac{\sigma_s(\vx')}{\bar{\sigma}(\vx')} L_s(\vx', \vomega) + \frac{\sigma_n(\vx')}{\bar{\sigma}(\vx')} L(\vx', \vomega) } \d s \end{equation*} デルタトラッキングではサーフェスに到達した場合など関与媒質の領域を抜けると処理を終えるので、今ここで興味があるのは右辺第二項目です。これを次のように $ L_{i, v}(\vx, \vomega) $ で表しましょう。 \begin{equation*} L_{i, v} \Prta{\vx, \vomega} = \int_0^{\norma{\vx_S - \vx}} p_d(s) \Brk{ \frac{1}{\bar{\sigma}(\vx')} L_e^V \Prta{\vx', \vomega} + \frac{\sigma_s(\vx')}{\bar{\sigma}(\vx')} L_s(\vx', \vomega) + \frac{\sigma_n(\vx')}{\bar{\sigma}(\vx')} L(\vx', \vomega) } \d s \end{equation*} 次に、関与媒質中の発光、散乱、偽の散乱を確率的に評価するとして確率 $ P_e(\vx') $, $ P_s(\vx') $, $ P_n(\vx') $ を導入します。 \begin{equation} \begin{split} L_{i, v} \Prta{\vx, \vomega} &= \int_0^{\norma{\vx_S - \vx}} p_d(s) \Bigl[ P_e(\vx') \frac{1}{\bar{\sigma}(\vx') P_e(\vx')} L_e^V \Prta{\vx', \vomega} + P_s(\vx') \frac{\sigma_s(\vx')}{\bar{\sigma}(\vx') P_s(\vx')} L_s(\vx', \vomega) + \hphantom{A} \\ & \quad \hphantom{\int_0^{\norma{\vx_S - \vx}} p_d(s) \Bigl[} P_n(\vx') \frac{\sigma_n(\vx')}{\bar{\sigma}(\vx') P_n(\vx')} L(\vx', \vomega) \Bigr] \d s \end{split} \label{eq:probabilistic_evaluation} \end{equation} そして確率的な評価を積分として表すためにヘヴィサイド関数を使用します。ヘヴィサイド関数は条件が真のとき1、偽のときは0を返します。これによって確率変数 $ X $ の確率的な評価は次のように表されます。 \begin{equation*} PX = \int_0^1 \mathcal{H}[y < P] X \d y \end{equation*} したがって、式 $ \eqref{eq:probabilistic_evaluation} $ は次のように表されます。 \begin{equation} \newcommand\Heaviside[1]{\mathcal{H}\Brk{#1}} \begin{split} L_{i, v} \Prta{\vx, \vomega} &= \int_0^{\norma{\vx_S - \vx}} p_d(s) \times \hphantom{A} \\ & \quad \hspace{10mm} \Bigl[ \int_0^1 \Heaviside{u_e < P_e(\vx')} \frac{1}{\bar{\sigma}(\vx') P_e(\vx')} L_e^V \Prta{\vx', \vomega} \d u_e + \hphantom{A} \\ & \quad \hspace{10mm} \hphantom{\Bigl[} \int_0^1 \Heaviside{u_s < P_s(\vx')} \frac{\sigma_s(\vx')}{\bar{\sigma}(\vx') P_s(\vx')} L_s(\vx', \vomega) \d u_s + \hphantom{A} \\ & \quad \hspace{10mm} \hphantom{\Bigl[} \int_0^1 \Heaviside{u_n < P_n(\vx')} \frac{\sigma_n(\vx')}{\bar{\sigma}(\vx') P_n(\vx')} L(\vx', \vomega) \d u_n \Bigr] \d s \end{split} \label{eq:integral_formulation_of_tracking_algorithm} \end{equation} これがデルタトラッキングを含む各種トラッキングアルゴリズムの積分形式の定式化となります。各積分を単純にモンテカルロ積分で置き換えることで様々なアルゴリズムに直接変換されます。その場合の経路のウェイトを図1に示しています。
重要な点として、$ \sigma_n(\vx) $ や確率 $ P_e(\vx') $, $ P_s(\vx') $, $ P_n(\vx') $ の選択によらず、それらのアルゴリズムからは正しい値を推定できることが挙げられます。例えば、$ \sigma_n(\vx) $ を負の値として考えることさえ可能です。その場合 $ \sigma_n(\vx) $ に物理的な解釈を与えることはもはや困難になりますが、依然として正しい答えを求めることは可能です。物理的解釈を考えることは直感的な理解のためには重要ですが、一方で視野が狭くなりがちです。積分定式化によって物理の視点から切り離して純粋な数学的な問題に落とし込むことで、様々な視点が開かれます。

結果として、関与媒質の領域を抜けた場合も考慮に含めたモンテカルロ推定量は次のようになります。(推定量がこのような形をとる理由については「関与媒質中の自由行程サンプリングと透過率の推定」を参照してください。) \begin{equation*} \newcommand\Heaviside[1]{\mathcal{H}\Brk{#1}} L_i \Prta{\vx, \vomega} \approx \begin{cases} \displaystyle \int_0^1 \Heaviside{u_e < P_e(\vx')} \frac{1}{\bar{\sigma}(\vx') P_e(\vx')} L_e^V \Prta{\vx', \vomega} \d u_e + \hphantom{A} \\ \displaystyle \hspace{10mm} \int_0^1 \Heaviside{u_s < P_s(\vx')} \frac{\sigma_s(\vx')}{\bar{\sigma}(\vx') P_s(\vx')} L_s(\vx', \vomega) \d u_s + \hphantom{A} & s < \norma{\vx_S - \vx} \\ \displaystyle \hspace{20mm} \int_0^1 \Heaviside{u_n < P_n(\vx')} \frac{\sigma_n(\vx')}{\bar{\sigma}(\vx') P_n(\vx')} L(\vx', \vomega) \d u_n \\ \displaystyle \vphantom{\frac{T}{P_S}} L_o \Prta{\vx_S, \vomega} & s \ge \norma{\vx_S - \vx} \end{cases} \end{equation*} 自由行程のサンプリングが物体表面までの距離を超えた場合( $ s \ge \norma{\vx_S - \vx} $ )は、物体表面からの放射輝度が評価されますが、この評価を行う確率は $ \sigma_n(\vx) $ に応じて変わります。例えばデルタトラッキングの場合、この評価を行う確率が低下することになります($ \sigma_n(\vx) $ は常に正の値なため)。そのため物体表面からの輝度の期待値が偽の粒子を含まない媒質の場合に比べて小さくなりそう、と思うかもしれません。しかし、偽の粒子との衝突が評価された場合にはこのモンテカルロ推定量を同じ方向 $ \vomega $ で再帰的に評価することになるので、同じ方向の物体表面からの放射輝度を評価する機会が再度生まれることで期待値としては正しくなります。

デルタトラッキングの積分定式化

ウェイト付きデルタトラッキング

参考文献

1. [Kutz2017] Peter Kutz, Ralf Habel, Yining Karl Li, Jan Novák - "Spectral and Decomposition Tracking for Rendering Heterogeneous Volumes", 2017