アルゴリズム解説

有理ベジエ曲線・曲面の最近点／交点探索
— 厳密法と近似法の具体的な仕組み

本稿は、有理ベジエ曲線と有理ベジエ曲面の最近点および交点を求めるツールで実装した 2つの区間絞り込み方式 —「厳密法（多項式係数＋凸包によるベジエクリッピング）」と 「近似法（サンプリング）」— の内部処理を、実装した関数名とともに具体的に解説する。

1問題設定

$n$次の有理ベジエ曲線 $C(t)$（制御点 $P_0,\dots,P_n$、重み $w_0,\dots,w_n$）と、 $u,v$ 方向とも $n$次の有理ベジエ曲面 $S(u,v)$（制御点 $P_{ij}$、重み $w_{ij}$）が与えられたとき、 次の2つを求めたい。

  • 最近点：$\lVert C(t)-S(u,v)\rVert$ を最小にするパラメータ $(t,u,v)$
  • 交点：$C(t)=S(u,v)$ となるすべてのパラメータの組（複数存在しうる）

どちらも、点 $Q=S(u,v)$ に対する曲線上の最近点条件（点から曲線への垂線条件）と、 点 $P=C(t)$ に対する曲面上の最近点条件を同時に満たす停留点探索に帰着する。 これは次の3本の方程式を $(t,u,v)$ について解くことと等価である。

$$ g(t)=C'(t)\cdot\bigl(C(t)-S(u,v)\bigr)=0,\qquad q_u(u,v)=S_u(u,v)\cdot\bigl(S(u,v)-C(t)\bigr)=0,\qquad q_v(u,v)=S_v(u,v)\cdot\bigl(S(u,v)-C(t)\bigr)=0 $$

交点では $C(t)-S(u,v)=\mathbf{0}$ となるため、この3条件は自動的に成り立つ（内積がゼロベクトルとの積になる）。 つまり交点は必ずこの停留点探索でも「解」として現れる、という点が両方式に共通する土台になっている。

2共通の反復構造

3変数を同時に解くのは煩雑なため、両方式とも次のように t の絞り込み と u, v の絞り込み を交互に繰り返す（ガウス・ザイデル型の反復）。

曲面点 Q=S(u,v) を固定して g(t) を解き [t_min, t_max] を絞る 曲線点 P=C(t) を固定して q_u=0 を解き [u_min, u_max] を絞る 新しい u を用い q_v=0 を解いて [v_min, v_max] を絞る P = C(t) Q = S(u,v) を更新して繰り返す
図1：t → u → v の順に片方ずつ固定しながら区間を絞り込む反復ループ。区間の絞り込み方法だけが厳密法／近似法で異なる。

実装では computeClosestExact() と computeClosestSampling() が この反復ループそのものであり、内部の「絞り込み1回分」の処理だけが2方式で入れ替わっている。

3近似法（サンプリング）approx

実装：computeClosestSampling() / gFuncGen() / quFuncGen()

近似法は、現在の区間 $[t_{\min},t_{\max}]$ を等間隔に $N$ 点（既定 $N=9$）サンプリングし、 各点で $g(t)$ を de Casteljau 法による通常の点・接線評価（積の法則）で数値的に求める。 最も絶対値が小さいサンプル点の両隣を新しい区間として採用する、という単純な縮小を行う。 

// 曲線側 1 反復ぶんの絞り込み（擬似コード）
ts = linspace(t_min, t_max, 9)
values = ts.map(t => gFuncGen(t, Q, curveNet))
i* = argmin |values[i]|
t_min, t_max = ts[i*-1], ts[i*+1]   // 両隣に縮小
t = ts[i*]

曲面側も同様に、現在の矩形 $[u_{\min},u_{\max}]\times[v_{\min},v_{\max}]$ を $5\times5$ グリッドでサンプリングし、 $q_u^2+q_v^2$ が最小のセルの周囲へ矩形を縮小する。

毎反復、区間はサンプル間隔の分だけ縮む（曲線側で約 $1/4$、曲面側で約 $1/2$）ため、 反復を重ねれば幾何級数的に収束する。ただし次の弱点がある。

  • サンプリング間隔より細かい符号変化（根の存在）を 見落としうる（エイリアシング）。
  • 区間内に複数の根がある場合、どちらに収束するかはサンプル位置に依存し制御できない。
  • 根が存在しないことを保証する手段がない（単に「これ以上絞れなかった」としか分からない）。

計算コストは低く実装も単純だが、上記の理由から探索の「健全性」（根を取りこぼさないという保証）が無い。

4厳密法（ベジエクリッピング）exact

実装：computeClosestExact() / buildDistNumer() / hullZeroBracket() / getSubCP()

厳密法は、$g(t)$ や $q_u(u,v),q_v(u,v)$ をサンプリングではなく厳密なベルンシュタイン多項式係数として構成し、 その制御多角形の凸包を使って「根を絶対に取りこぼさない」区間絞り込みを行う。 これが原論文でいう Bezier Clipping 法の核心部分である。

4.1　有理ベジエの基本演算

重み付き制御点 $Pw_i=w_iP_i$（ベクトル、次数 $n$）と重み列 $Wt_i=w_i$（スカラー、次数 $n$）の 2本のベルンシュタイン多項式として曲線を表す（同次座標表現）。

$$ C(t)=\frac{N(t)}{D(t)},\quad N(t)=\text{deCasteljau}(Pw,t),\quad D(t)=\text{deCasteljau}(Wt,t) $$

微分もベルンシュタイン基底のまま厳密に扱える。次数 $d$ の制御点列 $c_0,\dots,c_d$ の導関数は、 次数 $d-1$ の制御点列 $d(c_{i+1}-c_i)$ で 厳密に 表される（数値微分ではない）。 これと商の微分公式を組み合わせて $C'(t)$ を得る（bezDerivCP() / evalRationalDeriv()）。

4.2　距離関数を厳密な多項式として構成する

$C(t)=N/D$ を代入して整理すると、

$$ g(t)=\frac{\bigl(N'D-ND'\bigr)\cdot\bigl(N-QD\bigr)}{D(t)^3} $$

重みはすべて正なので $D(t)^3>0$ が常に成り立ち、符号は分子だけで決まる。そこで分子

$$ \tilde g(t) = \bigl(N'D-ND'\bigr)\cdot\bigl(N-QD\bigr) $$

を厳密なベルンシュタイン多項式として構成する。次数は $(2n-1)+n=3n-1$ で、原論文の記述と一致する。 構成には次の「ベジエ関数の積の公式」（畳み込み）を使う。次数 $p,q$ の多項式 $f,g$ の積は次数 $p+q$ の多項式になり、 その制御点は

$$ h_k=\sum_{i=\max(0,k-q)}^{\min(p,k)} \frac{\binom{p}{i}\binom{q}{k-i}}{\binom{p+q}{k}}\, f_i\, g_{k-i} $$

で厳密に求まる（bezMul()）。これをベクトル×スカラーの3成分に適用し（bezMulVecScalar()）、 差を取り（bezSubArr()）、最後に内積として3成分の積を足し合わせる（bezDotVec()）ことで $\tilde g(t)$ の係数列を一切近似なしに得る（buildDistNumer()）。

4.3　凸包によるゼロ交差区間の抽出

ベジエ曲線の重要な性質として、$\tilde g(t)\ (t\in[0,1])$ のグラフは、制御点集合 $\{(i/(3n{-}1),\ \tilde g_i)\}$ の凸包の内部に必ず収まる。したがって、この凸包が ゼロ水平線 $y=0$ と交差する $t$ の範囲を求めれば、それは区間内に存在しうるすべての根を 取りこぼしなく包含する。これが「クリッピング」の名の由来であり、 近似法の「サンプル点の最小値の周辺を取る」という発見的な絞り込みとは異なり、数学的に根の存在範囲を保証する。

y=0 t_lo t_hi 根はこの帯の中に必ず存在する
図2：制御多角形（青い折れ線）の凸包（薄い塗り）が y=0 と交差する区間 [t_lo, t_hi]（琥珀色の帯）が、根の存在を保証する新しい絞り込み区間になる。

実装（hullZeroBracket()）では、モノトーンチェーン法で凸包を求め、 その各辺について符号が変わる箇所の $y=0$ との交点を計算し、最小値と最大値を新しい $[t_{\min},t_{\max}]$ とする。 凸包が一度もゼロと交差しない場合は、その区間に根は存在しないと確定的に判定できる （この場合は絶対値最小の制御点近傍を暫定的な代替値として使う）。

4.4　部分区間への正確な再分割

次の反復では、狭まった $[t_{\min},t_{\max}]$ 上での $\tilde g$ を改めて構成し直す必要がある。 これは de Casteljau分割を2回 適用することで、サンプリングではなく厳密に、 元の制御点数のまま部分区間の制御点列を得られる（getSubCP() / subdivide()）。 実装では毎回、直前の縮小結果からではなく元の（最初の）制御点列から改めて分割し直しており、 反復を重ねても丸め誤差が蓄積しない設計にしている。

4.5　曲面への拡張：u成分とv成分を逐次抽出

2変数 $(u,v)$ を同時にクリッピングするのは複雑になるため、曲面が $u,v$ 方向にテンソル積構造を持つことを利用する。 $v$ を現在値に固定して制御網を $v$ 方向に de Casteljau で簡約すると、$S(u, v_{\text{fixed}})$ の 厳密な $u$ 方向 1変数ベルンシュタイン多項式（reduceAtV()）が得られる。 これは曲線の場合とまったく同じ4.2〜4.4の手続き（buildDistNumer() → hullZeroBracket() → getSubCP()）を そのまま適用できることを意味する。$u$ が求まったら、今度は $u$ 方向に簡約（reduceAtU()）して 同じ手続きで $v$ を絞り込む。これは原論文が式(4)から $u$成分の区間を、式(5)から $v$成分の区間を それぞれ抽出すると説明している構造と対応する。

5頑健性のための工夫

実装：coarseGlobalSeedMulti() / refineSurfaceOnly() / localPolish() / computeAllIntersections()

4章の手続きだけでは実用上いくつかの落とし穴があることが、総当たり探索との比較検証で判明した。 それぞれに対処を加えている。

5.1　複数の停留点問題 → マルチスタート

交互反復は「垂直条件を満たす停留点」には収束するが、その停留点が必ずしも大域最小（または真の交点）とは限らない。 近傍に複数の停留点が存在しうるためである。対策として、粗いグローバルグリッドで距離が小さい点を複数 （既定14点、パラメータ空間で互いに十分離れたもの）拾い、それぞれを独立の初期値として反復を実行し、 最終的に最良の結果を採用する。

5.2　境界極値問題 → 端点の明示チェック

曲線の端点（$t=0$ または $t=1$）で距離が最小になる場合、その点は $g(t)=0$ を満たすとは限らない （区間内部の停留条件ではなく境界での最小だから）。根探索だけでは原理的に見つからない。 そこで $t$ を $0,1$ に固定し、曲面側だけを同じ方式で絞り込む処理（refineSurfaceOnly()）を 別途実行し、通常の反復結果より優れていれば採用する。

5.3　局所的なずれ → 直接パターンサーチによる磨き上げ

5.1・5.2の対策を入れても、交互反復方式そのものの収束特性により、真の最小点のごく近傍から 出発してもわずかにずれた別の停留点に収束してしまう例が実測で確認された （交点近傍から開始しても距離が $0.015$ 程度で頭打ちになるケース）。垂直条件を介した間接的な反復ではなく、 実際のユークリッド距離 $\lVert C(t)-S(u,v)\rVert$ を目的関数として直接下げる座標降下パターンサーチ （localPolish()）を最終段に追加し、収束後の結果をさらに追い込む。

5.4　複数交点への対応

最近点は一般に一意だが、交点は曲線が曲面を何度も貫けば複数存在しうる。5.1のマルチスタートで得た 候補それぞれを独立に磨き上げたのち、距離が閾値（$0.01$）未満のものだけを抽出し、パラメータ空間で 近接する解同士を重複除去して、別々の交点として列挙する（computeAllIntersections()）。

これらの追加処理により、厳密法・近似法のどちらを選んでも、既定の複数次数（2〜5次）・意図的に 複数回交差させた形状で、総当たり探索と一致する結果が得られることを確認している。

62方式の比較

厳密法（ベジエクリッピング）近似法（サンプリング）
根の取りこぼし 凸包が区間内の根の存在を数学的に保証／棄却する サンプル間隔より細かい符号変化を見落とす可能性がある
1反復あたりの計算 多項式の積・和（畳み込み）と凸包計算がやや重い 評価点での単純な点・接線評価の繰り返しで軽量
区間の縮小率 凸包の形状に応じて変動（きれいに1点へ収束しやすい） サンプル数で決まる固定比（曲線約1/4、曲面約1/2）
原論文との対応 式(1)〜(5)のベジエクリッピング法にほぼ忠実 同じ枠組みを単純化した簡易近似
実装の複雑さ ベジエ積・凸包・分割など多くの補助関数が必要 既存の評価関数だけで実装できる

7参考文献

  • 西田友是,「有理ベジエ曲線・曲面間の最短距離検出法」, Visual Computing / グラフィクスとCAD合同シンポジウム 資料（本ツールの元になった解説論文）
  • T. Nishita, T. Sederberg, M. Kakimoto, "Ray Tracing Trimmed Rational Surface Patches," Computer Graphics, Vol.24, No.4, pp.337-345, 1990.
  • T. W. Sederberg, T. Nishita, "Curve intersection using Bézier clipping," Computer-Aided Design, Vol.22, No.9, 1990.（ベジエクリッピング法の基本文献）