1use std::collections::HashMap;
2
3use crate::canonicalize::{canonicalize_multi_polygon, canonicalize_polygon};
4use crate::error::Result;
5use crate::noding::node_lines;
6use crate::precision::PrecisionModel;
7use crate::predicates::{
8 point_in_polygon, point_in_ring, polygon_area, signed_area_coords, PointLocation,
9};
10use crate::types::{BBox, Coord, LineString, LinearRing, MultiPolygon, Polygon};
11
12pub fn intersection(
13 subject: &MultiPolygon,
14 clip: &MultiPolygon,
15 precision: PrecisionModel,
16) -> Result<MultiPolygon> {
17 overlay(subject, clip, precision, OverlayOperation::Intersection)
18}
19
20pub fn difference(
21 subject: &MultiPolygon,
22 clip: &MultiPolygon,
23 precision: PrecisionModel,
24) -> Result<MultiPolygon> {
25 overlay(subject, clip, precision, OverlayOperation::Difference)
26}
27
28#[derive(Debug, Clone, Copy, PartialEq, Eq)]
29enum OverlayOperation {
30 Intersection,
31 Difference,
32}
33
34fn overlay(
35 subject: &MultiPolygon,
36 clip: &MultiPolygon,
37 precision: PrecisionModel,
38 operation: OverlayOperation,
39) -> Result<MultiPolygon> {
40 if subject.is_empty() {
41 return Ok(MultiPolygon::empty());
42 }
43
44 if clip.is_empty() {
45 return Ok(match operation {
46 OverlayOperation::Intersection => MultiPolygon::empty(),
47 OverlayOperation::Difference => canonicalize_multi_polygon(subject, precision),
48 });
49 }
50
51 let faces = overlay_faces(subject, clip, precision);
52 let mut selected = Vec::new();
53 let mut unselected = Vec::new();
54 for face in faces {
55 let points = representative_points(&face, precision);
56 if points.is_empty() {
57 continue;
58 };
59
60 let is_selected = match operation {
61 OverlayOperation::Intersection => points.iter().any(|point| {
62 multi_polygon_contains_point(subject, *point, precision)
63 && multi_polygon_contains_point(clip, *point, precision)
64 }),
65 OverlayOperation::Difference => points.iter().any(|point| {
66 multi_polygon_contains_point(subject, *point, precision)
67 && !multi_polygon_contains_point(clip, *point, precision)
68 }),
69 };
70
71 if is_selected {
72 selected.push(face);
73 } else {
74 unselected.push(face);
75 }
76 }
77
78 Ok(assemble_selected_faces(selected, &unselected, precision))
79}
80
81fn overlay_faces(
82 subject: &MultiPolygon,
83 clip: &MultiPolygon,
84 precision: PrecisionModel,
85) -> Vec<Polygon> {
86 let lines = multi_polygon_lines(subject)
87 .into_iter()
88 .chain(multi_polygon_lines(clip))
89 .collect::<Vec<_>>();
90 let noded = node_lines(&lines, precision);
91 let arrangement = Arrangement::from_lines(&noded.lines, precision);
92 arrangement.polygonize_faces(precision)
93}
94
95fn multi_polygon_lines(multi_polygon: &MultiPolygon) -> Vec<LineString> {
96 multi_polygon
97 .polygons
98 .iter()
99 .flat_map(|polygon| {
100 std::iter::once(&polygon.exterior)
101 .chain(polygon.holes.iter())
102 .filter(|ring| ring.coords.len() >= 4)
103 .map(|ring| LineString::new(ring.coords.clone()))
104 })
105 .collect()
106}
107
108fn assemble_selected_faces(
109 selected: Vec<Polygon>,
110 unselected: &[Polygon],
111 precision: PrecisionModel,
112) -> MultiPolygon {
113 let mut output = merge_selected_faces(&selected, precision);
114 if output.is_empty() {
115 output = selected;
116 }
117 output = remove_repeated_exterior_loops(output, precision);
118 output = attach_nested_shells_as_holes(output, precision);
119
120 for hole_face in unselected {
121 let Some(point) = representative_point(hole_face, precision) else {
122 continue;
123 };
124 let Some(output_index) = containing_polygon_index(&output, point, precision) else {
125 continue;
126 };
127
128 if polygon_area(hole_face) >= polygon_area(&output[output_index]) {
129 continue;
130 }
131
132 if !ring_strictly_inside_polygon(&hole_face.exterior, &output[output_index], precision) {
133 continue;
134 }
135
136 output[output_index].holes.push(hole_face.exterior.clone());
137 }
138
139 canonicalize_multi_polygon(&MultiPolygon::new(output), precision)
140}
141
142fn remove_repeated_exterior_loops(
143 polygons: Vec<Polygon>,
144 precision: PrecisionModel,
145) -> Vec<Polygon> {
146 polygons
147 .into_iter()
148 .map(|mut polygon| {
149 polygon.exterior = remove_same_orientation_repeated_loops(&polygon.exterior, precision);
150 polygon
151 })
152 .collect()
153}
154
155fn remove_same_orientation_repeated_loops(
156 ring: &LinearRing,
157 precision: PrecisionModel,
158) -> LinearRing {
159 let mut out = Vec::<Coord>::new();
160 let mut indexes = HashMap::<CoordKey, usize>::new();
161
162 for coord in ring.coords.iter().copied() {
163 let coord = precision.snap_coord(coord);
164 if out
165 .last()
166 .is_some_and(|previous| precision.same_coord(*previous, coord))
167 {
168 continue;
169 }
170 if !out.is_empty() && precision.same_coord(out[0], coord) {
171 continue;
172 }
173
174 let key = CoordKey::new(coord);
175 if let Some(previous_index) = indexes.get(&key).copied() {
176 let mut loop_coords = out[previous_index..].to_vec();
177 loop_coords.push(out[previous_index]);
178 if signed_area_coords(&loop_coords) > precision.epsilon() {
179 for removed in out.drain(previous_index + 1..) {
180 indexes.remove(&CoordKey::new(removed));
181 }
182 continue;
183 }
184 }
185
186 indexes.insert(key, out.len());
187 out.push(coord);
188 }
189
190 if out.len() < 3 {
191 return LinearRing { coords: Vec::new() };
192 }
193
194 out.push(out[0]);
195 LinearRing::new(out)
196}
197
198fn ring_strictly_inside_polygon(
199 ring: &LinearRing,
200 polygon: &Polygon,
201 precision: PrecisionModel,
202) -> bool {
203 ring.coords
204 .iter()
205 .take(ring.coords.len().saturating_sub(1))
206 .all(|coord| {
207 matches!(
208 point_in_polygon(*coord, polygon, precision),
209 PointLocation::Interior
210 )
211 })
212}
213
214fn attach_nested_shells_as_holes(
215 mut polygons: Vec<Polygon>,
216 precision: PrecisionModel,
217) -> Vec<Polygon> {
218 if polygons.len() <= 1 {
219 return polygons;
220 }
221
222 let mut remove = vec![false; polygons.len()];
223 let mut hole_assignments = Vec::<(usize, LinearRing)>::new();
224
225 for inner_index in 0..polygons.len() {
226 let inner_area = polygon_area(&polygons[inner_index]);
227 let Some(point) = representative_point(&polygons[inner_index], precision) else {
228 continue;
229 };
230
231 let Some(outer_index) = polygons
232 .iter()
233 .enumerate()
234 .filter(|(outer_index, outer)| {
235 *outer_index != inner_index
236 && polygon_area(outer) > inner_area
237 && matches!(
238 point_in_ring(point, &outer.exterior, precision),
239 PointLocation::Interior
240 )
241 })
242 .min_by(|(_, left), (_, right)| {
243 polygon_area(left)
244 .partial_cmp(&polygon_area(right))
245 .unwrap_or(std::cmp::Ordering::Equal)
246 })
247 .map(|(outer_index, _)| outer_index)
248 else {
249 continue;
250 };
251
252 remove[inner_index] = true;
253 hole_assignments.push((outer_index, polygons[inner_index].exterior.clone()));
254 }
255
256 for (outer_index, hole) in hole_assignments {
257 if !remove[outer_index] {
258 polygons[outer_index].holes.push(hole);
259 }
260 }
261
262 polygons
263 .into_iter()
264 .enumerate()
265 .filter_map(|(index, polygon)| (!remove[index]).then_some(polygon))
266 .collect()
267}
268
269fn merge_selected_faces(selected: &[Polygon], precision: PrecisionModel) -> Vec<Polygon> {
270 if selected.len() <= 1 {
271 return selected.to_vec();
272 }
273
274 let mut edge_counts = HashMap::<EdgeKey, BoundaryEdge>::new();
275 for face in selected {
276 for (start, end) in face.exterior.segments() {
277 let key = EdgeKey::new(start, end);
278 edge_counts
279 .entry(key)
280 .and_modify(|edge| edge.count += 1)
281 .or_insert(BoundaryEdge {
282 start,
283 end,
284 count: 1,
285 });
286 }
287 }
288
289 let boundary_lines = edge_counts
290 .into_values()
291 .filter(|edge| edge.count == 1)
292 .map(|edge| LineString::new(vec![edge.start, edge.end]))
293 .collect::<Vec<_>>();
294
295 if boundary_lines.is_empty() {
296 return Vec::new();
297 }
298
299 Arrangement::from_lines(&boundary_lines, precision).polygonize_faces(precision)
300}
301
302#[derive(Debug, Clone, Copy)]
303struct BoundaryEdge {
304 start: Coord,
305 end: Coord,
306 count: usize,
307}
308
309#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
310struct EdgeKey(CoordKey, CoordKey);
311
312impl EdgeKey {
313 fn new(start: Coord, end: Coord) -> Self {
314 let start = CoordKey::new(start);
315 let end = CoordKey::new(end);
316 if start <= end {
317 Self(start, end)
318 } else {
319 Self(end, start)
320 }
321 }
322}
323
324#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
325struct CoordKey(u64, u64);
326
327impl CoordKey {
328 fn new(coord: Coord) -> Self {
329 Self(
330 normalize_zero(coord.x).to_bits(),
331 normalize_zero(coord.y).to_bits(),
332 )
333 }
334}
335
336fn containing_polygon_index(
337 polygons: &[Polygon],
338 point: Coord,
339 precision: PrecisionModel,
340) -> Option<usize> {
341 polygons
342 .iter()
343 .enumerate()
344 .filter(|(_, polygon)| {
345 matches!(
346 point_in_polygon(point, polygon, precision),
347 PointLocation::Interior
348 )
349 })
350 .min_by(|(_, left), (_, right)| {
351 polygon_area(left)
352 .partial_cmp(&polygon_area(right))
353 .unwrap_or(std::cmp::Ordering::Equal)
354 })
355 .map(|(index, _)| index)
356}
357
358fn multi_polygon_contains_point(
359 multi_polygon: &MultiPolygon,
360 point: Coord,
361 precision: PrecisionModel,
362) -> bool {
363 multi_polygon
364 .polygons
365 .iter()
366 .any(|polygon| polygon_contains_point(polygon, point, precision))
367}
368
369fn polygon_bbox(polygon: &Polygon) -> Option<BBox> {
370 BBox::from_coords(&polygon.exterior.coords)
371}
372
373#[derive(Debug, Clone)]
374struct Arrangement {
375 vertices: Vec<Coord>,
376 adjacency: Vec<Vec<usize>>,
377}
378
379impl Arrangement {
380 fn from_lines(lines: &[LineString], precision: PrecisionModel) -> Self {
381 let mut vertices = Vec::new();
382 let mut adjacency: Vec<Vec<usize>> = Vec::new();
383
384 for line in lines {
385 for (start, end) in line.segments() {
386 let start_index = vertex_index(&mut vertices, &mut adjacency, start, precision);
387 let end_index = vertex_index(&mut vertices, &mut adjacency, end, precision);
388 if start_index == end_index {
389 continue;
390 }
391
392 push_unique_neighbor(&mut adjacency[start_index], end_index);
393 push_unique_neighbor(&mut adjacency[end_index], start_index);
394 }
395 }
396
397 for index in 0..vertices.len() {
398 let origin = vertices[index];
399 adjacency[index].sort_by(|left, right| {
400 edge_angle(origin, vertices[*left])
401 .partial_cmp(&edge_angle(origin, vertices[*right]))
402 .unwrap_or(std::cmp::Ordering::Equal)
403 });
404 }
405
406 Self {
407 vertices,
408 adjacency,
409 }
410 }
411
412 fn polygonize_faces(&self, precision: PrecisionModel) -> Vec<Polygon> {
413 let mut visited = Vec::<(usize, usize)>::new();
414 let mut polygons = Vec::new();
415 let mut reversed_polygons = Vec::new();
416
417 for start in 0..self.vertices.len() {
418 for &end in &self.adjacency[start] {
419 if directed_edge_seen(&visited, start, end) {
420 continue;
421 }
422
423 let ring = self.walk_face(start, end, &mut visited);
424 if ring.len() < 3 {
425 continue;
426 }
427
428 let mut coords = ring
429 .iter()
430 .map(|index| self.vertices[*index])
431 .collect::<Vec<_>>();
432 coords.push(coords[0]);
433
434 let signed_area = signed_area_coords(&coords);
435 if signed_area.abs() <= precision.epsilon() {
436 continue;
437 }
438
439 if signed_area < 0.0 {
440 coords.reverse();
441 let polygon = canonicalize_polygon(
442 &Polygon::new(LinearRing::new(coords), Vec::new()),
443 precision,
444 );
445 if !polygon.is_empty()
446 && !reversed_polygons
447 .iter()
448 .any(|existing| existing == &polygon)
449 {
450 reversed_polygons.push(polygon);
451 }
452 continue;
453 }
454
455 let polygon = canonicalize_polygon(
456 &Polygon::new(LinearRing::new(coords), Vec::new()),
457 precision,
458 );
459 if !polygon.is_empty() && !polygons.iter().any(|existing| existing == &polygon) {
460 polygons.push(polygon);
461 }
462 }
463 }
464
465 if polygons.is_empty() {
466 reversed_polygons
467 } else {
468 polygons
469 }
470 }
471
472 fn walk_face(&self, start: usize, end: usize, visited: &mut Vec<(usize, usize)>) -> Vec<usize> {
473 let mut ring = Vec::new();
474 let mut current_start = start;
475 let mut current_end = end;
476 let max_steps = self
477 .adjacency
478 .iter()
479 .map(Vec::len)
480 .sum::<usize>()
481 .saturating_add(1);
482
483 for _ in 0..max_steps {
484 if directed_edge_seen(visited, current_start, current_end) {
485 break;
486 }
487
488 visited.push((current_start, current_end));
489 ring.push(current_start);
490
491 let Some(next) = self.next_face_vertex(current_start, current_end) else {
492 break;
493 };
494 current_start = current_end;
495 current_end = next;
496
497 if current_start == start && current_end == end {
498 break;
499 }
500 }
501
502 ring
503 }
504
505 fn next_face_vertex(&self, previous: usize, current: usize) -> Option<usize> {
506 let neighbors = self.adjacency.get(current)?;
507 let reverse_index = neighbors
508 .iter()
509 .position(|neighbor| *neighbor == previous)?;
510 Some(neighbors[(reverse_index + neighbors.len() - 1) % neighbors.len()])
511 }
512}
513
514fn vertex_index(
515 vertices: &mut Vec<Coord>,
516 adjacency: &mut Vec<Vec<usize>>,
517 coord: Coord,
518 precision: PrecisionModel,
519) -> usize {
520 if let Some(index) = vertices
521 .iter()
522 .position(|existing| precision.same_coord(*existing, coord))
523 {
524 return index;
525 }
526
527 vertices.push(precision.snap_coord(coord));
528 adjacency.push(Vec::new());
529 vertices.len() - 1
530}
531
532fn push_unique_neighbor(neighbors: &mut Vec<usize>, neighbor: usize) {
533 if !neighbors.contains(&neighbor) {
534 neighbors.push(neighbor);
535 }
536}
537
538fn edge_angle(origin: Coord, target: Coord) -> f64 {
539 (target.y - origin.y).atan2(target.x - origin.x)
540}
541
542fn directed_edge_seen(visited: &[(usize, usize)], start: usize, end: usize) -> bool {
543 visited.iter().any(|edge| edge.0 == start && edge.1 == end)
544}
545
546fn normalize_zero(value: f64) -> f64 {
547 if value == 0.0 {
548 0.0
549 } else {
550 value
551 }
552}
553
554fn representative_point(polygon: &Polygon, precision: PrecisionModel) -> Option<Coord> {
555 representative_points(polygon, precision).into_iter().next()
556}
557
558fn representative_points(polygon: &Polygon, precision: PrecisionModel) -> Vec<Coord> {
559 let mut points = Vec::new();
560
561 if let Some(point) = triangle_fan_representative_point(polygon, precision) {
562 push_unique_point(&mut points, point, precision);
563 }
564
565 if let Some(point) = polygon_centroid(polygon) {
566 if matches!(
567 point_in_polygon(point, polygon, precision),
568 PointLocation::Interior
569 ) {
570 push_unique_point(&mut points, point, precision);
571 }
572 }
573
574 if let Some(bbox) = polygon_bbox(polygon) {
575 let point = Coord::new(
576 (bbox.min.x + bbox.max.x) * 0.5,
577 (bbox.min.y + bbox.max.y) * 0.5,
578 );
579 if matches!(
580 point_in_polygon(point, polygon, precision),
581 PointLocation::Interior
582 ) {
583 push_unique_point(&mut points, point, precision);
584 }
585 }
586
587 for point in edge_probe_representative_points(polygon, precision) {
588 push_unique_point(&mut points, point, precision);
589 }
590
591 points
592}
593
594fn edge_probe_representative_points(polygon: &Polygon, precision: PrecisionModel) -> Vec<Coord> {
595 let mut points = Vec::new();
596 let extent = polygon_bbox(polygon)
597 .map(|bbox| {
598 let width = bbox.max.x - bbox.min.x;
599 let height = bbox.max.y - bbox.min.y;
600 (width * width + height * height).sqrt()
601 })
602 .unwrap_or(1.0);
603 let offsets = [
604 (extent * 1.0e-9).max(precision.epsilon() * 10.0),
605 (extent * 1.0e-7).max(precision.epsilon() * 10.0),
606 (extent * 1.0e-5).max(precision.epsilon() * 10.0),
607 ];
608
609 for (start, end) in polygon.exterior.segments() {
610 let dx = end.x - start.x;
611 let dy = end.y - start.y;
612 let length = (dx * dx + dy * dy).sqrt();
613 if length <= precision.epsilon() {
614 continue;
615 }
616
617 let midpoint = Coord::new((start.x + end.x) * 0.5, (start.y + end.y) * 0.5);
618 let normal = Coord::new(-dy / length, dx / length);
619 for offset in offsets {
620 for direction in [1.0, -1.0] {
621 let point = Coord::new(
622 midpoint.x + normal.x * offset * direction,
623 midpoint.y + normal.y * offset * direction,
624 );
625 if matches!(
626 point_in_polygon(point, polygon, precision),
627 PointLocation::Interior
628 ) {
629 push_unique_point(&mut points, point, precision);
630 }
631 }
632 }
633 }
634
635 points
636}
637
638fn push_unique_point(points: &mut Vec<Coord>, point: Coord, precision: PrecisionModel) {
639 if points
640 .iter()
641 .all(|existing| !precision.same_coord(*existing, point))
642 {
643 points.push(point);
644 }
645}
646
647fn polygon_centroid(polygon: &Polygon) -> Option<Coord> {
648 let coords = &polygon.exterior.coords;
649 if coords.len() < 4 {
650 return None;
651 }
652
653 let mut twice_area = 0.0;
654 let mut centroid_x = 0.0;
655 let mut centroid_y = 0.0;
656 for pair in coords.windows(2) {
657 let cross = pair[0].x * pair[1].y - pair[1].x * pair[0].y;
658 twice_area += cross;
659 centroid_x += (pair[0].x + pair[1].x) * cross;
660 centroid_y += (pair[0].y + pair[1].y) * cross;
661 }
662
663 if twice_area.abs() <= f64::EPSILON {
664 return None;
665 }
666
667 Some(Coord::new(
668 centroid_x / (3.0 * twice_area),
669 centroid_y / (3.0 * twice_area),
670 ))
671}
672
673fn triangle_fan_representative_point(
674 polygon: &Polygon,
675 precision: PrecisionModel,
676) -> Option<Coord> {
677 let coords = &polygon.exterior.coords;
678 if coords.len() < 4 {
679 return None;
680 }
681
682 let anchor = coords[0];
683 for pair in coords[1..coords.len() - 1].windows(2) {
684 let centroid = Coord::new(
685 (anchor.x + pair[0].x + pair[1].x) / 3.0,
686 (anchor.y + pair[0].y + pair[1].y) / 3.0,
687 );
688 if matches!(
689 point_in_polygon(centroid, polygon, precision),
690 PointLocation::Interior
691 ) {
692 return Some(centroid);
693 }
694 }
695
696 None
697}
698
699fn polygon_contains_point(polygon: &Polygon, point: Coord, precision: PrecisionModel) -> bool {
700 matches!(
701 point_in_polygon(point, polygon, precision),
702 PointLocation::Interior | PointLocation::Boundary
703 )
704}