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.
This commit is contained in:
@@ -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
|
||||
@@ -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
|
||||
Reference in New Issue
Block a user