空间数据结构和算法 (scipy.spatial)#
scipy.spatial 可以计算一组点的三角剖分、Voronoi图和凸包，通过利用 Qhull 库实现。
此外，它还包含用于最近邻点查询的 KDTree 实现，以及在各种度量中进行距离计算的实用工具。
Delaunay 三角剖分#
Delaunay 三角剖分是将一组点划分为一组不重叠的三角形，使得没有任何点位于任何三角形的外接圆内。实际上，这种三角剖分倾向于避免小角度的三角形。
可以使用 scipy.spatial 计算 Delaunay 三角剖分，如下所示：
>>> from scipy.spatial import Delaunay
>>> import numpy as np
>>> points = np.array([[0, 0], [0, 1.1], [1, 0], [1, 1]])
>>> tri = Delaunay(points)
我们可以将其可视化：
>>> import matplotlib.pyplot as plt
>>> plt.triplot(points[:,0], points[:,1], tri.simplices)
>>> plt.plot(points[:,0], points[:,1], 'o')
并添加一些进一步的装饰：
>>> for j, p in enumerate(points):
... plt.text(p[0]-0.03, p[1]+0.03, j, ha='right') # 标注点
>>> for j, s in enumerate(tri.simplices):
... p = points[s].mean(axis=0)
... plt.text(p[0], p[1], '#%d' % j, ha='center') # 标注三角形
>>> plt.xlim(-0.5, 1.5); plt.ylim(-0.5, 1.5)
>>> plt.show()
三角剖分的结构以以下方式编码：simplices 属性包含点在
构成三角形的``points``数组。例如：
>>> i = 1
>>> tri.simplices[i,:]
array([3, 1, 0], dtype=int32)
>>> points[tri.simplices[i,:]]
array([[ 1. , 1. ],
[ 0. , 1.1],
[ 0. , 0. ]])
此外，还可以找到相邻的三角形：
>>> tri.neighbors[i]
array([-1, 0, -1], dtype=int32)
这告诉我们，这个三角形有一个相邻的三角形 #0，但没有其他相邻三角形。此外，它告诉我们相邻的三角形 0 位于当前三角形的顶点 1 的对面：
>>> points[tri.simplices[i, 1]]
array([ 0. , 1.1])
确实，从图中我们可以看到情况确实如此。
Qhull 还可以对高维点集进行单纯形剖分（例如，在 3-D 中分割成四面体）。
共面点#
需要注意的是，并非所有点都必然出现在三角剖分的顶点中，这是由于在形成三角剖分时数值精度问题所致。考虑上述情况，但增加一个重复点：
>>> points = np.array([[0, 0], [0, 1], [1, 0], [1, 1], [1, 1]])
>>> tri = Delaunay(points)
>>> np.unique(tri.simplices.ravel())
array([0, 1, 2, 3], dtype=int32)
观察到点 #4（重复点）没有出现在三角剖分的顶点中。这一情况已被记录：
>>> tri.coplanar
array([[4, 0, 3]], dtype=int32)
这意味着点 4 位于三角形 0 和顶点 3 附近，但未包含在三角剖分中。
请注意，这种退化不仅可能因为重复点而发生，还可能因为更复杂的几何原因，甚至在乍看之下表现良好的点集中也会发生。
然而，Qhull 有一个 “QJ” 选项，它会随机扰动输入数据，直到退化问题得到解决：
>>> tri = Delaunay(points, qhull_options="QJ Pp")
>>> points[tri.simplices]
array([[[1, 0],
[1, 1],
[0, 0]],
[[1, 1],
[1, 1],
[1, 0]],
[[1, 1],
[0, 1],
[0, 0]],
[[0, 1],
[1, 1],
[1, 0]]])
[1, 1]]])
出现了两个新的三角形。然而，我们看到它们是退化的，并且面积为零。
凸包#
凸包是包含给定点集中所有点的最小凸对象。
这些可以通过 scipy.spatial 中的 Qhull 包装器计算如下：
>>> from scipy.spatial import ConvexHull
>>> rng = np.random.default_rng()
>>> points = rng.random((30, 2)) # 30 个二维随机点
>>> hull = ConvexHull(points)
凸包表示为一组 N 个一维单纯形，在二维中意味着线段。存储方案与上述讨论的 Delaunay 三角剖分的单纯形完全相同。
我们可以展示上述结果：
>>> import matplotlib.pyplot as plt
>>> plt.plot(points[:,0], points[:,1], 'o')
>>> for simplex in hull.simplices:
... plt.plot(points[simplex,0], points[simplex,1], 'k-')
>>> plt.show()
同样的结果可以通过 scipy.spatial.convex_hull_plot_2d 实现。
Voronoi 图#
Voronoi 图是将空间划分为给定点集的最近邻域的细分。
使用 scipy.spatial 有两种方法可以处理这个对象。首先，可以使用 KDTree 来回答“哪个点离这个点最近”的问题，并以此定义区域：
>>> from scipy.spatial import KDTree
>>> points = np.array([[0, 0], [0, 1], [0, 2], [1, 0], [1, 1], [1, 2],
... [2, 0], [2, 1], [2, 2]])
>>> tree = KDTree(points)
>>> tree.query([0.1, 0.1])
(0.14142135623730953, 0)
因此，点 (0.1, 0.1) 属于区域 0。用颜色表示：
>>> x = np.linspace(-0.5, 2.5, 31)
>>> y = np.linspace(-0.5, 2.5, 33)
>>> xx, yy = np.meshgrid(x, y)
>>> xy = np.c_[xx.ravel(), yy.ravel()]
>>> import matplotlib.pyplot as plt
>>> dx_half, dy_half = np.diff(x[:2])[0] / 2., np.diff(y[:2])[0] / 2.
>>> x_edges = np.concatenate((x - dx_half, [x[-1] + dx_half]))
>>> y_edges = np.concatenate((y - dy_half, [y[-1] + dy_half]))
>>> plt.pcolormesh(x_edges, y_edges, tree.query(xy)[1].reshape(33, 31), shading='flat')
>>> plt.plot(points[:,0], points[:,1], 'ko')
>>> plt.show()
然而，这并没有给出作为几何对象的Voronoi图。
通过`scipy.spatial`中的Qhull包装器，可以再次获得基于线和点的表示：
>>> from scipy.spatial import Voronoi
>>> vor = Voronoi(points)
>>> vor.vertices
array([[0.5, 0.5],
[0.5, 1.5],
[1.5, 0.5],
[1.5, 1.5]])
Voronoi顶点表示形成Voronoi区域多边形边缘的点集。在这种情况下，有9个不同的区域：
>>> vor.regions
[[], [-1, 0], [-1, 1], [1, -1, 0], [3, -1, 2], [-1, 3], [-1, 2], [0, 1, 3, 2], [2, -1, 0], [3, -1, 1]]
负值``-1``再次表示一个无穷远点。实际上，只有区域``[0, 1, 3, 2]``是有界的。请注意，由于与上述Delaunay三角剖分中类似的数值精度问题，Voronoi区域的数目可能少于输入点的数目。
分隔区域的脊线（二维中的线）被描述为与凸包片段类似的单纯形集合：
>>> vor.ridge_vertices
[[-1, 0], [-1, 0], [-1, 1], [-1, 1], [0, 1], [-1, 3], [-1, 2], [2, 3], [-1, 3], [-1, 2], [1, 3], [0, 2]]
这些数字表示构成线段的Voronoi顶点的索引。``-1``再次表示一个无穷远点——在12条线中，只有4条是有界的线段，而其他线段延伸至无穷远。
Voronoi脊线垂直于连接输入点的直线。 输入点。每条脊线对应的两点也被记录下来：
>>> vor.ridge_points
array([[0, 3],
[0, 1],
[2, 5],
[2, 1],
[1, 4],
[7, 8],
[7, 6],
[7, 4],
[8, 5],
[6, 3],
[4, 5],
[4, 3]], dtype=int32)
这些信息结合起来，足以构建完整的图表。
我们可以按如下方式绘制它。首先，绘制点和Voronoi顶点：
>>> plt.plot(points[:, 0], points[:, 1], 'o')
>>> plt.plot(vor.vertices[:, 0], vor.vertices[:, 1], '*')
>>> plt.xlim(-1, 3); plt.ylim(-1, 3)
绘制有限线段的方法与凸包相同，但现在我们必须处理无限边：
>>> for simplex in vor.ridge_vertices:
... simplex = np.asarray(simplex)
... if np.all(simplex >= 0):
... plt.plot(vor.vertices[simplex, 0], vor.vertices[simplex, 1], 'k-')
延伸到无穷远的脊线需要更多关注：
>>> center = points.mean(axis=0)
>>> for pointidx, simplex in zip(vor.ridge_points, vor.ridge_vertices):
... simplex = np.asarray(simplex)
... if np.any(simplex < 0):
... i = simplex[simplex >= 0][0] # 有限端Voronoi顶点
... t = points[pointidx[1]] - points[pointidx[0]] # 切线
... t = t / np.linalg.norm(t)
... n = np.array([-t[1], t[0]]) # 法线
... midpoint = points[pointidx].mean(axis=0)
... far_point = vor.vertices[i] + np.sign(np.dot(midpoint - center, n)) * n * 100
... plt.plot([vor.vertices[i,0], far_point[0]],
... [vor.vertices[i,1], far_point[1]], 'k--')
>>> plt.show()
此图也可以使用`scipy.spatial.voronoi_plot_2d`创建。
Voronoi图可以用来创建有趣的生成艺术。尝试调整这个``mandala``函数的设置，创造你自己的作品！
>>> import numpy as np
>>> from scipy import spatial
>>> import matplotlib.pyplot as plt
>>> def mandala(n_iter, n_points, radius):
... """使用Voronoi细分创建曼陀罗图形。
...
... 参数
... ----------
... n_iter : int
... 迭代次数，即生成等距点的次数。
... n_points : int
... 每次迭代绘制的点数。
... radius : scalar
... 径向扩展因子。
...
... 返回
... -------
... fig : matplotlib.Figure实例
...
... 注释
... -----
... 此代码改编自Audrey Roy Greenfeld [1]_ 和Carlos Focil-Espinosa [2]_ 的工作，他们使用Python代码创建了美丽的曼陀罗。
... 该代码基于Antonio Sánchez Chinchón的R代码 [3]_。
...
... 参考文献
... ----------
... .. [1] https://www.codemakesmehappy.com/2019/09/voronoi-mandalas.html
...
... .. [2] https://github.com/CarlosFocil/mandalapy
...
... .. [3] https://github.com/aschinchon/mandalas
...
... """
... fig = plt.figure(figsize=(10, 10))
... ax = fig.add_subplot(111)
... ax.set_axis_off()
... ax.set_aspect('equal', adjustable='box')
...
... angles = np.linspace(0, 2*np.pi * (1 - 1/n_points), num=n_points) + np.pi/2
... # 从单个中心点开始，迭代添加点
... xy = np.array([[0, 0]])
... for k in range(n_iter):
... t1 = np.array([])
... t2 = np.array([])
... # 在每次迭代中，围绕每个现有点添加`n_points`个新点
... for i in range(xy.shape[0]):
... t1 = np.append(t1, xy[i, 0] + radius**k * np.cos(angles))
... t2 = np.append(t2, xy[i, 1] + radius**k * np.sin(angles))
...
… xy = np.column_stack((t1, t2)) … … # 通过Voronoi图创建曼陀罗图形 … spatial.voronoi_plot_2d(spatial.Voronoi(xy), ax=ax) … … return fig
>>> # 修改以下参数以获得不同的图形
>>> n_iter = 3
>>> n_points = 6
>>> radius = 4
>>> fig = mandala(n_iter, n_points, radius)
>>> plt.show()
