简介本资源是一套基于Bowyer-Watson算法实现Delaunay三角剖分与Voronoi图生成的Matlab完整仿真方案面向本科及硕士阶段的科研学习者适用于计算几何、空间分析、地理信息系统、路径规划及图像处理等方向的基础建模需求。压缩包共含4个文件3张结果可视化PNG图1个核心m脚本总大小仅92KB轻量易部署适配Matlab 2014a/2019a环境内含可直接运行的代码及对应效果图便于理解算法流程与几何结构关系。目前已有132人下载学习资源聚焦经典计算几何问题提供从点集输入、增量构网、外接圆判据到对偶图转换的全流程实现代码注释清晰、逻辑分层明确特别适合初学者掌握Delaunay/Voronoi内在关联与Matlab向量化编程技巧。 做网格生成和空间分析这几年Delaunay三角网络和Voronoi泰森多边形一直是我绕不开的两个东西。最近整理旧项目时翻到一套用MATLAB写的代码基于Bowyer-Watson算法从零构建这两个结构代码虽然不算长但把计算几何里很多经典细节都踩了一遍。如果你也在做点集剖分、GIS邻近分析、区域划分或者只是想在MATLAB里彻底搞懂Delaunay/Voronoi的底层逻辑这篇应该能帮你省下不少折腾时间。1. 项目背景与核心思路Delaunay和Voronoi到底在算什么1.1 这两个几何结构为什么总是一起出现Delaunay三角网络和Voronoi泰森多边形本质上是一对对偶结构。Delaunay三角网满足一个很关键的性质每个三角形的外接圆内不包含其他点称为空外接圆性质。而Voronoi多边形是把平面划分成若干区域使得每个区域内的点到该区域控制点的距离最近。这两个结构的对偶关系体现在Voronoi图的每条边恰好对应Delaunay三角网的一条边而且是垂直平分的关系Voronoi图的每个顶点恰好是某个Delaunay三角形的外心。这个对偶关系在实际项目里非常有用。比如做无线基站覆盖分析你想知道每个基站管辖哪片区域直接算Voronoi做地形建模、有限元网格划分你需要一个不产生细长三角形的剖分直接上Delaunay。很多场景下算出了Delaunay基本就等于拿到了Voronoi反过来也一样。这也是为什么我在项目里把这两个东西放在一起实现一次构建两边受益。1.2 Bowyer-Watson算法的三步核心Bowyer-Watson算法是逐点插入法里最经典的一种核心思路可以拆成三步先构造一个足够大的超三角形把点集中所有点都包含进去保证初始状态下存在一个可以合法插入点的基础三角网。逐个插入点P找到所有外接圆包含P的三角形这些三角形被称为坏三角形。把它们全部删除留下一个多边形的空洞这个空洞叫影响空腔。用P和空腔边界上的所有边逐条连接生成新的三角形恢复三角网结构。循环执行第2、3步直到所有点插入完毕最后再把所有包含超三角形顶点的三角形删除。整个过程看起来不复杂但每一步都有隐含的边界情况超三角形取多大才不会影响最终结果坏三角形的判断用什么精度空腔边界如何提取才算干净去重和退化怎么处理。第四章里我会逐段拆MATLAB代码把这些问题都摊开来说。1.3 从Delaunay到Voronoi的映射关系因为Delaunay和Voronoi对偶所以生成Voronoi其实不需要重新跑算法。我采用的做法是先完整构建Delaunay三角网然后对每个三角形求外心。这样一来每个Delaunay三角形对应一个Voronoi顶点每条Delaunay内部边对应一条Voronoi边每个Delaunay点对应一个Voronoi胞元。实际操作时只需要遍历三角形列表计算外心再根据三角形之间的邻接关系把共享同一条边的两个三角形的外心连起来就是一条Voronoi边。构建每个点的胞元时找出所有以该点为顶点的三角形把这些三角形的外心按角度排序并连线就得到了该点对应的泰森多边形。这个映射关系看起来很直接但真正让代码“跑得稳”的细节全在外心和邻接关系的计算上下文会详细讲。2. 动手写代码前必须先做的三个设计决策2.1 点集和三角形用什么数据结构最顺手写MATLAB代码前首先要把数据结构定下来否则代码一长就会乱。我的做法是点集用N×2的矩阵P存储每行是一个点的x、y坐标坐标编号从1到N。三角形用一个Ntri×3的矩阵triList存储每一行存三个点的索引代表一个三角形。注意这里存的是索引不是坐标因为索引做去重、排序、邻接查找时都很快。为什么不用cell数组存点坐标因为索引操作可以避免每次创建新数组时的拷贝开销尤其在逐点插入过程中坏三角形的删除和新三角形的添加会频繁发生。我用triList矩阵配合逻辑索引能够用一行代码过滤掉坏三角形非常符合MATLAB的向量化风格。另外还有一个用途最终删除超三角形时只需要检查三角形三个顶点索引是否大于N就能快速筛掉比比较坐标靠谱得多。2.2 超三角形选多大才算“足够大”超三角形是整个Bowyer-Watson算法的启动条件。如果它不够大会有两种情况要么是某些点落在超三角形外部导致初始剖分不合法要么是三角形外接圆覆盖范围不够插入点时会漏掉部分坏三角形。这两种情况都会让最终三角网出现撕裂或重复三角形。我习惯这么取超三角形顶点先计算点集的包围盒取中心点c和包围盒对角线长度d然后以c为中心构造一个边长约为10倍d的等边三角形。把三个顶点放在离点群足够远的位置确保所有点包括它们的极端外接圆都被包含在内。需要提醒的是等边三角形优于直角或钝角三角形因为它的外接圆半径和边长比值较小能减少后续计算外接圆时的浮点误差。超三角形也不是越大越好。太大了会让初始三角形外接圆半径变得非常大判断坏三角形时就会出现大数吃小数的问题导致外心坐标误差明显增大。我实测下来10倍对角线是个比较稳的经验值既不干扰最终结果也不至于让数值精度失控。2.3 外接圆判断的几何式与浮点精度Bowyer-Watson的坏三角形判断核心是判断一个点P是否落在三角形ABC的外接圆内部。求外心坐标最直接的方法是解线性方程组也可以用下列公式function [cx, cy, r2] circumcircle(A, B, C) % A、B、C是两行或三行坐标点这里按2D点处理 d 2 * (A(1)*(B(2)-C(2)) B(1)*(C(2)-A(2)) C(1)*(A(2)-B(2))); if abs(d) 1e-12 cx inf; cy inf; r2 inf; return; end ax A(1); ay A(2); bx B(1); by B(2); cx_ C(1); cy_ C(2); a2 ax*ax ay*ay; b2 bx*bx by*by; c2 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; cx ux; cy uy; r2 (ax-ux)^2 (ay-uy)^2; end判断P是否在圆内时不能直接用dist(P, center) r因为浮点运算中边界情况极多。我统一用平方距离比较dist2 r2 tol其中tol取一个与点集尺度相关的值比如点集包围盒对角线长度的1e-9倍。引入容差后位于圆上的点也能稳定地进入删除流程避免因为0.0000001的误差导致插入失败这是我在实际调试中踩过的最典型的坑。3. MATLAB核心实现Bowyer-Watson算法逐段拆解3.1 初始化读入点集、去重、构造超三角形写了一个主函数delaunay_voronoi_demo.m整个流程从随机生成测试点开始。实际项目里可以直接替换成自己的坐标点。% 生成随机测试点稍微加点簇状分布更接近真实场景 rng(42); N 60; P [randn(N,1)*0.3 0.5, randn(N,1)*0.3 0.5]; P [P; rand(N,1)*3 1, rand(N,1)*3 1]; P unique(round(P, 10), rows, stable); N size(P, 1);这里unique去重很关键。如果你的数据里存在重复点Bowyer-Watson会在插入重复点时产生零面积三角形进而让外接圆计算出现分母接近0的情况后面所有判断都会失效。round(P,10)是为了避免浮点运算导致“同一坐标但差一位小数”的点没被当成重复点。构造超三角形我单独抽了一个子函数保持主流程清晰function T superTriangle(P) xmin min(P(:,1)); xmax max(P(:,1)); ymin min(P(:,2)); ymax max(P(:,2)); cx (xmin xmax) / 2; cy (ymin ymax) / 2; d sqrt((xmax - xmin)^2 (ymax - ymin)^2); R 10 * d; T [cx, cy 2*R; cx - sqrt(3)*R, cy - R; cx sqrt(3)*R, cy - R]; end用中心点和固定半径构造等边三角形是为了让三个顶点关于点集中心对称分布从源头上减小外心计算的病态程度。这一步不需要太精确但要记好最终结果里包含这三个顶点的三角形都要删掉。3.2 主循环逐点插入、找坏三角形、重建局部三角网主循环是整个算法的核心我把它拆成了几个清晰的步骤避免在一大坨代码里绕晕。% 初始化把超三角形顶点也加入坐标矩阵统一编号 T superTriangle(P); P_all [P; T]; triList [N1, N2, N3]; % 初始三角形 for i 1:N p P(i,:); % 1. 找出所有外接圆包含p的三角形 [cx, cy, r2] arrayfun((t) triangleCircum(P_all, triList(t,:)), 1:size(triList,1)); dist2 (cx - p(1)).^2 (cy - p(2)).^2; tol 1e-9 * (max(P_all(:,1)) - min(P_all(:,1))); badIdx find(dist2 r2 tol); if isempty(badIdx) error(当前点未落在任何三角形外接圆内超三角形可能不够大); end badTri triList(badIdx, :); triList(badIdx, :) []; % 2. 提取影响空腔边界统计每条边出现的次数 edges [badTri(:,1) badTri(:,2); badTri(:,2) badTri(:,3); badTri(:,3) badTri(:,1)]; edges sort(edges, 2); [uniqueEdges, ~, ic] unique(edges, rows, stable); edgeCount accumarray(ic, 1); boundaryEdges uniqueEdges(edgeCount 1, :); % 3. 连接新点与边界边生成新三角形 newTris [boundaryEdges(:,1), boundaryEdges(:,2), repmat(i, size(boundaryEdges,1), 1)]; triList [triList; newTris]; end这段代码里最值得讲的是边界提取。删除坏三角形后留下的空腔边界本质上是这些坏三角形里出现了奇数次的边。出现两次的边就是两个坏三角形共享的内部边删除三角形后它不应该再存在只出现一次的边才是空腔的真正外边界。我用sort把边的顶点排序再用unique配合accumarray统计每条边的出现次数一行代码就把内部边和边界边分离开了。这个方法比逐个邻接遍历高效得多而且在MATLAB里属于标准操作。循环结束后还差一步删除所有包含超三角形顶点的三角形。triList(any(triList N, 2), :) [];超过N的索引全部是超三角形的三个顶点之一直接删掉即可。到这一步一个不依赖MATLAB任何内置剖分函数的Delaunay三角网就构建完成了。3.3 去重和边界检查别让重复三角形毁了结果很多人写Bowyer-Watson时会遇到最后结果里出现重复三角形或者三角形三条边共线的情况这通常不是算法逻辑问题而是点集去重不彻底、容差设置不合适导致的。我在循环前已经对原始点做了unique去重但插入过程中因为容差原因仍然有可能生成面积接近0的三角形。所以在三角形生成后我加了一个兜底检查过滤掉面积小于某阈值的三角形。A P_all(triList(:,1), :); B P_all(triList(:,2), :); C P_all(triList(:,3), :); area2 abs((B(:,1)-A(:,1)).*(C(:,2)-A(:,2)) - (B(:,2)-A(:,2)).*(C(:,1)-A(:,1))); triList(area2 1e-12, :) [];这个过滤不能随便删因为面积过小的三角形虽然不影响后续Voronoi的生成但会在可视化时画出一堆肉眼几乎看不见的碎边影响图形输出的美观度。更重要的是它会降低后续邻接查找的效率。3.4 生成Voronoi从三角形外心到泰森多边形Delaunay三角网输出得到后Voronoi就省力了。先给每个三角形算外心再聚合到每个原始点上。% 计算每个Delaunay三角形外心 m size(triList, 1); vorVertices zeros(m, 2); for t 1:m [vorVertices(t,1), vorVertices(t,2)] triangleCircum(P_all, triList(t,:)); end如果只需要Voronoi边可以直接遍历每条内部Delaunay边找到共享边的两个三角形把两个外心连线。但我在项目里更多是需要“每个点的泰森多边形”所以我按点聚合所有邻接三角形的外心vorCells cell(N, 1); for i 1:N triIdx find(any(triList i, 2)); if isempty(triIdx) continue; end cellVerts vorVertices(triIdx, :); % 按角度排序保证多边形顶点顺序正确 angles atan2(cellVerts(:,2) - mean(cellVerts(:,2)), cellVerts(:,1) - mean(cellVerts(:,1))); [~, order] sort(angles); vorCells{i} cellVerts(order, :); end这里有一点要注意Voronoi胞元在凸包边界处是无限区域外心可能跑到无穷远。我在可视化时会把边界外的点裁剪到画布范围内MATLAB的axis命令虽然可以显示但裁剪操作最好自己做否则连线会画出很长的斜线。第四章可视化部分我会演示如何既保留边界射线的效果又不破坏画面。4. 实操记录从零跑通整个流程4.1 完整的函数组织方式写代码的时候我把功能拆成三个文件主脚本、核心算法函数、工具函数这样调试和维护都方便。主脚本负责生成点、调用算法、画图核心算法函数delaunayByBowyerWatson负责构建三角网工具函数circumcircle和superTriangle放在同一个文件末尾作为嵌套函数或者单独拆出来都行。如果你的点集有几个万级别逐点插入配for循环在MATLAB里会明显变慢。这时候有两个方向一是把坏三角形查找向量化像我在第三章里那样一次性计算所有三角形外接圆而不是再套一层for二是考虑加入网格空间索引减少每次插入时扫描的三角形数量。对于几千点的规模来说向量化已经足够我的实测结果是在普通笔记本上60个点可以瞬间完成2000个点大约需要0.3到1秒完全可以用在项目调试中。4.2 可视化把Delaunay和Voronoi画在一张图上画图是检验算法正确性最直观的手段。我习惯把Delaunay三角网用浅色线画出来Voronoi多边形用另一条颜色画在同一张图上能很清楚地看到对偶关系。figure; hold on; % 画Delaunay三角网 triplot(triList, P_all(:,1), P_all(:,2), Color, [0.6 0.6 0.6], LineWidth, 0.5); % 画Voronoi for i 1:N if ~isempty(vorCells{i}) patch(vorCells{i}(:,1), vorCells{i}(:,2), w, EdgeColor, [0.8 0.2 0.2], LineWidth, 1.2); end end plot(P(:,1), P(:,2), ko, MarkerFaceColor, k, MarkerSize, 4); axis equal;一个非常容易踩的坑是axis equal。如果没有等比例坐标轴Voronoi的多边形形状会被拉伸看起来像是算法错误实际上只是显示比例问题。我第一次排查了半天最后发现是坐标轴没等比例白白浪费了一个小时。4.3 对比验证用MATLAB内置函数当“裁判”自己写的算法对不对最直接的办法是拿MATLAB自带的delaunayTriangulation做交叉验证。虽然项目目的是练习底层实现但验证环节不能省。DT delaunayTriangulation(P); triBuiltin DT.ConnectivityList; % 比较三角形数量是否一致 fprintf(自定义 Delaunay 三角形数量: %d\n, size(triList, 1)); fprintf(内置 Delaunay 三角形数量: %d\n, size(triBuiltin, 1));注意三角形数量一致不代表两个剖分完全一样顶点索引顺序和三角形集合可能有微小差异这是Delaunay剖分不唯一导致的正常情况。更实际的验证方式是检查空外接圆性质和边界数量所有三角形外接圆内不含其它点凸包边界上的边数量等于凸包顶点数。我写了一个自动化检查脚本每个点插入后都断言一下当前三角网的三角形数量满足欧拉公式这样能尽早暴露问题而不是等最后画图才发现。5. 典型问题与排查实录5.1 外接圆误判浮点容差到底怎么设置才合理这是整个项目里我踩得最深的坑。一开始我把容差写成了固定值1e-10结果点集坐标范围比较大的时候比如坐标在10000量级1e-10完全被浮点精度淹没导致一部分恰好在圆上的点没有被判定为坏三角形最后生成Delaunay网格时出现明显的孔洞。后来我把容差改成相对值与点集包围盒对角线长度挂钩tol 1e-9 * diagLen。diagLen是点集范围的量级这样做的好处是无论你的数据是单位正方形内的小规模坐标还是几十万坐标值的地理坐标容差都能跟着缩放。建议读者在实现时一定不要用固定绝对值容差这算是一条通用经验。5.2 四点共圆退化情况下会出现非唯一剖分当四个点恰好共圆时Delaunay剖分并不唯一。Bowyer-Watson在这种情况下会随机选择两条对角线中的一条不能说算法错误但可能和你预期的结果不一致而且有时会在相邻位置产生一个极其扁平的三角形看起来像是bug。我处理这类问题的方式是加入微小扰动或者接受非唯一剖分。如果点集来自实际数据精确保存了原始坐标那就接受算法给出的结果如果点集是人工生成的比如规则网格点我会在预处理阶段给每个点加一个约1e-8量级的随机扰动打破共圆状态保证剖分结果的稳定性。这个技巧在有限元网格生成里也常用不只是为了“好看”。5.3 凸包边界上的Voronoi射线怎么处理理论上凸包边界边的Voronoi边是一条射线延伸到无穷远。但在屏幕上画图或者做区域统计时不能真的画到无穷远需要做一个裁剪。我的做法是直接取所有Voronoi胞元顶点如果某个顶点坐标超过画布范围的1.5倍就把它拉回到画布边缘然后再连线。这样做虽然牺牲了一点点几何精确度但对视觉呈现和后续的区域面积统计影响很小。如果你需要精确的无限区域Voronoi建议用polyshape和intersect做裁剪或者直接参考CGAL等成熟库的思路。5.4 算法速度慢当点集规模达到5000以上时怎么办逐点插入的Bowyer-Watson算法如果不加任何优化复杂度在点集规模增大后会明显变差。我在测试2000个点的时候还能接受到5000个点时每次插入都要扫描全部现有三角形耗时开始变得明显。如果遇到大数据集我的建议不是在MATLAB里硬堆代码优化而是换一个思路先用delaunayTriangulation内置函数跑通业务流程再用自己写的版本做小规模验证和教学。或者对算法加入空间索引比如把平面分成格网插入点P时只检查P所在格网及其邻接格网内的三角形这样能把扫描范围缩小一个量级。这两种方法我都试过后者的编码量确实比较大所以除非是性能敏感场景否则没必要一开始就做。一路实现下来我最大的感受是Bowyer-Watson本质上不复杂真正考验人的地方全在边界处理和数据结构的组织上。超三角形大小、容差设置、边界边提取方式、Voronoi顶点裁剪每一个细节都可能让最终结果从“看起来合理”变成“严格正确”。如果你打算在MATLAB里复现这套流程建议把本文第三章的核心代码完整跑一遍再结合第五章的踩坑清单做排查。等你自己亲手调通一遍这两个计算几何里的经典结构就会真正变成你工具箱里的东西而不是只会调包的黑盒。本文还有配套的精品资源点击获取 MATLAB实现Delaunay三角网与Voronoi图：Bowyer-Watson算法详解
简介本资源是一套基于Bowyer-Watson算法实现Delaunay三角剖分与Voronoi图生成的Matlab完整仿真方案面向本科及硕士阶段的科研学习者适用于计算几何、空间分析、地理信息系统、路径规划及图像处理等方向的基础建模需求。压缩包共含4个文件3张结果可视化PNG图1个核心m脚本总大小仅92KB轻量易部署适配Matlab 2014a/2019a环境内含可直接运行的代码及对应效果图便于理解算法流程与几何结构关系。目前已有132人下载学习资源聚焦经典计算几何问题提供从点集输入、增量构网、外接圆判据到对偶图转换的全流程实现代码注释清晰、逻辑分层明确特别适合初学者掌握Delaunay/Voronoi内在关联与Matlab向量化编程技巧。 做网格生成和空间分析这几年Delaunay三角网络和Voronoi泰森多边形一直是我绕不开的两个东西。最近整理旧项目时翻到一套用MATLAB写的代码基于Bowyer-Watson算法从零构建这两个结构代码虽然不算长但把计算几何里很多经典细节都踩了一遍。如果你也在做点集剖分、GIS邻近分析、区域划分或者只是想在MATLAB里彻底搞懂Delaunay/Voronoi的底层逻辑这篇应该能帮你省下不少折腾时间。1. 项目背景与核心思路Delaunay和Voronoi到底在算什么1.1 这两个几何结构为什么总是一起出现Delaunay三角网络和Voronoi泰森多边形本质上是一对对偶结构。Delaunay三角网满足一个很关键的性质每个三角形的外接圆内不包含其他点称为空外接圆性质。而Voronoi多边形是把平面划分成若干区域使得每个区域内的点到该区域控制点的距离最近。这两个结构的对偶关系体现在Voronoi图的每条边恰好对应Delaunay三角网的一条边而且是垂直平分的关系Voronoi图的每个顶点恰好是某个Delaunay三角形的外心。这个对偶关系在实际项目里非常有用。比如做无线基站覆盖分析你想知道每个基站管辖哪片区域直接算Voronoi做地形建模、有限元网格划分你需要一个不产生细长三角形的剖分直接上Delaunay。很多场景下算出了Delaunay基本就等于拿到了Voronoi反过来也一样。这也是为什么我在项目里把这两个东西放在一起实现一次构建两边受益。1.2 Bowyer-Watson算法的三步核心Bowyer-Watson算法是逐点插入法里最经典的一种核心思路可以拆成三步先构造一个足够大的超三角形把点集中所有点都包含进去保证初始状态下存在一个可以合法插入点的基础三角网。逐个插入点P找到所有外接圆包含P的三角形这些三角形被称为坏三角形。把它们全部删除留下一个多边形的空洞这个空洞叫影响空腔。用P和空腔边界上的所有边逐条连接生成新的三角形恢复三角网结构。循环执行第2、3步直到所有点插入完毕最后再把所有包含超三角形顶点的三角形删除。整个过程看起来不复杂但每一步都有隐含的边界情况超三角形取多大才不会影响最终结果坏三角形的判断用什么精度空腔边界如何提取才算干净去重和退化怎么处理。第四章里我会逐段拆MATLAB代码把这些问题都摊开来说。1.3 从Delaunay到Voronoi的映射关系因为Delaunay和Voronoi对偶所以生成Voronoi其实不需要重新跑算法。我采用的做法是先完整构建Delaunay三角网然后对每个三角形求外心。这样一来每个Delaunay三角形对应一个Voronoi顶点每条Delaunay内部边对应一条Voronoi边每个Delaunay点对应一个Voronoi胞元。实际操作时只需要遍历三角形列表计算外心再根据三角形之间的邻接关系把共享同一条边的两个三角形的外心连起来就是一条Voronoi边。构建每个点的胞元时找出所有以该点为顶点的三角形把这些三角形的外心按角度排序并连线就得到了该点对应的泰森多边形。这个映射关系看起来很直接但真正让代码“跑得稳”的细节全在外心和邻接关系的计算上下文会详细讲。2. 动手写代码前必须先做的三个设计决策2.1 点集和三角形用什么数据结构最顺手写MATLAB代码前首先要把数据结构定下来否则代码一长就会乱。我的做法是点集用N×2的矩阵P存储每行是一个点的x、y坐标坐标编号从1到N。三角形用一个Ntri×3的矩阵triList存储每一行存三个点的索引代表一个三角形。注意这里存的是索引不是坐标因为索引做去重、排序、邻接查找时都很快。为什么不用cell数组存点坐标因为索引操作可以避免每次创建新数组时的拷贝开销尤其在逐点插入过程中坏三角形的删除和新三角形的添加会频繁发生。我用triList矩阵配合逻辑索引能够用一行代码过滤掉坏三角形非常符合MATLAB的向量化风格。另外还有一个用途最终删除超三角形时只需要检查三角形三个顶点索引是否大于N就能快速筛掉比比较坐标靠谱得多。2.2 超三角形选多大才算“足够大”超三角形是整个Bowyer-Watson算法的启动条件。如果它不够大会有两种情况要么是某些点落在超三角形外部导致初始剖分不合法要么是三角形外接圆覆盖范围不够插入点时会漏掉部分坏三角形。这两种情况都会让最终三角网出现撕裂或重复三角形。我习惯这么取超三角形顶点先计算点集的包围盒取中心点c和包围盒对角线长度d然后以c为中心构造一个边长约为10倍d的等边三角形。把三个顶点放在离点群足够远的位置确保所有点包括它们的极端外接圆都被包含在内。需要提醒的是等边三角形优于直角或钝角三角形因为它的外接圆半径和边长比值较小能减少后续计算外接圆时的浮点误差。超三角形也不是越大越好。太大了会让初始三角形外接圆半径变得非常大判断坏三角形时就会出现大数吃小数的问题导致外心坐标误差明显增大。我实测下来10倍对角线是个比较稳的经验值既不干扰最终结果也不至于让数值精度失控。2.3 外接圆判断的几何式与浮点精度Bowyer-Watson的坏三角形判断核心是判断一个点P是否落在三角形ABC的外接圆内部。求外心坐标最直接的方法是解线性方程组也可以用下列公式function [cx, cy, r2] circumcircle(A, B, C) % A、B、C是两行或三行坐标点这里按2D点处理 d 2 * (A(1)*(B(2)-C(2)) B(1)*(C(2)-A(2)) C(1)*(A(2)-B(2))); if abs(d) 1e-12 cx inf; cy inf; r2 inf; return; end ax A(1); ay A(2); bx B(1); by B(2); cx_ C(1); cy_ C(2); a2 ax*ax ay*ay; b2 bx*bx by*by; c2 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; cx ux; cy uy; r2 (ax-ux)^2 (ay-uy)^2; end判断P是否在圆内时不能直接用dist(P, center) r因为浮点运算中边界情况极多。我统一用平方距离比较dist2 r2 tol其中tol取一个与点集尺度相关的值比如点集包围盒对角线长度的1e-9倍。引入容差后位于圆上的点也能稳定地进入删除流程避免因为0.0000001的误差导致插入失败这是我在实际调试中踩过的最典型的坑。3. MATLAB核心实现Bowyer-Watson算法逐段拆解3.1 初始化读入点集、去重、构造超三角形写了一个主函数delaunay_voronoi_demo.m整个流程从随机生成测试点开始。实际项目里可以直接替换成自己的坐标点。% 生成随机测试点稍微加点簇状分布更接近真实场景 rng(42); N 60; P [randn(N,1)*0.3 0.5, randn(N,1)*0.3 0.5]; P [P; rand(N,1)*3 1, rand(N,1)*3 1]; P unique(round(P, 10), rows, stable); N size(P, 1);这里unique去重很关键。如果你的数据里存在重复点Bowyer-Watson会在插入重复点时产生零面积三角形进而让外接圆计算出现分母接近0的情况后面所有判断都会失效。round(P,10)是为了避免浮点运算导致“同一坐标但差一位小数”的点没被当成重复点。构造超三角形我单独抽了一个子函数保持主流程清晰function T superTriangle(P) xmin min(P(:,1)); xmax max(P(:,1)); ymin min(P(:,2)); ymax max(P(:,2)); cx (xmin xmax) / 2; cy (ymin ymax) / 2; d sqrt((xmax - xmin)^2 (ymax - ymin)^2); R 10 * d; T [cx, cy 2*R; cx - sqrt(3)*R, cy - R; cx sqrt(3)*R, cy - R]; end用中心点和固定半径构造等边三角形是为了让三个顶点关于点集中心对称分布从源头上减小外心计算的病态程度。这一步不需要太精确但要记好最终结果里包含这三个顶点的三角形都要删掉。3.2 主循环逐点插入、找坏三角形、重建局部三角网主循环是整个算法的核心我把它拆成了几个清晰的步骤避免在一大坨代码里绕晕。% 初始化把超三角形顶点也加入坐标矩阵统一编号 T superTriangle(P); P_all [P; T]; triList [N1, N2, N3]; % 初始三角形 for i 1:N p P(i,:); % 1. 找出所有外接圆包含p的三角形 [cx, cy, r2] arrayfun((t) triangleCircum(P_all, triList(t,:)), 1:size(triList,1)); dist2 (cx - p(1)).^2 (cy - p(2)).^2; tol 1e-9 * (max(P_all(:,1)) - min(P_all(:,1))); badIdx find(dist2 r2 tol); if isempty(badIdx) error(当前点未落在任何三角形外接圆内超三角形可能不够大); end badTri triList(badIdx, :); triList(badIdx, :) []; % 2. 提取影响空腔边界统计每条边出现的次数 edges [badTri(:,1) badTri(:,2); badTri(:,2) badTri(:,3); badTri(:,3) badTri(:,1)]; edges sort(edges, 2); [uniqueEdges, ~, ic] unique(edges, rows, stable); edgeCount accumarray(ic, 1); boundaryEdges uniqueEdges(edgeCount 1, :); % 3. 连接新点与边界边生成新三角形 newTris [boundaryEdges(:,1), boundaryEdges(:,2), repmat(i, size(boundaryEdges,1), 1)]; triList [triList; newTris]; end这段代码里最值得讲的是边界提取。删除坏三角形后留下的空腔边界本质上是这些坏三角形里出现了奇数次的边。出现两次的边就是两个坏三角形共享的内部边删除三角形后它不应该再存在只出现一次的边才是空腔的真正外边界。我用sort把边的顶点排序再用unique配合accumarray统计每条边的出现次数一行代码就把内部边和边界边分离开了。这个方法比逐个邻接遍历高效得多而且在MATLAB里属于标准操作。循环结束后还差一步删除所有包含超三角形顶点的三角形。triList(any(triList N, 2), :) [];超过N的索引全部是超三角形的三个顶点之一直接删掉即可。到这一步一个不依赖MATLAB任何内置剖分函数的Delaunay三角网就构建完成了。3.3 去重和边界检查别让重复三角形毁了结果很多人写Bowyer-Watson时会遇到最后结果里出现重复三角形或者三角形三条边共线的情况这通常不是算法逻辑问题而是点集去重不彻底、容差设置不合适导致的。我在循环前已经对原始点做了unique去重但插入过程中因为容差原因仍然有可能生成面积接近0的三角形。所以在三角形生成后我加了一个兜底检查过滤掉面积小于某阈值的三角形。A P_all(triList(:,1), :); B P_all(triList(:,2), :); C P_all(triList(:,3), :); area2 abs((B(:,1)-A(:,1)).*(C(:,2)-A(:,2)) - (B(:,2)-A(:,2)).*(C(:,1)-A(:,1))); triList(area2 1e-12, :) [];这个过滤不能随便删因为面积过小的三角形虽然不影响后续Voronoi的生成但会在可视化时画出一堆肉眼几乎看不见的碎边影响图形输出的美观度。更重要的是它会降低后续邻接查找的效率。3.4 生成Voronoi从三角形外心到泰森多边形Delaunay三角网输出得到后Voronoi就省力了。先给每个三角形算外心再聚合到每个原始点上。% 计算每个Delaunay三角形外心 m size(triList, 1); vorVertices zeros(m, 2); for t 1:m [vorVertices(t,1), vorVertices(t,2)] triangleCircum(P_all, triList(t,:)); end如果只需要Voronoi边可以直接遍历每条内部Delaunay边找到共享边的两个三角形把两个外心连线。但我在项目里更多是需要“每个点的泰森多边形”所以我按点聚合所有邻接三角形的外心vorCells cell(N, 1); for i 1:N triIdx find(any(triList i, 2)); if isempty(triIdx) continue; end cellVerts vorVertices(triIdx, :); % 按角度排序保证多边形顶点顺序正确 angles atan2(cellVerts(:,2) - mean(cellVerts(:,2)), cellVerts(:,1) - mean(cellVerts(:,1))); [~, order] sort(angles); vorCells{i} cellVerts(order, :); end这里有一点要注意Voronoi胞元在凸包边界处是无限区域外心可能跑到无穷远。我在可视化时会把边界外的点裁剪到画布范围内MATLAB的axis命令虽然可以显示但裁剪操作最好自己做否则连线会画出很长的斜线。第四章可视化部分我会演示如何既保留边界射线的效果又不破坏画面。4. 实操记录从零跑通整个流程4.1 完整的函数组织方式写代码的时候我把功能拆成三个文件主脚本、核心算法函数、工具函数这样调试和维护都方便。主脚本负责生成点、调用算法、画图核心算法函数delaunayByBowyerWatson负责构建三角网工具函数circumcircle和superTriangle放在同一个文件末尾作为嵌套函数或者单独拆出来都行。如果你的点集有几个万级别逐点插入配for循环在MATLAB里会明显变慢。这时候有两个方向一是把坏三角形查找向量化像我在第三章里那样一次性计算所有三角形外接圆而不是再套一层for二是考虑加入网格空间索引减少每次插入时扫描的三角形数量。对于几千点的规模来说向量化已经足够我的实测结果是在普通笔记本上60个点可以瞬间完成2000个点大约需要0.3到1秒完全可以用在项目调试中。4.2 可视化把Delaunay和Voronoi画在一张图上画图是检验算法正确性最直观的手段。我习惯把Delaunay三角网用浅色线画出来Voronoi多边形用另一条颜色画在同一张图上能很清楚地看到对偶关系。figure; hold on; % 画Delaunay三角网 triplot(triList, P_all(:,1), P_all(:,2), Color, [0.6 0.6 0.6], LineWidth, 0.5); % 画Voronoi for i 1:N if ~isempty(vorCells{i}) patch(vorCells{i}(:,1), vorCells{i}(:,2), w, EdgeColor, [0.8 0.2 0.2], LineWidth, 1.2); end end plot(P(:,1), P(:,2), ko, MarkerFaceColor, k, MarkerSize, 4); axis equal;一个非常容易踩的坑是axis equal。如果没有等比例坐标轴Voronoi的多边形形状会被拉伸看起来像是算法错误实际上只是显示比例问题。我第一次排查了半天最后发现是坐标轴没等比例白白浪费了一个小时。4.3 对比验证用MATLAB内置函数当“裁判”自己写的算法对不对最直接的办法是拿MATLAB自带的delaunayTriangulation做交叉验证。虽然项目目的是练习底层实现但验证环节不能省。DT delaunayTriangulation(P); triBuiltin DT.ConnectivityList; % 比较三角形数量是否一致 fprintf(自定义 Delaunay 三角形数量: %d\n, size(triList, 1)); fprintf(内置 Delaunay 三角形数量: %d\n, size(triBuiltin, 1));注意三角形数量一致不代表两个剖分完全一样顶点索引顺序和三角形集合可能有微小差异这是Delaunay剖分不唯一导致的正常情况。更实际的验证方式是检查空外接圆性质和边界数量所有三角形外接圆内不含其它点凸包边界上的边数量等于凸包顶点数。我写了一个自动化检查脚本每个点插入后都断言一下当前三角网的三角形数量满足欧拉公式这样能尽早暴露问题而不是等最后画图才发现。5. 典型问题与排查实录5.1 外接圆误判浮点容差到底怎么设置才合理这是整个项目里我踩得最深的坑。一开始我把容差写成了固定值1e-10结果点集坐标范围比较大的时候比如坐标在10000量级1e-10完全被浮点精度淹没导致一部分恰好在圆上的点没有被判定为坏三角形最后生成Delaunay网格时出现明显的孔洞。后来我把容差改成相对值与点集包围盒对角线长度挂钩tol 1e-9 * diagLen。diagLen是点集范围的量级这样做的好处是无论你的数据是单位正方形内的小规模坐标还是几十万坐标值的地理坐标容差都能跟着缩放。建议读者在实现时一定不要用固定绝对值容差这算是一条通用经验。5.2 四点共圆退化情况下会出现非唯一剖分当四个点恰好共圆时Delaunay剖分并不唯一。Bowyer-Watson在这种情况下会随机选择两条对角线中的一条不能说算法错误但可能和你预期的结果不一致而且有时会在相邻位置产生一个极其扁平的三角形看起来像是bug。我处理这类问题的方式是加入微小扰动或者接受非唯一剖分。如果点集来自实际数据精确保存了原始坐标那就接受算法给出的结果如果点集是人工生成的比如规则网格点我会在预处理阶段给每个点加一个约1e-8量级的随机扰动打破共圆状态保证剖分结果的稳定性。这个技巧在有限元网格生成里也常用不只是为了“好看”。5.3 凸包边界上的Voronoi射线怎么处理理论上凸包边界边的Voronoi边是一条射线延伸到无穷远。但在屏幕上画图或者做区域统计时不能真的画到无穷远需要做一个裁剪。我的做法是直接取所有Voronoi胞元顶点如果某个顶点坐标超过画布范围的1.5倍就把它拉回到画布边缘然后再连线。这样做虽然牺牲了一点点几何精确度但对视觉呈现和后续的区域面积统计影响很小。如果你需要精确的无限区域Voronoi建议用polyshape和intersect做裁剪或者直接参考CGAL等成熟库的思路。5.4 算法速度慢当点集规模达到5000以上时怎么办逐点插入的Bowyer-Watson算法如果不加任何优化复杂度在点集规模增大后会明显变差。我在测试2000个点的时候还能接受到5000个点时每次插入都要扫描全部现有三角形耗时开始变得明显。如果遇到大数据集我的建议不是在MATLAB里硬堆代码优化而是换一个思路先用delaunayTriangulation内置函数跑通业务流程再用自己写的版本做小规模验证和教学。或者对算法加入空间索引比如把平面分成格网插入点P时只检查P所在格网及其邻接格网内的三角形这样能把扫描范围缩小一个量级。这两种方法我都试过后者的编码量确实比较大所以除非是性能敏感场景否则没必要一开始就做。一路实现下来我最大的感受是Bowyer-Watson本质上不复杂真正考验人的地方全在边界处理和数据结构的组织上。超三角形大小、容差设置、边界边提取方式、Voronoi顶点裁剪每一个细节都可能让最终结果从“看起来合理”变成“严格正确”。如果你打算在MATLAB里复现这套流程建议把本文第三章的核心代码完整跑一遍再结合第五章的踩坑清单做排查。等你自己亲手调通一遍这两个计算几何里的经典结构就会真正变成你工具箱里的东西而不是只会调包的黑盒。本文还有配套的精品资源点击获取