SciPy - 空间算法

SciPy - 空间数据 (scipy.spatial)

scipy.spatial 模块提供了处理空间数据结构和几何计算的算法。它利用强大的 Qhull 库来完成诸如 Delaunay 三角剖分（Delaunay triangulations）、Voronoi 图（Voronoi diagrams）和凸包（convex hulls）等任务。此外，它还包括用于高效最近邻搜索（nearest-neighbor searches）的 KDTree 实现以及各种距离度量（distance metrics）的实用工具。

Delaunay 三角剖分

什么是 Delaunay 三角剖分？

在计算几何中，点集 P 的 Delaunay 三角剖分 DT(P) 是一种三角剖分，使得 P 中没有点位于 DT(P) 中任何三角形的外接圆（circumcircle）内部。此属性倾向于最大化三角剖分中所有三角形的最小角，从而避免出现“瘦长”的三角形。它们广泛应用于网格生成（mesh generation）、地理信息系统（geographic information systems）和表面重建（surface reconstruction）。

我们可以使用 scipy.spatial.Delaunay 计算 Delaunay 三角剖分。

import numpy as np
from scipy.spatial import Delaunay
import matplotlib.pyplot as plt

# Define a set of 2D points
points = np.array([
    [0, 4],
    [2, 1.1],
    [1, 3],
    [1, 2],
    [3, 3],
    [0.5, 0.5]
])

# Compute the Delaunay triangulation
tri = Delaunay(points)

# The 'simplices' attribute gives the indices of points forming each triangle
print("Simplices (triangles defined by point indices):")
print(tri.simplices)

# Plotting the triangulation
plt.figure(figsize=(7, 5))
plt.triplot(points[:, 0], points[:, 1], tri.simplices.copy(), color='blue')
plt.plot(points[:, 0], points[:, 1], 'o', color='red', markersize=5)
plt.title('Delaunay Triangulation')
plt.xlabel('X-coordinate')
plt.ylabel('Y-coordinate')
plt.grid(True)
# To display the plot (if not in an interactive environment like Jupyter):
# plt.show()
# For this tutorial, we'll describe the plot:
# This code will generate a plot displaying the red points and blue lines
# forming the Delaunay triangulation that connects these points.

tri.simplices 的输出将是一个数组，其中每行包含三个整数，表示在 points 数组中构成一个三角形的点的索引。该图将显示这些点通过线连接形成三角剖分。

三角剖分中的共面点

识别共面点

处理 3D Delaunay 三角剖分（或更高维度）时，某些点可能与现有单形（simplices）共面（或共超平面，co-hyperplanar），但它们本身并非三角剖分的一部分。Delaunay 对象的 coplanar 属性有助于识别这些点。它是一个数组 [point_idx, facet_idx, furthest_facet_idx]，指示 point_idx 与 facet_idx 共面，而 furthest_facet_idx 是沿该面的法线方向距离该点最远的面（facet）。

对于 2D 三角剖分，共面点通常意味着，如果它们位于现有边上，则是共线点，或者这些点会创建退化（扁平）单形。

import numpy as np
from scipy.spatial import Delaunay

# Points, including a point (1,1) that is part of a line segment if others exist
# And one point (index 4) intentionally made redundant with point (index 3)
points_coplanar_test = np.array([
    [0, 0],  # 0
    [0, 1],  # 1
    [1, 0],  # 2
    [1, 1],  # 3
    [1, 1]   # 4 (duplicate of point 3)
])

# Qhull (used by Delaunay) typically handles duplicate points by ignoring them or specific options.
# Let's use the 'QJ' option to jog points slightly if they are degenerate.
# Without 'QJ', Qhull might error or merge duplicates.
# For simplicity, let's make point 4 slightly off to see coplanar behavior.
points_coplanar_test_slight_diff = np.array([
    [0, 0],        # 0
    [2, 0],        # 1
    [1, 1.5],      # 2 (apex)
    [1, 0.00001]   # 3 (very close to the base 0-1)
])

tri_coplanar = Delaunay(points_coplanar_test_slight_diff)
print("Coplanar points/facets info (point_idx, nearest_facet_idx, furthest_facet_idx):")
print(tri_coplanar.coplanar)
# The output structure is (coplanar_point_index, simplex_index, vertex_index_within_simplex)
# or related to facets depending on Qhull version and context.

tri.coplanar 的输出将是一个数组，指示哪些点与三角剖分中的哪些面（2D 中是三角形，3D 中是四面体）共面。例如，array([[3, 0, 1]], dtype=int32) 可能意味着索引为 3 的点与面 0 共面，并且该面的顶点 1 是相关的（解释可能取决于 Qhull 的具体细节）。确切的输出可能因点配置和 Qhull 选项而异。

实际上，重复的点通常会预先过滤。coplanar 属性在 3D 及以上维度中更重要，因为一个点可能位于四面体的一个面上，但不是其顶点之一。

凸包

什么是凸包？

点集 X 的凸包（convex hull）是包含 X 的最小凸集。想象一下用一根橡皮筋拉伸围绕所有点；橡皮筋形成的形状就是凸包。它是计算几何中的一个基本概念，在碰撞检测、图像处理和模式识别中有应用。

scipy.spatial.ConvexHull 计算凸包。

import numpy as np
from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt

# Generate some random 2D points
np.random.seed(42) # for reproducibility
points_hull = np.random.rand(30, 2) # 30 random points in 2D

# Compute the convex hull
hull = ConvexHull(points_hull)

print("Vertices forming the convex hull (indices):")
print(hull.vertices) # Indices of points forming the hull, in order
print("\nEdges (simplices) of the convex hull:")
print(hull.simplices) # Pairs of vertex indices forming the edges of the hull

# Plotting the convex hull
plt.figure(figsize=(7, 5))
plt.plot(points_hull[:, 0], points_hull[:, 1], 'o', label='All Points', markersize=4)
# Plot the hull by connecting vertices in order
for simplex in hull.simplices:
    plt.plot(points_hull[simplex, 0], points_hull[simplex, 1], 'k-')
# Highlight hull vertices
plt.plot(points_hull[hull.vertices,0], points_hull[hull.vertices,1], 'ro', markersize=6, label='Hull Vertices')

plt.title('Convex Hull of Random Points')
plt.xlabel('X-coordinate')
plt.ylabel('Y-coordinate')
plt.legend()
plt.grid(True)
# plt.show()
# This code will display a scatter plot of the 30 random points.
# The convex hull will be drawn as a black polygon enclosing the outermost points.
# The vertices of this polygon will be highlighted with red circles.

hull.vertices 属性给出构成凸包顶点的点索引，按逆时针顺序排列。hull.simplices 提供凸包的边。

其他空间工具

scipy.spatial 还包括：

• Voronoi: 计算 Voronoi 图。
• KDTree 和 cKDTree: 用于快速最近邻查找。
• distance: 一个模块，包含计算点或点集之间各种距离度量（欧几里得距离 Euclidean, 曼哈顿距离 Manhattan 等）的函数。

这些工具使得 scipy.spatial 成为一个强大的模块，用于几何分析和点云数据处理。有关高级用法和更多示例，请参阅 SciPy Spatial 文档。