mathematicskit.geometry#
Convex hull via scipy.spatial.ConvexHull, with a hand-rolled Graham
scan pedagogical comparison in 2D; Delaunay triangulation and its
Voronoi dual via scipy.spatial; hand-rolled segment intersection and
point-in-polygon tests; polygon area/centroid via the shoelace formula;
and curvature, arc length, and the Frenet-Serret frame for parametric
curves.
mathematicskit.geometry: computational geometry, built directly on scipy.spatial.
Convex hull via scipy.spatial.ConvexHull, with a hand-rolled Graham
scan kept as a pedagogical comparison in 2D; Delaunay triangulation via
scipy.spatial.Delaunay and its dual Voronoi diagram via
scipy.spatial.Voronoi; line/segment intersection and point-in-
polygon tests (hand-rolled – no direct scipy equivalent); polygon area
and centroid via the shoelace formula; and curvature, arc length, and
the Frenet-Serret frame for parametric plane/space curves (derivatives
via numpy.gradient, arc length via scipy.integrate, the frame
itself assembled by mathematicskit’s own code); Heron’s formula and
Pick’s lattice-point counts; Euler’s polyhedron formula on 3D hulls;
Gaussian and mean curvature of parametric surfaces; the smallest
enclosing circle (Welzl); the Hausdorff distance; Bézier curves via de
Casteljau; Douglas-Peucker line simplification; and the closest pair of
points by divide and conquer.
- class mathematicskit.geometry.CircleResult(center, radius)[source]#
Bases:
objectContainer for a circle in the plane.
- class mathematicskit.geometry.ClosestPairResult(indices, distance)[source]#
Bases:
objectContainer for the closest pair of points in a point set.
- class mathematicskit.geometry.ConvexHullResult(points, vertices, simplices, volume, area, method='')[source]#
Bases:
objectContainer for a convex hull.
- Parameters:
- class mathematicskit.geometry.CurveFrameResult(t, position, tangent, normal, binormal=None, curvature=<factory>, torsion=None, arc_length=<factory>)[source]#
Bases:
objectContainer for a parametric curve’s Frenet-Serret frame.
- Parameters:
- class mathematicskit.geometry.LatticePolygonResult(interior, boundary, area)[source]#
Bases:
objectContainer for the lattice-point counts of a polygon with integer vertices.
- class mathematicskit.geometry.PolyhedronResult(vertices, edges, faces)[source]#
Bases:
objectContainer for the vertex, edge, and face counts of a convex polyhedron.
- class mathematicskit.geometry.SurfaceCurvatureResult(u, v, points, gaussian, mean, area_element)[source]#
Bases:
objectContainer for the curvatures of a parametric surface \(\mathbf{r}(u, v)\) on a grid.
- Parameters:
- class mathematicskit.geometry.TriangulationResult(points, simplices)[source]#
Bases:
objectContainer for a Delaunay triangulation.
- class mathematicskit.geometry.VoronoiResult(points, vertices, regions, ridge_points)[source]#
Bases:
objectContainer for a Voronoi diagram (the Delaunay triangulation’s dual).
- regions: list#
Each input point’s Voronoi region, as indices into vertices (
-1marks an unbounded region).
- mathematicskit.geometry.bezier_curve(control_points, t)[source]#
Points \(\mathbf{B}(t) = \sum_i \binom{n}{i} t^i (1-t)^{n-i} \mathbf{p}_i\) on a Bézier curve.
- Parameters:
control_points (array_like, shape (n + 1, d))
t (array_like) – Parameter values in \([0, 1]\).
- Return type:
- Returns:
ndarray, shape (len(t), d)
Examples
>>> bezier_curve([[0, 0], [1, 2], [2, 0]], [0.0, 0.5, 1.0]).tolist() [[0.0, 0.0], [1.0, 1.0], [2.0, 0.0]]
- mathematicskit.geometry.circle(radius=1.0)[source]#
A circle of the given radius in the xy-plane: \(\gamma(t) = (r\cos t, r\sin t)\).
Known curvature: \(\kappa = 1/r\) everywhere.
- Parameters:
radius (
float)- Returns:
callable –
curve(t) -> ndarrayof shape(len(t), 2).
- mathematicskit.geometry.closest_pair(points)[source]#
The two closest points in a planar point set, in \(O(n \log n)\) time.
Sorts by \(x\), splits at the median, solves both halves, and then checks only points within the current best distance \(\delta\) of the dividing line. Sorted by \(y\), each such point needs to be compared with at most seven of its successors.
- Parameters:
points (array_like, shape (n, 2)) – At least two points.
- Return type:
- Returns:
ClosestPairResult
Examples
>>> result = closest_pair([[0, 0], [5, 5], [1, 1], [9, 0], [5.5, 5]]) >>> result.indices, result.distance ((1, 4), 0.5)
- mathematicskit.geometry.convex_hull(points)[source]#
Convex hull of a point set, via
scipy.spatial.ConvexHull(Qhull).Works in any dimension; the 2D case returns hull vertices in counterclockwise order. See O’Rourke, Computational Geometry in C, 2nd ed., Ch. 4.
- Parameters:
points (
ndarray)- Return type:
- Returns:
ConvexHullResult
Examples
>>> import numpy as np >>> points = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.5, 0.5]]) >>> result = convex_hull(points) >>> sorted(int(v) for v in result.vertices) # the interior point (0.5, 0.5) is excluded [0, 1, 2, 3] >>> round(result.volume, 6) # area of the unit square 1.0
- mathematicskit.geometry.cylinder_surface(radius=1.0)[source]#
The cylinder \((R\cos v, R\sin v, u)\): curved in space, yet \(K = 0\), since it unrolls flat.
Examples
>>> np.round(cylinder_surface(1.0)(2.0, 0.0), 12).tolist() [1.0, 0.0, 2.0]
- mathematicskit.geometry.de_casteljau(control_points, t)[source]#
All intermediate point sets of de Casteljau’s algorithm at parameter
t.Each level replaces consecutive points \(p_i, p_{i+1}\) by \((1-t)p_i + t\,p_{i+1}\); the single point of the last level is the curve point \(\mathbf{B}(t)\).
- Parameters:
- Return type:
- Returns:
list of ndarray – Level
khas shape(n + 1 - k, d).
Examples
>>> levels = de_casteljau([[0, 0], [1, 2], [2, 0]], 0.5) >>> levels[-1].tolist() [[1.0, 1.0]]
- mathematicskit.geometry.delaunay_triangulation(points)[source]#
Delaunay triangulation of a point set, via
scipy.spatial.Delaunay(Qhull).The triangulation maximizing the minimum angle among all possible triangulations of the point set – equivalently, no point lies inside any triangle’s circumcircle. See de Berg et al., Computational Geometry, 3rd ed., Ch. 9.1-9.3.
- Parameters:
points (
ndarray)- Return type:
- Returns:
TriangulationResult
Examples
>>> import numpy as np >>> points = np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]) >>> result = delaunay_triangulation(points) >>> result.simplices.shape[1] # each simplex is a triangle: 3 vertices 3
- mathematicskit.geometry.douglas_peucker(points, epsilon)[source]#
Simplify a polyline so that no removed point lies farther than
epsilonfrom the result.Keeps the two endpoints, finds the interior point farthest from the segment joining them, and recurses on both halves if that distance exceeds
epsilon; otherwise drops every interior point.- Parameters:
- Return type:
- Returns:
ndarray, shape (m, 2) – The retained vertices, a subset of
pointsin their original order.
Examples
>>> douglas_peucker([[0, 0], [1, 0.1], [2, -0.1], [3, 5], [4, 6], [5, 7]], 1.0).tolist() [[0.0, 0.0], [2.0, -0.1], [3.0, 5.0], [5.0, 7.0]]
- mathematicskit.geometry.ellipse(a=2.0, b=1.0)[source]#
An ellipse in the xy-plane: \(\gamma(t) = (a\cos t, b\sin t)\).
Curvature varies with
t(maximal at the ends of the major axis, minimal at the ends of the minor axis) – a useful non-constant- curvature test case, unlikecircle().
- mathematicskit.geometry.frenet_serret_frame(curve, t)[source]#
The Frenet-Serret frame (tangent, normal, binormal), curvature, and arc length of a parametric curve.
Given \(\gamma(t)\), the unit tangent is \(T = \gamma'/\|\gamma'\|\), curvature \(\kappa = \|\gamma' \times \gamma''\|/\|\gamma'\|^3\) (3D) or the signed 2D analogue, the (principal) normal \(N = T'/\|T'\|\), and (3D only) binormal \(B = T \times N\) with torsion \(\tau\) from the rate of change of \(B\). First and second derivatives are estimated via
numpy.gradient(); arc length viascipy.integrate.cumulative_trapezoid()on \(\|\gamma'\|\). See do Carmo, Differential Geometry of Curves and Surfaces, 2nd ed., Ch. 1.3-1.5.- Parameters:
- Return type:
- Returns:
CurveFrameResult
Examples
>>> import numpy as np >>> # A circle of radius 2 has constant curvature 1/2. >>> circle = lambda t: np.column_stack([2.0 * np.cos(t), 2.0 * np.sin(t)]) >>> t = np.linspace(0.0, 2.0 * np.pi, 400, endpoint=False) >>> result = frenet_serret_frame(circle, t) >>> bool(np.allclose(result.curvature[5:-5], 0.5, atol=1e-2)) True
- mathematicskit.geometry.graham_scan(points)[source]#
Convex hull of a planar point set via the Graham scan.
Sorts points by angle around the lowest (then leftmost) point, and sweeps them, popping the hull-in-progress whenever the last three points make a clockwise (non-left) turn – kept hand-rolled specifically to compare against
convex_hull()’s Qhull-based result, not as the primary 2D convex-hull API. See de Berg et al., Computational Geometry, 3rd ed., Ch. 1.1.- Parameters:
points (
ndarray)- Return type:
- Returns:
ConvexHullResult –
verticesholds the hull’s point indices in counterclockwise order.simplicesis left empty, since the Graham scan produces only the ordered vertex cycle, not a facet list.volumeis the enclosed area via the shoelace formula on those vertices – matchingconvex_hull(), where Qhull likewise reports the enclosed area asvolumein 2D.areais left at0.0: Qhull uses it for the hull’s perimeter, which this scan does not compute.
Examples
>>> import numpy as np >>> points = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.5, 0.5]]) >>> result = graham_scan(points) >>> sorted(int(v) for v in result.vertices) [0, 1, 2, 3] >>> round(result.volume, 6) 1.0
- mathematicskit.geometry.hausdorff_distance(a, b)[source]#
The Hausdorff distance \(d_H(A, B) = \max\{\sup_{a}\inf_{b}\|a-b\|,\; \sup_{b}\inf_{a}\|a-b\|\}\).
The smallest \(r\) such that each set lies within distance \(r\) of the other. Unlike the distance between closest points, it is zero only when the (closed) sets coincide.
- Parameters:
- Return type:
- Returns:
float
Examples
>>> hausdorff_distance([[0, 0], [1, 0]], [[0, 0], [1, 0], [1, 3]]) 3.0
- mathematicskit.geometry.helix(radius=1.0, pitch=1.0)[source]#
A circular helix: \(\gamma(t) = (r\cos t, r\sin t, pt)\).
Known curvature \(\kappa = r/(r^2+p^2)\) and torsion \(\tau = p/(r^2+p^2)\), both constant.
- mathematicskit.geometry.heron_area(a, b, c)[source]#
Area of a triangle from its three side lengths, by Heron’s formula.
With semi-perimeter \(s = (a+b+c)/2\), \(A = \sqrt{s(s-a)(s-b)(s-c)}\). Evaluated in William Kahan’s numerically stable arrangement, which stays accurate for needle-like triangles where the textbook form loses digits: with \(a \ge b \ge c\),
\[A = \tfrac14\sqrt{(a+(b+c))(c-(a-b))(c+(a-b))(a+(b-c))}.\]See W. Kahan, “Miscalculating Area and Angles of a Needle-like Triangle” (lecture notes, University of California, Berkeley, 2014).
- Parameters:
- Return type:
- Returns:
float
Examples
>>> heron_area(3.0, 4.0, 5.0) 6.0
- mathematicskit.geometry.lattice_point_counts(vertices)[source]#
Count the lattice points inside and on the boundary of a simple polygon with integer vertices.
Boundary points are counted edge by edge with \(\gcd(|\Delta x|, |\Delta y|)\); interior points are counted directly by testing every lattice point in the bounding box with
point_in_polygon(). Pick’s theorem says the area equals \(I + B/2 - 1\); the result reports both sides so the theorem can be checked.- Parameters:
- Return type:
- Returns:
LatticePolygonResult
Examples
>>> result = lattice_point_counts([[0, 0], [4, 0], [4, 3], [0, 3]]) >>> result.interior, result.boundary, result.area, result.pick_area (6, 14, 12.0, 12.0)
- mathematicskit.geometry.min_enclosing_circle(points, seed=0)[source]#
The smallest circle containing every point, by Welzl’s move-to-front algorithm.
Processes the points in random order. Whenever a point falls outside the current circle, the optimal circle must pass through it, and the search restarts on the points seen so far with that point fixed on the boundary. The optimum is determined by at most three points.
- Parameters:
- Return type:
- Returns:
CircleResult
Examples
>>> result = min_enclosing_circle([[0, 0], [2, 0], [1, 0.5]]) >>> result.center.tolist(), result.radius ([1.0, 0.0], 1.0)
- mathematicskit.geometry.point_in_polygon(point, polygon)[source]#
Whether point lies inside polygon, via the ray-casting (crossing-number) algorithm.
Casts a ray from point in the \(+x\) direction and counts how many polygon edges it crosses; an odd count means the point is inside (the standard Jordan-curve-theorem argument). See O’Rourke, Computational Geometry in C, 2nd ed., Ch. 7.4.
- Parameters:
- Return type:
- Returns:
bool
Examples
>>> import numpy as np >>> square = np.array([[0.0, 0.0], [4.0, 0.0], [4.0, 4.0], [0.0, 4.0]]) >>> point_in_polygon(np.array([2.0, 2.0]), square) True >>> point_in_polygon(np.array([5.0, 2.0]), square) False
- mathematicskit.geometry.polygon_area(vertices)[source]#
Signed area of a simple polygon via the shoelace formula.
\(A = \tfrac12 \left|\sum_{i} (x_i y_{i+1} - x_{i+1} y_i)\right|\), indices modulo
n(the polygon implicitly closes back to its first vertex). Positive for counterclockwise-ordered vertices, negative for clockwise (the absolute value is returned here).- Parameters:
vertices (
ndarray) – Polygon vertices in order (either winding direction).- Return type:
- Returns:
float
Examples
>>> polygon_area(np.array([[0.0, 0.0], [4.0, 0.0], [4.0, 3.0], [0.0, 3.0]])) 12.0
- mathematicskit.geometry.polygon_centroid(vertices)[source]#
Centroid (center of mass, assuming uniform density) of a simple polygon.
\(C_x = \dfrac{1}{6A}\sum_i (x_i+x_{i+1})(x_iy_{i+1}-x_{i+1}y_i)\), similarly for \(C_y\) – the shoelace-formula-weighted average, not the plain average of the vertex coordinates (which is only correct for a regular polygon). See Bourke (1988).
Examples
>>> polygon_centroid(np.array([[0.0, 0.0], [4.0, 0.0], [4.0, 3.0], [0.0, 3.0]])) array([2. , 1.5])
- mathematicskit.geometry.polyhedron_counts(points, decimals=9)[source]#
Count the vertices, edges, and faces of the convex hull of 3D points.
Triangles whose supporting planes agree (to
decimalsplaces) are merged into one face; an edge is counted only where two different faces meet.- Parameters:
- Return type:
- Returns:
PolyhedronResult
Examples
>>> import itertools >>> cube = np.array(list(itertools.product((0, 1), repeat=3)), dtype=float) >>> result = polyhedron_counts(cube) >>> result.vertices, result.edges, result.faces, result.euler_characteristic (8, 12, 6, 2)
- mathematicskit.geometry.segment_intersection(p1, p2, p3, p4)[source]#
Intersection point of segments
p1-p2andp3-p4, if they cross.Writes both segments parametrically (\(p_1 + t(p_2-p_1)\), \(p_3 + u(p_4-p_3)\)) and solves the resulting \(2\times2\) linear system for \(t, u \in [0, 1]\) via Cramer’s rule; returns
Noneif the segments are parallel or don’t overlap within their endpoints. See de Berg et al., Computational Geometry, 3rd ed., Ch. 2.- Parameters:
- Return type:
- Returns:
ndarray, shape (2,), or None
Examples
>>> import numpy as np >>> p = segment_intersection(np.array([0.0, 0.0]), np.array([2.0, 2.0]), np.array([0.0, 2.0]), np.array([2.0, 0.0])) >>> np.allclose(p, [1.0, 1.0]) True >>> segment_intersection(np.array([0.0, 0.0]), np.array([1.0, 0.0]), np.array([0.0, 1.0]), np.array([1.0, 1.0])) is None True
- mathematicskit.geometry.sphere_surface(radius=1.0)[source]#
The sphere \(\mathbf{r}(u, v) = R(\sin u \cos v, \sin u \sin v, \cos u)\); \(K = 1/R^2\).
Examples
>>> np.round(sphere_surface(1.0)(np.pi / 2, 0.0), 12).tolist() [1.0, 0.0, 0.0]
- mathematicskit.geometry.surface_curvature(surface, u, v)[source]#
Gaussian curvature \(K\) and mean curvature \(H\) of \(\mathbf{r}(u, v)\) on a grid.
\[K = \frac{LN - M^2}{EG - F^2}, \qquad H = \frac{EN - 2FM + GL}{2(EG - F^2)},\]with \(E = \mathbf{r}_u\cdot\mathbf{r}_u\), \(F = \mathbf{r}_u\cdot\mathbf{r}_v\), \(G = \mathbf{r}_v\cdot\mathbf{r}_v\), and \(L, M, N\) the components of \(\mathbf{r}_{uu}, \mathbf{r}_{uv}, \mathbf{r}_{vv}\) along the unit normal. Derivatives are second-order accurate, including at the grid’s edges (
edge_order=2).- Parameters:
surface (
Callable) –surface(U, V)returns an array of shape(3, *U.shape).u (array_like) – Increasing, equally spaced 1D parameter grids.
v (array_like) – Increasing, equally spaced 1D parameter grids.
- Return type:
- Returns:
SurfaceCurvatureResult
Examples
>>> u, v = np.linspace(0.5, 2.5, 81), np.linspace(0, 2 * np.pi, 81) >>> result = surface_curvature(sphere_surface(2.0), u, v) >>> round(float(np.median(result.gaussian)), 4) # 1 / R^2 0.25
- mathematicskit.geometry.torus_surface(major=2.0, minor=1.0)[source]#
The torus \(((R + r\cos u)\cos v, (R + r\cos u)\sin v, r\sin u)\); \(K = \cos u / (r(R + r\cos u))\).
- Parameters:
- Return type:
- Returns:
callable
Examples
>>> np.round(torus_surface(2.0, 1.0)(0.0, 0.0), 12).tolist() [3.0, 0.0, 0.0]
- mathematicskit.geometry.voronoi_diagram(points)[source]#
Voronoi diagram of a point set, via
scipy.spatial.Voronoi(Qhull).The dual of the Delaunay triangulation: partitions the plane (or higher-dimensional space) into regions, one per input point, consisting of all locations closer to that point than to any other. Voronoi vertices are exactly the circumcenters of the dual Delaunay triangles. See de Berg et al., Computational Geometry, 3rd ed., Ch. 7.
- Parameters:
points (
ndarray)- Return type:
- Returns:
VoronoiResult
Examples
>>> import numpy as np >>> points = np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]) >>> result = voronoi_diagram(points) >>> len(result.regions) >= 4 True