[自由研究] WebGL2で陰関数モデルをレンダリングしてみた の補足

N.Mです。この記事は、以前投稿したWebGL2での陰関数レンダリングについての補足記事です。

内容が多かったので、記事を2つに分けちゃいます。まだご覧になられていない方は先に以下の記事をご覧ください。
[自由研究]WebGL2で陰関数モデルをレンダリングしてみた

導入

以前の記事では、陰関数レンダリングをする際に、WebGL2を使用した実装を中心にまとめました。その中で、シェーダ内の細かい数式やアルゴリズムについては、ほとんど触れなかったため、今回その部分を解説しようと思います。

この記事で触れること

• マーチングキューブレイのためのシェーダの詳細（数式やアルゴリズム）

この記事で触れないこと

• 前回の記事と内容のかぶる話
  全体で何をやろうとしているのか、マーチングキューブレイ（と勝手に名づけた手法）がそもそも何なのかなどについては、以前の記事をご覧ください。

頂点シェーダのvRayDirection

mcrVshader.vert.js

vRayDirection = (mvMatrixTranspose * mvMatrix * vec4(position, 1.0)).xyz;

頂点シェーダに出てきた上記の式は、格子立方体座標系で見た際に、カメラからのレイの方向がどうなるかを計算しています。結果をvRayDirectionとしてフラグメントシェーダに入力します。

mvMatrixやmvMatrixTransposeが何なのかなどを説明するには、モデル変換、ビュー変換の理解が必要です。
（プロジェクション変換はこの数式では不要。）

モデル変換、ビュー変換

頂点シェーダでは、頂点バッファオブジェクトなどからの情報をもとに頂点の位置を決めますが、頂点位置の計算処理はシェーダに自前で書く必要があります。頂点バッファオブジェクトからのデータをそのまま使っただけでは、遠近感のある立体的なレンダリングを行うことはできません。

そこで、WebGLに限らずグラフィックスでは以下の3つの変換がよく使用されます。この話はwgld.orgさんの記事でも解説されているので、合わせてご覧ください。（記事を書く際の参考にもしています。）

• モデル変換
  モデルに対して、どのような移動、回転、スケールをかけたいかを表す変換です。
  同じモデルに対して異なるモデル変換をかけて描画することで、違った位置、向き、大きさのモデルを効率的に複数描画することができます。

• ビュー変換
  カメラ中心を原点に、カメラの奥行方向をz軸に、カメラの上方向をy軸に設定した座標系（ビュー座標系）にする変換です。
  この変換により、カメラから見てモデルがどの位置にあるかが求まります。

• プロジェクション変換
  ビュー変換で求まった座標をWebGLなどの描画環境に合わせた座標にする変換です。
  手前のものが大きく、奥のものが小さく映るような遠近感を出すための透視投影変換を主に行います。WebGLではxyzの各成分が-1~+1に収まる頂点を描画しますが、その範囲に収めるよう変換するクリッピング処理もここで行います。

これらの変換は4×4の行列として表現され、頂点位置を示すベクトルにかけることで変換が行われます。また順番も重要で モデル変換→ビュー変換→プロジェクション変換 の順でかける必要があります。頂点位置を列ベクトル¹として表した場合は、頂点シェーダで計算される頂点位置は以下のようになります。

（プロジェクション変換行列）×（ビュー変換行列）×（モデル変換行列）×（頂点位置）

vRayDirectionの計算

では本題に戻ります。今回の頂点シェーダに出てくるmvMatrixは（ビュー変換行列）×（モデル変換行列）のことです。またmvMatrixTransposeはmvMatrixの行列を転置したものです。（転置とは行列のm行n列の成分をn行m列の成分に移動させる行列の処理です。）

vRayDirectionを求めるために以下のことをやっています。

1. ビュー座標系でのレイの方向を求める。レイが頂点positionを通る場合、カメラの中心が原点であるため、mvMatrix * positionで求まります。
2. 格子立方体座標系のxyz各軸のベクトルが、ビュー座標系だとどの方向になるかを求めるため、これらのベクトルにmvMatrixをかける。
3. 1のビュー座標系でのレイの方向に対し、2のモデル変換・ビュー変換後の各軸ベクトルとの内積をとり、これをvRayDirectionの各成分とする。
   • vRayDirectionのx成分 = （1のレイ方向ベクトル）・（2のモデル・ビュー変換後のx軸ベクトル）
   • vRayDirectionのy成分 = （1のレイ方向ベクトル）・（2のモデル・ビュー変換後のy軸ベクトル）
   • vRayDirectionのz成分 = （1のレイ方向ベクトル）・（2のモデル・ビュー変換後のz軸ベクトル）

2についてですが、モデル・ビュー変換後の各軸ベクトル $\boldsymbol{u}_x, \boldsymbol{u}_y, \boldsymbol{u}_z$ はmvMatrixを $\boldsymbol{VM}$とすると

\displaylines{
\boldsymbol{u}_x = \boldsymbol{VM} \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \end{pmatrix} \\
\boldsymbol{u}_y = \boldsymbol{VM} \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \end{pmatrix} \\
\boldsymbol{u}_z = \boldsymbol{VM} \begin{pmatrix} 0 \\ 0 \\ 1 \\ 0 \end{pmatrix} \\
}

となり、$\boldsymbol{u}_x$はmvMatrixの1列目、$\boldsymbol{u}_y$はmvMatrixの2列目、$\boldsymbol{u}_z$はmvMatrixの3列目となります。

3の計算はmvMatrixの転置を利用して以下のように表現できます。mvMatrixの4列目を$\boldsymbol{u}_w$、モデル・ビュー変換前の頂点positionは$\boldsymbol{p}_b$、変換後を$\boldsymbol{p}_a$とします（つまり、$\boldsymbol{p}_a = \boldsymbol{VM} \boldsymbol{p}_b$となります）。

(\boldsymbol{VM})^T \boldsymbol{VM} \boldsymbol{p}_b = \begin{pmatrix} \boldsymbol{u}^T_x \\ \boldsymbol{u}^T_y \\ \boldsymbol{u}^T_z \\ \boldsymbol{u}^T_w \end{pmatrix} \boldsymbol{p}_a
= \begin{pmatrix} \boldsymbol{u}_x \cdot \boldsymbol{p}_a \\ \boldsymbol{u}_y \cdot \boldsymbol{p}_a \\ \boldsymbol{u}_z \cdot \boldsymbol{p}_a \\ \boldsymbol{u}_w \cdot \boldsymbol{p}_a \end{pmatrix}

この式の左部分がシェーダで使用している式で、右部分の上から3つの成分がvRayDirectionの各成分の内積になります。² 第4成分は不要なので、シェーダでは.xyzでvec3にしています。

vRayDirectionの正規化

vRayDirectionの正規化は次のフラグメントシェーダで行っています。

行列には線形性があるため、以下の2つが一致します。

• 頂点のベクトルを線形補間してから行列をかけた結果
• 頂点のベクトルに行列をかけたものを線形補間した結果

フラグメントシェーダには、ポリゴンを構成する頂点のデータを線形補間した結果が入力として送られます。この場合、線形性から、以下の等式が成り立ちます。（$\boldsymbol{p}_0$, $\boldsymbol{p}_1$, $\boldsymbol{p}_2$はポリゴンの頂点位置、$w_0$, $w_1$, $w_2$を頂点を線形補間する際のウェイトとします。）

\displaylines{
w_0(\boldsymbol{VM})^T \boldsymbol{VM} \boldsymbol{p}_0 + w_1(\boldsymbol{VM})^T \boldsymbol{VM} \boldsymbol{p}_1 + w_2(\boldsymbol{VM})^T \boldsymbol{VM} \boldsymbol{p}_2 \\
= (\boldsymbol{VM})^T \boldsymbol{VM}(w_0\boldsymbol{p}_0 + w_1\boldsymbol{p}_1 + w_2\boldsymbol{p}_2)
}

この$w_0\boldsymbol{p}_0 + w_1\boldsymbol{p}_1 + w_2\boldsymbol{p}_2$は該当ピクセルに映るポリゴン上の点の位置となります。

頂点シェーダで正規化してしまうと、正規化後のベクトルを線形補間することになります。線形補間する前に正規化をかけてしまうと上記の線形性が崩れてしまいます。要は

\displaylines{
w_0 * normalize((\boldsymbol{VM})^T \boldsymbol{VM} \boldsymbol{p}_0) + w_1 * normalize((\boldsymbol{VM})^T \boldsymbol{VM} \boldsymbol{p}_1) \\
+ w_2 * normalize((\boldsymbol{VM})^T \boldsymbol{VM} \boldsymbol{p}_2) \\
\not= normalize((\boldsymbol{VM})^T \boldsymbol{VM}(w_0\boldsymbol{p}_0 + w_1\boldsymbol{p}_1 + w_2\boldsymbol{p}_2))
}

ということです。このため、正規化はフラグメントシェーダで行う必要があります。

フラグメントシェーダのmaxT

mcrFshader.frag.js

vec3 rayOrigin = clamp(vRayOrigin, 0.0, 1.0);
vec3 rayDirection = normalize(vRayDirection);

// infinityにはあらかじめ適当な大きい値 (10000.0) を入れています。
float maxT = infinity;
for (int idx = 0; idx < 3; idx++) {
  float signDir = step(0.0, rayDirection[idx]);

  // epsilonには0除算を回避するために、適当な小さな値 (0.0001) を入れています。
  float candidateT = mix(rayOrigin[idx], 1.0 - rayOrigin[idx], signDir) / (abs(rayDirection[idx]) + epsilon);

  // rayDirectionの成分が0.0に近い場合はcandidateTの値をinfinityとして扱い、無効にしています。
  maxT = min(maxT, mix(infinity, candidateT, step(epsilon, abs(rayDirection[idx]))));
}

// 基本的にはないはずですが、maxTがinfinityのまま変わらなかった場合や負の値になってしまった場合の対応を書いています。
maxT = mix(0.0, maxT, step(epsilon, abs(infinity - maxT)));
maxT = max(maxT, 0.0);

フラグメントシェーダのこの部分ではmaxTの計算をしています。このmaxTは、レイが立方体にあたった位置（rayOrigin³）からrayDirectionの方向へ伸ばし、立方体の外に出るまでの距離にあたります。イメージとしては以下の画像のような感じです。

やってることとしてはrayDirectionを正規化した後、以下の条件すべてを満たすTの最大値をもとめ、これをレイが立方体を抜けない上限maxTとしています。各条件に対し、上限となるTが計算されますが、そのなかで最も小さいものがすべての条件を満たすTの最大値です。

• 0.0 ≦ rayOrigin.x + T * rayDirection.x ≦ 1.0
• 0.0 ≦ rayOrigin.y + T * rayDirection.y ≦ 1.0
• 0.0 ≦ rayOrigin.z + T * rayDirection.z ≦ 1.0

rayDirectionの成分の正負によって式が変わるため、正ならsignDirが1.0, 負ならsignDirが0.0になるようにstep関数で分岐しています。このsignDirの値で式が変わるように線形補間をするためのmix関数を使用しています。signDirは0.0か1.0のどちらかであるため、mix関数のどちらか片方のみが使用されます。

フラグメントシェーダでのレイマーチング

mcrFshader.frag.js

vec3 currentRayPos = rayOrigin;

// divideNumには最初の分割数である8を入れています。
float invDivideNum = 1.0 / float(divideNum);
float oneStepT = maxT * invDivideNum;
float currentT = 0.0;

float hitFlag = 0.0;
// iterNumには繰り返し回数である24を入れています。
for (int idx = 0; idx < iterNum; idx++) {
  // currentValueが現在のレイの位置での陰関数の値で、その正負 (currentSign) で内外判定を行います。
  float currentValue = trilinearInterpolation(currentRayPos, funcValueMinusZ, funcValuePlusZ) - isosurfaceValue;
  float currentSign = step(0.0, currentValue);
  hitFlag = max(hitFlag, 1.0 - currentSign);
  currentT -= mix(oneStepT, 0.0, currentSign);
  oneStepT *= mix(invDivideNum, 1.0, currentSign);
  currentT = clamp(currentT + oneStepT, 0.0, maxT);

  // 計算精度の問題か、currentTの範囲を制限してもレイの位置が立方体の外に出てしまうことがあるようだったので、さらにclampで制限しています。
  currentRayPos = clamp(rayOrigin + currentT * rayDirection, 0.0, 1.0);
}

このコードでは少しずつレイを進め、レイと陰関数モデルとの衝突点がどこかを調べています。最初の移動量oneStepTは立方体内を通るレイの長さmaxTの1/8としています。現在のレイの位置が陰関数モデルの内部にある場合は、currentSignが0.0になります。この場合はhitFlagを1.0とし、いったんレイの位置を前ステップの位置に戻して、レイの移動量をさらに1/8にしています。

こうすることで、少ない繰り返し回数で衝突位置の誤差を小さくすることができます。今回の場合だと、24回の処理の繰り返しで最悪のケース⁴でも誤差を立方体を通るレイの長さの$(\frac{1}{8})^{\frac{24}{3}} = \frac{1}{512}$に抑えることができます。移動量を立方体を通るレイの長さの1/512にして、レイが立方体の外に出るまで調べようとした場合は512回の繰り返し処理が必要になるため、だいぶ処理を軽くできます。

フラグメントシェーダでの法線の計算

mcrFshader.frag.js

vec3 surfaceNormal;
surfaceNormal.x = trilinearInterpolation(currentRayPos, funcNormalXMinusZ, funcNormalXPlusZ);
surfaceNormal.y = trilinearInterpolation(currentRayPos, funcNormalYMinusZ, funcNormalYPlusZ);
surfaceNormal.z = trilinearInterpolation(currentRayPos, funcNormalZMinusZ, funcNormalZPlusZ);
surfaceNormal = normalize(surfaceNormal);

この行の通り、各成分での関数の偏微分値を立方体の頂点ごとにあらかじめ計算し、トリリニア補間しています。関数値のトリリニア補間の式を$x$で偏微分すると$y$と$z$のバイリニア補間のような式に、$y$で偏微分すると$z$と$x$のバイリニア補間のような式に、$z$で偏微分すると$x$と$y$のバイリニア補間のような式になります。このため、偏微分値の計算を頂点ごとのトリリニア補間で計算しても問題なさそうと判断しています。

上の図のように各頂点の関数値が設定されている場合、トリリニア補間の式は以下のようになります。（0 ≦ $x, y, z$ ≦ 1です。）

\displaylines{
f(x, y, z) = (1-x)(1-y)(1-z)v_{000} + x(1-y)(1-z)v_{100} \\
+ (1-x)y(1-z)v_{010} + xy(1-z)v_{110} \\
+ (1-x)(1-y)zv_{001} + x(1-y)zv_{101} \\
+ (1-x)yzv_{011} + xyzv_{111}
}

これを例えば$x$で偏微分すると、以下のようになり、$y$と$z$でのバイリニア補間をしたような式になります。（似たような式になるので、$y$と$z$の偏微分式は割愛します。）

\displaylines{
\frac{\partial f}{\partial x} = (1-y)(1-z)(v_{100} - v_{000}) + y(1-z)(v_{110} - v_{010}) \\
+ (1-y)z(v_{101} - v_{001}) + yz(v_{111} - v_{011})
}

このトリリニア補間の式を偏微分した結果を、法線の成分に使用することも考えられます。しかし、隣接する立方体で頂点での関数値が異なるため、立方体の境界で偏微分の値も不連続となります。結果として、下の画像のような格子の境界がはっきり見えるような結果になってしまいます。

まとめ

前回の記事で説明しなかった、シェーダ内の詳細のアルゴリズムや数式について解説してみました。マーチングキューブレイの追実装をしなくとも、この記事がシェーダで何かを実装する際の参考になればうれしいです。

1. 変換を4×4行列として表現しているので、頂点位置も4成分必要です。グラフィックスではxyzの座標に、wの成分を加えた同次座標を使用し、透視投影はこの同次座標の性質を利用しています。詳細の話はこの記事では割愛しますが、wgld.orgさんではここで解説されています。 ↩

2. positionには第4成分として1を入れていますが、モデル変換行列、ビュー変換行列の性質上、$\boldsymbol{u}_x$, $\boldsymbol{u}_y$, $\boldsymbol{u}_z$の第4成分が0であるので、第4成分の積は無視できます。 ↩

3. rayOriginでclampを入れているのは、頂点シェーダで線形補間されるvRayOriginが0.0~1.0の範囲をわずかに超えてしまうのを回避するためです。これがなかった場合、一部の状況で描画がおかしくなってしまいました。 ↩

4. レイが立方体を抜けるぎりぎりの位置で衝突するようなケースの場合。1/8だけ進める処理を8回、次に1/64だけ進める処理を8回、最後に1/512だけ進める処理を8回行うようなケースです。 ↩