境界面を守って死ぬメッシュ生成器 — 制約付きドロネーと Steiner 点

中点分割では侵入点が 4 から 5 へ戻り得ます。止めるのは基準点と t の規則です。

STL を一つ tetgen に食わせたら、プログラムがそのまま落ちました。ログの最終行は A segment and a facet intersect at point。形状を調べても壊れた箇所は見つかりません。面は閉じていて、自己交差もなく、別の商用メッシャは問題なく通します。ファイルの問題ではありませんでした。

Thingi10k データセットの有効なモデル 4408 個のうち約 8.5% で、tetgen は同じ壊れ方をします。この記事では、その 8.5% がどこから来るのかを追いかけます。原因は二系統です。一つは浮動小数点、もう一つは浮動小数点とまったく関係ありません。後者のほうが面白いです。

三角形分割はあなたの境界を覚えていない

ドロネー三角形分割(空外接円条件を満たす分割)は点集合の性質です。入力が点群なら完璧に動きます。問題は、CFD が点群を渡してくれないことにあります。

私たちが渡すのは PLC(piecewise-linear complex、区分線形複体 — 頂点・線分・多角形面が互いに整合して噛み合った複体)です。翼の表面、シリンダ壁、流入パッチ。これらはメッシュの中に辺と面そのままの姿で生き残らなければなりません。壁境界条件を近似された面に課すわけにはいきません。

ところが頂点だけ入れてドロネーを回すと、その辺が生き残る保証はまったくありません。別の辺がその上を横切ってしまいます。こうして消えたものを missing segment と呼びます。

TetWild や Quartet のように境界を近似する方式は、この問題を回避します。代わりに境界条件を投影して移す必要があり、その投影は損失があり、全単射(bijective)である保証もありません。壁面摩擦や熱流束をその上に課せば、精度はその分だけ飛びます。

侵入点のスライダを線分へ寄せると、直径円が先に赤くなり、さらに寄せて緑の線分が切れます。

Steiner points inserted: 0

intruder clearance = 58

スライダが変えるのは距離だけです。円が赤いあいだ線分はまだ生きていて、黄色い辺が横切ったときに missing になります。split once は基準点 $r$ を選んで、その位置で切ります。

消えた線分には必ず侵入点がある

線分を復元する方法は割ることだけです。ところが、どこで割るかがアルゴリズムの収束を決めます。やみくもに半分で割ると終わらない場合が出てきます。

判定基準は一つです。線分 $e = \langle v_1, v_2 \rangle$ の直径円(両端点を直径とする最小外接円) $D$ の中に他の頂点が一つもなければ、 $e$ は strongly Delaunay です。

\|x - \tfrac{v_1+v_2}{2}\| \;\ge\; \tfrac{1}{2}\|v_2 - v_1\| \quad \forall x \in V \setminus \{v_1, v_2\}

$V$ は PLC の頂点集合、左辺は直径円の中心までの距離です。この条件が成り立てば、 $e$ は必ずドロネー三角形分割の中にあります。

向きを反転して読むのが要点です。線分が消えたなら、直径円の中に必ず頂点がある。 その頂点たちを侵入点(encroaching point)と呼びます。逆は成り立ちません。侵入点があっても線分が無事に生き残ることはあります。もっと大きく膨らんだ外接円が空であれば、それで済むからです。

どこで割るかが収束を決める

基準点 $r$ は、侵入点のうち $(v_1, v_2, r)$ を通る外接円の半径が最大の点です。最も深く侵入した点ではなく、最も広く邪魔をする点を選びます。

分割位置は $r$ までの距離で決めます。 $R_1 = \|r - v_1\|$ 、 $R_2 = \|r - v_2\|$ 、 $L = \|v_2 - v_1\|$ とします。

t = \begin{cases} L/2, & R_1 > L/2 \;\text{かつ}\; R_2 > L/2 \\ \min(R_1,\, L - R_2), & \text{その他} \end{cases}

$t$ は $v_1$ から測った分割点までの距離です。二番目の場合は、 $v_1$ または $v_2$ を中心に $r$ へ届く円と線分が交わる位置です。この規則があってはじめて、反復が有限回で終わることが証明されます。

ここで一つ目の落とし穴が生じます。円と線分の交点なので、 $t$ は無理数になりえます。入力座標がすべて double でも、分割点の座標は double では表現されません。丸めた瞬間にその点は元の線分からわずかに外れ、アルゴリズムのすべての証明が無効になります。

迂回路は、座標を保存しないことです。分割点を $t v_1 + (1-t) v_2$ という式のまま持ち回れば、その点は定義上つねに線分の上にあります。こうした点を implicit point、そのうち二点の線形結合であるものを LNC(linear combination)と呼びます。orient3d と inSphere 述語を LNC を受け取れるよう拡張すれば、浮動小数点ハードウェアを使いながら符号を間違えずに済みます。

出力の $t/L$ を規則に当てると中点ではありません。0.453、0.516、0.698、0.315 です。四回とも、コードが切る範囲 0.12 $L$ から 0.88 $L$ の中にあります。step 2 の 0.698 は半分からかなり離れています。基準点 $r$ は最も深く入った点ではなく、 $(v_1, v_2, r)$ の外接半径がいちばん大きい点 #8、#15、#11、#4 です。侵入点数は 14、4、5、3 と行き来します。4 から 5 へ戻ったのは、新しい部分線分がそれまでなかった直径円の侵入を作ったからです。単調減少がなくても、分割 4 回、Steiner 点 4 個、部分線分 5 個で、すべての部分線分がドロネー辺になります。

空外接円で点を入れる Bowyer–Watsonは、点だけ見ればこの辺集合を作ります。境界の線分はその集合に自動では入りません。だからループが消えた線分を探して、もう一度切ります。

3D には答えが存在しない多面体がある#

2D では、交差しない線分の集合でありさえすれば、制約付きドロネー三角形分割は必ず存在します。3D では違います。

三角柱の上面を $\theta$ だけねじってみましょう。側面の四角形三つが平面性を失います。各四角形は対角線一本で割らなければならず、割り方は二通りあります。どちらを選ぶかは符号一つが決めます。

\mathrm{orient3d}(a,b,c,d) = \big((b-a) \times (c-a)\big) \cdot (d-a)

$(a,b,c)$ は外を向く三角形、 $d$ は残りの頂点です。値が負なら $d$ が内側にあり、その対角線は凸(convex)な辺、正なら凹(reflex)な辺です。

三面すべてを凹側に選ぶと Schönhardt の多面体(1928)になります。閉じていて、単純で、自己交差もありません。そして自分の六つの頂点だけでは、四面体に絶対に分割できません。凹の対角線を選ぶと候補 15 個が一度に失格になり、内側の点一つが 8 個の四面体でその穴を埋めます。

ここで、CDT 定理がなぜ線分の話しかしないのかが見えてきます。すべての線分が strongly Delaunay であれば、CDT は存在します。多面体をどう作り直すかではなく、線分をどこで割るかだけを決めればよいという意味です。分割は純粋に位相操作です。入力形状は 1 mm も動きません。

面の復元 — 掘り抜いて詰め直す

線分がすべて復元できたら、次は面です。PLC の面 $f$ をメッシュの辺が一本でも貫通すれば、その面は missing です。

手順はこうです。 $f$ を貫く辺に接する四面体を全部集めてキャビティを作ります。 $f$ の平面でキャビティを上下の半分 $C_1$ , $C_2$ に割ります。平面の上に載る頂点は両方に入れます。各半分の頂点で局所的なドロネー分割 $D_i$ を計算し、キャビティの内側に入る四面体だけを選んで詰めます。 $D_i$ は凸ですが $C_i$ は凹でありうるので、全部が使われるわけではありません。

問題は、 $C_i$ の境界三角形が $D_i$ に現れないときです。このときは反対側の四面体を足してキャビティを広げ、再計算します。この拡張が二つ目の失敗地点です。

Python で数える侵入点と分割回数#

2D で線分一本を復元する全体ループです。外接円判定 → 消えた線分の検出 → 侵入点の収集 → 基準点の選択 → 分割、そして反復。

import numpy as np
from itertools import combinations
 
def circumcircle(a, b, c):
    """三点の外接円(中心、半径)。共線なら (None, None)。"""
    (ax, ay), (bx, by), (cx, cy) = a, b, c
    d = 2.0 * (ax*(by-cy) + bx*(cy-ay) + cx*(ay-by))
    if abs(d) < 1e-12:
        return None, None
    ux = ((ax*ax+ay*ay)*(by-cy) + (bx*bx+by*by)*(cy-ay) + (cx*cx+cy*cy)*(ay-by)) / d
    uy = ((ax*ax+ay*ay)*(cx-bx) + (bx*bx+by*by)*(ax-cx) + (cx*cx+cy*cy)*(bx-ax)) / d
    ctr = np.array([ux, uy])
    return ctr, float(np.linalg.norm(ctr - np.asarray(a)))
 
def delaunay_edges(pts):
    """空外接円条件を通過した三角形の辺の集合。"""
    n = len(pts)
    edges = set()
    for i, j, k in combinations(range(n), 3):
        ctr, rad = circumcircle(pts[i], pts[j], pts[k])
        if ctr is None:
            continue
        rest = [m for m in range(n) if m not in (i, j, k)]
        if rest and np.linalg.norm(pts[rest] - ctr, axis=1).min() < rad - 1e-9:
            continue                      # 外接円の中に点があればドロネーではない
        edges |= {(i, j), (j, k), (i, k)}
    return {(min(a, b), max(a, b)) for a, b in edges}
 
def encroaching(pts, i1, i2):
    """線分の直径円の中に入った頂点 = 侵入点。"""
    mid = 0.5 * (pts[i1] + pts[i2])
    rad = 0.5 * float(np.linalg.norm(pts[i2] - pts[i1]))
    return [k for k in range(len(pts))
            if k not in (i1, i2) and np.linalg.norm(pts[k] - mid) < rad - 1e-9]
 
def recover_segment(points, chain, max_split=16):
    pts = [np.asarray(p, dtype=float) for p in points]
    for step in range(max_split):
        E = delaunay_edges(np.array(pts))
        gone = [s for s in range(len(chain)-1)
                if (min(chain[s], chain[s+1]), max(chain[s], chain[s+1])) not in E]
        if not gone:
            return np.array(pts), chain, step
        s = gone[0]
        i1, i2 = chain[s], chain[s+1]
        vd = encroaching(np.array(pts), i1, i2)
        v1, v2 = pts[i1], pts[i2]
        L = float(np.linalg.norm(v2 - v1)); u = (v2 - v1) / L
        r = max(vd, key=lambda k: circumcircle(v1, v2, pts[k])[1] or 0.0)
        R1 = float(np.linalg.norm(pts[r] - v1))
        R2 = float(np.linalg.norm(pts[r] - v2))
        t = L/2 if (R1 > L/2 and R2 > L/2) else (R1 if R1 <= R2 else L - R2)
        pts.append(v1 + u * float(np.clip(t, 0.12*L, 0.88*L)))
        chain = chain[:s+1] + [len(pts)-1] + chain[s+1:]
        print(f"  step {step}: 侵入点 {len(vd)}個, 基準点 #{r}, t/L = {t/L:.3f}")
    raise RuntimeError("分割回数の上限を超過")
 
rng = np.random.default_rng(20260729)
P = [np.array([0.0, 0.0]), np.array([10.0, 0.0])]
P += [rng.uniform([1.0, -3.0], [9.0, 3.0]) for _ in range(14)]
 
pts, chain, nsplit = recover_segment(P, [0, 1])
E = delaunay_edges(pts)
ok = all((min(chain[s], chain[s+1]), max(chain[s], chain[s+1])) in E
         for s in range(len(chain)-1))
print(f"Steiner 点 {len(pts)-len(P)}個, 分割 {nsplit}回, 部分線分 {len(chain)-1}個")
print(f"すべての部分線分がドロネー辺か: {ok}")

実行結果です。

  step 0: 侵入点 14個, 基準点 #8, t/L = 0.453
  step 1: 侵入点 4個, 基準点 #15, t/L = 0.516
  step 2: 侵入点 5個, 基準点 #11, t/L = 0.698
  step 3: 侵入点 3個, 基準点 #4, t/L = 0.315
Steiner 点 4個, 分割 4回, 部分線分 5個
すべての部分線分がドロネー辺か: True

twist θ = 18°

reflex split で $\theta$ を上げると、候補 15 個のうち有効な四面体が 0 になります。convex split は同じ形状を分け、Steiner point は 8 個で埋めます。黄色い対角線は固体の外へ出た辺です。

注目したいのは step 1 です。侵入点が 14 個から 4 個へ落ちたあと、step 2 で再び 5 個に増えます。新しく作った部分線分が、それまでなかった侵入関係を生むからです。それでも $t$ 規則のおかげで、単調に減らなくても有限回で終わります。半分で割るだけの実装が無限ループになる理由はここにあります。

残る二つの失敗 — 丸めと理論

tetgen が落ちる 8.5% には、二つの原因が混ざっています。8.5% なら $4408 \times 0.085 \approx 375$ 件です。シングルコア 5 時間を 4408 で割ると、モデルあたり約 4.1 秒です。

片方は丸めです。Steiner 点の座標を double にスナップした瞬間、入力 PLC がわずかに変形します。Shewchuk のフィルタ付き述語を使っても無駄です。述語が正確でも、入力がすでに間違っているからです。LNC の式として持ち回れば、この系統は消えます。

もう片方は数値と無関係です。キャビティの拡張は、二つの半分の分割の内部が互いに重ならないことを暗黙に仮定しています。ところが拡張の過程で、復元中の面の平面を越えて反対側の四面体を引き込む場合が存在します。二つの分割が交差してしまいます。著者らがテストした 4408 個のモデルのうち、ちょうど 2 個で発生しました。無限精度で計算しても失敗します。アルゴリズムそのものの穴です。

正確な数値型(CORE ライブラリ)で実装し直せば一つ目は解決しますが、二つ目は残ります。そして速度が実用範囲を外れます。中規模のファイル一つに数時間かかります。座標を $t \in (0,1)$ 一つで有理数パラメータ化し、間接述語を使うほうが、実際に使える折衷案です。4408 個のモデル全部がシングルコアで約 5 時間で処理されます。

中点で割ると復元ループが終わらない

$t/L$ が 0.5 でないのは、そこが終了条件だからです。出力は 0.453、0.516、0.698、0.315 で、侵入点は 14 から 4 へ落ちたあと 5 に戻りました。半分で切るだけでは、この増加を切れません。基準点 $r$ と $t$ の規則は、増えても有限回で止まるためにあります。

座標を double にスナップすると、線分上にあった点が外へずれます。述語が正確でも、入力がすでに違っていれば役に立ちません。丸めを直しても 2 個は空洞が平面を越えて重なります。点を別の順で覆う前進フロントは境界を先に敷いて中を埋めますが、Schönhardt のような穴は Steiner 点なしでは閉じません。面に壁関数を課すなら、作った点を式のまま持ち、その 2 個も含めて越えます。