このページではMathJaxによる数式表示を行っています。対応しているブラウザーであれば少し待てば数式が表示されます。
(画面左下にMathJaxのタイプセットステータスが表示されます。)
環境によっては表示が崩れているかもしれません。その際はブラウザーを変更してみてください。

双方向パストレーシング実装

双方向パストレーシング」ではアルゴリズムの理論面を紹介しました。このページでは実装に関して擬似コードを使って解説します。

サブパスの生成

光源 $ \mathrm{light} $ を選択、選ぶ確率 $ P_\mathrm{light} $
$ \mathrm{light} $ から位置 $ \by_0 $ をサンプル、面積に関する確率密度 $ p_A(\by_0) $
光源点の情報を取得: $ \by_0 $ における法線 $ \vec{n}_{\by_0} $、EDF $ L_e^{(1)}(\by_0, \cdot) $
1点目のウェイトを計算 $ \displaystyle \alpha_1^L = \frac{L_e^{(0)}(\by_0)}{P_\mathrm{light} \cdot p_A(\by_0)} $
サブパスに頂点を追加 $ \mathrm{lightSubPath.append} ( \by_0, \vec{n}_{\by_0}, \Vec{0}, L_e^{(1)}(\by_0, \cdot), \alpha_1^L, p_A(\by_0), 1 ) $
\rh{30}
光の出射方向 $ \vomega_{\by_0 \by_1} $ をサンプル、方向に関する確率密度 $ p_\sigma(\vomega_{\by_0 \by_1}) $
cos項を計算 $ \mathrm{cosLast} \LAR \abs{\vomega_{\by_0 \by_1} \cdot \vec{n}_{\by_0}} $
2点目のウェイトを計算 $ \displaystyle \alpha_2^L = \frac{L_e^{(1)}(\by_0, \vomega_{\by_0 \by_1}) \cdot \mathrm{cosLast}}{p_\sigma(\vomega_{\by_0 \by_1})} \alpha_1^L $
レイを生成 $ \mathrm{ray} \LAR \{ \mathrm{org:\;} \by_0,\; \mathrm{dir:\;} \vomega_{\by_0 \by_1} \} $
続きのサブパスを生成 $ \mathrm{generateSubPath} ( \alpha_2^L, \mathrm{ray}, p_\sigma(\vomega_{\by_0 \by_1}), \mathrm{cosLast}, \mathrm{lightSubPath} ) $

まずは光線サブパスの生成手順について説明します。

1行目 光源の数は往往にして複数個になりますが、基本的には視線サブパス1本に対して光線サブパスも1本とするので、ある一つの光源 $ \mathrm{light} $ を確率 $ P_\mathrm{light} $ で選択します。それぞれの光源を等しい確率で選ぶ場合には $ P_\mathrm{light} $ は 1 / (光源数) となりますが、各光源の出力に応じて重点的サンプリングを行う場合には、それを考慮した値にする必要があります。
2行目 選んだ光源 $ \mathrm{light} $ の発光面上の位置 $ \by_0 $ をサンプリングします。面積測度に関する確率密度 $ p_A(\by_0) $ も得ます。
3行目 光源上のサンプル点 $ \by_0 $ における各種情報を得ます。EDFとは「無限遠光源の扱い方」でも示したように、放射発散度と放射輝度の関係を表す量です(便宜的に使用している名称で一般的ではありません)。
4行目 1点目のウェイト $ \alpha_1^L $ を計算します。基本的には理論編のページに示した通りですが、光源を選ぶ確率 $ P_\mathrm{light} $ も分母に含まれていることに注意してください。分母の積は、全ての光源をまとまったひとつの光源と考えて、その面上の位置をサンプリングする確率密度と捉えることもできます。それどころか、シーン中の全ての物体表面をひとつの光源として考えることもできます。その場合、非光源面は「輝度ゼロで光っている」と捉え直すことで、重点的サンプリングの観点で、本当に光っている部分だけをサンプリングしているということになります。
5行目 光線サブパスの1点目の情報を記録します。サブパスの頂点情報に記録するものとして、位置や法線、トレースしてきた方向、BSDF、ウェイト、サブパスを生成する確率密度、この頂点を作る時のロシアンルーレットの確率などが挙げられます。トレースしてきた方向というものは1点目には存在せず、後の計算で使われることもないので適当な方向(ここではゼロベクトル $ \Vec{0} $)を記録します。BSDFは、サブパス接続時に接続辺の両端のBSDFや確率密度の値を評価するために記録します。そのため何らかの入射・出射方向に関するBSDFの評価値ではなく、関数そのものの情報を記録します。ただし、1点目に関してはBSDFではなく、それに相当するEDF $ L_e^{(1)}(\by_0, \cdot) $ を記録します。ロシアンルーレットの確率も基本的には1点目で考えることが無く、1点目は100%生成されるので $ 1 $ を記録します。BSDFが異方性を持つ場合には接平面ベクトル2つも記録する必要があり、状況に応じて必要な情報を全て記録します。
7行目 光源上のサンプル位置 $ \by_0 $ から光の出射方向 $ \vomega_{\by_0 \by_1} $ をサンプリングします。立体角測度に関する確率密度 $ p_\sigma(\vomega_{\by_0 \by_1}) $ も得ます。 位置のサンプルも含めて、無限遠光源の場合の確率密度は特殊な扱い方を必要とします。「無限遠光源の扱い方」をご参照ください。
8行目 後の計算で使用するcos項を計算します。
9行目 2点目のウェイト $ \alpha_2^L $ を計算します。
10行目 サンプリングした光源上の位置 $ \by_0 $ と、光の出射方向 $ \vomega_{\by_0 \by_1} $ からレイを生成します。
11行目 光線・視線サブパス共通の処理に続きます。
レンズ上の位置 $ \bz_0 $ をサンプル、面積に関する確率密度 $ p_A(\bz_0) $
レンズ上の点の情報を取得: $ \bz_0 $ における法線 $ \vec{n}_{\bz_0} $、IDF $ W_e^{(1)}(\bz_0, \cdot) $
1点目のウェイトを計算 $ \displaystyle \alpha_1^E = \frac{W_e^{(0)}(\bz_0)}{p_A(\bz_0)} $
サブパスに頂点を追加 $ \mathrm{eyeSubPath.append} ( \bz_0, \vec{n}_{\bz_0}, \Vec{0}, W_e^{(1)}(\bz_0, \cdot), \alpha_1^E, p_A(\bz_0), 1 ) $
\rh{30}
光の入射方向 $ \vomega_{\bz_0 \bz_1} $ をサンプル、方向に関する確率密度 $ p_\sigma(\vomega_{\bz_0 \bz_1}) $
cos項を計算 $ \mathrm{cosLast} \LAR \abs{\vomega_{\bz_0 \bz_1} \cdot \vec{n}_{\bz_0}} $
2点目のウェイトを計算 $ \displaystyle \alpha_2^E = \frac{W_e^{(1)}(\bz_0, \vomega_{\bz_0 \bz_1}) \cdot \mathrm{cosLast}}{p_\sigma(\vomega_{\bz_0 \bz_1})} \alpha_1^E $
レイを生成 $ \mathrm{ray} \LAR \{ \mathrm{org:\;} \bz_0,\; \mathrm{dir:\;} \vomega_{\bz_0 \bz_1} \} $
続きのサブパスを生成 $ \mathrm{generateSubPath} ( \alpha_2^E, \mathrm{ray}, p_\sigma(\vomega_{\bz_0 \bz_1}), \mathrm{cosLast}, \mathrm{eyeSubPath}) $

続いて視線サブパスの生成手順について説明します。おおよそ光線サブパスと同じ手順です。

1行目 レンズ(もしくはセンサー)面上の位置 $ \bz_0 $ をサンプリングします。面積測度に関する確率密度 $ p_A(\bz_0) $ も得ます。
2行目 レンズ面上のサンプル点 $ \bz_0 $ における各種情報を得ます。IDF(Importance Distribution Function)はEDFと同様に、センサー応答の位置に関する成分と、方向も含めたセンサー応答の関係を表す量です(これも便宜的に定義した名称で一般的ではありませんので注意)。
3行目 1点目のウェイト $ \alpha_1^E $ を計算します。
4行目 視線サブパスの1点目の情報を記録します。光線サブパスの1点目と同様に、トレースしてきた方向には適当な値、ロシアンルーレットの確率には $ 1 $ を記録します。BSDFの代わりとしてIDF $ W_e^{(1)}(\bz_0, \cdot) $ を記録します。
6行目 レンズ面上のサンプル位置 $ \bz_0 $ から光の入射方向 $ \vomega_{\bz_0 \bz_1} $ をサンプリングします。立体角測度に関する確率密度 $ p_\sigma(\vomega_{\bz_0 \bz_1}) $ も得ます。この確率密度の求め方については「被写界深度」をご参照ください。
7行目 後の計算で使用するcos項を計算します。
8行目 2点目のウェイト $ \alpha_2^E $ を計算します。
9行目 サンプリングしたレンズ面上の位置 $ \bz_0 $ と、光の入射方向 $ \vomega_{\bz_0 \bz_1} $ からレイを生成します。
10行目 光線・視線サブパス共通の処理に続きます。
function $ \mathrm{generateSubPath} $( パスのウェイト $ \alpha $, レイ $ \mathrm{ray} $, 立体角測度の確率密度 $ p_{\sigma, \mathrm{last}} $, 最後のcos項 $ \mathrm{cosLast} $, サブパス $ \mathrm{subPath} $ ) {
	最後のロシアンルーレットの確率 $ P_{\mathrm{RR}} \LAR 1 $
	while レイ $ \mathrm{ray} $ をトレースして衝突がある {
		衝突点の情報を取得: 位置 $ \vx $, 法線 $ \vec{n} $, トレースしてきた方向 $ \vomega_\mathrm{traced} \LAR -\mathrm{ray.dir} = $ ($ \vomega_o $ or $ \vomega_i $), BSDF $ f_s(\cdot) $
		\rh{30}
		最後の頂点との距離の2乗 $ \mathrm{dist2} \LAR \abs{\vx - \mathrm{subPath.lastVertex.position}}^2 $
		現在の点 $ \vx $ をサンプリングする面積に関する確率密度 $ \displaystyle p_{A, \vx} \LAR \frac{p_{\sigma, \mathrm{last}} \cdot \abs{\vomega_\mathrm{traced} \cdot \vec{n}}}{\mathrm{dist2}} $
		サブパスに頂点を追加 $ \mathrm{subPath.append} ( \vx, \vec{n}, \vomega_\mathrm{traced}, f_s(\cdot), \alpha, p_{A, \vx}, P_{RR} ) $
		\rh{30}
		if 視線サブパストレース中に現在の点が光源 or 光線サブパストレース中に現在の点がレンズの場合 {
			$ \ldots $
		}
		\rh{30}
		次のトレース方向をサンプル $ \vomega_\mathrm{sampled} = $ ($\vomega_i $ or $ \vomega_o $)、方向に関する確率密度 $ p_\sigma(\vomega_\mathrm{sampled}) $, 逆方向の確率密度 $ p_\sigma(\vomega_\mathrm{traced}) $
		ウェイトを更新 $ \displaystyle \alpha \LAR \alpha \times \frac{f_s(\vx, \vomega_i, \vomega_o) \abs{\vomega_\mathrm{sampled} \! \cdot \vec{n}}}{p_\sigma(\vomega_\mathrm{sampled})} $
		\rh{30}
		ロシアンルーレットの確率を決める $ P_{RR} \LAR \ldots $
		if $ [0, 1) $ の一様乱数 $ \xi $ に関して $ \xi \ge P_{RR} $
			break
		\rh{30}
		新たなレイを生成 $ \mathrm{ray} \LAR \{ \mathrm{org:\;} \vx,\; \mathrm{dir:\;} \vomega_\mathrm{sampled} \} $
		ウェイトにロシアンルーレットの確率を反映 $ \displaystyle \alpha \LAR \alpha \times \frac{1}{P_{RR}} $
		\rh{30}
		最後から2番目の頂点を逆方向からサンプリングする場合の値を記録
		面積に関する確率密度 $ \displaystyle \mathrm{subPath.secondLastVertex.revAreaPDF} \LAR \frac{p_\sigma(\vomega_\mathrm{traced}) \cdot \mathrm{cosLast}}{\mathrm{dist2}} $
		ロシアンルーレットの確率 $ \displaystyle \mathrm{subPath.secondLastVertex.revRRProb} \LAR \ldots $
		\rh{30}
		最後の確率密度・cos項を更新 $ p_{\sigma, \mathrm{last}} \LAR p_\sigma(\vomega_\mathrm{sampled}),\; \mathrm{cosLast} \LAR \abs{\vomega_\mathrm{sampled} \cdot \vec{n}} $
	}
}

光線・視線サブパスの構築は、いずれも3点目以降はBSDFのサンプリングとなりほぼ共通の処理となるため、ひとつの関数にまとめることが考えられます。おおまかな流れとしてはパストレーシングと似ています。

2行目 光線・視線サブパスともに2点目のトレース前にはまだロシアンルーレットは使わないと思うので、確率としては 1 になります。
3行目 経路をBSDFのサンプリングとレイトレースを繰り返しながら延長します。
4行目 衝突点 $ \vx $ の各種情報を得ます。レイトレースしてきた方向 $ \vomega_\mathrm{traced} $ は、光線サブパスの場合は光の入射方向 $ \vomega_i $、視線サブパスの場合は光の出射方向 $ \vomega_o $ となります。
6行目 後の計算で使うために、現在の衝突点 $ \vx $ とサブパスの最後の頂点間の距離の2乗を計算します。いずれかの点が無限遠光源上の場合、特殊な処理が必要となります。
7行目 レイと法線からなるcos項、距離の2乗を用いて現在の衝突点 $ \vx $ をサンプリングする確率密度を立体角測度から面積測度に変換します。
8行目 サブパスの頂点情報を記録します。1点目と同様に後のサブパス接続時に必要となる量を記録します。
10-12行目 光線サブパスをトレース中にレンズ面にヒットした場合、視線サブパスをトレース中に光源面にヒットした場合は、サブパスだけで完全な光輸送経路を構築できたことになるので、MISウェイトの計算と寄与の蓄積が行えます。ここの処理については、後ほどMISウェイト計算について説明した後に改めて説明します。
14行目 次のトレース方向をBSDFからサンプリングします。サンプリングされる方向 $ \vomega_\mathrm{sampled} $ は、光線サブパスの場合は光の出射方向 $ \vomega_o $、視線サブパスの場合は光の入射方向 $ \vomega_i $ となります。サンプリング方向に関する立体角測度の確率密度 $ p_\sigma(\vomega_\mathrm{sampled}) $ を求めますが、加えてトレースしてきた方向に関する立体角測度の確率密度 $ p_\sigma(\vomega_\mathrm{traced}) $ も求めます。
例として光線サブパスのトレースを考えます。サンプリング時に与えられる「トレースしてきた方向」は光の入射方向 $ \vomega_i $となり、出射方向 $ \vomega_o $ をサンプリング、確率密度を求めますが、同時に「トレースしてきた方向」を $ \vomega_o $ とした場合に $ \vomega_i $ をサンプリングする確率密度も求めます。この逆方向の確率密度は、「逆方向を順方向として」トレースするときの確率密度と一致する必要があります。つまり、光線サブパストレース時に求める逆方向の確率密度は、視線サブパストレース時に2つの方向を入れ替えて与えた場合に得られる順方向の確率密度に一致する必要があります。(何を言っているかわかりにくいと思いますが、MISウェイト計算で使用する各戦略の確率密度を矛盾なく導くことを考えればわかると思います。)
逆方向にサブパスを構築する際の確率密度をあらかじめ求めておくことで、後のサブパス接続時に行うMISウェイトの計算が効率的になります。
15行目 サブパスのウェイト $ \alpha $ を更新します。
17-19行目 ロシアンルーレットの確率 $ P_{RR} $ を決定してサブパスの延長を確率的に打ち切ります。
21行目 新たなレイを生成します。
22行目 ロシアンルーレットの確率 $ P_{RR} $ をサブパスのウェイトに反映します。
24-26行目 ロシアンルーレットで打ち切られることなく、サンプリングした方向へのサブパスの延長が確定したことで、最後から2番目に記録されているサブパスの頂点を逆方向からサンプリングする場合の確率密度も確定します。サンプリング時に求めておいた逆方向の確率密度 $ p_\sigma(\vomega_\mathrm{traced}) $ を面積測度に変換して記録します。また、逆方向にトレースした場合のロシアンルーレットの確率も記録します。逆方向のロシアンルーレットの確率も確率密度と同じく、順方向・逆方向で辻褄を合わせる必要があります。
28行目 確率密度・cos項を更新します。

サブパスの接続とMISウェイトの計算

for $ t \LAR 1 $ to 視線サブパスの頂点数 {
	視線サブパスの $ t $ 頂点目を取得 $ \mathrm{eyeVtx} \LAR \mathrm{eyeSubPath.getNthVertex(} t \mathrm{)} $
	for $ s \LAR 1 $ to 光線サブパスの頂点数 {
		光線サブパスの $ s $ 頂点目を取得 $ \mathrm{lightVtx} \LAR \mathrm{lightSubPath.getNthVertex(} s \mathrm{)} $
		\rh{30}
		接続辺のベクトルと距離の2乗を計算
		$ \vomega_\mathrm{connect} \LAR \mathrm{normalize(lightVtx.position - eyeVtx.position)} $, $ \mathrm{dist2_{connect}} \LAR \abs{\mathrm{lightVtx.position - eyeVtx.position}}^2 $
		光線サブパスの端点のcos項を計算 $ \mathrm{cosLightEnd} = \abs{\vomega_\mathrm{connect} \cdot \mathrm{lightVtx.normal}} $
		視線サブパスの端点のcos項を計算 $ \mathrm{cosEyeEnd} = \abs{\vomega_\mathrm{connect} \cdot \mathrm{eyeVtx.normal}} $
		接続辺の幾何項を計算 $ \displaystyle G_\mathrm{connect} = \frac{\mathrm{cosLightEnd} \cdot \mathrm{cosEyeEnd}}{\mathrm{dist2_{connect}}} $
		\rh{30}
		光線サブパスの頂点のBSDFを取得 $ f_{s, \mathrm{lightEnd}}(\cdot) \LAR \mathrm{lightVtx.function} $
		BSDFの値を計算 $ f_{s, \mathrm{lightEnd}} = f_{s, \mathrm{lightEnd}}(\mathrm{lightVtx.position}, \mathrm{lightVtx.tracedDir}, -\vomega_\mathrm{connect}) $
		$ \mathrm{lightVtx.tracedDir} $ 方向からトレースして $ -\vomega_\mathrm{connect} $ 方向にサンプリングする場合と、逆方向の場合の確率密度を計算
		$ p_{\sigma, \mathrm{extendLight, 1st}} \LAR p_\sigma(-\vomega_\mathrm{connect}) $, $ p_{\sigma, \mathrm{extendEye, 2nd}} \LAR p_\sigma(\mathrm{lightVtx.tracedDir}) $
		\rh{30}
		視線サブパスの頂点のBSDFを取得 $ f_{s, \mathrm{eyeEnd}}(\cdot) \LAR \mathrm{eyeVtx.function} $
		BSDFの値を計算 $ f_{s, \mathrm{eyeEnd}} = f_{s, \mathrm{eyeEnd}}(\mathrm{eyeVtx.position}, \vomega_\mathrm{connect}, \mathrm{eyeVtx.tracedDir}) $
		$ \mathrm{eyeVtx.tracedDir} $ 方向からトレースして $ \vomega_\mathrm{connect} $ 方向にサンプリングする場合と、逆方向の場合の確率密度を計算
		$ p_{\sigma, \mathrm{extendEye, 1st}} \LAR p_\sigma(\vomega_\mathrm{connect}) $, $ p_{\sigma, \mathrm{extendLight, 2nd}} \LAR p_\sigma(\mathrm{eyeVtx.tracedDir}) $
		\rh{30}
		放射束の伝達度合いを計算 $ c_{s, t} \LAR f_{s, \mathrm{lightEnd}} \, G_\mathrm{connect} \, f_{s, \mathrm{eyeEnd}} $
		\rh{20}
        if $ \mathrm{lightVtx} $ と $ \mathrm{eyeVtx} $ の間に遮蔽物がある or $ c_{s, t} $ がゼロ
        	continue
        \rh{40}
        光線サブパスの1度目の延長における確率密度の測度を変換、ロシアンルーレットの確率も求める
        $ \displaystyle p_{A, \mathrm{extendLight, 1st}} = \frac{p_{\sigma, \mathrm{extendLight, 1st}} \cdot \mathrm{cosEyeEnd}}{\mathrm{dist2_{connect}}} $
        $ P_{RR, \mathrm{extendLight, 1st}} \LAR \ldots $
        if $ t > 1 $
        	視線サブパスの接続頂点の1つ前の頂点を取得 $ \mathrm{eyeVtx2nd} \LAR \mathrm{eyeSubPath.getNthVertex(} t - 1 \mathrm{)} $
        	光線サブパスの2度目の延長における確率密度の測度を変換、ロシアンルーレットの確率も求める
        	$ \vomega_{\mathrm{2nd}} \LAR \mathrm{normalize(eyeVtx2nd.position - eyeVtx.position)} $, $ \mathrm{dist2} \LAR \abs{\mathrm{eyeVtx2nd.position - eyeVtx.position}}^2 $
        	$ \displaystyle p_{A, \mathrm{extendLight, 2nd}} = \frac{p_{\sigma, \mathrm{extendLight, 2nd}} \cdot \abs{\vomega_\mathrm{2nd} \cdot \mathrm{eyeVtx2nd.normal}}}{\mathrm{dist2}} $
        	$ P_{RR, \mathrm{extendLight, 2nd}} \LAR \ldots $
        \rh{40}
        視線サブパスの1度目の延長における確率密度の測度を変換、ロシアンルーレットの確率も求める
        $ \displaystyle p_{A, \mathrm{extendEye, 1st}} = \frac{p_{\sigma, \mathrm{extendEye, 1st}} \cdot \mathrm{cosLightEnd}}{\mathrm{dist2_{connect}}} $
        $ P_{RR, \mathrm{extendEye, 1st}} \LAR \ldots $
        if $ s > 1 $
        	光線サブパスの接続頂点の1つ前の頂点を取得 $ \mathrm{lightVtx2nd} \LAR \mathrm{lightSubPath.getNthVertex(} s - 1 \mathrm{)} $
        	視線サブパスの2度目の延長における確率密度の測度を変換、ロシアンルーレットの確率も求める
        	$ \vomega_{\mathrm{2nd}} \LAR \mathrm{normalize(lightVtx2nd.position - lightVtx.position)} $, $ \mathrm{dist2} \LAR \abs{\mathrm{lightVtx2nd.position - lightVtx.position}}^2 $
        	$ \displaystyle p_{A, \mathrm{extendEye, 2nd}} = \frac{p_{\sigma, \mathrm{extendEye, 2nd}} \cdot \abs{\vomega_\mathrm{2nd} \cdot \mathrm{lightVtx2nd.normal}}}{\mathrm{dist2}} $
        	$ P_{RR, \mathrm{extendEye, 2nd}} \LAR \ldots $
        \rh{40}
        MISウェイトを計算 $ \begin{aligned}[t] w_{s, t} \LAR \mathrm{calculateMISWeight} (& \mathrm{lightSubPath}, \mathrm{eyeSubPath}, s, t, \\& p_{A, \mathrm{extendLight, 1st}}, P_{RR, \mathrm{extendLight, 1st}}, p_{A, \mathrm{extendLight, 2nd}}, P_{RR, \mathrm{extendLight, 2nd}}, \\& p_{A, \mathrm{extendEye, 1st}}, P_{RR, \mathrm{extendEye, 1st}}, p_{A, \mathrm{extendEye, 2nd}}, P_{RR, \mathrm{extendEye, 2nd}} ) \end{aligned} $
        寄与を計算 $ C_{s, t} \LAR w_{s, t} \times \mathrm{lightVtx.alpha} \times c_{s, t} \times \mathrm{eyeVtx.alpha} $
        \rh{30}
        if $ t > 1 $
        	処理中のピクセル位置 $ j $ に寄与を蓄積する $ I_{eye, j} \LAR I_{eye, j} + C_{s, t} $
        else
        	視線サブパスの接続頂点 $ \mathrm{eyeVtx} $ と接続辺のベクトル $ \vomega_\mathrm{connect} $ からピクセル位置 $ j^* $ を逆算する
        	逆算したピクセル位置 $ j^* $ に寄与を蓄積する $ I_{light, j^*} \LAR I_{light, j^*} + C_{s, t} $
	}
}

光線サブパス・視線サブパス、両方を構築し終わった後、それぞれの頂点間の接続・MISウェイトの計算を行い、寄与の蓄積を行なっていきます。

1-2行目 視線サブパスの各頂点分ループします。サブパスの $ t $ 頂点目 $ \mathrm{eyeVtx} $ を取得します。
3-4行目 光線サブパスの各頂点分ループします。サブパスの $ s $ 頂点目 $ \mathrm{lightVtx} $ を取得します。
6-7行目 接続頂点間の方向 $ \vomega_\mathrm{connect} $ と距離の2乗 $ \mathrm{dist2_{connect}} $ を計算します。
8-10行目 それぞれのサブパスの接続頂点におけるcos項と距離2乗を用いて接続の幾何項 $ G_\mathrm{connect} $ を計算します。
12-15行目 光線サブパスの頂点におけるBSDFの値と、光線サブパスを延長したと考える場合の確率密度 $ p_{\sigma, \mathrm{extendLight, 1st}} $ を計算します。確率密度はサブパス頂点に記録されているトレース方向 $ \mathrm{lightVtx.tracedDir} $ と、接続辺の方向 $ -\vomega_\mathrm{connect} $ から求めることができます。同時に2つの方向を入れ替えた場合の確率密度も計算しておきます。この確率密度は、サブパス接続後に視線サブパスの2度目の延長を考える場合の確率密度に相当します。なお、$ s = 1 $ の場合はBSDFではなく、EDFとなります。
17-20行目 視線サブパスの頂点におけるBSDFの値と、視線サブパスを延長したと考える場合の確率密度 $ p_{\sigma, \mathrm{extendEye, 1st}} $ を計算します。確率密度はサブパス頂点に記録されているトレース方向 $ \mathrm{eyeVtx.tracedDir} $ と、接続辺の方向 $ \vomega_\mathrm{connect} $ から求めることができます。同時に2つの方向を入れ替えた場合の確率密度も計算しておきます。この確率密度は、サブパス接続後に光線サブパスの2度目の延長を考える場合の確率密度に相当します。なお、$ t = 1 $ の場合はBSDFではなく、IDFとなります。
22行目 サブパス間の放射束の伝達度合いを表す $ c_{s, t} $ を計算します。
24-25行目 サブパス頂点間に遮蔽物がある場合や、$ c_{s, t} $ がゼロとなる場合には、この接続の処理をスキップします。
27-35行目 上で求めた光線サブパスを延長する場合の確率密度2つを面積測度に変換、2度の延長それぞれのロシアンルーレットの確率も求めます。
37-45行目 上で求めた視線サブパスを延長する場合の確率密度2つを面積測度に変換、2度の延長それぞれのロシアンルーレットの確率も求めます。
47行目 2つのサブパス、現在のそれぞれの頂点インデックス $ s, t $、そして接続に依存して決まる確率密度を使ってMISウェイトを計算します。MISウェイト計算には理論編で示したように、各戦略の確率密度の比だけが必要となりますが、その計算に必要な値は各頂点に記録されている値と、ここで追加で渡した確率密度で全て揃います。
48行目 MISウェイトをかけた寄与を計算します。
50-54行目 $ t > 1 $ の戦略の場合には素直に、処理中のピクセル位置に寄与を蓄積するだけですが、 $ t \le 1 $ の戦略の場合にはまずピクセル位置を現在の接続から逆算する必要があります。
function $ \mathrm{calculateMISWeight} $( 光線・視線サブパス $ \mathrm{lightSubPath} $, $ \mathrm{eyeSubPath} $, それぞれのサブパスの頂点数 $ \mathrm{numLightVtx} $, $ \mathrm{numEyeVtx} $
									光線サブパスの1, 2度目の延長に関する確率密度 $ p_{A, \mathrm{extendLight, 1st}}, P_{RR, \mathrm{extendLight, 1st}}, p_{A, \mathrm{extendLight, 2nd}}, P_{RR, \mathrm{extendLight, 2nd}} $
									視線サブパスの1, 2度目の延長に関する確率密度 $ p_{A, \mathrm{extendEye, 1st}}, P_{RR, \mathrm{extendEye, 1st}}, p_{A, \mathrm{extendEye, 2nd}}, P_{RR, \mathrm{extendEye, 2nd}} $ ) {
	MISウェイトの逆数を初期化 $ \mathrm{recMISWeight} \LAR 1 $
	\rh{30}
	光線サブパスの延長(視線サブパスの短縮)を考える
	if $ \mathrm{numEyeVtx} > 0 $
		視線サブパスの接続頂点を取得 $ \mathrm{eyeEndVtx} \LAR \mathrm{eyeSubPath.getNthVertex(} \mathrm{numEyeVtx} \mathrm{)} $
		1度目の延長における確率密度の比を計算 $ \displaystyle \mathrm{densityRatio} \LAR \frac{p_{A, \mathrm{extendLight, 1st}} P_{RR, \mathrm{extendLight, 1st}}}{\mathrm{eyeEndVtx.areaPDF} \times \mathrm{eyeEndVtx.RRProb}} $
		if 短縮した頂点の確率密度にデルタ関数が含まれていない
			$ \mathrm{recMISWeight} \LAR \mathrm{recMISWeight} + \mathrm{densityRatio}^2 $
		\rh{30}
		if $ \mathrm{numEyeVtx} - 1 > 0 $
			視線サブパスの接続頂点の1つ前の頂点を取得 $ \mathrm{eyeNextEndVtx} \LAR \mathrm{eyeSubPath.getNthVertex(} \mathrm{numEyeVtx} - 1 \mathrm{)} $
			2度目の延長における確率密度の比を計算 $ \displaystyle \mathrm{densityRatio} \LAR \mathrm{densityRatio} \times \frac{p_{A, \mathrm{extendLight, 2nd}} P_{RR, \mathrm{extendLight, 2nd}}}{\mathrm{eyeNextEndVtx.areaPDF} \times \mathrm{eyeNextEndVtx.RRProb}} $
			if 今回と前回、短縮した頂点のいずれの確率密度にもデルタ関数が含まれていない
				$ \mathrm{recMISWeight} \LAR \mathrm{recMISWeight} + \mathrm{densityRatio}^2 $
			\rh{30}
			for $ t \LAR \mathrm{numEyeVtx} - 2 $ to $ 1 $
				視線サブパスの短縮頂点を取得 $ \mathrm{newLightVtx} \LAR \mathrm{eyeSubPath.getNthVertex(} t \mathrm{)} $
				延長における確率密度の比を計算 $ \displaystyle \mathrm{densityRatio} \LAR \mathrm{densityRatio} \times \frac{\mathrm{newLightVtx.revAreaPDF} \times \mathrm{newLightVtx.revRRProb}}{\mathrm{newLightVtx.areaPDF} \times \mathrm{newLightVtx.RRProb}} $
				if 今回と前回、短縮した頂点のいずれの確率密度にもデルタ関数が含まれていない
					$ \mathrm{recMISWeight} \LAR \mathrm{recMISWeight} + \mathrm{densityRatio}^2 $
	\rh{40}
	視線サブパスの延長(光線サブパスの短縮)を考える
	if $ \mathrm{numLightVtx} > 0 $
		光線サブパスの接続頂点を取得 $ \mathrm{lightEndVtx} \LAR \mathrm{lightSubPath.getNthVertex(} \mathrm{numLightVtx} \mathrm{)} $
		1度目の延長における確率密度の比を計算 $ \displaystyle \mathrm{densityRatio} \LAR \frac{p_{A, \mathrm{extendEye, 1st}} P_{RR, \mathrm{extendEye, 1st}}}{\mathrm{lightEndVtx.areaPDF} \times \mathrm{lightEndVtx.RRProb}} $
		if 短縮した頂点の確率密度にデルタ関数が含まれていない
			$ \mathrm{recMISWeight} \LAR \mathrm{recMISWeight} + \mathrm{densityRatio}^2 $
		\rh{30}
		if $ \mathrm{numLightVtx} - 1 > 0 $
			光線サブパスの接続頂点の1つ前の頂点を取得 $ \mathrm{lightNextEndVtx} \LAR \mathrm{lightSubPath.getNthVertex(} \mathrm{numLightVtx} - 1 \mathrm{)} $
			2度目の延長における確率密度の比を計算 $ \displaystyle \mathrm{densityRatio} \LAR \mathrm{densityRatio} \times \frac{p_{A, \mathrm{extendEye, 2nd}} P_{RR, \mathrm{extendEye, 2nd}}}{\mathrm{lightNextEndVtx.areaPDF} \times \mathrm{lightNextEndVtx.RRProb}} $
			if 今回と前回、短縮した頂点のいずれの確率密度にもデルタ関数が含まれていない
				$ \mathrm{recMISWeight} \LAR \mathrm{recMISWeight} + \mathrm{densityRatio}^2 $
			\rh{30}
			for $ s \LAR \mathrm{numLightVtx} - 2 $ to $ 1 $
				光線サブパスの短縮頂点を取得 $ \mathrm{newEyeVtx} \LAR \mathrm{lightSubPath.getNthVertex(} s \mathrm{)} $
				延長における確率密度の比を計算 $ \displaystyle \mathrm{densityRatio} \LAR \mathrm{densityRatio} \times \frac{\mathrm{newEyeVtx.revAreaPDF} \times \mathrm{newEyeVtx.revRRProb}}{\mathrm{newEyeVtx.areaPDF} \times \mathrm{newEyeVtx.RRProb}} $
				if 今回と前回、短縮した頂点のいずれの確率密度にもデルタ関数が含まれていない
					$ \mathrm{recMISWeight} \LAR \mathrm{recMISWeight} + \mathrm{densityRatio}^2 $
	\rh{30}
	return $ 1 \; / \; \mathrm{recMISWeight} $
}
4行目 MISウェイトの逆数を $ 1 $ で初期化します。これは現在の接続に対応する戦略の確率密度に相当します。
6-23行目 光線サブパスの延長(視線サブパスの短縮)によって考えられる戦略の確率密度の比を加算していきます。基本的に視線サブパスの頂点を除去し終わるまで行いますが、$ t = 0 $ の戦略を考えない場合は、視線サブパスの1点目だけを残して終わります。
7-11行目 光線サブパスを1段階延長する際の確率密度の比を計算・加算を行います。
8行目 視線サブパスの除去される頂点 $ \mathrm{eyeEndVtx} $ を取得します。
9行目 延長に相当する項である分子は接続に依存した量である $ p_{A, \mathrm{extendLight, 1st}} $ と $ P_{RR, \mathrm{extendLight, 1st}} $ になります。一方で短縮に相当する項である分母は、接続には依存せず、除去される頂点に記録されている確率密度、確率からなります。
10-11行目 除去される頂点がデルタ関数にそってサンプリング(つまり完全スペキュラー面のサンプリング)されていた場合には、その確率密度に無限大を含むことになり、除去後に相当する戦略の確率密度は除去前に比べて無限に小さくなるので、その確率密度の(比の)蓄積を無視する必要があります。そうでない場合のみ、確率密度の比を蓄積します。パワーヒューリスティックを用いるため、厳密には比の2乗を加算します。
13-17行目 光線サブパスの2度目の延長に関する確率密度の比・加算を行います。
14行目 視線サブパスの除去される頂点 $ \mathrm{eyeNextEndVtx} $ を取得します。
15行目 延長に相当する項である分子はまだ接続に依存した量である $ p_{A, \mathrm{extendLight, 2nd}} $ と $ P_{RR, \mathrm{extendLight, 2nd}} $ になります。一方で短縮に相当する項である分母は、接続には依存せず、除去される頂点に記録されている確率密度、確率からなります。現在の接続に対する比を計算しなければならないので、計算済みの比に対してさらに乗算します。
16-17行目 確率密度の比を行うにあたっては、今回除去する頂点に加えて、前回除去した頂点、2つともがデルタ関数に無関係である必要があります。
19-23行目 光線サブパスの3度目以降の延長に関する確率密度の比・加算を、視線サブパスの頂点を除去し終わるまで行います。ここからは比の計算を全て頂点に記録されている情報だけで行えます。
25-42行目 視線サブパスの延長(光線サブパスの短縮)によって考えられる戦略の確率密度の比を加算していきます。手順としては光線サブパスの延長と全く同じなので説明は省略します。
44行目 ここまでで、MISウェイトの逆数を計算できたので、最後に逆数をとって返します。

$ s = 0 $ または $ t = 0 $ の戦略の寄与の蓄積

if 光線サブパストレース中
	レンズ面のサンプリングによって現在の点 $ \vx $ を得る確率密度を計算 $ p_{A, \mathrm{extendEye, 1st}} $	
	レンズ面における光の入射方向のサンプリングによって一つ前の点を得る確率密度を計算 $ \displaystyle p_{A, \mathrm{extendEye, 2nd}} \LAR \frac{p_\sigma(\vomega_\mathrm{traced}) \cdot \mathrm{cosLast}}{\mathrm{dist2}} $
	$ \begin{aligned}[t] w_{s, 0} \LAR \mathrm{calculateMISWeight} (& \mathrm{subPath}, \mathrm{null}, \mathrm{subPath.numVertices}, 0, \\& 0, 0, 0, 0, \\& p_{A, \mathrm{extendEye, 1st}}, 1, p_{A, \mathrm{extendEye, 2nd}}, 1 ) \end{aligned} $
	寄与を計算 $ C_{s, 0} \LAR w_{s, 0} \times \alpha \times W_e(\vx, \vomega_\mathrm{traced}) $
	現在のレイからピクセル位置 $ j^* $ を逆算して寄与を蓄積する $ I_{light, j^*} \LAR I_{light, j^*} + C_{s, 0} $
else
	ヒットしている光源が選ばれる確率を計算 $ P_\mathrm{light} $
	光源面のサンプリングによって現在の点 $ \vx $ を得る確率密度を計算 $ p_{A, \mathrm{extendLight, 1st}} $	
	光源が選ばれる確率を乗算 $ p_{A, \mathrm{extendLight, 1st}} \LAR P_\mathrm{light} p_{A, \mathrm{extendLight, 1st}} $
	光源面における光の出射方向のサンプリングによって一つ前の点を得る確率密度を計算 $ \displaystyle p_{A, \mathrm{extendLight, 2nd}} \LAR \frac{p_\sigma(\vomega_\mathrm{traced}) \cdot \mathrm{cosLast}}{\mathrm{dist2}} $
	$ \begin{aligned}[t] w_{0, t} \LAR \mathrm{calculateMISWeight} (& \mathrm{null}, \mathrm{subPath}, 0, \mathrm{subPath.numVertices}, \\& p_{A, \mathrm{extendLight, 1st}}, 1, p_{A, \mathrm{extendLight, 2nd}}, 1, \\& 0, 0, 0, 0 ) \end{aligned} $
	寄与を計算 $ C_{0, t} \LAR w_{0, t} \times L_e(\vx, \vomega_\mathrm{traced}) \times \alpha $
	処理中のピクセル位置 $ j $ に寄与を蓄積する $ I_{eye, j} \LAR I_{eye, j} + C_{0, t} $

サブパストレース中に完全な光輸送経路が完成した場合は、寄与の計算にもう一方のサブパスの情報を必要としません。

1-6行目 光線サブパストレース中にレンズ(センサー)面にヒットした場合の処理です。現在処理中のピクセルとは無関係なピクセル位置にレイがヒットすることになるため、ピクセル位置の逆算が必要となります。
2行目 現在の点 $ \vx $ が視線サブパスのトレース、つまりレンズ面のサンプリングによって生成される確率密度を計算します。
3行目 一つ前の点がレンズ面から光の入射方向をサンプリングすることによって得られる確率密度を計算します。面積測度に変換する必要がありますが、一つ前の点におけるcos項と距離の2乗から計算できます。
4行目 MISウェイトを計算します。最初の2つの頂点生成にはロシアンルーレットを使っていないので対応する引数には確率 $ 1 $ を渡しています。
5-6行目 MISウェイトとセンサー応答、そして現在のサブパスのウェイトを用いて寄与を計算します。そして、現在のレイからセンサー上にあるピクセル位置を逆算して寄与を蓄積します。
7-14行目 視線サブパストレース中に光源面にヒットした場合の処理です。基本的には光線サブパスの場合と同じ処理ですが、ヒットした光源自体が選ばれる確率を考慮する必要があります。一方、寄与の蓄積先は単純に現在処理中のピクセル位置になります。

コードまとめ

light imageとeye imageを初期化 $ I_{light, j} \LAR 0,\; I_{eye, j} \LAR 0 $ for 全ての画素 $ j $
for $ 0 $ to サンプル数 $ N - 1 $ {
	for 画素 $ j $ in 全てのピクセル {
		光線サブパス $ \mathrm{lightSubPath} $, 視線サブパス $ \mathrm{eyeSubPath} $ を初期化
		\rh{30}
		// アルゴリズム1 .光線サブパスの構築
		...
		\rh{30}
		// アルゴリズム2. 視線サブパスの構築
		...
		\rh{30}
		// アルゴリズム4. サブパスの接続
		...
	}
}
light imageとeye imageを加算 $ I_j \LAR I_{light, j} + I_{eye, j} $ for 全ての画素 $ j $
画素値をサンプル数で割る $ I_j \LAR I_j \; / \; N $ for 全ての画素 $ j $

ここまでに説明した各コード片を組み合わせて双方向パストレーシングのアルゴリズムとなります。プログラムの最も外枠としては、(light imageがありますが)基本的には各ピクセル独立に処理して最後にサンプル数で平均をとって結果出力、というパストレーシングと同じ処理になります。

1行目 ほとんどの戦略の寄与の蓄積先となるeye imageに加えて $ t \le 1 $ の戦略用にlight imageの初期化を行います。
16行目 eye imageとlight imageの画像を合成します。
17行目 サンプル数 $ N $ で結果を割って、結果画像を求めます。

デバッグ

双方向パストレーシングはここまでに示した通り、(疑似コードだとそんなに大したものに見えないかもしれないですが)結構複雑な実装になります。そのため、いきなりはアルゴリズム全体としての検証が難しいかもしれません。そこで双方向パストレーシングのデバッグとして次に示すような、部分的な検証を重ねていくことをお勧めします。

リファレンスレンダラーの確立

双方向パストレーシングに限らず、新たな光輸送アルゴリズムを実装した際には正当性の検証として、既存のアルゴリズムの結果(十分時間をかけてレンダリングして分散を抑えたもの)と一致するかを見ることがよく行われます。この検証を意味あるものとする上で、既存のアルゴリズムの実装が本当に正しいのか?ということに注意する必要があります。正当さの証明は難しい問題ですが、いくつか試せることがあります。
双方向パストレーシングに対するリファレンスレンダラーとしてはパストレーシングが使われることが多いと思われます。しかし、収束した画像を作る上で純粋なパストレーシングだと非常に時間がかかってしまうため、Next Event Estimationを使ったりすると思いますが、その実装が純粋なパストレーシングと一致するかどうかは予め十分に検証しておきましょう。

ライトトレーシングの検証

レンズ面のサンプリングやBSDFのサンプリングは、パストレーシング、双方向パストレーシングで共通して使われる一方で、光源からのレイの生成や、$ t \le 1 $ の戦略におけるセンサー上の位置の逆算などパストレーシングには無かった機能の実装が必要になります。これら新機能の実装とMISウェイト計算の検証を同時に行うとバグの特定が難しくなるので、まずは前者の正当性の検証を行うべきです。$ t \le 1 $ の戦略は、実はライトトレーシングという名前のアルゴリズムになります。まずはライトトレーシングが正しく実装できるか、既存アルゴリズムと一致した結果が得られるか、が重要になります。注意点として、ライトトレーシングにおけるNext Event Estimation (パストレーシングとは逆に、レンズ上の1点を明示的にサンプリングして接続)では、完全スペキュラー面をレンダリングできない、などアルゴリズムの特性上比較が不可能な部分もあります。

特定の経路長の検証

双方向パストレーシングではある経路長に関して様々な戦略を考えますが、これらの戦略は得意不得意はあっても(サンプリング不可能な場合を除いて)同じ結果に収束するはずです。特定の戦略・特定の経路長のみで寄与を蓄積させるようにして、各戦略から生成される結果が一致するかを検証します。注意点としてはこの際にはMISウェイトをかけずに結果を出力するようにします。逆に、MISウェイトをかけた結果とこれらを比較すれば、MISの効果を実感することができるでしょう。

拡散反射面のみのシーン

スペキュラー面が絡んだMISウェイトの計算はコード実装上は特殊処理を行う必要があり、間違った実装になりやすいため、まずは拡散反射面のみで正当性の検証を行った上で、スペキュラー面のテストを行うことをお勧めします。

法線補間は使わない

法線補間を使うことで、少ないポリゴン数でも滑らかなライティングに見せることができますが、パストレーシングと違って双方向パストレーシングでは光源からのトレースも行う関係上、パストレーシングと一致した結果を得るためには少々特殊な考慮が必要となります。この問題については「双方向アルゴリズムにおける法線補間の扱い方」で紹介したいと思いますが、まず、アルゴリズムの本質的な正しさの検証を行う上では法線補間を使わないシーンを使うべきでしょう。

さいごに

パストレーシングに比べると双方向からトレースを行う関係もあって疑似コードが随分と巨大になってしまいましたが、逆に言えば双方向のところはほぼ同じことを行なっているため、その分を省いて考えれば少しは気が楽になります。理論をしっかりと理解した上で読んでいただければ個々の部分はそれほど難しいことはしていないとわかると思います。

参考文献


Today : 00000071 Total : 00061010
Copyright © 2017 shocker.