1use robust::{orient2d, Coord as RobustCoord};
2
3use crate::precision::PrecisionModel;
4use crate::types::{Coord, LinearRing, Polygon};
5
6#[derive(Debug, Clone, Copy, PartialEq, Eq)]
7pub enum PointLocation {
8 Exterior,
9 Boundary,
10 Interior,
11}
12
13pub fn orientation(a: Coord, b: Coord, c: Coord) -> f64 {
14 orient2d(
15 RobustCoord { x: a.x, y: a.y },
16 RobustCoord { x: b.x, y: b.y },
17 RobustCoord { x: c.x, y: c.y },
18 )
19}
20
21pub fn signed_ring_area(ring: &LinearRing) -> f64 {
22 signed_area_coords(&ring.coords)
23}
24
25pub fn signed_area_coords(coords: &[Coord]) -> f64 {
26 coords
27 .windows(2)
28 .map(|pair| pair[0].x * pair[1].y - pair[1].x * pair[0].y)
29 .sum::<f64>()
30 * 0.5
31}
32
33pub fn polygon_area(polygon: &Polygon) -> f64 {
34 if polygon.is_empty() {
35 return 0.0;
36 }
37 let holes = polygon
38 .holes
39 .iter()
40 .map(|ring| signed_ring_area(ring).abs())
41 .sum::<f64>();
42 signed_ring_area(&polygon.exterior).abs() - holes
43}
44
45pub fn is_ring_ccw(ring: &LinearRing) -> bool {
46 signed_ring_area(ring) > 0.0
47}
48
49pub fn point_on_segment(point: Coord, a: Coord, b: Coord, precision: PrecisionModel) -> bool {
50 if orientation(a, b, point).abs() > precision.epsilon() {
51 return false;
52 }
53
54 point.x >= a.x.min(b.x) - precision.epsilon()
55 && point.x <= a.x.max(b.x) + precision.epsilon()
56 && point.y >= a.y.min(b.y) - precision.epsilon()
57 && point.y <= a.y.max(b.y) + precision.epsilon()
58}
59
60pub fn point_in_ring(point: Coord, ring: &LinearRing, precision: PrecisionModel) -> PointLocation {
61 if ring.coords.len() < 4 {
62 return PointLocation::Exterior;
63 }
64
65 let mut inside = false;
66 for (a, b) in ring.segments() {
67 if point_on_segment(point, a, b, precision) {
68 return PointLocation::Boundary;
69 }
70
71 let crosses = (a.y > point.y) != (b.y > point.y);
72 if crosses {
73 let x_at_y = (b.x - a.x) * (point.y - a.y) / (b.y - a.y) + a.x;
74 if x_at_y > point.x {
75 inside = !inside;
76 }
77 }
78 }
79
80 if inside {
81 PointLocation::Interior
82 } else {
83 PointLocation::Exterior
84 }
85}
86
87pub fn point_in_polygon(
88 point: Coord,
89 polygon: &Polygon,
90 precision: PrecisionModel,
91) -> PointLocation {
92 match point_in_ring(point, &polygon.exterior, precision) {
93 PointLocation::Exterior => PointLocation::Exterior,
94 PointLocation::Boundary => PointLocation::Boundary,
95 PointLocation::Interior => {
96 for hole in &polygon.holes {
97 match point_in_ring(point, hole, precision) {
98 PointLocation::Interior => return PointLocation::Exterior,
99 PointLocation::Boundary => return PointLocation::Boundary,
100 PointLocation::Exterior => {}
101 }
102 }
103 PointLocation::Interior
104 }
105 }
106}
107
108#[derive(Debug, Clone, PartialEq)]
109pub enum SegmentIntersection {
110 Point(Coord),
111 Overlap(Coord, Coord),
112}
113
114pub fn segment_intersection(
115 a1: Coord,
116 a2: Coord,
117 b1: Coord,
118 b2: Coord,
119 precision: PrecisionModel,
120) -> Option<SegmentIntersection> {
121 let o1 = orientation(a1, a2, b1);
122 let o2 = orientation(a1, a2, b2);
123 let o3 = orientation(b1, b2, a1);
124 let o4 = orientation(b1, b2, a2);
125 let eps = precision.epsilon();
126
127 if o1.abs() <= eps && o2.abs() <= eps && o3.abs() <= eps && o4.abs() <= eps {
128 return collinear_overlap(a1, a2, b1, b2, precision);
129 }
130
131 if (o1 > eps && o2 > eps) || (o1 < -eps && o2 < -eps) {
132 return None;
133 }
134 if (o3 > eps && o4 > eps) || (o3 < -eps && o4 < -eps) {
135 return None;
136 }
137
138 let denom = (a1.x - a2.x) * (b1.y - b2.y) - (a1.y - a2.y) * (b1.x - b2.x);
139 if denom.abs() <= eps {
140 return None;
141 }
142
143 let a_cross = a1.x * a2.y - a1.y * a2.x;
144 let b_cross = b1.x * b2.y - b1.y * b2.x;
145 let x = (a_cross * (b1.x - b2.x) - (a1.x - a2.x) * b_cross) / denom;
146 let y = (a_cross * (b1.y - b2.y) - (a1.y - a2.y) * b_cross) / denom;
147 let point = precision.snap_coord(Coord::new(x, y));
148
149 if point_on_segment(point, a1, a2, precision) && point_on_segment(point, b1, b2, precision) {
150 Some(SegmentIntersection::Point(point))
151 } else {
152 None
153 }
154}
155
156fn collinear_overlap(
157 a1: Coord,
158 a2: Coord,
159 b1: Coord,
160 b2: Coord,
161 precision: PrecisionModel,
162) -> Option<SegmentIntersection> {
163 let use_x = (a2.x - a1.x).abs() >= (a2.y - a1.y).abs();
164 let mut points = [a1, a2, b1, b2];
165 points.sort_by(|left, right| {
166 let l = if use_x { left.x } else { left.y };
167 let r = if use_x { right.x } else { right.y };
168 l.partial_cmp(&r).unwrap_or(std::cmp::Ordering::Equal)
169 });
170
171 let start = points[1];
172 let end = points[2];
173 if point_on_segment(start, a1, a2, precision)
174 && point_on_segment(start, b1, b2, precision)
175 && point_on_segment(end, a1, a2, precision)
176 && point_on_segment(end, b1, b2, precision)
177 {
178 if precision.same_coord(start, end) {
179 Some(SegmentIntersection::Point(precision.snap_coord(start)))
180 } else {
181 Some(SegmentIntersection::Overlap(
182 precision.snap_coord(start),
183 precision.snap_coord(end),
184 ))
185 }
186 } else {
187 None
188 }
189}
190
191pub fn is_convex_ring(ring: &LinearRing, precision: PrecisionModel) -> bool {
192 let coords = &ring.coords;
193 if coords.len() < 4 {
194 return false;
195 }
196 let mut sign = 0i8;
197 for i in 0..coords.len() - 1 {
198 let a = coords[i];
199 let b = coords[(i + 1) % (coords.len() - 1)];
200 let c = coords[(i + 2) % (coords.len() - 1)];
201 let o = orientation(a, b, c);
202 if o.abs() <= precision.epsilon() {
203 continue;
204 }
205 let current = if o > 0.0 { 1 } else { -1 };
206 if sign == 0 {
207 sign = current;
208 } else if sign != current {
209 return false;
210 }
211 }
212 true
213}