10.scipy Spatial

1import numpy as np
2import pylab as pl
3from scipy import spatial

1import matplotlib as mpl
2mpl.rcParams['font.sans-serif'] = ['SimHei']

1x = np.sort(np.random.rand(100))
2idx = np.searchsorted(x, 0.5)
3print (x[idx], x[idx - 1]) #距离0.5最近的数是这两个数中的一个

1from scipy import spatial
2np.random.seed(42)
3N = 100
4points = np.random.uniform(-1, 1, (N, 2))
5kd = spatial.cKDTree(points)
6
7targets = np.array([(0, 0), (0.5, 0.5), (-0.5, 0.5), (0.5, -0.5), (-0.5, -0.5)])
8dist, idx = kd.query(targets, 3)

1r = 0.2
2idx2 = kd.query_ball_point(targets, r)
3idx2

1idx3 = kd.query_pairs(0.1) - kd.query_pairs(0.08)
2idx3

 1#%figonly=用cKDTree寻找近旁点
 2x, y = points.T
 3colors = "r", "b", "g", "y", "k"
 4
 5fig, (ax1, ax2, ax3) = pl.subplots(1, 3, figsize=(12, 4))
 6
 7for ax in ax1, ax2, ax3:
 8    ax.set_aspect("equal")
 9    ax.plot(x, y, "o", markersize=4)
10    
11for ax in ax1, ax2:
12    for i in range(len(targets)):
13        c = colors[i]
14        tx, ty = targets[i]
15        ax.plot([tx], [ty], "*", markersize=10, color=c)
16
17for i in range(len(targets)):
18    nx, ny = points[idx[i]].T
19    ax1.plot(nx, ny, "o", markersize=10, markerfacecolor="None", 
20             markeredgecolor=colors[i], markeredgewidth=1)
21    
22    nx, ny = points[idx2[i]].T
23    ax2.plot(nx, ny, "o", markersize=10, markerfacecolor="None", 
24             markeredgecolor=colors[i], markeredgewidth=1)
25    
26    ax2.add_artist(pl.Circle(targets[i], r, fill=None, linestyle="dashed"))
27    
28for pidx1, pidx2 in idx3:
29    sx, sy = points[pidx1]
30    ex, ey = points[pidx2]
31    ax3.plot([sx, ex], [sy, ey], "r", linewidth=2, alpha=0.6)
32    
33ax1.set_xlabel(u"搜索最近的3个近旁点")
34ax2.set_xlabel(u"搜索距离在0.2之内的所有近旁点")
35ax3.set_xlabel(u"搜索所有距离在0.08到0.1之间的点对");

1from scipy.spatial import distance
2dist1 = distance.squareform(distance.pdist(points))
3dist2 = distance.cdist(points, targets)
4print(dist1.shape)
5print(dist2.shape)

1print (dist[:, 0]) # cKDTree.query()返回的与targets最近的距离
2print (np.min(dist2, axis=0))

1dist1[np.diag_indices(len(points))] = np.inf
2nearest_pair = np.unravel_index(np.argmin(dist1), dist1.shape)
3print (nearest_pair, dist1[nearest_pair])

1dist, idx = kd.query(points, 2)
2print (idx[np.argmin(dist[:, 1])], np.min(dist[:, 1]))

1N = 1000000
2start = np.random.uniform(0, 100, N)
3span = np.random.uniform(0.01, 1, N)
4span = np.clip(span, 2, 100)
5end = start + span

1def naive_count_at(start, end, time):
2    mask = (start < time) & (end > time)
3    return np.sum(mask)

 1#%figonly=使用二维K-d树搜索指定区间的在线用户
 2def _():
 3    N = 100
 4    start = np.random.uniform(0, 100, N)
 5    span = np.random.normal(40, 10, N)
 6    span = np.clip(span, 2, 100)
 7    end = start + span
 8
 9    time = 40
10    
11    fig, ax = pl.subplots(figsize=(8, 6))
12    ax.scatter(start, end)
13    mask = (start < time) & (end > time)
14    start2, end2 = start[mask], end[mask]
15    ax.scatter(start2, end2, marker="x", color="red")
16    rect = pl.Rectangle((-20, 40), 60, 120, alpha=0.3)
17    ax.add_patch(rect)
18    ax.axhline(time, color="k", ls="--")
19    ax.axvline(time, color="k", ls="--")
20    ax.set_xlabel("Start")
21    ax.set_ylabel("End")
22    ax.set_xlim(-20, 120)
23    ax.set_ylim(-20, 160)
24    ax.plot([0, 120], [0, 120])
25
26_()

 1class KdSearch(object):
 2    def __init__(self, start, end, leafsize=10):
 3        self.tree = spatial.cKDTree(np.c_[start, end], leafsize=leafsize)
 4        self.max_time = np.max(end)
 5      
 6    def count_at(self, time):
 7        max_time = self.max_time
 8        to_search = spatial.cKDTree([[time - max_time, time + max_time]])
 9        return self.tree.count_neighbors(to_search, max_time, p=np.inf)
10    
11naive_count_at(start, end, 40) == KdSearch(start, end).count_at(40)

1np.random.seed(42)
2points2d = np.random.rand(10, 2)
3ch2d = spatial.ConvexHull(points2d)
4print(ch2d.simplices)
5print(ch2d.vertices)

1#%fig=二维平面上的凸包
2poly = pl.Polygon(points2d[ch2d.vertices], fill=None, lw=2, color="r", alpha=0.5)
3ax = pl.subplot(aspect="equal")
4pl.plot(points2d[:, 0], points2d[:, 1], "go")
5for i, pos in enumerate(points2d):
6    pl.text(pos[0], pos[1], str(i), color="blue")
7ax.add_artist(poly);

1np.random.seed(42)
2points3d = np.random.rand(40, 3)
3ch3d = spatial.ConvexHull(points3d)
4ch3d.simplices.shape

1points2d = np.array([[0.2, 0.1], [0.5, 0.5], [0.8, 0.1],
2                     [0.5, 0.8], [0.3, 0.6], [0.7, 0.6], [0.5, 0.35]])
3vo = spatial.Voronoi(points2d)

1print(vo.vertices); print(vo.regions); print(vo.ridge_vertices)

1bound = np.array([[-100, -100], [-100,  100], 
2                  [ 100,  100], [ 100, -100]])
3vo2 = spatial.Voronoi(np.vstack((points2d, bound)))

 1#%figonly=沃罗诺伊图将空间分割为多个区域
 2fig, (ax1, ax2) = pl.subplots(1, 2, figsize=(9, 4.5))
 3ax1.set_aspect("equal")
 4ax2.set_aspect("equal")
 5spatial.voronoi_plot_2d(vo, ax=ax1)
 6for i, v in enumerate(vo.vertices):
 7    ax1.text(v[0], v[1], str(i), color="red")
 8    
 9for i, p in enumerate(points2d):
10    ax1.text(p[0], p[1], str(i), color="blue")
11
12n = len(points2d)
13color = pl.cm.rainbow(np.linspace(0, 1, n))
14for i in range(n):
15    idx = vo2.point_region[i]
16    region = vo2.regions[idx]
17    poly = pl.Polygon(vo2.vertices[region], facecolor=color[i], alpha=0.5, zorder=0)
18    ax2.add_artist(poly)
19ax2.scatter(points2d[:, 0], points2d[:, 1], s=40, c=color, linewidths=2, edgecolors="k")
20ax2.plot(vo2.vertices[:, 0], vo2.vertices[:, 1], "ro", ms=6)
21
22for ax in ax1, ax2:
23    ax.set_xlim(0, 1)
24    ax.set_ylim(0, 1)

1print (vo.point_region)
2print (vo.regions[6])

1x = np.array([46.445, 263.251, 174.176, 280.899, 280.899, 
2              189.358, 135.521, 29.638, 101.907, 226.665])
3y = np.array([287.865, 250.891, 287.865, 160.975, 54.252,
4              160.975, 232.404, 179.187, 35.765, 71.361])
5points2d = np.c_[x, y]
6dy = spatial.Delaunay(points2d)
7vo = spatial.Voronoi(points2d) 
8print(dy.simplices)
9print(vo.vertices)

 1#%fig=德劳内三角形的外接圆与圆心
 2cx, cy = vo.vertices.T
 3
 4ax = pl.subplot(aspect="equal")
 5spatial.delaunay_plot_2d(dy, ax=ax)
 6ax.plot(cx, cy, "r.")
 7for i, (cx, cy) in enumerate(vo.vertices):
 8    px, py = points2d[dy.simplices[i, 0]]
 9    radius = np.hypot(cx - px, cy - py)
10    circle = pl.Circle((cx, cy), radius, fill=False, ls="dotted")
11    ax.add_artist(circle)
12ax.set_xlim(0, 300)
13ax.set_ylim(0, 300);