番外編

ちょっと前に「WebGL テッセレーションビューワ」を作ったということで、使用した技術等を中心として 2回に渡り投稿しました。

［前編］: WebGL でテッセレーションして、曲面パッチ間が滑らかに繋がった物体を表示（１）
［後編］: WebGL でテッセレーションして、曲面パッチ間が滑らかに繋がった物体を表示（２）

そのうち後編では、独自高次元曲面パッチ（便宜上「Beziex パッチ」と呼称）に関し、「パッチを "分割して出来たポリゴン" の頂点」の頂点位置を求める式について述べています。
ちなみに 3D モデルでは頂点位置と共に "法線" も重要なのですが、後編の中では法線の計算式を割愛しました。

実は割愛したのは、「説明のボリュームが大きくなりそう」というのが一番の理由です。しかし、やっぱり説明して欲しいというご要望を頂いたので、番外編として投稿することにしました。
なお本稿は、上記 2回の記事をご覧になっていることを前提に書いておりますので、まず先にそちらの方を読んでいただければと思います。

デモ＆ソース

毎度で申し訳ありませんが、本稿にもデモとソース（github）のリンクを貼っておきます。

デモ
Beziex_We : (https://beziex.netlify.com/beziex_we/)

デモの使い方
デモと使い方 : (beziex/Beziex_We [github] のところの README.md 内)

ソースコード（github）
beziex/Beziex_We [github] : (https://github.com/beziex/Beziex_We)

法線に関するおさらい

前編にて曲面パッチモデルでは、テッセレーション後のポリゴン頂点に対して、法線データを事前に用意する必要は無いことを述べました。曲面の u方向と v方向の接線（偏微分値）の外積が法線になるからです。

また後編では、「Beziex パッチは 2枚の曲面を混ぜ合わせたものになっており、この 2枚を曲面 A と曲面 B と呼ぶことにすると、曲面 A、B はそれぞれ下図のようなコントロールポイントを持っている」ということを書いています。

なお後編の頂点位置を求める式では、この 2曲面それぞれの位置を計算し、その後これらを混ぜ合わせていたと思います。今回、u、v 方向の偏微分値を求める方法もこれと同じで、各曲面の偏微分値を求め、後で混ぜ合わせます。

というわけで、まずは曲面 A の偏微分値を求めてみましょう。なお偏微分値といっても u方向と v方向の 2種類があって、Beziex パッチでは両者の計算方法が異なります。
そこで先に（曲面 A の）u 方向の偏微分値の式について述べたいと思います。

u 方向の偏微分値の式（曲面 A）

下図（曲面 A）の 4本の黒い折線について考えます。

コントロールポイントが 4個ある方が「3次ベジェ曲線」、7個ある方が「6次ベジェ曲線」ですが、各ベジェ曲線の "導関数" は下記のようになります。

\frac{d}{du}Q_{a0}(u) = 3((1-u)^2(P_{a10}-P_{a00})+2u(1-u)(P_{a20}-P_{a10})+u^2(P_{a30}-P_{a20}))\\

\frac{d}{du}Q_{a1}(u) = 6((1-u)^5(P_{a11}-P_{a01})+5u(1-u)^4(P_{a21}-P_{a11})+10u^2(1-u)^3(P_{a31}-P_{a21})\\
+10u^3(1-u)^2(P_{a41}-P_{a31})+5u^4(1-u)(P_{a51}-P_{a41})+u^5(P_{a61}-P_{a51}))

\frac{d}{du}Q_{a2}(u) = 6((1-u)^5(P_{a12}-P_{a02})+5u(1-u)^4(P_{a22}-P_{a12})+10u^2(1-u)^3(P_{a32}-P_{a22})\\
+10u^3(1-u)^2(P_{a42}-P_{a32})+5u^4(1-u)(P_{a52}-P_{a42})+u^5(P_{a62}-P_{a52}))

\frac{d}{du}Q_{a3}(u) = 3((1-u)^2(P_{a13}-P_{a03})+2u(1-u)(P_{a23}-P_{a13})+u^2(P_{a33}-P_{a23}))\\

ただし本プログラムではシェーダ内での計算量削減のため、

S_{ai0}(u) = 3(P_{a(i+1)0}-P_{ai0})\hspace{50px}(i = 0～2)\\
S_{ai1}(u) = 6(P_{a(i+1)1}-P_{ai1})\hspace{50px}(i = 0～5)\\
S_{ai2}(u) = 6(P_{a(i+1)2}-P_{ai2})\hspace{50px}(i = 0～5)\\
S_{ai3}(u) = 3(P_{a(i+1)3}-P_{ai3})\hspace{50px}(i = 0～2)\\

として、$S_{aij}(u)$ を CPU 側（TypeScript コード）で事前計算しています。それが下記のソースコードです。

BxGlShader.convVtfInfo

private convVtfInfo(hPosBez0: BxBezier3Line3|null, hPosBez1: BxBezier6Line3|null, hPosBez2: BxBezier6Line3|null, hPosBez3: BxBezier3Line3|null,
    vPosBez0: BxBezier3Line3|null, vPosBez1: BxBezier6Line3|null, vPosBez2: BxBezier6Line3|null, vPosBez3: BxBezier3Line3|null,
    hDiffBez0: BxBezier2Line3|null, hDiffBez1: BxBezier5Line3|null, hDiffBez2: BxBezier5Line3|null, hDiffBez3: BxBezier2Line3|null,
    vDiffBez0: BxBezier2Line3|null, vDiffBez1: BxBezier5Line3|null, vDiffBez2: BxBezier5Line3|null, vDiffBez3: BxBezier2Line3|null,
    surfaceOfs: number, vtfInfo: number[]): void {
                                        :
                                        :
    this.fromBezier2(<BxBezier2Line3>hDiffBez0, surfaceOfs, 44, vtfInfo);
    this.fromBezier5(<BxBezier5Line3>hDiffBez1, surfaceOfs, 47, vtfInfo);
    this.fromBezier5(<BxBezier5Line3>hDiffBez2, surfaceOfs, 53, vtfInfo);
    this.fromBezier2(<BxBezier2Line3>hDiffBez3, surfaceOfs, 59, vtfInfo);
                                        :
                                        :
}

さらに GPU 側（シェーダ）では、上記の値を用いて次のように計算を行います。

（shader.vert）getDiffH_U

void getDiffH_U( vec3 hDiffBez0_U[3], vec3 hDiffBez1_U[6], vec3 hDiffBez2_U[6], vec3 hDiffBez3_U[3], float u, float v, out vec3 diffH_U )
{
    vec3  bezV[ 4 ];

    bezV[ 0 ] = posBez2( hDiffBez0_U, u );
    bezV[ 1 ] = posBez5( hDiffBez1_U, u );
    bezV[ 2 ] = posBez5( hDiffBez2_U, u );
    bezV[ 3 ] = posBez2( hDiffBez3_U, u );
                    :
                    :
}

これで、各ベジェ曲線の導関数での「特定の u 値」におけるベクトル値を求めることが出来ました。さらにこの 4個のベクトル値を 3次ベジェ曲線（導関数では無い）で混ぜ合わせることで、「曲面 A での u 方向の偏微分値」を得られます。

\frac{\partial}{\partial u}R_{a}(u,v) = (1-v)^3(\frac{d}{du}Q_{a0}(u))+3v(1-v)^2(\frac{d}{du}Q_{a1}(u))\\
+3v^2(1-v)(\frac{d}{du}Q_{a2}(u))
+v^3(\frac{d}{du}Q_{a3}(u))\\

それを GPU 側で行っているのが、同じく getDiffH_U 関数の次の箇所となるわけです。

（shader.vert）getDiffH_U

void getDiffH_U( vec3 hDiffBez0_U[3], vec3 hDiffBez1_U[6], vec3 hDiffBez2_U[6], vec3 hDiffBez3_U[3], float u, float v, out vec3 diffH_U )
{
                    :
                    :
    diffH_U = posBez3( bezV, v );
}

v 方向の偏微分値の式（曲面 A）

次に、（曲面 A の）v 方向の偏微分値の式について述べます。
まず v 方向の偏微分値を取り出したい場合、u 方向は微分の必要はありません。そのため前項の図にある 4本のベジェ曲線における「特定の u 値」の情報は、そのまま頂点位置として求めます。

Q_{a0}(u) = (1-u)^3P_{a00}+3u(1-u)^2P_{a10}+3u^2(1-u)P_{a20}+u^3P_{a30}\\

Q_{a1}(u) = (1-u)^6P_{a01}+6u(1-u)^5P_{a11}+15u^2(1-u)^4P_{a21}
+20u^3(1-u)^3P_{a31}\\
+15u^4(1-u)^2P_{a41}+6u^5(1-u)P_{a51}+u^6P_{a61}

Q_{a2}(u) = (1-u)^6P_{a02}+6u(1-u)^5P_{a12}+15u^2(1-u)^4P_{a22}
+20u^3(1-u)^3P_{a32}\\
+15u^4(1-u)^2P_{a42}+6u^5(1-u)P_{a52}+u^6P_{a62}

Q_{a3}(u) = (1-u)^3P_{a03}+3u(1-u)^2P_{a13}+3u^2(1-u)P_{a23}+u^3P_{a33}\\

その後、（普通に考えれば）$Q_{a0}(u)$ ～ $Q_{a3}(u)$ を 3次ベジェ曲線のパラメータとして、本ベジェ曲線の（偏）導関数を作った後、「特定の v 値」の偏微分値を求めることになります。

\frac{\partial}{\partial v}R_{a}(u,v) = 3((1-v)^2(Q_{a1}(u)-Q_{a0}(u))+2v(1-v)(Q_{a2}(u)-Q_{a1}(u))\\
+v^2(Q_{a3}(u)-Q_{a2}(u)))\\

しかし本プログラムでは計算量を減らすため、上の式を変形して直接 3次ベジェ曲線のパラメータから偏微分値を取り出すことにしました。

\frac{\partial}{\partial v}R_{a}(u,v) = 3(-(1-v)^2Q_{a0}(u) + (1-v)(1-3v)Q_{a1}(u)\\
+ (1-v)(2-3v)Q_{a2}(u) + v^2Q_{a3}(u))\\

この（曲面 A の）v 方向の偏微分値の式を、バーテックスシェーダ内に実装したものが、次のコードです。

（shader.vert）getDiffV_U

vec3 diffBez3( vec3 bez[4], float t )
{
    float  t2  = t * t;
    float  mt  = 1.0 - t;
    float  mt2 = mt * mt;
    float  d3t = 3.0 * t;

    return 3.0 * ( -mt2                * bez[ 0 ]
                 +  mt * ( 1.0 - d3t ) * bez[ 1 ]
                 +  t  * ( 2.0 - d3t ) * bez[ 2 ]
                 +  t2                 * bez[ 3 ] );
}

void getDiffV_U( vec3 vPosBez0_U[4], vec3 vPosBez1_U[7], vec3 vPosBez2_U[7], vec3 vPosBez3_U[4], float u, float v, out vec3 diffV_U )
{
    vec3  bezU[ 4 ];

    bezU[ 0 ] = posBez3( vPosBez0_U, v );
    bezU[ 1 ] = posBez6( vPosBez1_U, v );
    bezU[ 2 ] = posBez6( vPosBez2_U, v );
    bezU[ 3 ] = posBez3( vPosBez3_U, v );

    diffV_U = diffBez3( bezU, u );
}

曲面 A と曲面 B の偏微分値の混ぜ合わせ

これまでは曲面 A の偏導関数について述べてきましたが、曲面 B も同様に行うことが出来ます。（曲面 B は、曲面 A の $u$ と $v$ を入れ替えたものなので）

次に曲面 A と曲面 B の偏微分値の混ぜ合わせるわけですが、偏微分値には u方向と v方向があるので、以下のように混ぜ合わせます。

「曲面 A の u方向偏微分値」と「曲面 B の u方向偏微分値」を混ぜて、最終的な u方向偏微分値を求める
「曲面 A の v方向偏微分値」と「曲面 B の v方向偏微分値」を混ぜて、最終的な v方向偏微分値を求める

ここでは u方向について述べましょう。（v方向も同様です）

後編で示した通り、グレゴリーパッチ（Brown パッチ）を応用した混ぜ合わせの式は以下の通りでしたね。

S(u,v) = \frac{u(1-u)R_{a}(u,v)+v(1-v)R_{b}(u,v)}{u(1-u)+v(1-v)} \\

今回は $S(u,v)$ の u方向偏微分値を求めたいので、下記のようになります。（式が長くなったので分割しました）

T_{u0}(u,v) = (u(1-u)(\frac{\partial}{\partial u}R_{a}(u,v)) + (-2u+1)R_{a}(u,v) \\
+ v(1-v)(\frac{\partial}{\partial u}R_{b}(u,v)))(u(1-u) + v(1-v)) \\

T_{u1}(u,v) = (u(1-u)R_{a}(u,v)+v(1-v)R_{b}(u,v))(-2u+1) \\

\frac{\partial}{\partial u}S(u,v) = \frac{T_{u0}(u,v) - T_{u1}(u,v)}{(u(1-u) + v(1-v))^2} \\

すなわち $\frac{\partial}{\partial ｕ}S(u,v)$ が最終的な「u 方向の偏微分値」ということになります。また同様に $\frac{\partial}{\partial v}S(u,v)$ を求めれば、それが「v 方向の偏微分値」ということです。
そしてこの $\frac{\partial}{\partial ｕ}S(u,v)$ と $\frac{\partial}{\partial v}S(u,v)$ の外積から、法線が算定されるわけです。
ちなみに偏微分値を求める式の場合でも、四隅が特異点になるので、四隅のみ別の方法で求めます。

なおこの混ぜ合わせもバーテックスシェーダ内で行っており、コードは以下のようになっています。

（shader.vert）getDiffMix_U

vec3 getDiffMixMain_U( vec3 diffH_U, vec3 diffV_U, float u, float v, vec3 posH, vec3 posV )
{
    float  pu = u * ( 1.0 - u );
    float  du = -2.0 * u + 1.0;
    float  pv = v * ( 1.0 - v );

    float  posDenom   = pu + pv;
    float  diffDenomU = du;
    
    vec3  posNumer = pu * posH + pv * posV;

    vec3  diffNumerU0 = pu * diffH_U + du * posH;
    vec3  diffNumerU1 = pv * diffV_U;
    vec3  diffNumerU  = diffNumerU0 + diffNumerU1;

    vec3  numerU = diffNumerU * posDenom - posNumer * diffDenomU;

    return numerU / ( posDenom * posDenom );
}

void getDiffMix_U( vec3 diffH_U, vec3 diffV_U, float u, float v, vec3 posH, vec3 posV, out vec3 diffU )
{
    diffU = ( ( kEps > u || u > ( 1.0 - kEps ) ) && ( kEps > v || v > ( 1.0 - kEps ) ) ) ? diffH_U : getDiffMixMain_U( diffH_U, diffV_U, u, v, posH, posV );
}

偏微分値の長さが 0 の場合は？

ここでもし u 方向または v 方向の偏微分値の片方、もしくは両方の長さが 0 だった時を考えましょう。この場合当然、外積の値も 0 になります。

すなわち法線の長さが 0 になるので、正常なシェーディング状態ではなくなります。それでは困ってしまいますね。本プログラムでは特定の条件でそれが起こりうるので、次項からその条件と対策について述べたいと思います。
ちなみに「u 方向偏微分ベクトルと v 方向偏微分ベクトルの成す角度」が 0 度または 180 度の時も外積（法線長）が 0 になりますが、ここでは考えません。

クローズ／オープン境界曲線とスムース／クリース境界曲線

まず下図をご覧ください。

3D モデルには、穴が開いているものとそうでないものがありますが、上記は穴が開いているモデルです。
このような穴開き曲面パッチモデルには必ず上図の青線のように、以下に示したタイプの境界曲線があります。

境界曲線の片側にしか曲面パッチが存在しない

それに対し、緑の境界曲線は両側に曲面パッチが存在します。ここでは、青い境界曲線のようなタイプを「オープン境界曲線」、緑の境界曲線のようなタイプを「クローズ境界曲線」と呼ぶことにしましょう。
つまり穴が開いていない曲面パッチモデルには、クローズ境界曲線しか無いことになります。

またクローズ境界曲線ではさらにタイプが分かれます。下図を見てください。

青い境界曲線の両側の曲面パッチは、境界曲線を境の折れ曲がっていることが分かります。つまり C0 接続です。それに対し、緑の境界曲線では滑らかに繋がっており、G1 接続であることが分かります。
本稿では、青い境界曲線のようなタイプを「クリース境界曲線」、緑の境界曲線のようなタイプを「スムース境界曲線」と呼ぶことにしました。

ちなみに 3D モデルでは、このように「一部分を角ばらせる」ことや「穴が開いている」ことが必要となる場合も多いと思います。
そして、この「Beziex パッチモデルデータ」を作成した（「ポリゴン → 独自拡張の Catmull–Clark サブディビジョン → Beziex パッチ」に変換する）未公開コンバータは、このような「一部分を角ばらせる」「穴が開いている」にも対応しています。

しかしこのコンバータではアルゴリズムの関係上、「オープン境界曲線」や「クリース境界曲線」上の（ある方向の）偏微分ベクトルの長さが 0 になります。
つまりこの時は法線長が 0 になるので、対策を取らないといけないわけです。

境界横断導関数（CBD）

境界曲線上の偏微分値（偏微分ベクトル）には、境界曲線自身の（偏）微分値の他に

境界曲線と交差する方向の偏微分値

があります。これを「境界横断導関数（Cross Boundary Derivative）」、略して CBD と呼びます。
下図で言えば、青い線がそれぞれの「境界曲線上の点」から延びる CBD です。

Beziex パッチの偏微分値には、u 方向と v 方向の 2つがあることを前述しましたが、もし境界曲線が u 方向に延びている場合は、CBD は v 方向の偏微分値ということになります。境界曲線が v 方向に延びている時は、その逆です。

実は（Beziex パッチの）「オープン境界曲線」「クリース境界曲線」では、CBD のベクトル長が 0 になるのです。よって対策は、この CBD に対して行うことになります。

3次ベジェ曲線端点の（偏）微分ベクトル長が 0 になる条件

前項の図を例に、「"上の境界曲線" と "下の境界曲線" 上の同じ u 値の点」を上下の CBD を使って結ぶと、下図のような 3次ベジェ曲線となります（図では制御点を結んだ折線で表示）

またこの3次ベジェ曲線を単独で図にすると、下記のようになりますね。

さらに上の図のベジェ曲線に対し、P0 と P1 を同位置にしたものが下図です。

実はこの場合、P0 上での微分ベクトル長が 0 になります。つまり P0 が「（前述したパッチの）下の境界曲線上の点」とすると、（その境界曲線上の）CBD の長さが 0 であるということを意味します。

微分ベクトル長の補正（概念）

ここまで、「CBD の長さが 0 であるがゆえに、法線長も 0 になってしまう」ということについて、しつこく述べてきました。逆に言えば、もし CBD が（適した方向に）長さを持つことが出来れば、正常な法線が生成されることになるわけです。

そうであれば、「長さ 0 の CBD」を補正して「わずかに長さを持つ CBD」にしてやればどうでしょう？
ただ（厳密に言えば）微妙に形状が変わってしまうことになりますが、「わずかに長さを持つ」程度であれば、ほとんど誤差レベルです。（元々コンバータによって、Catmull–Clark サブディビジョンから "近似" 変換したデータですし）

しかし（わずかな長さであっても）CBD のベクトル方向は、適した向きにあるべきです。この "適した向き" とは、どこなのでしょうか？
ここで再度、前項で示した図を見てみましょう。

この場合、P0 上にて微分ベクトルが「わずかに長さを持つ」とは、「P0 と P1 をわずかに離す」ということになります。

そして結論から言うと、適した方向とは

P0（元々の P1 の位置）から P2 に向かうベクトルの方向

なのです。すなわち P1 は、P0 と P2 を結ぶ直線上にあれば良いことになります。（ただし P1 は P0 のごく近傍）

微分ベクトル長の補正（証明）

前項の方向が適したものである理由を、簡単に証明してみましょう。まずこの 3次ベジェ曲線の "導関数" は、以下の通りです。

\frac{d}{dt}P(t) = 3((1-t)^2(P_{1}-P_{0})+2t(1-t)(P_{2}-P_{1})+t^2(P_{3}-P_{2}))\\

ここで、$P_{0}$ と $P_{１}$ が同位置にあるとすると、下記のようになります。

\frac{d}{dt}P(t) = 3(2t(1-t)(P_{2}-P_{1})+t^2(P_{3}-P_{2}))\\
= 3t(2(1-t)(P_{2}-P_{1})+t(P_{3}-P_{2}))

また上式を "$3t$" と "$2(1-t)(P_{2}-P_{1})+t(P_{3}-P_{2})$" に分けてみましょう。

Q_{a}(t) = 3t\\

Q_{b}(t) = 2(1-t)(P_{2}-P_{1})+t(P_{3}-P_{2})\\

\frac{d}{dt}P(t) = Q_{a}(t)Q_{b}(t)\\

ここで $Q_{a}(t)$ と $Q_{b}(t)$ は、$ｔ=0$ では、

Q_{a}(0) = 0\\

Q_{b}(0) = 2(P_{2}-P_{1})\\

となります。さらに $P_{0}=P_{1}$ なので、

Q_{b}(t) = 2(P_{2}-P_{0})\\

です。すなわち $P_{0}$ と $P_{１}$ が同位置の場合、

$t=0$ での微分ベクトルは、$2(P_{2}-P_{0})$ ベクトルを 0 倍したもの

と考えられます。つまり（長さは 0 ですが）方向は「$P_{0}$ から $P_{２}$ に向かうベクトル」ということになりますね。

微分ベクトル長の補正（実装）

実はこの補正処理、https://github.com/beziex/Beziex_We の初版では（未公開の）コンバータ内で行っていました。しかし今回、本稿を寄稿するに当たり、

本プログラムに補正処理を追加
初版では「（コンバータで）補正済みの gzjson ファイル」を置いていたが、「補正前の gzjson ファイル」に変更

というコミットを行ったわけです。（ちなみに前編、後編に記載したコードは、今回のコミットで変わってないはず）

なおこの補正処理は最初に 1回だけやれば良いものなので、CPU 側で行っています。具体的には、BxCmCorrectSeparator という TypeScript で記述したクラスで実行するようにしました。（ちょっと長いので、本稿にコードは載せません）

終わりに

やっぱり「法線（偏微分値）の計算式」の説明は長くなってしまいましたね。でも、何とか発表出来て良かったです。

他にも今回、UI 関連では

Vue.js + TypeScript + WebGL (not three.js)

という、あまりネットにも載ってない組み合わせで実装しているので、要望があれば寄稿したいな、と思っています。
でも本プログラムは、初心者向けでは無いからな～（説明しづらいコードなのです）