沃罗诺伊图算法实战：从暴力法到KDTree与扫描线

1. 从帝国边界到算法落地沃罗诺伊图到底在解决什么问题第一次看到帝国边界划分这个说法很多人会以为这是历史或地理话题。其实它描述的是一个非常经典的几何问题假设地图上有若干个权力中心城市、据点、基站、门店每个中心都想控制离自己最近的那片区域那么最终形成的势力范围图就是沃罗诺伊图Voronoi Diagram。每个中心对应一个细胞Cell细胞内的任意一点到该中心的距离都小于到其他任何中心的距离。这个问题的核心价值在于它把最近邻归属这件事从模糊的直觉变成了可计算、可绘制、可验证的几何结构。你只要给出若干种子点算法就能自动帮你把整个平面切分干净不留缝隙、不重叠。听起来简单但真正动手实现时会发现边界到底怎么画、无穷远的区域怎么处理、退化情况多个点共线、重合怎么兜底全是坑。这篇是2/2也就是下半部分。上半部分我们聊了沃罗诺伊图的数学定义、对偶关系和德劳内三角剖分互为对偶以及它在自然界和城市格局里的直观体现。这一篇要解决的是工程落地怎么把一张沃罗诺伊图真正算出来、画出来、用起来。我会从最朴素的暴力法讲起一路讲到分治法、扫描线法再给出可直接运行的代码最后聊聊它在选址、路径规划、生成艺术里的真实用法。适合谁看如果你已经理解了沃罗诺伊图是什么但卡在怎么算、怎么画、怎么用这一步这篇就是为你写的。哪怕你只是想在项目里画一张好看的分区图或者想搞明白游戏里那种自然裂纹、细胞纹理是怎么生成的也能直接抄作业。提示本文所有代码基于 Python依赖 numpy 和 matplotlib不需要任何额外几何库就能跑通基础版本。进阶部分会提到 scipy.spatial 的 Voronoi 类但核心逻辑我会手写一遍方便你理解原理。2. 暴力法先跑通为什么逐点比较是最该先写的版本2.1 暴力法的核心逻辑与像素化实现很多人一上来就想找现成的库结果调出来的图边界是断的、颜色是乱的反而更懵。我的建议永远是先用最笨的办法跑通一遍。暴力法的思路简单到不能再简单——把整个画布切成一个个像素点对每个像素计算它到所有种子点的距离谁最近就归谁。用生活化的类比想象你在一张方格纸上画了很多个点然后拿一支笔一格一格地涂色每格涂成离它最近的那个点的颜色。涂完沃罗诺伊图就出来了。这个方法的时间复杂度是 O(像素数 × 种子点数)听起来很慢但对于几百个种子点、一张 800×800 的图现代计算机几秒钟就能跑完完全够用。import numpy as np import matplotlib.pyplot as plt def voronoi_bruteforce(seeds, width800, height800): # seeds: (N, 2) 的数组每行是一个种子点的 (x, y) xs np.arange(width) ys np.arange(height) xx, yy np.meshgrid(xs, ys) # 生成所有像素坐标 # 展平成 (H*W, 2) points np.stack([xx.ravel(), yy.ravel()], axis1).astype(float) # 计算每个像素到每个种子的距离形状 (H*W, N) dists np.linalg.norm(points[:, None, :] - seeds[None, :, :], axis2) # 取最近种子的索引 labels np.argmin(dists, axis1).reshape(height, width) return labels # 随机生成 20 个种子点 np.random.seed(42) seeds np.random.rand(20, 2) * 800 labels voronoi_bruteforce(seeds) plt.figure(figsize(8, 8)) plt.imshow(labels, cmaptab20, originlower) plt.scatter(seeds[:, 0], seeds[:, 1], cblack, s30, markerx) plt.title(Voronoi by brute force) plt.axis(off) plt.show()这段代码跑出来的图颜色块就是各个帝国的领地黑叉是权力中心。你会发现边界是锯齿状的因为像素是离散的。想要边界更平滑把分辨率调高就行代价是计算量按平方增长。2.2 暴力法的三个隐藏坑第一个坑是坐标顺序。numpy 的 meshgrid 默认返回的是 (行, 列)对应到图像是 (y, x)而种子点通常按 (x, y) 给。如果你不统一画出来的点会整体转置看起来对但就是怪。我习惯在生成 points 时显式写成[xx.ravel(), yy.ravel()]保证和 seeds 的 (x, y) 一致。第二个坑是距离度量。默认用欧氏距离没问题但如果你做的是城市街区划分曼哈顿距离L1可能更符合实际因为道路是横平竖直的。改法很简单把np.linalg.norm换成np.abs(...).sum(axis2)即可。不同距离度量会得到完全不同的图这一点在选址场景里很关键。第三个坑是内存爆炸。points[:, None, :] - seeds[None, :, :]这一步会生成一个 (H*W, N, 2) 的中间数组。800×800 的图配 100 个种子就是 800×800×100×2 ≈ 1.28 亿个浮点数约 1GB 内存。种子一多就崩。解决办法是分块计算或者直接用 scipy 的 KDTree 查询最近邻把复杂度从 O(N) 降到 O(log N)。注意暴力法适合验证理解和做小规模 Demo生产环境请务必换用 KDTree 或直接调用成熟的 Voronoi 库。我见过有人拿暴力法去处理 4K 分辨率的图结果内存直接打满进程被系统杀掉。3. 从 O(N) 到 O(log N)KDTree 加速与 scipy 实战3.1 KDTree 为什么能加速最近邻查询暴力法慢在每个像素都要和所有种子比一遍。KDTree 的思路是先把种子点组织成一棵二叉树按坐标轴交替切分空间。查询某个像素的最近种子时从根节点往下走边走边排除掉那些肯定不可能更近的子树。这样平均只需要比较 log(N) 个种子而不是 N 个。用生活类比你要在一个城市里找最近的便利店。暴力法是挨家挨户问距离KDTree 是先看你在哪个区再在区里找再在街道里找层层缩小范围。种子点分布越均匀KDTree 的效果越好如果所有点挤在一条线上树会退化成链表加速效果打折。from scipy.spatial import cKDTree def voronoi_kdtree(seeds, width800, height800): tree cKDTree(seeds) xs np.arange(width) ys np.arange(height) xx, yy np.meshgrid(xs, ys) points np.stack([xx.ravel(), yy.ravel()], axis1).astype(float) # query 返回 (距离, 索引) dists, labels tree.query(points, k1) return labels.reshape(height, width)实测下来同样 800×800、100 个种子暴力法要 3 到 5 秒KDTree 只要 0.2 秒左右快了十几倍。种子越多差距越明显。3.2 scipy.spatial.Voronoi 的正确打开方式如果你不需要像素级的归属图而是想要矢量边界一堆线段和多边形那就该用scipy.spatial.Voronoi。它直接返回顶点、边、区域等信息画出来的边界是光滑的直线段不是锯齿。from scipy.spatial import Voronoi, voronoi_plot_2d points np.random.rand(30, 2) vor Voronoi(points) fig, ax plt.subplots(figsize(8, 8)) voronoi_plot_2d(vor, axax, show_pointsTrue, show_verticesFalse, line_colorsorange) ax.set_xlim(0, 1) ax.set_ylim(0, 1) plt.show()但这里有个大坑scipy 的 Voronoi 默认不处理无穷远的区域。位于点集凸包上的种子它的细胞会延伸到无穷远scipy 用 -1 表示这些区域的顶点索引。你直接画会发现边缘的细胞是开口的图不封闭。解决办法是手动在点集外围加一圈虚拟点把无穷远区域框回来或者用vor.ridge_points和vor.ridge_vertices自己重建边界。我常用的技巧是在原始点集的包围盒外扩一圈均匀撒 8 到 16 个虚拟点算完 Voronoi 后只保留原始点对应的细胞。这样边缘就封闭了而且虚拟点离得够远不会影响内部结构。3.3 两种方案的选型对照方案输出形式适用场景速度边界质量暴力法像素标签图教学、小图、自定义距离慢锯齿KDTree像素标签图大图、自定义距离、批量查询快锯齿scipy.Voronoi矢量顶点与边精确几何、路径规划、面积计算快光滑选型逻辑很直接要图就用前两个要几何结构就用第三个。如果你既要图又要精确面积可以先用 scipy 算出多边形再用 shapely 求面积最后用 matplotlib 填充颜色。4. 分治与扫描线理解工业级算法的底层思路4.1 分治法把大问题切成两半工业级的沃罗诺伊图算法主流是分治法和扫描线法。分治法的思路是把所有种子点按 x 坐标排序从中间切成左右两半分别递归求各自的沃罗诺伊图最后把两张图缝合起来。缝合的关键是找一条分割链merge curve它由左右两侧种子点之间的垂直平分线片段组成。为什么分治能到 O(N log N)因为每次递归把问题规模减半递归深度是 log N每层的缝合操作是线性的。这个复杂度已经是最优的了因为任何算法都至少要读一遍所有点。分治法的难点全在缝合这一步。你要从无穷远开始沿着分割链往下走每走一步都要判断当前是哪个左侧点和哪个右侧点在争夺边界。实现起来代码量不小但理解了它对偶的德劳内三角剖分就会清晰很多——德劳内三角剖分有成熟的分治实现而沃罗诺伊图就是它的对偶两者可以互相转换。4.2 扫描线法一条线扫过去事件驱动扫描线法Fortune 算法更巧妙。想象一条水平线从上往下扫扫过的区域里每个种子点都生长出一个抛物线形状的波前。这些抛物线的下包络beach line就是当前已确定边界的轮廓。当扫描线遇到新种子点或者两条抛物线相交形成一个沃罗诺伊顶点就触发一个事件更新数据结构。Fortune 算法用平衡二叉树维护 beach line用优先队列维护事件整体也是 O(N log N)。它的优势是在线处理点可以一个一个来不需要一次性拿到全部。这在流式数据场景里很有用比如实时更新的基站选址。不过说实话除非你要自己写几何库否则没必要手撸 Fortune 算法。它的边界情况三个点共线、四个点共圆处理起来极其繁琐一个符号错误就全盘皆输。我的建议是理解思想用现成库。scipy、CGAL、shapely 都提供了可靠的实现。4.3 退化情况算法工程师的噩梦不管用哪种算法退化情况都是必须面对的。常见的退化有三类点重合两个种子点坐标完全相同。这时它们的细胞面积都是零算法可能除零或死循环。处理办法是去重或者给重合点加一个极小的随机偏移。多点共线三个以上种子点在同一直线上。这时会出现退化边垂直平分线互相平行没有交点。需要特殊分支处理。四点共圆四个种子点恰好在同一个圆上。这时沃罗诺伊图会出现一个四岔路口四条边交于一点。浮点误差会让这个点分裂成两个很近的点导致图出现微小裂缝。提示在实际项目里我习惯在输入种子点后先做一次抖动jitter给每个坐标加上 1e-6 量级的随机扰动。这能规避 90% 以上的退化问题代价是边界位置有极微小的偏移肉眼完全看不出来。5. 沃罗诺伊图的真实应用选址、路径与生成艺术5.1 选址与势力范围划分回到帝国边界这个比喻。现实中连锁门店的配送范围、外卖平台的骑手分区、通信基站的信号覆盖本质上都是沃罗诺伊问题。给定若干服务点每个点负责离它最近的区域这样总配送距离最短。但纯沃罗诺伊有个问题它假设所有中心权重相同。现实中大城市的门店覆盖范围应该比小城市大。解决办法是加权沃罗诺伊图Weighted Voronoi也叫乘权沃罗诺伊或加权沃罗诺伊。加权版本里距离变成d(p, s) - w_s权重大的点能抢到更多区域。实现时只需在距离计算里减去权重其余逻辑不变。def weighted_voronoi(seeds, weights, width800, height800): xs np.arange(width) ys np.arange(height) xx, yy np.meshgrid(xs, ys) points np.stack([xx.ravel(), yy.ravel()], axis1).astype(float) dists np.linalg.norm(points[:, None, :] - seeds[None, :, :], axis2) dists dists - weights[None, :] # 加权重 labels np.argmin(dists, axis1).reshape(height, width) return labels这个改动虽小但效果立竿见影。权重可以按人口、销售额、订单量来设让分区更符合业务实际。5.2 路径规划与避障在机器人路径规划里沃罗诺伊图有一个经典用法广义沃罗诺伊图Generalized Voronoi Diagram。把障碍物的顶点和边当作种子生成的沃罗诺伊边就是离所有障碍物都尽量远的路径。机器人沿着这些边走天然就能避开障碍因为每条边到最近障碍的距离是局部最大的。这个思路在无人机航线、自动驾驶换道决策里都有应用。它的优点是路径安全裕度高缺点是路径往往比较绕不是最短路径。实际工程里通常先用沃罗诺伊图生成一条安全走廊再用优化算法在走廊内拉直。5.3 生成艺术与程序化纹理游戏和设计领域沃罗诺伊图是生成自然纹理的利器。细胞裂纹、干裂土地、长颈鹿斑纹、玻璃碎片都能用它做出来。做法是随机撒种子点算沃罗诺伊图然后给每个细胞填充略有差异的颜色再叠加噪声和边缘高光。更进一步可以用劳埃德松弛Lloyds Relaxation让细胞变得均匀。做法是反复把每个种子点移到它当前细胞的质心迭代几十次后细胞会变得像蜂巢一样规整。这个技巧在生成有机但均匀的纹理时特别好用。def lloyd_relaxation(seeds, iterations10, width800, height800): for _ in range(iterations): labels voronoi_kdtree(seeds, width, height) new_seeds [] for i in range(len(seeds)): mask (labels i) if mask.sum() 0: new_seeds.append(seeds[i]) continue ys, xs np.where(mask) new_seeds.append([xs.mean(), ys.mean()]) seeds np.array(new_seeds) return seeds跑完松弛再画出来的图就非常舒服没有特别大或特别小的细胞视觉上很平衡。6. 手写一个完整的沃罗诺伊可视化工具6.1 需求拆解与模块划分把前面所有东西串起来我们做一个完整的小工具功能包括随机或手动指定种子点、支持加权、支持劳埃德松弛、输出彩色分区图并标注种子点。模块分三块种子管理、计算核心、可视化。种子管理负责生成、去重、抖动。计算核心用 KDTree 做最近邻查询支持权重。可视化用 matplotlib 的 imshow 加散点颜色用 tab20 或 hsv 色图循环。6.2 完整代码与参数说明import numpy as np import matplotlib.pyplot as plt from scipy.spatial import cKDTree class VoronoiTool: def __init__(self, width800, height800, seed42): self.width width self.height height self.rng np.random.default_rng(seed) self.seeds None self.weights None def random_seeds(self, n, jitter1e-6): pts self.rng.random((n, 2)) * [self.width, self.height] pts self.rng.normal(0, jitter, pts.shape) # 抖动防退化 self.seeds pts self.weights np.zeros(n) return self def set_weights(self, weights): self.weights np.asarray(weights, dtypefloat) return self def relax(self, iterations10): for _ in range(iterations): labels self._compute() new_seeds [] for i in range(len(self.seeds)): mask (labels i) if mask.sum() 0: new_seeds.append(self.seeds[i]) continue ys, xs np.where(mask) new_seeds.append([xs.mean(), ys.mean()]) self.seeds np.array(new_seeds) return self def _compute(self): tree cKDTree(self.seeds) xs np.arange(self.width) ys np.arange(self.height) xx, yy np.meshgrid(xs, ys) points np.stack([xx.ravel(), yy.ravel()], axis1).astype(float) # 加权查询多个候选再按权重修正 dists, idx tree.query(points, kmin(len(self.seeds), 8)) if dists.ndim 1: dists dists[:, None] idx idx[:, None] adjusted dists - self.weights[idx] best np.argmin(adjusted, axis1) labels idx[np.arange(len(best)), best] return labels.reshape(self.height, self.width) def plot(self, cmaptab20): labels self._compute() plt.figure(figsize(8, 8)) plt.imshow(labels, cmapcmap, originlower, interpolationnearest) plt.scatter(self.seeds[:, 0], self.seeds[:, 1], cblack, s25, markerx, linewidths1.5) plt.axis(off) plt.tight_layout() plt.show() # 使用示例 tool VoronoiTool(width600, height600) tool.random_seeds(25).relax(iterations8) tool.plot()参数说明width/height决定分辨率越大越精细但越慢jitter是防退化抖动幅度1e-6 足够relax的迭代次数一般 5 到 15 次太多会过度规整失去自然感k8是加权查询的候选数权重差异大时可以调大。6.3 实测效果与调参心得实测下来600×600 配 25 个种子松弛 8 次整个流程不到 2 秒。松弛后的细胞大小非常均匀适合做规划感强的图不松弛则更自然适合做有机纹理。调参上我踩过的坑权重不能设得太大。如果某个点的权重远超其他点它会吞掉大半个图其他点被挤到边缘看起来像 bug。经验是权重差异控制在距离尺度的 20% 以内比如图宽 600权重差不超过 120。另一个心得是颜色映射的选择。tab20 适合 20 个以内的细胞超过 20 个颜色会重复相邻细胞可能同色边界就糊了。细胞多的时候改用hsv或nipy_spectral或者干脆用随机颜色字典保证相邻不同色。7. 那些文档不会告诉你的踩坑记录7.1 边界不封闭的排查链路第一次用 scipy.Voronoi 画图我发现边缘的细胞是开口的图看起来像被啃了一口。排查过程是这样的先打印vor.regions发现有些区域的顶点列表里出现了 -1再查文档确认 -1 代表无穷远顶点然后尝试手动加虚拟点发现虚拟点太近会挤压内部结构太远又不起作用。最终的解决方案是计算原始点集的包围盒在盒外按 1.5 倍盒宽的距离均匀撒 12 个虚拟点。这样既封闭了边缘又不影响内部。这个1.5 倍是我试出来的经验值太小会干扰太大则虚拟点自身的细胞会异常大虽然不影响原始点但画图时如果不过滤会很难看。7.2 浮点误差导致的裂缝有一次做面积统计发现所有细胞面积加起来比总面积少了 0.3%。查了半天发现是四点共圆导致的微小裂缝。浮点计算里本该交于一点的四条边实际交成了两个相距 1e-10 的点中间留了一条极细的缝。解决办法有两个一是加抖动从源头避免共圆二是后处理时做焊接把距离小于阈值的顶点合并。我一般两个都用抖动防患于未然焊接兜底。7.3 性能优化的三个层次性能优化我总结了三层第一层是换算法暴力换 KDTree第二层是降分辨率先低分辨率算归属再对边界像素做精细判定第三层是并行化把画布切成若干块多进程同时算。第二层特别实用。比如 4K 图先用 400×400 算一遍找到所有边界附近的像素只对这些像素做高精度距离计算。这样能省 90% 以上的计算量视觉上几乎看不出差别。注意并行化时要注意随机种子的管理每个进程用独立的种子否则抖动会重复退化问题反而更严重。8. 从沃罗诺伊图延伸出去的两个方向沃罗诺伊图本身已经够用了但如果你想把这件事做得更深有两个方向值得探索。第一个是德劳内三角剖分。它和沃罗诺伊图互为对偶沃罗诺伊图的每个顶点对应德劳内三角形的一个外接圆圆心每条边对应一条垂直平分线。理解了这层关系你就能在两者之间自由转换。德劳内三角剖分在网格生成、地形建模里用得极多而且它的算法比如 Bowyer-Watson比沃罗诺伊更容易实现。第二个是约束沃罗诺伊图。普通沃罗诺伊图的边界是直线但现实中边界可能是河流、山脉、行政线。约束沃罗诺伊图允许你指定某些边必须存在算法会在满足约束的前提下最小化偏离。这在真实地图分区里非常有用实现上通常要借助计算几何库。我在实际项目里的体会是先把无约束版本吃透再上约束。无约束版本的所有坑退化、边界、性能你都会在约束版本里再遇到一遍而且更难排查。基础打牢了后面都是水到渠成。最后分享一个小技巧如果你只是想快速看效果不想写代码很多在线工具和绘图软件都内置了沃罗诺伊功能。但如果你想真正掌控参数、做批量处理、集成到自己的系统里手写一遍绝对值得。我前后写了三版每一版都对最近邻归属这件事理解更深一层。第一版只会暴力第二版学会用树第三版才搞明白加权和松弛的配合。这个过程没有捷径但每一步的收获都很实在。