Delaunay 三角剖分 — 空外接圆条件与 Bowyer

要保住空外接圆，只删坏三角形；符号靠近0时，交给精确谓词。

1934 年，Boris Delaunay 为纪念恩师 Georgy Voronoi，定义了一种三角剖分。规则只有一行：任何三角形的外接圆内部都不能有其他点。正是这一个条件，成了今天几乎所有非结构网格生成器的核心。

一个空外接圆就决定网格好不好。Bowyer–Watson 逐点插入，并保住这条性质。浮点数会悄悄背叛这条规则。

用滑块把顶点 B 向上拉，三角形 A–B–C 的外接圆变大。

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.

点 D 进入圆内的那一刻，对角线 A–C 变得非法，Delaunay 解翻成 B–D。两个 min-angle 里较大的那个，属于被选中的对角线。

由一个空圆定义的网格

把同一组点连成三角形的方法有无数种。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$ 为基准的相对坐标，再附上它到 $d$ 的距离平方。这个值的符号取决于三角形的方向(orientation)：顺时针时符号会翻转，所以实际代码会乘上三角形的方向符号来归一化。

前面把 B 拉高、D 进入圆内的那一刻，就是这个行列式改符号的地方。

Bowyer–Watson — 落点之后再把洞填回去#

现在来构建网格。覆盖同一组点的另一种做法是前进前沿。Bowyer–Watson 是一种逐点插入、每步都保持 Delaunay 性质的增量(incremental)算法，流程分五步。

从一个大到能包住所有点的超级三角形开始。
对新点 $p$ ，找出所有外接圆包含 $p$ 的三角形 — 称为坏三角形。
坏三角形总是构成一个星形(star-shaped)的洞。它的边界，是恰好属于一个坏三角形的那些边。
删掉所有坏三角形，把 $p$ 与边界的每条边相连，形成新三角形。
放完所有点后，移除任何与超级三角形共享顶点的三角形。

关键在第 3 步的性质。插入 $p$ 后违反空外接圆条件的三角形，只有「外接圆包含 $p$ 的三角形」，而它们总是构成一个连通的洞。只把那个洞填回去，整个网格就重新成为 Delaunay。

每加一个点就重织的画面，放在 60 个点的输出正下方。

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

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，没有任何黄点会落进另一个三角形的圆里。

三角形数量并非偶然。在平面 Delaunay 中，三角形个数固定为 $2N - 2 - h$ ( $N$ 是点数， $h$ 是凸包上的点数)。当 $N = 60$ 、凸包上有 8 个点时，恰好得到 110 个。网格还没画出来，就能用个数先验算一遍。

把 110 用手再数一遍。2 乘 60 是 120。120 减去 2 得 118，118 再减去凸包上的 8 个点，得到 110。漏掉凸包项就会得到 118，代码却印出 110，看起来像是 bug。少掉的 8 就是 h。种子 7 的单位正方形、60 个点，给出的就是这个 8。

in-circle 只在行列式乘上方向符号后大于 1e-12 时，才把点算作内部。四点共圆时行列式为 0，检验不会通过。两条对角线都不是内部，算法就留下先铺下的那条。0 和 1e-12 之间极小的正数也当成外部。正方网格的四个角总是落在这段里。

超级三角形的 r 是到重心最远距离再加 1。顶点在 (cx−3r, cy−r)、(cx+3r, cy−r)、(cx, cy+3r)。底边宽 6r，从下端到上端高 4r。样本点落在半径 r−1 的圆盘里，超级顶点大约在中心外 3r。60 个点在插入之前就已经在里面。最后删掉与超级顶点共享的三角形，剩下 110 个。

当点几乎都落在同一个圆上

优雅的规则有陷阱。当四个点恰好落在同一个圆上(cocircular)时，in-circle 行列式为 0。任选哪条对角线都是 Delaunay，剖分不唯一。最坏的情况是正方网格 — 每个正方形的四角都共圆，于是最规则的输入反而是最含糊的输入。

真正的麻烦出在那个 0 附近，由浮点数引起。当行列式达到机器精度极限，它的符号会随机翻转，拓扑(topology)陷入矛盾，代码要么卡进死循环，要么崩溃。所以 Triangle、CGAL 这类实战库不用普通 double，而用自适应精确谓词(adaptive exact predicates，Shewchuk)。它们只精确计算到所需位数，保证符号永远不出错。

到了三维，陷阱又多一个。空外接球条件并不禁止近乎扁平的四面体(sliver)。所以三维网格要用 Delaunay refinement(Ruppert、Chew)插入 Steiner 点来提升质量。

最后一个实务陷阱是边界。原始 Delaunay 网格并不尊重你画的几何边界。要强制保留某条边界边，就需要上面说的约束 Delaunay。「为什么网格直接穿墙而过」就是漏掉这一步的信号。

在正方网格上照单相信 in-circle 符号时#

四个角落在同一个圆上，行列式就是 0。1e-12 的门槛不把这个 0 当成内部。两条对角线都是 Delaunay，代码却只留下先出现的那条。整张正方网格都是这种平局，于是最规则的输入先从拓扑上晃起来。

用三角形个数先筛比较便宜。N=60 时得到 118，就是漏了凸包项 h。得到 110，才和这个种子的 h=8 相符。个数对了墙还是被穿透，该怀疑的是没做边界恢复，不是空外接圆。

Delaunay 三角剖分 — 空外接圆条件与 Bowyer–Watson 插点法