CUDAを用いたシンプルなパストレーシング
gdep
CUDAを用いたシンプルなパストレーシング
岩崎 慶 先生
和歌山大学 システム工学部 准教授、プロメテックCGリサーチ 研究員
1. はじめに
CG分野において、画像を生成する手法の一つにレイトレーシングがあります。
レイトレーシングは、例えば視点から半直線の光線（レイ）を追跡（トレース）することで、視点に到達する光のエネルギー（輝度）を求めることができます。
図1に、視点からレイを放射して、シーン中の物体との交点を求めるレイトレーシングの例を示します。
レイトレーシングで必要となる、レイと物体（通常は三角形メッシュで表現されます）との交差計算は、レイごとに独立して計算することができるため、GPUによる並列処理に適しています。
そのため、近年では、NVIDIA の RTX Technology や Microsoft の DirectX Ray Tracing など、GPUを用いた高速なレイトレーシングが実装・開発されており、次世代の3Dグラフィクス技術として注目を浴びています。
赤と緑の画素を通るレイは独立なため、並列に処理できる。
本稿では、物体を球のみに限定し、NVIDIAの並列プログラミングモデルであるCUDAを用いて、画像を生成する手法について解説します。
2. レイトレーシング
まず基本的なレイトレーシングを、CUDAを用いて実装してみます。
はじめに、各画素の輝度値（RGB3成分単精度）を格納するためのバッファを用意します。
今回、GPU側でレイトレーシングを実行して各画素の輝度値を求め、その結果をCPU側でビットマップファイルとして保存することを考えます。
そのため、GPUとCPUの両方からアクセスできるメモリであるUnified Memoryを利用します。
Unified Memoryの確保は、cudaMallocManaged関数を用います。
変数widthとheightはそれぞれバッファの横と縦の画素数で、device_bufferがUnified Memoryへのポインタです。
次に、レイが各画素内のどの位置を通るかを決めるために、乱数を用意します。
CUDAを用いて乱数を計算する手法の一つに、cuRANDがあります。cuRANDを用いて一様乱数を生成するには curand_uniform関数を用います。
curand_uniform関数は、curandState型の変数を引数として入力する必要があるため、まずcurandState用のメモリを確保します。
図2↓にアンチエイリアシング（anti-aliasing）を考慮したレイトレーシング用のカーネル関数 render_aa を示します。
カーネル関数 render_aa の引数について説明します。
- pixels は、各画素の値を格納するバッファへのポインタで、先程説明した device_buffer を渡します。
- scene は、球で構成されたシーンへのポインタを表します。
- rand_state は、乱数生成用の状態を格納した配列へのポインタ、eye は視点の3次元座標、ns は1画素あたりのレイのサンプル数を表します。
カーネル関数 render_aa では、各画素に対して1つのスレッドを割り当てて、ns本のレイを追跡します。
77行目から79行目で画素とスレッドの対応づけを行なっています。
x, y がそれぞれ左から x番目、下から y番目の画素を表しており、79行目はpixelsの範囲外へアクセスしないようにしています。
81行目の intersection は構造体で、レイと物体（今回のプログラムでは球）との交点の情報（3次元位置，法線，色）を保持します。
83行目の ray はレイのためのクラスで、レイの始点と単位方向ベクトルをメンバ変数を保持します。
86行目の変数 m_p は1画素の大きさを表します。
89行目から100行目の forループで、ns本のレイを作成します。
90行目から92行目で画素中の点を通る方向ベクトルを計算しています。
前述の通り、画素内の位置をの乱数で決めるため、curand_uniform関数を用いています。
レイの方向を単位ベクトルとして渡すため、94行目の normalize で正規化してレイr を初期化しています。
95行目の intersect関数は、レイと交点情報を引数としてとり、レイがシーン中の物体と交差した場合はtrueを返し、交差しない場合はfalseを返す（sceneクラスの）メンバ関数です。レイr が物体と交差した場合、レイの始点に最も近い交点の情報を変数 isect に記録します。
96行目から98行目では、交点の法線（isect.m_n）を変数Lに足しこみます。
ここで、変数L を輝度値として各画素に格納したいのですが、法線の xyz成分を取りうるため、範囲に収まるように1を足して0.5倍しています。
101行目から103行目で、交点の法線情報を累積した値（変数L）をレイの数 ns で割ることで輝度を平均化し、画素のRGB値として格納します。
図3左に render_aa の実行結果を示します。
このシーンは、7個の球から構成されています。左右、上下、奥の平面にみえる物体は、半径の大きい球の一部です。
右上図3(a)は anti-aliasing 処理をしない場合、右下図3(b)は anti-aliasing 処理をした場合。右上図と比較して、エイリアシングが緩和されています。
図3右の2枚の画像は、render_aa の実行結果を拡大した図です。図3右上の画像は、anti-aliasing 処理を無効にして作成した画像です。
レイの通る画素内の位置を、乱数で決定するのではなく、画素の中心を通るようにして作成した画像です。（図2の 93行目をアンコメントすると作成できます）。
図3右上の画像に比べて、図3右下の画像はエリアシングが緩和されていることがわかります。
3. パストレーシング
物理則に基づいて画像を生成する場合、光源から放たれた光が、シーン中の物体表面で反射などを繰り返して、スクリーンの各画素を通って視点に到達する光の輝度を積分する必要があります。
通常、この積分は解析的に求めることができないため、モンテカルロ法を用いて推定されることが多いです。
モンテカルロ法を用いて視点に到達する光の輝度を求める手法の一つに、パストレーシングがあります。
パストレーシングは、アルゴリズムが比較的単純で、サンプル数を増やせば物理的に正しい画像に収束することが知られています。
そのため、Sony Pictures Imageworks の Arnold、 Pixar の RenderMan、 Weta DigitalのManuka、Disney の Hyperion Renderer など、プロダクションレベルでのレンダラ（画像生成プログラム）で使用されています。
パストレーシングは、光源から視点へ向かう光のエネルギーを、光の経路（パス）をサンプリングすることによって求めます。
最もシンプルなパストレーシングでは、視点からレイを追跡し、光源と交差するまで経路を追跡します。
図4にパストレーシングの例を示します。
視点 x0 から緑色の画素を通るレイを生成し、レイトレーシングによりレイと物体表面との交点 x1 を計算します。
交点 x1 が光源上にないため、あらたな方向をサンプリングし、交点 x1 をレイの始点として再びレイトレーシングし交点 x2 を計算します。
視点から光源上の点への経路をサンプリング
交点 x2 も光源上にないため、再び方向をサンプリングしレイトレーシングにより交点 x3 を計算します。
交点 x3 は光源上の点のため、ここで経路追跡を終了します。
経路 x-=x0x1x2x3 に沿って視点 x0 に到達する光のエネルギーは、光源上の点 x3 の放射輝度に、交点 x1x2 での（余弦項で重み付けされた）反射率を乗算することで計算されます。
各画素の色は、サンプリングされた経路のエネルギーを、その経路のサンプリング確率密度で割った値の平均値で求められます。
図5に、パストレーシングを用いて画像生成するカーネル関数 render を示します。
画素とスレッドとの対応など、基本となる箇所は前述の render_aa 関数と類似しているため、パストレーシングに関する箇所のみ抜粋しています。
234行目で、各交点での反射率をサンプリング方向の確率密度で割った値を格納する変数 tp を初期化しておきます。
236行目で、視点からのレイがシーン中の物体と交差するか判定します。
238行目で、交点が光源上にあるかどうかを is_emissive 関数で判定します。
視点からのレイが直接光源と交差した場合は、光源の放射輝度を変数Lに足しこみます。光源と交差しない場合、経路追跡をはじめます。
本来、光源と交差するまで経路を追跡するか、ロシアンルーレットにより経路追跡を終了すべきですが、今回は簡単のため、forループにより経路の長さが10を超えたら経路追跡を終了します。
241行目から262行目が経路追跡処理に相当します。
今回は物体表面の材質が全て光沢反射すると仮定します（拡散反射の場合は245行目と246行目をアンコメントし、248行目から250行目をコメントアウトすると実行できます）。
249行目で、交点における光の反射分布に応じて方向をサンプリングし、変数dに格納しています。
またこの時、サンプリングした方向の確率密度を変数 pdf_w に格納します。
250行目で、eval_ggx 関数によって計算される反射率と、dot（d, isect.m_n）で計算される余弦項の積を、確率密度 pdf_w で割った値を、変数 tp に乗算します。
交点（isect.m_p）とサンプリングした方向 d を用いてレイ r を初期化して、レイトレーシングを行います。
255行目において、レイが交差した点が光源上にあるか判定し、光源上であれば変数 tp と光源の放射輝度 Le を乗算して変数 L に足しこみ、経路追跡を終了します。 交点が光源上でなければ経路追跡を継続します。
図6に、パストレーシングによる画像生成例を示します。
図6左は、シーン中の物体表面の材質が全て拡散反射すると仮定して生成した画像で、図6右は、光沢反射を仮定して生成した画像です。
どちらの画像も、左右の赤色と青色の球体で反射した間接光によるカラーブリーディングや映り込みがみてとれます。
図6：カーネル関数 render の実行例．左図は物体表面の材質をすべて拡散反射する材質として生成した画像
右図はすべて光沢反射として生成した画像。どちらもサンプル数を4096としている。
4. おわりに
本稿では、CUDAを用いてシンプルなレイトレーシングとパストレーシングの実装例を紹介しました。
レイトレーシングやパストレーシングは、レイごとに独立した処理であるため、GPUによる並列処理が比較的しやすいアルゴリズムといえます。
今回はCUDAでの簡易実装に焦点をあて、自前でレイと球との交差判定を計算しましたが、交差判定部分を RTX Technologies や DirectX Ray Tracing の API に置き換えることも可能です。
パストレーシングでは、視点側から経路を生成しますが、視点側と光源側の双方から経路を生成する手法として双方向パストレーシングがあります。
我々の研究グループでは、双方向パストレーシングの新しい手法を提案しています（図7参照）。
こちらの手法は現時点でCPUによる並列処理を行なっていますが、 今後CUDAを利用して高速化することも考えられます。
【参考サイト】
本稿で紹介したCUDAのプログラムは以下のサイトによりご覧いただけます。
https://github.com/kiwasaki/simpleptgpu/
また、我々の研究グループが開発した新しい双方向パストレーシングの技術論文と簡易版プログラムは以下のサイトによりご覧いただけます。
http://web.wakayama-u.ac.jp/~iwasaki/project/risbpt/
http://web.wakayama-u.ac.jp/~iwasaki/project/tsrbpt/
https://github.com/kiwasaki/simple_ris_bpt/
図7：双方向パストレーシングによる画像
著者紹介
岩崎 慶 先生
和歌山大学 システム工学部 准教授
プロメテックCGリサーチ 研究員
1999年 東京大学理学部情報科学科卒業
2001年 同大学大学院新領域創成科学研究科博士前期課程修了
2004年 同大学大学院新領域創成科学研究科博士後期課程修了
同年 和歌山大学システム工学部情報通信システム学科助手
2007年 同講師
2009年 同准教授．科学博士．主としてコンピュータグラフィクスに関する研究に従事
