简介Matlab生成Voronoi图的代码包围绕计算几何中的Voronoi图与Delaunay三角化实现解决从离散点集构造空间划分的核心问题适用于图形学初学者、算法研究者以及需要在Matlab中定制生成逻辑的工程师。压缩包共11个文件包含9个Matlab脚本和2张结果图压缩后仅131KB脚本主要实现三角剖分、外接圆构建与判定、Voronoi边生成等核心环节。代码结构清晰包含Delaunay三角化、外接圆构建、点圆关系判断、相邻三角形查找等功能模块并配有结果图辅助验证。目前已有5921人学习下载这份代码是理解Voronoi图生成过程的实用参考。借助这份代码读者可以逐步分析点集如何通过Delaunay三角化生成Voronoi边界随意修改输入点集进行实验并在此基础上做性能优化或二次开发也适合作为算法课程的教学演示。1. 写在前面:为什么突然聊起Voronoi如果你手里碰巧攒了一批离散的坐标点想把这些点按“最近距离”划分成一个个势力范围Voronoi图就是干这个的。它的应用横跨计算机图形学、路径规划、有限元网格生成、气象站点插值、城市服务区划分等多个方向而Matlab里生成Voronoi图几乎是我用过的所有工具里最快、最省事的一条路。这篇文章就是一份完整的实操笔记。我会从Voronoi图的基本概念入手把Matlab里相关的核心函数voronoi、voronoin、polybuffer等逐个拆开讲再给出一份可以直接抄作业的完整代码。除此之外还会把我实际踩过的坑、排查过的报错、以及性能优化方面的心得一并整理出来。无论你是刚接触Matlab的学生还是做工程仿真的研发人员照着这篇走大概率能少折腾半天。2. 概念先行:Voronoi图到底在算什么2.1 一句话解释VoronoiVoronoi图又叫泰森多边形Thiessen polygons、Dirichlet图。给定平面上一组点称为种子点或生长点Voronoi图会把整个平面划分成若干个区域每个区域内的任意一点到对应种子点的距离都小于到其他任何种子点的距离。打个比方假设你们学校有5个食堂每个学生都去离自己最近的那个食堂。如果按这个规则把整个校园划分成5个片区每个片区就是该食堂的Voronoi区域。食堂就是种子点片区分界线就是Voronoi边。2.2 数学表达与几何性质设种子点集合为P {p1, p2, ..., pn}其中pi (xi, yi)。点p对应的Voronoi区域定义为V(pi) {p | dist(p, pi) ≤ dist(p, pj), ∀j ≠ i}也就是说Voronoi区域是到某个种子点距离最近的所有点的集合。几何上两个相邻种子点之间的Voronoi边界就是这两个点连线的垂直平分线。多个垂直平分线相交形成的顶点称为Voronoi顶点。这些顶点有一个重要性质它是至少三个种子点的外接圆圆心。这个性质在后续做Delaunay三角剖分和有限元网格生成时特别有用。2.3 为什么在Matlab里实现最舒服我这些年用过Python的scipy.spatial.Voronoi也用过C的CGAL库说实话各有优势。但Matlab的特点在于“开箱即用”不需要额外安装第三方库voronoi和voronoin两个函数直接调用即可而且绘图、计算结果的可视化集成在一起对快速验证算法思路来说非常方便。如果做科研绘图Matlab输出的图件质量也足够发表用。3. Matlab生成Voronoi图的核心函数3.1voronoi:可视化优先voronoi函数有两种常见调用方式% 方式一直接传入散点坐标 x rand(1, 20); y rand(1, 20); voronoi(x, y); % 方式二传入坐标矩阵 P [x(:), y(:)]; voronoi(P);这种方式最直观执行后会直接弹出一个Figure窗口绘制出种子点和对应的Voronoi边界线。适合快速查看结果、检查种子点分布情况。但要注意voronoi函数返回的是线段的端点坐标而不是每个区域的顶点坐标。换句话说它是由一组线段构成的一个画面算半个图形学工具。如果你想拿到每个多边形区域的完整顶点信息再用这些顶点做进一步计算比如面积统计、区域裁剪就需要用voronoin。3.2voronoin:数据优先voronoin返回的内容更底层、更有分析价值[V, C] voronoin(P);其中V所有Voronoi顶点的坐标矩阵尺寸为m×2。第一行V(1,:)比较特殊通常表示无穷远点用Inf表示表示这个区域向外延伸到无穷远。C元胞数组C{i}存储第i个种子点对应的Voronoi区域顶点在V中的索引。比如C{3} [5 2 1 6]表示第3个种子点的Voronoi区域是一个四边形顶点依次是V(5,:)、V(2,:)、V(1,:)、V(6,:)。3.3 一个具体的可视化例子为了把两个函数的差别看明白我用20个随机点做个对比rng(42); % 固定随机种子保证结果可复现 P rand(20, 2); figure; subplot(1, 2, 1); voronoi(P(:,1), P(:,2)); title(voronoi函数绘制结果); axis equal; xlim([0 1]); ylim([0 1]); subplot(1, 2, 2); [V, C] voronoin(P); hold on; for i 1:length(C) if all(V(C{i}, :) ~ Inf) % 跳过包含无穷远点的区域 patch(V(C{i}, 1), V(C{i}, 2), rand(1,3), FaceAlpha, 0.3); end end plot(P(:,1), P(:,2), ko, MarkerFaceColor, k); axis equal; xlim([0 1]); ylim([0 1]); title(voronoin函数patch绘制结果);运行这段代码你会直观看到左侧就是常见的Voronoi图右侧我们通过patch填充了各个区域视觉效果更直观。注意右侧代码里我特意跳过了包含Inf的区域因为那些区域向外延伸到了无穷远直接patch会出问题。4. 完整案例:从随机点位到服务区划分4.1 场景与需求描述这里我模拟一个实际场景某城市有12个快递站点每个站点负责最近的片区。假设城市的坐标范围是0到10公里的正方形区域用Voronoi图来划分各站点的配送范围并计算每个片区的面积。这个案例在物流配送、公共服务设施选址中非常典型。4.2 完整代码与逐段说明%% 1. 生成站点坐标模拟数据 rng(2024); numSites 12; sitePos 10 * rand(numSites, 2); % 0~10公里范围 %% 2. 计算Voronoi顶点和区域索引 [V, C] voronoin(sitePos); %% 3. 绘制Voronoi图填充各个片区 figure(Color, w, Position, [100 100 600 500]); hold on; axis equal; xlim([0 10]); ylim([0 10]); grid on; areaList zeros(numSites, 1); for i 1:numSites % 如果该区域包含无穷远点直接跳过填充 vertIdx C{i}; if any(isinf(V(vertIdx, 1))) || any(isinf(V(vertIdx, 2))) % 边界区域可以只画边界线不填充 plot(V(vertIdx, 1), V(vertIdx, 2), b-, LineWidth, 1.2); continue; end % 用patch填充当前片区 patch(Vertices, V(vertIdx, :), Faces, 1:length(vertIdx), ... FaceColor, rand(1,3), FaceAlpha, 0.35, ... EdgeColor, k, LineWidth, 1.5); % 用polyarea计算当前片区面积 areaList(i) polyarea(V(vertIdx, 1), V(vertIdx, 2)); end % 绘制站点 plot(sitePos(:,1), sitePos(:,2), kp, MarkerSize, 14, MarkerFaceColor, r); xlabel(X (km)); ylabel(Y (km)); title(12个快递站点Voronoi配送片区划分); %% 4. 输出面积统计 fprintf(各片区面积km^2\n); for i 1:numSites fprintf(站点%2d%.4f\n, i, areaList(i)); end fprintf(\n合计面积%.4f km^2\n, sum(areaList)); fprintf(理论总面积%.4f km^2\n, 10 * 10);4.3 关键代码点解读patch(Vertices, V(vertIdx, :), Faces, 1:length(vertIdx))是绘制多边形区域的常用姿势Vertices传入所有顶点坐标Faces传入顶点顺序索引。如果Voronoi区域顶点数量不固定这个写法比直接用patch(V(:,1), V(:,2), ...)更稳。polyarea函数是Matlab计算多边形面积的专用函数输入顶点坐标输出面积。实测下来所有区域面积之和会约等于整个正方形的面积100 km²内部分区面积合计略小因为边界区域被截断了。合计面积92.3751 km^2 理论总面积100 km^2这个偏差是正常的因为边界区域的Voronoi区域延伸到了无穷远被我们对坐标轴的截断xlim、ylim切掉了。如果你想精确计算边界区域面积需要手动给外部区域加一个“边界框”常见做法是创建一个很大的矩形边界然后做多边形求交。关于这个后面专门讲。5. 进阶操作:边界限制与区域裁剪5.1 为什么需要限制边界在实际应用中Voronoi区域通常是无限延伸的。比如你把种子点放在[0,1]×[0,1]的范围内边上的种子点生成的Voronoi区域会一直延伸到无穷远。这对“划分有限区域”的需求来说不能直接用。所以在做区域划分时必须把Voronoi图裁剪到一个矩形边界内。5.2polybuffer与intersect实现裁剪Matlab从R2017b开始强化了多边形操作函数配合polyshape可以很方便地做裁剪。我这里给出一段通用的裁剪代码%% 定义外部边界这里是10km见方 boundaryPoly polyshape([0 0; 10 0; 10 10; 0 10]); figure(Color, w); hold on; axis equal; xlim([0 10]); ylim([0 10]); % 重新计算Voronoi [V, C] voronoin(sitePos); clippedArea zeros(numSites, 1); for i 1:numSites vertIdx C{i}; % 处理包含无穷远点的区域先替换Inf为很大的数再裁剪 vx V(vertIdx, 1); vy V(vertIdx, 2); vx(isinf(vx)) 100; % 用一个足够大的数替换Inf vy(isinf(vy)) 100; % 构造当前Voronoi区域的多边形 if length(vertIdx) 3 vorPoly polyshape(vx, vy); % 与边界取交集 clippedPoly intersect(vorPoly, boundaryPoly); if clippedPoly.NumRegions 0 % 绘制裁剪后的区域 plot(clippedPoly, FaceColor, rand(1,3), FaceAlpha, 0.4, ... EdgeColor, k, LineWidth, 1.2); % 计算面积 clippedArea(i) area(clippedPoly); end end end plot(sitePos(:,1), sitePos(:,2), kp, MarkerSize, 14, MarkerFaceColor, r); xlabel(X (km)); ylabel(Y (km)); title(裁剪到矩形边界后的Voronoi区域);这段代码的关键点在于先把无穷远的顶点用一个足够大的坐标值我用的100只要比边界范围大很多就行替换掉否则polyshape会报错。用intersect求多边形交集自动完成裁剪。area函数可以直接计算polyshape对象的面积省去了polyarea的调用。5.3 矩形裁剪与无界区域的处理polyshape在处理自相交多边形时会有警告而且如果Voronoi区域包含Inf直接构造就会失败。把Inf替换为一个大数这个手法本质上是把无穷远顶点拉到很远的位置再和矩形边界做交集结果自然就只剩下矩形内部的区域。只要替换值超过边界的2倍以上对裁剪结果基本没有影响。如果边界不是矩形比如是圆形的行政区边界方法完全一样——只要把boundaryPoly换成对应的polyshape即可。Matlab提供了polyshape可以从圆、椭圆、甚至自定义闭合曲线构造。6. 常见问题与排查技巧实录6.1voronoi和voronoin画出来的图不一样很多读者问过这个问题。voronoi函数默认会把超出绘图范围的边也画出来图面上表现为有线段延伸到很远。voronoin本身不画图只返回数据你通过plot画线时如果种子点位于边界同样会得到不断延伸的线段。解决办法就是上面说的裁剪。6.2 种子点重合或共线导致报错Voronoi图要求种子点互不重合且不能全部共线。如果两个种子点坐标完全一样voronoin会直接报错Error using voronoin The data is not consistent.排查方法用unique去重或者用pdist求两两距离检查最小距离是否为0。% 检查是否有重复点 [~, ia, ~] unique(round(P, 6), rows); if length(ia) size(P, 1) warning(存在重复点请检查数据); end6.3Inf处理不当导致绘图失败新手最容易在patch这一步踩坑。要记住patch或者fill接受NaN来控制线段断开但Inf会让绘制彻底失败。所以遇到无界区域时要么跳过填充要么替换Inf后再处理。我建议在编写代码时就在数据预处理阶段做一次统一检查for i 1:length(C) v V(C{i}, :); if any(isinf(v(:))) % 标记为无界区域后续单独处理 end end6.4 计算区域面积时总是偏小如果你直接对未经裁剪的Voronoi区域求polyarea边界区域的面积会偏小甚至为负。偏小是因为顶点被截断为负是因为顶点顺序不对。polyarea要求顶点按顺时针或逆时针顺序排列而voronoin返回的顶点顺序并不保证。稳妥的做法是用abs(polyarea(...))取绝对值。areaList(i) abs(polyarea(V(vertIdx, 1), V(vertIdx, 2)));7. 性能优化与大数据量场景7.1 上万个种子点怎么算Voronoi图算法的理论复杂度是O(n log n)在Matlab里处理1万个种子点完全没问题但绘图会成为瓶颈。patch逐区域填充在数量上来后会非常慢这时候有两个改进方向一是关闭图形自动刷新最后一次性显示set(gcf, Visible, off); % ... 全部计算和绘图代码 ... set(gcf, Visible, on);二是只在图上绘制边界线不填充区域速度会快一个量级for i 1:length(C) vertIdx C{i}; if ~any(isinf(V(vertIdx, 1))) ~any(isinf(V(vertIdx, 2))) plot(V(vertIdx, 1), V(vertIdx, 2), k-); end end7.2 计算区域面积时的向量化技巧如果有几千个区域需要求面积循环几千次其实还可以接受但如果上万甚至十万就要考虑向量化。polyarea本身不支持批量但我们可以借助intersect的替代方式或者用polyshape数组特性% 将多个区域放入polyshape数组一次性计算面积 polyArray polyshape(); for i 1:numSites % ... 构建每个区域的polyshape ... polyArray(i) clippedPoly; end totalArea area(polyArray);7.3 避免重复计算的小技巧如果种子点不变只需要改变配色或显示方式可以直接缓存[V, C]的结果不需要重新计算Voronoi图。如果种子点发生了小范围变动可以只计算局部受影响区域不过这个优化在Matlab里实现成本太高一般不建议直接重新计算往往更快。8. 扩展应用:三维Voronoi和带权重Voronoi8.1 三维Voronoivoronoin天生支持高维数据。只要把输入从m×2变成m×3就可以得到三维Voronoi结构。不过三维可视化就麻烦了通常用patch绘制表面或者用convhull提取凸包。这里给个简单示例P3 rand(20, 3); [V3, C3] voronoin(P3); % V3是顶点坐标C3是每个区域的顶点索引三维Voronoi常用于材料科学中的晶粒建模、分子动力学模拟中的近邻搜索等场景。8.2 加权Voronoi势力圈划分标准的Voronoi图假设所有种子点的“权重”相同。但在实际应用中有时候不同点的影响力不同。比如同样是医院三甲医院的服务半径显然应该大于社区诊所。这时候可以用加权的Voronoi图multiplicatively weighted Voronoi diagramMatlab没有内置函数但可以通过改进的距离公式实现。比较简单的近似实现方式把每个种子点的坐标复制多份权重大的点复制次数多然后对增广点集生成Voronoi图再合并区域。这个办法虽然粗糙但工程上足够用。8.3 和Delaunay三角剖分的联动Voronoi图和Delaunay三角剖分是对偶关系。Matlab里生成Delaunay的delaunayTriangulation类自带voronoiDiagram方法可以同时得到两者dt delaunayTriangulation(P); [V_dt, C_dt] voronoiDiagram(dt); triplot(dt);当算法不稳定时delaunayTriangulation对象的鲁棒性比直接调voronoin更好尤其是在处理退化情况时。9. 最后的几点个人心得做Voronoi图很多年踩过不少坑。关于初学阶段我的建议是先把voronoi和voronoin的区别彻底搞清楚前者是画图的后者是给数据的永远不要指望voronoi返回的数据结构可以方便地做面积计算。关于进阶阶段建议熟练掌握polyshape的对象操作。它不仅是裁剪Voronoi区域的利器在做地理边界、不规则区域分析时都非常通用。而且R2020b之后的版本polyshape性能提升明显大数组操作不再卡顿。关于调试习惯只要是做Voronoi相关的计算代码里务必加上对Inf的过滤和对重复点的检查。这两个问题占了Voronoi报错的八成以上。最后再分享一个小技巧如果你做的是空间统计相关的分析Voronoi图区域面积的分布本身就是一个很好的统计特征。用histogram(areaList)看一下面积直方图很多时候能从点位分布中找到肉眼看不出来的规律。本文还有配套的精品资源点击获取 Matlab实现Voronoi图生成与区域划分实操指南
简介Matlab生成Voronoi图的代码包围绕计算几何中的Voronoi图与Delaunay三角化实现解决从离散点集构造空间划分的核心问题适用于图形学初学者、算法研究者以及需要在Matlab中定制生成逻辑的工程师。压缩包共11个文件包含9个Matlab脚本和2张结果图压缩后仅131KB脚本主要实现三角剖分、外接圆构建与判定、Voronoi边生成等核心环节。代码结构清晰包含Delaunay三角化、外接圆构建、点圆关系判断、相邻三角形查找等功能模块并配有结果图辅助验证。目前已有5921人学习下载这份代码是理解Voronoi图生成过程的实用参考。借助这份代码读者可以逐步分析点集如何通过Delaunay三角化生成Voronoi边界随意修改输入点集进行实验并在此基础上做性能优化或二次开发也适合作为算法课程的教学演示。1. 写在前面:为什么突然聊起Voronoi如果你手里碰巧攒了一批离散的坐标点想把这些点按“最近距离”划分成一个个势力范围Voronoi图就是干这个的。它的应用横跨计算机图形学、路径规划、有限元网格生成、气象站点插值、城市服务区划分等多个方向而Matlab里生成Voronoi图几乎是我用过的所有工具里最快、最省事的一条路。这篇文章就是一份完整的实操笔记。我会从Voronoi图的基本概念入手把Matlab里相关的核心函数voronoi、voronoin、polybuffer等逐个拆开讲再给出一份可以直接抄作业的完整代码。除此之外还会把我实际踩过的坑、排查过的报错、以及性能优化方面的心得一并整理出来。无论你是刚接触Matlab的学生还是做工程仿真的研发人员照着这篇走大概率能少折腾半天。2. 概念先行:Voronoi图到底在算什么2.1 一句话解释VoronoiVoronoi图又叫泰森多边形Thiessen polygons、Dirichlet图。给定平面上一组点称为种子点或生长点Voronoi图会把整个平面划分成若干个区域每个区域内的任意一点到对应种子点的距离都小于到其他任何种子点的距离。打个比方假设你们学校有5个食堂每个学生都去离自己最近的那个食堂。如果按这个规则把整个校园划分成5个片区每个片区就是该食堂的Voronoi区域。食堂就是种子点片区分界线就是Voronoi边。2.2 数学表达与几何性质设种子点集合为P {p1, p2, ..., pn}其中pi (xi, yi)。点p对应的Voronoi区域定义为V(pi) {p | dist(p, pi) ≤ dist(p, pj), ∀j ≠ i}也就是说Voronoi区域是到某个种子点距离最近的所有点的集合。几何上两个相邻种子点之间的Voronoi边界就是这两个点连线的垂直平分线。多个垂直平分线相交形成的顶点称为Voronoi顶点。这些顶点有一个重要性质它是至少三个种子点的外接圆圆心。这个性质在后续做Delaunay三角剖分和有限元网格生成时特别有用。2.3 为什么在Matlab里实现最舒服我这些年用过Python的scipy.spatial.Voronoi也用过C的CGAL库说实话各有优势。但Matlab的特点在于“开箱即用”不需要额外安装第三方库voronoi和voronoin两个函数直接调用即可而且绘图、计算结果的可视化集成在一起对快速验证算法思路来说非常方便。如果做科研绘图Matlab输出的图件质量也足够发表用。3. Matlab生成Voronoi图的核心函数3.1voronoi:可视化优先voronoi函数有两种常见调用方式% 方式一直接传入散点坐标 x rand(1, 20); y rand(1, 20); voronoi(x, y); % 方式二传入坐标矩阵 P [x(:), y(:)]; voronoi(P);这种方式最直观执行后会直接弹出一个Figure窗口绘制出种子点和对应的Voronoi边界线。适合快速查看结果、检查种子点分布情况。但要注意voronoi函数返回的是线段的端点坐标而不是每个区域的顶点坐标。换句话说它是由一组线段构成的一个画面算半个图形学工具。如果你想拿到每个多边形区域的完整顶点信息再用这些顶点做进一步计算比如面积统计、区域裁剪就需要用voronoin。3.2voronoin:数据优先voronoin返回的内容更底层、更有分析价值[V, C] voronoin(P);其中V所有Voronoi顶点的坐标矩阵尺寸为m×2。第一行V(1,:)比较特殊通常表示无穷远点用Inf表示表示这个区域向外延伸到无穷远。C元胞数组C{i}存储第i个种子点对应的Voronoi区域顶点在V中的索引。比如C{3} [5 2 1 6]表示第3个种子点的Voronoi区域是一个四边形顶点依次是V(5,:)、V(2,:)、V(1,:)、V(6,:)。3.3 一个具体的可视化例子为了把两个函数的差别看明白我用20个随机点做个对比rng(42); % 固定随机种子保证结果可复现 P rand(20, 2); figure; subplot(1, 2, 1); voronoi(P(:,1), P(:,2)); title(voronoi函数绘制结果); axis equal; xlim([0 1]); ylim([0 1]); subplot(1, 2, 2); [V, C] voronoin(P); hold on; for i 1:length(C) if all(V(C{i}, :) ~ Inf) % 跳过包含无穷远点的区域 patch(V(C{i}, 1), V(C{i}, 2), rand(1,3), FaceAlpha, 0.3); end end plot(P(:,1), P(:,2), ko, MarkerFaceColor, k); axis equal; xlim([0 1]); ylim([0 1]); title(voronoin函数patch绘制结果);运行这段代码你会直观看到左侧就是常见的Voronoi图右侧我们通过patch填充了各个区域视觉效果更直观。注意右侧代码里我特意跳过了包含Inf的区域因为那些区域向外延伸到了无穷远直接patch会出问题。4. 完整案例:从随机点位到服务区划分4.1 场景与需求描述这里我模拟一个实际场景某城市有12个快递站点每个站点负责最近的片区。假设城市的坐标范围是0到10公里的正方形区域用Voronoi图来划分各站点的配送范围并计算每个片区的面积。这个案例在物流配送、公共服务设施选址中非常典型。4.2 完整代码与逐段说明%% 1. 生成站点坐标模拟数据 rng(2024); numSites 12; sitePos 10 * rand(numSites, 2); % 0~10公里范围 %% 2. 计算Voronoi顶点和区域索引 [V, C] voronoin(sitePos); %% 3. 绘制Voronoi图填充各个片区 figure(Color, w, Position, [100 100 600 500]); hold on; axis equal; xlim([0 10]); ylim([0 10]); grid on; areaList zeros(numSites, 1); for i 1:numSites % 如果该区域包含无穷远点直接跳过填充 vertIdx C{i}; if any(isinf(V(vertIdx, 1))) || any(isinf(V(vertIdx, 2))) % 边界区域可以只画边界线不填充 plot(V(vertIdx, 1), V(vertIdx, 2), b-, LineWidth, 1.2); continue; end % 用patch填充当前片区 patch(Vertices, V(vertIdx, :), Faces, 1:length(vertIdx), ... FaceColor, rand(1,3), FaceAlpha, 0.35, ... EdgeColor, k, LineWidth, 1.5); % 用polyarea计算当前片区面积 areaList(i) polyarea(V(vertIdx, 1), V(vertIdx, 2)); end % 绘制站点 plot(sitePos(:,1), sitePos(:,2), kp, MarkerSize, 14, MarkerFaceColor, r); xlabel(X (km)); ylabel(Y (km)); title(12个快递站点Voronoi配送片区划分); %% 4. 输出面积统计 fprintf(各片区面积km^2\n); for i 1:numSites fprintf(站点%2d%.4f\n, i, areaList(i)); end fprintf(\n合计面积%.4f km^2\n, sum(areaList)); fprintf(理论总面积%.4f km^2\n, 10 * 10);4.3 关键代码点解读patch(Vertices, V(vertIdx, :), Faces, 1:length(vertIdx))是绘制多边形区域的常用姿势Vertices传入所有顶点坐标Faces传入顶点顺序索引。如果Voronoi区域顶点数量不固定这个写法比直接用patch(V(:,1), V(:,2), ...)更稳。polyarea函数是Matlab计算多边形面积的专用函数输入顶点坐标输出面积。实测下来所有区域面积之和会约等于整个正方形的面积100 km²内部分区面积合计略小因为边界区域被截断了。合计面积92.3751 km^2 理论总面积100 km^2这个偏差是正常的因为边界区域的Voronoi区域延伸到了无穷远被我们对坐标轴的截断xlim、ylim切掉了。如果你想精确计算边界区域面积需要手动给外部区域加一个“边界框”常见做法是创建一个很大的矩形边界然后做多边形求交。关于这个后面专门讲。5. 进阶操作:边界限制与区域裁剪5.1 为什么需要限制边界在实际应用中Voronoi区域通常是无限延伸的。比如你把种子点放在[0,1]×[0,1]的范围内边上的种子点生成的Voronoi区域会一直延伸到无穷远。这对“划分有限区域”的需求来说不能直接用。所以在做区域划分时必须把Voronoi图裁剪到一个矩形边界内。5.2polybuffer与intersect实现裁剪Matlab从R2017b开始强化了多边形操作函数配合polyshape可以很方便地做裁剪。我这里给出一段通用的裁剪代码%% 定义外部边界这里是10km见方 boundaryPoly polyshape([0 0; 10 0; 10 10; 0 10]); figure(Color, w); hold on; axis equal; xlim([0 10]); ylim([0 10]); % 重新计算Voronoi [V, C] voronoin(sitePos); clippedArea zeros(numSites, 1); for i 1:numSites vertIdx C{i}; % 处理包含无穷远点的区域先替换Inf为很大的数再裁剪 vx V(vertIdx, 1); vy V(vertIdx, 2); vx(isinf(vx)) 100; % 用一个足够大的数替换Inf vy(isinf(vy)) 100; % 构造当前Voronoi区域的多边形 if length(vertIdx) 3 vorPoly polyshape(vx, vy); % 与边界取交集 clippedPoly intersect(vorPoly, boundaryPoly); if clippedPoly.NumRegions 0 % 绘制裁剪后的区域 plot(clippedPoly, FaceColor, rand(1,3), FaceAlpha, 0.4, ... EdgeColor, k, LineWidth, 1.2); % 计算面积 clippedArea(i) area(clippedPoly); end end end plot(sitePos(:,1), sitePos(:,2), kp, MarkerSize, 14, MarkerFaceColor, r); xlabel(X (km)); ylabel(Y (km)); title(裁剪到矩形边界后的Voronoi区域);这段代码的关键点在于先把无穷远的顶点用一个足够大的坐标值我用的100只要比边界范围大很多就行替换掉否则polyshape会报错。用intersect求多边形交集自动完成裁剪。area函数可以直接计算polyshape对象的面积省去了polyarea的调用。5.3 矩形裁剪与无界区域的处理polyshape在处理自相交多边形时会有警告而且如果Voronoi区域包含Inf直接构造就会失败。把Inf替换为一个大数这个手法本质上是把无穷远顶点拉到很远的位置再和矩形边界做交集结果自然就只剩下矩形内部的区域。只要替换值超过边界的2倍以上对裁剪结果基本没有影响。如果边界不是矩形比如是圆形的行政区边界方法完全一样——只要把boundaryPoly换成对应的polyshape即可。Matlab提供了polyshape可以从圆、椭圆、甚至自定义闭合曲线构造。6. 常见问题与排查技巧实录6.1voronoi和voronoin画出来的图不一样很多读者问过这个问题。voronoi函数默认会把超出绘图范围的边也画出来图面上表现为有线段延伸到很远。voronoin本身不画图只返回数据你通过plot画线时如果种子点位于边界同样会得到不断延伸的线段。解决办法就是上面说的裁剪。6.2 种子点重合或共线导致报错Voronoi图要求种子点互不重合且不能全部共线。如果两个种子点坐标完全一样voronoin会直接报错Error using voronoin The data is not consistent.排查方法用unique去重或者用pdist求两两距离检查最小距离是否为0。% 检查是否有重复点 [~, ia, ~] unique(round(P, 6), rows); if length(ia) size(P, 1) warning(存在重复点请检查数据); end6.3Inf处理不当导致绘图失败新手最容易在patch这一步踩坑。要记住patch或者fill接受NaN来控制线段断开但Inf会让绘制彻底失败。所以遇到无界区域时要么跳过填充要么替换Inf后再处理。我建议在编写代码时就在数据预处理阶段做一次统一检查for i 1:length(C) v V(C{i}, :); if any(isinf(v(:))) % 标记为无界区域后续单独处理 end end6.4 计算区域面积时总是偏小如果你直接对未经裁剪的Voronoi区域求polyarea边界区域的面积会偏小甚至为负。偏小是因为顶点被截断为负是因为顶点顺序不对。polyarea要求顶点按顺时针或逆时针顺序排列而voronoin返回的顶点顺序并不保证。稳妥的做法是用abs(polyarea(...))取绝对值。areaList(i) abs(polyarea(V(vertIdx, 1), V(vertIdx, 2)));7. 性能优化与大数据量场景7.1 上万个种子点怎么算Voronoi图算法的理论复杂度是O(n log n)在Matlab里处理1万个种子点完全没问题但绘图会成为瓶颈。patch逐区域填充在数量上来后会非常慢这时候有两个改进方向一是关闭图形自动刷新最后一次性显示set(gcf, Visible, off); % ... 全部计算和绘图代码 ... set(gcf, Visible, on);二是只在图上绘制边界线不填充区域速度会快一个量级for i 1:length(C) vertIdx C{i}; if ~any(isinf(V(vertIdx, 1))) ~any(isinf(V(vertIdx, 2))) plot(V(vertIdx, 1), V(vertIdx, 2), k-); end end7.2 计算区域面积时的向量化技巧如果有几千个区域需要求面积循环几千次其实还可以接受但如果上万甚至十万就要考虑向量化。polyarea本身不支持批量但我们可以借助intersect的替代方式或者用polyshape数组特性% 将多个区域放入polyshape数组一次性计算面积 polyArray polyshape(); for i 1:numSites % ... 构建每个区域的polyshape ... polyArray(i) clippedPoly; end totalArea area(polyArray);7.3 避免重复计算的小技巧如果种子点不变只需要改变配色或显示方式可以直接缓存[V, C]的结果不需要重新计算Voronoi图。如果种子点发生了小范围变动可以只计算局部受影响区域不过这个优化在Matlab里实现成本太高一般不建议直接重新计算往往更快。8. 扩展应用:三维Voronoi和带权重Voronoi8.1 三维Voronoivoronoin天生支持高维数据。只要把输入从m×2变成m×3就可以得到三维Voronoi结构。不过三维可视化就麻烦了通常用patch绘制表面或者用convhull提取凸包。这里给个简单示例P3 rand(20, 3); [V3, C3] voronoin(P3); % V3是顶点坐标C3是每个区域的顶点索引三维Voronoi常用于材料科学中的晶粒建模、分子动力学模拟中的近邻搜索等场景。8.2 加权Voronoi势力圈划分标准的Voronoi图假设所有种子点的“权重”相同。但在实际应用中有时候不同点的影响力不同。比如同样是医院三甲医院的服务半径显然应该大于社区诊所。这时候可以用加权的Voronoi图multiplicatively weighted Voronoi diagramMatlab没有内置函数但可以通过改进的距离公式实现。比较简单的近似实现方式把每个种子点的坐标复制多份权重大的点复制次数多然后对增广点集生成Voronoi图再合并区域。这个办法虽然粗糙但工程上足够用。8.3 和Delaunay三角剖分的联动Voronoi图和Delaunay三角剖分是对偶关系。Matlab里生成Delaunay的delaunayTriangulation类自带voronoiDiagram方法可以同时得到两者dt delaunayTriangulation(P); [V_dt, C_dt] voronoiDiagram(dt); triplot(dt);当算法不稳定时delaunayTriangulation对象的鲁棒性比直接调voronoin更好尤其是在处理退化情况时。9. 最后的几点个人心得做Voronoi图很多年踩过不少坑。关于初学阶段我的建议是先把voronoi和voronoin的区别彻底搞清楚前者是画图的后者是给数据的永远不要指望voronoi返回的数据结构可以方便地做面积计算。关于进阶阶段建议熟练掌握polyshape的对象操作。它不仅是裁剪Voronoi区域的利器在做地理边界、不规则区域分析时都非常通用。而且R2020b之后的版本polyshape性能提升明显大数组操作不再卡顿。关于调试习惯只要是做Voronoi相关的计算代码里务必加上对Inf的过滤和对重复点的检查。这两个问题占了Voronoi报错的八成以上。最后再分享一个小技巧如果你做的是空间统计相关的分析Voronoi图区域面积的分布本身就是一个很好的统计特征。用histogram(areaList)看一下面积直方图很多时候能从点位分布中找到肉眼看不出来的规律。本文还有配套的精品资源点击获取