1use std::collections::{HashMap, HashSet};
2use std::f64::consts::PI;
3
4use crate::canonicalize::canonicalize_polygon;
5use crate::error::{GeometryError, Result};
6use crate::noding::node_lines;
7use crate::precision::PrecisionModel;
8use crate::predicates::{
9 is_ring_ccw, orientation, point_in_ring, segment_intersection, signed_ring_area, PointLocation,
10 SegmentIntersection,
11};
12use crate::types::{Coord, LineString, LinearRing, MultiPolygon, Polygon};
13
14#[derive(Debug, Clone, Copy, PartialEq, Eq)]
15pub enum JoinStyle {
16 Round,
17 Mitre,
18 Bevel,
19}
20
21#[derive(Debug, Clone, Copy, PartialEq, Eq)]
22pub enum CapStyle {
23 Round,
24 Flat,
25 Square,
26}
27
28#[derive(Debug, Clone, Copy, PartialEq)]
29pub struct BufferOptions {
30 pub quadrant_segments: u32,
31 pub join_style: JoinStyle,
32 pub cap_style: CapStyle,
33 pub mitre_limit: f64,
34 pub simplify_factor: f64,
35}
36
37impl Default for BufferOptions {
38 fn default() -> Self {
39 Self {
40 quadrant_segments: 8,
41 join_style: JoinStyle::Round,
42 cap_style: CapStyle::Round,
43 mitre_limit: 5.0,
44 simplify_factor: 0.01,
45 }
46 }
47}
48
49pub fn buffer_polygon(
50 polygon: &Polygon,
51 distance: f64,
52 options: BufferOptions,
53 precision: PrecisionModel,
54) -> Result<MultiPolygon> {
55 if polygon.is_empty() || distance == 0.0 {
56 return Ok(MultiPolygon::new(vec![canonicalize_polygon(
57 polygon, precision,
58 )]));
59 }
60
61 let polygon = precision.snap_polygon(polygon);
62 let shell_coords = remove_repeated_points(&polygon.exterior.coords, precision);
63 if shell_coords.len() < 4 {
64 return Ok(MultiPolygon::empty());
65 }
66 if distance < 0.0 && is_eroded_completely(&shell_coords, distance) {
67 return Ok(MultiPolygon::empty());
68 }
69
70 let offset_distance = distance.abs();
71 let shell_side = if distance < 0.0 {
72 Position::Right
73 } else {
74 Position::Left
75 };
76 let exterior_curve = ring_curve(
77 &shell_coords,
78 shell_side,
79 offset_distance,
80 options,
81 precision,
82 )?;
83
84 let mut curves = vec![exterior_curve.clone()];
85 let mut direct_holes = Vec::new();
86 for hole in &polygon.holes {
87 let hole_coords = remove_repeated_points(&hole.coords, precision);
88 if hole_coords.len() < 4 {
89 continue;
90 }
91 if distance > 0.0 && is_eroded_completely(&hole_coords, -distance) {
92 continue;
93 }
94 let hole_curve = ring_curve(
95 &hole_coords,
96 shell_side.opposite(),
97 offset_distance,
98 options,
99 precision,
100 )?;
101 direct_holes.push(LinearRing::new(hole_curve.clone()));
102 curves.push(hole_curve);
103 }
104
105 let polygonized = polygonize_buffer_curves(&curves, precision);
106 if !polygonized.is_empty() {
107 return Ok(MultiPolygon::new(polygonized));
108 }
109
110 let direct = canonicalize_polygon(
111 &Polygon::new(LinearRing::new(exterior_curve), direct_holes),
112 precision,
113 );
114 if direct.is_empty() {
115 Ok(MultiPolygon::empty())
116 } else {
117 Ok(MultiPolygon::new(vec![direct]))
118 }
119}
120
121pub fn line_buffer(
122 line: &LineString,
123 distance: f64,
124 options: BufferOptions,
125 precision: PrecisionModel,
126) -> Result<MultiPolygon> {
127 if line.coords.len() < 2 || distance <= 0.0 {
128 return Ok(MultiPolygon::empty());
129 }
130
131 let snapped = precision.snap_line(line);
132 let coords = remove_repeated_points(&snapped.coords, precision);
133 if coords.len() < 2 {
134 let Some(center) = snapped.coords.first().copied() else {
135 return Ok(MultiPolygon::empty());
136 };
137 let polygon = canonicalize_polygon(
138 &Polygon::new(
139 LinearRing::new(point_buffer_ring(center, distance, options, precision)),
140 Vec::new(),
141 ),
142 precision,
143 );
144 return if polygon.is_empty() {
145 Ok(MultiPolygon::empty())
146 } else {
147 Ok(MultiPolygon::new(vec![polygon]))
148 };
149 }
150
151 let coords = if coords.len() == 2 && options.join_style == JoinStyle::Round {
152 single_segment_line_buffer_ring(&coords, distance, options, precision)?
153 } else {
154 line_buffer_curve(&coords, distance, options, precision)?
155 };
156 let polygonized = polygonize_buffer_curves(std::slice::from_ref(&coords), precision);
157 if !polygonized.is_empty() {
158 return Ok(MultiPolygon::new(polygonized));
159 }
160
161 let polygon = canonicalize_polygon(
162 &Polygon::new(LinearRing::new(coords), Vec::new()),
163 precision,
164 );
165 if polygon.is_empty() {
166 Ok(MultiPolygon::empty())
167 } else {
168 Ok(MultiPolygon::new(vec![polygon]))
169 }
170}
171
172#[derive(Debug, Clone, Copy, PartialEq, Eq)]
173enum Position {
174 Left,
175 Right,
176}
177
178impl Position {
179 const fn opposite(self) -> Self {
180 match self {
181 Self::Left => Self::Right,
182 Self::Right => Self::Left,
183 }
184 }
185}
186
187const CLOCKWISE: i8 = -1;
188const COUNTERCLOCKWISE: i8 = 1;
189const COLLINEAR: i8 = 0;
190const OFFSET_SEGMENT_SEPARATION_FACTOR: f64 = 1.0e-3;
191const INSIDE_TURN_VERTEX_SNAP_DISTANCE_FACTOR: f64 = 1.0e-3;
192const CURVE_VERTEX_SNAP_DISTANCE_FACTOR: f64 = 1.0e-6;
193const MAX_CLOSING_SEG_LEN_FACTOR: f64 = 80.0;
194const NUM_SIMPLIFY_PTS_TO_CHECK: usize = 10;
195
196fn ring_curve(
197 coords: &[Coord],
198 mut side: Position,
199 distance: f64,
200 options: BufferOptions,
201 precision: PrecisionModel,
202) -> Result<Vec<Coord>> {
203 if coords.len() <= 2 {
204 return line_buffer_curve(coords, distance, options, precision);
205 }
206
207 if coords.len() >= 4 && is_ccw_jts(coords) {
208 side = side.opposite();
209 }
210
211 let simplify_tolerance = simplify_tolerance(distance, options, side == Position::Right);
212 let simplified = simplify_buffer_input_line(coords, simplify_tolerance);
213 if simplified.len() < 4 {
214 return Err(GeometryError::InvalidGeometry(
215 "buffer simplification produced an invalid ring".to_owned(),
216 ));
217 }
218
219 let n = simplified.len() - 1;
220 let mut generator = OffsetSegmentGenerator::new(distance, options, precision);
221 generator.init_side_segments(simplified[n - 1], simplified[0], side)?;
222 for i in 1..=n {
223 generator.add_next_segment(simplified[i], i != 1)?;
224 }
225 generator.close_ring();
226 Ok(generator.coordinates())
227}
228
229fn line_buffer_curve(
230 coords: &[Coord],
231 distance: f64,
232 options: BufferOptions,
233 precision: PrecisionModel,
234) -> Result<Vec<Coord>> {
235 let dist_tol = simplify_tolerance(distance, options, false);
236 let simp1 = simplify_buffer_input_line(coords, dist_tol);
237 if simp1.len() < 2 {
238 return Ok(Vec::new());
239 }
240
241 let mut generator = OffsetSegmentGenerator::new(distance, options, precision);
242 let n1 = simp1.len() - 1;
243 generator.init_side_segments(simp1[0], simp1[1], Position::Left)?;
244 for point in simp1.iter().take(n1 + 1).skip(2) {
245 generator.add_next_segment(*point, true)?;
246 }
247 generator.add_last_segment();
248 generator.add_line_end_cap(simp1[n1 - 1], simp1[n1])?;
249
250 let simp2 = simplify_buffer_input_line(coords, -dist_tol);
251 let n2 = simp2.len() - 1;
252 generator.init_side_segments(simp2[n2], simp2[n2 - 1], Position::Left)?;
253 if n2 >= 2 {
254 for i in (0..=n2 - 2).rev() {
255 generator.add_next_segment(simp2[i], true)?;
256 }
257 }
258 generator.add_last_segment();
259 generator.add_line_end_cap(simp2[1], simp2[0])?;
260 generator.close_ring();
261
262 Ok(generator.coordinates())
263}
264
265fn single_segment_line_buffer_ring(
266 coords: &[Coord],
267 distance: f64,
268 options: BufferOptions,
269 precision: PrecisionModel,
270) -> Result<Vec<Coord>> {
271 let start = coords[0];
272 let end = coords[1];
273 let dx = end.x - start.x;
274 let dy = end.y - start.y;
275 let length = (dx * dx + dy * dy).sqrt();
276 if length <= precision.epsilon() {
277 return Ok(Vec::new());
278 }
279
280 let unit = Coord::new(dx / length, dy / length);
281 let normal = Coord::new(-unit.y, unit.x);
282 let cap_steps = options.quadrant_segments.max(1) as usize * 2;
283 let angle_step = PI / cap_steps as f64;
284 let mut ring = Vec::with_capacity(2 * cap_steps + 3);
285
286 ring.push(offset_point(start, normal, distance, precision));
287 ring.push(offset_point(end, normal, distance, precision));
288
289 for step in 1..=cap_steps {
290 let rotated = rotate_vector(normal, -angle_step * step as f64);
291 ring.push(offset_point(end, rotated, distance, precision));
292 }
293
294 let opposite = Coord::new(-normal.x, -normal.y);
295 ring.push(offset_point(start, opposite, distance, precision));
296
297 for step in 1..=cap_steps {
298 let rotated = rotate_vector(opposite, -angle_step * step as f64);
299 ring.push(offset_point(start, rotated, distance, precision));
300 }
301
302 Ok(ring)
303}
304
305fn point_buffer_ring(
306 center: Coord,
307 distance: f64,
308 options: BufferOptions,
309 precision: PrecisionModel,
310) -> Vec<Coord> {
311 let steps = options.quadrant_segments.max(1) as usize * 4;
312 (0..=steps)
313 .map(|step| {
314 let angle = 2.0 * PI * step as f64 / steps as f64;
315 precision.snap_coord(Coord::new(
316 center.x + distance * angle.cos(),
317 center.y + distance * angle.sin(),
318 ))
319 })
320 .collect()
321}
322
323fn offset_point(point: Coord, unit: Coord, distance: f64, precision: PrecisionModel) -> Coord {
324 precision.snap_coord(Coord::new(
325 point.x + unit.x * distance,
326 point.y + unit.y * distance,
327 ))
328}
329
330fn rotate_vector(vector: Coord, angle: f64) -> Coord {
331 let cos = angle.cos();
332 let sin = angle.sin();
333 Coord::new(
334 vector.x * cos - vector.y * sin,
335 vector.x * sin + vector.y * cos,
336 )
337}
338
339fn simplify_tolerance(distance: f64, options: BufferOptions, right_side: bool) -> f64 {
340 let tolerance = distance * options.simplify_factor.max(0.0);
341 if right_side {
342 -tolerance
343 } else {
344 tolerance
345 }
346}
347
348#[derive(Debug, Clone, Copy)]
349struct Segment {
350 p0: Coord,
351 p1: Coord,
352}
353
354impl Segment {
355 const fn new(p0: Coord, p1: Coord) -> Self {
356 Self { p0, p1 }
357 }
358}
359
360struct OffsetSegmentGenerator {
361 distance: f64,
362 options: BufferOptions,
363 precision: PrecisionModel,
364 fillet_angle_quantum: f64,
365 closing_seg_length_factor: f64,
366 minimum_vertex_distance: f64,
367 points: Vec<Coord>,
368 s0: Coord,
369 s1: Coord,
370 s2: Coord,
371 offset0: Segment,
372 offset1: Segment,
373 side: Position,
374}
375
376impl OffsetSegmentGenerator {
377 fn new(distance: f64, options: BufferOptions, precision: PrecisionModel) -> Self {
378 let quadrant_segments = options.quadrant_segments.max(1) as f64;
379 let closing_seg_length_factor =
380 if options.quadrant_segments >= 8 && options.join_style == JoinStyle::Round {
381 MAX_CLOSING_SEG_LEN_FACTOR
382 } else {
383 1.0
384 };
385
386 Self {
387 distance,
388 options,
389 precision,
390 fillet_angle_quantum: PI / 2.0 / quadrant_segments,
391 closing_seg_length_factor,
392 minimum_vertex_distance: distance * CURVE_VERTEX_SNAP_DISTANCE_FACTOR,
393 points: Vec::new(),
394 s0: Coord::default(),
395 s1: Coord::default(),
396 s2: Coord::default(),
397 offset0: Segment::new(Coord::default(), Coord::default()),
398 offset1: Segment::new(Coord::default(), Coord::default()),
399 side: Position::Left,
400 }
401 }
402
403 fn coordinates(self) -> Vec<Coord> {
404 self.points
405 }
406
407 fn init_side_segments(&mut self, s1: Coord, s2: Coord, side: Position) -> Result<()> {
408 self.s1 = s1;
409 self.s2 = s2;
410 self.side = side;
411 self.offset1 = compute_offset_segment(Segment::new(s1, s2), side, self.distance)?;
412 Ok(())
413 }
414
415 fn add_next_segment(&mut self, point: Coord, add_start_point: bool) -> Result<()> {
416 self.s0 = self.s1;
417 self.s1 = self.s2;
418 self.s2 = point;
419 self.offset0 =
420 compute_offset_segment(Segment::new(self.s0, self.s1), self.side, self.distance)?;
421 self.offset1 =
422 compute_offset_segment(Segment::new(self.s1, self.s2), self.side, self.distance)?;
423
424 if same_coord_exact(self.s1, self.s2) {
425 return Ok(());
426 }
427
428 let orient = orientation_index(self.s0, self.s1, self.s2);
429 let outside_turn = (orient == CLOCKWISE && self.side == Position::Left)
430 || (orient == COUNTERCLOCKWISE && self.side == Position::Right);
431
432 if orient == COLLINEAR {
433 self.add_collinear(add_start_point);
434 } else if outside_turn {
435 self.add_outside_turn(orient, add_start_point);
436 } else {
437 self.add_inside_turn();
438 }
439 Ok(())
440 }
441
442 fn add_last_segment(&mut self) {
443 self.add_point(self.offset1.p1);
444 }
445
446 fn add_line_end_cap(&mut self, p0: Coord, p1: Coord) -> Result<()> {
447 let segment = Segment::new(p0, p1);
448 let offset_left = compute_offset_segment(segment, Position::Left, self.distance)?;
449 let offset_right = compute_offset_segment(segment, Position::Right, self.distance)?;
450 let angle = (p1.y - p0.y).atan2(p1.x - p0.x);
451
452 match self.options.cap_style {
453 CapStyle::Round => {
454 self.add_point(offset_left.p1);
455 self.add_directed_fillet(
456 p1,
457 angle + PI / 2.0,
458 angle - PI / 2.0,
459 CLOCKWISE,
460 self.distance,
461 );
462 self.add_point(offset_right.p1);
463 }
464 CapStyle::Flat => {
465 self.add_point(offset_left.p1);
466 self.add_point(offset_right.p1);
467 }
468 CapStyle::Square => {
469 let square_offset = Coord::new(
470 self.distance.abs() * angle.cos(),
471 self.distance.abs() * angle.sin(),
472 );
473 self.add_point(Coord::new(
474 offset_left.p1.x + square_offset.x,
475 offset_left.p1.y + square_offset.y,
476 ));
477 self.add_point(Coord::new(
478 offset_right.p1.x + square_offset.x,
479 offset_right.p1.y + square_offset.y,
480 ));
481 }
482 }
483
484 Ok(())
485 }
486
487 fn add_outside_turn(&mut self, orient: i8, add_start_point: bool) {
488 if self.offset0.p1.distance(self.offset1.p0)
489 < self.distance * OFFSET_SEGMENT_SEPARATION_FACTOR
490 {
491 self.add_point(self.offset0.p1);
492 return;
493 }
494
495 match self.options.join_style {
496 JoinStyle::Mitre => self.add_mitre_join(),
497 JoinStyle::Bevel => {
498 self.add_point(self.offset0.p1);
499 self.add_point(self.offset1.p0);
500 }
501 JoinStyle::Round => {
502 if add_start_point {
503 self.add_point(self.offset0.p1);
504 }
505 self.add_corner_fillet(
506 self.s1,
507 self.offset0.p1,
508 self.offset1.p0,
509 orient,
510 self.distance,
511 );
512 self.add_point(self.offset1.p0);
513 }
514 }
515 }
516
517 fn add_inside_turn(&mut self) {
518 if let Some(point) = segment_intersection_point(self.offset0, self.offset1, self.precision)
519 {
520 self.add_point(point);
521 return;
522 }
523
524 self.add_point(self.offset0.p1);
525 if self.offset0.p1.distance(self.offset1.p0)
526 >= self.distance * INSIDE_TURN_VERTEX_SNAP_DISTANCE_FACTOR
527 {
528 if self.closing_seg_length_factor > 0.0 {
529 self.add_point(Coord::new(
530 (self.closing_seg_length_factor * self.offset0.p1.x + self.s1.x)
531 / (self.closing_seg_length_factor + 1.0),
532 (self.closing_seg_length_factor * self.offset0.p1.y + self.s1.y)
533 / (self.closing_seg_length_factor + 1.0),
534 ));
535 self.add_point(Coord::new(
536 (self.closing_seg_length_factor * self.offset1.p0.x + self.s1.x)
537 / (self.closing_seg_length_factor + 1.0),
538 (self.closing_seg_length_factor * self.offset1.p0.y + self.s1.y)
539 / (self.closing_seg_length_factor + 1.0),
540 ));
541 } else {
542 self.add_point(self.s1);
543 }
544 self.add_point(self.offset1.p0);
545 }
546 }
547
548 fn add_collinear(&mut self, add_start_point: bool) {
549 let dx0 = self.s1.x - self.s0.x;
550 let dy0 = self.s1.y - self.s0.y;
551 let dx1 = self.s2.x - self.s1.x;
552 let dy1 = self.s2.y - self.s1.y;
553 if dx0 * dx1 + dy0 * dy1 >= 0.0 {
554 return;
555 }
556
557 if self.options.join_style == JoinStyle::Bevel
558 || self.options.join_style == JoinStyle::Mitre
559 {
560 if add_start_point {
561 self.add_point(self.offset0.p1);
562 }
563 self.add_point(self.offset1.p0);
564 return;
565 }
566
567 self.add_corner_fillet(
568 self.s1,
569 self.offset0.p1,
570 self.offset1.p0,
571 CLOCKWISE,
572 self.distance,
573 );
574 }
575
576 fn add_mitre_join(&mut self) {
577 if let Some(point) = line_intersection_centered(self.offset0, self.offset1, self.precision)
578 {
579 let mitre_ratio = if self.distance <= 0.0 {
580 1.0
581 } else {
582 point.distance(self.s1) / self.distance.abs()
583 };
584 if mitre_ratio <= self.options.mitre_limit {
585 self.add_point(point);
586 return;
587 }
588 }
589
590 self.add_point(self.offset0.p1);
591 self.add_point(self.offset1.p0);
592 }
593
594 fn add_corner_fillet(
595 &mut self,
596 center: Coord,
597 start: Coord,
598 end: Coord,
599 direction: i8,
600 radius: f64,
601 ) {
602 let dx0 = start.x - center.x;
603 let dy0 = start.y - center.y;
604 let mut start_angle = dy0.atan2(dx0);
605 let dx1 = end.x - center.x;
606 let dy1 = end.y - center.y;
607 let end_angle = dy1.atan2(dx1);
608
609 if direction == CLOCKWISE {
610 if start_angle <= end_angle {
611 start_angle += 2.0 * PI;
612 }
613 } else if start_angle >= end_angle {
614 start_angle -= 2.0 * PI;
615 }
616
617 self.add_point(start);
618 self.add_directed_fillet(center, start_angle, end_angle, direction, radius);
619 self.add_point(end);
620 }
621
622 fn add_directed_fillet(
623 &mut self,
624 center: Coord,
625 start_angle: f64,
626 end_angle: f64,
627 direction: i8,
628 radius: f64,
629 ) {
630 let direction_factor = if direction == CLOCKWISE { -1.0 } else { 1.0 };
631 let total_angle = (start_angle - end_angle).abs();
632 let n_segments = (total_angle / self.fillet_angle_quantum + 0.5).trunc() as usize;
633 if n_segments < 1 {
634 return;
635 }
636
637 let angle_inc = total_angle / n_segments as f64;
638 for i in 0..n_segments {
639 let angle = start_angle + direction_factor * i as f64 * angle_inc;
640 self.add_point(Coord::new(
641 center.x + radius * angle.cos(),
642 center.y + radius * angle.sin(),
643 ));
644 }
645 }
646
647 fn add_point(&mut self, point: Coord) {
648 let point = self.precision.snap_coord(point);
649 if self
650 .points
651 .last()
652 .is_some_and(|last| point.distance(*last) < self.minimum_vertex_distance)
653 {
654 return;
655 }
656 self.points.push(point);
657 }
658
659 fn close_ring(&mut self) {
660 let Some(first) = self.points.first().copied() else {
661 return;
662 };
663 if self
664 .points
665 .last()
666 .is_none_or(|last| !same_coord_exact(first, *last))
667 {
668 self.points.push(first);
669 }
670 }
671}
672
673fn compute_offset_segment(segment: Segment, side: Position, distance: f64) -> Result<Segment> {
674 let side_sign = if side == Position::Left { 1.0 } else { -1.0 };
675 let dx = segment.p1.x - segment.p0.x;
676 let dy = segment.p1.y - segment.p0.y;
677 let length = (dx * dx + dy * dy).sqrt();
678 if length == 0.0 {
679 return Err(GeometryError::InvalidGeometry(
680 "buffer encountered a zero-length segment".to_owned(),
681 ));
682 }
683 let ux = side_sign * distance * dx / length;
684 let uy = side_sign * distance * dy / length;
685 Ok(Segment::new(
686 Coord::new(segment.p0.x - uy, segment.p0.y + ux),
687 Coord::new(segment.p1.x - uy, segment.p1.y + ux),
688 ))
689}
690
691fn segment_intersection_point(
692 left: Segment,
693 right: Segment,
694 precision: PrecisionModel,
695) -> Option<Coord> {
696 match segment_intersection(left.p0, left.p1, right.p0, right.p1, precision) {
697 Some(SegmentIntersection::Point(point)) => Some(point),
698 _ => None,
699 }
700}
701
702fn line_intersection_centered(
703 left: Segment,
704 right: Segment,
705 precision: PrecisionModel,
706) -> Option<Coord> {
707 let min_x = left.p0.x.min(left.p1.x).max(right.p0.x.min(right.p1.x));
708 let max_x = left.p0.x.max(left.p1.x).min(right.p0.x.max(right.p1.x));
709 let min_y = left.p0.y.min(left.p1.y).max(right.p0.y.min(right.p1.y));
710 let max_y = left.p0.y.max(left.p1.y).min(right.p0.y.max(right.p1.y));
711 let mid_x = (min_x + max_x) / 2.0;
712 let mid_y = (min_y + max_y) / 2.0;
713
714 let p1 = Coord::new(left.p0.x - mid_x, left.p0.y - mid_y);
715 let p2 = Coord::new(left.p1.x - mid_x, left.p1.y - mid_y);
716 let q1 = Coord::new(right.p0.x - mid_x, right.p0.y - mid_y);
717 let q2 = Coord::new(right.p1.x - mid_x, right.p1.y - mid_y);
718
719 let px = p1.y - p2.y;
720 let py = p2.x - p1.x;
721 let pw = p1.x * p2.y - p2.x * p1.y;
722 let qx = q1.y - q2.y;
723 let qy = q2.x - q1.x;
724 let qw = q1.x * q2.y - q2.x * q1.y;
725 let x = py * qw - qy * pw;
726 let y = qx * pw - px * qw;
727 let w = px * qy - qx * py;
728
729 if w == 0.0 {
730 return None;
731 }
732
733 let x_int = x / w;
734 let y_int = y / w;
735 if !x_int.is_finite() || !y_int.is_finite() {
736 return None;
737 }
738
739 Some(precision.snap_coord(Coord::new(x_int + mid_x, y_int + mid_y)))
740}
741
742fn simplify_buffer_input_line(input: &[Coord], distance_tol: f64) -> Vec<Coord> {
743 let distance_tol_abs = distance_tol.abs();
744 let angle_orientation = if distance_tol < 0.0 {
745 CLOCKWISE
746 } else {
747 COUNTERCLOCKWISE
748 };
749 let mut is_deleted = vec![false; input.len()];
750
751 loop {
752 let changed =
753 delete_shallow_concavities(input, &mut is_deleted, distance_tol_abs, angle_orientation);
754 if !changed {
755 break;
756 }
757 }
758
759 input
760 .iter()
761 .zip(is_deleted)
762 .filter_map(|(coord, deleted)| (!deleted).then_some(*coord))
763 .collect()
764}
765
766fn delete_shallow_concavities(
767 input: &[Coord],
768 is_deleted: &mut [bool],
769 distance_tol: f64,
770 angle_orientation: i8,
771) -> bool {
772 let mut index = 1;
773 let mut mid_index = find_next_non_deleted_index(index, is_deleted);
774 let mut last_index = find_next_non_deleted_index(mid_index, is_deleted);
775 let mut changed = false;
776
777 while last_index < input.len() {
778 let mut middle_deleted = false;
779 if is_deletable(
780 input,
781 index,
782 mid_index,
783 last_index,
784 distance_tol,
785 angle_orientation,
786 ) {
787 is_deleted[mid_index] = true;
788 middle_deleted = true;
789 changed = true;
790 }
791
792 if middle_deleted {
793 index = last_index;
794 } else {
795 index = mid_index;
796 }
797 mid_index = find_next_non_deleted_index(index, is_deleted);
798 last_index = find_next_non_deleted_index(mid_index, is_deleted);
799 }
800
801 changed
802}
803
804fn find_next_non_deleted_index(index: usize, is_deleted: &[bool]) -> usize {
805 let mut next = index + 1;
806 while next < is_deleted.len() && is_deleted[next] {
807 next += 1;
808 }
809 next
810}
811
812fn is_deletable(
813 input: &[Coord],
814 i0: usize,
815 i1: usize,
816 i2: usize,
817 distance_tol: f64,
818 angle_orientation: i8,
819) -> bool {
820 let p0 = input[i0];
821 let p1 = input[i1];
822 let p2 = input[i2];
823 orientation_index(p0, p1, p2) == angle_orientation
824 && point_to_segment_distance(p1, p0, p2) < distance_tol
825 && is_shallow_sampled(input, p0, p1, i0, i2, distance_tol)
826}
827
828fn is_shallow_sampled(
829 input: &[Coord],
830 p0: Coord,
831 p1: Coord,
832 i0: usize,
833 i2: usize,
834 distance_tol: f64,
835) -> bool {
836 let mut increment = (i2 - i0) / NUM_SIMPLIFY_PTS_TO_CHECK;
837 if increment == 0 {
838 increment = 1;
839 }
840
841 let mut i = i0;
842 while i < i2 {
843 if point_to_segment_distance(p1, p0, input[i]) >= distance_tol {
844 return false;
845 }
846 i += increment;
847 }
848 true
849}
850
851fn point_to_segment_distance(point: Coord, a: Coord, b: Coord) -> f64 {
852 if same_coord_exact(a, b) {
853 return point.distance(a);
854 }
855
856 let len2 = b.distance_squared(a);
857 let r = ((point.x - a.x) * (b.x - a.x) + (point.y - a.y) * (b.y - a.y)) / len2;
858 if r <= 0.0 {
859 return point.distance(a);
860 }
861 if r >= 1.0 {
862 return point.distance(b);
863 }
864
865 let s = ((a.y - point.y) * (b.x - a.x) - (a.x - point.x) * (b.y - a.y)) / len2;
866 s.abs() * len2.sqrt()
867}
868
869fn is_eroded_completely(ring: &[Coord], buffer_distance: f64) -> bool {
870 if ring.len() < 4 {
871 return buffer_distance < 0.0;
872 }
873 let Some(bbox) = crate::types::BBox::from_coords(ring) else {
874 return true;
875 };
876 let min_dimension = (bbox.max.y - bbox.min.y).min(bbox.max.x - bbox.min.x);
877 buffer_distance < 0.0 && 2.0 * buffer_distance.abs() > min_dimension
878}
879
880#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
881struct CoordKey(u64, u64);
882
883#[derive(Debug, Clone)]
884struct HalfEdge {
885 from: usize,
886 to: usize,
887 reverse: usize,
888 angle: f64,
889 visited: bool,
890}
891
892fn polygonize_buffer_curves(curves: &[Vec<Coord>], precision: PrecisionModel) -> Vec<Polygon> {
893 let lines = curves
894 .iter()
895 .filter(|curve| curve.len() >= 4)
896 .map(|curve| LineString::new(curve.clone()))
897 .collect::<Vec<_>>();
898 let noded = node_lines(&lines, precision);
899 if noded.lines.is_empty() {
900 return Vec::new();
901 }
902
903 let mut vertices = Vec::<Coord>::new();
904 let mut vertex_ids = HashMap::<CoordKey, usize>::new();
905 let mut undirected_edges = HashSet::<(usize, usize)>::new();
906 let mut edges = Vec::<HalfEdge>::new();
907 let mut outgoing = Vec::<Vec<usize>>::new();
908
909 for line in noded.lines {
910 if line.coords.len() != 2 || precision.same_coord(line.coords[0], line.coords[1]) {
911 continue;
912 }
913 let from = vertex_id(
914 line.coords[0],
915 &mut vertices,
916 &mut vertex_ids,
917 &mut outgoing,
918 );
919 let to = vertex_id(
920 line.coords[1],
921 &mut vertices,
922 &mut vertex_ids,
923 &mut outgoing,
924 );
925 let edge_key = if from < to { (from, to) } else { (to, from) };
926 if !undirected_edges.insert(edge_key) {
927 continue;
928 }
929
930 let forward = edges.len();
931 let reverse = forward + 1;
932 edges.push(HalfEdge {
933 from,
934 to,
935 reverse,
936 angle: angle_between(vertices[from], vertices[to]),
937 visited: false,
938 });
939 edges.push(HalfEdge {
940 from: to,
941 to: from,
942 reverse: forward,
943 angle: angle_between(vertices[to], vertices[from]),
944 visited: false,
945 });
946 outgoing[from].push(forward);
947 outgoing[to].push(reverse);
948 }
949
950 for edge_ids in &mut outgoing {
951 edge_ids.sort_by(|left, right| {
952 edges[*left]
953 .angle
954 .partial_cmp(&edges[*right].angle)
955 .unwrap_or(std::cmp::Ordering::Equal)
956 });
957 }
958
959 let mut shells = Vec::<LinearRing>::new();
960 let mut holes = Vec::<LinearRing>::new();
961 for edge_id in 0..edges.len() {
962 if edges[edge_id].visited {
963 continue;
964 }
965 let Some(ring) = trace_face(edge_id, &mut edges, &outgoing, &vertices, precision) else {
966 continue;
967 };
968 let area = signed_ring_area(&ring);
969 if area > precision.epsilon() {
970 shells.push(ring);
971 } else if area < -precision.epsilon() {
972 holes.push(ring);
973 }
974 }
975 assemble_polygons(shells, holes, precision)
976}
977
978fn vertex_id(
979 coord: Coord,
980 vertices: &mut Vec<Coord>,
981 vertex_ids: &mut HashMap<CoordKey, usize>,
982 outgoing: &mut Vec<Vec<usize>>,
983) -> usize {
984 let key = coord_key(coord);
985 if let Some(id) = vertex_ids.get(&key) {
986 return *id;
987 }
988
989 let id = vertices.len();
990 vertices.push(coord);
991 vertex_ids.insert(key, id);
992 outgoing.push(Vec::new());
993 id
994}
995
996fn trace_face(
997 start_edge: usize,
998 edges: &mut [HalfEdge],
999 outgoing: &[Vec<usize>],
1000 vertices: &[Coord],
1001 precision: PrecisionModel,
1002) -> Option<LinearRing> {
1003 let mut edge_id = start_edge;
1004 let mut coords = Vec::new();
1005
1006 loop {
1007 if edges[edge_id].visited {
1008 if edge_id == start_edge {
1009 break;
1010 }
1011 return None;
1012 }
1013
1014 edges[edge_id].visited = true;
1015 coords.push(vertices[edges[edge_id].from]);
1016
1017 let reverse = edges[edge_id].reverse;
1018 let destination = edges[edge_id].to;
1019 let edge_ids = &outgoing[destination];
1020 let reverse_index = edge_ids.iter().position(|id| *id == reverse)?;
1021 let next_index = if reverse_index == 0 {
1022 edge_ids.len() - 1
1023 } else {
1024 reverse_index - 1
1025 };
1026 edge_id = edge_ids[next_index];
1027
1028 if edge_id == start_edge {
1029 break;
1030 }
1031 }
1032
1033 if coords.len() < 3 {
1034 return None;
1035 }
1036 coords.push(coords[0]);
1037 let ring = LinearRing::new(coords);
1038 (signed_ring_area(&ring).abs() > precision.epsilon()).then_some(ring)
1039}
1040
1041fn assemble_polygons(
1042 shells: Vec<LinearRing>,
1043 holes: Vec<LinearRing>,
1044 precision: PrecisionModel,
1045) -> Vec<Polygon> {
1046 let mut polygons = shells
1047 .into_iter()
1048 .map(|shell| Polygon::new(shell, Vec::new()))
1049 .collect::<Vec<_>>();
1050
1051 for hole in holes {
1052 let Some(sample) = hole.coords.first().copied() else {
1053 continue;
1054 };
1055 let mut target_index = None;
1056 let mut target_area = f64::INFINITY;
1057 for (index, polygon) in polygons.iter().enumerate() {
1058 if matches!(
1059 point_in_ring(sample, &polygon.exterior, precision),
1060 PointLocation::Interior
1061 ) {
1062 let area = signed_ring_area(&polygon.exterior).abs();
1063 if area < target_area {
1064 target_area = area;
1065 target_index = Some(index);
1066 }
1067 }
1068 }
1069 if let Some(index) = target_index {
1070 polygons[index].holes.push(hole);
1071 }
1072 }
1073
1074 polygons
1075 .into_iter()
1076 .filter(|polygon| !polygon.is_empty())
1077 .collect()
1078}
1079
1080fn angle_between(from: Coord, to: Coord) -> f64 {
1081 (to.y - from.y).atan2(to.x - from.x)
1082}
1083
1084fn coord_key(coord: Coord) -> CoordKey {
1085 CoordKey(
1086 normalize_zero(coord.x).to_bits(),
1087 normalize_zero(coord.y).to_bits(),
1088 )
1089}
1090
1091fn normalize_zero(value: f64) -> f64 {
1092 if value == 0.0 {
1093 0.0
1094 } else {
1095 value
1096 }
1097}
1098
1099fn remove_repeated_points(coords: &[Coord], precision: PrecisionModel) -> Vec<Coord> {
1100 let mut out = Vec::with_capacity(coords.len());
1101 for coord in coords {
1102 let coord = precision.snap_coord(*coord);
1103 if out
1104 .last()
1105 .is_none_or(|last| !same_coord_exact(*last, coord))
1106 {
1107 out.push(coord);
1108 }
1109 }
1110 out
1111}
1112
1113fn is_ccw_jts(ring: &[Coord]) -> bool {
1114 if ring.len() < 4 {
1115 return is_ring_ccw(&LinearRing::new(ring.to_vec()));
1116 }
1117
1118 let n = ring.len() - 1;
1119 let mut high_point = ring[0];
1120 let mut high_index = 0;
1121 for (index, point) in ring.iter().enumerate().take(n + 1).skip(1) {
1122 if point.y > high_point.y {
1123 high_point = *point;
1124 high_index = index;
1125 }
1126 }
1127
1128 let mut prev_index = high_index;
1129 loop {
1130 prev_index = if prev_index == 0 { n } else { prev_index - 1 };
1131 if !same_coord_exact(ring[prev_index], high_point) || prev_index == high_index {
1132 break;
1133 }
1134 }
1135
1136 let mut next_index = high_index;
1137 loop {
1138 next_index = (next_index + 1) % n;
1139 if !same_coord_exact(ring[next_index], high_point) || next_index == high_index {
1140 break;
1141 }
1142 }
1143
1144 let prev = ring[prev_index];
1145 let next = ring[next_index];
1146 if same_coord_exact(prev, high_point)
1147 || same_coord_exact(next, high_point)
1148 || same_coord_exact(prev, next)
1149 {
1150 return false;
1151 }
1152
1153 let disc = orientation_index(prev, high_point, next);
1154 if disc == COLLINEAR {
1155 prev.x > next.x
1156 } else {
1157 disc > 0
1158 }
1159}
1160
1161fn orientation_index(a: Coord, b: Coord, c: Coord) -> i8 {
1162 let value = orientation(a, b, c);
1163 if value > 0.0 {
1164 COUNTERCLOCKWISE
1165 } else if value < 0.0 {
1166 CLOCKWISE
1167 } else {
1168 COLLINEAR
1169 }
1170}
1171
1172fn same_coord_exact(a: Coord, b: Coord) -> bool {
1173 a.x == b.x && a.y == b.y
1174}