この記事はレイトレアドベントカレンダー2023 ⧉の記事として書かれました。
@Shocker_0x15です。最近Tessellation-Free Displacement Mappingという手法を実装したのですが、そこで使われていた三角形上の関数へのアフィン演算の適用が面白かったので記事を書いてみます。

精度保証付き演算

floatやdoubleといった浮動小数点数は有限の精度しか持たないため、ごく一部の数値を除いて近似的にしか表現できません。そのため長い計算プロセスを経ると誤差が蓄積・増大してしまい計算結果の信頼性に問題が出ます。そこで特定の数値のかわりに、表現したい数値を確実に含む値の範囲を何かしらの形式で表現、それらの範囲に対する各種演算も併せて定義することで通常の数値計算の結果に加えて結果の有りうる誤差も同時に評価する手法が開発されています。このような手法を精度保証付き演算と呼びます。

精度保証付き演算の理論の詳細や実装についてはこのページでは簡単な紹介に留めますが参考にした文献をページ末尾に載せておきます。また、アフィン演算(や区間演算)に関してはちょうどレイトレ合宿9 ⧉で@ykozw88 ⧉さんがセミナー ⧉を行っていたのでこちらも参考にしてみてください。

区間演算

精度保証付き演算の最もシンプルな実現方法として区間演算(Interval Arithmetic)と呼ばれる手法があります。区間演算では表現したい数値を含む下限と上限のペアを「区間」として表現します。例えば $ a $ という実数を表現する区間 $ A $ は次のように表されます。
$$
A = \left[ \underline{A}, \overline{A} \right]
$$
ここで $ \underline{A} $, $ \overline{A} $ は $ a $ を2つの間に含む計算機上で表現可能な実数値です。$ a \in \left[ \underline{A}, \overline{A} \right] $ とも表せます。そして区間同士の二項演算や区間に対する単項演算が、「入力区間に含まれるすべての実数値の組み合わせに演算を適用した結果すべてを含む区間」を返すように定義します。例えば区間 $ A $ と $ B $ の加算は次のように定義します。
$$
A + B = \left[ \underline{A} \underline{+} \underline{B}, \overline{A} \overline{+} \overline{B} \right]
$$
$ \underline{+} $, $ \overline{+} $ は下向き・上向きに丸めモードを設定した上で加算を実行することを表しています。C++であれば std::fesetround などで調べると良いでしょう。2つの区間の下限同士、上限同士の加算が入力区間中のあり得る組み合わせ全てに対して下限・上限となるのは明らかですね。一方減算は次のようになります。
$$
A - B = \left[ \underline{A} \underline{-} \overline{B}, \overline{A} \overline{-} \underline{B} \right]
$$
2つの区間の減算として最も小さくなるのは $ A $ の下限から $ B $ の上限を引いた場合ですね。上限の計算も同様に加算と比較して関係が反転しています。乗算や除算についてはもう少しややこしくなっていて、例えば乗算は次のようになります。
$$
A \cdot B = \left[ \min \left( \underline{A} \underline{\cdot} \underline{B}, \underline{A} \underline{\cdot} \overline{B}, \overline{A} \underline{\cdot} \underline{B}, \overline{A} \underline{\cdot} \overline{B} \right), \max \left( \underline{A} \overline{\cdot} \underline{B}, \underline{A} \overline{\cdot} \overline{B}, \overline{A} \overline{\cdot} \underline{B}, \overline{A} \overline{\cdot} \overline{B} \right) \right]
$$
符号の正負によって大小が反転するためすべての組み合わせに対して最小・最大を考える必要があります。加減乗除以外にも区間に対する様々な単項演算($ \sqrt{} $ など)を定義します。プログラミングの観点では演算子や関数のオーバーロードを使うことで通常の数値と同様の感覚で所望の計算を実装できます。

区間演算の欠点

区間演算は非常に簡単な原理なのですが、区間がすぐに大きくなりすぎるという欠点があります。数値間の誤差の相関を考慮に入れていないことが原因です。極端な例として同じ区間 $ A $ 同士の減算を考えてみます。
$$
A - A = \left[ \underline{A} \underline{-} \overline{A}, \overline{A} \overline{-} \underline{A} \right]
$$
意味的には同じ数値を表現している区間同士を減算すれば答えは完璧にゼロになって欲しいところですが、区間演算はそれを考慮に入れていないため無駄に区間が広がってしまいます。この欠点を解消する手法としてアフィン演算(Affine Arithmetic)があります。

アフィン演算

アフィン演算ではある数値を次のように表現します。
$$
[A] = a_0 + a_1 \epsilon_1 + \ldots + a_n \epsilon_n
$$
ここで $ \epsilon_i $ は数値の変動を表す $ -1 \leq \epsilon_i \leq 1 $ の範囲を持つ変数(しかし何かを解けば本当の数値がわかるような通常の意味での変数ではありません、ダミー変数とここでは呼びます)、$ a_i $ はそれらにかかる係数です。$ a_0 $ のみ単独で存在しています。ダミー変数と係数の線形結合として表すためアフィン形式(Affine form)と呼びます。ダミー変数の数はアフィン形式に対する演算の過程で増えることもあります。実装的には係数を保持します。

計算の最初に与える入力値それぞれ(が独立であるならば)には異なるダミー変数を割り当ててアフィン形式を初期化します。例えば二次元の関数 $ z=f(x, y) $ を考えるとき $ x $ と $ y $ それぞれに対応するアフィン形式を次のように定義します。

\begin{align}
[x] &= x_c + x_r \epsilon_x + 0 \epsilon_y \\
[y] &= y_c + 0 \epsilon_x + y_r \epsilon_y \\
\end{align}

ここで $ x_c $, $ y_c $ はそれぞれの値の変動範囲の中心値、$ x_r $, $ y_r $ は変動の「半径」を表します(ダミー変数 $ \epsilon $ の範囲が $ -1 \leq \epsilon \leq 1 $ であるため)。$ x $ の値はダミー変数 $ \epsilon_x $ にあわせて変動しますが、$ [y] $ の持つ $ \epsilon_x $ の係数はゼロであり、また $ [x] $ の持つ $ \epsilon_y $ の係数もゼロであるため、これら2つの数値のあいだには相関が一切無いこと意味しています。

アフィン形式同士の加算・減算やアフィン形式の定数倍は通常の文字式と全く同様(丸め誤差については一旦無視)にして計算することができます。したがって区間演算の欠点で述べたような同じ数同士の減算は次のように綺麗にゼロになってくれます。
$$
[x] - [x] = (x_c - x_c) + (x_r - x_r) \epsilon_x = 0
$$

乗算や減算については詳細は示しませんがもちろん演算が定義されます。ただしダミー変数の2次以上の項が出ないように線形近似を行います。その際の近似誤差を考慮するため計算結果には新たなダミー変数が付加されます。非線形演算以外にも丸め誤差を別途考慮する場合にまた別のダミー変数が付加されます。アフィン演算ではこのように計算の過程でダミー変数が増えていくことにより計算負荷が高くなりがちな欠点があるのですが、誤差の推定性能とトレードオフでダミー変数の増加を抑える方法も考えられています。

区間演算・アフィン演算の学習がてら実装したリポジトリのリンクを参考までに貼っておきます。
(比較的シンプルな実装にとどめているので見通しは良いかもしれませんが、あっている保証はありません。歴史ある分野なので本当に調べたい場合は既存の歴史あるオープンソースライブラリなどを見るべきです。)
https://github.com/shocker-0x15/AffineArithmetic

三角形上の関数

三角形自体は三次元空間中の図形として表されたりしますが、三頂点の位置が与えられていれば三角形上の位置は二次元のパラメターで表すことができます。具体的には三角形の重心座標(barycentric coordinates)やテクスチャー座標が使われるでしょう。

重心座標

もう一つよく使われるパラメターとして三角形の重心座標があります。重心座標表現では三角形の各頂点 $ A $, $ B $, $ C $ それぞれのウェイトとして3つの数値 $ b_A $, $ b_B $, $ b_C $ を使います。3つの重心座標の間には
$$
b_A + b_B + b_C = 1
$$
という条件があるため、重心座標は実質二次元のパラメターです。例えば三角形が含まれる平面中の位置 $ \boldsymbol{p} $ は三頂点の位置 $ \boldsymbol{p}_A $, $ \boldsymbol{p}_B $, $ \boldsymbol{p}_C $ を使って
$$
\boldsymbol{p}(b_A, b_B, b_C) = b_A \cdot \boldsymbol{p}_A + b_B \cdot \boldsymbol{p}_B + b_C \cdot \boldsymbol{p}_C
$$
として表します。点 $ \boldsymbol{p} $ が三角形内にある場合、重心座標は必ず $ 0 \leq b_A, b_B, b_C \leq 1 $ となります。重心座標はここに示した位置を始めとして、三頂点が持つ属性値を補間する用途に使われます。プログラム上で保持するのは二頂点目と三頂点目の値 $ b_B $ と $ b_C $ が多い気がします。

テクスチャー座標

三角形メッシュは多くの場合、物体表面のマテリアルや法線・高さなどを表すテクスチャーを貼り付けるため、各位置と合わせてテクスチャー座標 $ u $, $ v $ を持っています。テクスチャー座標は位置とは独立に設定できるため縮退(degenerated, 2つ以上の頂点が同じ座標を持っている状態)してしまうことがあるものの、縮退が無ければ三頂点間の補間として求めたあるテクスチャー座標と三角形中の位置は一対一対応であるため、テクスチャー座標は三角形上の位置を特定するための二次元のパラメターとみなすことができます。

テクスチャー座標の関数への変換

三角形上の位置や法線などは重心座標の関数として表される一方、テクスチャーは当然テクスチャー座標の関数として表されます。三角形上のパラメターの関数を考える場合、いずれかのパラメター表現に揃える必要があります。重心座標に関して表す場合はテクスチャー座標 $ \boldsymbol{t} = (u, v)^T $ を上で示した位置と同様、次のように三頂点のテクスチャー座標 $ \boldsymbol{t}_A $, $ \boldsymbol{t}_B $, $ \boldsymbol{t}_C $ の補間として表すことで全体を重心座標の関数として表現できます。
$$
\boldsymbol{t}(b_B, b_C) = b_A \cdot \boldsymbol{t}_A + b_B \cdot \boldsymbol{t}_B + b_C \cdot \boldsymbol{t}_C
$$
逆にテクスチャー座標に関する三角形上の関数を考える場合、テクスチャー座標から重心座標への変換が必要になります。まず上の重心座標からテクスチャー座標への変換を行列 $ \mathcal{T} _{bc}^{tc} $ を用いて表します。

\begin{pmatrix}
\boldsymbol{t} \\
1
\end{pmatrix}
=
\begin{pmatrix}
u \\
v \\
1
\end{pmatrix}
=
\mathcal{T} _{bc}^{tc}
\begin{pmatrix}
b_A \\
b_B \\
b_C
\end{pmatrix}

$ \mathcal{T} _{bc}^{tc} $ は具体的には次のような内容です。

\mathcal{T} _{bc}^{tc} = 
\begin{pmatrix}
u_A & u_B & u_C \\
v_A & v_B & v_C \\
1 & 1 & 1
\end{pmatrix}
=
\begin{pmatrix}
\boldsymbol{t}_A & \boldsymbol{t}_B & \boldsymbol{t}_C \\
1 & 1 & 1
\end{pmatrix}

$ \mathcal{T} _{bc}^{tc} $ の逆行列 $ \mathcal{T} _{tc}^{bc} = \left( \mathcal{T} _{bc}^{tc} \right)^{-1} $ を使えばテクスチャー座標から重心座標へ変換できます。すなわち次の式を考えます。

\begin{pmatrix}
b_A \\
b_B \\
b_C
\end{pmatrix}
=
\mathcal{T} _{tc}^{bc}
\begin{pmatrix}
u \\
v \\
1
\end{pmatrix}
=
\begin{pmatrix}
u_A & u_B & u_C \\
v_A & v_B & v_C \\
1 & 1 & 1
\end{pmatrix}^{-1}
\begin{pmatrix}
u \\
v \\
1
\end{pmatrix}

三角形上の関数へのアフィン演算の適用

ここでは例としてハイトマップによるディスプレイスメントを適用した三角形の位置を考えてみます。変位適用後の表面の位置 $ \boldsymbol{S} $ は次の式で表されます。
$$
\boldsymbol{S}(u, v) = \boldsymbol{p}(u, v) + h(u, v) \cdot \boldsymbol{\hat{n}}(u, v)
$$
ここで $ h(u, v) $ は $ u $, $ v $ に対応するハイトマップテクスチャーの値、$ \boldsymbol{p}(u, v) $ は重心座標補間した位置、$ \boldsymbol{\hat{n}}(u, v) $ は次の示すように重心座標補間した法線 $ \boldsymbol{n}(u, v) $ を正規化したものです。
$$
\boldsymbol{\hat{n}}(u, v) = \frac{\boldsymbol{n}(u, v)}{|| \boldsymbol{n}(u, v) ||}
$$
$ \boldsymbol{p} $ と $ \boldsymbol{n} $ は具体的には次の形になります。

\boldsymbol{p}(u, v)
=
\begin{pmatrix}
\boldsymbol{p}_A & \boldsymbol{p}_B & \boldsymbol{p}_C \\
\end{pmatrix}
\begin{pmatrix}
b_A \\
b_B \\
b_C
\end{pmatrix}
=
\begin{pmatrix}
\boldsymbol{p}_A & \boldsymbol{p}_B & \boldsymbol{p}_C \\
\end{pmatrix}
\mathcal{T} _{tc}^{bc}
\begin{pmatrix}
u \\
v \\
1
\end{pmatrix}

\boldsymbol{n}(u, v)
=
\begin{pmatrix}
\boldsymbol{n}_A & \boldsymbol{n}_B & \boldsymbol{n}_C \\
\end{pmatrix}
\begin{pmatrix}
b_A \\
b_B \\
b_C
\end{pmatrix}
=
\begin{pmatrix}
\boldsymbol{n}_A & \boldsymbol{n}_B & \boldsymbol{n}_C \\
\end{pmatrix}
\mathcal{T} _{tc}^{bc}
\begin{pmatrix}
u \\
v \\
1
\end{pmatrix}

純粋なテクスチャー座標に代わりに、テクスチャー座標の動く範囲を表した二次元のアフィン形式(実装上は三次元目に固定の1)を使えば $ p $, $ n $ に対応するアフィン形式が得られます。アフィン形式の法線ベクトルの正規化演算も定義することができます。あとはテクスチャー座標の範囲に対応するハイト値の(保守的な)範囲を何らかの方法で求めアフィン形式で表現すればそれらをまとめて $ \boldsymbol{S} $、つまり変位適用後の表面の位置の範囲を表す三次元のアフィン形式が得られます。アフィン形式のままだと範囲として少し扱いづらいため、これを区間に変換すればxyzそれぞれの範囲、つまりAABBを求めることができます。区間に変換することで誤差の範囲は広がってしまいますが、最初から全て区間演算で計算した場合に比べて誤差の増大を抑えることができます。

ところで精度保証付き演算というと計算誤差などの基本的には小さい範囲の入力を扱うイメージがありましたが、ここではガッツリ動く範囲を入力として与えています。理論的には意図せず生まれる誤差の範囲も意図的に与える入力の範囲も同じですね。

テクスチャー座標のアフィン形式

上述のアフィン形式の紹介ではスカラー値に関する形式を示しましたが、テクスチャー座標のようなベクトルに関してもアフィン形式を定義できます。各次元が互いに完全に相関しているのなら単に係数がベクトルになるだけです。したがってベクトル $ \boldsymbol{v} = (v_x, v_y, v_z)^T $ に対応するアフィン形式は次のようになります。

\left[
\begin{pmatrix}
v_x \\
v_y \\
v_z
\end{pmatrix}
\right]
=
\begin{pmatrix}
v_{x, 0} \\
v_{y, 0} \\
v_{z, 0}
\end{pmatrix}
+
\begin{pmatrix}
v_{x, 1} \\
v_{y, 1} \\
v_{z, 1}
\end{pmatrix}
\epsilon_1
+
\ldots
+
\begin{pmatrix}
v_{x, n} \\
v_{y, n} \\
v_{z, n}
\end{pmatrix}
\epsilon_n

アフィン演算は前述の通り推定性能と計算コストのトレードオフとしてダミー変数をどれだけ保持するかという観点がありました。ここでは専用のダミー変数はテクスチャー座標 $ u $, $ v $ に対応する $ \epsilon_u $, $ \epsilon_v $ のみとし、$ u $, $ v $ に無関係の計算過程で発生する近似はすべて単一のダミー変数 $ \epsilon_K $ で面倒を見ることにします。すなわち次のようなアフィン形式を扱います。

\left[
\begin{pmatrix}
v_x \\
v_y \\
v_z
\end{pmatrix}
\right]
=
\begin{pmatrix}
v_{x, c} \\
v_{y, c} \\
v_{z, c}
\end{pmatrix}
+
\begin{pmatrix}
v_{x, u} \\
v_{y, u} \\
v_{z, u}
\end{pmatrix}
\epsilon_u
+
\begin{pmatrix}
v_{x, v} \\
v_{y, v} \\
v_{z, v}
\end{pmatrix}
\epsilon_v
+
\begin{pmatrix}
v_{x, K} \\
v_{y, K} \\
v_{z, K}
\end{pmatrix}
\epsilon_K

入力となるテクスチャー座標のアフィン形式の作り方の例として三角形中のテクセルに対応する範囲と、三角形全体の範囲を以下に示します。

三角形中のテクセル

テクセルの形状はテクスチャー空間内ではuv軸に沿った矩形(上図中左)なので簡単です。テクセルの中心 $ (u_c, v_c) $ を基準として $ \epsilon_u $, $ \epsilon_v $ の係数にu軸v軸に沿ったテクセルの幅 $ w_t $ の半分の大きさのベクトルを設定するだけです。

\left[
\begin{pmatrix}
u \\
v \\
1
\end{pmatrix}
\right]
=
\begin{pmatrix}
u_c \\
v_c \\
1
\end{pmatrix}
+
\begin{pmatrix}
0.5 \cdot w_t \\
0 \\
0
\end{pmatrix}
\epsilon_u
+
\begin{pmatrix}
0 \\
0.5 \cdot w_t \\
0
\end{pmatrix}
\epsilon_v

三角形全体

三角形全体に対応するテクスチャー座標のアフィン形式を表現したい場合、三角形全体を囲むテクスチャー空間中での長方形を考えるのが一つの手ではありますが、三角形外のテクスチャー座標の範囲も多く含んでしまうため値の範囲を過大評価してしまいます。代わりに三角形を3つの(重複する)平行四辺形の和(上図中右)として考えることで厳密に三角形に沿ったアフィン形式を考えることが提案されています。例として頂点 $ A $ をひとつの端点として持つ平行四辺形を考えるとアフィン形式は次のようになります。

\left[
\begin{pmatrix}
u \\
v \\
1
\end{pmatrix}
\right]
=
\begin{pmatrix}
0.5 \cdot u_A + 0.25 \cdot u_B + 0.25 \cdot u_C \\
0.5 \cdot v_A + 0.25 \cdot v_B + 0.25 \cdot v_C \\
1
\end{pmatrix}
+
\begin{pmatrix}
0.25 \cdot (u_B - u_A) \\
0.25 \cdot (v_B - v_A) \\
0
\end{pmatrix}
\epsilon_u
+
\begin{pmatrix}
0.25 \cdot (u_C - u_A) \\
0.25 \cdot (v_C - v_A) \\
0
\end{pmatrix}
\epsilon_v

平行四辺形の中心を基準として $ \epsilon_u $, $ \epsilon_v $ の係数に平行四辺形を構成する2つの方向に対応するベクトルを適切な大きさで設定します。

Tessellation-Free Displacement Mapping

最後にこの記事を書くきっかけとなったTessellation-Free Displacement Mapping (TFDM)を極々簡単に紹介します。TFDMはハイトマップを適用した高密度なジオメトリを事前テッセレーションを行うことなく省メモリ(かつそこそこ高速)にレイトレーシングでレンダリングすることを実現する手法です。

TFDMでは事前計算として、ハイトマップの各テクセルごとの最小値最大値(バイリニアサンプリングなどを行っているとテクセル内でも値が変動します)を記録したテクスチャーを作り、さらに2x2ピクセルごとに再帰的に最小値最大値でダウンサンプリングを行うことでハイトマップの最小値最大値を階層的に記録したminmaxミップマップと呼ばれるテクスチャーを作ります。

レイトレーシングのAS構築のためにはディスプレイスメントマッピングを適用する三角形ごとに変位後のサーフェス位置を囲むオブジェクト空間のAABBが必要です。上述の三角形全体に対応するテクスチャー座標のアフィン形式に基づく計算を行い三角形ごとのAABBを求めます。

レイトレーシングの段階ではレイが前述のAABBにヒットするとIntersection Shaderが起動されます。Intersection Shader内ではminmaxミップマップの最も粗いレベルからトラバーサルを始めて、テクセルの範囲とテクセルが持つ値の範囲を使って上述の通りテクセルに対応するAABBを計算、レイとAABBの交差判定を行います。レイがヒットする場合はミップレベルを一段下って再帰的にAABBの計算・テクセルとの交差判定を行っていきます。レイが目的のミップレベルに到達し次第実際の変位サーフェスとの交差判定を行います。ミップマップは暗黙的に四分木構造となっているためレイがあるテクセルの処理を終えた段階で次に訪れるテクセルを簡単に求めることができます。

このようにTFDMでは暗黙的な四分木構造であるminmaxミップマップとオンザフライなAABB計算、レイとの交差判定を組み合わせることで事前テッセレーションを行うことなく非常に省メモリに高密度ジオメトリのレイトレーシングによるレンダリングを実現しています。

なお、TFDMは2021年に発表された比較的最近の手法でしたがすでに同じ著者からRMIP (Rectangular Minmax Image Pyramid)という手法が発表されていたり、先日のアドベントカレンダーの記事でもありましたが@ShinjiOgakiさんによる非線形レイトレーシングなど目下アクティブな領域のようです。

以下にTFDMの実装例のリンクとサンプル画像を貼っておきます。
https://github.com/shocker-0x15/GfxExp

参考文献

区間演算の実装について (1)
http://verifiedby.me/adiary/070
Affine Arithmeticについて
http://verifiedby.me/kv/nv-boost/affine/affine.pdf
アフィン演算
https://speakerdeck.com/ykozw/ahuinyan-suan
Affine Arithmetic and its Applications to Computer Graphics
https://graphics.stanford.edu/~comba/papers/aa-93-09-sibgrapi-paper.pdf
Approximating implicit curves on triangulations with affine arithmetic
https://webdoc.sub.gwdg.de/ebook/serien/e/IMPA_A/720.pdf
Implementation and improvements of affine arithmetic
https://www.jstage.jst.go.jp/article/nolta/6/3/6_341/_article
Tessellation-Free Displacement Mapping for Ray Tracing
https://perso.telecom-paristech.fr/boubek/papers/TFDM/