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