ドロネー三角形分割 — 空の外接円条件とBowyer–Watson点挿入

1934年、Boris Delaunayは恩師Georgy Voronoiに敬意を表して、ある三角形分割を定義しました。規則はたった一行です。どの三角形の外接円の内部にも、ほかの点が入ってはいけない。この一つの条件が、今日のほぼすべての非構造格子生成器の心臓になりました。

この記事では、その「空の外接円」条件がなぜ良い格子を生むのかを見ます。そして点を一つずつ挿入して格子を育てるBowyer–Watsonアルゴリズムを、Python 40行で実装します。最後に、浮動小数点がこの優雅な規則をどう静かに裏切るかまで扱います。

一つの空円で定義される格子#

同じ点集合を三角形でつなぐ方法は無数にあります。Delaunayはその中のただ一つを選びます。すべての三角形の外接円(三頂点を通る円)が空である分割です。

三角形 $T = \{a, b, c\}$ の外接円を $C(T)$ とすると、条件はこう書けます。

\forall\, p \notin T : \quad p \notin \mathrm{int}\, C(T)

ここで $\mathrm{int}\,C(T)$ は外接円の内部、 $p$ は格子のほかのすべての点です。どの点もほかの三角形の外接円の内側にあってはいけません。

なぜこの規則なのでしょうか。Delaunay分割は、同じ点集合のあらゆる三角形分割の中で最小内角を最大化します。細長い三角形(sliver)を可能なかぎり避けるという意味です。細い三角形は、有限体積法・有限要素法で勾配再構成と行列の条件数を壊します。だから格子生成器はDelaunayを既定値にします。

おまけに、Delaunay三角形分割の双対(dual)がそのままVoronoi図です。一方の構造が得られれば、もう一方が無料でついてきます。

InCircle行列式 — 円の内外を一発で判定#

「点が外接円の内側にあるか」は、外接円の中心や半径を求めずに一発で答えられます。反時計回りに並んだ三角形 $(a, b, c)$ に対し、点 $d$ が外接円の内部にある必要十分条件は、次の行列式が正であることです。

\det \begin{pmatrix} 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{pmatrix} > 0

各行は、一つの頂点を $d$ 基準の相対座標に移し、その座標と距離の二乗を並べたものです。この値の符号は三角形の向き(orientation)に依存します。時計回りだと符号が反転するので、実際のコードでは三角形の向きの符号を掛けて正規化します。

下の図で、スライダーを使って頂点Bを上に引き上げてみましょう。

vertex B height (drag B up toward the diagonal)

The circle is the circumcircle of triangle A–B–C. Slide B upward: when D crosses into the circle the diagonal A–C becomes illegal and the Delaunay solution flips to B–D. Notice the larger of the two min-angle numbers always belongs to the chosen diagonal — Delaunay maximizes the smallest angle.

Bが上がるほど、三角形A–B–Cの外接円が大きくなります。点Dが円の内側に入った瞬間、対角線A–Cは不正になり、Delaunay解はB–Dに反転します。二つのmin-angleの数値のうち、大きいほうが常に選ばれた対角線に付くことも確認してみてください。

Bowyer–Watson — 点を落として穴を埋め直す#

いよいよ格子を作ります。Bowyer–Watsonは、点を一つずつ挿入しながら毎ステップでDelaunay性を保つ逐次(incremental)アルゴリズムです。流れは五段階です。

すべての点を包めるほど大きなスーパー三角形から始めます。
新しい点 $p$ について、外接円が $p$ を含む三角形をすべて探します — これらを悪い三角形と呼びます。
悪い三角形は常に一つの星形(star-shaped)の穴を作ります。その境界は、悪い三角形ちょうど一つだけに属する辺の集合です。
悪い三角形をすべて消し、 $p$ を境界の各辺とつないで新しい三角形を作ります。
すべての点を入れたあと、スーパー三角形の頂点を共有する三角形を取り除きます。

要は3の性質です。 $p$ の挿入で空の外接円条件を破る三角形は「外接円が $p$ を含む三角形」だけで、それらは常に連結した穴を作ります。その穴だけを埋め直せば、全体がふたたびDelaunayになります。

下のキャンバスをクリックして、点を自分で追加してみましょう。

show circumcircles

Each thin cyan ring is the circumcircle of one triangle. In a valid Delaunay mesh no vertex (yellow dot) ever sits inside another triangle's ring — that is the empty-circumcircle property. Add points and watch the mesh repair itself.

点を追加するたびに格子が自分で組み直されます。細いシアン色の円が各三角形の外接円です。show circumcirclesをオンにしておくと、どの黄色い点もほかの三角形の円に入らないという空の外接円性が、目で見えます。

Lawsonの辺反転 — 対角線一本の局所的な選択#

Bowyer–Watsonが穴をまるごと埋め直すのに対し、Lawsonの方法は辺を一本ずつ直します。新しい点を含む三角形を三つに分けたあと、影響を受けた各辺にin-circle検査をかけます。反対側の点が外接円の内側にあれば、その共有辺(対角線)を反転します。不正な辺がなくなるまで繰り返します。

一度の反転は、関わる二つの三角形の最小内角を必ず大きくします。だからループが終われば、最小内角が最大化されたDelaunay格子に自動的になります。上の図で見たA–C ↔ B–Dの切り替えが、まさにこの一度の反転です。

境界線を強制的に残す必要があるとき(constrained Delaunay)も、同じ道具を使います。復元したい辺ABが横切る三角形の帯(pipe)を見つけ、その中を対角線交換で分割し直せば、ABが格子辺として復元されます。

Python — ランダムな60点をDelaunayに#

逐次挿入の全体をnumpyに移します。外部ライブラリなしで動きます。

import numpy as np
 
def in_circle(a, b, c, d):
    """点 d が三角形 a,b,c の外接円に厳密に入れば True。
    三角形の向きの符号で正規化し、winding に依存しないようにする。"""
    ax, ay = a[0] - d[0], a[1] - d[1]
    bx, by = b[0] - d[0], b[1] - d[1]
    cx, cy = c[0] - d[0], c[1] - d[1]
    det = np.linalg.det(np.array([
        [ax, ay, ax * ax + ay * ay],
        [bx, by, bx * bx + by * by],
        [cx, cy, cx * cx + cy * cy],
    ]))
    orient = (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])
    return det * np.sign(orient) > 1e-12
 
def super_triangle(points):
    pts = np.asarray(points)
    cx, cy = pts.mean(axis=0)
    r = np.max(np.linalg.norm(pts - [cx, cy], axis=1)) + 1.0
    return [(cx - 3 * r, cy - r), (cx + 3 * r, cy - r), (cx, cy + 3 * r)]
 
def bowyer_watson(points):
    st = super_triangle(points)
    tris = [(st[0], st[1], st[2])]
    for p in points:
        bad = [t for t in tris if in_circle(t[0], t[1], t[2], p)]
        edge_n = {}                              # 各辺の出現回数を集計
        for t in bad:
            for e in [(t[0], t[1]), (t[1], t[2]), (t[2], t[0])]:
                key = frozenset((e[0], e[1]))
                edge_n[key] = edge_n.get(key, 0) + 1
        boundary = [tuple(k) for k, n in edge_n.items() if n == 1]
        tris = [t for t in tris if t not in bad]
        for e in boundary:                       # 点 p で穴を埋め直す
            tris.append((e[0], e[1], p))
    s = set(st)
    return [t for t in tris if not (set(t) & s)]
 
rng = np.random.default_rng(7)
pts = [tuple(p) for p in rng.random((60, 2))]
mesh = bowyer_watson(pts)
print(f"points = {len(pts)}   triangles = {len(mesh)}")

出力はおおよそ次の通りです。

points = 60   triangles = 110

三角形の数は偶然ではありません。平面Delaunayでは三角形の個数は $2N - 2 - h$ で決まります( $N$ は点数、 $h$ は凸包上の点数)。 $N = 60$ で凸包上の点が8個なら、ちょうど110個になります。格子が描かれる前に、個数だけで検算できます。

点がほぼ一つの円の上に並ぶとき#

優雅な規則には罠があります。四点がちょうど一つの円の上にある(cocircular)とき、in-circle行列式は0になります。どちらの対角線を選んでもDelaunayで、分割は一意ではありません。最悪なのは正方格子です — すべての正方形の四隅がcocircularなので、最も規則的な入力が最も曖昧な入力になります。

本当の問題は、その0の近くで浮動小数点が起こします。行列式がマシン精度の限界に達すると、符号がランダムに反転し、位相(topology)が矛盾に陥って、コードが無限ループに入るか落ちます。だからTriangleやCGALのような実戦ライブラリは、単なる double の代わりに適応的厳密述語(adaptive exact predicates、Shewchuk)を使います。必要なだけの桁を正確に計算し、符号だけは絶対に間違えないようにします。

三次元に行くと、罠がもう一つ増えます。空の外接球条件は、ほぼ平らな四面体(sliver)を防げません。だから3DメッシュではDelaunay refinement(Ruppert・Chew)でSteiner点を追加し、品質を上げます。

最後の実務的な罠は境界です。生のDelaunay格子は、あなたが描いた形状境界を尊重しません。境界辺を強制的に残すには、上で見たconstrained Delaunayが別途必要です。「なぜ格子が壁を突き抜けるのか」は、これを抜かしたという合図です。

格子を敷く前に書き留める三行#

Delaunay = 空の外接円 = 最小内角の最大化。三つの言い方は同じもので、すべてsliverを避けるという同じ目的です。
Bowyer–Watsonは「悪い三角形を消して、穴を点で埋め直す」の五行で終わります。鍵は、穴が常に星形だという定理です。
cocircularな入力と浮動小数点はin-circleの符号を壊します。正方格子が疑わしいなら、まず厳密述語を用意してください。