ボクセルに対しレイトレーシングを行う

グリッド上に並んだボクセルに対してレイトレを行う必要が出てきたので、今後同じようなことをする際に混乱しないためのメモです。

記法

• 実数 x に対し、\operatorname{floor}(x) で x の切り捨て（x を超えない最大の整数）を表します。
• 実数 x に対し、\operatorname{fract}(x) で x 小数部（正確には x-\operatorname{floor}(x)）を表します。

ボクセルのローカル座標系に持っていく

そのままでは計算し辛過ぎるので、とりあえずボクセルのローカル座標系にレイを持っていきます。ローカル座標系では、

• 各ボクセルの大きさは1
• ボクセルはXYZ座標軸に平行に並ぶ
• ボクセルの角は整数座標にある

とします。早い話がマイクラと同じです。

ここから、レイの始点を \bm p=(p_x,p_y,p_z)、正規化されたレイの進行方向を \bm d=(d_x,d_y,d_z) としてレイトレを行います。レイは単位時間あたりに距離1だけ進むものとして、時刻 t におけるレイの位置を \bm p+t\bm d で表すことにします。

また、説明のために時間の単位を秒にしておきます。つまり、レイは1秒に距離1だけ進み、t 秒後には座標 \bm p+t\bm d にいることになります。

やりたいことは次の通りです。

• 進めたレイが新たに交差するボクセル全てについて、交差する順に次の情報を得る：
  • 進入開始時刻
  • 進入するボクセルの法線ベクトル

座標軸ごとに分離する

次に、ボクセルの境界面を座標軸ごとに分離して考えます。

• X軸に垂直な方向に広がる平面（YZ平面）が間隔1で並んだもの
• Y軸に垂直な方向に広がる平面（XZ平面）が間隔1で並んだもの
• Z軸に垂直な方向に広がる平面（XY平面）が間隔1で並んだもの

これらを全て重ねると無限に広がるグリッドの境界が形成されます。そこで、ひとまず間隔1で並んだX軸に垂直な平面に対してだけ判定を行うことを考えてみます。面倒を避けるため、以下では d_x>0 の場合を考えます。

1秒経つごとにレイのX座標は d_x だけ変化するので、レイは 1/d_x 秒ごとに1回ずつ平面と交わることになります。レイが最初に平面と（新たに）交わるまでに、X軸方向に 1-\operatorname{fract}(p_x) 進む必要があるので、レイが平面と交差する時刻は

\frac{1-\operatorname{fract}(p_x)}{d_x}+N_x\frac1{d_x} = \left\{\frac{1-\operatorname{fract}(p_x)}{d_x},\frac{2-\operatorname{fract}(p_x)}{d_x},\dots\right\}

と列挙できます。ここで、N_x は0以上の整数です。同様に、Y軸とZ軸に垂直な平面と交差する時刻についても、0以上の整数 N_y, N_z を用いてそれぞれ次のように列挙できます。ただし、同様に d_y,d_z>0 を仮定します。

\frac{1-\operatorname{fract}(p_y)}{d_y}+N_y\frac1{d_y}\\
\frac{1-\operatorname{fract}(p_z)}{d_z}+N_z\frac1{d_z}

これらの時刻において、レイが何らかのボクセルに進入することが分かりました。

衝突イベントを処理する

望む情報を得るためには、上記の衝突時刻を小さい順に処理していく必要があります。これは次のような処理によって実現できます。

進入時刻を得る

現在の時刻 t と、各座標軸について「あと何秒で次の平面と交わるか」を表す変数 s_x, s_y, s_z を用意し、次にように初期化します。

t &= 0 \\
s_x &= (1-\operatorname{fract}(p_x))/d_x \\
s_y &= (1-\operatorname{fract}(p_y))/d_y \\
s_z &= (1-\operatorname{fract}(p_z))/d_z

次に、s_x, s_y, s_z のうち最も小さいものについて、その方向からレイがボクセルに進入したことを報告し、時刻 t をその値だけ進め、s_x, s_y, s_z の値を更新します。これを満足するまで^[1]繰り返します。

t, s_x, s_y, s_z の更新は、例えば s_x が最小の場合は次のように行います。

1. t \gets t+s_x（時間を進める）
2. s_y \gets s_y-s_x（残り時間を減らす）
3. s_z \gets s_z-s_x（残り時間を減らす）
4. s_x \gets 1/d_x（次の衝突時刻を設定する）

進入したボクセルを特定する

次に、各反復において進入するボクセルを特定することを考えます。

中心が (x+0.5,y+0.5,z+0.5) にあるボクセルを整数の組 (x,y,z) で表すことにします。レイが現在所属しているボクセルの座標を表す変数 g_x, g_y, g_z を用意し、次のように初期化しておきます。

g_x &= \operatorname{floor}(p_x) \\
g_y &= \operatorname{floor}(p_y) \\
g_z &= \operatorname{floor}(p_z)

後は、イベントが発生する度に g_x, g_y, g_z の値を変化させることで次に進入するボクセルを特定できます。

例えば、s_x が一番小さかった場合は、X軸方向にボクセルの境界を跨いだことになるので g_x の値を1増やします。更新後の (g_x,g_y,g_z) がレイの進入するボクセルです。

衝突法線を得る

最後に各イベントについて進入方向に対応するボクセルの衝突法線を得ることを考えますが、これはレイの跨いだ平面の法線に対応するので、

• s_x が最小のときは (-1,0,0)
• s_y が最小のときは (0,-1,0)
• s_z が最小のときは (0,0,-1)

と分かります。

ゼロ除算の回避

d_x=0 のようなケースが発生した場合、上記の手順通りに計算するとゼロ除算が発生してしまうため、フラグを立ててその軸に関する計算を回避するか、s_x=10^9 のように巨大な値を設定しておくことで問題を回避できます。

おわり

以上で必要な情報が揃ったため、後は好きなようにレイトレすることができます。実装の際は、進行方向について d_x,d_y,d_z>0 が仮定されていることだけ注意してください。正にしておいて後で結果を修正するか、各種の変分に符号を掛けることで対応できます。