用三角形覆盖点云的两种哲学

对角线由外接圆的符号翻转，阵面上下一个三角形由形状与尺寸的乘积来决定。

1934年，鲍里斯·德劳内从以老师格奥尔基·沃罗诺伊命名的图里抽出一个对偶。把点连成三角形，让任何一个外接圆里都不再有别的点。九十年过去了。网格生成器还在用这条规则运转。

种子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$ 在圆内，为负表示在圆外，为零表示恰好落在圆周上。没有除法，没有平方根。只看符号，因此对浮点误差相对稳健。这个符号在一次插入里怎么被打破，写在空外接圆与Bowyer–Watson插点里。

Bowyer–Watson：挖出空腔，再填回去#

反复做劳森翻转也能得到Delaunay网格。更直的一条路，是阿德里安·鲍耶与戴维·沃森1981年在同一份期刊上并排发表的增量插入。

把新点 $p$ 丢进已经是Delaunay的网格，凡是外接圆包含 $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是为了包住点的余量，和三角形尖不尖无关。

代码里的超级三角形，底边是尺度的两倍，高也是尺度的两倍，所以面积是尺度平方的两倍。把点云跨度取成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与推进阵面