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


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

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

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

## この記事で触れること
- **マーチングキューブレイのためのシェーダの詳細（数式やアルゴリズム）**

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

# 頂点シェーダの`vRayDirection`
```javascript: mcrVshader.vert.js
vRayDirection = (mvMatrixTranspose * mvMatrix * vec4(position, 1.0)).xyz;
```
頂点シェーダに出てきた上記の式は、格子立方体座標系で見た際に、カメラからのレイの方向がどうなるかを計算しています。結果を`vRayDirection`としてフラグメントシェーダに入力します。

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

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

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

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

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

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

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

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

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

## `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軸ベクトル）

![CubeViewCoord.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152791/59672e61-b7b3-9d14-4395-e97f1afb811a.png)


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

```math:
\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$となります）。

```math:
(\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`の各成分の内積になります。[^2] 第4成分は不要なので、シェーダでは`.xyz`で`vec3`にしています。

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

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

行列には線形性があるため、以下の2つが一致します。
- 頂点のベクトルを線形補間してから行列をかけた結果
- 頂点のベクトルに行列をかけたものを線形補間した結果

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

```math:
\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$は該当ピクセルに映るポリゴン上の点の位置となります。

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

```math:
\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`
```javascript: 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`[^3]）から`rayDirection`の方向へ伸ばし、立方体の外に出るまでの距離にあたります。イメージとしては以下の画像のような感じです。

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

![MaxTImage.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152791/2792aeed-6cf9-c7ec-cc1b-68fbc1607842.png)

やってることとしては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`関数のどちらか片方のみが使用されます。

# フラグメントシェーダでのレイマーチング
```javascript: 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回の処理の繰り返しで最悪のケース[^4]でも誤差を立方体を通るレイの長さの$(\frac{1}{8})^{\frac{24}{3}} = \frac{1}{512}$に抑えることができます。移動量を立方体を通るレイの長さの1/512にして、レイが立方体の外に出るまで調べようとした場合は512回の繰り返し処理が必要になるため、だいぶ処理を軽くできます。

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

# フラグメントシェーダでの法線の計算
```javascript: 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$のバイリニア補間のような式になります。このため、偏微分値の計算を頂点ごとのトリリニア補間で計算しても問題なさそうと判断しています。

![TrilinearInterpolation.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152791/2a3fbf1a-d26f-6fd0-f592-474da140d55c.png)

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

```math:
\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$の偏微分式は割愛します。）
```math:
\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})
}
```

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

![WrongNormalResult.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152791/d00c583b-d931-40af-e147-61e02526079a.png)

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