2
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

WebGL でテッセレーションして、曲面パッチ間が滑らかに繋がった物体を表示(番外編:法線について)

Last updated at Posted at 2019-09-01

番外編

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

[前編]: WebGL でテッセレーションして、曲面パッチ間が滑らかに繋がった物体を表示(1)
[後編]: WebGL でテッセレーションして、曲面パッチ間が滑らかに繋がった物体を表示(2)

そのうち後編では、独自高次元曲面パッチ(便宜上「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 はそれぞれ下図のようなコントロールポイントを持っている」ということを書いています。
beziex-patch-uv.png
なお後編の頂点位置を求める式では、この 2曲面それぞれの位置を計算し、その後これらを混ぜ合わせていたと思います。今回、u、v 方向の偏微分値を求める方法もこれと同じで、各曲面の偏微分値を求め、後で混ぜ合わせます。

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

u 方向の偏微分値の式(曲面 A)

下図(曲面 A)の 4本の黒い折線について考えます。
beziex-patch-u-mod.png
コントロールポイントが 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 u}S(u,v)$ が最終的な「u 方向の偏微分値」ということになります。また同様に $\frac{\partial}{\partial v}S(u,v)$ を求めれば、それが「v 方向の偏微分値」ということです。
そしてこの $\frac{\partial}{\partial u}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 になりますが、ここでは考えません。

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

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

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

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

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

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

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

境界横断導関数(CBD)

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

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

があります。これを「境界横断導関数(Cross Boundary Derivative)」、略して CBD と呼びます。
下図で言えば、青い線がそれぞれの「境界曲線上の点」から延びる CBD です。
beziex-patch-cbd.png
Beziex パッチの偏微分値には、u 方向と v 方向の 2つがあることを前述しましたが、もし境界曲線が u 方向に延びている場合は、CBD は v 方向の偏微分値ということになります。境界曲線が v 方向に延びている時は、その逆です。

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

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

前項の図を例に、「"上の境界曲線" と "下の境界曲線" 上の同じ u 値の点」を上下の CBD を使って結ぶと、下図のような 3次ベジェ曲線となります(図では制御点を結んだ折線で表示)
cbd-bezier-line.png
またこの3次ベジェ曲線を単独で図にすると、下記のようになりますね。
bezier-line-normal.png
さらに上の図のベジェ曲線に対し、P0 と P1 を同位置にしたものが下図です。
bezier-line-zero-length.png
実はこの場合、P0 上での微分ベクトル長が 0 になります。つまり P0 が「(前述したパッチの)下の境界曲線上の点」とすると、(その境界曲線上の)CBD の長さが 0 であるということを意味します。

微分ベクトル長の補正(概念)

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

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

しかし(わずかな長さであっても)CBD のベクトル方向は、適した向きにあるべきです。この "適した向き" とは、どこなのでしょうか?
ここで再度、前項で示した図を見てみましょう。
bezier-line-zero-length.png
この場合、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_{1}$ が同位置にあるとすると、下記のようになります。

\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)$ は、$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_{1}$ が同位置の場合、

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

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

微分ベクトル長の補正(実装)

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

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

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

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

終わりに

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

他にも今回、UI 関連では

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

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

2
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?