平面と NURBS 曲面:交差・最短最遠距離・境界箱
— 有理ベジエパッチへの分解と凸包性・ベジエクリッピング —
nurbs_plane_distance.html で使っている計算方法を、式の導出から順に説明します。第 1〜10 節は平面からの距離の極値と曲面全体の境界箱、第 11〜13 節は NURBS のセグメント分解・セグメントごとの包含箱・平面との交差(断面)です。共通するポイントは「有理関数のまま扱わず、同次座標の多項式(ベジエ関数)と凸包性で処理する」ことです。
1有理ベジエ曲面と同次座標
NURBS 曲面は、ノット挿入によって「多項式区間ごとの有理ベジエパッチ」に分解できます(コードの surfaceToBezier)。以下では 1 枚のパッチだけを考えます。
以下、曲面のパラメータを \(u,v\)(\(0\le u,v\le1\))とします。次数 \(p_u \times p_v\)、制御点 \(\mathbf P_{ij}\)、重み \(w_{ij}>0\) の有理ベジエ曲面は
$$ \mathbf S(u,v)=\frac{\displaystyle\sum_{i=0}^{p_u}\sum_{j=0}^{p_v} w_{ij}\,\mathbf P_{ij}\,B_i^{p_u}(u)\,B_j^{p_v}(v)} {\displaystyle\sum_{i=0}^{p_u}\sum_{j=0}^{p_v} w_{ij}\,B_i^{p_u}(u)\,B_j^{p_v}(v)}, \qquad B_i^{p}(u)=\binom{p}{i}u^i(1-u)^{p-i} $$です。分母を \(W(u,v)\) と書きます。\(w_{ij}>0\) なので \(W(u,v)>0\) です。
同次座標 \(\mathbf H_{ij}=(w_{ij}x_{ij},\,w_{ij}y_{ij},\,w_{ij}z_{ij},\,w_{ij})\) を使うと、分子と分母はどちらも普通の(多項式の)ベジエ曲面になります。分割(de Casteljau 法)も同次座標のまま線形補間すれば正しく行えます。
2平面からの符号付き距離は有理関数になる
平面を \(\Pi:\ \mathbf n\cdot\mathbf p=d\)(\(\mathbf n\) は単位法線)とします。曲面上の点の符号付き距離は
$$ f(u,v)=\mathbf n\cdot\mathbf S(u,v)-d . $$\(\mathbf n\cdot\mathbf S\) に第 1 節の式を代入すると、内積は和の中に入れられるので
$$ \mathbf n\cdot\mathbf S(u,v)=\frac{A(u,v)}{W(u,v)},\qquad A(u,v)=\sum_{i,j} \underbrace{w_{ij}\,(\mathbf n\cdot\mathbf P_{ij})}_{a_{ij}}\,B_i(u)B_j(v),\qquad W(u,v)=\sum_{i,j} w_{ij}\,B_i(u)B_j(v) $$となります。\(A\) と \(W\) はそれぞれ制御値 \(a_{ij}\)、\(w_{ij}\) をもつスカラーのベジエ関数です(\(a_{ij}\) は同次座標の第 1〜3 成分と \(\mathbf n\) の内積で、コードでは n.x*h[0]+n.y*h[1]+n.z*h[2])。したがって
3最短・最遠距離は g の最小値・最大値だけで決まる
パッチ全体(NURBS 曲面全体)での \(g\) の最小値を \(g_{\min}\)、最大値を \(g_{\max}\) とします。\(g\) は連続なので、曲面上の点の高さは区間 \([g_{\min},g_{\max}]\) をすべて取ります。距離は \(|g-d|\) なので、次のように場合分けできます。
| 平面の位置 | 最短距離 | 最遠距離 |
|---|---|---|
| \(d\le g_{\min}\)(曲面の下) | \(g_{\min}-d\)(最小点) | \(\max(g_{\max}-d,\ d-g_{\min})\) |
| \(d\ge g_{\max}\)(曲面の上) | \(d-g_{\max}\)(最大点) | |
\(g_{\min}| \(0\)(交線 \(g=d\) 上の全点) | |
4極値の候補:隅・境界・内部
パラメータ領域は正方形 \([0,1]^2\) です。閉じた領域上の滑らかな関数の最小・最大は、次のどれかで生じます。
境界上の極値は、内部の条件 \(\nabla g=\mathbf 0\) を満たすとは限りません(境界に沿った方向の微分だけが 0 になる)。そのため ① と ② を別に調べる必要があります。
5商の微分から多項式の方程式へ(有理関数を消す)
内部の停留点
商の微分公式より
$$ \frac{\partial g}{\partial u}=\frac{A_u W - W_u A}{W^2},\qquad \frac{\partial g}{\partial v}=\frac{A_v W - W_v A}{W^2}. $$\(W>0\) なので分母は決して 0 になりません。したがって、停留点の条件は分子だけで書けます。
次数は次のようになります(\(A,W\) は \(p_u\times p_v\) 次、\(A_u,W_u\) は \((p_u-1)\times p_v\) 次)。
| u 方向の次数 | v 方向の次数 | |
|---|---|---|
| \(G_u=A_uW-W_uA\) | \((p_u-1)+p_u=2p_u-1\) | \(p_v+p_v=2p_v\) |
| \(G_v=A_vW-W_vA\) | \(2p_u\) | \(2p_v-1\) |
実際には最高次の項が打ち消し合い、\(G_u\) の u 方向の真の次数は \(2p_u-2\) です(曲線の場合の \(a'w-w'a\) が \(2n-2\) 次になるのと同じ理由)。コードでは積の公式をそのまま使って \(2p_u-1\) 次で表しています。これは次数を 1 つ上げた表現(次数上げ)と同じ多項式なので、結果に影響はありません。
境界の有理曲線
たとえば \(v=0\) の境界は \(a(u)/w(u)\)(\(a(u)=\sum_i a_{i0}B_i(u)\)、\(w(u)=\sum_i w_{i0}B_i(u)\))という 1 変数の有理関数です。同様に
$$ \frac{d}{du}\frac{a}{w}=\frac{a'w-w'a}{w^2}=0\iff a'(u)\,w(u)-w'(u)\,a(u)=0 . $$これは \(2p_u-1\) 次の 1 変数ベジエ関数の根を求める問題で、有理ベジエ曲線の境界箱と全く同じ処理です。
非有理(すべての重みが等しい)場合との比較
重みがすべて 1 なら \(W\equiv1\)、\(W_u=W_v=0\) なので \(G_u=A_u\)、\(G_v=A_v\) になり、前回のベジエ曲面版の式(\(\partial g/\partial u=\mathbf n\cdot\mathbf S_u\))に一致します。有理版はその自然な拡張です。
6ベジエ関数の微分と積の公式
\(G_u, G_v\) の Bernstein 係数は、次の 2 つの公式だけで計算できます。
微分
$$ \frac{d}{du}\sum_{i=0}^{p} c_i B_i^{p}(u)=\sum_{i=0}^{p-1} p\,(c_{i+1}-c_i)\,B_i^{p-1}(u) $$曲面の u 方向の微分は、v 方向の各列ごとにこれを適用します(コードの nDs, nDt)。
積
\(p\) 次と \(q\) 次の 1 変数ベジエ関数の積は \(p+q\) 次で、係数は
$$ r_k=\sum_{i+j=k}\frac{\binom{p}{i}\binom{q}{j}}{\binom{p+q}{k}}\,c_i\,e_j \qquad(k=0,\dots,p+q). $$テンソル積曲面(\(m\times n\) 次と \(p\times q\) 次)では、各方向に独立に同じ規則を使います。
$$ r_{kl}=\sum_{i+i'=k}\ \sum_{j+j'=l} \frac{\binom{m}{i}\binom{p}{i'}}{\binom{m+p}{k}}\, \frac{\binom{n}{j}\binom{q}{j'}}{\binom{n+q}{l}}\; c_{ij}\,e_{i'j'} $$これらは \(B_i^p B_j^q=\frac{\binom pi\binom qj}{\binom{p+q}{i+j}}B_{i+j}^{p+q}\) という恒等式から直ちに導けます(コードの bz1Prod, nProd)。
7有理ベジエの凸包性(枝刈りの根拠)
\(g\) は次のように書き直せます。
$$ g(u,v)=\frac{\sum a_{ij}B_iB_j}{\sum w_{ij}B_iB_j} =\sum_{i,j}\lambda_{ij}(u,v)\,\frac{a_{ij}}{w_{ij}},\qquad \lambda_{ij}=\frac{w_{ij}B_i(u)B_j(v)}{W(u,v)} . $$\(w_{ij}>0\)、\(B\ge0\) なので \(\lambda_{ij}\ge0\)、さらに \(\sum\lambda_{ij}=1\) です。したがって \(g\) は 比 \(a_{ij}/w_{ij}=\mathbf n\cdot\mathbf P_{ij}\)(制御点の高さ)の凸結合になります。
アルゴリズムではこの不等式を分枝限定に使います。ある部分パッチの \([\min a/w,\ \max a/w]\) が、それまでに見つかった最小値・最大値をどちらも更新しえないなら、その部分パッチはそれ以上調べる必要がありません。
注意:\(G_u\)、\(G_v\) の凸包(次節)は「値の範囲」ではなく「0 になりうる範囲」を判定するためのものです。\(G_u\)、\(G_v\) は多項式なので、普通のベジエ関数の凸包性がそのまま使えます。
8ベジエクリッピングで停留点を求める
連立方程式 \(G_u=G_v=0\) の解(内部停留点)を、パラメータ領域を絞り込みながら探します。
1 方向のクリッピング
\(G_u\) の係数を \(c_{kl}\)(u 方向 \(k=0..K\)、v 方向 \(l=0..L\))とします。点群 \(\{(k/K,\ c_{kl})\}\)(すべての \(l\) について)を平面に置くと、\(G_u(u,v)\) のグラフはどの \(v\) についてもこの点群の凸包の内側にあります。したがって
v 方向は \(G_v\) について同じことを行います。どちらかが交わらなければ、その部分パッチには停留点がありません。
再帰の手順
- 枝刈り:制御値の比 \(a/w\) の範囲が現在の最小・最大を更新しえなければ終了(第 7 節)。
- クリッピング:\(G_u\) から u の範囲、\(G_v\) から v の範囲を求める。どちらかが空なら終了。
- 縮小できた場合(面積比 < 0.64):その範囲の部分パッチを同次座標で切り出し、手順 1 から繰り返す。
- 縮小できない場合(解が複数ありうる):4 分割して、それぞれに手順 1 から適用する。
- 十分小さくなった場合(幅 < 10⁻⁴):ニュートン法で \(G_u=G_v=0\) を精密化し、その点で \(g\) を評価して候補にする。
ニュートン法の候補は「曲面上の点」なので、評価した値は必ず本当の値域の中にあります。精密化がうまくいかない場合でも、最小・最大を誤って広げることはありません。
9具体例:四分円(有理 2 次曲線)
単位円の四分円は、重み付きの 2 次有理ベジエ曲線で正確に表せます。
$$ \mathbf P_0=(1,0),\ w_0=1;\qquad \mathbf P_1=(1,1),\ w_1=\tfrac{\sqrt2}{2};\qquad \mathbf P_2=(0,1),\ w_2=1 . $$方向 \(\mathbf n=(1,1)/\sqrt2\)(45° 方向)の高さ \(g=\mathbf n\cdot\mathbf C(t)\) の最大を求めます。真の答えは円上の 45° の点で、\(g=1\) です。
| 量 | Bernstein 係数 |
|---|---|
| \(\mathbf n\cdot\mathbf P_i\) | \(\tfrac1{\sqrt2},\ \sqrt2,\ \tfrac1{\sqrt2}\) |
| \(a_i=w_i(\mathbf n\cdot\mathbf P_i)\) | \(\tfrac1{\sqrt2},\ 1,\ \tfrac1{\sqrt2}\) |
| \(w_i\) | \(1,\ \tfrac{\sqrt2}{2},\ 1\) |
| \(a'\)(1 次) | \(2-\sqrt2,\ \sqrt2-2\) |
| \(w'\)(1 次) | \(\sqrt2-2,\ 2-\sqrt2\) |
| \(G=a'w-w'a\)(3 次) | \(1,\ \tfrac13,\ -\tfrac13,\ -1\) |
たとえば \(G\) の最初の係数は \(a'_0w_0-w'_0a_0=(2-\sqrt2)\cdot1-(\sqrt2-2)\cdot\frac1{\sqrt2}=1\) です。
\(t=1/2\) で
$$ A(\tfrac12)=\tfrac14\cdot\tfrac1{\sqrt2}+\tfrac12\cdot1+\tfrac14\cdot\tfrac1{\sqrt2}=\tfrac12+\tfrac{\sqrt2}{4},\qquad W(\tfrac12)=\tfrac14+\tfrac12\cdot\tfrac{\sqrt2}{2}+\tfrac14=\tfrac12+\tfrac{\sqrt2}{4}, $$より \(g=A/W=1\) となり、真の答えと一致します。端点では \(g=1/\sqrt2\approx0.707\) なので、最大は内部の停留点です。一方、制御点の高さの最大は \(\mathbf n\cdot\mathbf P_1=\sqrt2\approx1.414\) です。制御点だけから作る範囲は、実際の曲線よりかなり広くなることがわかります。
10境界箱(バウンディングボックス)
軸平行な境界箱の 6 つの面は、座標平面からの最近点・最遠点そのものです。
| 基準平面 | 法線 \(\mathbf n\) | \(g=\mathbf n\cdot\mathbf S\) | 得られる範囲 |
|---|---|---|---|
| yz 平面 | \((1,0,0)\) | \(x(u,v)=X/W\) | \([x_{\min},x_{\max}]\) |
| xz 平面 | \((0,1,0)\) | \(y(u,v)=Y/W\) | \([y_{\min},y_{\max}]\) |
| xy 平面 | \((0,0,1)\) | \(z(u,v)=Z/W\) | \([z_{\min},z_{\max}]\) |
ここで \(X,Y,Z\) は同次座標の第 1〜3 成分のベジエ関数です。たとえば x 方向の停留点方程式は
$$ X_uW-W_uX=0,\qquad X_vW-W_vX=0 $$で、第 5〜8 節の手順が \(\mathbf n\) を座標軸にしただけでそのまま使えます。全パッチについて最小・最大をとれば、NURBS 曲面全体を囲む最小の直方体
$$ [x_{\min},x_{\max}]\times[y_{\min},y_{\max}]\times[z_{\min},z_{\max}] $$が得られます。
11NURBS 曲面の有理ベジエセグメントへの分解
ここまでの議論はすべて「1 枚の有理ベジエパッチ」についてのものでした。NURBS 曲面は、ノットベクトル \(U\)、\(V\) の各区間(スパン)ごとに異なる多項式で表されるので、まずスパンごとの有理ベジエパッチ(セグメント)に分解します。
ノット挿入(Boehm の方法)
次数 \(p\) の B スプラインにノット \(\bar u\in[u_k,u_{k+1})\) を 1 つ挿入すると、形を変えずに制御点が 1 つ増えます。新しい制御点は
$$ \mathbf Q_i=\alpha_i\,\mathbf H_i+(1-\alpha_i)\,\mathbf H_{i-1},\qquad \alpha_i=\begin{cases}1 & i\le k-p\\[2pt] \dfrac{\bar u-u_i}{u_{i+p}-u_i} & k-p+1\le i\le k\\[6pt] 0 & i\ge k+1\end{cases} $$です。NURBS では、この線形補間を同次座標 \(\mathbf H_i=(w_ix_i,w_iy_i,w_iz_i,w_i)\) で行うのが正しいやり方です(デカルト座標で補間すると形が変わります)。コードの insertKnotH がこれに当たります。
各内部ノットの多重度が次数 \(p\) に達するまで挿入を繰り返すと、隣り合うスパンの間で曲線が「切り離され」、各スパンの \(p+1\) 個の制御点がそのままベジエ曲線の制御点になります(refineToBezierH)。曲面では u 方向・v 方向に順に行います(surfaceToBezier)。
各セグメントのパラメータは局所的に \(0\le u,v\le1\) に取り直されています。元の NURBS のパラメータ \(\hat u\) との関係は、スパン \([\hat u_a,\hat u_{a+1}]\) について \(\hat u=\hat u_a+u(\hat u_{a+1}-\hat u_a)\) です。画面の「ボトル本体 [3,3]」は、u 方向 3 番目・v 方向 3 番目のスパンのセグメントという意味です。
| 例:ペットボトル | u 方向(円, 2 次) | v 方向(母線) | セグメント数 |
|---|---|---|---|
| 本体 | 4 スパン | 12 点・3 次 → 9 スパン | 36 |
| キャップ側面 | 4 | 2 点・1 次 → 1 | 4 |
| キャップ天面 | 4 | 3 点・2 次 → 1 | 4 |
| リング | 4 | 4 点・1 次 → 3 | 12 |
| 合計 | 56(画面のスライダー #1〜#56) | ||
12各セグメントの包含箱(境界箱)
セグメント \(k\)(パッチ 1 枚)の包含箱は、第 10 節の計算をそのパッチだけに対して行えば得られます。法線を \(\mathbf e_x,\mathbf e_y,\mathbf e_z\) とした 3 回の極値計算で
$$ \mathrm{Box}_k=[x^{k}_{\min},x^{k}_{\max}]\times[y^{k}_{\min},y^{k}_{\max}]\times[z^{k}_{\min},z^{k}_{\max}], \qquad x^{k}_{\min}=\min_{(u,v)\in[0,1]^2}\frac{X_k(u,v)}{W_k(u,v)}\ \ \text{など}. $$各値は、そのパッチの ①4 隅、②4 本の境界の有理曲線の停留点(\(X'W-W'X=0\) の根)、③内部停留点(\(X_uW-W_uX=0,\ X_vW-W_vX=0\))のいずれかで生じます。画面ではこの 6 つの接点を ● で示し、種別(隅・境界・内部)も表示しています(コードの segBBox)。
全体の境界箱との関係
曲面全体はセグメントの和集合なので、全体の最小・最大は各セグメントの最小・最大の最小・最大です。
$$ x_{\min}=\min_k x^{k}_{\min},\qquad x_{\max}=\max_k x^{k}_{\max}\qquad(y,z\ \text{も同様}) $$つまり全体の境界箱 = セグメントの包含箱の和集合を囲む箱です。プログラムでは全体の箱(getSurfaceBBox)とセグメントの箱(segBBox)を別々に計算していますが、9 種類のサンプルすべてで両者が完全に一致することを確認しています。
セグメントの制御点の箱との比較
セグメントの制御点 \(\mathbf P_{ij}=\mathbf H_{ij}/w_{ij}\) の最小・最大で作る箱も、凸包性(第 7 節)により必ずセグメントを含みます。ただし次の点で厳密な包含箱より大きくなります。
- 端点以外の制御点は一般に曲面上にない(曲面から離れた位置にある)。
- 重みが大きい制御点ほど曲面がそこへ引き寄せられるが、届くわけではない。
画面の情報欄に出る「体積比(曲面/制御点)」はこの差を表します。分解直後のセグメントは小さく平坦に近いので、全体の場合より比が 100% に近くなる傾向があります。
13平面との交差(断面線)の計算
平面 \(\Pi:\mathbf n\cdot\mathbf p=d\) と曲面の交線は、符号付き距離 \(f(u,v)=A/W-d\) の零点集合 \(\{(u,v)\mid f(u,v)=0\}\) です。\(W>0\) なので、これは多項式の零点集合
$$ F(u,v)=A(u,v)-d\,W(u,v)=\sum_{i,j}w_{ij}\,(\mathbf n\cdot\mathbf P_{ij}-d)\,B_i(u)B_j(v)=0 $$と同じです。\(F\) の制御値 \(w_{ij}(\mathbf n\cdot\mathbf P_{ij}-d)\) の符号は、\(w_{ij}>0\) なので制御点 \(\mathbf P_{ij}\) が平面のどちら側にあるかだけで決まります。
① 凸包による棄却(枝刈り)
第 7 節と同じ理由で、\(f\) の値はパッチの制御点の符号付き距離 \(\mathbf n\cdot\mathbf P_{ij}-d\) の凸結合です。したがって
この判定は制御点の数だけの内積で済み、非常に安価です(コードの sliceRecurse の冒頭)。有理パッチでも、デカルト座標の制御点 \(\mathbf H_{ij}/w_{ij}\) を使えば非有理と全く同じ判定になります。
② 同次座標での 4 分割再帰
制御点が平面の両側にまたがる場合は、交差の可能性があるので、de Casteljau 法で \(u=\tfrac12,\ v=\tfrac12\) の 4 つの部分パッチに分け、それぞれに ① から適用します。分割は同次座標で線形補間します(コードの bisect1D, split4)。部分パッチは元の曲面の一部を正確に表す有理ベジエパッチなので、① の判定は何度でも正しく使えます。
③ 末端での線分抽出
最大深さ \(D\)(画面の「再帰分割深度」、計測平面の交線は \(\max(D,6)\))に達した部分パッチでは、4 隅の点を使います。ベジエパッチの 4 隅の制御点は曲面上の点そのものなので、隅の符号付き距離 \(f_0,\dots,f_3\) は正確な値です。
- 4 隅 \(\mathbf A,\mathbf B,\mathbf C,\mathbf D\) の \(f\) を計算する。
- 符号が変わる辺(例:\(f_A>0>f_B\))では、線形補間で交点を求める。 $$ \mathbf p=\mathbf A+\lambda(\mathbf B-\mathbf A),\qquad \lambda=\frac{f_A}{f_A-f_B} $$
- 得られた交点を 2 つずつ結んで、断面の線分にする。
末端パッチは十分小さく平坦に近いので、隅を結ぶ直線で曲面を近似した誤差は小さくなります。深さを 1 増やすと末端の大きさは半分になり、近似誤差はおよそ \(1/4\) になります(曲面の 2 次の項による誤差)。
- 交点は「隅を結ぶ線分上」で求めているので、厳密には曲面上の点ではありません(深さとともに急速に収束)。
- 末端パッチの内側だけで閉じる小さなループ(4 隅の符号がすべて同じなのに中で交差する場合)は、線分として出力されません。① の判定で棄却はされないので、深さを上げれば検出できます。
- 4 辺すべてで符号が変わる鞍点状のセルでは、交点の結び方に 2 通りの可能性があります。末端パッチを十分小さくすればこの曖昧さの影響は小さくなります。
④ 距離の極値・包含箱との関係
第 3 節の結果を使うと、交差の有無は厳密に判定できます。
$$ \text{パッチ } k \text{ が平面と交差する}\iff g^k_{\min}\le d\le g^k_{\max}\qquad(g^k=\mathbf n\cdot\mathbf S_k) $$| 判定法 | 使う量 | 性質 |
|---|---|---|
制御点の凸包(① / sliceRecurse) | \(\min,\max\ \mathbf n\cdot\mathbf P_{ij}\) | 安価・保守的(交差しないのに「交差しうる」と判定することがある)。分割すると鋭くなる。 |
| セグメントの包含箱(第 12 節) | 箱の 8 頂点の \(\mathbf n\cdot\mathbf p\) の最小・最大 | 座標軸方向には厳密。斜めの平面では箱の角の分だけ保守的。 |
| 法線方向の厳密な極値(第 3〜8 節) | \(g^k_{\min},g^k_{\max}\) | 厳密(必要十分)。法線ごとに 1 回の計算で、平行な平面すべてに使える。 |
プログラムでは、平行な複数の断面(等距離線)の配置範囲を \([g_{\min},g_{\max}]\) から決めています( 平面と NURBS 曲面の交差、および境界箱の計算について、従来から使われている方法と今回の方法を比べます。黄色の行が nurbs_plane_distance.html で使っている方法です。 この節は、計算幾何・CAD 分野の標準的な文献知識に基づく整理です。個々の文献を新たに検索・確認したものではありません。 平面は陰関数 \(\mathbf n\cdot\mathbf p=d\) で書けるので、曲面の式を代入すると交線はパラメータ空間の多項式 の零点集合になります(第 13 節)。一般の曲面どうしの交差より扱いやすく、多くの手法はこの性質を使います。 曲線の場合、厳密な方法は広く使われています。たとえば 3 次ベジエでは、導関数の 2 次方程式を解いて箱を求め、フォントや SVG の描画に利用されています。曲面の場合、実用システムでは制御点の箱か細分割した箱が主流です。光線追跡や干渉判定の階層構造(BVH)では、「必ず包含する」ことが重要で、少し大きい分には困らないからです。buildPlanes)。\(g_{\min}14アルゴリズム全体とコードの対応
// 法線 n に対する g = A/W の全体最小・最大(d に依存しない)
for 各有理ベジエパッチ:
a_ij = w_ij (n·P_ij), w_ij // A, W の係数
① 4 隅で g を評価
② 4 本の境界曲線: bz1Roots( a'w − w'a ) の根で g を評価
③ 内部: rec(部分パッチ)
枝刈り: [min a/w, max a/w] が最小・最大を更新しえない → return
Gu = A_u W − W_u A, Gv = A_v W − W_v A // nProd で係数計算
u 範囲 ← Gu の凸包 ∩ {y=0}, v 範囲 ← Gv の凸包 ∩ {y=0}
空 → return
十分小さい → ニュートン法で精密化し g を評価
十分縮んだ → その範囲の部分パッチで rec
それ以外 → 4 分割して rec
// 平面 (n, d) に対する結果
最遠 = max(gmax − d, d − gmin)
最短 = gmin − d / d − gmax / 0(交差:交線を描く)
// 境界箱
n = (1,0,0), (0,1,0), (0,0,1) で同じ計算 → [xmin,xmax]×[ymin,ymax]×[zmin,zmax]
// セグメント k の包含箱
上の ①②③ をパッチ k だけに対して 3 方向で実行 → Box_k (全体の箱 = Box_k の和集合)
// 平面 (n, d) との交差(断面線)
sliceRecurse(パッチ, 深さ):
制御点 P_ij がすべて平面の片側 → return // 凸包による棄却
深さ = D → 4 隅の f の符号変化から線形補間で線分を出力
それ以外 → 同次座標で 4 分割して sliceRecurse
処理 nurbs_plane_distance.html の関数 関数名の S/s は u 方向、T/t は v 方向を表します(例:nDs = ∂/∂u、nClipT = v 方向のクリップ)。NURBS → 有理ベジエパッチ(ノット挿入) insertKnotH, refineToBezierH, surfaceToBezier, buildBezierPatchesベジエ関数の微分・積・部分区間 bz1Deriv, bz1Prod, nDs, nDt, nProd, nSubS, nSubT停留点方程式 \(G_u,G_v\) stationaryNets凸包と y=0 の交差(クリッピング) hullZeroRange, nClipS, nClipT1 変数の根(境界曲線) bz1Roots1 パッチの極値探索(①②③) ratPatchExtrema全パッチ・法線ごとのキャッシュ ratExtremaAlong, getExtrema最短・最遠の場合分け planeDistanceResult平面との交差:凸包棄却・4 分割再帰 sliceRecurse, split4, bisect1D平面との交差:末端の線分抽出 leafExtract断面の配置(等距離線)・計測平面の交線 buildPlanes, getCrossSegs曲面全体の境界箱 getSurfaceBBoxセグメントごとの包含箱 segBBox, setSegIdx15従来法との比較
15.1 平面と NURBS 曲面の交差
手法 概要 長所 短所 出力 ポリゴン化(テッセレーション) 曲面を三角形メッシュにしてから三角形と平面を交差させる。3D プリンタのスライサー(STL の断面化)がこの方式。 単純で高速。実装が容易。 精度がメッシュの細かさで決まる。曲面上の厳密な点ではない。 線分 格子評価+マーチングスクエア \((u,v)\) の一様格子で \(F\) の符号を調べ、セル内を補間する。 実装が容易。点列化もしやすい。 格子より小さいループを見落とす。平坦部でも同じ密度で計算するので無駄が多い。 線分 再帰分割法(Lane–Riesenfeld 型) 制御点の凸包(または箱)で交差しない部分を棄却し、残りを de Casteljau で分割する。 交差しない部分を確実に除外でき、交線近くだけ細かくなる(適応的)。 末端は近似。出力は順序のない線分。 順序のない線分 追跡法(マーチング) 開始点から接線方向 \(\mathbf t=\mathbf n\times\mathbf N(u,v)\)(\(\mathbf N\) は曲面法線)に予測し、ニュートン法で \(F=0\) に補正して進む。 高精度で、順序付きの点列が得られる。刻み幅を曲率に合わせられる。 すべての開始点(特に閉ループ)を見つける保証が要る。平面と曲面が接する特異点で不安定。 順序付き点列 ハイブリッド法 分割法で開始点とループを確実に検出し、追跡法で精密な点列を作る。 確実さと精度を両立。CAD カーネルの標準的な構成。 実装が複雑。 順序付き点列 ループ検出条件 パッチ内に閉ループがないことを法線の向き(ガウス写像)の範囲から保証し、分割をそこで止める(Sederberg ら、Hohmeyer など)。 見落としがないことを理論的に保証。 判定自体の計算が必要。 (判定のみ) 区間演算・ベジエクリッピング \(F\) の Bernstein 係数の凸包と 0 の交差で \((u,v)\) 範囲を絞り込む。ベジエクリッピングは Nishita, Sederberg, Kakimoto(SIGGRAPH 1990)で、有理パッチの光線追跡として導入された。 収束が速い(根の近くで 2 次収束)。 交線は「曲線」なので、点の解法をそのまま適用するには工夫が要る。 点(範囲) sliceRecurse)は再帰分割法そのものです。棄却には制御点の凸包の符号判定(第 13 節 ①)を使い、末端では 4 隅を補間しています(格子評価と同じ考え方)。ベジエクリッピングは距離の極値計算(第 8 節)で使っており、交差の計算には使っていません。改善の方向としては、次のものがあります。
15.2 境界箱(包含箱)
手法 概要 包含の保証 箱の大きさ 計算量 制御点の箱 凸包性から制御点の最小・最大を使う(画面の点線)。 あり(\(w>0\)) 大きすぎることが多い 最小 サンプリング/メッシュの箱 評価した点の最小・最大を使う。 なし(極値を外すと小さすぎる) ほぼ正確 点数に比例 細分割した制御点の箱 分割を繰り返して制御点の箱を縮める。 あり 分割するほど正確(誤差は分割幅の 2 乗に比例) 中 区間演算・アフィン演算 式を区間で評価する。 あり 依存関係があると大きくなる 中 極値の厳密計算 導関数 = 0(有理なら \(X'W-W'X=0\))を解き、隅・境界・内部の候補を調べる(第 4〜10、12 節)。 あり(数値誤差の範囲で) 正確(最小の箱) やや大きい(今回のサンプルでは 1 方向あたり数 ms)
15.3 用途ごとの使い分け
用途 推奨する方法 理由 粗い絞り込み(干渉判定の前段など) 制御点の箱 最も安価で、必ず包含する。 厳密な寸法、最短・最遠距離 極値の厳密計算(今回の方法) 最小の箱と接点が得られる。 断面の点列 分割法と追跡法のハイブリッド 見落としがなく、高精度な順序付き点列が得られる。