平面から NURBS 曲面への最短・最遠距離と境界箱
— 有理関数の極値をベジエクリッピングで厳密に求める —
nurbs_plane_distance.html で使っている計算方法を、式の導出から順に説明します。ポイントは「有理関数のまま扱わず、分子だけの多項式（ベジエ関数）の問題に直す」ことです。
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}] $$が得られます。
11アルゴリズム全体とコードの対応
// 法線 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]
| 処理 | nurbs_plane_distance.html の関数 |
|---|---|
関数名の S/s は u 方向、T/t は v 方向を表します（例：nDs = ∂/∂u、nClipT = v 方向のクリップ）。 | |
| NURBS → 有理ベジエパッチ | 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 |
| 交線（最短距離 0 のとき） | sliceRecurse（元のスライサー） |
| 境界箱 | getSurfaceBBox |