1use std::collections::{HashMap, HashSet};
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 vertex_indexes = HashMap::<CoordKey, usize>::new();
383 let mut adjacency: Vec<Vec<usize>> = Vec::new();
384
385 for line in lines {
386 for (start, end) in line.segments() {
387 let start_index = vertex_index(
388 &mut vertices,
389 &mut vertex_indexes,
390 &mut adjacency,
391 start,
392 precision,
393 );
394 let end_index = vertex_index(
395 &mut vertices,
396 &mut vertex_indexes,
397 &mut adjacency,
398 end,
399 precision,
400 );
401 if start_index == end_index {
402 continue;
403 }
404
405 push_unique_neighbor(&mut adjacency[start_index], end_index);
406 push_unique_neighbor(&mut adjacency[end_index], start_index);
407 }
408 }
409
410 for index in 0..vertices.len() {
411 let origin = vertices[index];
412 adjacency[index].sort_by(|left, right| {
413 edge_angle(origin, vertices[*left])
414 .partial_cmp(&edge_angle(origin, vertices[*right]))
415 .unwrap_or(std::cmp::Ordering::Equal)
416 });
417 }
418
419 Self {
420 vertices,
421 adjacency,
422 }
423 }
424
425 fn polygonize_faces(&self, precision: PrecisionModel) -> Vec<Polygon> {
426 let mut visited = HashSet::<(usize, usize)>::new();
427 let mut polygons = Vec::new();
428 let mut reversed_polygons = Vec::new();
429
430 for start in 0..self.vertices.len() {
431 for &end in &self.adjacency[start] {
432 if directed_edge_seen(&visited, start, end) {
433 continue;
434 }
435
436 let ring = self.walk_face(start, end, &mut visited);
437 if ring.len() < 3 {
438 continue;
439 }
440
441 let mut coords = ring
442 .iter()
443 .map(|index| self.vertices[*index])
444 .collect::<Vec<_>>();
445 coords.push(coords[0]);
446
447 let signed_area = signed_area_coords(&coords);
448 if signed_area.abs() <= precision.epsilon() {
449 continue;
450 }
451
452 if signed_area < 0.0 {
453 coords.reverse();
454 let polygon = canonicalize_polygon(
455 &Polygon::new(LinearRing::new(coords), Vec::new()),
456 precision,
457 );
458 if !polygon.is_empty()
459 && !reversed_polygons
460 .iter()
461 .any(|existing| existing == &polygon)
462 {
463 reversed_polygons.push(polygon);
464 }
465 continue;
466 }
467
468 let polygon = canonicalize_polygon(
469 &Polygon::new(LinearRing::new(coords), Vec::new()),
470 precision,
471 );
472 if !polygon.is_empty() && !polygons.iter().any(|existing| existing == &polygon) {
473 polygons.push(polygon);
474 }
475 }
476 }
477
478 if polygons.is_empty() {
479 reversed_polygons
480 } else {
481 polygons
482 }
483 }
484
485 fn walk_face(
486 &self,
487 start: usize,
488 end: usize,
489 visited: &mut HashSet<(usize, usize)>,
490 ) -> Vec<usize> {
491 let mut ring = Vec::new();
492 let mut current_start = start;
493 let mut current_end = end;
494 let max_steps = self
495 .adjacency
496 .iter()
497 .map(Vec::len)
498 .sum::<usize>()
499 .saturating_add(1);
500
501 for _ in 0..max_steps {
502 if directed_edge_seen(visited, current_start, current_end) {
503 break;
504 }
505
506 visited.insert((current_start, current_end));
507 ring.push(current_start);
508
509 let Some(next) = self.next_face_vertex(current_start, current_end) else {
510 break;
511 };
512 current_start = current_end;
513 current_end = next;
514
515 if current_start == start && current_end == end {
516 break;
517 }
518 }
519
520 ring
521 }
522
523 fn next_face_vertex(&self, previous: usize, current: usize) -> Option<usize> {
524 let neighbors = self.adjacency.get(current)?;
525 let reverse_index = neighbors
526 .iter()
527 .position(|neighbor| *neighbor == previous)?;
528 Some(neighbors[(reverse_index + neighbors.len() - 1) % neighbors.len()])
529 }
530}
531
532fn vertex_index(
533 vertices: &mut Vec<Coord>,
534 vertex_indexes: &mut HashMap<CoordKey, usize>,
535 adjacency: &mut Vec<Vec<usize>>,
536 coord: Coord,
537 precision: PrecisionModel,
538) -> usize {
539 let snapped = precision.snap_coord(coord);
540 let key = CoordKey::new(snapped);
541 if let Some(index) = vertex_indexes.get(&key).copied() {
542 return index;
543 }
544
545 if let Some(index) = vertices
546 .iter()
547 .position(|existing| precision.same_coord(*existing, snapped))
548 {
549 vertex_indexes.insert(key, index);
550 return index;
551 }
552
553 vertices.push(snapped);
554 adjacency.push(Vec::new());
555 vertex_indexes.insert(key, vertices.len() - 1);
556 vertices.len() - 1
557}
558
559fn push_unique_neighbor(neighbors: &mut Vec<usize>, neighbor: usize) {
560 if !neighbors.contains(&neighbor) {
561 neighbors.push(neighbor);
562 }
563}
564
565fn edge_angle(origin: Coord, target: Coord) -> f64 {
566 (target.y - origin.y).atan2(target.x - origin.x)
567}
568
569fn directed_edge_seen(visited: &HashSet<(usize, usize)>, start: usize, end: usize) -> bool {
570 visited.contains(&(start, end))
571}
572
573fn normalize_zero(value: f64) -> f64 {
574 if value == 0.0 {
575 0.0
576 } else {
577 value
578 }
579}
580
581fn representative_point(polygon: &Polygon, precision: PrecisionModel) -> Option<Coord> {
582 representative_points(polygon, precision).into_iter().next()
583}
584
585fn representative_points(polygon: &Polygon, precision: PrecisionModel) -> Vec<Coord> {
586 let mut points = Vec::new();
587
588 if let Some(point) = triangle_fan_representative_point(polygon, precision) {
589 push_unique_point(&mut points, point, precision);
590 }
591
592 if let Some(point) = polygon_centroid(polygon) {
593 if matches!(
594 point_in_polygon(point, polygon, precision),
595 PointLocation::Interior
596 ) {
597 push_unique_point(&mut points, point, precision);
598 }
599 }
600
601 if let Some(bbox) = polygon_bbox(polygon) {
602 let point = Coord::new(
603 (bbox.min.x + bbox.max.x) * 0.5,
604 (bbox.min.y + bbox.max.y) * 0.5,
605 );
606 if matches!(
607 point_in_polygon(point, polygon, precision),
608 PointLocation::Interior
609 ) {
610 push_unique_point(&mut points, point, precision);
611 }
612 }
613
614 for point in edge_probe_representative_points(polygon, precision) {
615 push_unique_point(&mut points, point, precision);
616 }
617
618 points
619}
620
621fn edge_probe_representative_points(polygon: &Polygon, precision: PrecisionModel) -> Vec<Coord> {
622 let mut points = Vec::new();
623 let extent = polygon_bbox(polygon)
624 .map(|bbox| {
625 let width = bbox.max.x - bbox.min.x;
626 let height = bbox.max.y - bbox.min.y;
627 (width * width + height * height).sqrt()
628 })
629 .unwrap_or(1.0);
630 let offsets = [
631 (extent * 1.0e-9).max(precision.epsilon() * 10.0),
632 (extent * 1.0e-7).max(precision.epsilon() * 10.0),
633 (extent * 1.0e-5).max(precision.epsilon() * 10.0),
634 ];
635
636 for (start, end) in polygon.exterior.segments() {
637 let dx = end.x - start.x;
638 let dy = end.y - start.y;
639 let length = (dx * dx + dy * dy).sqrt();
640 if length <= precision.epsilon() {
641 continue;
642 }
643
644 let midpoint = Coord::new((start.x + end.x) * 0.5, (start.y + end.y) * 0.5);
645 let normal = Coord::new(-dy / length, dx / length);
646 for offset in offsets {
647 for direction in [1.0, -1.0] {
648 let point = Coord::new(
649 midpoint.x + normal.x * offset * direction,
650 midpoint.y + normal.y * offset * direction,
651 );
652 if matches!(
653 point_in_polygon(point, polygon, precision),
654 PointLocation::Interior
655 ) {
656 push_unique_point(&mut points, point, precision);
657 }
658 }
659 }
660 }
661
662 points
663}
664
665fn push_unique_point(points: &mut Vec<Coord>, point: Coord, precision: PrecisionModel) {
666 if points
667 .iter()
668 .all(|existing| !precision.same_coord(*existing, point))
669 {
670 points.push(point);
671 }
672}
673
674fn polygon_centroid(polygon: &Polygon) -> Option<Coord> {
675 let coords = &polygon.exterior.coords;
676 if coords.len() < 4 {
677 return None;
678 }
679
680 let mut twice_area = 0.0;
681 let mut centroid_x = 0.0;
682 let mut centroid_y = 0.0;
683 for pair in coords.windows(2) {
684 let cross = pair[0].x * pair[1].y - pair[1].x * pair[0].y;
685 twice_area += cross;
686 centroid_x += (pair[0].x + pair[1].x) * cross;
687 centroid_y += (pair[0].y + pair[1].y) * cross;
688 }
689
690 if twice_area.abs() <= f64::EPSILON {
691 return None;
692 }
693
694 Some(Coord::new(
695 centroid_x / (3.0 * twice_area),
696 centroid_y / (3.0 * twice_area),
697 ))
698}
699
700fn triangle_fan_representative_point(
701 polygon: &Polygon,
702 precision: PrecisionModel,
703) -> Option<Coord> {
704 let coords = &polygon.exterior.coords;
705 if coords.len() < 4 {
706 return None;
707 }
708
709 let anchor = coords[0];
710 for pair in coords[1..coords.len() - 1].windows(2) {
711 let centroid = Coord::new(
712 (anchor.x + pair[0].x + pair[1].x) / 3.0,
713 (anchor.y + pair[0].y + pair[1].y) / 3.0,
714 );
715 if matches!(
716 point_in_polygon(centroid, polygon, precision),
717 PointLocation::Interior
718 ) {
719 return Some(centroid);
720 }
721 }
722
723 None
724}
725
726fn polygon_contains_point(polygon: &Polygon, point: Coord, precision: PrecisionModel) -> bool {
727 matches!(
728 point_in_polygon(point, polygon, precision),
729 PointLocation::Interior | PointLocation::Boundary
730 )
731}