7.2 二维Delaunay三角化 (Delaunay_triangulation_2)
上一节: 7.1 二维三角化基础 下一节: 7.3 约束三角化 相关章节: 11.1 Voronoi图 (Delaunay对偶关系)
0. 直观理解：为什么需要Delaunay三角化？
生活类比：气泡堆积
想象你在吹一堆气泡，每个气泡的中心是一个数据点。当气泡相互挤压时，它们之间的接触线就形成了Delaunay三角化的边。这就是为什么Delaunay三角化被称为”对偶”于Voronoi图的原因——就像气泡的外壁和内壁的关系。
气泡堆积示意图：
O O O O
\ / \ /
\ / \ /
\ / \ /
---O--- vs ---O---
/ \ /|\
/ \ / | \
/ \ / | \
O O O | O
/ \
/ \
O O
不规则三角化 Delaunay三角化
(细长三角形多) (三角形更规则)
为什么Delaunay三角化很重要？
问题场景：假设你正在构建一个地形模型，需要将地形采样点连接成三角网格。
- 不好的三角化：产生大量细长三角形，导致地形表面出现不合理的陡坡
- Delaunay三角化：最大化最小角，避免细长三角形，地形更平滑
核心优势：
- 最优三角形质量：在所有可能的三角化中，Delaunay三角化最大化最小内角
- 最近邻关系：自动建立点集的自然邻接关系
- 对偶Voronoi图：同时获得Voronoi划分，用于空间分析
1. 理论基础
1.1 Delaunay三角化的数学定义
Delaunay三角化是计算几何中最优的三角化方式，由俄罗斯数学家Boris Delaunay于1934年提出。
形式化定义：给定平面点集 ，Delaunay三角化 是满足以下条件的三角化：
空外接圆性质（Empty Circumcircle Property）：对于 中的每个三角形，其外接圆内部不包含任何其他输入点。
数学表达：对于 中的任意三角形 ，满足：
存在性与唯一性定理：
- 存在性：对于任意有限点集 ，Delaunay三角化存在
- 唯一性：如果点集处于一般位置（general position），即不存在四点共圆，则Delaunay三角化唯一
对于非一般位置的点集（存在共圆点），Delaunay三角化不唯一，但Delaunay图（Delaunay graph）唯一。
1.2 空外接圆性质详解
空外接圆性质是Delaunay三角化的核心判定标准。让我们通过可视化来理解：
Delaunay三角形的外接圆：
p_i
/ \\
/ \\
/ O \\ O = 外接圆心（到三个顶点等距）
/ | \\ 圆内无其他点
/ | \\
p_j---+----p_k
半径R
非Delaunay情况：
p_i
/ \\
/ \\ p_l ← 点在外接圆内！
/ O /\\
/ | / \\
/ |/ \\
p_j---+----p_k
违反空外接圆性质，需要边翻转：将边(p_j, p_k)翻转为边(p_i, p_l)
几何解释：空外接圆性质确保三角形尽可能”均匀”，避免细长三角形。直观上，如果某点落在三角形的外接圆内，将该点与三角形顶点连接通常会产生更规则的三角形。
1.3 最大化最小角性质
Delaunay三角化具有最优的三角形质量：
定理（Delaunay, 1934）：在所有可能的三角化中，Delaunay三角化最大化最小内角。
角度比较示例：
非Delaunay三角化： Delaunay三角化：
p_3 p_3
/\\ /|\\
/ \\ / | \\
/ \\ / | \\
/ 5° \\ / | \\
p_1------p_2 p_1---+----p_2
\\ / \\ | /
\\ / \\ | /
\\ / \\|/
\\ / p_4
p_4
最小角 = 5° 最小角 = 30°
（细长三角形） （三角形更规则）
应用意义：
- 有限元分析：避免细长三角形导致的数值不稳定
- 地形建模：防止不合理的陡坡和地形扭曲
- 计算机图形学：更均匀的网格带来更好的渲染效果
- 插值计算：规则三角形的插值误差更小
1.4 对偶Voronoi图关系
Delaunay三角化与Voronoi图互为对偶图（Dual Graph）：
Voronoi图定义：对于点集 ，点 的Voronoi单元为：
对偶关系：
| Delaunay三角化 | Voronoi图 |
|---|---|
| 顶点（输入点） | Voronoi单元 |
| Delaunay边 | Voronoi边（单元边界） |
| 三角形外接圆心 | Voronoi顶点 |
| 三角形 | Voronoi顶点关联的区域 |
对偶关系可视化：
Delaunay三角化（实线） + Voronoi图（虚线）：
p_1
/|\
/ | \\ p_1, p_2, p_3, p_4 是Delaunay顶点
/ | \\ v_1, v_2 是Voronoi顶点（外接圆心）
/ v_1--...
/ /|\\ 实线 = Delaunay边
p_4-/-|-\\-p_2 虚线 = Voronoi边
\\/ v_2 \\\
\\ | / v_1 是三角形(p_1,p_2,p_4)的外接圆心
\\ | / v_2 是三角形(p_2,p_3,p_4)的外接圆心
\\|/
p_3
对偶性质：
- 每条Delaunay边垂直于对应的Voronoi边
- Delaunay边连接的两个点是相邻Voronoi单元的生成点
应用意义：
- 最近邻查询：Voronoi图天然支持最近邻搜索
- 设施选址：Voronoi单元表示最近服务区域
- 自然邻域插值：基于Voronoi单元的权重插值
- 形态分析：骨架提取、形状描述
1.5 退化情况处理
当点集不满足一般位置假设时，需要特殊处理：
共圆点（Cocircular Points）：多个点位于同一圆上
共圆点情况：
p_1 -------- p_2
/ \\
/ O 圆心 \\
/ \\
p_4 ------------ p_3
所有四点共圆，Delaunay三角化不唯一：
方案1: 三角形(p_1,p_2,p_3) + 三角形(p_1,p_3,p_4)
方案2: 三角形(p_1,p_2,p_4) + 三角形(p_2,p_3,p_4)
两种方案都满足空外接圆性质！
CGAL使用**符号扰动（Symbolic Perturbation）**技术处理退化：
- 通过模拟扰动打破平局
- 确保结果的一致性和唯一性
- 保持拓扑正确性
共线点（Collinear Points）：所有点位于同一直线上
共线点情况：
p_1 -- p_2 -- p_3 -- p_4 -- p_5
退化为1维：没有有限三角形
Delaunay三角化仅包含连接相邻点的边
CGAL使用"无限远顶点"技术统一处理
2. 算法详解
2.1 增量插入算法（Incremental Insertion）
CGAL采用增量算法构建Delaunay三角化，逐个插入点并维护Delaunay性质。
算法框架：
// 伪代码：增量Delaunay三角化
Delaunay_triangulation_2 build_delaunay(Point_set P) {
// 1. 初始化：创建包含所有点的超大三角形
Triangle super_triangle = create_super_triangle(P);
Delaunay_triangulation_2 dt(super_triangle);
// 2. 随机打乱点序（避免最坏情况）
random_shuffle(P.begin(), P.end());
// 3. 逐个插入点
for (Point p : P) {
insert_point(dt, p);
}
// 4. 移除超三角形相关顶点
remove_super_triangle(dt, super_triangle);
return dt;
}单点插入过程：
void insert_point(Delaunay_triangulation_2& dt, Point p) {
// 步骤1: 定位 - 找到包含点p的三角形
Face_handle f = dt.locate(p);
// 步骤2: 分裂 - 将三角形分裂为3个
if (p在f内部) {
split_triangle(dt, f, p); // 1个三角形 -> 3个三角形
} else if (p在f的边上) {
split_edge(dt, f, p); // 2个三角形 -> 4个三角形
} else {
// p与现有点重合，跳过或更新
return;
}
// 步骤3: 恢复Delaunay性质 - 边翻转
restore_delaunay(dt, p);
}2.2 Bowyer-Watson算法步骤
Bowyer-Watson算法是增量插入的经典实现，使用**冲突区域（Conflict Region）**概念：
Bowyer-Watson算法可视化：
步骤1: 插入新点p，找到包含它的三角形
a
/\\
/ \\
/ \\ 新点p在三角形(a,b,c)内
/ p \\
/ O \\
/__________\\
b c
步骤2: 找到所有外接圆包含p的三角形（冲突区域）
a
/\\
/ \\
/ T1 \\ T1, T2, T3的外接圆都包含p
/ p \\
/ O------O \\
/ T2 T3 \\
b____________c
冲突区域 = {T1, T2, T3} 形成"星形多边形"
步骤3: 删除冲突区域内的所有三角形
a
/ \\
/ \\
/ p \\
/ O \\
/ / \\ \\
/ / \\ \\
b/___________\\c
步骤4: 将p与边界连接，创建新三角形
a
/|\\
/ | \\
/ | \\
/ T1| T2\\
/____p____ \\
/ T3 | T4 \\
/______|______\\
b c
新三角形: (p,a,b), (p,b,c), (p,c,a) 等
算法复杂度分析：
- 定位步骤： 使用历史DAG或 使用行走法
- 冲突检测：，d为冲突区域大小（平均常数）
- 边翻转：
- 总体平均： 每点，总体
2.3 边翻转（Edge Flip）机制
边翻转是维护Delaunay性质的核心操作。
可翻转条件：两个相邻三角形形成凸四边形
边翻转操作：
翻转前（非Delaunay）： 翻转后（Delaunay）：
a a
/\\ /|\\
/ \\ / | \\
/ T1 \\ / | \\
/ \\ / T1| T2 \\
b--------c --> b---p---c
\\ / \\ | /
\\ T2 / \\ | /
\\ / \\|/
\\/ p
p
边(b,c)违反Delaunay性质， 翻转为边(a,p)
因为p在三角形(a,b,c)的外接圆内 现在满足空外接圆性质
边翻转判定：
// 检查边是否需要翻转
bool needs_flip(Face_handle f, int i) {
// f是面，i是边索引
Face_handle n = f->neighbor(i); // 邻接面
if (is_infinite(n)) return false; // 边界边不翻转
Vertex_handle v_op = n->vertex(n->index(f)); // 对角顶点
// 核心测试：对角顶点是否在外接圆内
// side_of_oriented_circle返回:
// ON_POSITIVE_SIDE: 在圆外（满足Delaunay，不翻转）
// ON_NEGATIVE_SIDE: 在圆内（违反Delaunay，需要翻转）
// ON_ORIENTED_BOUNDARY: 在圆上（退化情况）
return side_of_oriented_circle(f, v_op->point()) == ON_NEGATIVE_SIDE;
}
// 执行边翻转
void flip(Face_handle f, int i) {
Face_handle n = f->neighbor(i);
int j = n->index(f); // f在n中的索引
// 重新连接四个顶点
// 翻转前: f = (a,b,c), n = (b,p,c) 共享边(b,c)
// 翻转后: f = (a,b,p), n = (a,p,c) 共享边(a,p)
// 更新顶点关联关系
// 更新邻接面指针
// 维护组合数据结构的一致性
}边翻转算法流程：
void restore_delaunay(Vertex_handle v) {
// 收集所有与新顶点v关联的边
std::stack<Edge> edges_to_check;
Face_circulator fc = incident_faces(v), done(fc);
do {
int i = fc->index(v);
edges_to_check.push(Edge(fc, i));
} while (++fc \!= done);
// 边翻转循环
while (\!edges_to_check.empty()) {
Edge e = edges_to_check.top();
edges_to_check.pop();
Face_handle f = e.first;
int i = e.second;
if (needs_flip(f, i)) {
Face_handle n = f->neighbor(i);
// 执行翻转
flip(f, i);
// 将新产生的边加入检查队列
// 翻转产生两条新边，需要检查它们是否满足Delaunay
int ni = n->index(f);
edges_to_check.push(Edge(f, (i+1)%3));
edges_to_check.push(Edge(n, (ni+1)%3));
}
}
}2.4 时间复杂度分析
| 操作 | 平均复杂度 | 最坏复杂度 | 说明 |
|---|---|---|---|
| 构造（n个点） | 随机增量算法 | ||
| 插入单个点 | 期望复杂度，使用历史DAG | ||
| 删除顶点 | 需要重建局部区域 | ||
| 点定位 | walking策略 | ||
| 边翻转 | 局部操作 | ||
| 最近邻查询 | 利用Delaunay性质 | ||
| Voronoi顶点计算 | 外接圆心计算 |
复杂度对比：
不同三角化算法复杂度比较：
算法 构造复杂度 空间复杂度
--------------------------------------------------------
暴力算法 O(n^3) O(n)
扫描线算法 O(n log n) O(n)
分治算法 O(n log n) O(n)
增量插入(无随机) O(n^2) O(n)
增量插入(随机) O(n log n) O(n) <- CGAL默认
Delaunay层次结构 O(n log n) O(n) <- 最快，常数小
3. 实现细节
3.1 CGAL的Delaunay_triangulation_2类
类定义：
namespace CGAL {
template <class Traits_, class Tds_ = Default>
class Delaunay_triangulation_2
: public Triangulation_2<Traits_, Tds_> {
public:
// 类型定义
typedef Traits_ Geom_traits;
typedef Tds_ Triangulation_data_structure;
typedef typename Traits_::Point_2 Point;
typedef typename Traits_::Segment_2 Segment;
typedef typename Traits_::Triangle_2 Triangle;
// 句柄和迭代器（继承自Triangulation_2）
typedef typename Tds_::Vertex_handle Vertex_handle;
typedef typename Tds_::Face_handle Face_handle;
typedef typename Tds_::Edge Edge;
typedef typename Tds_::Vertex_iterator Vertex_iterator;
typedef typename Tds_::Face_iterator Face_iterator;
typedef typename Tds_::Edge_iterator Edge_iterator;
typedef typename Tds_::Face_circulator Face_circulator;
typedef typename Tds_::Vertex_circulator Vertex_circulator;
// 构造与析构
Delaunay_triangulation_2(const Geom_traits& traits = Geom_traits());
Delaunay_triangulation_2(InputIterator first, InputIterator last);
// 插入操作
Vertex_handle insert(const Point& p, Face_handle hint = Face_handle());
Vertex_handle insert(const Point& p, Locate_type lt, Face_handle loc, int li);
int insert(InputIterator first, InputIterator last);
// 删除操作
void remove(Vertex_handle v);
void remove(InputIterator first, InputIterator last);
// 点定位
Face_handle locate(const Point& p, Face_handle hint = Face_handle()) const;
Face_handle locate(const Point& p, Locate_type& lt, int& li,
Face_handle hint = Face_handle()) const;
// Delaunay特定查询
Vertex_handle nearest_vertex(const Point& p, Face_handle hint = Face_handle()) const;
bool is_delaunay(Face_handle f) const;
// Voronoi图相关
Point dual(Face_handle f) const; // Voronoi顶点（外接圆心）
Object dual(Edge e) const; // Voronoi边
Object dual(Edge_circulator ec) const;
Object dual(Edge_iterator ei) const;
// 几何谓词
Oriented_side side_of_oriented_circle(Face_handle f, const Point& p) const;
// 验证
bool is_valid(bool verbose = false, int level = 0) const;
};
} // namespace CGAL继承层次：
Triangulation_2<Traits, Tds> // 基类：通用三角化
└── Delaunay_triangulation_2 // Delaunay特化
└── Regular_triangulation_2 // 加权Delaunay推广
3.2 顶点/面/边的遍历
基本遍历模式：
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <iostream>
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using DT = CGAL::Delaunay_triangulation_2<K>;
using Point = K::Point_2;
using Vertex_handle = DT::Vertex_handle;
using Face_handle = DT::Face_handle;
using Edge = DT::Edge;
void traversal_examples(DT& dt) {
// ========== 顶点遍历 ==========
std::cout << "=== 所有顶点 ===" << std::endl;
for (auto vit = dt.vertices_begin(); vit \!= dt.vertices_end(); ++vit) {
Point p = vit->point();
std::cout << "顶点: (" << p.x() << ", " << p.y() << ")" << std::endl;
}
// 仅有限顶点（不包括无限远顶点）
for (auto vit = dt.finite_vertices_begin();
vit \!= dt.finite_vertices_end(); ++vit) {
// 处理有限顶点
}
// ========== 面遍历 ==========
std::cout << "\n=== 所有面 ===" << std::endl;
for (auto fit = dt.faces_begin(); fit \!= dt.faces_end(); ++fit) {
std::cout << "三角形: ";
for (int i = 0; i < 3; ++i) {
Point p = fit->vertex(i)->point();
std::cout << "(" << p.x() << ", " << p.y() << ") ";
}
std::cout << std::endl;
}
// 仅有限面（不包括与无限远顶点关联的面）
for (auto fit = dt.finite_faces_begin();
fit \!= dt.finite_faces_end(); ++fit) {
// 处理有限三角形
}
// ========== 边遍历 ==========
std::cout << "\n=== 所有边 ===" << std::endl;
for (auto eit = dt.edges_begin(); eit \!= dt.edges_end(); ++eit) {
Edge e = *eit; // Edge = pair<Face_handle, int>
Segment seg = dt.segment(e);
std::cout << "边: (" << seg.source().x() << ", " << seg.source().y()
<< ") -> (" << seg.target().x() << ", " << seg.target().y()
<< ")" << std::endl;
}
// 仅有限边
for (auto eit = dt.finite_edges_begin();
eit \!= dt.finite_edges_end(); ++eit) {
// 处理有限边
}
}循环器（Circulators）遍历邻域：
void circulator_examples(DT& dt, Vertex_handle v) {
// ========== 绕顶点的面循环器 ==========
std::cout << "=== 顶点 " << v->point() << " 的关联面 ===" << std::endl;
DT::Face_circulator fc = dt.incident_faces(v), done(fc);
if (fc \!= nullptr) { // 检查顶点是否关联面
do {
if (\!dt.is_infinite(fc)) { // 跳过无限面
std::cout << "三角形: ";
for (int i = 0; i < 3; ++i) {
Point p = fc->vertex(i)->point();
std::cout << "(" << p.x() << ", " << p.y() << ") ";
}
std::cout << std::endl;
}
++fc;
} while (fc \!= done);
}
// ========== 绕顶点的顶点循环器 ==========
std::cout << "\n=== 顶点 " << v->point() << " 的邻接顶点 ===" << std::endl;
DT::Vertex_circulator vc = dt.incident_vertices(v), vdone(vc);
if (vc \!= nullptr) {
do {
if (\!dt.is_infinite(vc)) { // 跳过无限远顶点
Point p = vc->point();
std::cout << "邻接顶点: (" << p.x() << ", " << p.y() << ")" << std::endl;
}
++vc;
} while (vc \!= vdone);
}
// ========== 绕边的顶点循环器 ==========
Edge e = ...; // 某条边
DT::Vertex_circulator ec = dt.incident_vertices(e);
// 遍历与边e关联的所有顶点（即共享该边的所有三角形）
}3.3 邻接关系查询
面-顶点关系：
void face_vertex_relationships(DT& dt) {
// 获取面的三个顶点
Face_handle f = ...;
Vertex_handle v0 = f->vertex(0);
Vertex_handle v1 = f->vertex(1);
Vertex_handle v2 = f->vertex(2);
// 获取顶点在面中的索引
int idx = f->index(v0); // 返回0
// 判断顶点是否属于面
bool has_vertex = f->has_vertex(v0); // true
// 获取面的对边顶点
Vertex_handle opposite = f->vertex(f->cw(0)); // 顺时针下一个
Vertex_handle opposite2 = f->vertex(f->ccw(0)); // 逆时针下一个
}面-面邻接关系：
void face_adjacency(DT& dt) {
Face_handle f = ...;
// 获取三个邻接面
Face_handle n0 = f->neighbor(0); // 对边0的邻接面
Face_handle n1 = f->neighbor(1); // 对边1的邻接面
Face_handle n2 = f->neighbor(2); // 对边2的邻接面
// 获取邻接面在本面中的索引
int idx_in_neighbor = f->index(n0); // n0在f中的索引
// 获取本面在邻接面中的索引
int idx_in_f = n0->index(f); // f在n0中的索引
// 遍历所有邻接面
std::cout << "面 " << f << " 的邻接面:" << std::endl;
for (int i = 0; i < 3; ++i) {
Face_handle ni = f->neighbor(i);
if (\!dt.is_infinite(ni)) {
std::cout << " 边 " << i << " 的邻接面: " << ni << std::endl;
} else {
std::cout << " 边 " << i << " 是边界边（邻接无限面）" << std::endl;
}
}
}高级邻接查询：
#include <vector>
#include <set>
// 获取顶点的所有邻接顶点
std::vector<Vertex_handle> get_adjacent_vertices(DT& dt, Vertex_handle v) {
std::vector<Vertex_handle> neighbors;
DT::Vertex_circulator vc = dt.incident_vertices(v), done(vc);
if (vc \!= nullptr) {
do {
if (\!dt.is_infinite(vc)) {
neighbors.push_back(vc);
}
++vc;
} while (vc \!= done);
}
return neighbors;
}
// 获取面的所有邻接面
std::vector<Face_handle> get_adjacent_faces(DT& dt, Face_handle f) {
std::vector<Face_handle> neighbors;
for (int i = 0; i < 3; ++i) {
Face_handle ni = f->neighbor(i);
if (\!dt.is_infinite(ni)) {
neighbors.push_back(ni);
}
}
return neighbors;
}
// 查找两个顶点的公共邻接顶点
std::vector<Vertex_handle> get_common_neighbors(DT& dt, Vertex_handle v1, Vertex_handle v2) {
std::set<Vertex_handle> n1(get_adjacent_vertices(dt, v1).begin(),
get_adjacent_vertices(dt, v1).end());
std::set<Vertex_handle> n2(get_adjacent_vertices(dt, v2).begin(),
get_adjacent_vertices(dt, v2).end());
std::vector<Vertex_handle> common;
std::set_intersection(n1.begin(), n1.end(), n2.begin(), n2.end(),
std::back_inserter(common));
return common;
}3.4 几何谓词（Predicate）
几何谓词是计算几何算法的基础，CGAL使用精确算术保证鲁棒性。
核心谓词：
// 1. 方向测试（Orientation Test）
// 判断三个点的方向：顺时针、逆时针或共线
Orientation orientation(const Point& p, const Point& q, const Point& r);
// 返回值: LEFT_TURN (逆时针), RIGHT_TURN (顺时针), COLLINEAR (共线)
// 2. 有向圆测试（Side of Oriented Circle）
// 判断点d相对于(a,b,c)的有向外接圆的位置
Oriented_side side_of_oriented_circle(const Point& a, const Point& b,
const Point& c, const Point& d);
// 返回值: ON_POSITIVE_SIDE (圆外), ON_NEGATIVE_SIDE (圆内),
// ON_ORIENTED_BOUNDARY (圆上)
// 3. 共线点上的顺序
// 判断共线点r是否在线段(p,q)上
bool collinear_are_ordered_along_line(const Point& p, const Point& r,
const Point& q);
// 4. 距离比较
Comparison_result compare_distance(const Point& p, const Point& q,
const Point& r);
// 比较 |p-q| 和 |p-r| 的大小谓词使用示例：
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <iostream>
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point = K::Point_2;
void predicate_examples() {
Point a(0, 0), b(1, 0), c(0, 1), d(0.5, 0.5);
// 方向测试
K::Orientation_2 orientation;
CGAL::Orientation orient = orientation(a, b, c);
std::cout << "orientation(a,b,c) = "
<< (orient == CGAL::LEFT_TURN ? "LEFT_TURN" :
orient == CGAL::RIGHT_TURN ? "RIGHT_TURN" : "COLLINEAR")
<< std::endl;
// 圆测试
K::Side_of_oriented_circle_2 circle_test;
CGAL::Oriented_side side = circle_test(a, b, c, d);
std::cout << "side_of_oriented_circle(a,b,c,d) = "
<< (side == CGAL::ON_POSITIVE_SIDE ? "OUTSIDE" :
side == CGAL::ON_NEGATIVE_SIDE ? "INSIDE" : "ON_CIRCLE")
<< std::endl;
// 距离比较
K::Compare_distance_2 compare_dist;
CGAL::Comparison_result cmp = compare_dist(a, b, c);
std::cout << "compare_distance(a,b,c) = "
<< (cmp == CGAL::SMALLER ? "|ab| < |ac|" :
cmp == CGAL::LARGER ? "|ab| > |ac|" : "|ab| = |ac|")
<< std::endl;
}构造器（Constructions）：
// 1. 外接圆心计算
K::Construct_circumcenter_2 circumcenter;
Point cc = circumcenter(a, b, c); // 三角形(a,b,c)的外接圆心
// 2. 中点计算
K::Construct_midpoint_2 midpoint;
Point mid = midpoint(a, b);
// 3. 线段构造
K::Construct_segment_2 segment;
Segment s = segment(a, b);
// 4. 三角形构造
K::Construct_triangle_2 triangle;
Triangle t = triangle(a, b, c);4. 使用示例
4.1 基本使用示例
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <iostream>
#include <vector>
#include <random>
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using DT = CGAL::Delaunay_triangulation_2<K>;
using Point = K::Point_2;
using Vertex_handle = DT::Vertex_handle;
using Face_handle = DT::Face_handle;
int main() {
// 创建Delaunay三角化对象
DT dt;
std::cout << "=== 基本Delaunay三角化示例 ===" << std::endl;
// ========== 插入点 ==========
std::cout << "\n--- 插入点 ---" << std::endl;
// 方法1: 逐个插入
Vertex_handle v1 = dt.insert(Point(0, 0));
Vertex_handle v2 = dt.insert(Point(1, 0));
Vertex_handle v3 = dt.insert(Point(0.5, 1));
std::cout << "逐个插入3个点" << std::endl;
std::cout << "顶点数: " << dt.number_of_vertices() << std::endl;
std::cout << "面数: " << dt.number_of_faces() << std::endl;
// 方法2: 批量插入（更高效）
std::vector<Point> points;
std::mt19937 gen(42);
std::uniform_real_distribution<> dis(0.0, 1.0);
for (int i = 0; i < 100; ++i) {
points.push_back(Point(dis(gen), dis(gen)));
}
dt.insert(points.begin(), points.end());
std::cout << "\n批量插入100个随机点" << std::endl;
std::cout << "总顶点数: " << dt.number_of_vertices() << std::endl;
std::cout << "总面数: " << dt.number_of_faces() << std::endl;
std::cout << "边数: " << dt.number_of_edges() << std::endl;
// ========== 基本查询 ==========
std::cout << "\n--- 基本查询 ---" << std::endl;
// 检查是否为空
std::cout << "三角化是否为空: " << (dt.is_empty() ? "是" : "否") << std::endl;
// 获取维度
std::cout << "维度: " << dt.dimension() << std::endl;
// 0: 只有顶点, 1: 共线, 2: 平面三角化
// 检查是否为Delaunay（验证）
std::cout << "是否满足Delaunay性质: "
<< (dt.is_valid() ? "是" : "否") << std::endl;
// ========== 遍历输出 ==========
std::cout << "\n--- 遍历有限面 ---" << std::endl;
int face_count = 0;
for (auto fit = dt.finite_faces_begin();
fit \!= dt.finite_faces_end(); ++fit) {
if (face_count < 5) { // 只输出前5个
std::cout << "面 " << ++face_count << ": ";
for (int i = 0; i < 3; ++i) {
Point p = fit->vertex(i)->point();
std::cout << "(" << p.x() << ", " << p.y() << ") ";
}
std::cout << std::endl;
}
}
// ========== 几何输出 ==========
std::cout << "\n--- 输出三角形 ---" << std::endl;
face_count = 0;
for (auto fit = dt.finite_faces_begin();
fit \!= dt.finite_faces_end(); ++fit) {
if (face_count < 3) {
Triangle tri = dt.triangle(fit);
std::cout << "三角形 " << ++face_count << ": ";
std::cout << "面积=" << tri.area() << std::endl;
}
}
return 0;
}4.2 顶点插入和删除
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <iostream>
#include <vector>
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using DT = CGAL::Delaunay_triangulation_2<K>;
using Point = K::Point_2;
using Vertex_handle = DT::Vertex_handle;
void insertion_deletion_example() {
DT dt;
std::cout << "=== 顶点插入和删除 ===" << std::endl;
// ========== 批量插入优化 ==========
std::cout << "\n--- 批量插入优化 ---" << std::endl;
std::vector<Point> points;
for (int i = 0; i < 20; ++i) {
for (int j = 0; j < 20; ++j) {
points.push_back(Point(i, j));
}
}
// 方法1: 直接插入（CGAL会自动进行空间排序优化）
dt.insert(points.begin(), points.end());
std::cout << "直接插入 " << points.size() << " 个点" << std::endl;
std::cout << "顶点数: " << dt.number_of_vertices() << std::endl;
// 方法2: 手动空间排序（当知道点分布时更高效）
DT dt2;
CGAL::spatial_sort(points.begin(), points.end());
dt2.insert(points.begin(), points.end());
std::cout << "空间排序后插入，结果相同但可能更快" << std::endl;
// ========== 带提示的插入 ==========
std::cout << "\n--- 带提示的插入 ---" << std::endl;
// 当知道新点大概位置时，提供提示面可加速定位
Point new_point(10.5, 10.5);
Face_handle hint = dt.locate(Point(10, 10)); // 先定位到附近
Vertex_handle v = dt.insert(new_point, hint);
std::cout << "插入点 (10.5, 10.5)，使用提示加速" << std::endl;
// ========== 顶点删除 ==========
std::cout << "\n--- 顶点删除 ---" << std::endl;
std::cout << "删除前顶点数: " << dt.number_of_vertices() << std::endl;
// 删除单个顶点
dt.remove(v);
std::cout << "删除后顶点数: " << dt.number_of_vertices() << std::endl;
// 批量删除
std::vector<Vertex_handle> to_remove;
int count = 0;
for (auto vit = dt.finite_vertices_begin();
vit \!= dt.finite_vertices_end() && count < 10; ++vit, ++count) {
to_remove.push_back(vit);
}
dt.remove(to_remove.begin(), to_remove.end());
std::cout << "批量删除10个顶点后: " << dt.number_of_vertices() << std::endl;
// ========== 移动顶点 ==========
std::cout << "\n--- 移动顶点 ---" << std::endl;
Vertex_handle vm = dt.insert(Point(5, 5));
std::cout << "原位置: (" << vm->point().x() << ", " << vm->point().y() << ")" << std::endl;
// 移动顶点到新位置（如果新位置很近，这比删除+插入更高效）
dt.move(vm, Point(5.5, 5.5));
std::cout << "新位置: (" << vm->point().x() << ", " << vm->point().y() << ")" << std::endl;
// 验证Delaunay性质仍然保持
std::cout << "移动后仍满足Delaunay: " << (dt.is_valid() ? "是" : "否") << std::endl;
}
int main() {
insertion_deletion_example();
return 0;
}4.3 最近邻查询
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <iostream>
#include <vector>
#include <cmath>
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using DT = CGAL::Delaunay_triangulation_2<K>;
using Point = K::Point_2;
using Vertex_handle = DT::Vertex_handle;
void nearest_neighbor_example() {
DT dt;
std::cout << "=== 最近邻查询 ===" << std::endl;
// 创建测试点集
std::vector<Point> points = {
Point(0, 0), Point(1, 0), Point(0, 1), Point(1, 1),
Point(0.5, 0.5), Point(2, 2), Point(-1, -1),
Point(0.5, 0), Point(0, 0.5), Point(0.25, 0.25)
};
dt.insert(points.begin(), points.end());
// ========== 单最近邻查询 ==========
std::cout << "\n--- 单最近邻查询 ---" << std::endl;
Point query(0.6, 0.6);
Vertex_handle nearest = dt.nearest_vertex(query);
std::cout << "查询点: (" << query.x() << ", " << query.y() << ")" << std::endl;
std::cout << "最近邻: (" << nearest->point().x()
<< ", " << nearest->point().y() << ")" << std::endl;
// 计算实际距离
double dist = std::sqrt(CGAL::squared_distance(query, nearest->point()));
std::cout << "欧氏距离: " << dist << std::endl;
// ========== k近邻查询 ==========
std::cout << "\n--- k近邻查询（手动实现）---" << std::endl;
// CGAL 2D Delaunay不直接支持k近邻，需要手动实现
auto k_nearest = [&dt](const Point& q, int k) {
std::vector<std::pair<Vertex_handle, double>> neighbors;
// 遍历所有顶点计算距离
for (auto vit = dt.finite_vertices_begin();
vit \!= dt.finite_vertices_end(); ++vit) {
double d = std::sqrt(CGAL::squared_distance(q, vit->point()));
neighbors.push_back({vit, d});
}
// 按距离排序
std::partial_sort(neighbors.begin(),
neighbors.begin() + std::min(k, (int)neighbors.size()),
neighbors.end(),
[](const auto& a, const auto& b) {
return a.second < b.second;
});
neighbors.resize(std::min(k, (int)neighbors.size()));
return neighbors;
};
auto knn = k_nearest(query, 3);
std::cout << "3个最近邻:" << std::endl;
for (const auto& [vh, dist] : knn) {
std::cout << " (" << vh->point().x() << ", " << vh->point().y()
<< ") 距离=" << dist << std::endl;
}
// ========== 范围查询 ==========
std::cout << "\n--- 范围查询 ---" << std::endl;
double radius = 0.5;
std::cout << "查询点 (" << query.x() << ", " << query.y() << ") 周围 "
<< radius << " 范围内的点:" << std::endl;
for (auto vit = dt.finite_vertices_begin();
vit \!= dt.finite_vertices_end(); ++vit) {
double d = std::sqrt(CGAL::squared_distance(query, vit->point()));
if (d <= radius) {
std::cout << " (" << vit->point().x() << ", " << vit->point().y()
<< ") 距离=" << d << std::endl;
}
}
// ========== 利用Delaunay性质加速范围查询 ==========
std::cout << "\n--- 利用Delaunay邻域加速 ---" << std::endl;
// 先找到最近邻，然后在其邻域内搜索
Vertex_handle seed = dt.nearest_vertex(query);
std::vector<Vertex_handle> local_candidates;
// 获取seed的邻接顶点
DT::Vertex_circulator vc = dt.incident_vertices(seed), done(vc);
if (vc \!= nullptr) {
do {
if (\!dt.is_infinite(vc)) {
local_candidates.push_back(vc);
}
++vc;
} while (vc \!= done);
}
std::cout << "最近邻的局部邻域有 " << local_candidates.size() << " 个顶点" << std::endl;
std::cout << "这些顶点是候选的近邻点" << std::endl;
}
int main() {
nearest_neighbor_example();
return 0;
}4.4 自然邻域插值
自然邻域插值（Natural Neighbor Interpolation）是一种基于Voronoi图的高精度插值方法。
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <iostream>
#include <vector>
#include <utility>
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using DT = CGAL::Delaunay_triangulation_2<K>;
using Point = K::Point_2;
using Vertex_handle = DT::Vertex_handle;
// 带高度信息的顶点
typedef std::pair<Vertex_handle, double> Vertex_value;
void natural_neighbor_interpolation_example() {
DT dt;
std::cout << "=== 自然邻域插值 ===" << std::endl;
// 创建采样点（带高度值）
std::vector<std::pair<Point, double>> samples = {
{{0, 0}, 0.0},
{{1, 0}, 1.0},
{{0, 1}, 1.0},
{{1, 1}, 2.0},
{{0.5, 0.5}, 1.0},
{{0.5, 0}, 0.5},
{{0, 0.5}, 0.5},
{{1, 0.5}, 1.5},
{{0.5, 1}, 1.5}
};
// 存储顶点与值的映射
std::map<Vertex_handle, double> values;
for (const auto& [p, v] : samples) {
Vertex_handle vh = dt.insert(p);
values[vh] = v;
}
std::cout << "插入 " << samples.size() << " 个采样点" << std::endl;
// ========== 自然邻域坐标 ==========
std::cout << "\n--- 计算自然邻域坐标 ---" << std::endl;
Point query(0.6, 0.6);
std::cout << "查询点: (" << query.x() << ", " << query.y() << ")" << std::endl;
// 计算自然邻域坐标
std::vector<std::pair<Vertex_handle, double>> coords;
double normalization_factor;
CGAL::natural_neighbor_coordinates_2(dt, query,
std::back_inserter(coords),
normalization_factor);
std::cout << "自然邻域坐标 (未归一化):" << std::endl;
for (const auto& [vh, coord] : coords) {
std::cout << " 顶点 (" << vh->point().x() << ", " << vh->point().y()
<< "): 坐标=" << coord << std::endl;
}
std::cout << "归一化因子: " << normalization_factor << std::endl;
// ========== 插值计算 ==========
std::cout << "\n--- 插值计算 ---" << std::endl;
// 使用自然邻域坐标进行加权平均
double interpolated_value = 0.0;
for (const auto& [vh, coord] : coords) {
interpolated_value += coord * values[vh];
}
interpolated_value /= normalization_factor;
std::cout << "插值结果: " << interpolated_value << std::endl;
// 与最近邻插值比较
Vertex_handle nearest = dt.nearest_vertex(query);
double nearest_value = values[nearest];
std::cout << "最近邻插值结果: " << nearest_value << std::endl;
// ========== 梯度估计 ==========
std::cout << "\n--- 梯度估计 ---" << std::endl;
// 使用自然邻域坐标估计梯度
// 梯度 = sum(coord_i * value_i * (p_i - query)) / normalization
double grad_x = 0.0, grad_y = 0.0;
for (const auto& [vh, coord] : coords) {
double dx = vh->point().x() - query.x();
double dy = vh->point().y() - query.y();
grad_x += coord * values[vh] * dx;
grad_y += coord * values[vh] * dy;
}
grad_x /= normalization_factor;
grad_y /= normalization_factor;
std::cout << "估计梯度: (" << grad_x << ", " << grad_y << ")" << std::endl;
// ========== 多个查询点插值 ==========
std::cout << "\n--- 网格插值 ---" << std::endl;
std::cout << "在规则网格上插值:" << std::endl;
for (double x = 0.2; x <= 0.8; x += 0.2) {
for (double y = 0.2; y <= 0.8; y += 0.2) {
Point p(x, y);
std::vector<std::pair<Vertex_handle, double>> nncoords;
double norm;
CGAL::natural_neighbor_coordinates_2(dt, p,
std::back_inserter(nncoords), norm);
double val = 0.0;
for (const auto& [vh, coord] : nncoords) {
val += coord * values[vh];
}
val /= norm;
std::cout << " (" << x << ", " << y << "): " << val << std::endl;
}
}
}
int main() {
natural_neighbor_interpolation_example();
return 0;
}5. 复杂度分析
5.1 时间/空间复杂度表格
| 操作 | 平均复杂度 | 最坏复杂度 | 空间复杂度 | 说明 |
|---|---|---|---|---|
| 构造（n个点） | 随机增量算法 | |||
| 插入单个点 | 期望复杂度 | |||
| 删除顶点 | 需要重建局部区域 | |||
| 移动顶点 | 小位移时高效 | |||
| 点定位 | walking策略 | |||
| 最近邻查询 | 利用Delaunay性质 | |||
| 遍历所有面 | 线性遍历 | |||
| Voronoi图计算 | 计算所有对偶 | |||
| 验证 | 检查Delaunay性质 |
5.2 与其他三角化对比
特性对比表：
特性 Delaunay 约束Delaunay 正则三角化
-------------------------------------------------------------------
构造复杂度 O(n log n) O(n log n) O(n log n)
插入复杂度 O(log n) O(log n) O(log n)
空外接圆性质 有 有（局部） 有（正交圆）
最大化最小角 是 是（约束内） 否
支持约束边 否 是 否
支持加权点 否 否 是
Voronoi图 标准 受限Voronoi 幂图
最近邻查询 O(log n) O(log n) O(log n)
空间复杂度 O(n) O(n) O(n + h)
h=隐藏点数
适用场景选择：
- Delaunay_triangulation_2：通用高质量三角化，需要最优三角形质量
- Constrained_Delaunay_triangulation_2：需要保留特征边（如地形中的山脊线）
- Regular_triangulation_2：处理加权点，如不同大小的球体填充
6. 可视化内容
6.1 空外接圆性质的ASCII图示
Delaunay三角化的空外接圆验证：
示例1: 满足Delaunay性质的三角化
p_1(0,2)
/\\
/ \\
/ \\
/ O \\ O = 外接圆心
/ /|\\ \\ R = 外接圆半径
/ R/ | \\R \\ 圆内无其他点
/ / | \\ \\
/__/____|____\\__\\
p_0(0,0) | p_2(2,0)
|
p_3(2,2) ← 在圆外（满足Delaunay）
验证：p_3到外接圆心距离 > R
示例2: 违反Delaunay性质的情况
p_1(0,2)
/\\
/ \\
/ \\
/ O \\
/ /|\\ \\
/ R/ | \\R \\
/ / | \\ \\
/__/____|____\\__\\
p_0(0,0) | p_2(2,0)
|
p_3(1,1) ← 在圆内（违反Delaunay）
验证：p_3到外接圆心距离 < R
需要边翻转：将边(p_0,p_2)翻转为边(p_1,p_3)
翻转后：
p_1(0,2)
/|\\
/ | \\
/ | \\
/ | \\
/ | \\
/ | \\
/ | \\
/_______|_______\\
p_0(0,0) p_3(1,1) p_2(2,0)
现在两个三角形都满足空外接圆性质！
6.2 边翻转过程的可视化
边翻转详细过程：
步骤1: 识别需要翻转的边
p_1 p_1
/\\ /|\\
/ \\ / | \\
/ T1 \\ / | \\
/ \\ --> / T1| T2 \\
/ O \\ / | \\
p_0---p---p_2 p_0--p_3--p_2
\\ / \\ | /
\\ T2 / \\ | /
\\ / \\ | /
\\/ \\|/
p_3 p_3
T1 = 三角形(p_0,p_1,p) T1' = 三角形(p_0,p_1,p_3)
T2 = 三角形(p_0,p,p_2) T2' = 三角形(p_1,p_2,p_3)
边(p_0,p)需要翻转为边(p_1,p_3)
步骤2: 检查翻转条件
翻转前四边形(p_0,p_1,p_2,p_3)必须是凸的：
凸四边形（可翻转）： 凹四边形（不可翻转）：
p_1 p_1
/ \\ / \\
/ \\ / \\
/ \\ / \\
p_0 p_2 p_0------p_2
\\ / \\ /
\\ / \\ /
\\ / \\/
p_3 p_3
凸的 凹的
(对角线在内部) (对角线在外部)
步骤3: 执行翻转
翻转前邻接关系：
- T1: 顶点(p_0,p_1,p), 邻接(T3,T2,T4)
- T2: 顶点(p_0,p,p_2), 邻接(T5,T1,T6)
翻转后邻接关系：
- T1': 顶点(p_0,p_1,p_3), 邻接(T3,T2',T6)
- T2': 顶点(p_1,p_2,p_3), 邻接(T1',T5,T4)
步骤4: 验证Delaunay性质
翻转后检查新边(p_1,p_3)是否满足Delaunay：
- 三角形(p_0,p_1,p_3)的外接圆不包含p_2
- 三角形(p_1,p_2,p_3)的外接圆不包含p_0
如果满足，翻转成功！
如果不满足，继续翻转其他边...
6.3 Delaunay vs 非Delaunay对比图
相同点集的不同三角化对比：
输入点集：
p_1 p_2
\\ /
\\ /
\\ /
p_3
/ \\
/ \\
/ \\
p_4 p_5
非Delaunay三角化（贪心算法）：
p_1-------p_2
\\ /
\\ /
\\ /
\\ /
p_3
/ \\
/ \\
/ \\
p_4-------p_5
问题：
- 三角形(p_1,p_2,p_3)非常细长
- 最小角约15度
- 三角形(p_3,p_4,p_5)同样细长
- 数值稳定性差
- 插值误差大
Delaunay三角化：
p_1-------p_2
\\ /
\\ /
\\ /
\\ /
p_3
/ \\
/ \\
/ \\
p_4-------p_5
优化：
- 边(p_1,p_2)被翻转为边(p_3,p_x)
- 三角形更规则
- 最小角约45度
- 数值稳定性好
- 插值精度高
角度统计对比：
+-------------------+------------+------------+
| 指标 | 非Delaunay | Delaunay |
+-------------------+------------+------------+
| 最小角 | 15.2° | 42.8° |
| 最大角 | 127.5° | 87.3° |
| 平均角 | 58.4° | 60.0° |
| 标准差 | 32.1° | 15.6° |
| 细长三角形数量 | 3 | 0 |
+-------------------+------------+------------+
结论：Delaunay三角化产生更规则的三角形，提高后续计算质量！
7. 总结
Delaunay_triangulation_2是CGAL中最重要、最常用的三角化类：
核心优势
- 最优三角形质量：最大化最小角，避免细长三角形
- 高效算法：随机增量算法，期望 构造时间
- 丰富的查询：点定位、最近邻、Voronoi图等
- 数值鲁棒性：精确谓词保证几何正确性
适用场景
- 地形建模：从采样点构建高质量地形网格
- 有限元分析：生成数值稳定的计算网格
- 计算机图形学：网格变形、纹理映射、 morphing
- 插值与逼近：自然邻域插值、曲面重建
- 最近邻搜索：高效的空间索引结构
使用建议
- 内核选择：默认使用
Exact_predicates_inexact_constructions_kernel - 批量插入：使用
insert(first, last)而非逐个插入 - 空间排序：对大规模数据先进行
spatial_sort - 验证调试：开发时使用
is_valid()检查三角化正确性
学习路径
- 掌握
Triangulation_2的基本概念和操作 - 理解Delaunay性质及其几何意义
- 学习边翻转算法和增量构造
- 探索Voronoi图和对偶关系
- 应用到实际问题：插值、网格生成、空间查询
Delaunay三角化是计算几何的基石，理解其原理和实现对于解决各类几何问题至关重要。CGAL的 Delaunay_triangulation_2 提供了高效、鲁棒、易用的实现，是学习和应用的理想工具。