NURBS曲面の包含箱・包含球と交差判定ベジエクリッピング法の考え方——距離関数のベジエ形式と、制御点の凸包による棄却——を、包含体の計算にも曲面の交差判定にも一貫して使う

NURBS曲面を有理ベジエパッチに分解すると、各パッチ上の点は制御点の凸結合になります。この凸包性と、距離関数がベジエ関数で書けることを使うのがベジエクリッピング法です。この文書では、同じ道具で「包含箱・包含球の計算」と「曲面どうしの交差判定」の両方が解けることを示します。さらに、包含体を前処理に使って、交差判定する曲面(パッチ)を絞り込む方法までを説明します。ページ内のデモで、クリッピングが両方の用途で同じように働くようすを確かめられます。

1全体の流れ — ひとつの道具で二つの仕事

曲面を直接扱うと、最大値の計算も交差の計算も難しくなります。そこで、まず曲面を有理ベジエパッチの集まりに分解します。次に、求めたい量を距離関数として書きます。距離関数は、ベジエ形式の制御係数で表せます。最後に、制御点の凸包を使って、答えが入り得ない範囲を捨てていきます。この「距離関数をベジエ形式で書き、凸包で捨てる」手順がベジエクリッピング法の考え方です。

NURBS曲面有理ベジエパッチ ベジエクリッピング 包含箱・包含球曲面の交差判定 箱・球交線 ノット+重みに分解(§2) 距離関数をベジエ形式で書き、制御点の凸包で残す範囲を絞る(§3・§4)同じ処理を両方で使う 最大値を追う(§5〜§7)帯に入る範囲を絞る(§8) 前処理(§9)交差判定するパッチを絞り込む
図1 処理の流れ。包含体の計算(上)と曲面の交差判定(下)は、どちらも中央の「ベジエクリッピング」の考え方で解きます。上で求めた包含体は、下の交差判定の前処理に使えます。

要点:ベジエクリッピング法は、包含体の計算にも交差判定にも使える

次の表のとおり、どの仕事でも「距離関数 fを作り、その制御係数の凸包からf が水準を超え得る範囲だけを残す」ことを繰り返します。違うのは、距離関数の作り方と「水準の置き方」だけです。

仕事距離関数 f残す範囲(凸包で判定)得られるもの
包含箱(③)軸方向の座標(x, y, z)=座標平面までの距離f が、見つけた最大値 L を超え得る範囲箱の6面と接点
包含球(④)球の中心からの距離(の2乗)f が、現在の半径(の2乗)を超え得る範囲球の中心・半径と接点
曲面の交差判定(§8)相手の曲面を包む平行2平面(FAT平面)までの距離f が、帯 [lo, hi] に入り得る範囲交線(の線分列)
交差の前処理(§9)包含体(箱・球)と、相手の帯・包含体との位置関係重なり得るパッチ交差判定する曲面(候補)

どの仕事でも、有理曲面では距離関数が分数になります。ただし、その符号や水準との大小は分子だけで判定できます(§4.2)。分子の制御係数は、重みを掛けた wij(fij − 水準) です。したがって、非有理のときと同じ次数のベジエ関数のまま扱えます。パッチの分割も、どちらの仕事でも斉次座標のde Casteljauです。

4つの包含体は、次のとおりです。番号は開発の順序でもあり、①②で骨組みを作り、③で厳密な最大値の求め方を加え、④で球の中心の最適化を加えます。

方式最大・最小の求め方性質
① 包含箱(制御点)制御点の x,y,z 座標の最大・最小一瞬で求まる。必ず曲面を含むが、有理円弧などでは緩い
② 包含球(制御点)①の箱の中心から制御点までの最大距離同上。回転体の有理円弧では特に大きくなる
③ 包含箱(厳密)曲面上の真の最大・最小を、パッチの分割で追い込む6面が曲面に接する(許容誤差以内)
④ 包含球(厳密)四隅の最小包含球から成長させる(または箱の中心から縮める)球面が曲面に接し、最適性の証明つき

2NURBS曲面を有理ベジエパッチに分解する

2.1 NURBS曲面と斉次座標

NURBS曲面は、制御点 Pij、重み wij、ノットベクトルから決まるB-スプライン基底 Ni,pNj,q を用いて次のように書けます。

S(u,v) = Σi,j Ni,p(u) Nj,q(v) wij PijΣi,j Ni,p(u) Nj,q(v) wij

制御点に重みを掛けた Pwij = (wx, wy, wz, w) を4次元の点と見なすと、分子と分母をまとめた4次元の曲面は、ふつうの(非有理な)B-スプライン曲面になります。3次元の曲面は、これを最後の成分で割った(射影した)ものです。この見方が以降の計算の土台です。

2.2 ノット挿入によるベジエ分解

内部ノットの多重度が次数と同じになるまでノットを挿入すると、各ノット区間が独立した1枚のベジエパッチとして取り出せるようになります。u 方向が p 次、v 方向が q 次なら、パッチは (p+1)×(q+1) 個の制御点を持ちます。ノット ū を1つ挿入するとき、斉次座標の制御点は次のように更新されます(Boehmのアルゴリズム)。

Qi = αi Pi + (1 − αi) Pi−1,  αi = ūuiui+pui  (0 ≤ αi ≤ 1)

係数 α は 0 以上 1 以下なので、新しい制御点の重みは元の重みの凸結合です。元の重みがすべて正なら、分解後の重みもすべて正のままです。次の節の凸包性が使えるのは、この正値性のおかげです。

2.3 円弧・回転体は有理2次ベジエで厳密に表せる

中心 C、軸ベクトル UV の円(楕円)弧 P(θ) = C + cosθ U + sinθ V(角度範囲 θ0θ1、Δθ = θ1θ0)は、3つの制御点と重みで厳密に表現できます。

P0 = P(θ0),  P2 = P(θ1),  P1 = C + cosθm U + sinθm Vcos(Δθ/2),  w = (1, cos(Δθ/2), 1)

ここで θm は角度範囲の中央です。母線の弧と回転の弧のテンソル積(重みは積)を取ると、ドーナツ・球・円柱などの回転体がNURBSとして正確に得られます。このとき中間の制御点 P1 は円弧の外側へ 1/cos(Δθ/2) 倍だけ張り出します。90°の弧なら √2 倍です。§5で見る簡易版の「緩さ」の正体はこれです。

形状(本アプリ)パッチ数パッチの種類
ドーナツ163×3の有理ベジエ(重み cos45°を含む)
楕円体8単位球のNURBSにアフィン変換(重みは不変)
円柱4次数混在(円周方向2次・有理/高さ方向1次)
12母線は非有理2次、円周方向は有理2次
ユタティーポット32双3次の非有理ベジエ(重みはすべて1)

補足 本アプリの曲面は、すでに有理ベジエパッチに分解された形で持っています。以降の計算は「パッチの配列」だけを入力にするので、ノット挿入で分解した一般のNURBS曲面もそのまま渡せます。

3有理ベジエ曲面と凸包性

m×n 次の有理ベジエパッチは、ベルンシュタイン基底 BiBj を使って、次のように書けます。

S(u,v) = Σi,j Rij(u,v) Pij,  Rij = wij Bi(u) Bj(v)W(u,v),  W = Σi,j wij Bi(u) Bj(v)

すべての wij が正なら Rij ≥ 0 かつ ΣRij = 1 です。つまり曲面上の点は、制御点の凸結合で、制御点の凸包に入ります。ここから、包含体の計算に必要な次の2つが得られます。

  • 上界 凸関数 φ(座標成分、中心からの距離など)に対して、φ(S) ≤ ΣRij φ(Pij) ≤ max φ(Pij)(イェンセンの不等式)。曲面上の φ の最大値は、制御点での最大値を超えません。
  • 下界 四隅の制御点 S(0,0)、S(1,0)、S(0,1)、S(1,1) は曲面上の点そのものです。そこでの φ の値は、曲面上の実在の値なので、最大値の下界になります。

分割すると、上界は曲面へ近づきます。パッチを t = 1/2 で分割するには、斉次座標 (wx, wy, wz, w) の上で通常のde Casteljauを行えば、得られる小パッチは元の曲面を厳密に再現する有理ベジエパッチになります。パッチが小さくなるほど、制御網は曲面に二次の速さで近づきます(パッチの大きさを h として誤差 O(h2))。

4距離関数とベジエクリッピング

4.1 軸方向の距離関数 → 箱の6面

座標軸 a(x, y, z)に垂直な平面までの符号付き距離は、単に座標成分です。

da(u,v) = ea · S(u,v) (= x, y, z

3つの軸それぞれについて、全パッチ上での最小値と最大値を求めれば、その6つの値が包含箱の6面です。一般の平面 n·Sd でも同じ扱いができます(交差判定の FAT 平面がこれです)。

4.2 有理式は「分子」で判定できる

有理曲面では S = X/W なので、平面までの距離は次のように分数になります。

f(u,v) = n·Sd = n·Xd WW = Σi,j wij (n·Pijd) Bi(u) Bj(v)W(u,v)

分母 W は正なので、f の符号や「f がしきい値 L を超えるか」は、分子だけで決まります。分子は係数 wij(n·Pijd) のベジエ関数で、非有理の場合と同じ次数のままです。重みが正なので、係数が0以下になる条件は n·Pijd と同じです。

4.3 ベジエクリッピングとは

ベジエクリッピング(Sederberg–Nishita)は、方程式 f(t) = 0 の根を求める方法です。まず、f をベジエ形式で書き、制御点 (i/n, fi) を作ります。次に、その凸包が横軸(f = 0)と交わる t の範囲を求め、根はその範囲にしかない、として区間を縮めます。縮めた区間で同じことを繰り返します。単純な根では、収束は二次であることが知られています。

t0, t1] ← 凸包 conv{ (i/n, fi) } と 直線 y = 0 が交わる t の範囲

凸包が横軸に触れない区間には根がないので、その区間は捨てられます。この「凸包で捨てる」判定が、あとの応用のすべての土台です。

4.4 同じ道具を、二つの目的に使う

「0」の代わりに別の水準を置くだけで、同じ処理が包含体の計算にも交差判定にも使えます。

  • 交差判定:水準は「帯 [lo, hi]」。相手の曲面をFAT平面で包み、こちらの曲面からその平面までの距離 f を考えます。f が [lo, hi] に入り得ない範囲は、相手に届かないので捨てます。残る範囲は、「f − hi ≤ 0 になり得る範囲」と「f − lo ≥ 0 になり得る範囲」の共通部分です。
  • 包含箱:水準は「これまでに見つけた最大値 L」。f を座標成分にして、fL > 0 になり得る範囲だけを残します。L は、曲面上の実在点(サブパッチの端点)での値で更新します。残る範囲が最大値の位置へ絞られ、L が最大値に収束します。
  • 包含球:f を「中心からの距離」に替えるだけ。距離の2乗は 2n 次のベジエ関数になり、水準は半径の2乗です(§7.4)。西田・中村(2024)は、この距離関数の導関数の制御点でクリップして、最遠点を求めています。

共通する部品

① 距離関数を、制御係数 wij(fij − 水準) のベジエ関数(有理の分子)として作る。② 制御点の凸包で、残す範囲を決める(残らなければ捨てる)。③ 残らない範囲を捨てるか、二分割(曲面は4分割)して、斉次座標のde Casteljauで厳密なサブパッチに切り出す。この3つは、包含体の計算と交差判定でまったく同じです。デモ4(§8)で、同じクリッピングが「最大値」と「交差」の両方に働くようすを確かめられます。

4.5 本実装での位置づけ

  • 曲面の交差判定(§8):区間を縮める本来のクリッピングです(clipToSlab)。ただし、凸包と軸の交点を直接求める代わりに、係数を3段階まで二分割して、係数が同符号でない区間の外包を取る近似実装です。
  • 包含体の計算(§6・§7):「凸包による棄却+二分割」(分枝限定)で実装しています。区間を縮める版は、水準 L でクリップするデモ4や、西田・中村(2024)の導関数によるクリップで実現できます。

5簡易版:制御点から求める(① ②)

5.1 方法

§3 の上界をそのまま答えにします。全パッチの制御点について、各座標の最小・最大を取れば①の包含箱です。②の包含球は、①の箱の中心 c から全制御点までの距離の最大値を半径にします(球は凸なので、制御点をすべて含めば凸包、つまり曲面も含みます)。

mina = min Pij,a,  maxa = max Pij,a,  c = (min + max)/2,  r = max |Pijc|

計算量はパッチ数×制御点数だけで、必ず曲面を含みます。ただし、制御点は曲面の外にあるので、その分だけ緩くなります。

5.2 なぜ緩くなるか

§2.3 のとおり、有理円弧の中間の制御点は円弧の外側へ 1/cos(Δθ/2) 倍だけ張り出します。ドーナツ(大半径 2.2、小半径 0.9)では、外周の角にあたるパッチに制御点 (3.1, 0.9, 3.1) があり、原点からの距離は √(2×3.1² + 0.9²) ≈ 4.475 です。曲面の実際の最遠点までの距離 3.1 に比べて約44%大きくなります。次のデモで、重みを動かして緩さを確かめてください。

デモ1 制御点の見積もりと真の値2次元の有理ベジエ曲線で、簡易版と厳密版を比べます。制御点はドラッグで動かせます。
曲線 簡易版(制御点) 厳密版(曲線に接する) 接点

150°の円弧は重みが cos75° ≈ 0.26。中間の制御点が大きく外へ飛び出すので、簡易版はとても緩くなります。重みを1にすると放物線になり、制御点は曲線に近づきます。

6厳密な包含箱(③):分割で最大値を追い込む

6.1 アルゴリズム

求めるのは、凸関数 φ(= ±x, ±y, ±z)の曲面上の最大値です。各サブパッチには2つの値を持たせます。

  • 上界 U サブパッチの制御点での φ の最大(§3の凸包性による)
  • 下界 L これまでに評価した曲面上の点(サブパッチの四隅)での φ の最大。全体で共有する
// 入力: パッチ集合, 凸関数 φ, 許容誤差 tol
全パッチをヒープ(キー = 上界 U)に入れ、四隅で L を更新
while ヒープが空でない:
    P ← U が最大のサブパッチ
    if U(P) ≤ L + tol: break           // これ以上 tol を超えて改善しない
    P を 4 分割(u,v を 1/2 で、斉次座標の de Casteljau)
    for 子 C in 4分割の結果:
        C の四隅で L を更新、C の制御点から U(C) を計算
        if U(C) > L + tol: ヒープに入れる    // そうでなければ捨てる(凸包による棄却)
返り値: max φ ∈ [L, L + tol]、L を与えた曲面上の点 = 接点

上界が最大のものから分割する(最良優先探索)ので、最大値に無関係な場所を細かく分割する無駄が少なくなります。終了時には、最大値が L(曲面上の実在点の値)と L + tol の間に閉じ込められます。箱の面には保証付きの上界を使うので、曲面は必ず箱に含まれ、接点との隙間は tol 以下です。

6.2 1次元で見る

次のデモは、同じ考え方を曲線で1ステップずつ見せます。右のグラフは、選んだ方向の距離関数 f(t) です。橙の横線が各区間の上界 U、赤の横線が下界 L です。L + tol 以下の区間は捨てられ、残った橙の区間が最大値の位置へ絞られていきます。

デモ2 分枝限定で最大値を追い込む「1ステップ」を押して、上界が最大の区間が二分割されるようすを見てください。
残っている区間(上界 U) 捨てた区間 下界 L(曲線上の実在点) 現在の最大の上界

左:曲線と、方向に垂直な2本の直線(赤=L、緑破線=現在の最大のU)。この2本の間に箱の面があります。

6.3 収束と計算量

制御点と曲面の差は、サブパッチの大きさ h の2乗に比例して小さくなるので、UL ≤ tol にするには、サブパッチをおよそ h ≈ √tol の大きさまで分割すれば足ります。最大点が1点なら、各段階で残る候補は数個で済みます。

一方、最大値が円周のように連続して並ぶ場合(ドーナツの外周の赤道、円柱の上下の縁など)は、円周に沿ったサブパッチがすべて同じ値の近くに残るので、分割数が増えます。従来の上界(制御点での距離の最大)では、ドーナツの最遠点探索(§7)で許容誤差(箱の対角線比)を 1e-4、1e-5、1e-6 と厳しくすると、分割数は 672、1712、6096 回と、おおよそ √10 倍ずつ増えました。この増加は、§7.4 の「2n次の分子による判定」で解消できます(ドーナツ・円柱・球では分割が0回になります)。

7厳密な包含球(④):中心を動かして半径を縮める

7.1 箱の中心は最適とは限らない

箱の中心を球の中心にして、距離関数 D(u,v) = |Sc| の最大値を求めれば、③と同じ分枝限定で半径 R0 が得られます(φ を中心からの距離にするだけです。距離は凸関数なので、§3の上界がそのまま使えます)。しかし箱の中心が最適な中心とは限りません。たとえば、150°の円弧は、弦を直径とする円が最小の円ですが、箱の中心はそこからずれています。

7.2 中心を動かして小さくする

中心 c に対する半径 R(c) = max|Sc| は、凸関数の最大なので c の凸関数です。局所的に小さくなる中心は、そのまま大域的な最適解になります。ただし折れ目のある関数なので、方向を変えて試すだけでは効率が悪く、次のように進めます。

// 支持点 S: 曲面上の実在の点の集合(最小包含球の下界を与える)
c ← ③の箱の中心;  R ← R(c);  S ← ③の接点 + c での最遠点
loop:
    ball ← MEB(S)                       // S の最小包含球(Welzl法)。半径 r(S)
    if R − r(S) ≤ tol: 終了             // 証明つきの終了条件
    t ← ball の中心                     // 中心の移動先の候補
    (R(t), 最遠点 p) ← 分枝限定で厳密に評価
    if R(t) < R: c ← t;  R ← R(t)       // 半径が小さくなる中心だけ採用
    p が ball の内側: 評価が粗いだけ → 精度を上げて再評価
    else: S に p を追加

なぜ終了できるか。S は曲面上の実在の点の集合なので、曲面を含むどの球も S を含みます。したがって S の最小包含球の半径 r(S) は、曲面の最小包含球の半径 R* を超えません。一方 R(c) は曲面を含む球の半径なので R* 以上です。つまり

r(S) ≤ R* ≤ R(c),  R(c) − R* ≤ R(c) − r(S) ≤ tol

となり、終了時に「これ以上、許容誤差以上は小さくならない」ことが保証されます。本アプリの各形状で実際に必要だった最遠点探索は、1〜8回程度でした。

7.3 最小包含球(Welzl法)

支持点の集合 S の最小包含球は、境界上の点(3次元では最大4点)で決まります。Welzl法は、点を1つずつ加えて「その点が現在の球の外なら、その点を境界に加えて再帰する」方法で、乱択順なら期待線形時間です。境界点から球を作る公式は次のとおりです(p0 を原点に取り、a, b, c を他の点との差とします)。

  • 2点:中点、半径は距離の半分。
  • 3点:3点を通る円の中心 o。平面内で 2a·o = |a|²、2b·o = |b|² を満たす点として求める。
  • 4点:中心の原点からの位置 o は 2a·o = |a|²、2b·o = |b|²、2c·o = |c|² の3×3連立方程式の解。
デモ3 中心を動かして半径を縮める「次の反復」で、支持点の追加・中心の移動・半径の縮小を1回ずつ見ます。
箱の中心の円(移動前) 現在の最良の円 支持点の最小包含円 r(S) 試した中心 支持点

150°の円弧は解析解が「弦を直径とする円」(半径 1.5 sin75° ≈ 1.449)で、1反復でそこへ到達します。90°の円弧は箱の中心がもう最適なので、まったく動きません。重みで偏らせた曲線は、数回の反復で少しずつ縮みます。

7.4 距離関数を 2n 次のベジエ関数として扱う(2n次の分子による判定)

§7.2 の最遠点探索では、サブパッチの上界に「幾何の制御点(m,n 次)での距離の最大」を使っていました。簡単ですが、距離が一定になる形状(ドーナツの赤道、円柱の縁、球)では分割が増えます。制御点が曲面の外へ張り出す分(有理円弧では 1/cos(Δθ/2) 倍)だけ上界が甘くなり、許容誤差を満たすまで 1/√tol に比例した分割が必要になるためです。

そこで、距離の2乗を、そのままベジエ関数として扱います。中心 c からの距離の2乗は、次のように分数になります。ここで ha,ij は斉次座標の成分(重み×座標)です。

D²(u,v) = |Sc|² = Σa Aa(u,vW(u,v,  Aa = Σi,j (ha,ijca wij) Bi(u) Bj(v)

ベルンシュタイン基底の積の公式 BimBi′m = [C(m,i)C(m,i′)/C(2m,i+i′)] Bi+i′2m を使うと、分子・分母はどちらも (2m, 2n) 次のベジエ関数になり、制御係数 NklDkl が得られます。分母の係数 Dkl は、重みが正なら正です。

D² = Σ Nkl Bk BlΣ Dkl Bk Bl  ⇒  max D² ≤ maxkl NklDkl

右辺の不等式は、D² が Nkl/Dkl の凸結合であることから従います。これが上界です。非有理(w = 1)なら Dkl = 1 で、Nkl は 2n 次の距離関数の制御係数そのものです(西田・中村(2024)の dk)。

  • 従来の上界を超えません。Nkl/Dklwijwi′j′(pij·pi′j′)(p = Pc)の加重平均で、p·p′ ≤ max|p|² だからです。
  • 距離が一定なら、上界が半径にちょうど一致します。円弧をその円の中心から測る場合、NR²D が恒等的に成り立ち、ベルンシュタイン係数は一意なので Nkl/Dkl = R² です。分割は要りません。
  • 条件は、重みが正であることだけです。

7.5 四隅からの成長

§7.2 は、箱の中心から出発して半径を縮める方式(上から近づく)でした。もう一つの方式は、四隅の最小包含球(下界)から成長させる方式です。

S ← 全パッチの四隅(曲面上の実在点);  ball ← MEB(S)     // r(S) は真の最小半径 R* の下界
loop:
    c ← ball の中心;  (R(c), 最遠点 p) ← 2n次の分子による判定で厳密に評価
    if R(c) − r(S) ≤ tol: 終了                 // 証明つき
    p が ball の内側: 評価が粗いだけ → 精度を上げる
    else: S に p を加える;  ball ← MEB(S)

西田・中村(2024)の「四端点を含む球から出発し、最遠点の方向へ中心を動かす」方式と、同じ枠組みです。厳密な包含箱を先に求める必要がなく、パッチ単位で見ると、多くのパッチで四隅の球がそのまま曲面を含みます(ティーポット32面のうち24面)。

8曲面の交差判定 — 同じクリッピングで交線を求める

ここが要点

§5〜§7 で最大値を追うために使ったのと同じ道具——距離関数のベジエ形式と、制御点の凸包による棄却——で、曲面どうしの交差も解けます。違うのは、水準が「これまでの最大値 L」ではなく、「相手の曲面を包む帯 [lo, hi]」になることだけです。

8.1 問題

2枚の曲面 S1(u,v) と S2(s,t) の交差は、S1(u,v) = S2(s,t)、つまり未知数4つ・条件3つの連立方程式です。解は1次元、つまり交線になります。ベジエパッチ1対ごとに、(u,v,s,t) の4次元の箱を、凸包で絞っていきます。

8.2 相手の曲面を包む「FAT平面」

曲面 Q について、代表となる平面 n·x = d を作ります。法線は四隅から、平面は制御点の重心を通します。次に、Q の制御点から、平面までの符号付き距離の範囲を求めます。

lo = minij (n·Qijd),  hi = maxij (n·Qijd)

重みが正なら凸包性(§3)より、Q は平行2平面の間(帯 [lo, hi])に必ず入ります。この帯をFAT平面と呼びます。曲線どうしの交点計算で使うFAT Lineの、曲面版です。

8.3 もう一方の曲面 P の距離関数と、帯へのクリップ

P から平面までの距離は、§4.2 と同じ形の分数です。

f(u,v) = n·P(u,v) − d = Σ wij (n·Pijd) Bi(u) Bj(v)W(u,v)

f が帯 [lo, hi] に入り得るのは、次の2つが同時に成り立つ範囲だけです。W > 0 なので、分子の係数だけで判定できます。

G1,ij = wij (n·Pijd − hi)  が ≤ 0 になり得る
G2,ij = wij (n·Pijd − lo)  が ≥ 0 になり得る

実装(clipToSlab)では、u 方向のクリップのために、v 方向の添字 j について G1 の最小と G2 の最大を取り、u の1変数ベジエ関数(包絡線)にします。その係数を最大3段階まで二分割し、係数が同符号でない区間(凸包が軸に触れる区間)だけを残して、その外包(に3%の余裕)を残す u の範囲とします。v 方向も同様です。残る範囲が元の面積の92%以下なら、その範囲を斉次座標のde Casteljauでサブパッチとして切り出します。これを、PQ の帯へ、QP の帯へ、と交互に行います。

8.4 全体のアルゴリズム

// パッチ対 (P, Q)。flatTol=平坦の許容, maxDepth=最大の再帰の深さ
ssiRec(P, Q, depth):
    if P,Q の制御点の境界箱が重ならない(余裕 flatTol): return      // 最も安い棄却
    Q の FAT平面 [lo,hi] で P をクリップ(残らなければ return)        // ベジエクリッピング
    P の FAT平面 [lo,hi] で Q をクリップ(残らなければ return)
    if depth ≥ maxDepth または P,Q がどちらも十分平坦(平面からのずれ < flatTol):
        P,Q を各2つの三角形で近似し、三角形どうしの交線を線分として出力    // 末端
        return
    大きい方(平面から遠い方)を 4分割して再帰                        // 有理de Casteljau

末端では、パッチは平面に近いので、三角形どうしの交線(線分)で置き換えます。得られる線分の列が交線です。有理曲面でも、分子の次数は非有理と同じなので、アルゴリズムは変わりません。重みをすべて1にしたコピーで同じ計算を走らせる「非有理とみなす判定」と比べると、重みが交線に与える影響を確かめられます。

8.5 同じクリッピングが「最大値」と「交差」の両方に働く(デモ4)

次のデモは、2次元の有理ベジエ曲線で、まったく同じクリッピング処理を2つの目的に使います。「最大値」は包含箱(③)の考え方で、直線 n·x の値の最大を追います。「交差」は、直線とその両側の帯 [lo, hi](相手の曲面のFAT平面に相当)に、曲線が入る範囲を絞ります。下段のグラフは、直前に処理した区間の「分子の制御点」と凸包です。凸包が水準の側に入る範囲だけが残されます。

デモ4 同じクリッピングを「最大値」と「交差」に区間を縮める本来のベジエクリッピングです。凸包が水準の側に入らない範囲を捨て、縮まらなければ二分割します。
目的
残っている区間 末端(帯に収まった区間) 捨てた部分 水準 L(最大値の下界)

「最大値」で 1ステップずつ進めると、凸包が0の軸の上に出る範囲(緑)へ区間が絞られ、水準 L が最大値へ近づきます。「帯に入る範囲」に切り替えると、同じ凸包の判定で、曲線が帯を横切る場所だけが残ります。区間が帯の中に完全に収まったら「十分平坦」として末端にします。曲面の交差判定(§8.4)の末端と同じ考え方です。

補足 デモでは、区間を縮める本来のクリッピングを、最大値の探索にも使っています。本アプリの包含体の計算(§6・§7)は、簡単のため「凸包による棄却+二分割」で実装しています。一方、曲面の交差判定(§8.3・§8.4)は、区間を縮めるクリッピングです。

8.6 包含体の計算と交差判定で共通の部品

部品包含体の計算(§5〜§7)曲面の交差判定(§8)
距離関数座標成分(軸方向)、中心からの距離(の2乗)相手のFAT平面までの符号付き距離
有理の扱い分子の係数 wij(fij − 水準)分子の係数 G1G2(同形)
水準見つけた最大値 L(球では半径)帯 [lo, hi]
捨てる条件制御点の凸包が水準の側に入らない制御点の凸包が帯に入らない
分割斉次座標のde Casteljauで4分割斉次座標のde Casteljauで4分割・切り出し
正の重みが必要必要(凸包性のため)必要(凸包性のため)

9包含体で前処理して、交差判定する曲面を抽出する

9.1 目的

曲面Aは、たくさんのベジエパッチから成ります。曲面Bとの交差判定(§8)を、Aの全パッチに対して行うのは無駄です。Bから遠いパッチは、そもそも交差しません。そこで、包含体で安く判定して、交差判定する曲面(候補)だけを §8 のクリッピングに渡します。包含体の計算もクリッピングの考え方で求めたものなので、包含体の計算と交差判定が、一つの流れとしてつながります。

9.2 階層包含体

曲面Aについて、包含体を2階層で持ちます。

  • 上位:NURBS全体の包含体(箱・球)。球は、全パッチの点集合の最小包含球です。パッチごとの球を包む外接球よりも、小さくなります。
  • 下位:ベジエパッチごとの包含体。パッチごとに、箱・球を、制御点版・厳密版のどちらでも求められます。

曲面Bは1枚のパッチなので、包含体は1つです。

9.3 候補パッチを抽出する手順

  1. Aの全体の包含体と、Bの包含体を比べます。箱なら各軸で区間が重なるか、球なら中心間の距離が半径の和以下か、で判定します。判定には ssiRec と同じ余裕 flatTol を見込みます。重ならなければ、Aの全パッチを一括で除外します。
  2. 重なるときは、Aの各パッチの包含体とBの包含体を比べ、重ならないパッチを除外します。
  3. さらに、BのFAT平面(§8.2)と、Aのパッチの包含体を比べます。包含体が帯 [lo, hi] の外にあれば、そのパッチはBに届かないので除外します。箱の帯への射影は、次のように求められます。
箱の射影 = m ± Σk |nk| hk,  m = n·(箱の中心) − d,  hk = 箱の半分の幅  (球は m ± 半径)
  1. 残った候補パッチだけを、§8.4 の ssiRec に渡します。
ティーポットと曲面Bの交差。候補パッチの境界を細いオレンジで表示
候補パッチの境界(細いオレンジ)。交線(赤)が通るパッチが、32枚のうち10枚に絞られています。淡い色のパッチは除外されました。
ティーポットの全体とパッチごとの包含箱(階層表示)
階層包含体(箱)。外側の大きな箱がNURBS全体、内側の箱がパッチごとです。除外されたパッチの箱は薄く表示されます。

9.4 安全性

包含体は曲面を必ず含み、Bの帯もBを必ず含むので、除外したパッチとBは交差しません。この判定で交線を見落とすことはなく、結果は前処理なしと同じになります。実際に、7形状×60構成×有理/非有理×8種=6720通りで、前処理ありの交線が前処理なしと完全に一致し、除外したパッチに交差が隠れていた例は0件でした。

Bの包含体は、Bを動かすたびに、現在の形で必ず計算し直します。古い包含体を使うと、交差を見落とすからです。Aの包含体は、形状を変えるまで再利用できます。

9.5 効果

Bを高さ −3〜+3 の範囲で動かし(61点、6.5角の傾いた平面)、Aのパッチのうち候補として残る割合を測りました(小さいほど絞れています)。

形状箱のみ箱+FAT平面球のみ球+FAT平面
ドーナツ69%34%100%66%
楕円体75%34%100%48%
円柱100%57%100%62%
69%22%100%27%
ユタティーポット61%18%100%31%
球(NURBS)77%43%100%52%
  • FAT平面を併用すると、絞り込みが大きく進みます。Bは大きな傾いた平面なので、その軸平行箱・球はとても大きく、多くのパッチと重なります。一方、Bの平行2平面は、Bにぴったり沿います。
  • 球は、大きな平面に対しては箱より弱いことが分かります。
  • 計算時間の短縮は小さいです。ssiRec の最初の判定(制御点の箱)が、すでに同じ役割を果たしているからです。ティーポットにBを交差させた例では、ssiRec の判定ノードが 176 から 154 に減りましたが、交線は前処理なしと同じ130線分でした。パッチ数が非常に多い場合や、Bが遠い場合(全体で一括除外)ほど、有利になります。
  • 前処理にかかる時間は小さく、ティーポット32パッチ分の箱(厳密)の計算が約3ms、Bの更新ごとの計算と判定が0.01ms程度です。

10実行例と数値

10.1 包含体:4方式の比較

本アプリ(nurbs_ssInt_bv.html)の組み込み形状で、4方式を比べた結果です。許容誤差は箱の対角線の 1e-5 倍です。

ティーポットの制御点による包含箱(ワイヤフレーム)
① 制御点による包含箱。制御点はティーポットの外にはみ出します。
ティーポットの厳密な包含箱(半透明)と接点
③ 厳密な包含箱。6つの面それぞれに曲面上の接点(●)があります。
ドーナツの制御点による包含球(半透明)
② ドーナツの制御点による包含球(R = 4.475)。有理円弧の制御点が外へ張り出すため、大きく余ります。
ドーナツの厳密な包含球(半透明)
④ 同じドーナツの厳密な包含球(R = 3.100)。外周の赤道に接しています。
ティーポットの厳密な包含球。点線が移動前の球、実線が移動後の球
④ ティーポットの厳密な包含球。○が箱の中心、●が移動後の中心で、線がその軌跡です。点線が移動前(箱の中心)の球で、実線が移動後の球です。半径は 3.006 から 2.914 へ、約3.1%小さくなりました。
形状③箱の体積
(厳密/制御点)
②球の半径
(制御点)
④箱の中心の
球の半径 R0
④移動後の
半径 R
縮小球の体積
(厳密/制御点)
最遠点
探索
ドーナツ100.0%4.4753.10013.10010.00%33.2%1
楕円体100.0%2.8722.00002.00000.00%33.8%1
円柱100.0%2.5422.05192.05190.00%52.6%1
69.8%1.9251.45781.41133.19%39.4%7
ユタティーポット98.6%3.0593.00632.91363.08%86.4%5
1枚の双3次ベジエ92.4%4.4344.41714.34771.57%94.3%3
傾き・重みを付けた1枚*98.7%4.6234.62124.60810.28%99.0%3

読み取れること ドーナツ・楕円体・円柱は軸に沿った形なので、制御点の箱がそのまま厳密な箱と一致します(軸方向の端の制御点が曲面上にある)。一方、球は制御点版が大きく、ドーナツと楕円体では厳密版の3倍ほどの体積になります。対称な形状(ドーナツ・楕円体・円柱)では箱の中心が最初から最適で、④は動きません。非対称な壺・ティーポットでは、中心を動かして3%ほど半径が縮みます。
* 傾き・うねり・重み(最大3、最小0.3など)を与えた双3次有理ベジエ1枚(検証用)。

10.2 ④の高速化:2n次の分子による判定と、四隅からの成長

最遠点探索のコストを、従来の上界(制御点での距離の最大)と、§7.4 の 2n次の分子による判定で比べました。数値は分割回数(箱と球の合計)です。

形状従来の上界
+箱の中心から
2n次の分子
+四隅から成長
ドーナツ17120
円柱10120
楕円体560
球(距離が一定)12万超(上限で打切り)0
1625382
ユタティーポット9578

距離が円周状(ドーナツの赤道、円柱の縁)や曲面全体(球)で一定になる形状は、分割が0回になりました。壺のように最大点が円周の内側にある場合は、減るものの0にはなりません。効いているのは主に 2n次の分子で、四隅からの成長は、厳密な包含箱の計算が要らなくなる利点です。4通りの組み合わせの半径は、許容誤差の範囲内で一致しました(乱数200曲面でも一致)。

10.3 交差の前処理

前処理の効果は §9.5 の表のとおりです。ティーポットにBを交差させた例では、32パッチのうち10パッチが候補になり、交線は前処理なしと同じ130線分でした。

10.4 検証

検証は、曲面上の密なサンプル点との比較(すべて箱・球に含まれ、箱の面は最大値との差が許容誤差の範囲内)、標本点に対する別方式の最小化との比較、乱数で作った有理ベジエ曲面200個(重みは約0.07〜14、次数混在)での包含・証明ギャップ・縮小の確認で行いました。計算時間は、許容誤差 1e-5 でどの形状も約35ms以下でした。

11実装の対応表と使い方

関数段階役割
bvBoxControl制御点の座標の最小・最大から包含箱を作る
bvSphereControl①の箱の中心から制御点までの最大距離で包含球を作る
bvFlatパッチを斉次座標の平坦な配列にする(高速な分割のため)
bvHalve / bvSplit4斉次座標のde Casteljauで、1方向に二分割/4分割する。任意次数に対応
BVHeap上界が最大のサブパッチを先に取り出す最大ヒープ
bvUpper / bvMaximizeWith③④上界の評価(差し替え可能)と、凸関数の曲面上の最大値を求める分枝限定
bvBoxExact±x, ±y, ±z の6回の最大化で、厳密な包含箱と接点を得る
bvUpperNum距離²を (2m,2n) 次の有理ベジエ関数と見て、制御係数の比で上から押さえる(§7.4)
bvBallFrom / bvMEB2〜4点から球を作る公式と、Welzl法による最小包含球
bvFarthest中心から曲面上の最遠点を求める(判定は 2n次の分子/従来を選択)
bvSphereGrow / bvSphereExact四隅から成長(§7.5)/箱の中心から縮める(§7.2)
bvVolume / bvHier階層包含体を1つ求める/全体とパッチごとの階層を求める(§9.2)
bvSlabOf / bvVolInSlab / bvCullCandidates前処理Bの平行2平面/包含体と帯の判定/交差判定する候補パッチの抽出(§9.3)
meanPlane / slabOf交差FAT平面(平均平面と帯 [lo, hi])を作る(§8.2)
clipToSlab交差相手の帯へのクリップ。分子の係数 G1G2 から残す範囲を求め、サブパッチを切り出す(§8.3)
split4RGrid / extractPatchUV交差有理de Casteljauによる4分割・サブパッチの切り出し
ssiRec / ssiCompute交差再帰による交線の計算(§8.4)/前処理つきで候補パッチだけを渡す
triTriSegment交差末端で、三角形どうしの交線を線分として取り出す
bvCompute / bvGetResultUI曲面ごとの結果キャッシュ。曲面Bの操作直後は約0.1秒待ってから再計算(表示用)

アプリでの使い方

  1. サイドバーの「★ 包含箱・包含球」で方式(①〜④)を1つ選びます。表示形式(ワイヤフレーム/半透明)と、対象(曲面A/曲面B)も選びます。
  2. 「表示する階層」で、NURBS全体/パッチごと/両方を切り替えられます。パッチごとの包含体は、そのパッチの箱・球です。
  3. 厳密版(③④)は「厳密版の許容誤差」スライダーで精度を変えられます。④では、距離判定(従来/2n次の分子)と出発点(箱の中心/四隅から成長)を選べ、「4通りを全形状で比較」で計算コストを比べられます。
  4. 「★ 交差の前処理」で、包含体による候補の絞り込みのON/OFF、使う包含体(箱・球、制御点・厳密)、Bの平行2平面(FAT平面)の併用を選べます。候補パッチは細いオレンジの輪郭で強調され、除外されたパッチは淡く表示されます。統計欄に、候補パッチの数と前処理の時間が出ます。
  5. 「前処理あり/なしを比較」で、9通りの前処理について、候補パッチの数・ssiRec の判定ノード数・交線の数・時間を比べられます。曲面Bを動かして交差する位置にしてから比べてください。
  6. 曲面Bの制御点や重みを動かすと、交差は毎回、Bの現在の形で計算し直されます。厳密版の包含体の表示は、操作が止まったあとに再計算されます。

12注意点

  • 重みは正であること。凸包性は、すべての重みが正(または同符号)のときに成り立ちます。負や0の重みを含む曲面には、この方法は使えません。
  • 許容誤差の意味。箱の面・球面には保証付きの上界を使うので、曲面は必ず含まれます。接点との隙間は、指定した許容誤差(箱の対角線に対する比)以下です。
  • 最大値が連続して並ぶ形状。円周状に同じ値が並ぶ場合、従来の上界では分割数が増えますが、§7.4 の 2n次の分子による判定で大きく減らせます(最大点が円周の内側にある場合は残ります)。暴走を防ぐため、1回の最大化には分割回数の上限(30万回)があり、上限に達した場合も「保証付きの上界」は有効です。
  • 退化したパッチ。球の極のように辺が1点に潰れたパッチも、重みが正なら問題なく扱えます。
  • 交差判定の精度。ssiRec の末端は、平面(三角形)による近似です。交線は線分の列で、精度は平坦の許容(flatTol)と最大の深さで決まります。前処理は、この結果を変えずに候補パッチを減らすだけです。
  • 前処理の効果。制御点の箱による判定は ssiRec の最初の判定と重なるので、時間短縮は小さくなります。Bの平行2平面(FAT平面)の併用が、絞り込みに効きます。
  • 軸に平行な箱のみ。向き付き箱(OBB)が必要な場合は、座標系を回転してから同じ計算をします。
  • トリム曲面。トリム(切り取り)された曲面には対応していません。切り取り前の曲面全体の包含体は、切り取り後の曲面も含みますが、その分だけ緩くなります。