3.10.3 沃罗诺伊图
沃罗诺伊图(Voronoi Diagram)是一种空间分割算法，它根据指定的N个胞点将空间分割为N个区域，每个区域中的所有坐标点都与该区域对应的胞点最近。
points2d = np.array([[0.2, 0.1], [0.5, 0.5], [0.8, 0.1],
[0.5, 0.8], [0.3, 0.6], [0.7, 0.6], [0.5, 0.35]])
vo = spatial.Voronoi(points2d)
vo.vertices vo.regions vo.ridge_vertices
-------------------- ----------------- -----------------
[[ 0.5 , 0.045 ], [[-1, 0, 1], [[-1, 0],
[ 0.245 , 0.351 ], [-1, 0, 2], [0, 1],
[ 0.755 , 0.351 ], [], [-1, 1],
[ 0.3375, 0.425 ], [6, 4, 3, 5], [0, 2],
[ 0.6625, 0.425 ], [5, -1, 1, 3], [-1, 2],
[ 0.45 , 0.65 ], [4, 2, 0, 1, 3], [3, 5],
[ 0.55 , 0.65 ]] [6, -1, 2, 4], [3, 4],
[6, -1, 5]] [4, 6],
[5, 6],
[1, 3],
[-1, 5],
[2, 4],
[-1, 6]]
使用voronoi_plot_2d()可以将沃罗诺伊图显示为图表，效果如图3-61(左)所示。图中蓝色小圆点为points2d指定的胞点，红色大圆点表示Voronoi.vertices中的点，图中为每个vertices点标注了下标。由虚线和实线将空间分为7个区域，以虚线为边的区域为无限大的区域，一直向外延伸，全部由实线构成的区域为有限区域。每个区域都以vertices中的点为顶点。
图3-61 沃罗诺伊图将空间分割为多个区域
Voronoi.regions是区域列表，其中每个区域由一个列表(忽略空列表)表示，列表中的整数为vertices中的序号，包含-1的区域为无限区域。例如[6, 4, 3, 5]为图中正中心的那块区域。
Voronoi.ridge_vertices是区域分割线列表，每条分割线由vertices中的两个序号构成，包含-1的分割线为图中的虚线，其长度为无限长。
如果希望将每块区域以不同颜色填充，但由于外围的区域是无限大的，因此无法使用matplotlib绘图，可以在外面添加4个点将整个区域围起来，这样每个points2d中的胞点都对应一个有限区域。在图3-61(右)中，黑圈点为points2d指定的胞点，将空间中与其最近的区域填充成胞点的颜色。
bound = np.array([[-100, -100], [-100, 100],
[ 100, 100], [ 100, -100]])
vo2 = spatial.Voronoi(np.vstack((points2d, bound)))
使用沃罗诺伊图可以解决最大空圆问题：找到一个半径最大的圆，使得其圆心在一组点的凸包区域之内，并且圆内没有任何点。根据图3-61可知最大空圆必定是以vertices为圆心，以与最近胞点的距离为半径的圆中的某一个。为了找到胞点与vertices之间的联系，可以使用Voronoi.point_region属性。point_region[i]是与第i个胞点(points2d[i])最近的区域的编号。例如由下面的输出可知，下标为5的蓝色点与下标为6的区域对应，而由Voronoi.regions[6]可知，该区域由编号为6、2、4的三个vertices点(图中红色的点)构成。因此vertices中与胞点points2d[5]最近的点的序号为[2, 4, 6]。
print vo.point_region
print vo.regions[6]
[0 3 1 7 4 6 5]
[6, -1, 2, 4]
下面是计算最大空圆的程序，效果如图3-62所示。程序中：❶使用pylab.Polygon.contains_point()判断用ConvexHull计算的凸包多边形是否包含vertices中的点，用以在❷处剔除圆心在凸包之外的圆。
图3-62 使用沃罗诺伊图计算最大空圆
❸vertice_point_map是一个字典，它的键为vertices点的下标，值为与其最近的几个points2d点的序号。整个字典使用point_region和regions构建，注意这里剔除了所有在凸包之外的vertices点。
❹对于vertice_point_map中的每一对点，找到距离最大的那对点，即可得出圆心坐标和圆的半径。
from collections import defaultdict
n = 50
np.random.seed(42)
points2d = np.random.rand(n, 2)
vo = spatial.Voronoi(points2d)
ch = spatial.ConvexHull(points2d)
poly = pl.Polygon(points2d[ch.vertices]) ❶
vs = vo.vertices
convexhull_mask = [poly.contains_point(p, radius=0) for p in vs] ❷
vertice_point_map = defaultdict(list) ❸
for index_point, index_region in enumerate(vo.point_region):
region = vo.regions[index_region]
if -1 in region: continue
for index_vertice in region:
if convexhull_mask[index_vertice]:
vertice_point_map[index_vertice].append(index_point)
def dist(p1, p2):
return ((p1-p2)**2).sum()**0.5
max_cicle = max((dist(points2d[pidxs[0]], vs[vidx]), vs[vidx]) ❹
for vidx, pidxs in vertice_point_map.iteritems())
r, center = max_cicle
print "r = ", r, ", center = ", center
r = 0.174278456762 , center = [ 0.46973363 0.59356531]