平面と 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])。したがって

符号付き距離 $$ f(u,v)=g(u,v)-d,\qquad g(u,v)=\frac{A(u,v)}{W(u,v)} $$ \(g\) は「点を法線方向へ射影した高さ」で、平面の位置 \(d\) を含みません。
平面 Π : n·p = d n 曲面 S(u,v) 最短 = gmin − d 最遠 = gmax − d
図1 平面が曲面の下にある場合。最短点は法線方向の高さ g が最小の点、最遠点は最大の点。

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\) 上の全点)
重要な帰結 問題は「有理関数 \(g=A/W\) の、曲面全体での最小値・最大値とその点を求める」ことに帰着します。\(g\) は \(d\) を含まないので、法線 \(\mathbf n\) を決めて一度計算すれば、平面を平行移動しても再計算は不要です(コードは法線ごとに結果をキャッシュしています)。交差する場合の交線は、元のスライサーの断面抽出で描いています。

4極値の候補:隅・境界・内部

パラメータ領域は正方形 \([0,1]^2\) です。閉じた領域上の滑らかな関数の最小・最大は、次のどれかで生じます。

(0,0)(1,1) ① 4 隅:値をそのまま評価 ② 境界の 4 本の有理曲線上の停留点 a′w − w′a = 0 の根(1 変数) ③ 内部の停留点 ∂g/∂u = ∂g/∂v = 0(2 変数の連立)
図2 パラメータ領域 [0,1]² における極値の候補。境界上の極値は境界曲線を 1 変数関数として調べれば足りる。

境界上の極値は、内部の条件 \(\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 になりません。したがって、停留点の条件は分子だけで書けます。

停留点方程式(多項式) $$ G_u(u,v)=A_u W-W_u A=0,\qquad G_v(u,v)=A_v W-W_v A=0 $$ \(G_u, G_v\) は有理関数ではなく多項式なので、ベジエ関数(Bernstein 多項式)として係数を正確に表せます。

次数は次のようになります(\(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_{i,j}\frac{a_{ij}}{w_{ij}}\ \le\ g(u,v)\ \le\ \max_{i,j}\frac{a_{ij}}{w_{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\) についてもこの点群の凸包の内側にあります。したがって

u 方向のクリップ範囲 凸包と直線 \(y=0\) が交わる区間 \([u_{lo},u_{hi}]\) の外では、\(G_u\neq0\) なので停留点は存在しない。

v 方向は \(G_v\) について同じことを行います。どちらかが交わらなければ、その部分パッチには停留点がありません。

再帰の手順

  1. 枝刈り:制御値の比 \(a/w\) の範囲が現在の最小・最大を更新しえなければ終了(第 7 節)。
  2. クリッピング:\(G_u\) から u の範囲、\(G_v\) から v の範囲を求める。どちらかが空なら終了。
  3. 縮小できた場合(面積比 < 0.64):その範囲の部分パッチを同次座標で切り出し、手順 1 から繰り返す。
  4. 縮小できない場合(解が複数ありうる):4 分割して、それぞれに手順 1 から適用する。
  5. 十分小さくなった場合(幅 < 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\) です。

0 t=1 1 1/3 −1/3 −1 t = 1/2
図3 G(t) の係数(橙)とその凸包。4 点が一直線 y = 1 − 2x 上にあるため凸包は線分(太い橙線)につぶれ、G(t) 自体もこの直線になる(3 次で表現しているが実際は 1 次。一般には第 5 節の打ち消しで 2 次になり、この例では左右対称性によりさらに 1 次まで下がる)。y = 0 との交点 t = 1/2 が 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}] $$

が得られます。

制御点の箱との関係 重みが正なら第 7 節の凸包性により、曲面は必ずデカルト座標の制御点 \(\mathbf P_{ij}\) の最小・最大で作る箱(画面の点線)の中に入ります。ただしこの範囲は一般に広すぎます(第 9 節の例では、45° 方向の高さの上限が実際の 1 に対して \(\sqrt2\))。厳密な境界箱は、曲面が実際に箱の 6 面に接する点(画面の ●)を求めて作るので、余分な隙間がありません。

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)。

u 方向のノットベクトル(例:円 U = [0,0,0, ¼,¼, ½,½, ¾,¾, 1,1,1]) 0 ¼ ½ ¾ 1 セグメント 1セグメント 2セグメント 3セグメント 4 局所 u∈[0,1]局所 u∈[0,1] 局所 u∈[0,1]局所 u∈[0,1] 内部ノット(橙)の多重度を p にすると、各スパンが独立した有理ベジエセグメントになる。
図4 ノット挿入によるスパンごとのセグメント分解(u 方向)。v 方向も同様に分け、u のスパン数 × v のスパン数 枚のパッチになる。

各セグメントのパラメータは局所的に \(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
キャップ側面42 点・1 次 → 14
キャップ天面43 点・2 次 → 14
リング44 点・1 次 → 312
合計56(画面のスライダー #1〜#56)
分解の利点 セグメントは「1 本の多項式(同次座標)」なので、第 1〜10 節の議論(凸包性、ベジエ関数の微分・積、de Casteljau 分割)がそのまま使えます。また、スパンの境目で分割しておくと、パッチごとの形が単純になり、凸包による判定が鋭くなります。

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% に近くなる傾向があります。

包含箱の使い道 セグメントの包含箱は、そのセグメントが空間のどこにあるかを示す最小の直方体です。光線・平面・他の物体との干渉判定で、「箱と交わらないセグメントは調べなくてよい」という前段の絞り込み(粗い判定)に使えます。平面との交差では、次節の \([g^k_{\min},g^k_{\max}]\) がさらに直接的な判定になります。

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\) の凸結合です。したがって

交差しない十分条件 すべての制御点が平面の同じ側にある(全部 \(>0\) または全部 \(<0\))なら、そのパッチ(部分パッチ)は平面と交差しない。

この判定は制御点の数だけの内積で済み、非常に安価です(コードの sliceRecurse の冒頭)。有理パッチでも、デカルト座標の制御点 \(\mathbf H_{ij}/w_{ij}\) を使えば非有理と全く同じ判定になります。

② 同次座標での 4 分割再帰

制御点が平面の両側にまたがる場合は、交差の可能性があるので、de Casteljau 法で \(u=\tfrac12,\ v=\tfrac12\) の 4 つの部分パッチに分け、それぞれに ① から適用します。分割は同次座標で線形補間します(コードの bisect1D, split4)。部分パッチは元の曲面の一部を正確に表す有理ベジエパッチなので、① の判定は何度でも正しく使えます。

u=0u=1 ■ 棄却:制御点が平面の片側に収まる ■ 末端:最大深さまで分割された部分パッチ — 交線 f(u,v)=0 交線から離れた領域は浅い段階で大きく棄却され、 交線の近くだけが細かく分割される(適応的)。 末端の数 ≈ 交線の長さ ÷ 末端の大きさ
図5 パラメータ空間での再帰分割の様子(模式図)。棄却されたパッチの内側は二度と調べない。

③ 末端での線分抽出

最大深さ \(D\)(画面の「再帰分割深度」、計測平面の交線は \(\max(D,6)\))に達した部分パッチでは、4 隅の点を使います。ベジエパッチの 4 隅の制御点は曲面上の点そのものなので、隅の符号付き距離 \(f_0,\dots,f_3\) は正確な値です。

  1. 4 隅 \(\mathbf A,\mathbf B,\mathbf C,\mathbf D\) の \(f\) を計算する。
  2. 符号が変わる辺(例:\(f_A>0>f_B\))では、線形補間で交点を求める。 $$ \mathbf p=\mathbf A+\lambda(\mathbf B-\mathbf A),\qquad \lambda=\frac{f_A}{f_A-f_B} $$
  3. 得られた交点を 2 つずつ結んで、断面の線分にする。

末端パッチは十分小さく平坦に近いので、隅を結ぶ直線で曲面を近似した誤差は小さくなります。深さを 1 増やすと末端の大きさは半分になり、近似誤差はおよそ \(1/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}]\) から決めています(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, nClipT
1 変数の根(境界曲線)bz1Roots
1 パッチの極値探索(①②③)ratPatchExtrema
全パッチ・法線ごとのキャッシュratExtremaAlong, getExtrema
最短・最遠の場合分けplaneDistanceResult
平面との交差:凸包棄却・4 分割再帰sliceRecurse, split4, bisect1D
平面との交差:末端の線分抽出leafExtract
断面の配置(等距離線)・計測平面の交線buildPlanes, getCrossSegs
曲面全体の境界箱getSurfaceBBox
セグメントごとの包含箱segBBox, setSegIdx
検証結果 同梱の 9 種類のサンプル形状(ペットボトル、ワイングラス、マグカップ、電球、スプーン、主翼と胴体、トーラス、球、重み付き波形)について、それぞれ 11 方向で計算した最小・最大を、各パッチ 61×61 点の密なサンプリングと比較しました。差は 10⁻¹⁵ 程度で、取りこぼしはありませんでした。中央の重みが大きい有理パッチでは、内部の停留点(③)による極大も正しく検出しています。また、全サンプルで「セグメントの包含箱の和集合」と「曲面全体の境界箱」が一致することを確認しています。

15従来法との比較

平面と NURBS 曲面の交差、および境界箱の計算について、従来から使われている方法と今回の方法を比べます。黄色の行が nurbs_plane_distance.html で使っている方法です。

この節は、計算幾何・CAD 分野の標準的な文献知識に基づく整理です。個々の文献を新たに検索・確認したものではありません。

15.1 平面と NURBS 曲面の交差

平面は陰関数 \(\mathbf n\cdot\mathbf p=d\) で書けるので、曲面の式を代入すると交線はパラメータ空間の多項式

$$ F(u,v)=\sum_{i,j}w_{ij}(\mathbf n\cdot\mathbf P_{ij}-d)\,B_i(u)B_j(v)=0 $$

の零点集合になります(第 13 節)。一般の曲面どうしの交差より扱いやすく、多くの手法はこの性質を使います。

手法概要長所短所出力
ポリゴン化(テッセレーション)曲面を三角形メッシュにしてから三角形と平面を交差させる。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)

曲線の場合、厳密な方法は広く使われています。たとえば 3 次ベジエでは、導関数の 2 次方程式を解いて箱を求め、フォントや SVG の描画に利用されています。曲面の場合、実用システムでは制御点の箱か細分割した箱が主流です。光線追跡や干渉判定の階層構造(BVH)では、「必ず包含する」ことが重要で、少し大きい分には困らないからです。

今回の方法の特徴(境界箱) 今回の方法は、法線方向 \(\mathbf n\) の支持関数(凸解析の用語) $$ h(\mathbf n)=\max_{(u,v)}\ \mathbf n\cdot\mathbf S(u,v) $$ を厳密に求めるものです。条件 \(G_u=G_v=0\) は「曲面の法線が \(\mathbf n\) と平行になる点」を求めることと同じで、方向 \(\mathbf n\) から見たシルエットの極点を求める古典的な問題にあたります。このため、軸平行の箱だけでなく、任意の向きの箱や k-DOP(複数方向の範囲で囲む凸多面体)も同じ処理で作れます。

15.3 用途ごとの使い分け

用途推奨する方法理由
粗い絞り込み(干渉判定の前段など)制御点の箱最も安価で、必ず包含する。
厳密な寸法、最短・最遠距離極値の厳密計算(今回の方法)最小の箱と接点が得られる。
断面の点列分割法と追跡法のハイブリッド見落としがなく、高精度な順序付き点列が得られる。