(画面左下にMathJaxのタイプセットステータスが表示されます。)
環境によっては表示が崩れているかもしれません。その際はブラウザーを変更してみてください。
レンダリング方程式 (Rendering Equation) 2
このページではレンダリング方程式の基本形とは異なるバージョンについて解説します。
3点形式のレンダリング方程式
「レンダリング方程式」の式 $ (1) $, $ (2) $ では、反射光を表す項が、微小立体角に関する積分として表されていました。これをシーン中の物体表面の3点 $ \vx'', \vx, \vx' $ を考え、微小面積に関する積分として書き直すことができます。
シーン中の微小面積と微小立体角の関係は次の式で表されます。 \begin{eqnarray} \d\vomega_i = \frac{\absb{\vomega_{\vx'' \vx} \cdot \vec{n}''}}{\norm{\vx'' - \vx}^2} V \Prta{\vx'' \LRAR \vx} \d A \Prta{\vx''} \label{eq:solidangle_to_area} \end{eqnarray} ここで $ V \Prta{\vx'' \LRAR \vx} $ は $ \vx'' $ と $ \vx $ の間に遮蔽物がない場合は1、ある場合には0となる可視関数を表しています。ある微小面積 $ \d A \Prta{\vx''} $ が点 $ \vx $ 上の単位半球上に占める面積(＝微小立体角)は距離の2乗に反比例します。また上の図からもわかるように、半球に対して面が傾いている場合には微小立体角は小さくなります。
式 $ \eqref{eq:solidangle_to_area} $ の関係式を微小立体角に関するレンダリング方程式に代入する事で、3点形式のレンダリング方程式を次のように書くことができます。 \begin{eqnarray} L_o \Prta{\vx \RAR \vx'} &=& L_e \Prta{\vx \RAR \vx'} + \int_{\mathcal{M}} f_s \Prta{\vx'' \RAR \vx \RAR \vx'} L_o \Prta{\vx'' \RAR \vx} \frac{\absb{\vomega_{\vx'' \vx} \cdot \vec{n}} \absb{\vomega_{\vx'' \vx} \cdot \vec{n}''}}{\norm{\vx'' - \vx}^2} V \Prta{\vx'' \LRAR \vx} \d A \Prta{\vx''} \nonumber \\ &=& L_e \Prta{\vx \RAR \vx'} + \int_{\mathcal{M}} f_s \Prta{\vx'' \RAR \vx \RAR \vx'} L_o \Prta{\vx'' \RAR \vx} G \Prta{\vx'' \! \LRAR \vx} \d A \Prta{\vx''} \label{eq:3point_LTE} \end{eqnarray} ここで $ \mathcal{M} $ (ManifoldのM？)はシーン中の物体表面全ての集合を表しています。また、$ G \Prta{\vx'' \! \LRAR \vx} $ は幾何項と呼ばれ、次の形をとります。 \begin{eqnarray*} G \Prta{\vx'' \! \LRAR \vx} = \frac{\absb{\vomega_{\vx'' \vx} \cdot \vec{n}} \absb{\vomega_{\vx'' \vx} \cdot \vec{n}''}}{\norm{\vx'' - \vx}^2} V \Prta{\vx'' \LRAR \vx} \end{eqnarray*}
ノイマン級数展開
モンテカルロ法では、積分したい関数を区間内でランダムにサンプリングすることで、積分値を推定することができます。ただし、そのためには任意のサンプルに対する被積分関数の値が評価できる必要があります。レンダリング方程式を解いて左辺の $ L_o $ を推定するために、反射光を表す積分項にモンテカルロ法を適用する事を考えます。しかし積分内には $ L_o $ が存在し、直接評価できないことがわかります。このままではモンテカルロ法を適用できないために、レンダリング方程式の変形を行います。
まず、表記を簡単にするために積分オペレーター $ T $ を定義します。 \begin{eqnarray*} < Tg > \Prta{\vx \RAR \vx'} := \int_{\mathcal{M}} f_s \Prta{\vx'' \RAR \vx \RAR \vx'} g \Prta{\vx'' \RAR \vx} G \Prta{\vx'' \LRAR \vx} \d A \Prta{\vx''} \end{eqnarray*} この積分オペレーターは言わずもがなレンダリング方程式の反射項の積分を表しています。通常、「関数」と言えば数値を受け取りますが、積分オペレーターは関数を受け取ります。これによってレンダリング方程式 $ \eqref{eq:3point_LTE} $ を次のように非常にシンプルな形に書き表す事ができます。 \begin{eqnarray*} L_o = L_e + T L_o \end{eqnarray*} この式の左辺を右辺に代入し、再帰的に計算を行う事で次の級数を得ます。 \begin{eqnarray} L_o = L_e + T L_e + T^2 L_e + T^3 L_e + \cdots = \sum_{i = 0}^{\infty} T^i L_e \label{eq:Neumann_series} \end{eqnarray} この級数をノイマン級数と呼びます。レンダリング方程式のような、左辺に「未知の関数」、右辺に「既知の関数」と「積分に含まれた左辺と同じ未知の関数」によって構成される積分方程式を第二種フレドホルム積分方程式と呼び(※)、その解の一つの表現としてこのノイマン級数があるらしいです。このノイマン級数 $ \eqref{eq:Neumann_series} $ は結果として、ある点を出る放射輝度は、光源から出た光が0, 1, 2, ...回の反射したものの和として表されることを意味しています。そして無限級数とは言え、右辺が既知の関数のみで構成されていることがわかります。
※ レンダリング方程式は第二種フレドホルム積分方程式である、と長らく思われていましたが、実は一般的にはそうではないという指摘をした論文もあるようです。自分は読んでいないので興味のある方は[Soler2022]を読んでみてください。とは言え、現実的にはこのノイマン級数展開に基づいた理論の下、様々なレンダリングアルゴリズムが開発され実用されているので実際は問題にならないことが多いのでしょう。多分。
経路積分による定式化
ノイマン級数 $ \eqref{eq:Neumann_series} $に従い、レンダリング方程式を各経路長(1〜∞)の寄与の和として書き直すことができます。 \begin{equation} L_o \Prta{\vx_1 \RAR \vx_0} = \sum_{k = 1}^{\infty} \underbrace{ \int_{\mathcal{M}} K \Prta{\vx_2, \vx_1, \vx_0} \int_{\mathcal{M}} K \Prta{\vx_3, \vx_2, \vx_1} \cdots \int_{\mathcal{M}} K \Prta{\vx_k, \vx_{k - 1}, \vx_{k - 2}} }_{k - 1個} L_e \Prta{\vx_k \RAR \vx_{k - 1}} \hspace{1mm} \underbrace{ \vphantom{\int_{\mathcal{M}}} \d A_{k} \d A_{k - 1} \ldots \d A_2 }_{k - 1個} \label{eq:path_integral_LTE} \end{equation} ここで、$ K \Prta{\vx''\!, \vx, \vx'} $ はBSDFと幾何項 $ G $ の積を表しており、次のようになっています。 \begin{equation*} K \Prta{\vx''\!, \vx, \vx'} := f_s \Prta{\vx''\! \RAR \vx \RAR \vx'} G \Prta{\vx''\! \LRAR \vx} \end{equation*} この変形によって任意の点 $ \vx_1 $ から $ \vx_0 $ に届く放射輝度が既知の関数 $ L_e $ の単なる多次元積分（とその総和）となりました。モンテカルロ法も適用する事ができそうです。
Measurement Contribution Function
シーン中の光の量を測定するカメラのセンサーのようなものを考えると、その測定値は次に示すmeasurement equation (計測方程式？)によって表されます。
\begin{eqnarray*}
I_j = \int_{\mathcal{M}} \int_{\mathcal{M}} W_e^j \Prta{\vx_1 \RAR \vx_0} L_o \Prta{\vx_1 \RAR \vx_0} G \Prta{\vx_1 \LRAR \vx_0} \d A \Prta{\vx_1} \d A \Prta{\vx_0}
\end{eqnarray*}
ここで $ I_j $ は $ j $ 番目の測定単位、CCDやCMOSセンサーで言うところの画素に相当し、$ W_e^j \Prta{\vx_1 \RAR \vx_0} $ は $ j $ 番目のセンサーの応答を表します。ちなみにセンサー応答はflux responsivityなどとも呼ばれ、単位は $ \Brka{\mathit{S} \cdot \mathrm{W}^{-1}} $ です。ここで $ S $ は任意のセンサー応答の単位で例えば電圧や電流、フィルムの感光度合いなどに対応します。センサーに対してシーン中の各点から入射する放射輝度に応答をかけ、さらにセンサーの領域全体を考慮する事で測定値を計算できます。外側の積分範囲もシーン中の全ての物体表面 $ \mathcal{M} $ として表されていますが、一般的に考えるとセンサー領域外では応答が0であると考えられるため、実質センサー領域に等しくなります。
これに式 $ \eqref{eq:path_integral_LTE} $ を代入すると次のようになります。
\begin{eqnarray*}
I_j &=& \int_{\mathcal{M}} \int_{\mathcal{M}} W_e^j \Prta{\vx_1 \RAR \vx_0} \times\\
&& \hspace{5mm} \sum_{k = 1}^{\infty}
\underbrace{
\int_{\mathcal{M}} K \Prta{\vx_2, \vx_1, \vx_0}
\int_{\mathcal{M}} K \Prta{\vx_3, \vx_2, \vx_1}
\cdots \int_{\mathcal{M}} K \Prta{\vx_k, \vx_{k - 1}, \vx_{k - 2}}
}_{k - 1個}
L_e \Prta{\vx_k \RAR \vx_{k - 1}} \hspace{1mm}
\underbrace{
\vphantom{\int_{\mathcal{S}^2}} \d A_{k} \d A_{k - 1} \ldots \d A_2
}_{k - 1個} \times\\
&& \hspace{10mm} G \Prta{\vx_1 \LRAR \vx_0} \d A \Prta{\vx_1} \d A \Prta{\vx_0}\\
&=& \sum_{k = 1}^{\infty} \underbrace{\int_{\mathcal{M}} \int_{\mathcal{M}} \cdots \int_{\mathcal{M}}}_{k + 1個}
W_e^j \Prta{\vx_1 \RAR \vx_0} G \Prta{\vx_1 \LRAR \vx_0} \prod_{i = 1}^{k - 1} \Brcb{K \Prta{\vx_{i + 1}, \vx_{i}, \vx_{i - 1}}}
L_e \Prta{\vx_k \RAR \vx_{k - 1}} \hspace{1mm} \underbrace{\vphantom{\int_{\mathcal{M}}} \d A_k \d A_{k - 1} \ldots \d A_0}_{k + 1個}
\end{eqnarray*}
多重積分を次の表記を使って簡単に書いてみます。
\begin{eqnarray*}
\int_{\Omega_k} f \Prta{\vx_k \vx_{k - 1} \ldots \vx_0} \d\mu_k \Prta{\vx_k \vx_{k - 1} \ldots \vx_0} :=
\underbrace{\int_{\mathcal{M}} \int_{\mathcal{M}} \cdots \int_{\mathcal{M}}}_{k + 1個}
f \Prta{\vx_k, \vx_{k - 1} \ldots \vx_0} \underbrace{\vphantom{\int_{\mathcal{M}}} \d A_k \d A_{k - 1} \ldots \d A_0}_{k + 1個}
\end{eqnarray*}
ここで $ \Omega_k $ は経路長 $ k $ の経路空間を表しています。測定値の式は次のように若干簡単化されます。
\begin{eqnarray*}
I_j &=& \sum_{k = 1}^{\infty} \int_{\Omega_k}
W_e^j \Prta{\vx_1 \RAR \vx_0} G \Prta{\vx_1 \LRAR \vx_0} \prod_{i = 1}^{k - 1} \Brcb{K \Prta{\vx_{i + 1}, \vx_{i}, \vx_{i - 1}}}
L_e \Prta{\vx_k \RAR \vx_{k - 1}} \hspace{1mm} \d\mu_k \Prta{\vx_k \vx_{k - 1} \ldots \vx_0}
\end{eqnarray*}
さらに経路(Path)の表記 $ \barx_k = \vx_k \vx_{k - 1} \ldots \vx_0 $ を導入すると、被積分関数を次のように表す事ができます。
\begin{eqnarray*}
f_j \Prta{\barx_k} &=& W_e^j \Prta{\vx_1 \RAR \vx_0} G \Prta{\vx_1 \LRAR \vx_0} \prod_{i = 1}^{k - 1} \Brc{f_s \Prta{\vx_{i + 1} \RAR \vx_{i} \RAR \vx_{i - 1}} G \Prta{\vx_{i + 1} \LRAR \vx_{i}}} L_e \Prta{\vx_k \RAR \vx_{k - 1}}
\end{eqnarray*}
この $ f_j(\barx) $ をmeasurement contribution function(明確な日本語訳は無さそう。計測寄与関数？)と呼びます。これによって次の式を得ます。
\begin{eqnarray*}
I_j &=& \sum_{k = 1}^{\infty} \int_{\Omega_k}
f_j \Prta{\barx_k} \hspace{1mm} \d\mu_k \Prta{\barx_k}
\end{eqnarray*}
最後に、あらゆる経路長を含む空間 $ \Omega $ を次のように定義する事でレンダリング方程式をただの積分問題へと導く事ができます。
\begin{equation*}
\Omega = \bigcup_{k = 1}^{\infty} \Omega_k
\end{equation*}
\begin{equation*}
I_j = \int_\Omega f_j \Prt{\barx} \d\mu \Prt{\barx}
\end{equation*}
$ \barx $ はシーン中での経路を表しており、$ f_j \Prta{\barx} $ によってその経路による寄与を表します。そしてそれらを全経路空間で積分する事で $ j $ 番目のセンサーの測定値を計算する事ができます。積分の総和だったものがただの積分になっていますが、数学的には総和は積分の一種として捉えられるそうです。お疲れさまでした。
参考文献
- [Jensen2002] Henrik Wann Jensen(著), 苗村 健(訳) - "フォトンマッピング―実写に迫るコンピュータグラフィックス", 2002, オーム社
- [Soler2022] Cyril Soler, Ronak Molazem, Kartic Subr - "A Theoretical Analysis of Compactness of the Light Transport Operator", 2022
- [Veach1998] Eric Veach, Leonidas J. Guibas - "ROBUST MONTE CARLO METHODS FOR LIGHT TRANSPORT SIMULATION", 1998