Spatial algorithms: KDTree, distances and Delaunay

Answer nearest-neighbour queries quickly, compute distance matrices without loops, triangulate point clouds, and align two sets of vectors.

Nearest neighbours with cKDTree

import numpy as np
from scipy.spatial import cKDTree, distance

rng = np.random.default_rng(0)
points = rng.random((20_000, 3))
queries = rng.random((500, 3))

tree = cKDTree(points)                       # build once, query many times

d, idx = tree.query(queries, k=5)            # 5 nearest per query
print(d.shape, idx.shape)                    # (500, 5) for both

within = tree.query_ball_point(queries[0], r=0.1)
counts = tree.query_ball_tree(cKDTree(queries), r=0.05)

# pairs closer than a threshold without building a full matrix
pairs = tree.query_pairs(r=0.02, output_type="ndarray")
print(within[:5], len(pairs))

• Building the tree costs O(n log n) and pays off as soon as you run more than a handful of queries.
• query with k=1 is the common case; workers=-1 uses all cores and can halve the wall time.
• Set distance_upper_bound to skip distant matches and speed up large queries.
• For a periodic domain, use the minimum-image convention on the distances or replicating points; a plain KDTree does not wrap around edges.

Distance matrices

a = rng.random((1000, 4))
b = rng.random((800, 4))

compact = distance.pdist(a, metric="euclidean")      # 1D, n*(n-1)/2 entries
square = distance.squareform(compact)                # (1000, 1000) dense
cross = distance.cdist(a, b, metric="cosine")        # (1000, 800)

print(compact.shape, square.shape, cross.shape)

Triangulations and alignment

from scipy.spatial import Delaunay, Voronoi, ConvexHull, Rotation, procrustes

pts = np.random.default_rng(1).random((60, 2))

tri = Delaunay(pts)
print(tri.simplices.shape)                   # triangles as index triples
mask = tri.find_simplex(np.array([[0.3, 0.4], [0.9, 0.1]]))
print(mask)

hull = ConvexHull(pts)
print(hull.volume, hull.area if hasattr(hull, "area") else hull.area)

# rigid alignment of two sets of vectors (Kabsch algorithm)
A = np.random.default_rng(2).random((10, 3))
B = A @ np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]], dtype=float) + 0.1
rot, rssd = Rotation.align_vectors(A, B)
print(rot.as_matrix().round(3), rssd)

result = procrustes(A, B)                    # rotation, scale and translation
print(result.scale)

These are the primitives behind collision detection, mesh generation, geographic clustering and molecular alignment. Each takes a point array of shape (n_points, n_dimensions) and returns indices into it, so your original data stays the source of truth.

FAQ

Related

Image processing with scipy.ndimage Linear algebra with scipy.linalg

Last refreshed 2026-09-18.