Spatial

from scipy.spatial import (
    KDTree, cKDTree,
    Voronoi, ConvexHull, Delaunay, QhullError,
    distance_matrix,
    SphericalVoronoi,
    Rectangle
)
from scipy.spatial.distance import (
    cdist, pdist, squareform,
    euclidean, cityblock, cosine, correlation,
    minkowski, chebyshev, hamming, jaccard, dice, rogerstanimoto,
    directed_hausdorff
)
import numpy as np

from scipy.spatial import KDTree

points = np.random.rand(1000, 3)   # N points in 3-D

tree = KDTree(points)
tree = KDTree(points, leafsize=10)  # leafsize controls tree depth (default 10)
tree = KDTree(points, compact_nodes=True, balanced_tree=True)

# cKDTree: same API but C-implemented (historically faster; now KDTree is C too)
from scipy.spatial import cKDTree
tree = cKDTree(points)   # identical API

# Single query point
query_pt = np.array([0.5, 0.5, 0.5])

# k nearest neighbors
dist, idx = tree.query(query_pt, k=1)      # 1 nearest
dist, idx = tree.query(query_pt, k=5)      # 5 nearest; dist/idx shape (5,)

# Multiple query points
query_pts = np.random.rand(100, 3)
dist, idx = tree.query(query_pts, k=3)     # dist/idx shape (100, 3)

# Options
dist, idx = tree.query(query_pts, k=5,
                        p=2,               # Minkowski p-norm (2=Euclidean, 1=Manhattan, inf=Chebyshev)
                        eps=0.0,           # approximate: at most (1+eps)*true_dist
                        distance_upper_bound=0.1,   # return only within this distance
                        workers=-1)        # parallel with all CPUs (SciPy 1.7+)

# Handling missing neighbors (when distance_upper_bound set)
# Missing: idx = tree.n (out of bounds), dist = inf
valid = dist < np.inf

# All neighbors within radius r
results = tree.query_ball_point(query_pt, r=0.1)   # returns list of indices
results = tree.query_ball_point(query_pts, r=0.1)  # list of lists

# Options
results = tree.query_ball_point(query_pt, r=0.1,
                                 p=2,              # p-norm
                                 eps=0.0,
                                 return_sorted=True,
                                 return_length=False)  # return count only

# Count without returning indices
count = tree.query_ball_point(query_pt, r=0.1, return_length=True)

# Pairs within distance r between two trees
tree2 = KDTree(np.random.rand(500, 3))
pairs = tree.query_ball_tree(tree2, r=0.05)   # pairs[i] = list of j in tree2

# Count neighbors within r for each query point
counts = tree.count_neighbors(tree2, r=0.1)          # single r
counts = tree.count_neighbors(tree2, r=[0.05, 0.1])  # multiple r (cumulative)

# Pairs within r in same tree
sparse_dist = tree.sparse_distance_matrix(tree2, max_distance=0.1)
# Returns dok_matrix; .toarray() for dense

# All-pairs distances within r
from scipy.spatial import cKDTree
pairs = list(tree.query_pairs(r=0.05))   # set of (i,j) pairs with i<j

from scipy.spatial.distance import cdist

A = np.random.rand(50, 3)
B = np.random.rand(30, 3)

D = cdist(A, B, metric='euclidean')   # shape (50, 30)
D = cdist(A, B, metric='minkowski', p=3)
D = cdist(A, B, metric='cosine')
D = cdist(A, B, metric='correlation')
D = cdist(A, B, metric='cityblock')   # L1 / Manhattan
D = cdist(A, B, metric='chebyshev')   # L-inf
D = cdist(A, B, metric='hamming')     # for binary / integer vectors
D = cdist(A, B, metric='jaccard')     # set distance
D = cdist(A, B, metric='mahalanobis', VI=np.linalg.inv(np.cov(A.T)))

# Custom metric
D = cdist(A, B, metric=lambda u, v: np.sum(np.abs(u - v)**3)**(1/3))

from scipy.spatial.distance import pdist, squareform

A = np.random.rand(50, 3)

# Condensed form (upper triangle, length n*(n-1)/2)
d = pdist(A)                        # Euclidean
d = pdist(A, metric='cosine')
d = pdist(A, metric='minkowski', p=1.5)

# Convert to/from square matrix
D_sq = squareform(d)      # condensed → (n×n) square matrix
d_back = squareform(D_sq) # square → condensed

from scipy.spatial.distance import euclidean, cityblock, cosine

euclidean([1, 0, 0], [0, 1, 0])   # sqrt(2)
cityblock([1, 2, 3], [4, 5, 6])   # 9
cosine([1, 0], [0, 1])            # 1.0 (orthogonal)

# Directed Hausdorff distance
from scipy.spatial.distance import directed_hausdorff

u = np.array([[0, 0], [1, 0], [1, 1]])
v = np.array([[0, 0], [0.5, 0.5]])
d, idx_u, idx_v = directed_hausdorff(u, v)

# Symmetric Hausdorff
haus = max(directed_hausdorff(u, v)[0], directed_hausdorff(v, u)[0])

from scipy.spatial import Voronoi, voronoi_plot_2d

points = np.random.rand(20, 2)
vor = Voronoi(points)

vor.points           # input points
vor.vertices         # Voronoi vertices
vor.ridge_points     # pairs of input points for each ridge
vor.ridge_vertices   # Voronoi vertex indices for each ridge (-1 = at infinity)
vor.regions          # list of vertex index lists per region (-1 = unbounded)
vor.point_region     # region index for each input point

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
voronoi_plot_2d(vor, ax=ax, show_vertices=True, line_colors='black')

# Incremental (add points later)
vor = Voronoi(points, incremental=True)
vor.add_points(new_points)

from scipy.spatial import Delaunay

points = np.random.rand(50, 2)
tri = Delaunay(points)

tri.simplices         # (nsimplex, ndim+1) array of vertex indices
tri.neighbors         # for each simplex: neighbor simplex indices (-1 = boundary)
tri.points            # input points
tri.convex_hull       # indices of boundary simplices

# Point location (find which triangle a point is in)
query = np.array([[0.5, 0.5], [0.1, 0.9]])
simplices = tri.find_simplex(query)    # -1 if outside hull

# Barycentric coordinates
b = tri.transform[simplices, :2] @ (query - tri.transform[simplices, 2]).T
# (more useful via tri.transform directly)

# Coplanarity check
in_hull = simplices >= 0

# Incremental
tri = Delaunay(points, incremental=True)
tri.add_points(new_points)

# Convex hull (as by-product)
from scipy.spatial import ConvexHull
hull = ConvexHull(points)
hull_vertices = hull.vertices          # indices of hull vertices (sorted CCW in 2D)
hull_simplices = hull.simplices        # facets (edges in 2D, triangles in 3D)
hull.volume                            # area in 2D, volume in 3D
hull.area                              # perimeter in 2D, surface area in 3D
hull.equations                         # facet hyperplane equations [normal, offset]

from scipy.spatial import ConvexHull

pts = np.random.rand(100, 3)  # 3-D points
hull = ConvexHull(pts)

hull.vertices   # indices of hull vertices
hull.simplices  # triangular facets (indices into pts)
hull.volume     # volume
hull.area       # surface area
hull.equations  # (nfacet, 4) — [A, b] s.t. A @ x + b <= 0 inside

# Point-in-hull test (Delaunay method)
tri = Delaunay(pts)
inside = tri.find_simplex(test_pts) >= 0

# Incremental hull
hull = ConvexHull(pts, incremental=True)
hull.add_points(new_pts)

from scipy.spatial import convex_hull_plot_2d
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
convex_hull_plot_2d(ConvexHull(pts[:, :2]), ax=ax)

from scipy.spatial import SphericalVoronoi

# Points on the unit sphere
n = 100
theta = np.random.uniform(0, np.pi, n)
phi = np.random.uniform(0, 2*np.pi, n)
pts = np.column_stack([np.sin(theta)*np.cos(phi),
                        np.sin(theta)*np.sin(phi),
                        np.cos(theta)])

sv = SphericalVoronoi(pts, radius=1.0, center=np.array([0, 0, 0]))
sv.sort_vertices_of_regions()     # sort region vertices in CCW order

sv.vertices              # coordinates of Voronoi vertices on sphere
sv.regions               # list of vertex indices per region
sv.calculate_areas()     # area of each Voronoi cell (sum ≈ 4π)

from scipy.spatial import distance_matrix

A = np.random.rand(10, 3)
B = np.random.rand(8, 3)

D = distance_matrix(A, B)    # shape (10, 8), L2 by default
D = distance_matrix(A, B, p=1)   # L1 norm
D = distance_matrix(A, B, p=np.inf)   # Chebyshev