extern crate alloc;
use geometry_coords::CoordinateScalar;
use geometry_cs::{CartesianFamily, CoordinateSystem};
use geometry_tag::{LinestringTag, PointTag, PolygonTag, SameAs, SegmentTag};
use geometry_trait::{
Geometry, Linestring as LinestringTrait, Point as PointTrait, Polygon as PolygonTrait,
Ring as RingTrait, Segment as SegmentTrait, segment_end, segment_start,
};
use crate::within::{WithinPoly, WithinStrategy};
pub use crate::reversal::Reversed;
pub trait IntersectsStrategy<A, B> {
fn intersects(&self, a: &A, b: &B) -> bool;
}
#[derive(Debug, Default, Clone, Copy)]
pub struct CartesianIntersects;
impl<A, B> IntersectsStrategy<A, B> for CartesianIntersects
where
A: Geometry,
B: Geometry,
A::Kind: IntersectsPairStrategy<B::Kind>,
<A::Kind as IntersectsPairStrategy<B::Kind>>::S: IntersectsStrategy<A, B>,
{
#[inline]
fn intersects(&self, a: &A, b: &B) -> bool {
<<A::Kind as IntersectsPairStrategy<B::Kind>>::S as Default>::default().intersects(a, b)
}
}
#[derive(Debug, Default, Clone, Copy)]
pub struct IxPointPoint;
#[derive(Debug, Default, Clone, Copy)]
pub struct IxPointSegment;
#[derive(Debug, Default, Clone, Copy)]
pub struct IxSegmentSegment;
#[derive(Debug, Default, Clone, Copy)]
pub struct IxLinestringSegment;
#[derive(Debug, Default, Clone, Copy)]
pub struct IxLinestringLinestring;
#[derive(Debug, Default, Clone, Copy)]
pub struct IxPointPolygon;
#[derive(Debug, Default, Clone, Copy)]
pub struct IxLinestringPolygon;
#[derive(Debug, Default, Clone, Copy)]
pub struct IxPolygonPolygon;
#[derive(Debug, Default, Clone, Copy)]
pub struct IxSegmentPoint;
#[derive(Debug, Default, Clone, Copy)]
pub struct IxSegmentLinestring;
#[derive(Debug, Default, Clone, Copy)]
pub struct IxPolygonPoint;
#[derive(Debug, Default, Clone, Copy)]
pub struct IxPolygonLinestring;
impl<A, B> IntersectsStrategy<A, B> for IxPointPoint
where
A: PointTrait,
B: PointTrait<Scalar = A::Scalar>,
<A::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
{
#[inline]
fn intersects(&self, a: &A, b: &B) -> bool {
points_equal::<A, B>(a, b, A::DIM)
}
}
impl<P, S> IntersectsStrategy<P, S> for IxPointSegment
where
P: PointTrait,
S: SegmentTrait<Point = P>,
P: geometry_trait::PointMut + Default,
P::Scalar: CoordinateScalar,
<P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
{
#[inline]
fn intersects(&self, p: &P, s: &S) -> bool {
point_on_segment(p, &segment_start(s), &segment_end(s))
}
}
impl<A, B, P> IntersectsStrategy<A, B> for IxSegmentSegment
where
A: SegmentTrait<Point = P>,
B: SegmentTrait<Point = P>,
P: PointTrait + geometry_trait::PointMut + Default,
P::Scalar: CoordinateScalar,
<P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
{
#[inline]
fn intersects(&self, a: &A, b: &B) -> bool {
segments_intersect(
&segment_start(a),
&segment_end(a),
&segment_start(b),
&segment_end(b),
)
}
}
impl<L, S, P> IntersectsStrategy<L, S> for IxLinestringSegment
where
L: LinestringTrait<Point = P>,
S: SegmentTrait<Point = P>,
P: PointTrait + geometry_trait::PointMut + Default,
P::Scalar: CoordinateScalar,
<P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
{
#[inline]
fn intersects(&self, ls: &L, s: &S) -> bool {
let s1 = segment_start(s);
let s2 = segment_end(s);
let mut it = ls.points();
let Some(mut prev) = it.next() else {
return false;
};
for curr in it {
if segments_intersect(prev, curr, &s1, &s2) {
return true;
}
prev = curr;
}
false
}
}
impl<A, B, P> IntersectsStrategy<A, B> for IxLinestringLinestring
where
A: LinestringTrait<Point = P>,
B: LinestringTrait<Point = P>,
P: PointTrait,
P::Scalar: CoordinateScalar,
<P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
{
#[inline]
fn intersects(&self, a: &A, b: &B) -> bool {
let mut ia = a.points();
let Some(mut pa) = ia.next() else {
return false;
};
for qa in ia {
let mut ib = b.points();
let Some(mut pb) = ib.next() else {
return false;
};
for qb in ib {
if segments_intersect(pa, qa, pb, qb) {
return true;
}
pb = qb;
}
pa = qa;
}
false
}
}
impl<P, G> IntersectsStrategy<P, G> for IxPointPolygon
where
P: PointTrait,
G: PolygonTrait<Point = P>,
P::Scalar: CoordinateScalar,
<P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
{
#[inline]
fn intersects(&self, p: &P, pg: &G) -> bool {
WithinPoly.covered_by(p, pg)
}
}
impl<L, G, P> IntersectsStrategy<L, G> for IxLinestringPolygon
where
L: LinestringTrait<Point = P>,
G: PolygonTrait<Point = P>,
P: PointTrait,
P::Scalar: CoordinateScalar,
<P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
{
fn intersects(&self, ls: &L, pg: &G) -> bool {
let Some(first) = ls.points().next() else {
return false;
};
if WithinPoly.covered_by(first, pg) {
return true;
}
if linestring_crosses_ring(ls, pg.exterior()) {
return true;
}
for hole in pg.interiors() {
if linestring_crosses_ring(ls, hole) {
return true;
}
}
false
}
}
impl<A, B, P> IntersectsStrategy<A, B> for IxPolygonPolygon
where
A: PolygonTrait<Point = P>,
B: PolygonTrait<Point = P>,
P: PointTrait,
P::Scalar: CoordinateScalar,
<P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
{
fn intersects(&self, a: &A, b: &B) -> bool {
if let Some(v) = a.exterior().points().next() {
if WithinPoly.covered_by(v, b) {
return true;
}
}
if let Some(v) = b.exterior().points().next() {
if WithinPoly.covered_by(v, a) {
return true;
}
}
if rings_cross(a.exterior(), b.exterior()) {
return true;
}
for hole in a.interiors() {
if rings_cross(hole, b.exterior()) {
return true;
}
}
for hole in b.interiors() {
if rings_cross(a.exterior(), hole) {
return true;
}
}
false
}
}
impl<S, P> IntersectsStrategy<S, P> for IxSegmentPoint
where
S: SegmentTrait<Point = P>,
P: PointTrait + geometry_trait::PointMut + Default,
P::Scalar: CoordinateScalar,
<P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
{
#[inline]
fn intersects(&self, s: &S, p: &P) -> bool {
IxPointSegment.intersects(p, s)
}
}
impl<S, L, P> IntersectsStrategy<S, L> for IxSegmentLinestring
where
S: SegmentTrait<Point = P>,
L: LinestringTrait<Point = P>,
P: PointTrait + geometry_trait::PointMut + Default,
P::Scalar: CoordinateScalar,
<P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
{
#[inline]
fn intersects(&self, s: &S, ls: &L) -> bool {
IxLinestringSegment.intersects(ls, s)
}
}
impl<G, P> IntersectsStrategy<G, P> for IxPolygonPoint
where
G: PolygonTrait<Point = P>,
P: PointTrait,
P::Scalar: CoordinateScalar,
<P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
{
#[inline]
fn intersects(&self, pg: &G, p: &P) -> bool {
IxPointPolygon.intersects(p, pg)
}
}
impl<G, L, P> IntersectsStrategy<G, L> for IxPolygonLinestring
where
G: PolygonTrait<Point = P>,
L: LinestringTrait<Point = P>,
P: PointTrait,
P::Scalar: CoordinateScalar,
<P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
{
#[inline]
fn intersects(&self, pg: &G, ls: &L) -> bool {
IxLinestringPolygon.intersects(ls, pg)
}
}
#[doc(hidden)]
pub trait IntersectsPairStrategy<K2> {
type S: Default;
}
impl IntersectsPairStrategy<PointTag> for PointTag {
type S = IxPointPoint;
}
impl IntersectsPairStrategy<SegmentTag> for PointTag {
type S = IxPointSegment;
}
impl IntersectsPairStrategy<SegmentTag> for SegmentTag {
type S = IxSegmentSegment;
}
impl IntersectsPairStrategy<SegmentTag> for LinestringTag {
type S = IxLinestringSegment;
}
impl IntersectsPairStrategy<LinestringTag> for LinestringTag {
type S = IxLinestringLinestring;
}
impl IntersectsPairStrategy<PolygonTag> for PointTag {
type S = IxPointPolygon;
}
impl IntersectsPairStrategy<PolygonTag> for LinestringTag {
type S = IxLinestringPolygon;
}
impl IntersectsPairStrategy<PolygonTag> for PolygonTag {
type S = IxPolygonPolygon;
}
impl IntersectsPairStrategy<PointTag> for SegmentTag {
type S = IxSegmentPoint;
}
impl IntersectsPairStrategy<LinestringTag> for SegmentTag {
type S = IxSegmentLinestring;
}
impl IntersectsPairStrategy<PointTag> for PolygonTag {
type S = IxPolygonPoint;
}
impl IntersectsPairStrategy<LinestringTag> for PolygonTag {
type S = IxPolygonLinestring;
}
#[inline]
fn points_equal<A, B>(a: &A, b: &B, dim: usize) -> bool
where
A: PointTrait,
B: PointTrait<Scalar = A::Scalar>,
{
let mut i = 0;
while i < dim {
let eq = match i {
0 => a.get::<0>() == b.get::<0>(),
1 => a.get::<1>() == b.get::<1>(),
2 => a.get::<2>() == b.get::<2>(),
3 => a.get::<3>() == b.get::<3>(),
_ => panic!("points_equal: dimension exceeds MAX_DIM (4)"),
};
if !eq {
return false;
}
i += 1;
}
true
}
fn point_on_segment<P>(p: &P, s1: &P, s2: &P) -> bool
where
P: PointTrait,
P::Scalar: CoordinateScalar,
{
let px = p.get::<0>();
let py = p.get::<1>();
let ax = s1.get::<0>();
let ay = s1.get::<1>();
let bx = s2.get::<0>();
let by = s2.get::<1>();
let cross = (bx - ax) * (py - ay) - (by - ay) * (px - ax);
if cross != P::Scalar::ZERO {
return false;
}
let (xlo, xhi) = if ax <= bx { (ax, bx) } else { (bx, ax) };
let (ylo, yhi) = if ay <= by { (ay, by) } else { (by, ay) };
xlo <= px && px <= xhi && ylo <= py && py <= yhi
}
fn segments_intersect<P>(p1: &P, p2: &P, p3: &P, p4: &P) -> bool
where
P: PointTrait,
P::Scalar: CoordinateScalar,
{
let x1 = p1.get::<0>();
let y1 = p1.get::<1>();
let x2 = p2.get::<0>();
let y2 = p2.get::<1>();
let x3 = p3.get::<0>();
let y3 = p3.get::<1>();
let x4 = p4.get::<0>();
let y4 = p4.get::<1>();
let d1 = side_sign((x3, y3), (x4, y4), (x1, y1));
let d2 = side_sign((x3, y3), (x4, y4), (x2, y2));
let d3 = side_sign((x1, y1), (x2, y2), (x3, y3));
let d4 = side_sign((x1, y1), (x2, y2), (x4, y4));
if ((d1 > 0 && d2 < 0) || (d1 < 0 && d2 > 0)) && ((d3 > 0 && d4 < 0) || (d3 < 0 && d4 > 0)) {
return true;
}
if d1 == 0 && point_on_segment(p1, p3, p4) {
return true;
}
if d2 == 0 && point_on_segment(p2, p3, p4) {
return true;
}
if d3 == 0 && point_on_segment(p3, p1, p2) {
return true;
}
if d4 == 0 && point_on_segment(p4, p1, p2) {
return true;
}
false
}
#[inline]
fn side_sign<T: CoordinateScalar>(a: (T, T), b: (T, T), c: (T, T)) -> i32 {
let v = (b.0 - a.0) * (c.1 - a.1) - (b.1 - a.1) * (c.0 - a.0);
if v > T::ZERO {
1
} else if v < T::ZERO {
-1
} else {
0
}
}
#[inline]
fn segment_bounds_disjoint<P>(p1: &P, p2: &P, p3: &P, p4: &P) -> bool
where
P: PointTrait,
{
let x1 = p1.get::<0>();
let y1 = p1.get::<1>();
let x2 = p2.get::<0>();
let y2 = p2.get::<1>();
let x3 = p3.get::<0>();
let y3 = p3.get::<1>();
let x4 = p4.get::<0>();
let y4 = p4.get::<1>();
(x1 < x3 && x1 < x4 && x2 < x3 && x2 < x4)
|| (x3 < x1 && x3 < x2 && x4 < x1 && x4 < x2)
|| (y1 < y3 && y1 < y4 && y2 < y3 && y2 < y4)
|| (y3 < y1 && y3 < y2 && y4 < y1 && y4 < y2)
}
fn linestring_crosses_ring<L, R, P>(ls: &L, r: &R) -> bool
where
L: LinestringTrait<Point = P>,
R: RingTrait<Point = P>,
P: PointTrait,
P::Scalar: CoordinateScalar,
{
let mut ils = ls.points();
let Some(mut pls) = ils.next() else {
return false;
};
for qls in ils {
if ring_edge_crosses_segment(r, pls, qls) {
return true;
}
pls = qls;
}
false
}
fn ring_edge_crosses_segment<R, P>(r: &R, pls: &P, qls: &P) -> bool
where
R: RingTrait<Point = P>,
P: PointTrait,
P::Scalar: CoordinateScalar,
{
let mut ir = r.points();
let Some(mut pr) = ir.next() else {
return false;
};
let first = pr;
let mut has_edge = false;
for qr in ir {
has_edge = true;
if !segment_bounds_disjoint(pls, qls, pr, qr) && segments_intersect(pls, qls, pr, qr) {
return true;
}
pr = qr;
}
if has_edge && pr.get::<0>() == first.get::<0>() && pr.get::<1>() == first.get::<1>() {
return false;
}
!segment_bounds_disjoint(pls, qls, pr, first) && segments_intersect(pls, qls, pr, first)
}
fn rings_cross<Ra, Rb, P>(a: &Ra, b: &Rb) -> bool
where
Ra: RingTrait<Point = P>,
Rb: RingTrait<Point = P>,
P: PointTrait,
P::Scalar: CoordinateScalar,
{
let edges_a = ring_edges(a);
let edges_b = ring_edges(b);
for (pa, qa) in &edges_a {
for (pb, qb) in &edges_b {
if segments_intersect(*pa, *qa, *pb, *qb) {
return true;
}
}
}
false
}
fn ring_edges<R>(r: &R) -> alloc::vec::Vec<(&R::Point, &R::Point)>
where
R: RingTrait,
R::Point: PointTrait,
{
let pts: alloc::vec::Vec<&R::Point> = r.points().collect();
let mut out = alloc::vec::Vec::with_capacity(pts.len());
if pts.len() < 2 {
return out;
}
for w in pts.windows(2) {
out.push((w[0], w[1]));
}
out.push((*pts.last().unwrap(), *pts.first().unwrap()));
out
}
#[cfg(test)]
mod tests {
use super::{CartesianIntersects, IntersectsStrategy, Reversed, linestring_crosses_ring};
use geometry_cs::Cartesian;
use geometry_model::{Point2D, Polygon, Ring, Segment, linestring, polygon};
type P = Point2D<f64, Cartesian>;
fn pt(x: f64, y: f64) -> P {
Point2D::new(x, y)
}
use geometry_model::Linestring;
type LS = Linestring<P>;
#[test]
fn ls_crosses_segment() {
let ls: LS = linestring![(1.0, 1.0), (3.0, 3.0), (2.0, 5.0)];
let s = Segment::new(pt(2.0, 0.0), pt(2.0, 6.0));
assert!(CartesianIntersects.intersects(&ls, &s));
}
#[test]
fn ls_touches_segment_endpoint() {
let ls: LS = linestring![(1.0, 1.0), (3.0, 3.0)];
let s = Segment::new(pt(1.0, 0.0), pt(1.0, 1.0));
assert!(CartesianIntersects.intersects(&ls, &s));
}
#[test]
fn ls_disjoint_from_segment() {
let ls: LS = linestring![(1.0, 1.0), (3.0, 3.0)];
let s = Segment::new(pt(3.0, 0.0), pt(4.0, 1.0));
assert!(!CartesianIntersects.intersects(&ls, &s));
}
#[test]
fn ls_crosses_ls() {
let a: LS = linestring![(0.0, 0.0), (2.0, 0.0), (3.0, 0.0)];
let b: LS = linestring![(0.0, 0.0), (1.0, 1.0), (2.0, 2.0)];
assert!(CartesianIntersects.intersects(&a, &b));
}
#[test]
fn ls_overlap_collinear() {
let a: LS = linestring![(0.0, 0.0), (2.0, 0.0), (3.0, 0.0)];
let b: LS = linestring![(1.0, 0.0), (4.0, 0.0), (5.0, 0.0)];
assert!(CartesianIntersects.intersects(&a, &b));
}
#[test]
fn ls_inside_polygon() {
let ls: LS = linestring![(1.0, 1.0), (2.0, 2.0)];
let p: Polygon<P> = polygon![[
(0.0, 0.0),
(10.0, 0.0),
(10.0, 10.0),
(0.0, 10.0),
(0.0, 0.0)
]];
assert!(CartesianIntersects.intersects(&ls, &p));
}
#[test]
fn ls_outside_polygon() {
let ls: LS = linestring![(11.0, 0.0), (12.0, 12.0)];
let p: Polygon<P> = polygon![[
(0.0, 0.0),
(10.0, 0.0),
(10.0, 10.0),
(0.0, 10.0),
(0.0, 0.0)
]];
assert!(!CartesianIntersects.intersects(&ls, &p));
}
#[test]
fn empty_linestring_has_no_ring_crossing() {
let ls = Linestring::<P>::new();
let ring: Ring<P> = Ring::from_vec(vec![
pt(0.0, 0.0),
pt(0.0, 1.0),
pt(1.0, 1.0),
pt(1.0, 0.0),
pt(0.0, 0.0),
]);
assert!(!linestring_crosses_ring(&ls, &ring));
}
#[test]
fn reversed_pair_compiles_and_agrees() {
let ls: LS = linestring![(1.0, 1.0), (2.0, 2.0)];
let p: Polygon<P> = polygon![[
(0.0, 0.0),
(10.0, 0.0),
(10.0, 10.0),
(0.0, 10.0),
(0.0, 0.0)
]];
let forward = CartesianIntersects.intersects(&ls, &p);
let reversed = Reversed(CartesianIntersects).intersects(&p, &ls);
assert_eq!(forward, reversed);
}
fn _accepts_readonly_point<A, B, S>(s: &S, a: &A, b: &B) -> bool
where
A: geometry_trait::Point,
B: geometry_trait::Point,
S: IntersectsStrategy<A, B>,
{
s.intersects(a, b)
}
use super::IxPointPoint;
#[test]
#[allow(
clippy::used_underscore_items,
reason = "the test exists to run the compile-time witness's body"
)]
fn readonly_point_witness_computes_equality() {
assert!(_accepts_readonly_point(
&IxPointPoint,
&pt(1.0, 1.0),
&pt(1.0, 1.0)
));
assert!(!_accepts_readonly_point(
&IxPointPoint,
&pt(1.0, 1.0),
&pt(2.0, 2.0)
));
}
}