求若干坐标点最近邻区间分界线:用Python计算Voronoi图

问题背景

Fig.1 Voronoi Map

在地理空间数据处理中我们可能遇到这样的问题:
如图1所示, 在某一区域中有N只狮子(Xn, n > 1),每只狮子都有自己的地盘,且每只狮子的制霸能力相同(即狮子和狮子之间画中垂线而统治)。
目标是计算每只狮子制霸的区域和边界。

此问题还可以推广到飞机场、地铁站管辖的区域等问题。

Voronoi图

此问题可以通过一个叫做Voronoi的数学模型解决。

Voronoi图是一组连续多边形组成,多边形的边界是由连接的垂直平分线组成。M 个在工平面上有区则的点。按照最近邻原则划分平面,每一个点与它最近邻的区域关联,与每个点相关联的区城(成多边用是唯一的,它由这些点的空网分布所决定。
Ref.: 刘耀林.土地信息系统:中国农业出版社,2011

关于Voronoi图的详细介绍请参考 百度百科

给定若干点计算Voronoi图的算法,一般是通过叫做 “Delaunay三角剖分”。关于 “Delaunay三角剖分” 请参考博文 voronoi图的和Delaunay三角剖分,这里也不再赘述。

这里主要介绍Voronoi图的Python实现。

Voronoi图的Python实现

可以通过Python的SciPy工具轻松计算Voronoi图:

文档点这里

测试环境

  • Python 3.6.1
  • SciPy 0.18.1

计算:

1
2
3
4
5
6
import numpy as np
from scipy.spatial import Voronoi
# 初始化种子点
points = np.array([[0, 0], [0, 1], [0, 2], [1, 0], [1, 1], [1, 2], [2, 0], [2, 1], [2, 2]])
# 计算Voronoi图
vor = Voronoi(points=points)

输入:

  • points : 种子点集。ndarray float数组。shape (点数量, 维度)。

可视化:

1
2
3
4
from scipy.spatial import voronoi_plot_2d
import matplotlib.pyplot as plt
voronoi_plot_2d(vor)
plt.show()

输出:

输入种子点[[0, 0], [0, 1], [0, 2], [1, 0], [1, 1], [1, 2], [2, 0], [2, 1], [2, 2]]后输出如下:

可视化结果

各输出元素


图像引用自: 《scipy Voronoi》

输出

Voronoi([[0, 0], [0, 1], [0, 2], [1, 0], [1, 1], [1, 2], [2, 0], [2, 1], [2, 2]])

  • max_bound: 种子点的最大边界
  • min_bound: 种子点的最小边界
  • ndim: 种子点维度
  • npoints: 种子点数量
  • point_region: 种子点所对应的区域下标
  • points: 种子点
  • regions: Voronoi区域
  • ridge_dict: 分界线字典
  • ridge_points: 分界线交点
  • ridge_vertices:
  • vertices: 分界线交点,顶点

关系

关系说明的时候从ridge_points开始,如下图

输出

  1. 在ridge_points中,每一个ridge_index对应两个point_index,这两个点的连线就是相邻点的连接线
  2. 在points里面,每一个point_index对应一个point_coords
1
2
3
4
5
ridge_points = vor.ridge_points
points = vor.points

for ridge_point in ridge_points:
plt.plot(points[ridge_point][:,0], points[ridge_point][:,1], "-", label="neighbor line")

相邻点的连接线

  1. 在point_region里面,每一个point_index对应一个region_index
  2. 在regions里面,每一个region_index对应n各vertice_index
  3. 在vertices里面,每一个vertice_index对应一个vertice_coords
1
2
3
4
5
6
7
8
9
10
11
point_region = vor.point_region
regions = vor.regions
vertices = vor.vertices
for point_index, region_index in enumerate(point_region):
point = points[point_index]
vertice_index = regions[region_index]
for i in vertice_index:
if i != -1: #注意除去vertice_index为0的情况,表示该种子点位于边界
print(vertices[i])
plt.plot(vertices[i], "o")
print("")
  1. 在ridge_vertices里面,每一个ridge_index对应两个vertice_index

题外话:通过Mashgrid计算K-最近邻(kNN)的判别边界

Voronoi可以认为是当k=1时的kNN模型,求其判别区域。
当k > 1时,一般是将区域分割成若干MashGrid,再分别对每个MashGrid归为其中一类进行可视化。
请参考:
Graph k-NN decision boundaries in Matplotlib

k-NN decision boundaries