フォトンマッピングの理論

このページではフォトンマッピングの理論について解説しています。前提知識としてパストレーシングのアルゴリズムや、モンテカルロ積分, Path Integral LTEなどについて知っているものとします。これらの話題については

memoRANDOM

で詳しく解説されているので、是非見て下さい。

TL;DR

$K$ をカーネル関数として $f(x) \approx \int_{\mathbb{R}} K(x - y)f(y)dy$ と関数近似できる
この近似をPath Integral LTEに適用する
これをモンテカルロ積分で解く

フォトンマッピングとは

Henrik Wann Jensenという方が90年代の終わりに提唱した光輸送アルゴリズムです。光輸送アルゴリズムには パストレーシング(Path Tracing) や 双方向パストレーシング(Bidirectional Path Tracing - BDPT) といった手法がありますが、フォトンマッピングは 集光模様(Caustics) や SDSパス(SDS Path) と呼ばれる光の経路のレンダリングに強いアルゴリズムです。実装すると以下のような水の底でゆらめく光の模様を綺麗にレンダリングすることができます。

SDSパスとはSpecular - Diffuse - Specularとなるような光の経路のことです。例えばガラスに覆われているようなDiffuse面はSDSパスを持ちます。

このような経路はパストレーシングや双方向パストレーシングなどの手法ではSpecular面に挟まれているが故にレンダリングが難しいです。特に光源に点光源、カメラにピンホールカメラを使っているような場合には、このようなパスをサンプリングする確率は0となり、レンダリングが不可能となります。

アルゴリズムの流れ

フォトンマッピングのアルゴリズムの基本的な流れは以下の通りです。

まず、光源から フォトン(Photon) と呼ばれる光のエネルギーを運ぶ仮想の粒子を飛ばし、パストレーシングと同じようにシーン中で反射を繰り返します。この時、Diffuse面に当たる度にフォトンを フォトンマップ(Photon Map) に格納しておきます。

次に視点側からレイトレーシングを行い、Diffuse面に当たったらフォトンマップから近くにあるフォトンをクエリし、反射放射輝度の推定を行います。フォトンには点の位置, 入射方向, Throughputなどの情報が格納されているので、これらを利用して放射輝度の推定を行います。

Specular面に当たった場合にはDiffuse面に当たるまで再帰的にレイを追跡し、Diffuse面でフォトンマップを用いた放射輝度推定を行います。

フォトンマッピングとは光源側から生成したパスをフォトンマップにキャッシュしておき、カメラ側から生成したパスと効率良く繋げる手法と言えます。つまり双方向パストレーシングのような双方向アルゴリズムの1つです。

また、パスをつなげる際には近くにあるフォトンを利用しているため、一定の範囲で "緩く" 接続していると言えます。これによってSDSパスをレンダリングすることが可能になっています。これはパストレーシングや双方向パストレーシングなどのような "点と点を結ぶ" 手法とは大きく異なり、フォトンマッピングを特徴付ける重要な点と言えます。

フォトンマッピングの理論

ここではフォトンマッピングの理論的な部分について

Hachisuka, Toshiya, Jacopo Pantaleoni, and Henrik Wann Jensen. "A path space extension for robust light transport simulation." ACM Transactions on Graphics (TOG) 31.6 (2012): 1-10.

をベースにしながら解説していきます。ここで述べる定式化はJensenのオリジナルのものとは全く異なりますが、フォトンマッピングのアルゴリズムをPath Integral LTEの形で書けるため理論的により明確です。モンテカルロ積分としての解釈も容易なので、後でThroughputの計算を導出する際にも役に立ちます。

カーネル関数近似

ある関数 $f(x)$ の値を、 $x$ 周辺の $f$ の値を使って推定することを考えます(関数補間)。

カーネル関数(Kernel Function) $K(x)$ とは以下のように積分して1になるような関数のことを言います。

\int_{\mathbb{R}} K(x)dx = 1

このような関数を用いると $f(x)$ は以下のように $x$ 周辺の値を使って近似することができます。

f(x) \approx \int_{\mathbb{R}} K(x - y)f(y)dy

これを カーネル関数近似 と呼ぶことにします。

例えば $K$ としてデルタ関数を用いれば

\int_{\mathbb{R}} \delta(x - y)f(y)dy = f(x)

となり, $f(x)$ と一致します。

もちろん実際にはデルタ関数を用いることは出来ないので、以下のように $0$ でピークを持つような関数を用いて推定を行います。この関数は平均0, 分散 $\sigma^2$ の正規分布の密度関数です。このようなカーネルをガウシアンカーネルと呼びます。

K(x) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{1}{2\sigma^2}x^2)

幅の広いカーネルを用いるほど関数近似の誤差は大きくなります。

上で説明したのは1次元の場合でしたが、多次元の場合は

K(x, y) = K(\|x - y\|)

のように点 $x$ と点 $y$ の間の距離を使うことで1次元のカーネルに帰着させて同様に関数近似が行えます。

f(x) \approx \int_{\mathbb{R^d}} K(x, y)f(y)dy

3点形式のLTEのカーネル関数近似

以下のような3点形式のLTEに対してカーネル関数近似を適用することを考えます。

$L_o(x_1 \to x_2)$ : $x_1$ から $x_2$ に出射する放射輝度
$\mathcal{M}$ : シーン上の点全体
$L_i(x_0 \to x_1)$ : $x_0$ から $x_1$ に入射する放射輝度
$f$ : BSDF
$G$ : 幾何項

L_o(x_1 \to x_2) = \int_{\mathcal{M}} L_i(x_0 \to x_1)f(x_0 \to x_1 \to x_2)G(x_0, x_1)dA(x_0)

今、 $x_1$ で $L_i$ , $f$ , $G$ の値を知ることができないとします。この場合 $x_1$ でカーネル関数近似を適用すると

L_o(x_1 \to x_2) = \int_{\mathcal{M}}\int_{\mathcal{M}} K(y, x_1)L_i(x_0 \to y)f(x_0 \to y \to x_2)G(x_0, y)dA(x_0)dA(y) \tag{1}

となり、 $x_1$ 周辺の $y$ の位置での $L_i$ や $f$ の値を使って推定することができます。これが3点形式のLTEにカーネル関数近似を適用した形です。

3点形式のLTEのカーネル関数近似(緩く接続した場合)

式(1)では $x_1$ を全て $y$ に置き換えてカーネル関数近似を行っていました。全て置き換えるのではなく、部分的に置き換える場合を考えます。

L_o(x_1 \to x_2) = \int_{\mathcal{M}}\int_{\mathcal{M}} K(y, x_1)L_i(x_0 \to y)f(x_0 \to y, x_1 \to x_2)G(x_0, y)dA(x_0)dA(y)

ここでは $f$ の中身を $x_0 \to y \to x_2$ から $x_0 \to y, x_1 \to x_2$ に変えています。図で表すと以下のような形です。

このパスは $y$ と $x_1$ を無理やり繋げた形になっています。これは物理的に正確な意味でのパスではありません(地面に沿って光が進んでいるということになる)。そのような意味で"緩く"接続されたパスと言えます。

何故このような緩い接続を考えるかと言うと、冒頭でも述べた通り、フォトンマッピングではこの緩い接続がSDSパスなどをレンダリングするために重要であり、これを理論的に説明するためにこのような緩いカーネル関数近似を考えます。

Path Integral LTEのカーネル関数近似

続いてカーネル関数近似をPath Integral LTEに適用することを考えます。簡単のため経路長3のパスのPath Integral LTEを考えます。

$I_j$ : $j$ 番目のピクセルの寄与
$\mathcal{M}$ : シーン上の点全体
$L_e$ : 光源から出る放射輝度
$G$ : 幾何項
$f$ : BSDF
$W_e^{j}$ : $j$ 番目のピクセルのセンサー応答

I_j = \int_{\mathcal{M}^4}L_e(x_0 \to x_1)G(x_0, x_1)f(x_0 \to x_1, \to x_2)G(x_1, x_2)f(x_1 \to x_2 \to x_3)G(x_2, x_3)W_e^{j}(x_2 \to x_3)dA(x_0)dA(x_1)dA(x_2)dA(x_3)

今, $x_2$ で緩い接続の意味でのカーネル関数近似を適用するとします。式は以下のようになります。

I_j = \int_{\mathcal{M}^5} L_e(x_0 \to x_1)G(x_0, x_1)f(x_0 \to x_1, \to y)G(x_1, y)K(y, x_2)f(x_1 \to y, x_2 \to x_3)G(x_2, x_3)W_e^{j}(x_2 \to x_3)dA(x_0)dA(x_1)dA(y)dA(x_2)dA(x_3)

この式は $x_2$ においてフォトンマップを用いて放射輝度推定を行っていることに対応します(フォトンの位置は $y$ )。

モンテカルロ積分

先程のPath Integral LTEの緩いカーネル関数近似をモンテカルロ積分することを考えます。

今、光源側からフォトンを飛ばして $x_0, x_1, y$ と経路を構築し、カメラ側からパストレーシングによって $x_3, x_2$ と経路を構築したとします。この時、先程の緩いカーネル関数近似を行ったPath Integral LTEのモンテカルロ積分は以下のように書けます。

$N$ : 経路生成のサンプル数
$p(x)$ : 点 $x$ を生成する確率密度
$p(x_i|x_j)$ : 点 $x_j$ から点 $x_i$ を生成する時の条件付き確率密度

I_j \approx \frac{1}{N}\sum_{i=1}^N \frac{L_e(x_0^i \to x_1^i)G(x_0^i, x_1^i)f(x_0^i \to x_1^i, \to y^i)G(x_1^i, y^i)K(y^i, x_2^i)f(x_1^i \to y^i, x_2^i \to x_3^i)G(x_2^i, x_3^i)W_e^{j}(x_2^i \to x_3^i)}{p(x_0^i)p(x_1^i|x_0^i)p(y^i|x_1^i)p(x_2^i|x_3^i)p(x_3^i)} \tag{2}

ここで上付きの $i$ は $i$ 番目のサンプルを表しています。この式は各 $i$ で毎回光源側とカメラ側からランダムに経路構築を行い、 $y$ と $x_2$ を緩く接続するようなアルゴリズムを表現しています。

フォトンのThroughputを

\beta_i = \frac{L_e(x_0^i \to x_1^i)G(x_0^i, x_1^i)f(x_0^i \to x_1^i, \to y^i)G(x_1^i, y^i)}{p(x_0^i)p(x_1^i|x_0^i)p(y^i|x_1^i)}

カメラ側からのレイのThroughputを

T_i = \frac{G(x_2^i, x_3^i)W_e^{j}(x_2^i \to x_3^i)}{p(x_2^i|x_3^i)p(x_3^i)}

とすると、式(2)は以下のように書けます。

I_j \approx \frac{1}{N}\sum_{i=1}^N\beta_i K(y^i, x_2^i)f(x_1^i \to y^i, x_2^i \to x_3^i)T_i

フォトンマッピング

式(2)は各 $i$ で光源側とカメラ側から毎回別々の経路を構築するような式になっていました。フォトンマッピングでは冒頭で述べた通り、フォトンマップを使うことで予め光源側からのパスの生成を行ってキャッシュしておきます。

フォトンマップ作成に使うフォトン数を $N_p$ 、カメラ側からの経路生成のサンプル数を $N$ とすると、式(2)は以下のように書き換えることができます。

I_j \approx \frac{1}{N}\sum_{i=1}^N \bigg(\frac{1}{N_p}\sum_{p=1}^{N_p}\beta_pK(y^p, x_2^i)f(x_1^p \to y^p, x_2^i \to x_3^i)\bigg)T_i \tag{3}

ここで $()$ で囲った部分がフォトンマップを用いた放射輝度推定に対応する部分です。フォトンマッピングの説明でよく 密度推定(Density Estimation) と言われている部分に対応します。(以上の議論の流れを見れば分かる通り、密度推定というよりは関数近似と呼ぶほうが適切なように思います)

$T_i$ の部分はカメラ側からのパストレーシングの部分に対応します。これを $()$ の項と掛け合わせることで $j$ 番目のピクセルの値が推定できます。

式(3)では $()$ 内の総和を全てのフォトンに対して行っていますが、実際にはカーネル関数は点 $x_2$ から離れたところでは0または非常に小さな値を取ります。例えばカーネル関数として

K(x, y) = \begin{cases} \frac{1}{\pi r^2} & (\|x-y\| \le r) \\ 0 & (otherwise) \end{cases}

という半径 $r$ の円上で定数となるようなカーネルを使うと、点 $x_2$ から半径 $r$ 以上の位置にあるフォトンに対して総和を取っても、カーネルの値が0になるため無駄です。そこで実際には $x_2$ の近傍にあるフォトンのみを用いて総和を取ります。このためフォトンマップを表現するデータ構造としては近傍探索が効率良く行えるOctreeやKd-treeなどのデータ構造が選ばれます。

式(3)がフォトンマッピングの理論的な部分を全て説明する式です。今回は経路長3の場合で考えましたが、他の経路長の場合でも同様の式を得ることができます。

次回は今回の内容を元に、フォトンマッピングの実装について説明します。

次回

References

Jensen, Henrik Wann. Realistic image synthesis using photon mapping. AK Peters/crc Press, 2001.
https://pbr-book.org/3ed-2018/Light_Transport_III_Bidirectional_Methods/Stochastic_Progressive_Photon_Mapping#
Jensen, Henrik Wann. "Global illumination using photon maps." Eurographics workshop on Rendering techniques. Springer, Vienna, 1996.
Veach, Eric. Robust Monte Carlo methods for light transport simulation. Stanford University, 1998.
Hachisuka, Toshiya, Jacopo Pantaleoni, and Henrik Wann Jensen. "A path space extension for robust light transport simulation." ACM Transactions on Graphics (TOG) 31.6 (2012): 1-10.