diff --git a/src/freeroute/geometry/__init__.py b/src/freeroute/geometry/__init__.py index 30944c1..b0d45f7 100644 --- a/src/freeroute/geometry/__init__.py +++ b/src/freeroute/geometry/__init__.py @@ -22,6 +22,8 @@ from .float_point import FloatPoint from .limits import CRIT_DOUBLE, CRIT_INT from .line import Line from .point import IntPoint, Point, RationalPoint, point, rational_point +from .polygon import Polygon +from .polygon_shape import PolygonShape from .side import Side, Signum from .simplex import Simplex from .tile import TileShape @@ -47,4 +49,6 @@ __all__ = [ "IntBox", "TileShape", "Simplex", + "Polygon", + "PolygonShape", ] diff --git a/src/freeroute/geometry/line.py b/src/freeroute/geometry/line.py index 61777a9..41ae3bd 100644 --- a/src/freeroute/geometry/line.py +++ b/src/freeroute/geometry/line.py @@ -198,6 +198,28 @@ class Line: result = self.side_of(line_1.intersection(line_2)) return result + def function_value_approx(self, x: float) -> float: + """Approximate ``y`` on this (non-vertical) line at abscissa ``x``.""" + p1 = self.a.to_float() + p2 = self.b.to_float() + dx = p2.x - p1.x + if dx == 0: + return 0.0 + dy = p2.y - p1.y + det = p1.x * p2.y - p2.x * p1.y + return (dy * x - det) / dx + + def function_in_y_value_approx(self, y: float) -> float: + """Approximate ``x`` on this (non-horizontal) line at ordinate ``y``.""" + p1 = self.a.to_float() + p2 = self.b.to_float() + dy = p2.y - p1.y + if dy == 0: + return 0.0 + dx = p2.x - p1.x + det = p1.x * p2.y - p2.x * p1.y + return (dx * y + det) / dy + # --- intersection (exact) ---------------------------------------------- def intersection(self, other: Line) -> Point: diff --git a/src/freeroute/geometry/polygon.py b/src/freeroute/geometry/polygon.py new file mode 100644 index 0000000..adb0795 --- /dev/null +++ b/src/freeroute/geometry/polygon.py @@ -0,0 +1,77 @@ +"""``Polygon`` — a normalized corner list (no duplicate or collinear corners). + +Ports ``geometry/planar/Polygon.java``. Used as the corner normalizer behind +:class:`~freeroute.geometry.polygon_shape.PolygonShape`. + +**Arithmetic model:** corner de-duplication and collinearity removal are +**exact** (exact ``Point.side_of``); the winding-number test is approximate +(sum of ``angle_approx``), exactly as upstream. +""" + +from __future__ import annotations + +import math + +from .point import Point +from .side import Side + +__all__ = ["Polygon"] + + +class Polygon: + """A list of corners with no two consecutive equal and no three collinear.""" + + __slots__ = ("corners",) + + def __init__(self, points) -> None: + corners: list[Point] = list(points) + corner_removed = True + while corner_removed: + corner_removed = False + if not corners: + break + # remove consecutive duplicate points + deduped: list[Point] = [corners[0]] + for pt in corners[1:]: + if pt != deduped[-1]: + deduped.append(pt) + if len(deduped) != len(corners): + corner_removed = True + corners = deduped + # remove a point collinear with its previous and next point + if len(corners) >= 3: + for i in range(1, len(corners) - 1): + if corners[i].side_of(corners[i - 1], corners[i + 1]) == Side.COLLINEAR: + del corners[i] + corner_removed = True + break + self.corners: list[Point] = corners + + def corner_array(self) -> list[Point]: + return list(self.corners) + + def revert_corners(self) -> Polygon: + return Polygon(list(reversed(self.corners))) + + def winding_number_after_closing(self) -> int: + """Winding number of the closed polygon: >0 counter-clockwise, <0 clockwise.""" + arr = self.corners + if len(arr) < 2: + return 0 + first_side = arr[1].difference_by(arr[0]) + prev_side = first_side + corner_count = len(arr) + if arr[0] == arr[corner_count - 1]: + corner_count -= 1 + angle_sum = 0.0 + for i in range(1, corner_count + 1): + if i == corner_count - 1: + next_side = arr[0].difference_by(arr[i]) + elif i == corner_count: + next_side = first_side + else: + next_side = arr[i + 1].difference_by(arr[i]) + angle_sum += prev_side.angle_approx(next_side) + prev_side = next_side + angle_sum /= 2.0 * math.pi + return round(angle_sum) diff --git a/src/freeroute/geometry/polygon_shape.py b/src/freeroute/geometry/polygon_shape.py new file mode 100644 index 0000000..5ee36a6 --- /dev/null +++ b/src/freeroute/geometry/polygon_shape.py @@ -0,0 +1,262 @@ +"""``PolygonShape`` — a simple (possibly concave) polygon and ``split_to_convex``. + +Ports ``geometry/planar/PolygonShape.java``, in particular the recursive +``split_to_convex`` that decomposes a simple polygon into convex +:class:`~freeroute.geometry.simplex.Simplex` tiles — the decomposition the +router's search tree relies on. + +**Arithmetic model:** corner orientation / convexity tests are **exact** +(exact ``Point.side_of``); the division-point search is **approximate** (it uses +``float`` line evaluations and rounds the split point to an integer corner, as +upstream). Point containment is exact and delegated to the convex tiles. + +The recursion starts the concave-corner search at corner 0 (deterministic); +upstream seeds a fixed PRNG to vary the start. The choice only affects *which* +valid decomposition is produced, never its correctness. +""" + +from __future__ import annotations + +from .float_point import FloatPoint +from .line import Line +from .point import Point +from .polygon import Polygon +from .side import Side +from .simplex import Simplex + +__all__ = ["PolygonShape"] + + +class PolygonShape: + """A simple polygon, normalized to a canonical counter-clockwise corner list.""" + + __slots__ = ("corners", "_convex_pieces") + + def __init__(self, points) -> None: + poly = Polygon(points) + if poly.winding_number_after_closing() < 0: + poly = poly.revert_corners() # make counter-clockwise + curr = poly.corner_array() + last = len(curr) - 1 + if last > 0 and curr[0] == curr[last]: + last -= 1 + if last >= 2 and curr[last].side_of(curr[last - 1], curr[0]) == Side.COLLINEAR: + last -= 1 + first = 0 + if last - first >= 2 and curr[0].side_of(curr[1], curr[last]) == Side.COLLINEAR: + first += 1 + # rotate so the corner with the lowest y (then lowest x) comes first + start = first + start_c = curr[start].to_float() + for i in range(start + 1, last + 1): + c = curr[i].to_float() + if c.y < start_c.y or (c.y == start_c.y and c.x < start_c.x): + start = i + start_c = c + self.corners: list[Point] = curr[start : last + 1] + curr[first:start] + self._convex_pieces: list[Simplex] | None = None + + def __repr__(self) -> str: + return f"PolygonShape(<{len(self.corners)} corners>)" + + # --- accessors ---------------------------------------------------------- + + def corner(self, no: int) -> Point: + return self.corners[no] + + def border_line_count(self) -> int: + return len(self.corners) + + def border_line(self, no: int) -> Line: + nxt = self.corners[0] if no == len(self.corners) - 1 else self.corners[no + 1] + return Line(self.corners[no], nxt) + + def dimension(self) -> int: + n = len(self.corners) + if n == 0: + return -1 + if n == 1: + return 0 + if n == 2: + return 1 + return 2 + + def is_empty(self) -> bool: + return len(self.corners) == 0 + + def area(self) -> float: + if self.dimension() <= 2 and len(self.corners) < 3: + return 0.0 + result = 0.0 + corners = [c.to_float() for c in self.corners] + n = len(corners) + prev_corner = corners[n - 2] + curr_corner = corners[n - 1] + for i in range(n): + next_corner = corners[i] + result += curr_corner.x * (next_corner.y - prev_corner.y) + prev_corner = curr_corner + curr_corner = next_corner + return 0.5 * abs(result) + + # --- convexity ---------------------------------------------------------- + + def is_convex(self) -> bool: + corners = self.corners + n = len(corners) + if n <= 2: + return True + prev_point = corners[n - 1] + curr_point = corners[0] + next_point = corners[1] + for ind in range(n): + if next_point.side_of(prev_point, curr_point) == Side.ON_THE_RIGHT: + return False + prev_point = curr_point + curr_point = next_point + next_point = corners[0] if ind == n - 2 else corners[(ind + 2) % n] + return True + + # --- containment (via convex decomposition) ---------------------------- + + def is_outside(self, point: Point) -> bool: + return all(piece.is_outside(point) for piece in self.split_to_convex()) + + def contains(self, point: Point) -> bool: + return not self.is_outside(point) + + def contains_inside(self, point: Point) -> bool: + # strictly inside some tile and not merely on a shared division edge + pieces = self.split_to_convex() + return any(p.contains_inside(point) for p in pieces) + + # --- the key operation: split into convex tiles ------------------------ + + def split_to_convex(self) -> list[Simplex]: + if self._convex_pieces is None: + pieces = _split_recu(self.corners) + if pieces is None: + self._convex_pieces = [] + else: + self._convex_pieces = [Simplex.from_corners(p) for p in pieces] + return self._convex_pieces + + +def _split_recu(corners: list[Point]) -> list[list[Point]] | None: + """Recursively divide a simple polygon into convex corner lists.""" + n = len(corners) + if n < 3: + return [corners] + + # find the first concave corner + concave_no = -1 + for start in range(n): + prev_corner = corners[start - 1] + curr_corner = corners[start] + next_corner = corners[(start + 1) % n] + if next_corner.side_of(prev_corner, curr_corner) == Side.ON_THE_RIGHT: + concave_no = start + break + + if concave_no < 0: + return [corners] # already convex + + projection, corner_no_after = _division_point(corners, concave_no) + if projection is None: + return None # self-intersecting or degenerate + + split_pt = projection.round() + + # first piece: concave_no .. corner_no_after (exclusive) + split point + count = corner_no_after - concave_no + if count < 0: + count += n + count += 1 + first_arr: list[Point] = [] + idx = concave_no + for _ in range(count - 1): + first_arr.append(corners[idx]) + idx = (idx + 1) % n + first_arr.append(split_pt) + + # last piece: split point + corner_no_after .. concave_no + count = concave_no - corner_no_after + if count < 0: + count += n + count += 2 + last_arr: list[Point] = [split_pt] + idx = corner_no_after + for _ in range(1, count): + last_arr.append(corners[idx]) + idx = (idx + 1) % n + + c1 = _split_recu(PolygonShape(first_arr).corners) + if c1 is None: + return None + c2 = _split_recu(PolygonShape(last_arr).corners) + if c2 is None: + return None + return c1 + c2 + + +def _division_point(corners: list[Point], concave_no: int) -> tuple[FloatPoint | None, int]: + """Minimal axis-parallel division line from a concave corner (approximate).""" + n = len(corners) + concave = corners[concave_no].to_float() + before = corners[concave_no - 1].to_float() + after = corners[(concave_no + 1) % n].to_float() + + search_right = before.y > concave.y or concave.y > after.y + search_left = before.y < concave.y or concave.y < after.y + search_up = before.x < concave.x or concave.x < after.x + search_down = before.x > concave.x or concave.x > after.x + + min_dist = float("inf") + min_projection: FloatPoint | None = None + corner_no_after_min = 0 + + idx_after = (concave_no + 2) % n + corner_before = corners[idx_after - 1] if idx_after != 0 else corners[n - 1] + before_approx = corner_before.to_float() + + for _ in range(n - 2): + corner_after = corners[idx_after] + after_approx = corner_after.to_float() + + if before_approx.y != after_approx.y: # horizontal division + lo_y = min(before_approx.y, after_approx.y) + hi_y = max(before_approx.y, after_approx.y) + if lo_y <= concave.y <= hi_y: + line = Line(corner_before, corner_after) + x_intersect = line.function_in_y_value_approx(concave.y) + dist = abs(x_intersect - concave.x) + ok = dist < min_dist and ( + (search_right and x_intersect > concave.x and concave.y <= after_approx.y) + or (search_left and x_intersect < concave.x and concave.y >= after_approx.y) + ) + if ok: + min_dist = dist + corner_no_after_min = idx_after + min_projection = FloatPoint(x_intersect, concave.y) + + if before_approx.x != after_approx.x: # vertical division + lo_x = min(before_approx.x, after_approx.x) + hi_x = max(before_approx.x, after_approx.x) + if lo_x <= concave.x <= hi_x: + line = Line(corner_before, corner_after) + y_intersect = line.function_value_approx(concave.x) + dist = abs(y_intersect - concave.y) + ok = dist < min_dist and ( + (search_up and y_intersect > concave.y and concave.x >= after_approx.x) + or (search_down and y_intersect < concave.y and concave.x <= after_approx.x) + ) + if ok: + min_dist = dist + corner_no_after_min = idx_after + min_projection = FloatPoint(concave.x, y_intersect) + + corner_before = corner_after + before_approx = after_approx + idx_after = 0 if idx_after == n - 1 else idx_after + 1 + + return min_projection, corner_no_after_min diff --git a/src/freeroute/geometry/simplex.py b/src/freeroute/geometry/simplex.py index 6fba858..844387a 100644 --- a/src/freeroute/geometry/simplex.py +++ b/src/freeroute/geometry/simplex.py @@ -52,6 +52,18 @@ class Simplex(TileShape): arr.sort() return Simplex(arr)._remove_redundant_lines() + @staticmethod + def from_corners(corners) -> Simplex: + """Build a normalized simplex from the corners of a convex polygon. + + Ports ``TileShape.get_instance(Point[])``: one directed border line per + edge (corner ``i`` to ``i+1``), then normalized. + """ + corners = list(corners) + n = len(corners) + lines = [Line(corners[j], corners[(j + 1) % n]) for j in range(n)] + return Simplex.get_instance(lines) + @staticmethod def empty() -> Simplex: return Simplex(()) diff --git a/tests/geometry/test_polygon.py b/tests/geometry/test_polygon.py new file mode 100644 index 0000000..eccbf2d --- /dev/null +++ b/tests/geometry/test_polygon.py @@ -0,0 +1,148 @@ +"""Invariant tests for Polygon / PolygonShape.split_to_convex. + +No external oracle exists, so these assert the decomposition invariants from the +phase-2 brief for several concave polygons: the convex tiles' areas sum to the +polygon area (no gaps, no overlap), a point is contained in the polygon iff it +is contained in some tile, and no point is strictly inside more than one tile +(interiors disjoint). + +Source: ``geometry/planar/{Polygon,PolygonShape}.java``. +""" + +from __future__ import annotations + +import pytest + +from freeroute.geometry import IntPoint, Side +from freeroute.geometry.polygon import Polygon +from freeroute.geometry.polygon_shape import PolygonShape + + +def pts(*coords): + return [IntPoint(x, y) for x, y in coords] + + +# --- Polygon normalization -------------------------------------------------- + + +def test_polygon_removes_consecutive_duplicates(): + poly = Polygon(pts((0, 0), (0, 0), (10, 0), (10, 10), (10, 10), (0, 10))) + assert len(poly.corners) == 4 + + +def test_polygon_removes_middle_collinear_corner(): + # (5, 0) is collinear with its neighbours (0,0) and (10,0) -> dropped. + # Raw Polygon only removes *middle* collinear corners (endpoint-collinear + # removal is PolygonShape's job), so the trailing (0,5) is kept. + poly = Polygon(pts((0, 0), (5, 0), (10, 0), (10, 10), (0, 10), (0, 5))) + coords = [(c.x, c.y) for c in poly.corners] + assert (5, 0) not in coords + assert coords == [(0, 0), (10, 0), (10, 10), (0, 10), (0, 5)] + + +def test_polygon_shape_removes_endpoint_collinear_corner(): + # PolygonShape additionally drops the endpoint-collinear (0,5). + ps = PolygonShape(pts((0, 0), (5, 0), (10, 0), (10, 10), (0, 10), (0, 5))) + coords = {(c.x, c.y) for c in ps.corners} + assert coords == {(0, 0), (10, 0), (10, 10), (0, 10)} + + +def test_winding_number_sign(): + ccw = Polygon(pts((0, 0), (10, 0), (10, 10), (0, 10))) + cw = Polygon(pts((0, 0), (0, 10), (10, 10), (10, 0))) + assert ccw.winding_number_after_closing() > 0 + assert cw.winding_number_after_closing() < 0 + + +# --- convexity -------------------------------------------------------------- + + +def test_is_convex(): + square = PolygonShape(pts((0, 0), (10, 0), (10, 10), (0, 10))) + assert square.is_convex() + ell = PolygonShape(pts((0, 0), (60, 0), (60, 20), (20, 20), (20, 60), (0, 60))) + assert not ell.is_convex() + + +# --- split_to_convex invariants --------------------------------------------- + +SHAPES = { + "L": pts((0, 0), (60, 0), (60, 20), (20, 20), (20, 60), (0, 60)), + "plus": pts( + (20, 0), + (40, 0), + (40, 20), + (60, 20), + (60, 40), + (40, 40), + (40, 60), + (20, 60), + (20, 40), + (0, 40), + (0, 20), + (20, 20), + ), + "staircase": pts((0, 0), (30, 0), (30, 10), (20, 10), (20, 20), (10, 20), (10, 30), (0, 30)), + "square": pts((0, 0), (40, 0), (40, 40), (0, 40)), +} + + +def _bounds(corners): + xs = [c.x for c in corners] + ys = [c.y for c in corners] + return min(xs) - 3, max(xs) + 3, min(ys) - 3, max(ys) + 3 + + +@pytest.mark.parametrize("name", list(SHAPES)) +def test_split_area_conserved(name): + ps = PolygonShape(SHAPES[name]) + tiles = ps.split_to_convex() + assert tiles, "expected at least one convex tile" + assert sum(t.area() for t in tiles) == pytest.approx(ps.area()) + + +@pytest.mark.parametrize("name", list(SHAPES)) +def test_split_union_equals_polygon(name): + ps = PolygonShape(SHAPES[name]) + tiles = ps.split_to_convex() + lo_x, hi_x, lo_y, hi_y = _bounds(ps.corners) + for x in range(lo_x, hi_x + 1): + for y in range(lo_y, hi_y + 1): + p = IntPoint(x, y) + in_poly = ps.contains(p) + in_some_tile = any(t.contains(p) for t in tiles) + assert in_poly == in_some_tile, (x, y) + + +@pytest.mark.parametrize("name", list(SHAPES)) +def test_split_interiors_are_disjoint(name): + ps = PolygonShape(SHAPES[name]) + tiles = ps.split_to_convex() + lo_x, hi_x, lo_y, hi_y = _bounds(ps.corners) + for x in range(lo_x, hi_x + 1): + for y in range(lo_y, hi_y + 1): + p = IntPoint(x, y) + strict = sum(1 for t in tiles if t.contains_inside(p)) + assert strict <= 1, (x, y, strict) + + +def test_convex_polygon_is_single_tile(): + ps = PolygonShape(SHAPES["square"]) + tiles = ps.split_to_convex() + assert len(tiles) == 1 + assert tiles[0].dimension() == 2 + + +def test_concave_corner_detected_by_side_of(): + # the reflex corner of the L is (20, 20): its next corner is on the right + ell = PolygonShape(SHAPES["L"]) + corners = ell.corners + reflex_found = False + n = len(corners) + for i in range(n): + prev_c = corners[i - 1] + curr = corners[i] + nxt = corners[(i + 1) % n] + if nxt.side_of(prev_c, curr) == Side.ON_THE_RIGHT: + reflex_found = True + assert reflex_found