简介面向计算几何学习与二次开发的Matlab Voronoi图生成代码包基于Delaunay三角化实现空间剖分完整覆盖从点集输入、三角网格构建到Voronoi边生成的全流程。Voronoi图将平面分割为距离最近生成点的区域是计算几何中的重要结构。整个压缩包共11个文件包含9个.m源码文件和2张结果验证图大小仅131KB代码结构紧凑、模块划分清晰。目前已有5921人学习适用于计算几何课程设计、算法复现或用来替换Matlab内置voronoi、delaunay函数以满足定制需求。源码按功能拆分为Delaunay三角化、外接圆判断、相邻三角形查找、线条构建等独立模块便于逐段理解与二次修改附带的结果图可直观对照不同实现方式的输出快速验证算法正确性。该算法在GIS最近设施搜索、生物医学图像分析、材料微观结构分析等领域也有典型应用场景对计算几何初学者和课题研究者都具有较好的参考价值。 给一个配送站覆盖划分的需求站点坐标倒是齐全用Matlab画Voronoi图也很快但对方真正想要的是能导出到规划报告里的面状图要求每个站点的“势力范围”必须落在指定行政区边界内。这才发现光会用voronoi函数远远不够。这篇文章就围绕Matlab生成Voronoi图代码的实际落地过程来写从基础函数、数据结构、裁剪方式到性能优化把我踩过的坑和验证过的写法都梳理一遍适合手里有站点坐标、想做区域划分或者看板可视化的读者直接参考。1. Voronoi图到底是什么从最近邻划分到泰森多边形1.1 通俗理解一个“势力范围”Voronoi图的概念其实不复杂。给平面上撒一堆站点把整个平面按照“离哪个站点最近”切分成若干多边形每个多边形内部任意一点到本站点都比到其他站点更近。这个多边形就是这一站的“势力范围”。气象学里叫泰森多边形地理信息系统里也叫Thiessen多边形底层都是同一套东西。我最早接触它是在做配送站覆盖分析的时候每个站负责周围一片订单范围不能重叠也不能留空。Voronoi图天然满足这两个要求所以算法选型几乎没有犹豫。但真正动手后才发现Matlab默认画出来的东西和“能用的面状图”之间还隔着一层数据处理。1.2 为什么值得单独研究一段生成代码Matlab里生成Voronoi图理论上两个命令就能搞定看似门槛很低。但默认接口返回的是线框图形拿不到可直接计算的多边形面积也拿不到边界裁剪后的有效面片。加上Voronoi图在凸包外侧会出现无限延伸的边界直接按默认方式绘制经常会得到超出预期范围的线条。所以这里说的“Matlab生成Voronoi图代码”不只是一个调用函数的问题而是从原始点坐标到最终可交付面状结果的一整套流程。下面我从最基础的函数调用讲起逐步深入到工程中真正需要的细节。2. 先把基础函数跑通voronoi与voronoin的最小实现2.1 最直接的voronoi用法如果只是临时看一眼点集的分布用自带函数就够了。rng(42); P rand(20, 2); figure; voronoi(P(:, 1), P(:, 2)); axis equal;这段代码会画出20个随机点的Voronoi线框图。注意voronoi函数接收的参数分别是横坐标数组和纵坐标数组也可以直接传一个N×2的矩阵进去。它会自动计算并绘制出线框方便是方便但它返回的图形句柄在版本之间还不完全一致有的版本返回Line对象有的版本返回句柄数组后期想精确控制样式会很别扭。所以我的建议是voronoi只适合快速预览真正需要拿数据做分析时换成voronoin。2.2 voronoin返回的V和C到底是什么voronoin是获取Voronoi几何数据的标准入口用法很简单[V, C] voronoin(P);P是N×2的站点坐标矩阵。返回的V是顶点坐标矩阵所有Voronoi多边形的顶点都汇总在V里但并不是每个站点对应一组独立顶点相邻多边形会共享公共顶点所以V的行数往往比站点数多得多。C是一个元胞数组长度和站点数一致。C{i}保存的是第i个站点对应多边形顶点的索引号。举例来说如果C{3}是[5 8 12 6]那第3个站点所在多边形的四个顶点就是V的第5、第8、第12、第6行逆序或顺序围成一圈。刚开始用这个结构时我总喜欢直接拿C去索引V但经常忘了C里存的不是坐标是行号。这个弯一旦绕过来后面所有绘制和计算就顺了。2.3 索引1就是无穷远点大部分坑都从这来voronoin返回的V第一行通常是一个特殊顶点[Inf, Inf]。为什么会有这样一行因为凸包外侧站点的Voronoi区域会向平面无穷远处延伸计算机无法用有限坐标表示于是Matlab约定用索引1指向无穷远点。这意味着当你遍历C{i}发现里面有1时就知道这个多边形不是封闭的它至少有一个方向通向无穷远。for i 1:length(C) if ismember(1, C{i}) fprintf(站点%d的多边形包含无穷远点\n, i); end end这个细节是很多初学者容易忽略的。直接拿V(C{i}, :)去画图一旦索引到1就会画出跑到图框外面去的线甚至让坐标轴范围变得不可控。我接下来说的绘制和裁剪很多操作都是围绕怎么处理这个特殊点展开的。3. 从线框图到可交付面状图着色、裁剪与无限边处理3.1 按自己的方式绘制多边形面片要画出面状图就不能再用voronoi的默认线框了得用patch按多边形绘制。一个最基础的面状图代码像这样figure; hold on; for i 1:length(C) if all(C{i} ~ 1) patch(V(C{i}, 1), V(C{i}, 2), w, ... EdgeColor, [0.2 0.4 0.7], LineWidth, 1.2); end end plot(P(:, 1), P(:, 2), r., MarkerSize, 12); axis equal;这里用all(C{i} ~ 1)过滤掉包含无穷远点的多边形就是为了避免画出无限边界。注意patch会自动把首尾顶点连起来形成闭合区域不需要手动重复终点。如果使用plot然后手动闭合很容易漏掉最后一段边导致图形有一个小缺口。从实际效果来看用patch绘制的另一个好处是可以同时设置填充色、边线颜色和透明度做展示图会直观很多。3.2 给每个多边形上色并映射数值实际项目中很少有人只需要看白色多边形更多时候需要把一个指标映射到颜色上比如每个站点的订单量、覆盖人口、温度观测值等。这时可以利用patch的颜色映射能力把数值挂到每个多边形上。values rand(length(C), 1); % 示例假设每个站点有个指标 figure; hold on; for i 1:length(C) if all(C{i} ~ 1) patch(V(C{i}, 1), V(C{i}, 2), values(i), ... EdgeColor, none); end end colormap(parula); colorbar; axis equal;这个写法的关键是把values(i)作为patch的第三个输入参数。Matlab会把它解释成每个面片的CData然后结合当前colormap自动映射出颜色。这样画出来的图每个区域颜色深浅直接对应数值高低效果比在循环里手动改颜色要自然得多。不过有一点要提醒如果站点数量很大循环里逐面片绘制依然会慢。后面我会专门说大数据量场景下的优化写法。3.3 把超出视野范围的线处理干净包含无穷远点的多边形虽然在循环里被跳过了但这样会导致凸包外侧的区域直接空出来图上就看到一片空白。有两种处理方式第一种直接设置坐标轴范围让空白区域不出现在视野内这种做法适合对图形完整性要求不高的场景。xlim([0 1]); ylim([0 1]);第二种把这些外侧多边形补成一个有限的闭合图形先取有限顶点再按边界盒的四角补上对应的角点使区域封闭在绘图范围内。这个思路做起来稍复杂但效果更完整。如果只是常规汇报展示我会优先用第一种因为实现成本低视觉上也不会觉得缺东西。4. 复合边界下的裁剪方案区域约束、polyshape与退化处理4.1 为什么默认生成的结果不是你想要的范围前面生成的都是覆盖整个平面的Voronoi图但实际业务中你可能只需要某个矩形区域、某个行政区范围或某个自然地块内部的划分结果。这时候仅仅靠视觉范围限制不够必须做真正的几何裁剪把每个Voronoi多边形与约束区域求交。以矩形区域为例假设约束范围是[0,1]×[0,1]最直接的做法是用polyshape对象做相交运算。Matlab从R2017b开始完整支持polyshape后续版本做多边形布尔运算非常方便比老版本的polybool接口好用得多也稳定得多。4.2 用polyshape实现区域裁剪下面这段代码演示了如何对有限多边形做裁剪。约束区域用polyshape定义然后对每个不包含无穷远点的Voronoi多边形执行intersect。rng(42); P rand(20, 2); [V, C] voronoin(P); bx [0 1 1 0]; by [0 0 1 1]; boundary polyshape(bx, by); figure; hold on; for i 1:length(C) if all(C{i} ~ 1) poly polyshape(V(C{i}, 1), V(C{i}, 2)); polyClipped intersect(poly, boundary); if polyClipped.NumRegions 0 plot(polyClipped, FaceColor, [0.4 0.6 0.8], ... EdgeColor, [0.1 0.2 0.4], LineWidth, 1.0); end end end plot(P(:, 1), P(:, 2), r., MarkerSize, 12); axis equal;intersect返回的是两个多边形的公共部分可能是一个polyshape对象也可能多个区域。用polyClipped.NumRegions判断区域数量是为了避免空结果直接画出来报错。polyshape处理带洞多边形或退化多边形时会自己修正很多边角问题但有些情况仍然需要手动过滤。4.3 退化多边形的过滤一个不关注就会翻车的细节裁剪最容易翻车的点不在代码逻辑而在几何退化。当某个Voronoi顶点刚好落在约束边界上或者两个多边形相交成一条线时intersect返回的结果面积可能极小甚至趋近于0但形式上仍然是一个合法polyshape。如果不做过滤后续计算面积、统计覆盖度时会出现一堆接近0的伪区域结果看着就很离谱。我的处理习惯是统一加一层“面积阈值过滤”minArea 1e-6; if polyClipped.NumRegions 0 area(polyClipped) minArea plot(polyClipped); end这个阈值要根据坐标尺度调整。如果坐标范围是几十千米1e-6可能太小建议用约束区域面积的百万分之一作为下限。这样既能滤掉退化区域又不会误删正常小区域。还有一个更隐蔽的坑当某个站点紧贴约束边界时它本身的Voronoi多边形在边界外部的那一部分在裁剪后可能变成一个极窄的多边形面积很小但在业务上看它仍然存在。所以过滤前先想清楚需求是“完全不要微小区域”还是“保留但标记出来”这会影响过滤策略。5. 上万点集时的性能实测与绘制策略5.1 从一千点到五万点的耗时变化很多Voronoi图的入门示例都只画几十个点但实际场景里可能遇到几千、几万个站点。我拿随机分布的点做了简单测试用tic/toc统计voronoin的计算耗时结果大致如下点数量计算耗时约1000.002秒1,0000.015秒5,0000.08秒10,0000.18秒50,0001.20秒Matlab底层用的是Qhull库二维Voronoi计算本身很快几万点也只是一两秒的事所以计算阶段通常不是瓶颈。真正卡住的地方在绘制阶段尤其是我前面写的那种逐面片patch的写法画几万个多边形可能要十几秒甚至几十秒。如果需要对很多个点集批量生成我建议先把核心计算写成函数避免在循环里反复重复绘图代码。计算和绘制分开定位问题也会更方便。5.2 绘制阶段比计算阶段更容易卡死逐个调用patch绘制大量多边形时每个patch都是一个图形对象几万个图形对象会把OpenGL渲染器直接拖垮。我踩过这个坑一万个点左右代码跑了快二十秒。解决办法是合并绘制。把有效多边形的顶点预先填充到NaN分隔的大数组里一次性调用patch或者fill。还有一种更省事的思路利用polyshape的数组特性把所有裁剪后的多边形对象放进一个数组然后统一plot。shapes polyshape.empty(); for i 1:length(C) if all(C{i} ~ 1) poly polyshape(V(C{i}, 1), V(C{i}, 2)); polyClipped intersect(poly, boundary); if polyClipped.NumRegions 0 area(polyClipped) 1e-6 shapes(end1) polyClipped; end end end plot(shapes);polyshape数组统一plot的效率远高于逐面片画图而且颜色、边界样式可以统一设置。如果仍然觉得慢就该考虑在画图时降低显示精度比如用reduce函数简化多边形顶点数视觉影响不大但渲染速度能再上一个台阶。6. 数值陷阱与工程化封装建议6.1 重复站点带来的Qhull退化Voronoi计算建立在点集几何关系上如果两个站点坐标完全重复或者距离非常非常近Qhull会判定输入退化轻则报警告重则直接报错。我遇到过坐标数据里混入了重复采集点导致后面所有结果全部错乱的情况。解决方式很简单计算前先去重[P, ~, ic] unique(P, rows);别忘了去重后站点数量变了如果原本还关联了其他业务数据要同步处理索引对应关系。对于距离极近但未完全重复的点可以按最小间距阈值做一次合并比如先对P做四舍五入到小数点后6位然后再unique能减少很多不必要的数值扰动。6.2 坐标尺度过大或过小的影响坐标数值非常大或非常小时双精度浮点数也会在几何计算中损失精度。比如经纬度坐标直接参与计算虽然能出结果但裁剪判断时容易出现边界不稳。更稳妥的做法是先对点集做中心化和平移缩放统一到单位坐标系下完成所有Voronoi计算和裁剪最后再把结果缩放回原坐标系。P0 P - mean(P); scale max(abs(P0), [], all); P_norm P0 / scale; % 计算完成后再乘回scale并加回mean这一步看起来多此一举但能省掉很多说不清的边界误差。尤其是做高精度区域计算时坐标归一化是值得养成的好习惯。6.3 把代码封装成可复用的文件最后说一点工程化的经验。Voronoi图生成的代码如果每次用都重新复制一遍很容易在某个版本里改坏一个参数。我现在的做法是封装成一个函数输入站点坐标和约束边界输出裁剪后的polyshape数组function shapes voronoiClip(P, boundary, minArea) [V, C] voronoin(P); shapes polyshape.empty(); for i 1:length(C) if all(C{i} ~ 1) poly polyshape(V(C{i}, 1), V(C{i}, 2)); poly intersect(poly, boundary); if poly.NumRegions 0 area(poly) minArea shapes(end1) poly; end end end end这样每次调用只需要一行代码核心逻辑也不会被业务代码干扰。如果在实际项目里发现某些站点特别靠近约束边界导致裁剪后面积极小我的习惯是把阈值逻辑做成可选参数不传就不过滤保留完整几何结果方便排查问题。这个小习惯帮我在好几次数据异常时快速定位了问题来源。本文还有配套的精品资源点击获取 
长河君 · 2026/9/2 16:31:59 · 阅读 1.2k Matlab生成Voronoi图并裁剪到行政区边界的完整实践
简介面向计算几何学习与二次开发的Matlab Voronoi图生成代码包基于Delaunay三角化实现空间剖分完整覆盖从点集输入、三角网格构建到Voronoi边生成的全流程。Voronoi图将平面分割为距离最近生成点的区域是计算几何中的重要结构。整个压缩包共11个文件包含9个.m源码文件和2张结果验证图大小仅131KB代码结构紧凑、模块划分清晰。目前已有5921人学习适用于计算几何课程设计、算法复现或用来替换Matlab内置voronoi、delaunay函数以满足定制需求。源码按功能拆分为Delaunay三角化、外接圆判断、相邻三角形查找、线条构建等独立模块便于逐段理解与二次修改附带的结果图可直观对照不同实现方式的输出快速验证算法正确性。该算法在GIS最近设施搜索、生物医学图像分析、材料微观结构分析等领域也有典型应用场景对计算几何初学者和课题研究者都具有较好的参考价值。 给一个配送站覆盖划分的需求站点坐标倒是齐全用Matlab画Voronoi图也很快但对方真正想要的是能导出到规划报告里的面状图要求每个站点的“势力范围”必须落在指定行政区边界内。这才发现光会用voronoi函数远远不够。这篇文章就围绕Matlab生成Voronoi图代码的实际落地过程来写从基础函数、数据结构、裁剪方式到性能优化把我踩过的坑和验证过的写法都梳理一遍适合手里有站点坐标、想做区域划分或者看板可视化的读者直接参考。1. Voronoi图到底是什么从最近邻划分到泰森多边形1.1 通俗理解一个“势力范围”Voronoi图的概念其实不复杂。给平面上撒一堆站点把整个平面按照“离哪个站点最近”切分成若干多边形每个多边形内部任意一点到本站点都比到其他站点更近。这个多边形就是这一站的“势力范围”。气象学里叫泰森多边形地理信息系统里也叫Thiessen多边形底层都是同一套东西。我最早接触它是在做配送站覆盖分析的时候每个站负责周围一片订单范围不能重叠也不能留空。Voronoi图天然满足这两个要求所以算法选型几乎没有犹豫。但真正动手后才发现Matlab默认画出来的东西和“能用的面状图”之间还隔着一层数据处理。1.2 为什么值得单独研究一段生成代码Matlab里生成Voronoi图理论上两个命令就能搞定看似门槛很低。但默认接口返回的是线框图形拿不到可直接计算的多边形面积也拿不到边界裁剪后的有效面片。加上Voronoi图在凸包外侧会出现无限延伸的边界直接按默认方式绘制经常会得到超出预期范围的线条。所以这里说的“Matlab生成Voronoi图代码”不只是一个调用函数的问题而是从原始点坐标到最终可交付面状结果的一整套流程。下面我从最基础的函数调用讲起逐步深入到工程中真正需要的细节。2. 先把基础函数跑通voronoi与voronoin的最小实现2.1 最直接的voronoi用法如果只是临时看一眼点集的分布用自带函数就够了。rng(42); P rand(20, 2); figure; voronoi(P(:, 1), P(:, 2)); axis equal;这段代码会画出20个随机点的Voronoi线框图。注意voronoi函数接收的参数分别是横坐标数组和纵坐标数组也可以直接传一个N×2的矩阵进去。它会自动计算并绘制出线框方便是方便但它返回的图形句柄在版本之间还不完全一致有的版本返回Line对象有的版本返回句柄数组后期想精确控制样式会很别扭。所以我的建议是voronoi只适合快速预览真正需要拿数据做分析时换成voronoin。2.2 voronoin返回的V和C到底是什么voronoin是获取Voronoi几何数据的标准入口用法很简单[V, C] voronoin(P);P是N×2的站点坐标矩阵。返回的V是顶点坐标矩阵所有Voronoi多边形的顶点都汇总在V里但并不是每个站点对应一组独立顶点相邻多边形会共享公共顶点所以V的行数往往比站点数多得多。C是一个元胞数组长度和站点数一致。C{i}保存的是第i个站点对应多边形顶点的索引号。举例来说如果C{3}是[5 8 12 6]那第3个站点所在多边形的四个顶点就是V的第5、第8、第12、第6行逆序或顺序围成一圈。刚开始用这个结构时我总喜欢直接拿C去索引V但经常忘了C里存的不是坐标是行号。这个弯一旦绕过来后面所有绘制和计算就顺了。2.3 索引1就是无穷远点大部分坑都从这来voronoin返回的V第一行通常是一个特殊顶点[Inf, Inf]。为什么会有这样一行因为凸包外侧站点的Voronoi区域会向平面无穷远处延伸计算机无法用有限坐标表示于是Matlab约定用索引1指向无穷远点。这意味着当你遍历C{i}发现里面有1时就知道这个多边形不是封闭的它至少有一个方向通向无穷远。for i 1:length(C) if ismember(1, C{i}) fprintf(站点%d的多边形包含无穷远点\n, i); end end这个细节是很多初学者容易忽略的。直接拿V(C{i}, :)去画图一旦索引到1就会画出跑到图框外面去的线甚至让坐标轴范围变得不可控。我接下来说的绘制和裁剪很多操作都是围绕怎么处理这个特殊点展开的。3. 从线框图到可交付面状图着色、裁剪与无限边处理3.1 按自己的方式绘制多边形面片要画出面状图就不能再用voronoi的默认线框了得用patch按多边形绘制。一个最基础的面状图代码像这样figure; hold on; for i 1:length(C) if all(C{i} ~ 1) patch(V(C{i}, 1), V(C{i}, 2), w, ... EdgeColor, [0.2 0.4 0.7], LineWidth, 1.2); end end plot(P(:, 1), P(:, 2), r., MarkerSize, 12); axis equal;这里用all(C{i} ~ 1)过滤掉包含无穷远点的多边形就是为了避免画出无限边界。注意patch会自动把首尾顶点连起来形成闭合区域不需要手动重复终点。如果使用plot然后手动闭合很容易漏掉最后一段边导致图形有一个小缺口。从实际效果来看用patch绘制的另一个好处是可以同时设置填充色、边线颜色和透明度做展示图会直观很多。3.2 给每个多边形上色并映射数值实际项目中很少有人只需要看白色多边形更多时候需要把一个指标映射到颜色上比如每个站点的订单量、覆盖人口、温度观测值等。这时可以利用patch的颜色映射能力把数值挂到每个多边形上。values rand(length(C), 1); % 示例假设每个站点有个指标 figure; hold on; for i 1:length(C) if all(C{i} ~ 1) patch(V(C{i}, 1), V(C{i}, 2), values(i), ... EdgeColor, none); end end colormap(parula); colorbar; axis equal;这个写法的关键是把values(i)作为patch的第三个输入参数。Matlab会把它解释成每个面片的CData然后结合当前colormap自动映射出颜色。这样画出来的图每个区域颜色深浅直接对应数值高低效果比在循环里手动改颜色要自然得多。不过有一点要提醒如果站点数量很大循环里逐面片绘制依然会慢。后面我会专门说大数据量场景下的优化写法。3.3 把超出视野范围的线处理干净包含无穷远点的多边形虽然在循环里被跳过了但这样会导致凸包外侧的区域直接空出来图上就看到一片空白。有两种处理方式第一种直接设置坐标轴范围让空白区域不出现在视野内这种做法适合对图形完整性要求不高的场景。xlim([0 1]); ylim([0 1]);第二种把这些外侧多边形补成一个有限的闭合图形先取有限顶点再按边界盒的四角补上对应的角点使区域封闭在绘图范围内。这个思路做起来稍复杂但效果更完整。如果只是常规汇报展示我会优先用第一种因为实现成本低视觉上也不会觉得缺东西。4. 复合边界下的裁剪方案区域约束、polyshape与退化处理4.1 为什么默认生成的结果不是你想要的范围前面生成的都是覆盖整个平面的Voronoi图但实际业务中你可能只需要某个矩形区域、某个行政区范围或某个自然地块内部的划分结果。这时候仅仅靠视觉范围限制不够必须做真正的几何裁剪把每个Voronoi多边形与约束区域求交。以矩形区域为例假设约束范围是[0,1]×[0,1]最直接的做法是用polyshape对象做相交运算。Matlab从R2017b开始完整支持polyshape后续版本做多边形布尔运算非常方便比老版本的polybool接口好用得多也稳定得多。4.2 用polyshape实现区域裁剪下面这段代码演示了如何对有限多边形做裁剪。约束区域用polyshape定义然后对每个不包含无穷远点的Voronoi多边形执行intersect。rng(42); P rand(20, 2); [V, C] voronoin(P); bx [0 1 1 0]; by [0 0 1 1]; boundary polyshape(bx, by); figure; hold on; for i 1:length(C) if all(C{i} ~ 1) poly polyshape(V(C{i}, 1), V(C{i}, 2)); polyClipped intersect(poly, boundary); if polyClipped.NumRegions 0 plot(polyClipped, FaceColor, [0.4 0.6 0.8], ... EdgeColor, [0.1 0.2 0.4], LineWidth, 1.0); end end end plot(P(:, 1), P(:, 2), r., MarkerSize, 12); axis equal;intersect返回的是两个多边形的公共部分可能是一个polyshape对象也可能多个区域。用polyClipped.NumRegions判断区域数量是为了避免空结果直接画出来报错。polyshape处理带洞多边形或退化多边形时会自己修正很多边角问题但有些情况仍然需要手动过滤。4.3 退化多边形的过滤一个不关注就会翻车的细节裁剪最容易翻车的点不在代码逻辑而在几何退化。当某个Voronoi顶点刚好落在约束边界上或者两个多边形相交成一条线时intersect返回的结果面积可能极小甚至趋近于0但形式上仍然是一个合法polyshape。如果不做过滤后续计算面积、统计覆盖度时会出现一堆接近0的伪区域结果看着就很离谱。我的处理习惯是统一加一层“面积阈值过滤”minArea 1e-6; if polyClipped.NumRegions 0 area(polyClipped) minArea plot(polyClipped); end这个阈值要根据坐标尺度调整。如果坐标范围是几十千米1e-6可能太小建议用约束区域面积的百万分之一作为下限。这样既能滤掉退化区域又不会误删正常小区域。还有一个更隐蔽的坑当某个站点紧贴约束边界时它本身的Voronoi多边形在边界外部的那一部分在裁剪后可能变成一个极窄的多边形面积很小但在业务上看它仍然存在。所以过滤前先想清楚需求是“完全不要微小区域”还是“保留但标记出来”这会影响过滤策略。5. 上万点集时的性能实测与绘制策略5.1 从一千点到五万点的耗时变化很多Voronoi图的入门示例都只画几十个点但实际场景里可能遇到几千、几万个站点。我拿随机分布的点做了简单测试用tic/toc统计voronoin的计算耗时结果大致如下点数量计算耗时约1000.002秒1,0000.015秒5,0000.08秒10,0000.18秒50,0001.20秒Matlab底层用的是Qhull库二维Voronoi计算本身很快几万点也只是一两秒的事所以计算阶段通常不是瓶颈。真正卡住的地方在绘制阶段尤其是我前面写的那种逐面片patch的写法画几万个多边形可能要十几秒甚至几十秒。如果需要对很多个点集批量生成我建议先把核心计算写成函数避免在循环里反复重复绘图代码。计算和绘制分开定位问题也会更方便。5.2 绘制阶段比计算阶段更容易卡死逐个调用patch绘制大量多边形时每个patch都是一个图形对象几万个图形对象会把OpenGL渲染器直接拖垮。我踩过这个坑一万个点左右代码跑了快二十秒。解决办法是合并绘制。把有效多边形的顶点预先填充到NaN分隔的大数组里一次性调用patch或者fill。还有一种更省事的思路利用polyshape的数组特性把所有裁剪后的多边形对象放进一个数组然后统一plot。shapes polyshape.empty(); for i 1:length(C) if all(C{i} ~ 1) poly polyshape(V(C{i}, 1), V(C{i}, 2)); polyClipped intersect(poly, boundary); if polyClipped.NumRegions 0 area(polyClipped) 1e-6 shapes(end1) polyClipped; end end end plot(shapes);polyshape数组统一plot的效率远高于逐面片画图而且颜色、边界样式可以统一设置。如果仍然觉得慢就该考虑在画图时降低显示精度比如用reduce函数简化多边形顶点数视觉影响不大但渲染速度能再上一个台阶。6. 数值陷阱与工程化封装建议6.1 重复站点带来的Qhull退化Voronoi计算建立在点集几何关系上如果两个站点坐标完全重复或者距离非常非常近Qhull会判定输入退化轻则报警告重则直接报错。我遇到过坐标数据里混入了重复采集点导致后面所有结果全部错乱的情况。解决方式很简单计算前先去重[P, ~, ic] unique(P, rows);别忘了去重后站点数量变了如果原本还关联了其他业务数据要同步处理索引对应关系。对于距离极近但未完全重复的点可以按最小间距阈值做一次合并比如先对P做四舍五入到小数点后6位然后再unique能减少很多不必要的数值扰动。6.2 坐标尺度过大或过小的影响坐标数值非常大或非常小时双精度浮点数也会在几何计算中损失精度。比如经纬度坐标直接参与计算虽然能出结果但裁剪判断时容易出现边界不稳。更稳妥的做法是先对点集做中心化和平移缩放统一到单位坐标系下完成所有Voronoi计算和裁剪最后再把结果缩放回原坐标系。P0 P - mean(P); scale max(abs(P0), [], all); P_norm P0 / scale; % 计算完成后再乘回scale并加回mean这一步看起来多此一举但能省掉很多说不清的边界误差。尤其是做高精度区域计算时坐标归一化是值得养成的好习惯。6.3 把代码封装成可复用的文件最后说一点工程化的经验。Voronoi图生成的代码如果每次用都重新复制一遍很容易在某个版本里改坏一个参数。我现在的做法是封装成一个函数输入站点坐标和约束边界输出裁剪后的polyshape数组function shapes voronoiClip(P, boundary, minArea) [V, C] voronoin(P); shapes polyshape.empty(); for i 1:length(C) if all(C{i} ~ 1) poly polyshape(V(C{i}, 1), V(C{i}, 2)); poly intersect(poly, boundary); if poly.NumRegions 0 area(poly) minArea shapes(end1) poly; end end end end这样每次调用只需要一行代码核心逻辑也不会被业务代码干扰。如果在实际项目里发现某些站点特别靠近约束边界导致裁剪后面积极小我的习惯是把阈值逻辑做成可选参数不传就不过滤保留完整几何结果方便排查问题。这个小习惯帮我在好几次数据异常时快速定位了问题来源。本文还有配套的精品资源点击获取