Voronoi 图与 Delaunay 三角剖分：对偶的力量

在计算几何中，有一对结构如同硬币的正反两面，它们彼此依存、相互转化——这就是 Voronoi 图和 Delaunay 三角剖分。Voronoi 图回答”离哪个点最近”这一根本问题，而 Delaunay 三角剖分则给出了”最优三角网格”这一工程核心需求。

我第一次接触这对结构是在做游戏地图生成的时候。当时需要把一个平面随机划分为若干区域，每个区域围绕一个”种子点”自然生长。用 Voronoi 图实现后，地图的效果令人惊艳——每个区域就像自然界的晶体结构一样规则而又有机。后来在做有限元分析时，又不得不面对 Delaunay 三角剖分。这两个概念看似各自独立，实际上通过对偶性紧密关联。

本文将从几何直觉出发，深入剖析这对结构的定义、性质、算法和工程实践。

一、Voronoi 图：空间的自然划分

1.1 定义与直觉

给定平面上 n 个不同的点 \(S = \{p_1, p_2, \ldots, p_n\}\)（称为生成点或站点），Voronoi 图将平面划分为 n 个区域，每个区域 \(V(p_i)\) 包含所有离 \(p_i\) 比离其他任何点都近的点：

\[V(p_i) = \{x \in \mathbb{R}^2 \mid d(x, p_i) \leq d(x, p_j), \forall j \neq i\}\]

其中 \(d\) 通常指欧几里得距离。

直觉上，可以把每个生成点想象成一个”领地中心”。如果所有领地同时以相同速度向外扩张，两个领地相遇的边界就是 Voronoi 边。这也正是”晶体生长模型”的几何抽象。

1.2 基本性质

Voronoi 图具有以下关键性质：

（一）每个 Voronoi 区域是凸多边形。 因为每个区域是若干半平面的交集，半平面的交集总是凸集。对于边界上的站点，区域可能是无界的凸多边形。

（二）Voronoi 边是两个站点的垂直平分线的一部分。 如果两个 Voronoi 区域 \(V(p_i)\) 和 \(V(p_j)\) 共享一条边，那么这条边一定落在 \(p_i\) 和 \(p_j\) 的垂直平分线上。

（三）Voronoi 顶点是三个（或更多）站点的等距点。 每个 Voronoi 顶点到产生它的三个站点距离相等，即它是这三个站点外接圆的圆心。

（四）复杂度。 对于 n 个站点，Voronoi 图最多有 \(2n - 5\) 个顶点和 \(3n - 6\) 条边（由欧拉公式可得）。这意味着 Voronoi 图的规模是 \(O(n)\) 的。

1.3 退化情况

在实际工程中，退化情况是必须处理的：

共线点：如果所有站点共线，Voronoi 图退化为一系列平行线，没有有限顶点。
共圆点：如果四个或更多站点共圆，对应的 Voronoi 顶点度数大于 3，此时需要特殊处理。
重合点：如果两个站点重合，Voronoi 图未定义。在工程中通常需要预处理去除重复点。

1.4 形式化：半平面交

从算法角度看，每个 Voronoi 区域可以表示为半平面交：

\[V(p_i) = \bigcap_{j \neq i} H(p_i, p_j)\]

其中 \(H(p_i, p_j) = \{x \mid d(x, p_i) \leq d(x, p_j)\}\) 是包含 \(p_i\) 的半平面。

朴素地对每个站点计算 \(n - 1\) 个半平面的交集，单个区域需要 \(O(n \log n)\) 时间，总计 \(O(n^2 \log n)\)。这远不是最优的。

# 半平面交的概念演示（非实际使用）
def voronoi_cell_naive(sites, index):
    """计算第 index 个站点的 Voronoi 区域（半平面交）"""
    pi = sites[index]
    halfplanes = []
    for j, pj in enumerate(sites):
        if j == index:
            continue
        # 垂直平分线的法向量和偏移
        mid = ((pi[0] + pj[0]) / 2, (pi[1] + pj[1]) / 2)
        normal = (pj[0] - pi[0], pj[1] - pi[1])
        halfplanes.append((mid, normal))
    return intersect_halfplanes(halfplanes)

二、Delaunay 三角剖分：最优三角网格

2.1 定义

给定平面上的点集 \(S\)，Delaunay 三角剖分（Delaunay Triangulation，简称 DT）是 \(S\) 的一个三角剖分，满足空圆性质：任何三角形的外接圆内部不包含 \(S\) 中的其他点。

更精确地说，对于 DT 中的每个三角形 \(\triangle p_i p_j p_k\)，其外接圆的开圆盘内不包含 \(S\) 中的任何点。

2.2 核心性质

（一）最大化最小角性质。 在所有可能的三角剖分中，Delaunay 三角剖分最大化了所有三角形角度中的最小角。这意味着 Delaunay 三角剖分尽可能避免出现”瘦长”三角形，这对数值计算至关重要。

（二）唯一性。 当没有四点共圆时，Delaunay 三角剖分是唯一的。如果存在四点共圆的情况，可以有多个合法的 Delaunay 三角剖分（通过不同的对角线翻转得到）。

（三）局部可验证性。 一个三角剖分是 Delaunay 的，当且仅当对于每条内部边，其相邻的两个三角形满足空圆性质。这意味着我们可以通过局部边翻转将任意三角剖分转化为 Delaunay 三角剖分。

（四）包含最近邻图。 如果 \(p_j\) 是 \(p_i\) 的最近邻，那么边 \(p_i p_j\) 一定出现在 Delaunay 三角剖分中。更进一步，DT 包含最小生成树（MST）、相对邻域图（RNG）和 Gabriel 图。

2.3 空圆判定

判断点 \(p_d\) 是否在三角形 \(\triangle p_a p_b p_c\) 的外接圆内，可以通过计算以下行列式的符号：

\[ \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 \]

当 \(p_a, p_b, p_c\) 按逆时针排列时，行列式值大于 0 表示 \(p_d\) 在外接圆内。

/* 外接圆判定（InCircle 测试）
 * 返回值 > 0：pd 在三角形 pa-pb-pc 的外接圆内
 * 返回值 = 0：pd 在外接圆上
 * 返回值 < 0：pd 在外接圆外
 * 要求 pa, pb, pc 按逆时针排列
 */
double incircle(double ax, double ay,
                double bx, double by,
                double cx, double cy,
                double dx, double dy)
{
    double adx = ax - dx, ady = ay - dy;
    double bdx = bx - dx, bdy = by - dy;
    double cdx = cx - dx, cdy = cy - dy;

    double ab_det = adx * bdy - bdx * ady;
    double bc_det = bdx * cdy - cdx * bdy;
    double ca_det = cdx * ady - adx * cdy;

    double a_lift = adx * adx + ady * ady;
    double b_lift = bdx * bdx + bdy * bdy;
    double c_lift = cdx * cdx + cdy * cdy;

    return a_lift * bc_det + b_lift * ca_det + c_lift * ab_det;
}

2.4 方向判定

三角剖分中另一个基本谓词是方向测试（orientation test）：

/* 方向判定
 * 返回值 > 0：pa, pb, pc 逆时针排列
 * 返回值 = 0：三点共线
 * 返回值 < 0：pa, pb, pc 顺时针排列
 */
double orient2d(double ax, double ay,
                double bx, double by,
                double cx, double cy)
{
    return (bx - ax) * (cy - ay) - (by - ay) * (cx - ax);
}

这两个谓词是所有 Delaunay 三角剖分算法的基石。在工程中，浮点误差可能导致错误的判定结果，后面会详细讨论这个问题。

三、对偶性：Voronoi 与 Delaunay 的桥梁

3.1 对偶关系

Voronoi 图和 Delaunay 三角剖分互为对偶图，这是本文的核心观点。对偶关系体现在以下几个层面：

Voronoi 图	Delaunay 三角剖分	对应关系
Voronoi 区域 \(V(p_i)\)	Delaunay 顶点 \(p_i\)	面 ↔︎ 顶点
Voronoi 边	Delaunay 边	边 ↔︎ 边
Voronoi 顶点	Delaunay 三角形	顶点 ↔︎ 面

（一）顶点-面对偶。 Voronoi 图的每个顶点对应 Delaunay 三角剖分中的一个三角形。Voronoi 顶点是该三角形外接圆的圆心。

（二）边-边对偶。 Voronoi 图的每条边垂直于对应的 Delaunay 边，且两者交于中点。换句话说，如果 Voronoi 边分隔区域 \(V(p_i)\) 和 \(V(p_j)\)，那么 Delaunay 三角剖分中存在边 \(p_i p_j\)，并且这两条边互相垂直。

（三）面-顶点对偶。 Voronoi 图的每个区域 \(V(p_i)\) 对应 Delaunay 三角剖分中的顶点 \(p_i\)。\(V(p_i)\) 的边数等于 \(p_i\) 在 Delaunay 三角剖分中的度数。

3.2 从 Voronoi 到 Delaunay

有了 Voronoi 图，构造 Delaunay 三角剖分非常简单：对于每对共享 Voronoi 边的相邻区域 \(V(p_i)\) 和 \(V(p_j)\)，连接 \(p_i\) 和 \(p_j\) 即可。

def voronoi_to_delaunay(voronoi):
    """从 Voronoi 图构造 Delaunay 三角剖分"""
    delaunay_edges = set()
    for edge in voronoi.edges:
        # 每条 Voronoi 边分隔两个区域
        pi = edge.left_site
        pj = edge.right_site
        if pi is not None and pj is not None:
            delaunay_edges.add((min(pi, pj), max(pi, pj)))
    return delaunay_edges

3.3 从 Delaunay 到 Voronoi

反过来也一样直接：对于 Delaunay 三角剖分中的每个三角形，计算其外接圆圆心，这些圆心就是 Voronoi 顶点。连接共享边的相邻三角形的外接圆圆心，就得到 Voronoi 边。

def delaunay_to_voronoi(triangulation):
    """从 Delaunay 三角剖分构造 Voronoi 图"""
    voronoi_vertices = {}
    voronoi_edges = []

    for tri in triangulation.triangles:
        # 外接圆圆心即为 Voronoi 顶点
        center = circumcenter(tri.a, tri.b, tri.c)
        voronoi_vertices[tri.id] = center

    for edge in triangulation.interior_edges:
        tri1 = edge.left_triangle
        tri2 = edge.right_triangle
        v1 = voronoi_vertices[tri1.id]
        v2 = voronoi_vertices[tri2.id]
        voronoi_edges.append((v1, v2))

    return voronoi_vertices, voronoi_edges

3.4 对偶性的意义

对偶性的实际意义在于：我们只需要一个高效算法。 只要能高效构造 Voronoi 图或 Delaunay 三角剖分其中之一，就可以在线性时间内转化为另一个。工程中通常先构造 Delaunay 三角剖分（因为增量算法更容易实现），再转化为 Voronoi 图。

四、Fortune 扫描线算法

4.1 算法思想

Steven Fortune 在 1986 年提出了一个优雅的 \(O(n \log n)\) 算法来构造 Voronoi 图。传统扫描线算法维护一条水平或垂直的扫描线，但 Fortune 算法的巧妙之处在于引入了海滩线（beach line）这一概念。

关键洞察：扫描线将平面分为”已处理”和”未处理”两部分。扫描线上方的站点已经被处理，其 Voronoi 区域的下半部分由海滩线描述。海滩线是若干抛物线弧的下包络线。

4.2 海滩线

每个已处理的站点 \(p_i\) 贡献一段抛物线弧。以站点 \(p_i\) 为焦点、扫描线为准线的抛物线上的点到 \(p_i\) 和扫描线的距离相等。海滩线是所有这些抛物线的下包络线——即在每个 x 坐标处取 y 值最小的抛物线弧。

海滩线上相邻抛物线弧的交点称为断点（breakpoint），断点的轨迹就是 Voronoi 边。

4.3 事件类型

Fortune 算法处理两种事件：

站点事件（Site Event）：当扫描线经过一个新的站点时触发。新站点会在海滩线上产生一个新的抛物线弧，将原有的某段弧一分为二。

圆事件（Circle Event）：当海滩线上某段抛物线弧缩小到零长度时触发。此时，该弧两侧的断点汇合，对应一个 Voronoi 顶点的生成。圆事件之所以叫”圆事件”，是因为三段相邻弧对应的三个站点共圆，圆的最低点就是事件触发的位置。

4.4 算法流程

Fortune-Voronoi(S):
  1. 初始化事件队列 Q（按 y 坐标排序），将所有站点作为站点事件加入 Q
  2. 初始化海滩线 T（平衡二叉搜索树，叶子是抛物线弧，内部节点是断点）
  3. while Q 非空:
       取出 Q 中 y 坐标最小的事件 e
       if e 是站点事件:
         处理站点事件(e, T, Q)
       else:
         处理圆事件(e, T, Q)
  4. 处理剩余的海滩线弧（无限远的 Voronoi 边）

4.5 站点事件处理

处理站点事件(p, T, Q):
  1. 在 T 中找到 p 正上方的抛物线弧 alpha
  2. 如果 alpha 有关联的圆事件，将其标记为无效（false alarm）
  3. 用三段弧替换 alpha：<..., alpha_l, p, alpha_r, ...>
     其中 alpha_l 和 alpha_r 是 alpha 对应站点的弧
  4. 创建两条新的半边（未来的 Voronoi 边）
  5. 检查新的三元组是否可能产生圆事件：
     - 检查 alpha_l 的左邻居、alpha_l、p 三段弧
     - 检查 p、alpha_r、alpha_r 的右邻居三段弧

4.6 圆事件处理

处理圆事件(gamma, T, Q):
  1. 令 alpha 是 gamma 左侧的弧，beta 是右侧的弧
  2. 从 T 中删除 gamma
  3. 删除 alpha 和 beta 关联的旧圆事件（如果有）
  4. 计算 Voronoi 顶点（三个站点的外接圆圆心）
  5. 完成两条 Voronoi 半边，开始一条新的半边
  6. 检查新形成的三元组是否可能产生新的圆事件

4.7 复杂度分析

时间复杂度：\(O(n \log n)\)。事件队列使用优先队列，每次操作 \(O(\log n)\)。事件总数为 \(O(n)\)（站点事件 \(n\) 个，圆事件最多 \(2n - 5\) 个）。
空间复杂度：\(O(n)\)。海滩线和事件队列的规模都是 \(O(n)\)。

Fortune 算法是理论最优的（\(\Omega(n \log n)\) 下界可以通过排序问题的规约证明），但实现复杂度较高。工程中很多项目选择实现 Delaunay 三角剖分（增量算法更简单），再转化为 Voronoi 图。

# Fortune 算法的骨架伪代码（Python 风格）
class FortuneAlgorithm:
    def __init__(self, sites):
        self.event_queue = PriorityQueue()  # 按 y 坐标排序
        self.beach_line = BalancedBST()     # 海滩线
        self.voronoi = VoronoiDiagram()
        for site in sites:
            self.event_queue.push(SiteEvent(site))

    def run(self):
        while not self.event_queue.empty():
            event = self.event_queue.pop()
            if isinstance(event, SiteEvent):
                self.handle_site_event(event)
            elif isinstance(event, CircleEvent) and event.valid:
                self.handle_circle_event(event)
        self.finalize_edges()

    def handle_site_event(self, event):
        site = event.site
        arc_above = self.beach_line.find_arc_above(site.x)
        if arc_above.circle_event:
            arc_above.circle_event.valid = False

        new_arc = Arc(site)
        self.beach_line.replace_arc(arc_above, new_arc)
        self.voronoi.create_half_edges(arc_above.site, site)
        self.check_circle_event(new_arc.prev)
        self.check_circle_event(new_arc.next)

    def handle_circle_event(self, event):
        arc = event.arc
        left = arc.prev
        right = arc.next
        vertex = circumcenter(left.site, arc.site, right.site)

        self.voronoi.add_vertex(vertex)
        self.voronoi.finish_edges(arc, vertex)

        self.beach_line.remove(arc)

        if left.circle_event:
            left.circle_event.valid = False
        if right.circle_event:
            right.circle_event.valid = False

        self.check_circle_event(left)
        self.check_circle_event(right)

五、增量构造 Delaunay 三角剖分

5.1 基本思想

增量构造算法是最直观也最常用的 Delaunay 三角剖分算法。其思想非常简单：从一个包含所有点的超级三角形开始，逐个插入点，每次插入后通过边翻转维护 Delaunay 性质。

5.2 算法步骤

Incremental-Delaunay(S):
  1. 创建一个足够大的超级三角形 T0，包含 S 中的所有点
  2. for each 点 p in S:
       a. 找到包含 p 的三角形 t
       b. 将 t 分裂为三个子三角形（连接 p 到 t 的三个顶点）
       c. 对新产生的边执行 edge flipping，直到满足 Delaunay 性质
  3. 删除与超级三角形顶点关联的所有三角形

5.3 点定位

步骤 2a 中的点定位是性能的关键。几种常见策略：

遍历搜索：从任意三角形开始，沿着”靠近目标点”的方向遍历。平均 \(O(\sqrt{n})\)，最坏 \(O(n)\)。
跳跃-行走（jump-and-walk）：随机选择若干三角形，从最近的开始行走。期望 \(O(n^{1/3})\)。
历史 DAG：记录三角形分裂的历史，形成一棵有向无环图。每次定位从根节点开始向下搜索。期望 \(O(\log n)\)，但需要额外空间。

5.4 边翻转（Edge Flipping）

边翻转是增量算法的核心操作。对于一条共享边 \(e\)，如果其两侧的四边形不满足 Delaunay 条件（即某个三角形的外接圆包含对面的顶点），则翻转这条边：删除 \(e\)，连接四边形的另一条对角线。

edge_flip(edge e):
  令 e 连接顶点 a 和 b
  令 c 和 d 分别是 e 两侧的对顶点
  if incircle(a, b, c, d) > 0:   // d 在三角形 abc 的外接圆内
    删除边 ab
    添加边 cd
    递归检查新三角形的邻边

边翻转的一个重要性质：从任意三角剖分开始，通过有限次边翻转一定可以到达 Delaunay 三角剖分。翻转次数的上界是 \(O(n^2)\)，但在增量插入的场景下，每次插入的期望翻转次数是 \(O(1)\)（在随机插入顺序下）。

5.5 复杂度

随机插入顺序：期望 \(O(n \log n)\) 时间。
最坏情况：\(O(n^2)\) 时间（当点按特定顺序排列时，点定位退化）。
使用历史 DAG 进行点定位，可以保证期望 \(O(n \log n)\)。

六、Bowyer-Watson 算法

6.1 算法概述

Bowyer-Watson 算法是另一种增量构造 Delaunay 三角剖分的方法，由 Adrian Bowyer 和 David Watson 在 1981 年独立提出。与上一节的”分裂+翻转”策略不同，Bowyer-Watson 算法采用”删除+重新三角化”的策略。

6.2 核心思想

当插入一个新点 \(p\) 时： 1. 找到所有外接圆包含 \(p\) 的三角形（称为”坏三角形”）。 2. 删除这些坏三角形，形成一个星形多边形（star-shaped polygon）空洞。 3. 将空洞的每条边界边与 \(p\) 相连，形成新的三角形。

6.3 详细步骤

Bowyer-Watson(S):
  1. 创建超级三角形 T0，包含所有点
  2. triangulation = {T0}
  3. for each 点 p in S:
       bad_triangles = {}
       for each 三角形 t in triangulation:
         if p 在 t 的外接圆内:
           bad_triangles.add(t)
       // 找到空洞的边界
       polygon = {}
       for each 三角形 t in bad_triangles:
         for each 边 e of t:
           if e 不被 bad_triangles 中的其他三角形共享:
             polygon.add(e)
       // 删除坏三角形
       for each t in bad_triangles:
         triangulation.remove(t)
       // 用新三角形填充空洞
       for each 边 e in polygon:
         new_tri = triangle(e, p)
         triangulation.add(new_tri)
  4. 删除与超级三角形关联的三角形
  return triangulation

6.4 优化

朴素实现中，对每个新点需要遍历所有三角形来找坏三角形，时间复杂度 \(O(n)\) 每步，总计 \(O(n^2)\)。优化方法：

从包含新点的三角形开始扩展搜索：先定位包含新点的三角形，然后沿着邻接关系扩展。因为坏三角形区域是连通的，所以只需要搜索局部区域。
维护三角形邻接信息：使用半边数据结构（half-edge data structure）可以在 \(O(1)\) 时间内访问相邻三角形。

6.5 Bowyer-Watson 与增量翻转的比较

特性	Bowyer-Watson	增量翻转
核心操作	删除+重建	分裂+翻转
实现难度	较简单	中等
通用性	容易推广到高维	二维特有
数据结构要求	三角形列表+邻接	半边或类似结构
实际性能	稍慢（内存分配多）	稍快（原地修改）

在我的经验中，Bowyer-Watson 更适合原型开发和高维场景，而增量翻转更适合需要极致性能的二维场景。

七、完整 C 实现：增量 Delaunay 三角剖分

下面给出一个完整的、可编译运行的 C 实现。这个实现采用 Bowyer-Watson 算法，因为它的逻辑更清晰，适合教学和理解。

7.1 数据结构

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#define MAX_POINTS    1024
#define MAX_TRIANGLES (MAX_POINTS * 4)
#define EPSILON       1e-10

typedef struct {
    double x, y;
} Point;

typedef struct {
    int v[3];           /* 三个顶点在点数组中的索引 */
    int adj[3];         /* 三个邻接三角形的索引，-1 表示无邻接 */
    double cx, cy, cr;  /* 外接圆圆心和半径 */
    int alive;          /* 是否存活 */
} Triangle;

typedef struct {
    Point    points[MAX_POINTS + 3];  /* +3 用于超级三角形 */
    int      npoints;
    Triangle tris[MAX_TRIANGLES];
    int      ntris;
} Delaunay;

7.2 几何谓词

static double orient2d_val(const Point *a, const Point *b, const Point *c)
{
    return (b->x - a->x) * (c->y - a->y)
         - (b->y - a->y) * (c->x - a->x);
}

/* 计算三角形外接圆 */
static void calc_circumcircle(Delaunay *dt, int ti)
{
    Triangle *t = &dt->tris[ti];
    const Point *a = &dt->points[t->v[0]];
    const Point *b = &dt->points[t->v[1]];
    const Point *c = &dt->points[t->v[2]];

    double ax = a->x, ay = a->y;
    double bx = b->x, by = b->y;
    double cx = c->x, cy = c->y;

    double D = 2.0 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by));
    if (fabs(D) < EPSILON) {
        t->cx = t->cy = 0.0;
        t->cr = -1.0;  /* 退化三角形 */
        return;
    }

    double a2 = ax * ax + ay * ay;
    double b2 = bx * bx + by * by;
    double c2 = cx * cx + cy * cy;

    t->cx = (a2 * (by - cy) + b2 * (cy - ay) + c2 * (ay - by)) / D;
    t->cy = (a2 * (cx - bx) + b2 * (ax - cx) + c2 * (bx - ax)) / D;

    double dx = ax - t->cx;
    double dy = ay - t->cy;
    t->cr = sqrt(dx * dx + dy * dy);
}

/* 判断点 p 是否在三角形 ti 的外接圆内 */
static int in_circumcircle(Delaunay *dt, int ti, const Point *p)
{
    Triangle *t = &dt->tris[ti];
    if (t->cr < 0) return 0;  /* 退化三角形 */
    double dx = p->x - t->cx;
    double dy = p->y - t->cy;
    return (dx * dx + dy * dy) < (t->cr * t->cr - EPSILON);
}

7.3 初始化与超级三角形

void delaunay_init(Delaunay *dt)
{
    dt->npoints = 0;
    dt->ntris = 0;
}

/* 创建包围所有点的超级三角形 */
static void create_super_triangle(Delaunay *dt,
                                  double xmin, double ymin,
                                  double xmax, double ymax)
{
    double dx = xmax - xmin;
    double dy = ymax - ymin;
    double dmax = (dx > dy) ? dx : dy;
    double midx = (xmin + xmax) / 2.0;
    double midy = (ymin + ymax) / 2.0;

    /* 超级三角形的三个顶点 */
    int s0 = dt->npoints;
    int s1 = dt->npoints + 1;
    int s2 = dt->npoints + 2;

    dt->points[s0].x = midx - 20 * dmax;
    dt->points[s0].y = midy - dmax;
    dt->points[s1].x = midx;
    dt->points[s1].y = midy + 20 * dmax;
    dt->points[s2].x = midx + 20 * dmax;
    dt->points[s2].y = midy - dmax;

    int ti = dt->ntris++;
    dt->tris[ti].v[0] = s0;
    dt->tris[ti].v[1] = s1;
    dt->tris[ti].v[2] = s2;
    dt->tris[ti].adj[0] = dt->tris[ti].adj[1] = dt->tris[ti].adj[2] = -1;
    dt->tris[ti].alive = 1;
    calc_circumcircle(dt, ti);
}

7.4 核心插入算法

/* 边结构，用于记录空洞边界 */
typedef struct {
    int v0, v1;    /* 边的两个顶点 */
    int tri_idx;   /* 边所属的三角形 */
    int edge_idx;  /* 边在三角形中的局部索引 */
} Edge;

/* 插入一个新点 */
void delaunay_insert(Delaunay *dt, double x, double y)
{
    int pi = dt->npoints++;
    dt->points[pi].x = x;
    dt->points[pi].y = y;

    /* 步骤 1：找到所有外接圆包含新点的三角形 */
    int bad[MAX_TRIANGLES];
    int nbad = 0;
    Point p = {x, y};

    for (int i = 0; i < dt->ntris; i++) {
        if (!dt->tris[i].alive) continue;
        if (in_circumcircle(dt, i, &p)) {
            bad[nbad++] = i;
        }
    }

    /* 步骤 2：找到空洞的边界多边形 */
    Edge polygon[MAX_TRIANGLES * 3];
    int npoly = 0;

    for (int i = 0; i < nbad; i++) {
        Triangle *t = &dt->tris[bad[i]];
        /* 对三角形的每条边检查 */
        for (int e = 0; e < 3; e++) {
            int v0 = t->v[e];
            int v1 = t->v[(e + 1) % 3];
            /* 检查这条边是否被其他坏三角形共享 */
            int shared = 0;
            for (int j = 0; j < nbad; j++) {
                if (j == i) continue;
                Triangle *t2 = &dt->tris[bad[j]];
                for (int e2 = 0; e2 < 3; e2++) {
                    int u0 = t2->v[e2];
                    int u1 = t2->v[(e2 + 1) % 3];
                    if ((v0 == u1 && v1 == u0) || (v0 == u0 && v1 == u1)) {
                        shared = 1;
                        break;
                    }
                }
                if (shared) break;
            }
            if (!shared) {
                polygon[npoly].v0 = v0;
                polygon[npoly].v1 = v1;
                polygon[npoly].tri_idx = bad[i];
                polygon[npoly].edge_idx = e;
                npoly++;
            }
        }
    }

    /* 步骤 3：删除坏三角形 */
    for (int i = 0; i < nbad; i++) {
        dt->tris[bad[i]].alive = 0;
    }

    /* 步骤 4：从边界多边形的每条边到新点创建新三角形 */
    int new_tris[MAX_TRIANGLES * 3];
    int nnew = 0;

    for (int i = 0; i < npoly; i++) {
        int ti = dt->ntris++;
        Triangle *nt = &dt->tris[ti];
        nt->v[0] = polygon[i].v0;
        nt->v[1] = polygon[i].v1;
        nt->v[2] = pi;
        nt->adj[0] = nt->adj[1] = nt->adj[2] = -1;
        nt->alive = 1;
        calc_circumcircle(dt, ti);
        new_tris[nnew++] = ti;
    }

    /* 步骤 5：设置新三角形之间的邻接关系 */
    for (int i = 0; i < nnew; i++) {
        for (int j = i + 1; j < nnew; j++) {
            Triangle *ti = &dt->tris[new_tris[i]];
            Triangle *tj = &dt->tris[new_tris[j]];
            for (int ei = 0; ei < 3; ei++) {
                for (int ej = 0; ej < 3; ej++) {
                    int vi0 = ti->v[ei], vi1 = ti->v[(ei + 1) % 3];
                    int vj0 = tj->v[ej], vj1 = tj->v[(ej + 1) % 3];
                    if ((vi0 == vj1 && vi1 == vj0) ||
                        (vi0 == vj0 && vi1 == vj1)) {
                        ti->adj[ei] = new_tris[j];
                        tj->adj[ej] = new_tris[i];
                    }
                }
            }
        }
    }
}

7.5 完成与清理

/* 移除与超级三角形关联的三角形 */
void delaunay_finalize(Delaunay *dt, int original_npoints)
{
    int super_start = original_npoints;
    for (int i = 0; i < dt->ntris; i++) {
        if (!dt->tris[i].alive) continue;
        Triangle *t = &dt->tris[i];
        for (int j = 0; j < 3; j++) {
            if (t->v[j] >= super_start) {
                t->alive = 0;
                break;
            }
        }
    }
}

/* 打印所有存活的三角形 */
void delaunay_print(Delaunay *dt)
{
    printf("Delaunay Triangulation:\n");
    int count = 0;
    for (int i = 0; i < dt->ntris; i++) {
        if (!dt->tris[i].alive) continue;
        Triangle *t = &dt->tris[i];
        printf("  Triangle %d: (%.2f, %.2f) - (%.2f, %.2f) - (%.2f, %.2f)\n",
               count++,
               dt->points[t->v[0]].x, dt->points[t->v[0]].y,
               dt->points[t->v[1]].x, dt->points[t->v[1]].y,
               dt->points[t->v[2]].x, dt->points[t->v[2]].y);
        printf("    Circumcenter: (%.4f, %.4f), R=%.4f\n",
               t->cx, t->cy, t->cr);
    }
    printf("Total: %d triangles\n", count);
}

7.6 主函数与测试

int main(void)
{
    Delaunay dt;
    delaunay_init(&dt);

    double test_points[][2] = {
        {0.0, 0.0}, {1.0, 0.0}, {0.5, 1.0},
        {0.0, 1.0}, {1.0, 1.0}, {0.5, 0.5},
        {0.25, 0.75}, {0.75, 0.25}
    };
    int n = sizeof(test_points) / sizeof(test_points[0]);

    double xmin = 1e18, ymin = 1e18, xmax = -1e18, ymax = -1e18;
    for (int i = 0; i < n; i++) {
        if (test_points[i][0] < xmin) xmin = test_points[i][0];
        if (test_points[i][1] < ymin) ymin = test_points[i][1];
        if (test_points[i][0] > xmax) xmax = test_points[i][0];
        if (test_points[i][1] > ymax) ymax = test_points[i][1];
    }

    /* 预加载点，创建超级三角形 */
    for (int i = 0; i < n; i++) {
        dt.points[i].x = test_points[i][0];
        dt.points[i].y = test_points[i][1];
    }
    dt.npoints = n;
    create_super_triangle(&dt, xmin, ymin, xmax, ymax);

    /* 逐点插入——调用前面定义的 delaunay_insert */
    for (int i = 0; i < n; i++)
        delaunay_insert(&dt, test_points[i][0], test_points[i][1]);

    delaunay_finalize(&dt, n);
    delaunay_print(&dt);
    return 0;
}

以上实现优先考虑清晰度而非性能。生产环境中需要考虑：使用精确算术库（如 Jonathan Shewchuk 的 predicates.c）、更高效的点定位、内存池分配等。

八、工程实践：踩坑与对策

8.1 工程陷阱一览

在实际项目中，Voronoi/Delaunay 的工程陷阱比算法本身更让人头疼。下表总结了我在多个项目中遇到的典型问题：

问题	症状	根因	对策
浮点误差导致拓扑错误	三角形重叠、空洞	orient2d 和 incircle 判定不一致	使用 Shewchuk 的精确谓词或自适应精度
四点共圆退化	非法翻转、死循环	未处理共圆情况下的二义性	使用符号扰动（Simulation of Simplicity）
超级三角形太小	边界三角形异常	超级三角形未完全包围点集	使用足够大的缩放因子（10x 以上）
共线点导致退化	外接圆半径无穷大	三点共线时外接圆不存在	预处理检测共线点，使用数值阈值
大规模点集性能差	插入一个点耗时数秒	朴素搜索坏三角形 \(O(n)\)	基于邻接关系的局部搜索 + 随机化点序
内存碎片	频繁分配释放，GC 压力大	三角形频繁创建销毁	使用内存池或空闲链表
Voronoi 无界区域	绘制异常、裁剪错误	边界 Voronoi 区域延伸到无穷远	用包围盒裁剪 Voronoi 边
数值漂移	长时间运行后精度下降	坐标值累积误差	定期重新计算外接圆参数

8.2 浮点精度：最大的敌人

浮点精度问题是 Delaunay 三角剖分工程化的头号难题。考虑以下场景：

/* 不鲁棒的 orient2d 实现 */
double orient2d_naive(double ax, double ay,
                      double bx, double by,
                      double cx, double cy)
{
    return (bx - ax) * (cy - ay) - (by - ay) * (cx - ax);
}

当三点几乎共线时，上述计算可能因浮点舍入返回错误的符号。这会导致”点在三角形内”和”点在三角形外”的判断不一致，进而产生拓扑错误（如三角形重叠、空洞）。

Jonathan Shewchuk 的经典论文《Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates》给出了解决方案。他的方法首先尝试快速的浮点计算，只有在结果接近零时才逐步提高精度，最终使用精确算术保证正确性。

/*
 * Shewchuk 自适应精度 orient2d 的简化思路：
 * 1. 先用普通浮点计算，得到近似值 det
 * 2. 计算误差界 err_bound
 * 3. 如果 |det| > err_bound，直接返回（快速路径）
 * 4. 否则，逐步使用更高精度的展开算术
 */
double orient2d_adaptive(double ax, double ay,
                         double bx, double by,
                         double cx, double cy)
{
    double detleft  = (ax - cx) * (by - cy);
    double detright = (ay - cy) * (bx - cx);
    double det = detleft - detright;

    double detsum;
    if (detleft > 0.0) {
        if (detright <= 0.0) return det;
        detsum = detleft + detright;
    } else if (detleft < 0.0) {
        if (detright >= 0.0) return det;
        detsum = -detleft - detright;
    } else {
        return det;
    }

    double err_bound = 3.3306690738754716e-16 * detsum;
    if ((det >= err_bound) || (-det >= err_bound))
        return det;

    /* 需要更高精度的计算（省略） */
    return orient2d_exact(ax, ay, bx, by, cx, cy);
}

8.3 符号扰动

四点共圆（或三点共线）是另一类退化情况。理论上可以使用符号扰动（Simulation of Simplicity）来处理：对每个点加上一个无穷小的扰动 \(\epsilon_i\)，使得没有四点共圆。实际实现中，这等价于定义一个一致的”打破平局”规则。

/* 符号扰动的简单实现：按点的索引打破平局 */
static int incircle_sos(Delaunay *dt, int ai, int bi, int ci, int di)
{
    double result = incircle_exact(
        dt->points[ai].x, dt->points[ai].y,
        dt->points[bi].x, dt->points[bi].y,
        dt->points[ci].x, dt->points[ci].y,
        dt->points[di].x, dt->points[di].y
    );
    if (fabs(result) > EPSILON) return result > 0;
    /* 共圆情况：按索引最小的点决定 */
    return di < ai || di < bi || di < ci;
}

8.4 性能优化策略

对于大规模点集（百万级），需要以下优化：

（一）随机化插入顺序。 将点随机打乱后插入，可以保证增量算法的期望时间为 \(O(n \log n)\)。这是最简单也最有效的优化。

/* Fisher-Yates 洗牌 */
void shuffle_points(int *indices, int n)
{
    for (int i = n - 1; i > 0; i--) {
        int j = rand() % (i + 1);
        int tmp = indices[i];
        indices[i] = indices[j];
        indices[j] = tmp;
    }
}

（二）BRIO 排序。 Biased Randomized Insertion Order（BRIO）是 Amenta、Choi 和 Rote 提出的一种改进排序方式。它将点按空间位置分层（类似希尔伯特曲线排序），每层内随机打乱。这样做的好处是利用了空间局部性，使点定位更快。

（三）并行化。 Delaunay 三角剖分可以通过分而治之并行化：将点集空间分割为若干区域，分别构造局部三角剖分，然后合并。合并步骤是难点，但在现代多核处理器上可以获得显著加速。

九、应用场景

9.1 最近设施查询

Voronoi 图最直接的应用就是最近设施查询。给定一组设施（如医院、加油站），Voronoi 图将空间划分为若干区域，查询点所在的区域对应的设施就是最近的设施。

class NearestFacilityQuery:
    """基于 Voronoi 图的最近设施查询"""

    def __init__(self, facilities):
        self.facilities = facilities
        # 预处理：构建 Voronoi 图 + 点定位结构
        self.voronoi = build_voronoi(facilities)
        self.locator = build_point_location(self.voronoi)

    def query(self, point):
        """查询离 point 最近的设施"""
        region = self.locator.locate(point)
        return self.facilities[region.site_index]

结合梯形分解或 Kirkpatrick 层次结构，单次查询可以做到 \(O(\log n)\)。这比暴力搜索的 \(O(n)\) 快得多，尤其是当查询量很大时。

9.2 网格生成（有限元法）

Delaunay 三角剖分在有限元分析中的应用可能是工业界最重要的应用。有限元法需要将计算域离散为三角形网格，网格质量直接影响数值解的精度和稳定性。

Delaunay 三角剖分的”最大化最小角”性质保证了网格中没有过于”瘦长”的三角形，这对有限元的刚度矩阵条件数至关重要。

Delaunay 细化算法（如 Ruppert 算法）可以在保持 Delaunay 性质的同时，通过插入新点来改善网格质量：

Ruppert-Refinement(DT, min_angle):
  while 存在角度小于 min_angle 的三角形 t:
    计算 t 的外接圆圆心 c
    if c 侵犯了某条约束边 e:
      在 e 的中点插入新点
    else:
      在 c 处插入新点
    更新 Delaunay 三角剖分

Ruppert 算法可以保证生成的三角形最小角不小于约 20.7 度（理论界），实际应用中可以放松到约 33 度。

9.3 自然邻域插值

自然邻域插值（Natural Neighbor Interpolation）是一种基于 Voronoi 图的空间插值方法。给定一组带值的采样点和一个查询点 \(q\)，插值值为：

\[f(q) = \sum_{i} w_i \cdot f(p_i)\]

其中权重 \(w_i\) 由 Voronoi 区域的”被偷面积”决定。将 \(q\) 插入 Voronoi 图后，\(q\) 的 Voronoi 区域会”偷取”相邻区域的一部分面积，偷取面积的比例就是对应站点的权重。

def natural_neighbor_interpolation(sites, values, query):
    """自然邻域插值"""
    # 构建原始 Voronoi 图
    original_voronoi = build_voronoi(sites)
    original_areas = [region.area for region in original_voronoi.regions]

    # 插入查询点后重建 Voronoi 图
    new_sites = sites + [query]
    new_voronoi = build_voronoi(new_sites)

    # 计算被偷面积作为权重
    query_region = new_voronoi.regions[-1]  # 查询点的区域
    weights = []
    for i, site in enumerate(sites):
        stolen_area = original_areas[i] - new_voronoi.regions[i].area
        if stolen_area > 0:
            weights.append((i, stolen_area))

    total_weight = sum(w for _, w in weights)
    result = sum(values[i] * w / total_weight for i, w in weights)
    return result

自然邻域插值的优点是： - 插值是连续的（C1 连续），过渡自然 - 不需要指定搜索半径或邻居数 - 在采样点处精确还原（忠实插值）

9.4 地理信息系统

在 GIS 中，Voronoi 图和 Delaunay 三角剖分有广泛应用：

Thiessen 多边形：气象站观测数据的空间分析，用 Voronoi 图划分每个气象站的代表区域。
TIN（不规则三角网）：基于 Delaunay 三角剖分的地形建模，从离散高程点构建地表模型。
可达性分析：城市规划中，用 Voronoi 图分析公共设施（学校、医院）的服务范围。
空间聚类：结合 Voronoi 图和密度估计进行空间数据聚类。

9.5 路径规划

Voronoi 图在机器人路径规划中也有重要应用。给定一组障碍物，计算障碍物边界点的 Voronoi 图，其边就是”离所有障碍物最远的路径”。沿着 Voronoi 边移动可以最大化安全距离。

def voronoi_path_planning(obstacles, start, goal):
    """基于 Voronoi 图的路径规划"""
    # 1. 提取障碍物边界点
    boundary_points = extract_boundary_points(obstacles)

    # 2. 构建 Voronoi 图
    voronoi = build_voronoi(boundary_points)

    # 3. 在 Voronoi 图上构建道路图
    road_map = build_road_map(voronoi)

    # 4. 将起点和终点连接到最近的 Voronoi 边
    road_map.connect(start)
    road_map.connect(goal)

    # 5. 在道路图上搜索最短路径
    path = dijkstra(road_map, start, goal)
    return path

十、高维推广

10.1 高维 Voronoi 图

Voronoi 图的定义可以自然推广到任意维度：

\[V(p_i) = \{x \in \mathbb{R}^d \mid d(x, p_i) \leq d(x, p_j), \forall j \neq i\}\]

然而，高维 Voronoi 图面临严峻的维度灾难。对于 \(d\) 维空间中的 \(n\) 个点，Voronoi 图的复杂度可以达到 \(\Theta(n^{\lceil d/2 \rceil})\)。例如：

维度	Voronoi 顶点上界	n=100 时的量级
2	\(O(n)\)	~200
3	\(O(n^2)\)	~10,000
4	\(O(n^2)\)	~10,000
5	\(O(n^3)\)	~1,000,000
10	\(O(n^5)\)	~\(10^{10}\)

这意味着在高维空间中，精确计算 Voronoi 图/Delaunay 三角剖分是不切实际的。

10.2 高维替代方案

面对高维情况，实际可行的替代方案包括：

（一）近似最近邻。 使用 LSH（局部敏感哈希）、随机投影树、HNSW 等近似最近邻算法，放弃精确的 Voronoi 结构。

（二）降维后构造。 使用 PCA、t-SNE 等方法将高维数据投影到低维空间，然后构造 Voronoi 图。但这会损失信息。

（三）利用抛物面提升。 有一个优美的数学事实：\(d\) 维 Delaunay 三角剖分等价于 \(d+1\) 维下凸包。将每个点 \((x_1, \ldots, x_d)\) 提升到抛物面 \((x_1, \ldots, x_d, x_1^2 + \cdots + x_d^2)\)，然后计算下凸包，其投影回原空间就是 Delaunay 三角剖分。

def delaunay_via_lifting(points_2d):
    """通过抛物面提升将 2D Delaunay 转化为 3D 凸包"""
    # 将 2D 点提升到 3D 抛物面
    lifted = []
    for x, y in points_2d:
        lifted.append((x, y, x*x + y*y))

    # 计算 3D 下凸包
    lower_hull = convex_hull_lower(lifted)

    # 下凸包的面投影回 2D 就是 Delaunay 三角形
    triangles = []
    for face in lower_hull.faces:
        tri = [(lifted[v][0], lifted[v][1]) for v in face.vertices]
        triangles.append(tri)

    return triangles

这个提升对应关系是 Voronoi/Delaunay 理论中最深刻的结果之一，它揭示了这两种结构与凸包之间的内在联系。

10.3 加权 Voronoi 图

在某些应用中，不同站点的”影响力”不同。加权 Voronoi 图（也称 Power Diagram）使用加权距离：

\[V_w(p_i) = \{x \mid d(x, p_i)^2 - w_i \leq d(x, p_j)^2 - w_j, \forall j \neq i\}\]

其中 \(w_i\) 是站点 \(p_i\) 的权重。加权 Voronoi 图的对偶是加权 Delaunay 三角剖分（Regular Triangulation）。

十一、算法变体与扩展

11.1 约束 Delaunay 三角剖分

在实际应用中，经常需要三角剖分尊重某些预定义的边（如地形边界、建筑轮廓）。约束 Delaunay 三角剖分（Constrained Delaunay Triangulation，CDT）在满足约束边的前提下尽可能保持 Delaunay 性质。

CDT 的构造方法：

先构造不带约束的 Delaunay 三角剖分。
逐条插入约束边。如果约束边与现有边交叉，删除交叉的边，重新三角化。
对非约束边执行边翻转，恢复 Delaunay 性质（约束边不参与翻转）。

def constrained_delaunay(points, constraints):
    """构造约束 Delaunay 三角剖分"""
    # 第一步：标准 Delaunay 三角剖分
    dt = delaunay_triangulate(points)

    # 第二步：插入约束边
    for edge in constraints:
        p, q = edge
        if not dt.has_edge(p, q):
            # 找到与 pq 相交的所有边
            crossing_edges = dt.find_crossing_edges(p, q)
            # 删除交叉边，形成两个多边形
            upper_polygon, lower_polygon = dt.remove_edges(crossing_edges)
            # 重新三角化两个多边形
            dt.triangulate_polygon(upper_polygon)
            dt.triangulate_polygon(lower_polygon)
            # 标记 pq 为约束边
            dt.mark_constrained(p, q)

    # 第三步：边翻转恢复 Delaunay 性质
    for edge in dt.non_constrained_edges():
        if not dt.is_locally_delaunay(edge):
            dt.flip(edge)

    return dt

11.2 动态 Delaunay 三角剖分

在某些应用中（如粒子模拟），点集会随时间变化。动态 Delaunay 三角剖分需要支持：

插入：增量算法自然支持。
删除：删除一个点后，需要重新三角化产生的空洞。可以使用”耳朵剪除”（ear removal）方法。
移动：等价于删除+插入，但可以优化——如果移动距离小，只需要局部更新。

11.3 分治算法

Guibas 和 Stolfi 在 1985 年提出了一种优雅的分治 Delaunay 三角剖分算法：

Divide-Conquer-Delaunay(S):
  if |S| <= 3:
    直接构造三角剖分
    return
  将 S 按 x 坐标分为左半 S_L 和右半 S_R
  DT_L = Divide-Conquer-Delaunay(S_L)
  DT_R = Divide-Conquer-Delaunay(S_R)
  return Merge(DT_L, DT_R)

合并步骤需要找到左右两部分的”公切线”并逐步缝合。这个算法保证最坏情况 \(O(n \log n)\)，不需要随机化。但实现复杂度较高，且缓存性能不如增量算法。

11.4 不同距离度量下的 Voronoi 图

除了欧几里得距离，还可以使用其他度量：

曼哈顿距离 Voronoi 图：区域边界由水平、垂直和 45 度线组成。用于 VLSI 设计中的线路规划。
切比雪夫距离 Voronoi 图：区域边界是轴对齐的多边形。
测地线距离 Voronoi 图：在曲面或带障碍物的环境中，使用最短路径距离。用于机器人导航。

十二、个人思考与总结

12.1 我的实践体会

在多年的工程实践中，我对 Voronoi/Delaunay 有以下几点深刻体会：

第一，永远不要自己从零写精确谓词。 Shewchuk 的 predicates.c 是经过数十年考验的工业级实现。自己写的浮点谓词迟早会在某个极端情况下出问题。我曾在一个项目中自信满满地用朴素 double 运算实现了 Delaunay 三角剖分，在测试中一切正常，上线后遇到了一组近乎共线的点，整个三角剖分崩溃了。

第二，理解对偶性比记住公式更重要。 很多人把 Voronoi 和 Delaunay 当作两个独立的概念来学习，实际上它们是同一个数学对象的两个视角。理解了对偶性，你就能在两者之间自由切换，选择最适合当前问题的表示。

第三，增量算法加随机化是工程中的首选。 虽然 Fortune 算法和分治算法在理论上更优雅，但增量 Bowyer-Watson 算法实现简单、调试容易、性能足够。加上随机化插入顺序后，期望时间复杂度就是最优的 \(O(n \log n)\)。

第四，不要忽视边界情况。 在我参与的每一个使用 Delaunay 三角剖分的项目中，80% 的 bug 都来自边界情况：共线点、共圆点、重合点、极端坐标值。花在处理边界情况上的时间往往超过核心算法本身。

第五，选择合适的库。 除非有特殊需求（如需要定制化数据结构），否则不要自己实现。CGAL（C++）、Triangle（C）、scipy.spatial.Delaunay（Python）都是成熟的选择。

12.2 算法选择指南

场景	推荐算法	理由
原型开发、小规模数据	Bowyer-Watson	实现简单，代码短
大规模静态数据	分治算法或 BRIO+增量	最坏情况保证，缓存友好
需要 Voronoi 图	Fortune 扫描线或 DT 转换	Fortune 直接输出 Voronoi
动态点集	增量+删除	支持在线更新
约束三角剖分	CDT 专用算法	需要尊重约束边
高维数据	放弃精确构造，使用 ANN	维度灾难不可避免

12.3 算法复杂度总结

算法	时间复杂度	空间复杂度	是否最优
Fortune 扫描线	\(O(n \log n)\)	\(O(n)\)	是
分治 Delaunay	\(O(n \log n)\)	\(O(n)\)	是
增量 Delaunay（随机化）	\(O(n \log n)\) 期望	\(O(n)\)	期望最优
Bowyer-Watson（随机化）	\(O(n \log n)\) 期望	\(O(n)\)	期望最优
增量 Delaunay（最坏）	\(O(n^2)\)	\(O(n)\)	否
半平面交（朴素）	\(O(n^2 \log n)\)	\(O(n)\)	否

12.4 进一步阅读

Voronoi 图和 Delaunay 三角剖分的文献浩如烟海。以下是我认为最值得阅读的资料：

de Berg, M., Cheong, O., van Kreveld, M., Overmars, M. 《Computational Geometry: Algorithms and Applications》（第三版）。第 7 章和第 9 章。计算几何的标准教材，讲解清晰。
Shewchuk, J. R. 《Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates》, Discrete & Computational Geometry, 1997。精确谓词的开创性工作，附带可直接使用的 C 代码。
Fortune, S. 《A Sweepline Algorithm for Voronoi Diagrams》, Algorithmica, 1987。Fortune 扫描线算法的原始论文。
Guibas, L., Stolfi, J. 《Primitives for the Manipulation of General Subdivisions and the Computation of Voronoi Diagrams》, ACM Transactions on Graphics, 1985。分治算法和四边结构的经典论文。
Bowyer, A. 《Computing Dirichlet Tessellations》, The Computer Journal, 1981。Bowyer-Watson 算法的原始论文之一。
Ruppert, J. 《A Delaunay Refinement Algorithm for Quality 2-Dimensional Mesh Generation》, Journal of Algorithms, 1995。Delaunay 细化算法的重要论文。
Shewchuk, J. R. 《Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator》。Triangle 库的技术报告，包含大量工程实践经验。
Aurenhammer, F. 《Voronoi Diagrams – A Survey of a Fundamental Geometric Data Structure》, ACM Computing Surveys, 1991。Voronoi 图的经典综述。

算法系列导航：上一篇：扫描线算法 | 下一篇：R-tree 与空间索引

相关阅读：凸包算法 | 最近点对

目录