SciPy - Spatial

SciPy - Spatial Data (scipy.spatial)

The scipy.spatial module provides algorithms for working with spatial data structures and geometric computations. It leverages the robust Qhull library for tasks like Delaunay triangulations, Voronoi diagrams, and convex hulls. Additionally, it includes KDTree implementations for efficient nearest-neighbor searches and utilities for various distance metrics.

Delaunay Triangulations

What are Delaunay Triangulations?

In computational geometry, a Delaunay triangulation for a set of points P in a plane is a triangulation DT(P) such that no point in P lies inside the circumcircle of any triangle in DT(P). This property tends to maximize the minimum angle of all triangles in the triangulation, avoiding ‘skinny’ triangles. They are widely used in mesh generation, geographic information systems, and surface reconstruction.

We can compute Delaunay triangulations using scipy.spatial.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.

The output for tri.simplices would be an array where each row contains three integers, representing the indices of the points in the points array that form a triangle. The plot would show these points connected by lines to form the triangulation.

Coplanar Points in Triangulation

Identifying Coplanar Points

When dealing with 3D Delaunay triangulations (or higher dimensions), some points might be coplanar (or co-hyperplanar) with existing simplices but not part of the triangulation itself. The coplanar attribute of the Delaunay object helps identify these. It’s an array of [point_idx, facet_idx, furthest_facet_idx], indicating that point_idx is coplanar with facet_idx, and furthest_facet_idx is the facet furthest from this point along the facet’s normal.

For 2D triangulations, coplanar points usually mean collinear points if they lie on an existing edge, or points that would create a degenerate (flat) simplex.

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.

The output of tri.coplanar will be an array indicating which points are coplanar with which facets (triangles in 2D, tetrahedra in 3D) of the triangulation. For example, array([[3, 0, 1]], dtype=int32) could mean point with index 3 is coplanar with facet 0, and vertex 1 of that facet is relevant (interpretation can depend on Qhull details). The exact output can vary based on point configuration and Qhull options.

In practice, duplicate points are often pre-filtered. The coplanar attribute is more significant in 3D+ where a point might lie on the face of a tetrahedron but not be one of its vertices.

Convex Hulls

What are Convex Hulls?

The convex hull of a set of points X is the smallest convex set that contains X. Imagine stretching a rubber band around all the points; the shape formed by the rubber band is the convex hull. It’s a fundamental concept in computational geometry with applications in collision detection, image processing, and pattern recognition.

scipy.spatial.ConvexHull computes the convex hull.

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.

The hull.vertices attribute gives the indices of the points that form the vertices of the convex hull, ordered counter-clockwise. hull.simplices provides the edges of the hull.

Other Spatial Utilities

scipy.spatial also includes:

• Voronoi: Computes Voronoi diagrams.
• KDTree and cKDTree: For fast nearest-neighbor lookups.
• distance: A module with functions to compute various distance metrics (Euclidean, Manhattan, etc.) between points or sets of points.

These tools make scipy.spatial a powerful module for geometric analysis and processing of point cloud data. For advanced usage and more examples, refer to the SciPy Spatial documentation.