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

パストレーシング実装

このページはレイトレ合宿3!!!のアドベントカレンダー3週目の記事として制作しました。
みなさんこんにちは!shockerです。このウェブページではMCRTまわりの理論などについて紹介してきましたが、今まで肝心のレンダリング手法・特に実装面に関しては何も書いていませんでした。理論面は別のページで紹介済みですし、良い機会なのでMCRTを始めるにはまずここから!といった感じのパストレーシングの実装面について解説します。
と言いつつ私が初めて実装したのはフォトンマッピングでしたが、最初はパストレがおすすめだと思います...

パストレーシング」ではアルゴリズムの理論面を紹介しましたが、実際にプログラムで書くとなるといくつか考える点が出てきます。このページではまず純粋なパストレーシングの実装について書いた後、実用的な実装について紹介します。
パストレーシングは初歩的な手法ですが、発展手法と基本的には全く同じ絵を描き出すことができます。異なるのは真値への収束スピードだけです。そして初歩的、言い換えれば素直な手法なのでレンダリングにおける様々なフィーチャーを実装しやすく、映画などで用いられるプロダクションレンダラーでも数多く採用されています。侮るなかれ。

純粋なパストレーシング

画像を初期化 $ I_j \LAR 0 $ for 全ての画素 $ j $
for $ 0 $ to サンプル数 $ N - 1 $ {
	for 画素 $ j $ in 全てのピクセル {
		画素 $ j $ 中の1点をサンプル $ \vx_p $、面積に関する確率密度 $ p_A(\vx_p) $
		レンズ中の1点をサンプル $ \vx_0 $、面積に関する確率密度 $ p_A(\vx_0) $
		$ \vx_p $, $ \vx_0 $からトレースするレイを求める $ \mathrm{ray} \LAR \{ \mathrm{org:\;} \vx_0,\; \mathrm{dir:\;} \vomega_{01} \} $、方向に関する確率密度も求める $ p_{\sigma}(\vomega_{01}) $
		経路のウェイトを初期化 $ \displaystyle \alpha \LAR \frac{W_e^j(\vx_0, \vomega_{01}) \abs{\vomega_{01} \! \cdot \vec{n}_0}}{p_A(\vx_0) p_\sigma (\vomega_{01})} $
		\rh{30}
		while レイ $ \mathrm{ray} $ をトレースして衝突がある {
			衝突点の情報を取得: 位置 $ \vx $, 法線 $ \vec{n} $, 光の出射方向 $ -\mathrm{ray.dir} = \vomega_o $, BSDF $ f_s(\cdot) $
			\rh{30}
			if 衝突点 $ \vx $ が発光している
				$ I_j \LAR I_j + \alpha \times L_e(\vx, \vomega_o)$
			\rh{30}
			光の入射方向をサンプル $ \vomega_i $、方向に関する確率密度 $ p_\sigma(\vomega_i) $
			ウェイトを更新 $ \displaystyle \alpha \LAR \alpha \times \frac{f_s(\vx, \vomega_i, \vomega_o) \cdot \abs{\vomega_i \! \cdot \vec{n}}}{p_\sigma(\vomega_i)} $
			\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_i \} $
			ウェイトにロシアンルーレットの確率を反映 $ \displaystyle \alpha \LAR \alpha \times \frac{1}{P_{RR}} $
		}
	}
}
画素値をサンプル数で割る $ I_j \LAR I_j \; / \; N $ for 全ての画素 $ j $

まずはNext Event Estimationも用いない純粋なパストレーシングについて紹介します。正直に言うとこの実装は光源がそれなりに大きい場合でも分散が非常に大きくなり実用的では無いのですが、レンダリング方程式とモンテカルロ法の結びつきを考える上では最も直感的な形となるので重要です。

それでは順を追って説明していきます。

1行目 画像をゼロで初期化しています。RGB画像なら各画素には3成分あるでしょう。ここは問題ないですね。
2行目 サンプル数分のループが記述されており、このサンプル数 $ N $ が多い程生成される画像のクオリティが高くなります。が、当然実行時間も長くなります。
3行目 各ピクセルに関するループです。特に言うことはないでしょう。
4行目 現在のピクセル中のある1点をサンプリングしています。理論のページで説明した通り、パストレーシングでは目(つまりセンサー側)から経路を構築していきます。確率密度 $ p_A(\vx_p) $ はピクセル中で一様サンプリングを行うのならば単純に 1 / (ピクセルの面積)になります。
5行目 レンズ中の1点をサンプリングしています。この確率密度も一様サンプリングを行うのならば 1 / (レンズの面積)になります。ピンホールカメラモデルの場合は $ \vx_0 $ は固定です。この場合のPDFはディラックのデルタ関数になってしまい、確率密度が無限の値になりますが、プログラム的にはとりあえず1.0にしておいたりします。
6行目 ピクセル中の位置 $ \vx_p $ とレンズ上の位置 $ \vx_0 $ が決まれば、レンズはスペキュラー透過なので屈折後の方向が一意に定まります。レンズ上の点を原点として屈折後の方向を持ったレイを生成します。この方向に関する確率密度 $ p_\sigma(\vomega_{01}) $を計算します。この辺りは「被写界深度」のページで解説しています。
7行目 経路のウェイトを初期化します。理論のページで説明した通り、分子にはセンサー応答関数などが含まれています。センサー応答関数の決め方は自由です。何か特定のレンズ・センサーの組み合わせを模倣するのでも無い限り、分子は定数でも良いでしょう。 ところで、理想上のピンホールカメラは穴の大きさが無限小なので通過する光もそれに比例して無限小です。それでは真っ暗な画像しか生成できないので、センサー応答が無限大だと考えます。どんな無限大でも良いわけではなく、「5行目のレンズ面のPDFになるデルタ関数の無限大」に比例した無限大です。面倒なことを言いましたが、結局同程度の無限大であれば、ウェイトの分母分子でキャンセルし合うことで定数になります。
9行目 ここからはロシアンルーレットで打ち切られない限り、経路を延長しながら各経路長における寄与を蓄積していきます。
10行目 衝突点が見つかった場合、その点を $ \vx $ とします。この点における光の出射方向 $ \vomega_o $ は、パストレーシングでは逆方向のトレースを行っているためレイとは逆方向になります。法線 $ \vec{n} $やBSDF $ f_s(\cdot) $ の情報も取得しておきます。
12 - 13行目 衝突した物体が発光物体だった場合、方向 $ \vomega_o $ に出ている発光放射輝度 $ L_e $を取得、ウェイトをかけて画素値に蓄積します。この $ L_e $ は適当な実装なら適当な定数にしていますが、厳密に定義したいという方は「光源の出力と放射輝度の関係」をご参照ください。
15行目 より長い経路の寄与を考えるために、2~3個の乱数を使って光の入射方向 $ \vomega_i $ をサンプリングします。その方向をサンプリングする確率密度 $ p_\sigma(\vomega_i) $も計算しておきます。全球や半球からの一様サンプリングでも良いですが、理論のページで紹介したようにBSDFやBSDFとcos項の積に対する重点的サンプリングを行うのが良いでしょう。完全スペキュラー面は適当なサンプリングでは計算できません(詳細後述)。
16行目 (位置 $ \vx $ と)入射・出射方向 $ \vomega_i $, $ \vomega_o $が決まればBSDFの値が求まります。 完全スペキュラー面においてBSDFの重点的サンプリングを行っている場合はこの値がデルタ関数を含むことになりますが、プログラム上ではこれもとりあえず定数で表現します(詳細後述)。 求めたBSDFの値や確率密度を使って経路のウェイトを更新しておきます。
18行目 無限にトレースし続けていると処理が終わらないので、理論のページで紹介したロシアンルーレットを導入します。トレースを継続する確率 $ P_{RR} $ には適当な定数を用いても良いですが、より良い選択肢として、ウェイトの更新時の係数だったり、ウェイトの初期値に対する現在のウェイトの比などがあります。RGBレンダリングの場合、ウェイトは3要素の値を持っているので、3つのうちの最大値を用いたりします。
19-20行目 一様乱数を生成して $ P_{RR} $ 以上の値になった場合はループを抜けて、このピクセル、このサンプルに関する処理を終了します。
22行目 新たなレイを生成します。原点は現在の衝突点、方向はサンプリングした入射方向になります。この生成方法だとプログラム的に問題があるので解決方法は後述します。
23行目 トレース処理を続ける場合は、まずロシアンルーレットの確率に応じてウェイトを修正する必要があります。この修正によって期待値を真値に一致させたまま処理を有限にできます。
27行目 各ピクセル、サンプル数分蓄積された寄与をその数で割ることで平均値を算出します。

BSDFに対する重点的サンプリング

完全拡散反射面

理論のページでは、完全拡散面の重点的サンプリングにはBRDFとcos項の積に比例したサンプリング $ p_\sigma (\vomega) = \cos{\theta} \; / \; \pi $ がよく使用されると書きました。そのサンプリングはどのようにして行うのでしょうか。1つの答えとして逆関数法(inverse transform sampling, inversion sampling, etc)を用います。

まずはサンプル分布を比例させたいPDFをそのパラメターに関して不定積分してCDFの式を求めます。

\begin{eqnarray*} P(\theta, \phi) = \int \int \frac{\cos{\theta}}{\pi} \sin{\theta} d\theta d\phi &=& \frac{1}{2\pi} \int \int \sin{2\theta} d\theta d\phi \\ &=& \frac{1}{2\pi} \int \Brcc{A - \frac{1}{2} \cos{2\theta}} d\phi \\ &=& \frac{\phi}{2\pi} \Brcc{A - \frac{1}{2} \cos{2\theta}} + B \\ \end{eqnarray*}

境界条件 $ P(0, 0) = 0 $ より $ B = 0 $ が求まり、もう一つの境界条件 $ P(\pi \; / \; 2, 2\pi) = 1 $ より $ A = 1 \; / \; 2 $ が求まります。したがってCDFは以下の式になります。

\begin{equation*} P(\theta, \phi) = \frac{\phi}{2\pi} \Brcc{\frac{1}{2} - \frac{1}{2} \cos{2\theta}} = \frac{\phi}{2\pi} \sin^2 \theta \end{equation*}

$ \theta $ と $ \phi $ を独立に考えられるので、それぞれのCDFは $ P(\theta) = \sin^2 \theta $, $ P(\phi) = \phi \; / \; 2\pi $ となります。逆関数法では求めたCDFの逆関数に一様乱数を入力することで、所望の分布を得ることができます。つまり、$ \theta $, $ \phi $ は次のように求められます。

\begin{equation*} \phi = 2\pi u_1 \hspace{5mm} \theta = \sin^{-1} \hspace{-4mm} \sqrt{u_2} \end{equation*}

ここで、$ u_1, u_2 $ は $ [0, 1) $ に分布する一様乱数です。このようにして求めた入射方向のPDFは当然ですが $ p_\sigma (\vomega) = \cos{\theta} \; / \; \pi $ になります。
注意点として、算出した方向は法線を鉛直上向きとする座標系内での方向になるので、着目している物体表面の向きに合わせて座標変換する必要があります。

完全スペキュラー面

完全なスペキュラー反射、完全鏡面反射のBSDFを式で表すと次のようになります。

\begin{equation*} f_s(\vx, \vomega_i, \vomega_o) = \rho \frac{\delta(\vomega_i - \mathrm{R}_{\vec{n}}(\vomega_o))}{\abs{\cos \theta_i}} \end{equation*}

ここで $ \rho $ は反射率、$ \theta_i $ は入射方向の天頂角です。$ \mathrm{R}_{\vec{n}}(\vomega) $ は点 $ \vx $ における法線に関して与えられた方向の鏡面反射方向を返す関数です。デルタ関数によって表されているので、BSDFが値を持つのは入射角・反射角が鏡面反射方向の関係に「完全に」一致した場合のみです。少しでもずれると一切の値を持たないので、適当なサンプリングによって完全スペキュラー面のBSDFが値を持つ確率密度もゼロです。ですのでBSDFにそった重点的サンプリングを用いる必要があります。というより入射方向は一意に決定されます。この一意に決定される方向のPDFは次のようにデルタ関数で表されます。

\begin{equation*} p_\sigma(\vx, \vomega_i, \vomega_o) = \delta(\vomega_i - \mathrm{R}_{\vec{n}}(\vomega_o)) \end{equation*}

スペキュラー面のBSDF・PDFが値を持つ場合には無限大の数値となるため、プログラム内ではとりあえず定数で持ちますが、BSDFとPDFそれぞれに「同じ」デルタ関数が含まれることになるので、キャンセルされて面の反射率とcosのみが残ります。ちなみに反射率は簡単な実装なら適当な定数で良いですが、よりリアルな表現を目指すのならフレネル反射を考慮することになります。

自己交叉の回避

アルゴリズム1の18行目で現在の衝突点と、サンプリングした光の入射方向から次にトレースするレイを生成していますが、プログラムでこれをそのまま実装すると浮動小数点の誤差により、新しいレイをトレースしても現在の衝突点とほぼ同じ点に衝突があったと判定されてしまうことがあります。この問題を回避するためによく使用される方法として、以下のようにレイの原点に法線方向のオフセットを加えることがあります。

\begin{equation*} \mathrm{ray} \LAR \{ \mathrm{org:\;} \vx + \epsilon \vec{n},\; \mathrm{dir:\;} \vomega_i \} \end{equation*}

ここで $ \epsilon $ はごくごく小さな値で、典型的には0.0001とかそんな値です。注意点として、透過方向に光の入射方向がサンプリングされていた場合は法線を逆向きに考える必要があります。わりとアドホックな手法に見えますが大体はこれで解決できます。
他に考えられる手法としては、レイが衝突した物体のIDを記録しておき、次のトレース時にはそれを無視するといったことも考えられます。ただしこれは三角形メッシュではない解析的な曲面などをレイトレースする場合には無効になります。

Next Event Estimation

理論のページでも紹介した通り、純粋なパストレーシングの実装ではランダムにトレースされた経路が光源に当たらない限り寄与が蓄積できず、光源が小さくなればなるほど結果の分散が大きくなってしまいます。Next Event Estimationにより結果の改善が期待できます。

光源上の1点をサンプル $ \vx_l $、面積に関する確率密度 $ p_A(\vx_l) $
シャドウレイ $ \mathrm{shadowRay} \LAR \{ \mathrm{org:\;} \vx,\; \mathrm{dir:\;} \vomega_l \LAR \mathrm{normalize}(\vx_l - \vx) \} $
if レイ $ \mathrm{shadowRay} $ をトレースして $ \vx_l $ との間に遮蔽物が無い場合
	$ \displaystyle I_j \LAR I_j + \alpha \times L_e(\vx_l \RAR \vx) \frac{f_s(\vx, \vomega_l, \vomega_o) G(\vx_l \LRAR \vx)}{p_A(\vx_l)} $
1行目 Next Event Estimationでは光源上の一点 $ \vx_l $ を明示的にサンプリングします。光源面中で一様サンプリングを行う場合は、確率密度 $ p_A(\vx_l) $ は単純に 1 / (光源面の面積)になります。
2行目 現在の点から光源上のサンプル点へと向かうシャドウレイ $ \mathrm{shadowRay} $ を生成します。
3 - 4行目 シャドウレイをトレースして間に遮蔽物が無ければ、寄与を計算して画素値に蓄積します。

スペキュラー面の取り扱い

上でも書きましたが、完全鏡面のBSDFが一意的な方向選択以外で値を持つことは無いので、Next Event Estimationでは完全鏡面を取り扱うことができません。なので完全鏡面に関しては純粋なパストレーシングと同じく、経路延長時に光源にヒットできれば寄与の蓄積を行います。

コードまとめ

画像を初期化 $ I_j \LAR 0 $ for 全ての画素 $ j $
for $ 0 $ to サンプル数 $ N - 1 $ {
	for 画素 $ j $ in 全てのピクセル {
		画素 $ j $ 中の1点をサンプル $ \vx_p $、面積に関する確率密度 $ p_A(\vx_p) $
		レンズ中の1点をサンプル $ \vx_0 $、面積に関する確率密度 $ p_A(\vx_0) $
		$ \vx_p $, $ \vx_0 $からトレースするレイを求める $ \mathrm{ray} \LAR \{ \mathrm{org:\;} \vx_0,\; \mathrm{dir:\;} \vomega_{01} \} $、方向に関する確率密度も求める $ p_{\sigma}(\vomega_{01}) $
		経路のウェイトを初期化 $ \displaystyle \alpha \LAR \frac{W_e^j(\vx_0, \vomega_{01}) \abs{\vomega_{01} \! \cdot \vec{n}_0}}{p_A(\vx_0) p_\sigma (\vomega_{01})} $
		\rh{30}
		while レイ $ \mathrm{ray} $ をトレースして衝突がある {
			衝突点の情報を取得: 位置 $ \vx $, 法線 $ \vec{n} $, 光の出射方向 $ -\mathrm{ray.dir} = \vomega_o $, BSDF $ f_s(\cdot) $
			\rh{30}
			if 完全スペキュラー面をトレースしてきた かつ 衝突点 $ \vx $ が発光している
				$ I_j \LAR I_j + \alpha \times L_e(\vx, \vomega_o)$
			if BSDFが完全スペキュラー以外の成分を持っている { 
				光源上の1点をサンプル $ \vx_l $、面積に関する確率密度 $ p_A(\vx_l) $
				シャドウレイ $ \mathrm{shadowRay} \LAR \{ \mathrm{org:\;} \vx,\; \mathrm{dir:\;} \vomega_l \LAR \mathrm{normalize}(\vx_l - \vx) \} $
				if レイ $ \mathrm{shadowRay} $ をトレースして $ \vx_l $ との間に遮蔽物が無い場合
					$ \displaystyle I_j \LAR I_j + \alpha \times L_e(\vx_l \RAR \vx) \frac{f_s(\vx, \vomega_l, \vomega_o) G(\vx_l \LRAR \vx)}{p_A(\vx_l)} $
			}
			\rh{30}
			光の入射方向をサンプル $ \vomega_i $、方向に関する確率密度 $ p_\sigma(\vomega_i) $
			ウェイトを更新 $ \displaystyle \alpha \LAR \alpha \times \frac{f_s(\vx, \vomega_i, \vomega_o) \cdot \abs{\vomega_i \! \cdot \vec{n}}}{p_\sigma(\vomega_i)} $
			\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_i \} $
			ウェイトにロシアンルーレットの確率を反映 $ \displaystyle \alpha \LAR \alpha \times \frac{1}{P_\mathrm{RR}} $
		}
	}
}
画素値をサンプル数で割る $ I_j \LAR I_j \; / \; N $ for 全ての画素 $ j $

さいごに

Next Event Estimationによってかなり実用的なアルゴリズムになりますが、まだまだ改善の余地はあります。Spherical Triangle (Rectangle)の概念によってNext Event Estimation自体を改善する手法[Arvo1995, Ureña2013]もありますし、Multiple Importance Samplingを用いてBSDFにそった重点的サンプリングと組み合わせる方法もあります。

参考文献

  • [Pharr2010] Matt Pharr, Greg Humphreys - "PHYSICALLY BASED RENDERING: From Theory to Implementation Second Edition", 2010, Morgan Kaufmann
  • [Arvo1995] James Arvo - "Stratified Sampling of Spherical Triangles", 1995
  • [Ureña2013] Carlos Ureña - "An Area-Preserving Parametrization for Spherical Rectangles", 2013

Today : 00000095 Total : 00054709
Copyright © 2017 shocker.