beam: convex hull monotone chain with tests

wizard-work/beam/convex_hull.py

@@ -0,0 +1,101 @@
+"""Convex hull via Andrew's monotone chain.
+
+Returns the hull vertices in counter-clockwise order starting from the
+leftmost-lowest point. Collinear points on hull edges are excluded by
+default (keep_collinear=False); set True to keep them.
+
+Complexity: O(n log n).
+"""
+from __future__ import annotations
+from typing import Iterable, List, Sequence, Tuple
+
+Point = Tuple[float, float]
+
+
+def _cross(o: Point, a: Point, b: Point) -> float:
+    return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])
+
+
+def convex_hull(points: Iterable[Point], keep_collinear: bool = False) -> List[Point]:
+    pts = sorted(set(points))
+    if len(pts) <= 1:
+        return list(pts)
+
+    # Build lower hull
+    lower: List[Point] = []
+    for p in pts:
+        while len(lower) >= 2:
+            c = _cross(lower[-2], lower[-1], p)
+            if keep_collinear:
+                if c < 0:
+                    lower.pop()
+                else:
+                    break
+            else:
+                if c <= 0:
+                    lower.pop()
+                else:
+                    break
+        lower.append(p)
+
+    # Build upper hull
+    upper: List[Point] = []
+    for p in reversed(pts):
+        while len(upper) >= 2:
+            c = _cross(upper[-2], upper[-1], p)
+            if keep_collinear:
+                if c < 0:
+                    upper.pop()
+                else:
+                    break
+            else:
+                if c <= 0:
+                    upper.pop()
+                else:
+                    break
+        upper.append(p)
+
+    # Concatenate, removing last point of each (it's repeated)
+    hull = lower[:-1] + upper[:-1]
+    return hull
+
+
+def hull_area(hull: Sequence[Point]) -> float:
+    """Shoelace area (absolute)."""
+    n = len(hull)
+    if n < 3:
+        return 0.0
+    s = 0.0
+    for i in range(n):
+        x1, y1 = hull[i]
+        x2, y2 = hull[(i + 1) % n]
+        s += x1 * y2 - x2 * y1
+    return abs(s) / 2.0
+
+
+def point_in_hull(p: Point, hull: Sequence[Point]) -> bool:
+    """Ray casting point-in-polygon (boundary counts as inside)."""
+    n = len(hull)
+    if n == 0:
+        return False
+    if n == 1:
+        return p == hull[0]
+    if n == 2:
+        # segment
+        return _cross(hull[0], hull[1], p) == 0 and min(hull[0][0], hull[1][0]) <= p[0] <= max(hull[0][0], hull[1][0]) and min(hull[0][1], hull[1][1]) <= p[1] <= max(hull[0][1], hull[1][1])
+    # boundary check
+    for i in range(n):
+        if _cross(hull[i], hull[(i + 1) % n], p) == 0:
+            a, b = hull[i], hull[(i + 1) % n]
+            if min(a[0], b[0]) <= p[0] <= max(a[0], b[0]) and min(a[1], b[1]) <= p[1] <= max(a[1], b[1]):
+                return True
+    # ray casting
+    inside = False
+    j = n - 1
+    for i in range(n):
+        xi, yi = hull[i]
+        xj, yj = hull[j]
+        if ((yi > p[1]) != (yj > p[1])) and (p[0] < (xj - xi) * (p[1] - yi) / (yj - yi) + xi):
+            inside = not inside
+        j = i
+    return inside

wizard-work/beam/test_convex_hull.py

@@ -0,0 +1,44 @@
+import sys, os
+sys.path.insert(0, os.path.dirname(__file__))
+from convex_hull import convex_hull, hull_area, point_in_hull
+
+def test_square():
+    pts = [(0,0),(1,0),(1,1),(0,1),(0.5,0.5)]
+    h = convex_hull(pts)
+    assert len(h) == 4
+    assert hull_area(h) == 1.0
+
+def test_triangle():
+    h = convex_hull([(0,0),(2,0),(1,2),(1,1)])
+    assert len(h) == 3
+    assert abs(hull_area(h)-2.0) < 1e-9
+
+def test_collinear_excluded():
+    h = convex_hull([(0,0),(1,0),(2,0),(2,2),(0,2)])
+    # collinear bottom edge middle point excluded
+    assert (1,0) not in h
+    assert len(h) == 4
+
+def test_keep_collinear():
+    h = convex_hull([(0,0),(1,0),(2,0),(2,2),(0,2)], keep_collinear=True)
+    assert (1,0) in h
+
+def test_duplicates():
+    h = convex_hull([(0,0),(0,0),(1,0),(0,1)])
+    assert len(h) == 3
+
+def test_single():
+    assert convex_hull([(5,5)]) == [(5,5)]
+
+def test_empty():
+    assert convex_hull([]) == []
+
+def test_point_in():
+    h = [(0,0),(2,0),(2,2),(0,2)]
+    assert point_in_hull((1,1), h)
+    assert point_in_hull((0,0), h)  # boundary
+    assert not point_in_hull((3,3), h)
+
+def test_area_triangle():
+    h = [(0,0),(4,0),(0,3)]
+    assert abs(hull_area(h)-6.0) < 1e-9