企业官网建设流程全解析

简介：一份基于C++实现的Delaunay三角网算法源码包，面向需要学习计算几何、网格生成或进行GIS/图形学开发的开发者与在校学生。资源围绕Delaunay三角剖分的核心原理展开，包含完整工程文件、头文件与实现代码，可帮助理解空圆检测、增量插入及局部优化等关键环节。压缩包共42个文件，以h/cpp源码为主，辅以Visual Studio工程配置、MFC界面资源、图标位图及编译生成的obj/exe等，整体体积2.66MB，目录结构清晰，覆盖开发调试各阶段，便于直接编译运行或对照学习。已有1881人学习下载，适合作为算法入门与实践的参考素材。通过阅读代码，读者可掌握从点集数据结构、三角网构建到动态插入与优化的完整实现思路，也能借助附带的可执行文件快速验证算法效果，并以此为基础扩展更多应用场景，这些内容对深化计算几何理解与工程实践均有帮助。 做地形建模或者点云网格化的时候，Delaunay三角网算法几乎是绕不开的一个基础工具。这东西用C++从头实现一遍，比翻十篇论文管用得多，它能把最小角最大化、空外接圆这些几何性质落到实处，让你对计算几何的理解完全上一个台阶。这篇文章我会把从零实现Delaunay三角网的过程完整拆开，从原理到代码再到调试经验，一条线捋下来，适合正在学计算几何的学生，也适合工作中要处理点云、网格、高程数据的程序员。无论你是想搞懂算法本身，还是单纯想抄一份能跑的C++实现，这篇都能给你省下不少时间。

1. Delaunay三角网核心概念：从几何直觉到工程价值

1.1 空外接圆与最大化最小角：两个等价定义

Delaunay三角网的定义有很多种说法，但最常用的是空外接圆性质：在三角网中，任意一个三角形的外接圆内部不包含其他顶点。注意是内部不包含，边界上可以共圆。这个性质听起来简单，但它直接决定了三角网的形状质量。

另一个等价定义是最小角最大化。在所有可能的三角剖分中，Delaunay三角网能使所有三角形的最小内角达到最大。也就是说，它能尽量避免出现那种特别扁长、尖细的三角形。为什么在乎这个？在地形建模里，一个细长三角形往往意味着极差的插值稳定性——稍微挪动一个顶点，高程结果就可能剧烈抖动；在有限元网格生成里，细长三角形会导致刚度矩阵条件数恶化，数值求解直接变慢甚至发散。

这两个定义表面上是两回事，实际上是一致的。证明思路不复杂：如果某个三角形的外接圆内包含了另一个顶点，那么交换这个四边形的一条对角线，生成的两个新三角形的最小角一定比原来大。反复进行这种边翻转，最终就能收敛到Delaunay三角剖分。我在项目里用这个性质做过一个简单的优化器：先随便剖分，再逐边检查翻转条件，三四轮迭代之后网格质量肉眼可见地提升。虽然效率不如主流算法，但对理解本质非常有帮助。

1.2 为什么值得用C++重写一遍

说实话，如果只是为了用，直接调CGAL或者OpenCV的Subdiv2D就行，没必要自己造轮子。但自己用C++实现一遍的价值在于三件事。

第一是性能把控。处理几十万个点云时，Python脚本光是在循环里判断外接圆就要卡半天，C++配合连续内存的vector和局部性良好的遍历，做同样的事能快两三个数量级。我这里说的不是“C++比Python快”这种空泛结论，而是Delaunay算法的核心瓶颈就是大量的几何判定和内存操作，C++在这两方面几乎没有解释型语言的额外开销。

第二是数据结构理解。Delaunay三角网看起来只是一堆三角形，但背后涉及邻接关系、边界提取、动态删除与重建，这些恰恰是面试高频题“如何设计高效的几何拓扑结构”的活教材。你把它实现一遍，再去看半边结构（Half-Edge）就非常轻松。

第三是可控性。CGAL功能强但学习曲线陡，OpenCV的Subdiv2D只能做矩形区域的剖分且不好扩展。自己实现的话，想加个约束边、想输出带邻接信息的网格，都能随心所欲。我后来做地形TIN简化时，就是在自己这套实现上改的，如果当初只调库，改动成本会高很多。

2. 算法选型与方案设计：为什么选Bowyer-Watson

2.1 三种主流算法对比

Delaunay三角网的主流构建算法有三类：分治法、扫描线法和逐点插入法。我一开始也纠结过用哪个，后来根据实际需求锁定了逐点插入法中的Bowyer-Watson算法。

分治法思路是把点集按x坐标分成左右两半，分别递归构建，然后合并两个子三角网。合并过程需要找到上下两条公切线，这个步骤写起来非常绕，我试过一次，调试了整整两天才跑通，最后发现还有边界的特殊情况没处理。扫描线法类似，只是改为从左到右扫描，处理“事件点”时的状态维护同样不省心。

Bowyer-Watson的思路就直白得多：每插入一个新点，找出所有外接圆包含该点的三角形（称为坏三角形），把它们删掉，形成一个空洞，然后把空洞的边界边与新点连接，生成新的三角形。整个循环可以在20行伪代码内描述清楚。对于绝大多数项目来说，点集在几万到几十万量级时，期望O(n log n)的性能完全够用。而且因为它是增量式的，你天然具备了动态插入点的能力，这是其他两类算法不具备的灵活性。

2.2 数据结构设计：三角形列表、邻接关系与超级三角形

实现之前我先想清楚一个关键问题：三角形的存储结构怎么定。

最朴素的做法是直接存三角形顶点坐标，每个三角形是一个独立单元。但这样存在问题：删除坏三角形时，我需要知道哪些边被两个坏三角形共享，哪些边是边界边。如果三角形之间没有邻接信息，我就得双重遍历，事情就麻烦了。

我用的方案是“顶点索引数组 + 三角形索引数组”。所有点存放在一个std::vector 里，三角形只是一个包含三个int索引的结构体。这样每个点只存一次，三角形通过索引引用顶点，内存紧凑且在插入新点时可以复用索引。

struct Point2D { double x, y; }; struct Triangle { int a, b, c; };

为了快速统计边被几个坏三角形共用，我临时用一个std::map<std::pair<int, int>, int>来计数，把无序的边对统一成有序的（比如a<b），这样同一条边无论从哪个三角形遍历到，都能对应到同一个键。

另一个容易忽略的点是“超级三角形”。逐点插入法在初始状态需要一个大三角形把所有点都包住，这个三角形必须在算法结束后删除。超级三角形如果选得太小，会漏掉一些点导致剖分不完整；选得太大，虽然不会出错，但会让外接圆判定的坐标数值变大，影响浮点精度。我的经验是：先求出点集的包围盒，然后让超级三角形的三个顶点距离包围盒边界至少为包围盒对角线长度的20倍。这样既保证包裹完整，又不至于让坐标值大到超过double的安全范围。

3. C++核心实现拆解

3.1 基础结构定义与工具函数

在动手写主逻辑前，我会先准备几个工具函数，它们决定了代码的整洁度和后续调试成本。

首先是点去重。实际工程中拿到的点云经常有大量重复点，如果不处理，后续几何判定会遇到除以零的情况。我用一个基于哈希的unordered_set配合格点化来去重：将坐标除以一个容差，取整数作为哈希键。容差一般取点云包围盒对角线长度的1e-9，这样既能识别真正的重复点，又不会误删间距极小的有效点。

其次是归一化。直接把原始坐标代入外接圆判定公式，如果坐标量级在十万级别，平方项会达到10^10，double做减法时精度会受影响。我先把所有点减去质心坐标，让数值分布在原点附近，再乘一个缩放系数让包围盒边长约为1。这样几何计算时参与运算的数的量级都在1附近，精度问题小很多。

3.2 外接圆判定：数学推导与关键代码

外接圆判定是整个算法的核心热路径，每插入一个点，都要对当前所有三角形跑一次。所以这里的实现质量直接决定整体性能。

我的实现思路是解线性方程组。设三角形三个顶点为A、B、C，外接圆圆心为O，那么O到A、B、C的距离相等，由此可以列出关于O坐标的二元一次方程组。把A作为基准点，令A为原点，则B、C、P（待测点）都变为相对A的向量，推导如下：

O·B = |B|^2 / 2 O·C = |C|^2 / 2

这是一个2x2线性方程组，用克莱姆法则求解O，然后比较OP的距离与OA的距离即可。

inline bool inCircumcircle(const std::vector<Point2D>& pts, int a, int b, int c, const Point2D& p) { double ax = pts[a].x, ay = pts[a].y; double bx = pts[b].x - ax, by = pts[b].y - ay; double cx = pts[c].x - ax, cy = pts[c].y - ay; double dx = p.x - ax, dy = p.y - ay; double denom = 2.0 * (bx * cy - by * cx); if (std::fabs(denom) < 1e-12) { // 三点共线，视为退化的三角形，直接忽略 return false; } double ox = (cy * (bx * bx + by * by) - by * (cx * cx + cy * cy)) / denom; double oy = (bx * (cx * cx + cy * cy) - cx * (bx * bx + by * by)) / denom; double r2 = ox * ox + oy * oy; double dist = (dx - ox) * (dx - ox) + (dy - oy) * (dy - oy); return dist <= r2; }

这里的denom是三角形面积的2倍，如果它接近零，说明三点共线，根本构不成外接圆。实际点集经过去重后很少出现这种情况，但保底判断必须有，否则会出现除零崩溃。数值上还有个细节：如果denom为负数，说明三角形顶点顺序是顺时针，这不影响外接圆判定结果，因为最后比较的是距离平方。

3.3 逐点插入主循环：定位、删除、重建

主循环的逻辑非常清晰，但细节在“边的收集”这一步。每个坏三角形有三条边，如果一条边被两个坏三角形同时拥有，那它是内部边，删除三角形后这条边应该被”消化”掉，不参与新三角形的构建；如果一条边只被一个坏三角形拥有，那它是空洞的边界边，需要与当前插入点连接生成新三角形。

我第一版实现时直接用双重循环去判断“这条边是否同时属于两个坏三角形”，倒也能跑，但点多了以后整个构建过程肉眼可见地变慢。原因很简单，每次插入点，判断边是否重复是O(m^2)的，m为坏三角形数量，在点集密集区域m可能达到几十上百，这个开销不可忽视。

更聪明的做法是使用我前面提到的map计数：

std::map<std::pair<int, int>, int> edgeCount; for (const auto& tri : badTriangles) { addEdge(edgeCount, tri.a, tri.b); addEdge(edgeCount, tri.b, tri.c); addEdge(edgeCount, tri.c, tri.a); }

addEdge里统一把较小索引放在pair的前面，这样同一条边从不同方向遍历时能对上。遍历完所有坏三角形后，edgeCount中计数为1的边就是边界边。

然后删除坏三角形，遍历所有计数为1的边，用边的两个端点加上当前插入点构造新三角形，加入三角形列表。这里有个小细节：新三角形的顶点索引要保证顺序一致，我统一按逆时针方向存储，这样后续计算面积、法向量时不用再判断方向。

3.4 结果清洗与边界处理

所有点都插入完毕后，超级三角形的三个顶点还参与在一些三角形中。算法最后一步是遍历三角形列表，删除所有包含超级三角形顶点索引的三角形。这个步骤容易踩的坑是：删除后列表会有空洞，但你先不要急着压缩，因为接下来可能还有边界三角形的处理需求。

实际项目中，超级三角形删除后，三角网的凸包边界可能不够规整。你拿到的点集可能本身有凹性，Delaunay三角网默认生成的是凸包剖分，凹进去的区域会被“三角化填满”。如果应用场景需要的是凹包，必须在算法后处理阶段执行边界裁剪。常用的做法是在生成完整三角网后，把网格数据和原始点集做一次点在多边形内判断，剔除那些重心在原始多边形外的三角形。这一步在纯凸点集场景可以跳过，但在地形建模时几乎必然要做，因为真实地形边界很少是凸的。

4. 实操中的典型问题与性能优化

4.1 退化输入与数值精度坑

我最开始测试算法时，用的是一组规则网格点。结果发现剖分结果出现很多重叠三角形，程序甚至崩溃过。排查了很久，定位到问题是输入点存在重复坐标。规则网格里很多点坐标完全一致，去重这一步漏了之后，外接圆判定里出现denom=0的情况，导致返回false，三角形被错误保留，网格拓扑就乱了。

处理后，我把去重容差设成相对值而非绝对值。如果点云坐标单位是米，容差1e-9米基本就是共点；但如果坐标是经纬度，1e-9度对应约0.1毫米，尺度完全不同。相对容差用包围盒对角线长度做基准，自动适配各种坐标系，实测下来稳很多。

另一个坑是共线点。当四个点近似共圆时，外接圆判定的结果对浮点误差非常敏感。比如一个正方形的四个顶点，任意三个点的外接圆都经过第四个点，这时判定结果可能是true也可能是false，取决于浮点舍入。这个问题没有一个绝对正确的解法，工程上常用做法是加一个极小的随机扰动打破对称性，或者在判定时使用容差：

return dist <= r2 + 1e-12 * std::max(1.0, r2);

这个容差让接近圆边的点倾向于被认为在圆外，能显著减少边翻转抖动。代价是最终剖分可能略偏离理论最优，但实际工程中这点偏差完全可接受。

4.2 性能瓶颈与空间索引优化

用Naive的Bowyer-Watson实现测过一组10万个随机点，构建时间是7秒多。这性能离“能用”还有距离，瓶颈在于每插入一个点都要遍历当前所有三角形判断外接圆。

优化的第一步是理解期望复杂度。虽然理论期望是O(n log n)，但前提是插入点与查找到的三角形在空间上局部相关。朴素的遍历方式让每次查找退化成O(n)，整体退化到O(n^2)。

我实现了一个简单有效的空间索引：把当前三角形按重心所在栅格分桶。每个栅格记录落在该区域的三角形索引列表。插入新点时，先计算它所在的栅格，只遍历附近5x5栅格内的三角形，找到包含该点的三角形作为起始搜索位置，然后利用三角形邻接关系做行走搜索（Walk）。

// 栅格分桶初始化 std::unordered_map<int, std::vector<int>> grid; // 插入新点时： Point2D p = pts[i]; int gx = static_cast<int>((p.x - minX) / cellSize); int gy = static_cast<int>((p.y - minY) / cellSize); for (int dx = -2; dx <= 2; dx++) { for (int dy = -2; dy <= 2; dy++) { auto it = grid.find(gridKey(gx + dx, gy + dy)); // 遍历候选三角形，找到包含p的 } }

加了栅格索引后，10万个点降到1.2秒左右。后续我又把当前插入点附近的栅格三角形按距离排序，优先遍历重心近的三角形，进一步缩短了行走搜索的起点寻址时间。这个优化的经验是：不要一上来就上复杂的四叉树或R树，栅格化简单高效，已经能覆盖绝大多数的性能需求。

4.3 可视化联调：用OpenCV画三角网的经验

算法写完后，验证结果的最好方式是可视化。我用了OpenCV作为绘图工具，把三角网直接画出来看效果。

代码思路很简单：遍历三角形，用cv::line把三条边画出来，点的坐标从归一化空间映射回原始坐标再画。这个过程中我犯过一个方向错误：OpenCV的y轴向下，而数学坐标系y轴向上，画出来的三角网上下颠倒。排查时我还以为是剖分逻辑错了，折腾了半小时才发现是坐标变换的问题。从那以后我统一在类内部维护一个坐标变换标记，画图时显式做一次y翻转。

OpenCV还有一个值得注意的点是cv::polylines的闭合参数。画三角形边界时如果用polylines画不闭合路径，会少一条边。我直接改用三次cv::line逐条画，反而更清晰。此外，调试时可以在每个三角形重心画一个圆点，配合txt标注三角形索引，能非常直观地看出坏三角形删除和新三角形生成的过程。

5. 扩展方向：从三角网到更大的应用图景

把Delaunay三角网做出来后，往任何方向扩展都顺理成章。

最直接的扩展是生成Voronoi图。Delaunay三角网和Voronoi图是对偶关系：三角网的每个三角形对应一个Voronoi顶点，每条边对应一条Voronoi边。你只需要遍历三角形，求外接圆圆心，然后连接相邻三角形的外心，就得到完整的Voronoi图。我在地理邻近分析项目里就用过这个，求“每个网点最近的设施点”简直不要太方便。

第二个方向是地形TIN生成。拿到LIDAR点云后，用Delaunay三角网生成TIN模型，然后按地形特征（坡度、粗糙度）做三角形折叠简化，可以在保持地形精度的前提下把三角形数量降到原来的十分之一。这部分工作里，Delaunay三角网的邻接信息帮了大忙，它可以快速判断哪些三角形的简化代价最低。

还有一个方向是跟OpenCV结合做边缘检测辅助。我试过先对图像提取特征点，然后对特征点做Delaunay剖分，用三角形密度分布来判断纹理密集区域。效果虽然不如深度学习方法，但在传统图像处理流程里算是一个轻量且可解释性强的补充手段。

个人经验是，别急着在第一次实现时就堆太多优化和扩展，先把最核心的Bowyer-Watson循环跑通，把外接圆判定的数学搞扎实，后面每一步扩展都能踩在前一步的稳定地基上。我自己的第一版代码里，主循环不过60行，但我在上面迭代了整整两周才把所有边界情况处理干净。算法这种东西，代码只是表象，边界情况的处理和取舍才是真正的功力所在。

本文还有配套的精品资源，点击获取