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

1934年、ソ連の数学者ボリス・ドロネーは、恩師ゲオルギー・ボロノイの名を冠した図から一つの双対を取り出しました。平面に散らばった点を三角形で結び、どの三角形の外接円の中にも他の点が入らないようにする方法です。それから90年、この一つの規則はほぼすべてのCFD格子生成器の心臓部に埋め込まれています。

この記事では、非構造三角形格子を作る二つの哲学を追います。一つは点を一つずつ挿入する Bowyer–Watson 増分Delaunayアルゴリズム。もう一つは境界から内側へ要素を貪欲に積み上げる 前進フロント法(AFT) です。外接円判定を純粋なPythonで書き、格子が育つ様子を直接見て、二つの方法をいつ使い分けるかまで整理します。

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

同じ点集合を三角形で結ぶ方法は無数にあります。問題は「どの三角形分割が良いか」です。CFDにおける良い格子とは、細く尖った三角形が少ない格子です。平たい三角形は数値微分の誤差を増やし、線形系の条件数を悪化させます。

Delaunay三角形分割はこの要求に正面から答えます。可能な分割の中で 最小内角を最大化 します。つまり、格子内で最も尖った三角形の角度をできるだけ大きくします。点の位置だけ与えれば答えが一意に決まるのも魅力です。

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

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

Delaunay三角形分割を規定する文はただ一つです。どの三角形の外接円も内部に他の頂点を含まない。 これを空の外接円条件(empty circumcircle)と呼びます。

判定は行列式一つで終わります。点 $d$ が三角形 $abc$ の外接円の内部にあるか調べるには、$abc$ を反時計回りに置いて次を計算します。

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

四点が一つの四角形を成すとき、対角線をどちら向きに引くかはこの判定一つで決まります。下で上側の頂点を上下に引いてみてください。共有辺を越えた瞬間に対角線が反転します。

頂点 T を下の三角形の外接円の中へ押し込むと、対角線が LR から BT へ変わります。この一度の反転がまさにローソン・フリップ(Lawson flip)です。局所的に空の外接円条件を回復する最小の演算です。

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

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

発想は単純です。すでにDelaunayである格子に新しい点 $p$ を入れるとき、$p$ を外接円の中に含む三角形はすべて無効になります。この無効な三角形を消すと、多角形状の空き、すなわち 空洞(cavity) が生まれます。あとは空洞の境界辺それぞれと $p$ を結んで新しい三角形を埋めるだけです。

手順は三段階に要約されます。

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

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

下のキャンバスをクリックして点を一つずつ打ってみてください。挿入のたびに格子がBowyer–Watsonで組み直されます。

外接円の表示をオンにすると、どの円も他の頂点を含みません。点をどれほど密に打ってもこの性質は壊れません。それが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行ほどに収まります。

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

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

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

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

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

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

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

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

ここで $\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より高くなります。特に境界層のように方向性の強い異方性格子を作るとき強みが大きいです。

代償は速度です。新しい要素ごとにフロントの他の辺と交差するか検査しなければなりません。素朴に書くと要素一つで $O(n)$、全体が $O(n^2)$ に膨らみ得ます。元文書が強調するように、背景格子や四分木(quadtree)のような空間探索データ構造でこの検査を局所化することが実戦の核心です。

いつ何を選ぶか

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

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

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

そこで現代の格子生成器は両者を混ぜます。AFTで要素の配置順序を決めつつDelaunay連結基準で点を結ぶフロント–Delaunay方式が代表的です。二つの哲学の長所だけを取る折衷案です。

覚えておくこと

• 空の外接円条件 がDelaunayのすべてです。判定は符号だけ見る行列式一つで終わり、Bowyer–Watsonは空洞を掘って再び埋める60行のカーネルで実装されます。
• AFT は境界から内側へ $\lambda = \alpha\delta$ が最大の要素を貪欲に積みます。境界忠実度と異方性に強いですが、交差検査のコストがボトルネックです。
• 実務の格子生成器は両者を フロント–Delaunay で結合し、順序はAFT、連結はDelaunayに分担させます。