Matlab实现Voronoi图：从原理到代码，解决空间划分与选址问题

简介面向计算几何与科学计算的Matlab Voronoi图生成代码包适合需要从底层理解Voronoi图与Delaunay三角化的学生、研究者或算法爱好者。资源打破内置函数的黑箱围绕点集剖分、外接圆判断、邻接关系构建等核心步骤提供了一套可逐行研读的自定义函数实现。压缩包共11个文件以9个.m源文件和2张.png结果图为主整体约131KB结构紧凑便于对照学习。目前已有5921人学习下载。代码将Voronoi图生成拆分为多个独立模块三角剖分、线段构建、相邻三角形查找、外接圆计算与空圆条件判断等使用者既能按步骤理解几何关系也能直接调用或修改这些函数用于实际项目。两幅结果图直观对比了不同实现路径的输出效果方便验证算法正确性可作为Matlab图形学课程设计或计算几何入门的重要参考。 写到Voronoi图我得先讲个真事。前阵子有个做充电桩选址的朋友找我说他们现在划分每个站点的服务范围还靠人工画圈地图上几十个点一个个手动画边界干到凌晨两点还没弄完。我说你早该用Voronoi图了这玩意儿就是干这个的。他把数据发我我写了一段Matlab不到一分钟服务范围全部自动切好还带面积统计。他当场愣住问我这到底是什么神仙技术。其实真不是神仙技术就是计算几何里一个经典工具Matlab里实现它比你想象中简单但坑也确实不少。这篇我就把Voronoi图的原理、Matlab里几个函数的选型思路、完整可跑的代码、还有我踩过的那些坑一次性讲清楚。1. Voronoi图是啥能帮我解决什么问题1.1 从奶茶店选址说起想象一条街上开了几家奶茶店顾客自然选离自己最近的那家。如果按“最近原则”把整条街划分成几块服务区每家店独占一块这块区域里任何一个位置到这家的距离都小于到其他家的距离——这个划分结果就是Voronoi图。在计算几何里叫Voronoi Diagram气象学里叫泰森多边形通信工程里管它叫Voronoi分区有些做地理分析的也叫它Dirichlet tessellation。名字一大堆本质都是同一个东西。数学定义也不复杂给定平面上的一组种子点也叫站点、生成点平面被分成若干个单元每个单元内部的所有点到对应种子点的欧氏距离都小于到其他任意种子点的距离。写出来就是长这样[ V(p_i) { x \in X \mid d(x, p_i) \le d(x, p_j), \forall j \ne i } ]通俗讲就是每个种子点圈一块地谁距离它最近谁就归它管。1.2 我见过最实用的五大应用场景这玩意儿能干活的地方太多了我简单列几个我自己实际碰过或跟人聊过的场景充电桩、快递站、基站选址给出所有站点的坐标自动划出每站的服务范围一眼看出哪些站点服务面积过大、哪些区域覆盖重叠。城市公共设施规划学校、医院、消防站的服务半径分析按Voronoi单元算覆盖人口密度比拿圆规画圈科学得多。材料科学晶粒结构模拟晶粒长大过程天然就是Voronoi图的演化过程很多论文里的微观结构示意图就是Voronoi图加随机扰动生成的。生态学和动物行为研究动物领地的划分模型每个个体占一块地盘跟Voronoi图几乎一模一样。最近邻查找与路径规划机器人导航里快速找最近目标点用Voronoi图做空间划分能大幅缩短搜索时间。适合谁学做科研数据处理的研究生、搞空间分析的工程师、学计算几何的本科生甚至是有兴趣的地理学爱好者。只要你手里有一堆坐标点想知道这些点在空间上怎么“瓜分地盘”Voronoi图就是你要的工具。2. Matlab里跟Voronoi相关的函数到底该选哪个Matlab原生支持Voronoi图这一点比Python需要装scipy.spatial.Voronoi要省心。但用哪个函数取决于你想干什么——画个图看看还是拿数据做计算分析。我在实际中用下来主要涉及三个函数各有定位。2.1 voronoi(x,y)画图最快数据有限这是最直接的画图函数输入两组坐标向量直接出图x rand(1,20) * 100; y rand(1,20) * 100; voronoi(x, y);三行代码图就出来了。这个函数还有带返回值的版本返回的是Voronoi边的端点坐标[vx, vy] voronoi(x, y);这里vx和vy是2×m的矩阵每一列代表一条边的两个端点坐标。注意它返回的只是线段不是完整的多边形区域。这意味着它适合快速可视化但不适合做面积计算、拓扑分析这类后续操作。2.2 voronoin(P)数据最全核心主力这是我最常用的函数。输入一个n×2的坐标矩阵P返回顶点矩阵V和单元索引元胞数组C[V, C] voronoin(P);V是k×2的顶点坐标矩阵每一行是一个Voronoi顶点的坐标。C是一个n×1的元胞数组C{i}存的是第i个种子点对应的Voronoi单元按顺序经过的顶点索引——注意是索引不是坐标要查坐标得用V(C{i}, 1)和V(C{i}, 2)。这里有一个很关键、新手必坑的细节如果某个种子点位于凸包边界上它的Voronoi单元会延伸到无穷远此时C{i}的第一个元素是1而V(1,:)在Matlab中被定义为无穷远点[Inf Inf]。后面算面积的时候只要发现C{i}的第一个元素是1基本可以判定这个单元是开放的不能直接算面积。2.3 polybnd_voronoi边界裁剪工程必备原生voronoin算出来的无限延伸区域在实际工程里很头疼——你研究的是一个100×100的研究区结果边界上的点到图外去了这不合理。所以我在做正经项目时90%的情况会用File Exchange上的第三方函数polybnd_voronoi。这个函数的核心能力是在给定边界多边形比如矩形研究区边界内裁剪Voronoi图让所有单元都被限制在边界之内没有无限延伸问题边界处的单元也能完整闭合。函数的调法大致是[vx, vy] polybnd_voronoi(x, y, bdry_xy);其中bdry_xy是边界多边形的顶点坐标按顺序排列成闭合多边形。返回的vx和vy是裁剪后的Voronoi边线段端点用它画图就是完整闭合的研究区划分。还有个进阶版本调用方式能直接返回裁剪后每个单元的多边形顶点方便我直接算面积、做统计。这一点在做定量分析时太重要了后面实操环节我详细讲。3. 从零开始的完整实操Voronoi图生成与数据提取3.1 第一步准备种子点数据做任何空间分析第一步都是先把种子点弄干净。我习惯用固定的随机种子保证每次跑出来的结果一样方便调参和复现% 生成30个随机种子点分布在100×100的研究区内 rng(42); % 固定随机种子保证结果可复现 num 30; x rand(num,1) * 100; y rand(num,1) * 100; P [x, y]; % 组装成voronoin需要的输入格式如果你的种子点来自真实业务数据比如充电桩坐标、基站经纬度投影坐标记得先检查有没有重复点、有没有超出研究区的异常点。这两类脏数据会在后面引发非常诡异的问题后面第4节详述。3.2 第二步画图与数据提取双轨并行先上最核心的代码直接能跑% 基础可视化 figure(Color,w); voronoi(x,y); axis equal; grid on; xlim([0 100]); ylim([0 100]); title(基础Voronoi图);注意axis equal这行不能省。如果不加坐标轴比例不一致图形会被拉伸Voronoi图看起来就会变形垂直的边可能看起来不是垂直的影响你对结果直观判断。但如果要做数据分析光画图不够必须拿数据。用voronoin提取拓扑信息% 提取Voronoi图顶点和单元拓扑 [V, C] voronoin(P); % 查看前三个单元的顶点索引 disp(C(1:3));C{i}里的索引顺序是逆时针环绕的这直接决定了后面用patch填色、用polyarea算面积时不会出现自相交的混乱图形。Matlab在这点上的设计是很贴心的但前提是你知道它有序。3.3 第三步在边界上裁剪得到真正闭合的单元原生voronoin画出来的整个Voronoi图四周是发散出去的开放区域。研究区通常限定在一个矩形或任意多边形范围内所以需要裁剪。我用得最多的方案是polybnd_voronoi给一个矩形边界就能把所有开放单元收拢% 定义矩形研究区边界逆时针闭合 bdry [0 0; 100 0; 100 100; 0 100; 0 0]; % 调用polybnd_voronoi计算裁剪后的Voronoi边 [vx, vy] polybnd_voronoi(x, y, bdry); % 绘制裁剪后的Voronoi图 figure(Color,w); plot(vx, vy, b-, LineWidth, 1.2); hold on; plot(x, y, ro, MarkerFaceColor, r, MarkerSize, 6); axis equal; xlim([0 100]); ylim([0 100]); grid on; title(边界裁剪后的Voronoi图);这一段完成后每个种子点都被一个闭合的多边形围住而且所有多边形都落在研究区内部。这个形态才是做土地利用分析、服务范围统计时想要的结果。3.4 第四步计算每个单元的面积与服务范围统计有了裁剪后的Voronoi图面积计算就有意义了。我刚说polybnd_voronoi的进阶版本能返回裁剪后的单元多边形这里用Matlab自带的polyarea就能逐块算面积% 如果polybnd_voronoi返回裁剪后的顶点和单元索引 % 假设返回格式为 [vx, vy, bvx, bvy]其中bvx/bvy是裁剪后单元多边形坐标 % 下面用原生顶点数据做面积统计示例注意判断开放单元 areas zeros(num, 1); for i 1:num if any(C{i} 1) % 开放单元面积记为NaN表示无法直接计算 areas(i) NaN; else % 闭合单元用polyarea计算面积 areas(i) polyarea(V(C{i}, 1), V(C{i}, 2)); end end % 查看面积统计 fprintf(平均单元面积: %.2f\n, mean(areas, omitnan)); fprintf(最大单元面积: %.2f\n, max(areas)); fprintf(最小单元面积: %.2f\n, min(areas));这段代码跑完每个种子点对应的服务面积一目了然。实际项目中我还会把面积数据关联回原始业务表比如结合每个充电桩的功率、周边人口密度算加权服务覆盖率。这是Voronoi图最有价值的应用场景之一——把几何划分和业务指标打通。4. 实测中踩过的坑Voronoi图常见问题与排查4.1 无限顶点问题边界种子点的开放单元这是新手最常踩的坑。当你用voronoin时只要种子点位于整组点集的凸包边界上它的Voronoi单元就会向外无限延伸反映在C{i}里就是第一个元素是1。我之前有个同事拿着一组24个基站的数据算覆盖面积没检查开放单元就直接调polyarea跑出来一堆荒谬的负数面积他还以为算法有问题排查了一下午。解决办法就是在循环里加开放单元判断——if any(C{i} 1)。如果是开放单元要么跳过要么用polybnd_voronoi裁剪后再算。我个人的建议是做可视化的时候用voronoin没问题做面积统计时必须裁剪或者至少把开放单元单独标注出来避免统计数据被污染。4.2 顶点顺序错乱导致的多边形异常正常情况下C{i}里的顶点索引是按逆时针顺序排列的直接拿去patch或polyarea都行。但如果你手动处理过顶点顺序或者从某些第三方库导入数据顺序一旦变乱画出来的多边形就是个打结的线团面积计算也会出错。排查方法很简单随便取一个单元把C{i}和V(C{i},:)打出来看看手工连线画一画。如果发现顺序不对用sort函数结合角度排序可以修复但更保险的做法是保持原生数据不动所有下游操作都在原始C和V基础上做。4.3 种子点重合或共线种子点完全重合会让Voronoi图退化Matlab可能直接报错或者生成一个面积为0的畸变单元。共线更隐蔽——如果好几个点精确落在同一条直线上Voronoi图的边会出现共线退化图形呈锯齿状运行时不报错但结果不美观。4.4 voronoi图形显示不完整坐标轴比例陷阱别看这个问题小我见过不止一个人卡在这里。如果画voronoi图时忘了axis equal图形会被拉伸Voronoi图的边看起来就不垂直、不平行甚至会让人觉得算错了。实际上数据没错是坐标轴比例失真。加了axis equal后图形比例恢复正常边界也严丝合缝。这是一个典型的“图形渲染问题”被误认为是“算法问题”的场景。4.5 第三方函数安装与版本兼容polybnd_voronoi是在File Exchange上的函数下载后要把m文件放到当前工作目录或MATLAB路径下。有朋友遇到过下载后运行时提示“未定义函数”原因往往是路径没加进MATLAB。用addpath命令手动加一下即可addpath(你存放函数的路径);同时要注意不同Matlab版本对函数语法的兼容性老代码在新版本尤其是R2023b之后偶发报错基本都能在MathWorks官方文档的Release Notes里查到原因。5. 进阶玩法从“画个图”到“空间统计出图”5.1 按面积/类别给Voronoi单元上色拿到每个单元的面积后下一步就是把结果可视化。我常用的方式是用patch函数上色颜色映射服务面积大小这样一眼就能看出哪些区域覆盖不足figure(Color,w); hold on; for i 1:num if ~any(C{i} 1) % 用面积作为颜色映射的输入 patch(Faces, C{i}, Vertices, V, ... FaceColor, [0.2 0.6 0.9], ... EdgeColor, k, LineWidth, 0.8); end end plot(x, y, ro, MarkerFaceColor, r, MarkerSize, 6); axis equal; grid on;如果想按面积渐变着色可以在循环里给每个patch的FaceColor设置一个和面积成比例的值再搭配colormap来映射。5.2 三维Voronoi能玩吗Matlab的voronoin其实支持三维输入返回的是三维顶点和三维单元。但说实话三维Voronoi图的可视化在Matlab里做得不算友好出图效果一般。我个人的建议是如果你只是想做三维空间划分可视化用Python的scipy.spatial.Voronoi加matplotlib或者ParaView这类专业软件效果会好很多。Matlab的生是在二维数据和二维可视化上不必硬上三维。5.3 权重Voronoi图当每个点的“势力”不一样时标准的Voronoi图假设每个种子点的权重相同服务范围只取决于距离。但现实中并非如此——一个大功率充电桩的服务范围肯定比小功率的大一家旗舰店的辐射半径肯定比普通店大。这时就要用加权Voronoi图。Matlab原生没有直接实现加权Voronoi图但可以自己写。最常用的近似方法是用一个迭代求解器将每个种子点赋予权重通过优化划分边界逼近加权Voronoi图另一个思路在图像域做逐像素分类对所有像素点比较“距离/权重”取最小值归类。实测下来逐像素方法在几百个种子点级别完全可跑就是把每个像素当作一个数据点比较距离代码几十行效果很直观适合科研成果展示。5.4 结合地图数据的落地实践最后再说一个真实项目中常见的做法如果你手上是经纬度坐标想叠加到地图底图上做服务区展示先把经纬度用projfwd或geodetic2ecef等工具投影到平面坐标中国区域常用UTM投影或高斯-克吕格投影再用投影坐标做Voronoi图。得出的多边形再反投影回经纬度就能叠加到地图引擎上展示。这一步别省——直接用经纬度做Voronoi在数学上是不成立的因为球面上的测地距离不等于平面欧氏距离算出来的Voronoi边在地图上会偏差得离谱。从我个人的经验看Voronoi图的价值从来不在“能画出来”而在它帮你把空间划分问题转化为可量化的几何计算。不管是算服务范围、做设施布局优化还是分析空间分布规律只要你手里有一组坐标点Voronoi图就是那把最趁手的工具。Matlab做这件事的性价比极高原生的两个函数加上一个裁剪函数就能覆盖从快速出图到精细分析的全部场景。如果你在实际项目里跑出了奇怪的结果优先按上面第4节列的几个方向排查开放单元、顶点顺序、坐标轴比例、重合点八成问题都出在这几处。本文还有配套的精品资源点击获取