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() > orientation_epsilon(a, b, point, precision) {
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 if segment_bboxes_disjoint(a1, a2, b1, b2, precision) {
122 return None;
123 }
124
125 let o1 = orientation(a1, a2, b1);
126 let o2 = orientation(a1, a2, b2);
127 let o3 = orientation(b1, b2, a1);
128 let o4 = orientation(b1, b2, a2);
129 let s1 = orientation_sign(o1, orientation_epsilon(a1, a2, b1, precision));
130 let s2 = orientation_sign(o2, orientation_epsilon(a1, a2, b2, precision));
131 let s3 = orientation_sign(o3, orientation_epsilon(b1, b2, a1, precision));
132 let s4 = orientation_sign(o4, orientation_epsilon(b1, b2, a2, precision));
133
134 if s1 == 0 && s2 == 0 && s3 == 0 && s4 == 0 {
135 return collinear_overlap(a1, a2, b1, b2, precision);
136 }
137
138 if s1 != 0 && s1 == s2 {
139 return None;
140 }
141 if s3 != 0 && s3 == s4 {
142 return None;
143 }
144
145 let r = Coord::new(a2.x - a1.x, a2.y - a1.y);
146 let s = Coord::new(b2.x - b1.x, b2.y - b1.y);
147 let denom = cross(r, s);
148 if denom.abs() <= cross_epsilon(r, s, precision) {
149 return None;
150 }
151
152 let q_minus_p = Coord::new(b1.x - a1.x, b1.y - a1.y);
153 let t = cross(q_minus_p, s) / denom;
154 let point = precision.snap_coord(Coord::new(a1.x + t * r.x, a1.y + t * r.y));
155
156 if point_on_segment(point, a1, a2, precision) && point_on_segment(point, b1, b2, precision) {
157 Some(SegmentIntersection::Point(point))
158 } else {
159 None
160 }
161}
162
163fn cross(a: Coord, b: Coord) -> f64 {
164 a.x * b.y - a.y * b.x
165}
166
167fn orientation_epsilon(a: Coord, b: Coord, c: Coord, precision: PrecisionModel) -> f64 {
168 let ab = l1_distance(a, b);
169 let ac = l1_distance(a, c);
170 precision.epsilon() * (ab + ac).max(f64::EPSILON)
171}
172
173fn cross_epsilon(a: Coord, b: Coord, precision: PrecisionModel) -> f64 {
174 precision.epsilon() * (l1_vector_length(a) + l1_vector_length(b)).max(f64::EPSILON)
175}
176
177fn l1_distance(a: Coord, b: Coord) -> f64 {
178 (a.x - b.x).abs() + (a.y - b.y).abs()
179}
180
181fn l1_vector_length(vector: Coord) -> f64 {
182 vector.x.abs() + vector.y.abs()
183}
184
185fn segment_bboxes_disjoint(
186 a1: Coord,
187 a2: Coord,
188 b1: Coord,
189 b2: Coord,
190 precision: PrecisionModel,
191) -> bool {
192 let eps = precision.epsilon();
193 a1.x.min(a2.x) > b1.x.max(b2.x) + eps
194 || b1.x.min(b2.x) > a1.x.max(a2.x) + eps
195 || a1.y.min(a2.y) > b1.y.max(b2.y) + eps
196 || b1.y.min(b2.y) > a1.y.max(a2.y) + eps
197}
198
199fn orientation_sign(value: f64, epsilon: f64) -> i8 {
200 if value > epsilon {
201 1
202 } else if value < -epsilon {
203 -1
204 } else {
205 0
206 }
207}
208
209fn collinear_overlap(
210 a1: Coord,
211 a2: Coord,
212 b1: Coord,
213 b2: Coord,
214 precision: PrecisionModel,
215) -> Option<SegmentIntersection> {
216 let use_x = (a2.x - a1.x).abs() >= (a2.y - a1.y).abs();
217 let mut points = [a1, a2, b1, b2];
218 points.sort_by(|left, right| {
219 let l = if use_x { left.x } else { left.y };
220 let r = if use_x { right.x } else { right.y };
221 l.partial_cmp(&r).unwrap_or(std::cmp::Ordering::Equal)
222 });
223
224 let start = points[1];
225 let end = points[2];
226 if point_on_segment(start, a1, a2, precision)
227 && point_on_segment(start, b1, b2, precision)
228 && point_on_segment(end, a1, a2, precision)
229 && point_on_segment(end, b1, b2, precision)
230 {
231 if precision.same_coord(start, end) {
232 Some(SegmentIntersection::Point(precision.snap_coord(start)))
233 } else {
234 Some(SegmentIntersection::Overlap(
235 precision.snap_coord(start),
236 precision.snap_coord(end),
237 ))
238 }
239 } else {
240 None
241 }
242}
243
244pub fn is_convex_ring(ring: &LinearRing, precision: PrecisionModel) -> bool {
245 let coords = &ring.coords;
246 if coords.len() < 4 {
247 return false;
248 }
249 let mut sign = 0i8;
250 for i in 0..coords.len() - 1 {
251 let a = coords[i];
252 let b = coords[(i + 1) % (coords.len() - 1)];
253 let c = coords[(i + 2) % (coords.len() - 1)];
254 let o = orientation(a, b, c);
255 if o.abs() <= orientation_epsilon(a, b, c, precision) {
256 continue;
257 }
258 let current = if o > 0.0 { 1 } else { -1 };
259 if sign == 0 {
260 sign = current;
261 } else if sign != current {
262 return false;
263 }
264 }
265 true
266}