本节介绍scipy.spatial模块中提供的K-d树、凸包、沃罗诺伊图、德劳内三角化等空间算法类的使用方法。
一、计算最近旁点
众所周知，对于一维的已排序的数值，可以使用二分法快速找到与指定数值最近的数值。在下面的例子中，在一个升序排序的随机数数组x中，使用numpy.searchsorted()搜索离0.5最近的数。排序算法的时间复杂度为O(NlogN)，而每次二分搜索的时间复杂度为O(logN)。
x = np.sort(np.random.rand(100))
idx = np.searchsorted(x, 0.5)
print x[idx], x[idx - 1] #距离0.5最近的数是这两个数中的一个
0.542258714465 0.492205345391类似的算法可以推广到N维空间，spatial模块提供的cKDTree类使用K-d树快速搜索N维空间中的最近点。在下面的例子中，用平面上随机的100个点创建cKDTree对象，并对其搜索与targets中每个点距离最近的3个点。cKDTree.query()返回两个数组dist和idx，dist[i, :]是距离targets[i]最近的3个点的距离，而idx[i, :]是这些最近点的下标：
from scipy import spatial
np.random.seed(42)
N = 100
points = np.random.uniform(-1, 1, (N, 2))
kd = spatial.cKDTree(points)
targets = np.array([(0, 0), (0.5, 0.5), (-0.5, 0.5), (0.5, -0.5), (-0.5, -0.5)])
dist, idx = kd.query(targets, 3)
dist idx
----------------------------------------- --------------
[[ 0.15188266, 0.21919416, 0.27647793], [[48, 73, 81],
[ 0.09595807, 0.15745334, 0.22855398], [37, 78, 43],
[ 0.05009422, 0.17583445, 0.1807312 ], [79, 22, 92],
[ 0.11180181, 0.16618122, 0.18127473], [35, 58, 6],
[ 0.19015485, 0.19060739, 0.19361173]] [83, 7, 42]]cKDTree.query_ball_point()搜索与指定点在一定距离之内的所有点，它只返回最近点的下标，由于每个目标点的近旁点数不一定相同，因此idx2数组中的每个元素都是一个列表：
r = 0.2
idx2 = kd.query_ball_point(targets, r)
idx2
array([[48], [37, 78], [79, 92, 22], [58, 35, 6], [7, 55, 83, 42]], dtype=object)cKDTree.query_pairs()找到points中距离小于指定值的每一对点，它返回的是一个下标对的集合对象。下面的程序使用集合差集运算找出所有距离在0.08到0.1之间的点对：
idx3 = kd.query_pairs(0.1) - kd.query_pairs(0.08)
idx3
{(1, 46), (3, 21), (3, 82), (3, 95), (5, 16), (9, 30),
(10, 87), (11, 42), (11, 97), (18, 41), (29, 74), (32, 51),
(37, 78), (39, 61), (41, 61), (50, 84), (55, 83), (73, 81)}在图1中，与target中的每个点(用五角星表示)最近的点用与其相同的颜色标识：
图1 用cKDTree寻找近旁点
cKDTree的所有搜索方法都有一个参数p，用于定义计算两点之间距离的函数。读者可以尝试使用不同的p参数，观察图1的变化：
·p = 1：绝对值之和作为距离
·p = 2：欧式距离
·p = np.inf：最大坐标差值作为距离
此外，cKDTree.query_ball_tree()可以在两棵K-d树之间搜索距离小于给定值的所有点对。
distance子模块中的pdist()计算一组点中每对点的距离，而cdist()计算两组点中每对点的距离。由于pdist()返回的是一个压缩之后的一维数组，需要用squareform()将其转换成二维数组。dist1[i, j]是points中下标为i和j的两个点的距离，dist2[i, j]是points[i]和targets[j]之间的距离。
from scipy.spatial import distance
dist1 = distance.squareform(distance.pdist(points))
dist2 = distance.cdist(points, targets)
dist1.shape dist2.shape
----------- -----------
(100, 100) (100, 5)下面使用np.min()在dist2中搜索points中与targets距离最近的点，其结果与cKDTree.query()的结果相同：
print dist[:, 0] # cKDTree.query()返回的与targets最近的距离
print np.min(dist2, axis=0)
[ 0.15188266 0.09595807 0.05009422 0.11180181 0.19015485]
[ 0.15188266 0.09595807 0.05009422 0.11180181 0.19015485]为了找到points中最近的点对，需要将dist1对角线上的元素填充为无穷大：
dist1[np.diag_indices(len(points))] = np.inf
nearest_pair = np.unravel_index(np.argmin(dist1), dist1.shape)
print nearest_pair, dist1[nearest_pair]
(22, 92) 0.00534621024816用cKDTree.query()可以快速找到这个距离最近的点对：在K-d树中搜索它自己包含的点，找到与每个点最近的两个点，其中距离最近的点就是它本身，距离为0，而距离第二近的点就是每个点的最近旁点，然后只需要找到这些距离中最小的那个即可：
dist, idx = kd.query(points, 2)
print idx[np.argmin(dist[:, 1])], np.min(dist[:, 1])
[22 92] 0.00534621024816让我们看一个使用K-d树提高搜索速度的实例。下面的start和end数组保存用户登录和离开网站的时间，对于任意指定的时刻time，计算该时刻在线用户的数量。
N = 1000000
start = np.random.uniform(0, 100, N)
span = np.random.uniform(0.01, 1, N)
span = np.clip(span, 2, 100)
end = start + span下面的naive_count_at()采用逐个比较的方法计算指定时间的在线用户数量：
def naive_count_at(start, end, time):
mask = (start＜time)＆(end＞time)
return np.sum(mask)图2显示了如何使用K-d树实现快速搜索。图中每点的横坐标为start，纵坐标为end。由于end大于start，因此所有的点都在y=x斜线的上方。图中阴影部分表示满足(start＜time)＆(end＞time)条件的区域。该区域中的点数为时刻time时的在线人数。
图2 使用二维K-d树搜索指定区间的在线用户
tree.count_neighbors(other, r, p=2)可以统计tree中到other的距离小于r的点数，其中p为计算距离的范数。距离按照如下公式计算：
当p=2时，上面的公式就是欧几里得空间中的向量长度。当p=∞时，该距离变为各个轴上的最大绝对值：
N∞(x)=｜｜x｜｜∞=max(｜x1｜,｜x2｜…,｜xn｜)
当使用p=∞时可以计算某正方形区域之内的点数。我们将该正方形的中心设置为(time - max_time, time + max_time)，正方形的边长为2 * max_time，即r = max_time。其中max_time为end中的最大值。下面的KdSearch类实现了该算法：
class KdSearch(object):
def __init__(self, start, end, leafsize=10):
self.tree = spatial.cKDTree(np.c_[start, end], leafsize=leafsize)
self.max_time = np.max(end)
def count_at(self, time):
max_time = self.max_time
to_search = spatial.cKDTree([[time - max_time, time + max_time]])
return self.tree.count_neighbors(to_search, max_time, p=np.inf)
naive_count_at(start, end, 40) == KdSearch(start, end).count_at(40)
True下面比较运算时间，由结果可知创建K-d树需要约0.5秒时间，K-d树的搜索速度则为线性搜索的17倍左右。
请读者研究点数N和leafsize参数与创建K-d树和搜索时间之间的关系。
%time ks = KdSearch(start, end, leafsize=100)
%timeit naive_search(start, end, 40)
%timeit ks.count_at(40)
Wall time: 484 ms
100 loops, best of 3: 3.85 ms per loop
1000 loops, best of 3: 221 μs per loop二、凸包
所谓凸包是指N维空间中的一个区域，该区域中任意两点之间的线段都完全被包含在该区域之中，二维平面上的凸多边形就是典型的凸包。ConvexHull可以快速计算包含N维空间中点的集合的最小凸包。下面先看一个二维的例子：points2d是一组二维平面上的随机点，ch2d是这些点的凸包对象。ConvexHull.simplices是凸包的每条边线的两个顶点在points2d中的下标，由于它的形状为(5, 2)，因此凸包由5条线段构成。对于二维的情况，ConvexHull.vertices是凸多边形的每个顶点在points2d中的下标，按逆时针方向的顺序排列。
np.random.seed(42)
points2d = np.random.rand(10, 2)
ch2d = spatial.ConvexHull(points2d)
ch2d.simplices ch2d.vertices
-------------- ---------------
[[2, 5], [5, 2, 6, 1, 0]
[2, 6],
[0, 5],
[1, 6],
[1, 0]]使用matplotlib中的Polygon对象可以绘制如图3所示的多边形。
图3 二维平面上的凸包
poly = pl.Polygon(points2d[ch2d.vertices], fill=None, lw=2, color="r", alpha=0.5)
ax = pl.subplot(aspect="equal")
pl.plot(points2d[:, 0], points2d[:, 1], "go")
for i, pos in enumerate(points2d):
pl.text(pos[0], pos[1], str(i), color="blue")
ax.add_artist(poly)三维空间中的凸包是一个凸多面体，每个面都是一个三角形。在下面的例子中，由simplices的形状可知，所得到的凸包由38个三角形面构成：
np.random.seed(42)
points3d = np.random.rand(40, 3)
ch3d = spatial.ConvexHull(points3d)
ch3d.simplices.shape
(38, 3)下面的程序用TVTK直观地显示凸包，图4中所有的绿色圆球表示points3d中的点，由红色线段构成的三角形面表示凸多面体上的面。没有和红色线段相连的点就是凸包之内的点：
图4 三维空间中的凸包
from scpy2 import vtk_convexhull, vtk_scene, vtk_scene_to_array
actors = vtk_convexhull(ch3d)
scene = vtk_scene(actors, viewangle=22)
%array_image vtk_scene_to_array(scene)
scene.close()如果读者希望采用三维交互界面，可以在Notebook中执行如下代码：
%gui qt
from scpy2 import ivtk_scene
ivtk_scene(actors)三、沃罗诺伊图
沃罗诺伊图(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()可以将沃罗诺伊图显示为图表，效果如图5(左)所示。图中蓝色小圆点为points2d指定的胞点，红色大圆点表示Voronoi.vertices中的点，图中为每个vertices点标注了下标。由虚线和实线将空间分为7个区域，以虚线为边的区域为无限大的区域，一直向外延伸，全部由实线构成的区域为有限区域。每个区域都以vertices中的点为顶点。
图5 沃罗诺伊图将空间分割为多个区域
Voronoi.regions是区域列表，其中每个区域由一个列表(忽略空列表)表示，列表中的整数为vertices中的序号，包含-1的区域为无限区域。例如[6, 4, 3, 5]为图中正中心的那块区域。
Voronoi.ridge_vertices是区域分割线列表，每条分割线由vertices中的两个序号构成，包含-1的分割线为图中的虚线，其长度为无限长。
如果希望将每块区域以不同颜色填充，但由于外围的区域是无限大的，因此无法使用matplotlib绘图，可以在外面添加4个点将整个区域围起来，这样每个points2d中的胞点都对应一个有限区域。在图5(右)中，黑圈点为points2d指定的胞点，将空间中与其最近的区域填充成胞点的颜色。
bound = np.array([[-100, -100], [-100, 100],
[ 100, 100], [ 100, -100]])
vo2 = spatial.Voronoi(np.vstack((points2d, bound)))使用沃罗诺伊图可以解决最大空圆问题：找到一个半径最大的圆，使得其圆心在一组点的凸包区域之内，并且圆内没有任何点。根据图5可知最大空圆必定是以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]下面是计算最大空圆的程序，效果如图6所示。程序中：?使用pylab.Polygon.contains_point()判断用ConvexHull计算的凸包多边形是否包含vertices中的点，用以在?处剔除圆心在凸包之外的圆。
图6 使用沃罗诺伊图计算最大空圆
?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]四、德劳内三角化
德劳内三角化算法对给定的点集合的凸包进行三角形分割，使得每个三角形的外接圆都不含任何点。下面的程序演示了德劳内三角化的效果：
Delaunay对象dy的simplices属性是每个三角形的顶点在points2d中的下标。可以通过三角形的顶点坐标计算其外接圆的圆心，不过这里使用同一点集合的沃罗诺伊图的vertices属性，vertices就是每个三角形的外接圆的圆心。
x = np.array([46.445, 263.251, 174.176, 280.899, 280.899,
189.358, 135.521, 29.638, 101.907, 226.665])
y = np.array([287.865, 250.891, 287.865, 160.975, 54.252,
160.975, 232.404, 179.187, 35.765, 71.361])
points2d = np.c_[x, y]
dy = spatial.Delaunay(points2d)
vo = spatial.Voronoi(points2d)
dy.simplices vo.vertices
------------ --------------------------------
[[8, 5, 7], [[ 104.58977484, 127.03566055],
[1, 5, 3], [ 235.1285 , 198.68143374],
[5, 6, 7], [ 107.83960707, 155.53682482],
[6, 0, 7], [ 71.22104881, 228.39479887],
[0, 6, 2], [ 110.3105 , 291.17642838],
[6, 1, 2], [ 201.40695449, 227.68436282],
[1, 6, 5], [ 201.61895891, 226.21958623],
[9, 5, 8], [ 152.96231864, 93.25060083],
[4, 9, 8], [ 205.40381294, -90.5480267 ],
[5, 9, 3], [ 235.1285 , 127.45701644],
[9, 4, 3]] [ 267.91709907, 107.6135 ]]下面是用delaunay_plot_2d()绘制德劳内三角化的结果，此外还绘制每个外接圆及其圆心。可以看到三角形的所有顶点都不在任何外接圆之内，效果如图7所示：
图7 德劳内三角形的外接圆与圆心?
cx, cy = vo.vertices.T
ax = pl.subplot(aspect="equal")
spatial.delaunay_plot_2d(dy, ax=ax)
ax.plot(cx, cy, "r.")
for i, (cx, cy) in enumerate(vo.vertices):
px, py = points2d[dy.simplices[i, 0]]
radius = np.hypot(cx - px, cy - py)
circle = pl.Circle((cx, cy), radius, fill=False, ls="dotted")
ax.add_artist(circle)
ax.set_xlim(0, 300)
ax.set_ylim(0, 300)