用Python生成Voronoi图：从算法原理到代码实战（附完整源码）
用Python生成Voronoi图：从算法原理到代码实战（附完整源码）
在数据可视化与计算几何领域，Voronoi图以其独特的空间分割特性成为算法工程师的必备工具。这种将平面划分为若干区域的数据结构，每个区域包含距离特定生成点最近的所有空间点，在游戏开发、地理信息系统、机器学习等场景中展现出惊人的实用性。本文将彻底拆解其数学原理，并手把手带您实现三种不同效率的Python实现方案。
1. Voronoi图的数学本质与核心算法
1.1 空间分割的数学定义
给定平面上的点集P={p₁,p₂,...,pₙ}，Voronoi图将空间划分为n个单元，满足：
V(pᵢ) = {x | d(x,pᵢ) ≤ d(x,pⱼ), ∀j≠i}其中d(x,y)表示距离度量函数。当采用欧式距离时，单元边界由直线段构成；使用曼哈顿距离则会产生矩形边界。下图展示了不同距离函数下的单元形态差异：
| 距离类型 | 公式 | 单元特征 |
|---|---|---|
| 欧式距离(L₂) | √((x₁-x₂)²+(y₁-y₂)²) | 凸多边形 |
| 曼哈顿距离(L₁) | x₁-x₂ | |
| 切比雪夫距离 | max( | x₁-x₂ |
1.2 Delaunay三角剖分的对偶关系
Delaunay三角剖分与Voronoi图存在精确的对偶性：
- 连接相邻Voronoi单元生成点构成Delaunay三角形
- 每个Delaunay三角形的外接圆不包含其他生成点（空圆特性）
- 三角剖分最大化最小内角，避免狭长三角形
这种对偶关系使得我们可以通过以下步骤高效构建Voronoi图：
- 计算点集的Delaunay三角剖分
- 找出每个三角形的外接圆圆心
- 连接相邻三角形的外接圆圆心形成Voronoi边
1.3 算法效率对比
不同构建方法的复杂度差异显著：
# 暴力法伪代码 O(n²) def brute_force_voronoi(points): for each pixel in plane: find closest point assign pixel to point's region # Fortune算法伪代码 O(n log n) def fortune_algorithm(points): initialize sweep line and event queue while events exist: handle site/circle event update beach line and voronoi edges2. 基础实现：纯Python暴力生成法
2.1 核心逻辑实现
以下代码展示了最直观的逐像素判定方法：
import numpy as np import matplotlib.pyplot as plt def naive_voronoi(width, height, points): """生成分辨率width×height的Voronoi图""" image = np.zeros((height, width, 3)) for y in range(height): for x in range(width): distances = [np.hypot(x-px, y-py) for px, py in points] closest = np.argmin(distances) image[y,x] = colors[closest] return image # 生成50个随机点 points = np.random.rand(50, 2) * [800, 600] colors = np.random.rand(len(points), 3) plt.imshow(naive_voronoi(800, 600, points)) plt.scatter(points[:,0], points[:,1], c='white', s=10) plt.show()注意：该方法在100×100分辨率下需1万次距离计算，实际项目应避免直接使用
2.2 性能优化技巧
即使采用暴力法，仍有优化空间：
- 距离计算优化：用平方距离代替实际距离避免开方运算
- 区域剪枝：利用空间索引快速排除明显不可能是最近邻的点
- 并行计算：使用multiprocessing分块处理图像区域
# 优化后的距离计算 def squared_dist(a, b): return (a[0]-b[0])**2 + (a[1]-b[1])**23. 科学计算库加速：Scipy实践方案
3.1 使用scipy.spatial模块
Scipy提供的Voronoi类基于Qhull实现，算法复杂度为O(n log n)：
from scipy.spatial import Voronoi, voronoi_plot_2d points = np.random.rand(30, 2) vor = Voronoi(points) fig = voronoi_plot_2d(vor) plt.scatter(points[:,0], points[:,1], c='red', s=20) plt.show()3.2 处理无限区域
Scipy生成的Voronoi对象可能包含无限延伸的边，需要特殊处理：
def adjust_bounds(ax, points): """自动调整绘图边界""" margin = 0.1 * points.ptp(axis=0) xy_min = points.min(axis=0) - margin xy_max = points.max(axis=0) + margin ax.set_xlim(xy_min[0], xy_max[0]) ax.set_ylim(xy_min[1], xy_max[1]) fig, ax = plt.subplots() voronoi_plot_2d(vor, ax=ax, show_vertices=False) adjust_bounds(ax, points)3.3 获取单元多边形
提取有限Voronoi单元的实际坐标：
def get_finite_polygons(vor): """返回有限Voronoi单元的顶点坐标""" finite_regions = [] for reg in vor.regions: if not reg or -1 in reg: continue finite_regions.append(vor.vertices[reg]) return finite_regions polygons = get_finite_polygons(vor)4. 高级应用：Lloyd松弛与艺术生成
4.1 Lloyd算法实现
通过迭代优化使Voronoi单元更均匀：
def lloyd_relaxation(points, bounds, iterations=5): """执行Lloyd松弛算法""" for _ in range(iterations): vor = Voronoi(points) new_points = [] for i, reg in enumerate(vor.regions): if not reg or -1 in reg: continue polygon = vor.vertices[reg] centroid = np.mean(polygon, axis=0) new_points.append(centroid) points = np.clip(new_points, bounds[0], bounds[1]) return points # 初始随机点 points = np.random.rand(20, 2) * [1, 1] bounds = np.array([[0,0], [1,1]]) # 执行5次松弛 relaxed_points = lloyd_relaxation(points, bounds)4.2 艺术化Voronoi图案
结合颜色映射创建艺术效果：
def artistic_voronoi(size, points, color_func): """生成带艺术效果的Voronoi图""" from matplotlib.colors import to_rgb img = np.zeros((size[1], size[0], 3)) for y in range(size[1]): for x in range(size[0]): dists = np.sum((points - [x,y])**2, axis=1) closest = np.argmin(dists) img[y,x] = to_rgb(color_func(points[closest])) return img def rainbow_color(point): """根据点坐标生成彩虹色""" return plt.cm.hsv((point[0]+point[1])%1) art = artistic_voronoi((800,600), relaxed_points*[800,600], rainbow_color) plt.imshow(art)5. 工程实践中的性能陷阱
5.1 内存优化策略
处理百万级点时需要特殊技巧：
- 分块处理：将空间划分为网格单独处理
- 稀疏表示：只存储边界信息而非完整像素矩阵
- GPU加速：使用CUDA实现并行距离计算
# 使用sparse矩阵示例 from scipy.sparse import dok_matrix def sparse_voronoi(size, points): """稀疏矩阵存储边界像素""" edge_map = dok_matrix((size[1], size[0]), dtype=bool) for i in range(len(points)-1): for j in range(i+1, len(points)): # 计算两点垂直平分线 mid = (points[i] + points[j]) / 2 normal = np.array([points[i][1]-points[j][1], points[j][0]-points[i][0]]) # 标记边界像素... return edge_map5.2 动态Voronoi图更新
实时应用中需要增量更新：
class DynamicVoronoi: def __init__(self, initial_points): self.vor = Voronoi(initial_points) def add_point(self, new_point): # 增量更新数据结构 pass def remove_point(self, index): # 处理点删除 pass在游戏开发中，我曾遇到需要每帧更新Voronoi图的情况。最终采用分区策略，将静态和动态点分开处理，将计算量降低了70%。