パストレーシング[Kajiya1986]はモンテカルロ法を用いてレンダリング方程式の解を求めるレンダリング手法のひとつであり、様々な発展手法のベースとも考えられる重要な手法です。 パストレーシングでは現実の光の向きとは逆に、視点から光輸送経路を考えます。
図1. パストレーシングによるレンダリング結果 |
現実では、光源から出た光がカメラのセンサーや人間の眼に入射することで、画像が生成されたり知覚されたりします。ただ、光源から出た光のうちセンサーに入射するのはごく一部であるため、パストレーシングでは光の輸送経路を古典的レイトレーシングと同様にセンサー側から逆追跡します。
「レンダリング方程式」で紹介したレンダリング方程式を見ればわかるように、反射・透過する光の量を計算するためにはあらゆる入射方向に関して積分を解く必要がありますが、一般的なシーン設定では解析的に解くのが困難であるためモンテカルロ法を用いて、入射方向を確率的にサンプリングします。また、入射方向のサンプリング数も基本的には1つに絞ることで、再帰的な反射・透過を扱いながらも計算量の爆発を抑えます。代償としてレンダリング画像には推定結果の分散が視覚的なノイズとして現れることになりますが、パストレーシングを何回も行い平均をとることで分散を低減させていきます。
このようにパストレーシングは、レイトレーシングをベースとしてモンテカルロ法を利用したレンダリングを行うため、モンテカルロ・レイトレーシング(Monte Carlo Ray Tracing, MCRT)と呼ばれる手法に属し、大域照明を解くMCRTの中では最も基本的な手法と考えられます。大域照明を計算するMCRTには様々な手法が存在します。それぞれの手法に得意不得意がありますが基本的にはどの手法もモンテカルロ積分のもと、様々な効果━━ソフトシャドウや光沢反射、被写界深度やモーションブラー、間接照明、コースティクスなど━━を描き出すことができます。さらに「ボリュームレンダリング方程式」まで考慮すればもっと多くの効果━━煙や雲、表面下散乱、炎など━━を描くことができます。このページでは物体表面間だけの光輸送を考える基本的なパストレーシングについて紹介します。
図1にパストレーシングによるレンダリング画像を載せています。直方体の側面などが間接照明によって色づいていることが確認できます。
センサーのあるピクセル上の点 $ \vx_0 $ に、ある方向 $ \vomega $ から入射する放射輝度 $ L_i(\vx_0, \vomega) $ を考えます。真空のシーンを考える場合、放射輝度は進む距離に対して不変であるため、$ L_i $ は $ \vx_0 $ から $ \vomega $ 方向にレイ(光線)をトレースして最初に見つかる点 $ \vx_1 $ を $ \vomega $ 方向に出る放射輝度 $ L_o(\vx_1, \vomega) $ に等しくなります(入射に関しては方向の表記が逆になります, 図2)。 \begin{equation*} L_i(\vx_0, \vomega) = L_o(\vx_1, \vomega) \end{equation*}
「レンダリング方程式1」で説明したように、$ L_o $ は $ \vx_1 $ が光源として出す光と、周囲から入ってきて $ \vomega $ 方向に反射した光の和で表されます。従って、$ \vx_1 $ が $ \vomega $ 方向に光を自ら発している場合はその寄与がまず蓄積されます。次に反射光成分を評価する必要があるのですが、これは $ \vx_1 $ から見たあらゆる方向に関する積分です。この項に対してモンテカルロ法を適用し、あるひとつの光の入射方向 $ \vomega_{i, 1} $ を確率的に選択します。例えば周囲の方向から一様分布でサンプリングする場合は、全立体角が $ 4\pi $ なので確率密度は $ 1 \;/\; 4\pi $ となり、$ L_o $ の推定値は次のようになります。 \begin{equation} \Abka{L_o(\vx_1, \vomega)} = L_e(\vx_1, \vomega) + \frac{ f_s(\vx_1, \vomega_{i, 1}, \vomega) L_i(\vx_1, \vomega_{i, 1}) \absb{\vomega_{i, 1} \! \cdot \vec{n}_1} }{\frac{1}{4\pi}} \label{eq:LTE_sphere_sampling} \end{equation} 確率的に選択された方向 $ \vomega_{i, 1} $ から入射する放射輝度を求める必要がありますが、$ L_i(\vx_1, \vomega_{i, 1}) $ は直接評価することができません。しかし、この $ L_i $ も同様に $ \vx_1 $ から $ \vomega_{i, 1} $ 方向にレイをトレースして最初に見つかる点から出射される放射輝度に等しくなるため、式 $ \eqref{eq:LTE_sphere_sampling} $ は次のように書き換えることができます。 \begin{eqnarray*} \Abka{L_o(\vx_1, \vomega)} &=& L_e(\vx_1, \vomega) + \frac{f_s(\vx_1, \vomega_{i, 1}, \vomega) \absb{\vomega_{i, 1} \! \cdot \vec{n}_1}}{\frac{1}{4\pi}} \Brc{L_e(\vx_2, \vomega_{i, 1}) + \frac{f_s(\vx_2, \vomega_{i, 2}, \vomega_{i, 1}) L_i(\vx_2, \vomega_{i, 2}) \absb{\vomega_{i, 2} \! \cdot \vec{n}_1}}{\frac{1}{4\pi}}} \\ &=& L_e(\vx_1, \vomega) + \frac{f_s(\vx_1, \vomega_{i, 1}, \vomega) \absb{\vomega_{i, 1} \! \cdot \vec{n}_1}}{\frac{1}{4\pi}} L_e(\vx_2, \vomega_{i, 1}) + \\ && \hspace{50mm} \frac{f_s(\vx_1, \vomega_{i, 1}, \vomega) \absb{\vomega_{i, 1} \! \cdot \vec{n}_1}}{\frac{1}{4\pi}} \frac{f_s(\vx_2, \vomega_{i, 2}, \vomega_{i, 1}) L_i(\vx_2, \vomega_{i, 2}) \absb{\vomega_{i, 2} \! \cdot \vec{n}_1}}{\frac{1}{4\pi}} \end{eqnarray*} 続く項も同様に再帰的に展開を行うことで、$ L_o(\vx_1, \vomega) $ の推定値を得ることができます。この式は一項目が0回反射、二項目が1回反射、...、$ n $ 項目が $ n - 1 $ 回反射を表しています。「レンダリング方程式2」で説明したノイマン級数展開や経路積分の話と同じですね。
さて、ここまででセンサー上の点 $ \vx_0 $ に方向 $ \vomega $ から届く放射輝度の計算方法を示しました。センサーは面積を持っていますし、様々な方向からの光に応答すると考えられるので $ L_i(\vx_0, \vomega) = L_o(\vx_1, \vomega) $ をそれらに関して積分することで、画素の値を求めることができます。この部分にも同様にモンテカルロ法を適用することで、$ j $ 番目のピクセル値 $ I_j $ は次のように推定されます。 \begin{equation*} \Abk{I_j} = \frac{W_e^j(\vx_0, \vomega) L_i(\vx_0, \vomega) \absb{\vomega \!\! \cdot \vec{n}_0}}{p_A(\vx_0) p_\sigma (\vomega)} \end{equation*} ここで、$ W_e^j(\vx_0, \vomega) $ は点 $ \vx_0 $、方向 $ \vomega $ に対するセンサー応答、$ p_A(\vx_0) $ はピクセル $ j $ 中の位置 $ \vx_0 $ をサンプルする面積測度に関する確率密度、$ p_\sigma (\vomega) $ は $ \vx_0 $ から方向 $ \vomega $ をサンプルする立体角測度に関する確率密度です。このピクセル値はスカラーで色を考慮していないため実際には例えばRGBの3つの値に関して評価することになります。その場合は輝度 $ L $、BSDF $ f_s $、センサー応答 $ W_e $ がRGBの値を持ちます。一方でRGBの値を全て同時に計算したい場合は同じ経路を共有することになるので確率密度やcos項といった他の要素はスカラー値のままです。
カメラを考える場合には、シーンとセンサーの間にはレンズ、ピンホールカメラの場合には光が通る小さな穴があります(図3)。そしてセンサーに到達する光は全てレンズや穴を通ってくるものと考えられます。レンズは理想的に完全鏡面透過ですし、ピンホールも理想的には大きさ0の穴です。すなわち、レンズ上の点が決定される(ピンホールの場合は必ず同じ一点)と、センサー上の点 $ \vx_0 $ はシーン中の1点に一意に対応します。そのような理由から、多くのレンダリングシステムでは、レンズをシーン中の物体と考えず、レンズ上の点を光輸送経路の端点 $ \vx_0 $ と捉えています。また、センサー応答もレンズ上の位置と入射方向に関して設定されます。つまり、センサーとレンズの間の光輸送に関しては隠蔽されています。レンズ上の点 $ \vx_0 $ が決まった後、ピクセル上の位置 $ \vx_p $ が決まれば $ \vx_0 $ からトレースされるレイが定まり、結果シーン中の最初の衝突点 $ \vx_1 $ も定まります。この最初の衝突点 $ \vx_1 $ を生成するPDFの求め方などについては「被写界深度」にて説明しています。
前節で説明した再帰的なレイのトレースと $ L_e $ の蓄積により、ピクセルに入射する放射輝度 $ L_o(\vx_1, \vomega) $ の推定が行えます。
ただ、このままでは無限に再帰的な展開が行われることになってしまいますので、どこかで処理を打ち切る必要があります。そこで、「モンテカルロ法その他」の「無限和への適用」で説明した手法を用います。
各点で光の入射方向を確率的に選んでトレースを行っていましたが、このトレース自体を確率的に継続するようにします。1回目の反射について式で表すと次のようになります。
\begin{equation*}
\Abka{L_o(\vx_1, \vomega)} = L_e(\vx_1, \vomega) +
\begin{cases}
\displaystyle \frac{1}{P_{RR}} \frac{f_s(\vx_1, \vomega_{i, 1}, \vomega) L_i(\vx_1, \vomega_{i, 1}) \absb{\vomega_{i, 1} \! \cdot \vec{n}_1}}{\frac{1}{4\pi}} & P_{RR} \mbox{の確率} \\
0 & 1 - P_{RR} \mbox{の確率}
\end{cases}
\end{equation*}
この処理をロシアンルーレット (Russian Roulette)と呼びます。これによって無限にトレースすることなく、結果の期待値は真値に一致したままにできます。
透過物質では反射する光の経路と透過する光の経路が考えられますが、パストレーシングでは、基本的にどちらか一方にのみレイをトレースします。両方のレイをトレースすると、長い経路長を考えるほどレイの本数が爆発的に増加してしまいます。この反射か透過どちらかを選択する処理もロシアンルーレットと呼ばれます。この場合も期待値を真値に一致したままにするために、選ばれた一方の寄与を選ぶ確率で割る必要があります。
ここまでの説明では、方向のサンプリングとして一様分布にそったものを例としてきましたが、既知の情報を有効活用することで結果の分散を低減することができます。理想は被積分関数である $ f_s(\vx, \vomega_i, \vomega) L_i(\vx, \vomega_i) \absb{\vomega_i \!\! \cdot \vec{n}} $ に完全に比例した入射方向のサンプリングですが、$ L_i $ が不明なのでこれは不可能です。そこで、よく使用される最適化として、BSDFやBSDFとcos項の積に沿ったサンプリングがあります。
\begin{equation*}
\Abka{L_o(\vx_1, \vomega)} = L_e(\vx_1, \vomega) + \frac{f_s(\vx_1, \vomega_{i, 1}, \vomega) L_i(\vx_1, \vomega_{i, 1}) \absb{\vomega_{i, 1} \! \cdot \vec{n}_1}}{p_\sigma (\vomega_{i, 1})}
\end{equation*}
例として完全拡散反射面($ f_s = \rho / \pi $, $ \rho $ は反射率)に関する重点的サンプリングを考えます。完全拡散面は入射したエネルギーが半球内のあらゆる方向に等しい放射輝度をもって反射します。反射の場合はあくまで「半球内」に等しく広がるため、透過側から光が入射したとしてもこの場合のBSDFの値は0となり、「大きな値ほど高い確率でサンプリングする」という重点的サンプリング的な観点から、全立体角 $ 4\pi $ からのサンプルは明らかに無駄です。では、完全拡散面の場合は半球内から一様にサンプリングするのが良いかというと、そうではなくまだ最適化の余地があります。
半球内の方向から入射する放射束 $ = L_i(\vx_1, \vomega_{i, 1}) \absb{\vomega_{i, 1} \! \cdot \vec{n}_1} $ が同じなのであれば、拡散面の場合、どの方向をサンプルしても寄与が同じとなるため、半球内の一様サンプリングが最適となります。しかし実際には $ L_i $ が未知であるため、入射放射束が一定ということは事前にはわかりません。ただ、cos項の存在があるので、面に対して垂直に近い入射方向ほど高い寄与を持つ可能性が高いことは言えます。
そこで、完全拡散面に対してはBRDFとcos項の積に比例した入射方向のサンプリング( $ p_\sigma (\vomega) = \cos{\theta} \; / \; \pi $ )がよく使用されます(図4)。
(a) 純粋なパストレーシング | (b) Next Event Estimation |
図5. Next Event Estimationの効果 |
画像を生成するためには、光源とセンサーを結ぶ完全な光輸送経路が必要です。しかし、ここまでに説明した手法では、視点からひたすらレイをトレースしてたまたま光源にヒットした時だけしか完全な経路を生成できません。さらに光源が小さくなればなるほどレイが光源にヒットできる確率が小さくなっていくため、結果の分散が増大してしまいます。また、前述のBSDFにそった重点的サンプリングなどを行うことにより、むしろ光源面にヒットできる確率が減ってしまう場合もあります。そこで、実用的なパストレーシングでは光源上の点を明示的にサンプルして、視線経路の端点と接続を行う手法Next Event Estimation (NEE)が使用されます。
図5にNEEを使うことによる結果の改善についてレンダリング画像を載せています。両者は同等の計算時間で得たレンダリング結果ですがNEEを使用することで劇的に分散が低減していることがわかります。ただし両者はもっと時間をかければ同じ結果に収束していきます。
視点からレイをトレースする過程で、衝突点ごとにNEEを実行します。
最初の衝突点 $ \vx_1 $ を例に考えます。まず光源上の点 $ \vx_l $ が確率的に1点選ばれ、その時点での視線経路の端点 $ \vx_1 $ から $ \vx_l $ に対してシャドウレイ(shadow ray)をトレースします。シャドウレイが障害物にあたらなかった場合は寄与を蓄積します。光源上の点のサンプリング方法にも色々ありますが、例えば単純に面の中から一様分布でサンプリングするような場合、確率密度は面積に関して表される(=単位は $ \Brk{\mathrm{m}^{-2}} $)ほうが直感的でしょう。この場合は次の式のように長さ2の経路の寄与を計算します。
\begin{equation}
\Abka{L_{o, 2}(\vx_1, \vomega)} = \frac{f_s(\vx_l \RAR \vx_1 \RAR \vx_0) L_e(\vx_l \RAR \vx_1) G(\vx_l \LRAR \vx_1)}{p_A(\vx_l)} \label{eq:NEE}
\end{equation}
ここで、$ p_A(\vx_l) $ は光源の面積中から $ \vx_l $ を選択する確率密度です。$ G(\vx_l \LRAR \vx_1) = \absb{\vomega_{1l} \cdot \vec{n}_1} \absb{\vomega_{l1} \cdot \vec{n}_l} \;/\; \norm{\vx_l - \vx_1}^2 $ は現在の点と光源上の点との間の幾何項です($ \vomega_{1l} $, $ \vomega_{l1} $ は2点を結ぶ方向ベクトル、$ \vec{n}_{1} $, $ \vec{n}_{l} $ はそれぞれの点における法線)。シャドウレイが障害物にあたった場合は、現在の点と光源上の点を結ぶ光輸送経路がそもそも存在しないことになるので寄与の蓄積は行いません。この寄与の計算は「レンダリング方程式2」で紹介した3点形式のレンダリング方程式に基づいています。ここでは確率密度を面積に関して表すほうが自然だったので幾何項 $ G $ を含む積分(積分変数の単位が面積)で考えましたが、もちろん同じ計算を幾何項の代わりにcos項を含む積分(積分変数の単位が立体角)に関して考えることも可能です。上記の面積測度に関する確率密度 $ p_A $ は $ p_\sigma(\vomega_{1l}) = p_A(\vx_l) \cdot \norm{\vx_l - \vx_1}^2 \;/\; \absb{\vomega_{l1} \cdot \vec{n}_l} $ として立体角測度に関する確率密度 $ p_\sigma $ に変換することができます。幾何項の代わりにcos項を含むモンテカルロ推定量を考えると次のようになります。
\begin{equation*}
\Abka{L_{o, 2}(\vx_1, \vomega)} = \frac{f_s(\vx_l \RAR \vx_1 \RAR \vx_0) L_e(\vx_l \RAR \vx_1) \absb{\vomega_{1l} \cdot \vec{n}_1}}{p_\sigma (\vomega_{1l})} = \frac{f_s(\vx_l \RAR \vx_1 \RAR \vx_0) L_e(\vx_l \RAR \vx_1) \absb{\vomega_{1l} \cdot \vec{n}_1} \absb{\vomega_{l1} \cdot \vec{n}_l}}{p_A(\vx_l) \norm{\vx_l - \vx_1}^2}
\end{equation*}
この式の右辺は結局式 $ \eqref{eq:NEE} $ と同じになることがわかります。
NEEを使う場合も確率的な入射方向のサンプリングを行います。ただし、光源にヒットできたときの寄与の蓄積は行わず、単に長い光輸送経路を生成する目的に使用されます。例として長さ3の経路の場合には、$ \vx_1 $ でNEEを行って光源上の点をサンプリングするのとは別に入射方向もサンプリングして次の衝突点 $ \vx_2 $ を決定した後、点 $ \vx_2 $ においてNEEを実行します。これによって、視線経路の端点と光源点の間に障害物が無い限り確実に各経路長ごと(長さ1は除く)の寄与を計算することができます。なお、パストレーシングでは一般的にセンサー上の点 $ \vx_0 $ ではNEE (経路長1)を行わないので、最初の衝突点 $ \vx_1 $ が光源だった場合は例外的に寄与(経路長1)の蓄積を行います。
このページではモンテカルロレイトレーシングの基礎となる手法であるパストレーシングの理論について紹介しました。基礎的な手法とはいっても時間さえかければレンダリング方程式の真値に収束していく結果が得られます。「パストレーシング実装」で実装について解説したいと思います。
先の内容も含みますが関連資料としてスライドも貼っておきます。