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 一句话解释Voronoi

Voronoi图，又叫泰森多边形（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区域');

这段代码的关键点在于：

1. 先把无穷远的顶点用一个足够大的坐标值（我用的100，只要比边界范围大很多就行）替换掉，否则polyshape会报错。
2. 用intersect求多边形交集，自动完成裁剪。
3. 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('存在重复点，请检查数据'); end

6.3Inf处理不当导致绘图失败

新手最容易在patch这一步踩坑。要记住：patch或者fill接受NaN来控制线段断开，但Inf会让绘制彻底失败。所以遇到无界区域时，要么跳过填充，要么替换Inf后再处理。

我建议在编写代码时，就在数据预处理阶段做一次统一检查：

for i = 1:length(C) v = V(C{i}, :); if any(isinf(v(:))) % 标记为无界区域，后续单独处理 end end

6.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 end

7.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和带权重Voronoi

8.1 三维Voronoi

voronoin天生支持高维数据。只要把输入从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 diagram），Matlab没有内置函数，但可以通过改进的距离公式实现。

比较简单的近似实现方式：把每个种子点的坐标复制多份，权重大的点复制次数多，然后对增广点集生成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)看一下面积直方图，很多时候能从点位分布中找到肉眼看不出来的规律。

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