From 55a869ef7d4966f65ab29ec8b4e6263a3a65623d Mon Sep 17 00:00:00 2001 From: Ryan Malloy Date: Sat, 11 Jul 2026 18:01:29 -0600 Subject: [PATCH] Add directions and lines with exact rational intersection Ports Direction/IntDirection (equivalence classes of vectors, gcd- normalized, exact angular compare) and Line. Line.side_of uses an exact integer determinant; Line.intersection returns an IntPoint when the crossing is integral and a RationalPoint otherwise, with the orthogonal and 45-degree fast paths from the source preserved. Parallel lines yield a point at infinity (z=0). BigIntDirection is folded into IntDirection since unbounded int always fits. Tests cover integral and rational intersections, parallel-line infinity, exactness beyond double precision (verified via exact collinearity of the result), and direction normalization/ordering. --- src/freeroute/geometry/direction.py | 152 +++++++++++++++++++++ src/freeroute/geometry/line.py | 197 ++++++++++++++++++++++++++++ tests/geometry/test_direction.py | 61 +++++++++ tests/geometry/test_line.py | 102 ++++++++++++++ 4 files changed, 512 insertions(+) create mode 100644 src/freeroute/geometry/direction.py create mode 100644 src/freeroute/geometry/line.py create mode 100644 tests/geometry/test_direction.py create mode 100644 tests/geometry/test_line.py diff --git a/src/freeroute/geometry/direction.py b/src/freeroute/geometry/direction.py new file mode 100644 index 0000000..4698370 --- /dev/null +++ b/src/freeroute/geometry/direction.py @@ -0,0 +1,152 @@ +"""Directions: ``Direction`` base and ``IntDirection``. + +Ports ``geometry/planar/{Direction,IntDirection}.java``. A direction is an +equivalence class of vectors pointing the same way — FreeRouting prefers these +over angles because angle arithmetic is inexact. + +**Arithmetic model: exact, Python ``int``.** ``BigIntDirection`` is folded away: +with unbounded ``int`` a normalized direction always fits in +:class:`IntDirection`. The angular :meth:`IntDirection.compare_to` uses an exact +integer determinant (upstream uses a ``double`` for speed). +""" + +from __future__ import annotations + +from .side import Side, Signum +from .vector import IntVector, Vector + +__all__ = ["Direction", "IntDirection"] + + +class Direction: + """Abstract base for plane directions.""" + + def get_vector(self) -> Vector: # pragma: no cover - overridden + raise NotImplementedError + + def is_orthogonal(self) -> bool: # pragma: no cover - overridden + raise NotImplementedError + + def is_diagonal(self) -> bool: # pragma: no cover - overridden + raise NotImplementedError + + def is_multiple_of_45_degree(self) -> bool: + return self.is_orthogonal() or self.is_diagonal() + + def turn_45_degree(self, factor: int) -> Direction: # pragma: no cover + raise NotImplementedError + + def opposite(self) -> Direction: # pragma: no cover - overridden + raise NotImplementedError + + def side_of(self, other: Direction) -> Side: + return self.get_vector().side_of(other.get_vector()) + + def projection(self, other: Direction) -> Signum: + return self.get_vector().projection(other.get_vector()) + + def compare_to(self, other: Direction) -> int: # pragma: no cover - overridden + raise NotImplementedError + + def equals(self, other: Direction | None) -> bool: + if self is other: + return True + if other is None: + return False + if self.side_of(other) != Side.COLLINEAR: + return False + # collinear — reject the opposite direction. + return self.get_vector().projection(other.get_vector()) == Signum.POSITIVE + + @staticmethod + def get_instance(vector: Vector) -> Direction: + return vector.to_normalized_direction() + + @staticmethod + def from_points(p_from, p_to) -> Direction | None: + if p_from == p_to: + return None + return Direction.get_instance(p_to.difference_by(p_from)) + + +class IntDirection(Direction): + """A direction represented by an integer vector (not necessarily reduced).""" + + __slots__ = ("x", "y") + + def __init__(self, x: int, y: int) -> None: + self.x = int(x) + self.y = int(y) + + def __repr__(self) -> str: + return f"IntDirection({self.x}, {self.y})" + + def get_vector(self) -> IntVector: + return IntVector(self.x, self.y) + + def is_orthogonal(self) -> bool: + return self.x == 0 or self.y == 0 + + def is_diagonal(self) -> bool: + return abs(self.x) == abs(self.y) + + def opposite(self) -> IntDirection: + return IntDirection(-self.x, -self.y) + + def turn_45_degree(self, factor: int) -> IntDirection: + n = factor % 8 + x, y = self.x, self.y + table = { + 0: (x, y), + 1: (x - y, x + y), + 2: (-y, x), + 3: (-x - y, x - y), + 4: (-x, -y), + 5: (y - x, -x - y), + 6: (y, -x), + 7: (x + y, y - x), + } + nx, ny = table[n] + return IntDirection(nx, ny) + + def compare_to(self, other: Direction) -> int: + """Angular comparison with the positive x-axis (exact). + + Returns +1 if ``self`` has a strictly larger angle than ``other``, 0 if + equal, -1 otherwise. Ports the half-plane split in + ``IntDirection.compareTo``. + """ + if not isinstance(other, IntDirection): + return -other.compare_to(self) + y, x = self.y, self.x + oy, ox = other.y, other.x + if y > 0: + if oy < 0: + return -1 + if oy == 0: + return 1 if ox > 0 else -1 + elif y < 0: + if oy >= 0: + return 1 + else: # y == 0 + if x > 0: + if oy != 0 or ox < 0: + return -1 + return 0 + # x < 0 + if oy > 0 or (oy == 0 and ox > 0): + return 1 + if oy < 0: + return -1 + return 0 + # same open horizontal half-plane: compare by exact determinant. + determinant = ox * y - oy * x + return _sign(determinant) + + +def _sign(value: int) -> int: + if value > 0: + return 1 + if value < 0: + return -1 + return 0 diff --git a/src/freeroute/geometry/line.py b/src/freeroute/geometry/line.py new file mode 100644 index 0000000..dfe8ea5 --- /dev/null +++ b/src/freeroute/geometry/line.py @@ -0,0 +1,197 @@ +"""``Line`` — a directed line through two integer points. + +Ports ``geometry/planar/Line.java`` (the subset needed for the router's +half-plane work). Like the upstream class, a ``Line`` is defined by two +:class:`~freeroute.geometry.point.IntPoint` endpoints. + +**Arithmetic model: exact.** :meth:`side_of` uses an exact integer determinant, +and :meth:`intersection` returns an :class:`~freeroute.geometry.point.IntPoint` +when the crossing is integral and a +:class:`~freeroute.geometry.point.RationalPoint` otherwise (``z == 0`` when the +lines are parallel — the "point at infinity"). :meth:`intersection_approx` is +the fast ``float`` path. +""" + +from __future__ import annotations + +import math + +from .direction import Direction +from .float_point import FloatPoint +from .point import IntPoint, Point, rational_point +from .side import Side +from .vector import IntVector, Vector + +__all__ = ["Line"] + + +class Line: + """A directed line ``a -> b`` with integer endpoints.""" + + __slots__ = ("a", "b", "_dir") + + def __init__(self, a: Point, b: Point) -> None: + self.a = a + self.b = b + self._dir: Direction | None = None + + @staticmethod + def from_coords(ax: int, ay: int, bx: int, by: int) -> Line: + return Line(IntPoint(ax, ay), IntPoint(bx, by)) + + @staticmethod + def from_point_direction(a: Point, direction: Direction) -> Line: + return Line(a, a.translate_by(direction.get_vector())) + + def __eq__(self, other: object) -> bool: + if not isinstance(other, Line): + return NotImplemented + if self.side_of(other.a) != Side.COLLINEAR: + return False + return self.direction().equals(other.direction()) + + def __hash__(self) -> int: + return hash((self.a, self.b)) + + def __repr__(self) -> str: + return f"Line({self.a!r}, {self.b!r})" + + def direction(self) -> Direction: + if self._dir is None: + self._dir = Direction.get_instance(self.b.difference_by(self.a)) + return self._dir + + def opposite(self) -> Line: + return Line(self.b, self.a) + + # --- orientation -------------------------------------------------------- + + def side_of(self, point: Point) -> Side: + """Side of this line relative to ``point`` (exact). + + ``ON_THE_LEFT`` if the line passes to the left of the point, + ``COLLINEAR`` if the point lies on the line. + """ + return point.side_of(self.a, self.b).negate() + + def side_of_float(self, point: FloatPoint, tolerance: float = 0.0) -> Side: + """Approximate side test for a :class:`FloatPoint` with a tolerance band.""" + ax, ay = self.a.x, self.a.y + det = (self.b.y - ay) * (point.x - ax) - (self.b.x - ax) * (point.y - ay) + if det - tolerance > 0: + return Side.ON_THE_LEFT + if det + tolerance < 0: + return Side.ON_THE_RIGHT + return Side.COLLINEAR + + def signed_distance(self, point: FloatPoint) -> float: + ax, ay = self.a.x, self.a.y + dx = self.b.x - ax + dy = self.b.y - ay + det = dy * (point.x - ax) - dx * (point.y - ay) + return det / math.sqrt(dx * dx + dy * dy) + + def is_orthogonal(self) -> bool: + return self.direction().is_orthogonal() + + def is_diagonal(self) -> bool: + return self.direction().is_diagonal() + + def is_multiple_of_45_degree(self) -> bool: + return self.direction().is_multiple_of_45_degree() + + def is_parallel(self, other: Line) -> bool: + return self.direction().side_of(other.direction()) == Side.COLLINEAR + + def is_perpendicular(self, other: Line) -> bool: + v1 = self.direction().get_vector() + v2 = other.direction().get_vector() + return v1.projection(v2).value == 0 + + def overlaps(self, other: Line) -> bool: + return self.side_of(other.a) == Side.COLLINEAR and self.side_of(other.b) == Side.COLLINEAR + + is_equal_or_opposite = overlaps + + # --- translation -------------------------------------------------------- + + def translate_by(self, vector: Vector) -> Line: + if vector.is_zero(): + return self + return Line(self.a.translate_by(vector), self.b.translate_by(vector)) + + # --- intersection (exact) ---------------------------------------------- + + def intersection(self, other: Line) -> Point: + """Exact intersection point; ``result.is_infinite()`` iff parallel.""" + a = self.a + b = self.b + oa = other.a + ob = other.b + delta_1 = b.difference_by(a) + delta_2 = ob.difference_by(oa) + if not isinstance(delta_1, IntVector) or not isinstance(delta_2, IntVector): + raise NotImplementedError("Line.intersection only supports integer lines") + + fast = _fast_intersection(a, oa, delta_1, delta_2) + if fast is not None: + return fast + + det_1 = a.determinant(b) + det_2 = oa.determinant(ob) + det = delta_2.determinant(delta_1) + is_x = det_1 * delta_2.x - det_2 * delta_1.x + is_y = det_1 * delta_2.y - det_2 * delta_1.y + if det == 0: + return rational_point(is_x, is_y, 0) + return rational_point(is_x, is_y, det) + + def intersection_approx(self, other: Line) -> FloatPoint: + """Fast ``float`` intersection; parallel lines give a huge coordinate.""" + a, b = self.a, self.b + oa, ob = other.a, other.b + d1x = b.x - a.x + d1y = b.y - a.y + d2x = ob.x - oa.x + d2y = ob.y - oa.y + det_1 = a.x * b.y - a.y * b.x + det_2 = oa.x * ob.y - oa.y * ob.x + det = d2x * d1y - d2y * d1x + if det == 0: + big = float(2**31 - 1) + return FloatPoint(big, big) + return FloatPoint( + (d2x * det_1 - d1x * det_2) / det, + (d2y * det_1 - d1y * det_2) / det, + ) + + +def _fast_intersection( + a: IntPoint, oa: IntPoint, delta_1: IntVector, delta_2: IntVector +) -> Point | None: + """Closed-form intersection for the orthogonal/45-degree combinations.""" + if delta_1.x == 0: # this line vertical + if delta_2.y == 0: # other horizontal + return IntPoint(a.x, oa.y) + if delta_2.x == delta_2.y: # other right diagonal + return IntPoint(a.x, oa.y + a.x - oa.x) + if delta_2.x == -delta_2.y: # other left diagonal + return IntPoint(a.x, oa.y + oa.x - a.x) + elif delta_1.y == 0: # this line horizontal + if delta_2.x == 0: # other vertical + return IntPoint(oa.x, a.y) + if delta_2.x == delta_2.y: # other right diagonal + return IntPoint(oa.x + a.y - oa.y, a.y) + if delta_2.x == -delta_2.y: # other left diagonal + return IntPoint(oa.x + oa.y - a.y, a.y) + elif delta_1.x == delta_1.y: # this right diagonal + if delta_2.x == 0: # other vertical + return IntPoint(oa.x, a.y + oa.x - a.x) + if delta_2.y == 0: # other horizontal + return IntPoint(a.x + oa.y - a.y, oa.y) + elif delta_1.x == -delta_1.y: # this left diagonal + if delta_2.x == 0: # other vertical + return IntPoint(oa.x, a.y + a.x - oa.x) + if delta_2.y == 0: # other horizontal + return IntPoint(a.x + a.y - oa.y, oa.y) + return None diff --git a/tests/geometry/test_direction.py b/tests/geometry/test_direction.py new file mode 100644 index 0000000..892003c --- /dev/null +++ b/tests/geometry/test_direction.py @@ -0,0 +1,61 @@ +"""Tests for directions. + +Oracles hand-derived from ``geometry/planar/{Direction,IntDirection}.java``. +""" + +from __future__ import annotations + +from freeroute.geometry import Direction, IntDirection, IntPoint, IntVector, Side + + +def test_normalized_direction_reduces_by_gcd(): + d = IntVector(2, 4).to_normalized_direction() + assert isinstance(d, IntDirection) + assert (d.x, d.y) == (1, 2) + + +def test_normalized_direction_beyond_java_int(): + # gcd reduction stays exact for coordinates far past Java's int range. + d = IntVector(6 * 10**12, 9 * 10**12).to_normalized_direction() + assert (d.x, d.y) == (2, 3) + + +def test_turn_45_degree_cycles(): + right = IntDirection(1, 0) + assert (right.turn_45_degree(1).x, right.turn_45_degree(1).y) == (1, 1) + assert (right.turn_45_degree(2).x, right.turn_45_degree(2).y) == (0, 1) + assert (right.turn_45_degree(4).x, right.turn_45_degree(4).y) == (-1, 0) + assert (right.turn_45_degree(8).x, right.turn_45_degree(8).y) == (1, 0) + + +def test_opposite(): + o = IntDirection(1, 2).opposite() + assert (o.x, o.y) == (-1, -2) + + +def test_direction_equals_ignores_magnitude_not_sign(): + a = IntDirection(1, 1) + assert a.equals(IntDirection(3, 3)) # same ray + assert not a.equals(IntDirection(-1, -1)) # opposite ray + assert not a.equals(IntDirection(1, 0)) + + +def test_from_points_returns_none_for_equal(): + assert Direction.from_points(IntPoint(1, 1), IntPoint(1, 1)) is None + d = Direction.from_points(IntPoint(0, 0), IntPoint(4, 4)) + assert (d.x, d.y) == (1, 1) + + +def test_compare_to_angular_order(): + right = IntDirection(1, 0) + up = IntDirection(0, 1) + down = IntDirection(0, -1) + # angle(right)=0, angle(up)=90, angle(down)=270 + assert right.compare_to(up) == -1 + assert up.compare_to(right) == 1 + assert up.compare_to(down) == -1 + assert right.compare_to(IntDirection(2, 0)) == 0 + + +def test_side_of_direction(): + assert IntDirection(0, 1).side_of(IntDirection(1, 0)) == Side.ON_THE_LEFT diff --git a/tests/geometry/test_line.py b/tests/geometry/test_line.py new file mode 100644 index 0000000..a4f8ed3 --- /dev/null +++ b/tests/geometry/test_line.py @@ -0,0 +1,102 @@ +"""Tests for lines, emphasizing exact intersection. + +Oracles hand-derived from ``geometry/planar/Line.java``. +""" + +from __future__ import annotations + +from freeroute.geometry import IntPoint, Line, RationalPoint, Side + + +def test_side_of_point(): + # Line.side_of reports where the *line* is relative to the point + # (point.side_of(line).negate() in the source): a point above a +x line has + # the line on its right. + line = Line.from_coords(0, 0, 2, 0) # along +x + assert line.side_of(IntPoint(1, 1)) == Side.ON_THE_RIGHT + assert line.side_of(IntPoint(1, -1)) == Side.ON_THE_LEFT + assert line.side_of(IntPoint(5, 0)) == Side.COLLINEAR + + +def test_integral_intersection_fast_paths(): + # vertical x horizontal + v = Line.from_coords(3, 0, 3, 9) + h = Line.from_coords(0, 4, 9, 4) + assert v.intersection(h) == IntPoint(3, 4) + # horizontal x right-diagonal (y=x): y=4 -> (4,4) + diag = Line.from_coords(0, 0, 5, 5) + assert h.intersection(diag) == IntPoint(4, 4) + + +def test_general_integral_intersection(): + # y = x/2 and y = 1 -> x = 2 + a = Line.from_coords(0, 0, 2, 1) + b = Line.from_coords(0, 1, 4, 1) + assert a.intersection(b) == IntPoint(2, 1) + + +def test_rational_intersection_exact(): + # y = 2x/3 and y = 1 -> x = 1.5 -> RationalPoint (3,2,2) + a = Line.from_coords(0, 0, 3, 2) + b = Line.from_coords(0, 1, 5, 1) + result = a.intersection(b) + assert isinstance(result, RationalPoint) + assert result == RationalPoint(3, 2, 2) + assert not result.is_infinite() + # float approximation agrees + fp = result.to_float() + assert abs(fp.x - 1.5) < 1e-9 and abs(fp.y - 1.0) < 1e-9 + + +def test_parallel_intersection_is_infinite(): + a = Line.from_coords(0, 0, 1, 0) + b = Line.from_coords(0, 5, 1, 5) + result = a.intersection(b) + assert result.is_infinite() + assert a.is_parallel(b) + + +def test_intersection_exact_beyond_double_precision(): + # Two nearly-parallel steep lines crossing at a non-integer point far out. + # Coordinates chosen so a double determinant would lose bits; int stays exact. + a = Line.from_coords(0, 0, 1000000, 999999) + b = Line.from_coords(0, 1, 1000000, 1000000) + result = a.intersection(b) + assert isinstance(result, RationalPoint) + # Verify exactness: the intersection lies on both lines (exact side test). + assert a.side_of(result) == Side.COLLINEAR + assert b.side_of(result) == Side.COLLINEAR + + +def test_intersection_approx_matches_exact(): + a = Line.from_coords(0, 0, 3, 2) + b = Line.from_coords(0, 1, 5, 1) + approx = a.intersection_approx(b) + assert abs(approx.x - 1.5) < 1e-9 + assert abs(approx.y - 1.0) < 1e-9 + + +def test_translate_and_opposite(): + line = Line.from_coords(0, 0, 4, 0) + assert line.opposite().direction().equals(line.direction()) is False + from freeroute.geometry import IntVector + + moved = line.translate_by(IntVector(0, 5)) + assert moved.side_of(IntPoint(2, 5)) == Side.COLLINEAR + + +def test_perpendicular_and_parallel_predicates(): + horizontal = Line.from_coords(0, 0, 1, 0) + vertical = Line.from_coords(0, 0, 0, 1) + assert horizontal.is_perpendicular(vertical) + assert not horizontal.is_parallel(vertical) + assert horizontal.is_parallel(Line.from_coords(3, 7, 10, 7)) + + +def test_line_equality_same_set_same_direction(): + a = Line.from_coords(0, 0, 2, 0) + b = Line.from_coords(-5, 0, 9, 0) # same line, same direction + c = Line.from_coords(9, 0, -5, 0) # same set, opposite direction + assert a == b + assert a != c + assert a.is_equal_or_opposite(c)