点群を三角形で覆う二つの哲学

対角線は外接円の符号でひっくり返り、フロントの次の三角形は形状と寸法の積が選ぶ。

1934年、ボリス・ドロネーは師ゲオルギー・ボロノイの名を冠した図から双対を一つ取り出しました。平面上の点を三角形で結び、どの外接円の中にも他の点が入らないようにする規則です。90年たちました。格子生成器は今もこの規則で動いています。

シード3で拾った点12個は、三角形17個になります。その出力を読む前に、四点が四角形になるとき対角線がどちらへ引かれるかを手で確認します。上側の頂点を上下にドラッグすると、共有辺を越えた瞬間に対角線が反転します。

Drag vertex T through the shared edge. When T enters the circumcircle of the lower triangle, the diagonal snaps from LR to BT — a single Lawson flip restoring the empty-circumcircle rule.

頂点Tが下側三角形の外接円に入る位置を見てください。対角線がLRからBTへ変わるその一回が、ローソン・フリップです。空の外接円を局所的に戻す最小の操作です。

同じ点集合の三角形分割は無数にあります。CFDが嫌うのは細く尖った三角形です。平たいと数値微分の誤差が膨らみます。線形系の条件数も一緒に悪化します。

Delaunay三角形分割は、可能な分割の中で最小内角を最大化します。いちばん尖った角を、点の配置が許す限り鈍くします。点の位置だけ与えれば答えは一つなので、同じ点をあとから結び直しても接続は揺れません。

前進フロントは順序が逆です。先に点を撒きません。境界線を出発点にして、内側へ三角形を一つずつ貼ります。新しい点はそのとき作ります。各瞬間で最良の要素を選ぶ局所的な貪欲法です。

空の外接円 — Delaunayを定義する性質#

規定する文は一つです。どの三角形の外接円も、内部に他の頂点を含みません。 空の外接円条件(empty circumcircle)です。点 $d$ が三角形 $abc$ の外接円の内部にあるかは、行列式一つで判定します。 $abc$ を反時計回りに置いて次を計算します。

\begin{vmatrix} a_x-d_x & a_y-d_y & (a_x-d_x)^2+(a_y-d_y)^2 \\ b_x-d_x & b_y-d_y & (b_x-d_x)^2+(b_y-d_y)^2 \\ c_x-d_x & c_y-d_y & (c_x-d_x)^2+(c_y-d_y)^2 \end{vmatrix} > 0

ここで $a,b,c$ は三角形の三頂点、 $d$ は検査する点です。行列式が正なら $d$ は円の内側、負なら外側、0なら円周上にちょうど乗ります。割り算も平方根もありません。符号だけ見ればよいので、浮動小数点誤差にも比較的安全です。この符号が挿入の一歩でどう崩れるかは、空の外接円とBowyer–Watsonの点挿入に別記してあります。

Bowyer–Watson: 空洞を掘り、再び埋める#

ローソン・フリップを繰り返してもDelaunay格子は作れます。よりまっすぐな道は、1981年にエイドリアン・バウアーとデイビッド・ワトソンが同じ学術誌へ並べて出した増分挿入です。

すでにDelaunayである格子に新しい点 $p$ を入れると、 $p$ を外接円の中に含む三角形はすべて無効になります。消すと多角形の穴、すなわち**空洞(cavity)**が残ります。空洞の境界辺それぞれと $p$ を結べば、穴は埋まります。

手順は三段階です。

新しい点 $p$ を外接円に含む三角形(悪い三角形)をすべて見つける。
悪い三角形が成す空洞の境界辺を取り出す。境界辺は悪い三角形一つだけに属する辺である。
悪い三角形を消し、境界辺ごとに $p$ を結んで三角形を新しく作る。

最初はすべての入力点を包む巨大なスーパー三角形一つから出発します。点をすべて入れたあと、スーパー三角形の頂点に掛かる三角形を取り除けば、最終的なDelaunay格子だけが残ります。

キャンバスをクリックして点を落とすと、挿入のたびに悪い三角形が消え、空洞が埋まり直されます。

Show circumcircles

Click on the canvas to insert a point. Each insertion rebuilds the Delaunay triangulation via Bowyer–Watson. Toggle the circumcircles: every one stays empty of other vertices — that is the Delaunay property.

外接円の表示をオンにしてください。どの円も他の頂点を含んでいないはずです。点をどれほど密に打っても、この性質は壊れません。それがDelaunayです。

Python — 増分Delaunayをゼロから書く#

外接円判定から挿入カーネルまでを、ライブラリなしの純粋なPythonへ移したコードです。

import numpy as np
 
def circumcircle(a, b, c):
    """三点の外接円の中心と半径の二乗を返す。一直線ならNone。"""
    ax, ay = a; bx, by = b; cx, cy = c
    d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    if abs(d) < 1e-12:
        return None
    a2, b2, c2 = ax * ax + ay * ay, bx * bx + by * by, cx * cx + cy * cy
    ux = (a2 * (by - cy) + b2 * (cy - ay) + c2 * (ay - by)) / d
    uy = (a2 * (cx - bx) + b2 * (ax - cx) + c2 * (bx - ax)) / d
    return (ux, uy), (ax - ux) ** 2 + (ay - uy) ** 2
 
def in_circumcircle(p, tri, pts):
    """点pが三角形triの外接円の内部ならTrue。"""
    cc = circumcircle(pts[tri[0]], pts[tri[1]], pts[tri[2]])
    if cc is None:
        return False
    (ux, uy), r2 = cc
    return (p[0] - ux) ** 2 + (p[1] - uy) ** 2 < r2 - 1e-12
 
def bowyer_watson(points):
    """増分挿入でDelaunay三角形分割を作る。"""
    pts = list(points)
    xs = [p[0] for p in pts]; ys = [p[1] for p in pts]
    dmax = 20 * max(max(xs) - min(xs), max(ys) - min(ys))
    cx, cy = (min(xs) + max(xs)) / 2, (min(ys) + max(ys)) / 2
    base = len(pts)                                    # スーパー三角形の頂点開始インデックス
    pts += [(cx - dmax, cy - dmax), (cx + dmax, cy - dmax), (cx, cy + dmax)]
    tris = [(base, base + 1, base + 2)]
 
    for pi in range(base):                             # 点を一つずつ挿入
        p = pts[pi]
        bad = [t for t in tris if in_circumcircle(p, t, pts)]
        edge_count = {}                                # 空洞境界 = 一度だけ現れた辺
        for t in bad:
            for e in [(t[0], t[1]), (t[1], t[2]), (t[2], t[0])]:
                key = tuple(sorted(e))
                edge_count[key] = edge_count.get(key, 0) + 1
        boundary = [e for e, n in edge_count.items() if n == 1]
        tris = [t for t in tris if t not in bad]
        tris += [(a, b, pi) for a, b in boundary]      # pを境界辺に結ぶ
 
    return [t for t in tris if all(i < base for i in t)]  # スーパー三角形を除去
 
if __name__ == "__main__":
    rng = np.random.default_rng(3)
    P = [tuple(xy) for xy in rng.random((12, 2))]
    T = bowyer_watson(P)
    print(f"点 {len(P)}個 -> 三角形 {len(T)}個")
    print("最初の三角形の頂点インデックス:", T[0])

実行すると点は12個、三角形は17個です。凸包の中の三角形の数はおよそ $2n - 2 - h$ 個です( $n$ は点の数、 $h$ は凸包上の点の数)。純粋な挿入カーネルは60行ほどに収まります。

12点をその式へ入れると $2\times 12-2=22$ です。出力が17なので、式がぴったり合うとき凸包上の点は $22-17=5$ 個です。スーパー三角形のスケールは点群の幅の20倍です。サンプルは単位正方形から引かれるので幅は1を超えず、三頂点はそのスケールだけ箱の中心から離れて、入力12点を内側に置きます。20は包含のための余裕であり、三角形が細いかどうかとは無関係です。

コードが置くスーパー三角形は、底辺がスケールの2倍、高さもスケールの2倍です。面積はスケールの二乗の2倍になります。点群の幅を1に置くとスケールは20、面積は800です。点が入っている正方形の面積1と比べると、スーパー三角形はその800倍です。最後にスーパー頂点に掛かる三角形を取り除くので、この800は結果に残りません。残るのは点12個から作った三角形17個だけです。

本文の行列式は符号だけを見ます。コードは、距離の二乗が半径の二乗から1e-12を引いた値より小さいときだけ内部と判定します。円周のすぐ上に乗った点は外側へ押し出されます。三点がほぼ一直線なら同じ1e-12で中心計算を放棄します。空洞側では、辺が何回現れたかだけを数えます。二回出た辺は穴の内側です。一回だけの辺が新しい三角形の底辺になり、その本数だけ三角形が増えます。

境界はどう守るか — 制約とパイプ

ここまでは点さえあればよい問題でした。実際のCFD形状は違います。翼型の表面やシリンダーの壁のように、必ず格子辺として存在しなければならない境界線があります。ところが純粋なDelaunayは、この境界辺を横切ってしまいます。

解法は制約Delaunay三角形分割(constrained DT)です。失われた境界線 $AB$ を戻すには、まず $AB$ が横切る三角形の帯、いわゆる**パイプ(pipe)**を見つけます。点 $A$ に付いた三角形から出発し、 $AB$ が貫く隣接三角形をたどって $B$ に達すれば、パイプが完成します。

パイプが決まれば、その内部を再三角形分割して $AB$ を辺として復元します。対角線の交換(swapping of diagonals)や分割統治で処理します。境界は戻りますが、その近くの三角形はもはや完全なDelaunayではないことがあります。形状の忠実さと格子品質を引き換えにします。

前進フロント — 貪欲に要素を積む

前進フロント法(AFT)は、この境界問題をひっくり返して解きます。境界がそのまま出発線です。

境界セグメントの一覧から基準セグメント $AB$ を一つ選びます。この $AB$ と最良の三角形を成す点 $C$ を探します。候補は二種類です。既存フロント上の点を使うか、新しい内部点 $I$ を作るか。品質の高い方を選んで三角形を作り、フロントを更新します。新しくできた辺はフロントに追加し、すでにフロントにあった辺は閉じたものとみなして除去します。

品質は二つの因子の積で測ります。

\lambda = \alpha \cdot \delta

ここで $\alpha$ は三角形の形状因子(正三角形なら1、平たいほど0に近づく)、 $\delta$ は要求された要素サイズと実際の辺長の一致度(両者が同じなら1)です。 $\lambda$ が大きいほど良い要素です。

形状因子はこう書けます。

import math
 
def shape_quality(a, b, c):
    """形状因子alpha。正三角形なら1、平たいほど0に近づく。"""
    l1 = math.dist(a, b); l2 = math.dist(b, c); l3 = math.dist(c, a)
    area = abs((b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])) / 2
    denom = l1 * l1 + l2 * l2 + l3 * l3
    return 4 * math.sqrt(3) * area / denom if denom > 0 else 0.0

フロントが空になると(閉じた境界が内側へすべて集まると)格子が完成します。各段階で局所最適を選ぶため、AFTの初期格子品質はたいていDelaunayより高くなります。特に境界層のように方向性の強い異方性格子で、差が大きくなります。

一辺の長さが等しい正三角形では、面積は辺の二乗に $\sqrt{3}/4$ を掛けた値で、三辺の二乗和は辺の二乗の三倍です。コードの $4\sqrt{3}$ を面積に掛けてその和で割ると、 $4\sqrt{3}\cdot(\sqrt{3}/4)/(3)=1$ になります。一辺が伸びて高さが落ちると分子は減り分母は増えるので、 $\alpha$ は0へ向かいます。寸法の一致 $\delta$ が1でも、積 $\lambda$ は一緒に小さくなります。平たい要素がフロントに残らない理由です。

代償は速度です。新しい要素ごとに、フロントの他の辺と交差するかを検査しなければなりません。素朴に書くと要素一つで $O(n)$ 、全体が $O(n^2)$ に膨らみ得ます。背景格子や四分木(quadtree)でこの検査を局所化しないと、フロントが長くなるほど交差検査が生成を食い潰します。

いつ何を選ぶか

二つの方法は競争というより分業です。

Delaunayは点分布が与えられたとき、一意で速い答えを返します。凸包の充填、散乱データの補間、事後の格子改善に向きます。一方、境界を自力では守れず、制約処理が別に必要です。

AFTは境界忠実度と初期要素品質に優れ、異方性への拡張も容易です。代わりに実装が難しく、データ構造の設計で性能が大きく左右されます。

そこで現代の格子生成器は両者を混ぜます。配置の順序はAFT、点の連結はDelaunay。フロント–Delaunayがその分担です。

翼の辺が切れたら挿入だけでは終わらせない

点12個が三角形17個で終わるのは、点がすでに領域の中にあるときの話です。翼面が対角線で切られると、その17個の一部は形状を代表しません。その帯をパイプとして再分割するか、切れた辺をフロントの基準セグメントとして残し、次の点 $C$ をそれに対して選ぶ必要があります。

$\lambda = \alpha\delta$ が大きい三角形を積むだけでは、切れた境界は戻りません。制約辺を戻すために点をどこへ足すかは、Steiner点で切れた境界を戻す手順に書いてあります。挿入カーネルの60行を先に合わせ、境界が切れる形状のときだけ、その手順へ移ります。

点群を三角形で覆う二つの哲学 — Bowyer–Watsonと前進フロント