1 use crate::prelude::*;
2 use crate::{
3 Coordinate, Line, LineString, MultiLineString, MultiPolygon, Point, Polygon, Triangle,
4 };
5 use num_traits::Float;
6 use std::cmp::Ordering;
7 use std::collections::BinaryHeap;
8
9 use rstar::{RTree, RTreeNum};
10
11 /// Store triangle information
12 // current is the candidate point for removal
13 #[derive(Debug)]
14 struct VScore<T, I>
15 where
16 T: Float,
17 {
18 left: usize,
19 current: usize,
20 right: usize,
21 area: T,
22 // `visvalingam_preserve` uses `intersector`, `visvalingam` does not
23 intersector: I,
24 }
25
26 // These impls give us a min-heap
27 impl<T, I> Ord for VScore<T, I>
28 where
29 T: Float,
30 {
cmp(&self, other: &VScore<T, I>) -> Ordering31 fn cmp(&self, other: &VScore<T, I>) -> Ordering {
32 other.area.partial_cmp(&self.area).unwrap()
33 }
34 }
35
36 impl<T, I> PartialOrd for VScore<T, I>
37 where
38 T: Float,
39 {
partial_cmp(&self, other: &VScore<T, I>) -> Option<Ordering>40 fn partial_cmp(&self, other: &VScore<T, I>) -> Option<Ordering> {
41 Some(self.cmp(other))
42 }
43 }
44
45 impl<T, I> Eq for VScore<T, I> where T: Float {}
46
47 impl<T, I> PartialEq for VScore<T, I>
48 where
49 T: Float,
50 {
eq(&self, other: &VScore<T, I>) -> bool where T: Float,51 fn eq(&self, other: &VScore<T, I>) -> bool
52 where
53 T: Float,
54 {
55 self.area == other.area
56 }
57 }
58
59 /// Geometries that can be simplified using the topology-preserving variant
60 #[derive(Debug, Clone, Copy)]
61 enum GeomType {
62 Line,
63 Ring,
64 }
65
66 /// Settings for Ring and Line geometries
67 // initial min: if we ever have fewer than these, stop immediately
68 // min_points: if we detect a self-intersection before point removal, and we only
69 // have min_points left, stop: since a self-intersection causes removal of the spatially previous
70 // point, THAT could lead to a further self-intersection without the possibility of removing
71 // more points, potentially leaving the geometry in an invalid state.
72 #[derive(Debug, Clone, Copy)]
73 struct GeomSettings {
74 initial_min: usize,
75 min_points: usize,
76 geomtype: GeomType,
77 }
78
79 /// Simplify a line using the [Visvalingam-Whyatt](http://www.tandfonline.com/doi/abs/10.1179/000870493786962263) algorithm
80 //
81 // This method returns the **indices** of the simplified line
82 // epsilon is the minimum triangle area
83 // The paper states that:
84 // If [the new triangle's] calculated area is less than that of the last point to be
85 // eliminated, use the latter's area instead.
86 // (This ensures that the current point cannot be eliminated
87 // without eliminating previously eliminated points)
88 // (Visvalingam and Whyatt 2013, p47)
89 // However, this does *not* apply if you're using a user-defined epsilon;
90 // It's OK to remove triangles with areas below the epsilon,
91 // then recalculate the new triangle area and push it onto the heap
92 // based on Huon Wilson's original implementation:
93 // https://github.com/huonw/isrustfastyet/blob/25e7a68ff26673a8556b170d3c9af52e1c818288/mem/line_simplify.rs
visvalingam_indices<T>(orig: &LineString<T>, epsilon: &T) -> Vec<usize> where T: Float,94 fn visvalingam_indices<T>(orig: &LineString<T>, epsilon: &T) -> Vec<usize>
95 where
96 T: Float,
97 {
98 // No need to continue without at least three points
99 if orig.0.len() < 3 {
100 return orig.0.iter().enumerate().map(|(idx, _)| idx).collect();
101 }
102
103 let max = orig.0.len();
104
105 // Adjacent retained points. Simulating the points in a
106 // linked list with indices into `orig`. Big number (larger than or equal to
107 // `max`) means no next element, and (0, 0) means deleted element.
108 let mut adjacent: Vec<_> = (0..orig.0.len())
109 .map(|i| {
110 if i == 0 {
111 (-1_i32, 1_i32)
112 } else {
113 ((i - 1) as i32, (i + 1) as i32)
114 }
115 })
116 .collect();
117
118 // Store all the triangles in a minimum priority queue, based on their area.
119 //
120 // Invalid triangles are *not* removed if / when points are removed; they're
121 // handled by skipping them as necessary in the main loop by checking the
122 // corresponding entry in adjacent for (0, 0) values
123
124 // Compute the initial triangles
125 let mut pq = orig
126 .triangles()
127 .enumerate()
128 .map(|(i, triangle)| VScore {
129 area: triangle.area().abs(),
130 current: i + 1,
131 left: i,
132 right: i + 2,
133 intersector: (),
134 })
135 .collect::<BinaryHeap<VScore<T, ()>>>();
136 // While there are still points for which the associated triangle
137 // has an area below the epsilon
138 while let Some(smallest) = pq.pop() {
139 // This triangle's area is above epsilon, so skip it
140 if smallest.area > *epsilon {
141 continue;
142 }
143 // This triangle's area is below epsilon: eliminate the associated point
144 let (left, right) = adjacent[smallest.current];
145 // A point in this triangle has been removed since this VScore
146 // was created, so skip it
147 if left as i32 != smallest.left as i32 || right as i32 != smallest.right as i32 {
148 continue;
149 }
150 // We've got a valid triangle, and its area is smaller than epsilon, so
151 // remove it from the simulated "linked list"
152 let (ll, _) = adjacent[left as usize];
153 let (_, rr) = adjacent[right as usize];
154 adjacent[left as usize] = (ll, right);
155 adjacent[right as usize] = (left, rr);
156 adjacent[smallest.current as usize] = (0, 0);
157
158 // Now recompute the adjacent triangle(s), using left and right adjacent points
159 let choices = [(ll, left, right), (left, right, rr)];
160 for &(ai, current_point, bi) in &choices {
161 if ai as usize >= max || bi as usize >= max {
162 // Out of bounds, i.e. we're on one edge
163 continue;
164 }
165 let area = Triangle(
166 orig.0[ai as usize],
167 orig.0[current_point as usize],
168 orig.0[bi as usize],
169 )
170 .area()
171 .abs();
172 pq.push(VScore {
173 area,
174 current: current_point as usize,
175 left: ai as usize,
176 right: bi as usize,
177 intersector: (),
178 });
179 }
180 }
181 // Filter out the points that have been deleted, returning remaining point indices
182 orig.0
183 .iter()
184 .enumerate()
185 .zip(adjacent.iter())
186 .filter_map(|(tup, adj)| if *adj != (0, 0) { Some(tup.0) } else { None })
187 .collect::<Vec<usize>>()
188 }
189
190 // Wrapper for visvalingam_indices, mapping indices back to points
visvalingam<T>(orig: &LineString<T>, epsilon: &T) -> Vec<Coordinate<T>> where T: Float,191 fn visvalingam<T>(orig: &LineString<T>, epsilon: &T) -> Vec<Coordinate<T>>
192 where
193 T: Float,
194 {
195 let subset = visvalingam_indices(orig, epsilon);
196 // filter orig using the indices
197 // using get would be more robust here, but the input subset is guaranteed to be valid in this case
198 orig.0
199 .iter()
200 .zip(subset.iter())
201 .map(|(_, s)| orig[*s])
202 .collect()
203 }
204
205 /// Wrap the actual VW function so the R* Tree can be shared.
206 // this ensures that shell and rings have access to all segments, so
207 // intersections between outer and inner rings are detected
vwp_wrapper<T>( geomtype: &GeomSettings, exterior: &LineString<T>, interiors: Option<&[LineString<T>]>, epsilon: &T, ) -> Vec<Vec<Coordinate<T>>> where T: Float + RTreeNum,208 fn vwp_wrapper<T>(
209 geomtype: &GeomSettings,
210 exterior: &LineString<T>,
211 interiors: Option<&[LineString<T>]>,
212 epsilon: &T,
213 ) -> Vec<Vec<Coordinate<T>>>
214 where
215 T: Float + RTreeNum,
216 {
217 let mut rings = vec![];
218 // Populate R* tree with exterior and interior samples, if any
219 let mut tree: RTree<Line<_>> = RTree::bulk_load(
220 exterior
221 .lines()
222 .chain(
223 interiors
224 .iter()
225 .flat_map(|ring| *ring)
226 .flat_map(|line_string| line_string.lines()),
227 )
228 .collect::<Vec<_>>(),
229 );
230
231 // Simplify shell
232 rings.push(visvalingam_preserve(
233 geomtype, &exterior, epsilon, &mut tree,
234 ));
235 // Simplify interior rings, if any
236 if let Some(interior_rings) = interiors {
237 for ring in interior_rings {
238 rings.push(visvalingam_preserve(geomtype, &ring, epsilon, &mut tree))
239 }
240 }
241 rings
242 }
243
244 /// Visvalingam-Whyatt with self-intersection detection to preserve topologies
245 /// this is a port of the technique at https://www.jasondavies.com/simplify/
visvalingam_preserve<T>( geomtype: &GeomSettings, orig: &LineString<T>, epsilon: &T, tree: &mut RTree<Line<T>>, ) -> Vec<Coordinate<T>> where T: Float + RTreeNum,246 fn visvalingam_preserve<T>(
247 geomtype: &GeomSettings,
248 orig: &LineString<T>,
249 epsilon: &T,
250 tree: &mut RTree<Line<T>>,
251 ) -> Vec<Coordinate<T>>
252 where
253 T: Float + RTreeNum,
254 {
255 if orig.0.len() < 3 {
256 return orig.0.to_vec();
257 }
258 let max = orig.0.len();
259 let mut counter = orig.0.len();
260
261 // Adjacent retained points. Simulating the points in a
262 // linked list with indices into `orig`. Big number (larger than or equal to
263 // `max`) means no next element, and (0, 0) means deleted element.
264 let mut adjacent: Vec<_> = (0..orig.0.len())
265 .map(|i| {
266 if i == 0 {
267 (-1_i32, 1_i32)
268 } else {
269 ((i - 1) as i32, (i + 1) as i32)
270 }
271 })
272 .collect();
273 // Store all the triangles in a minimum priority queue, based on their area.
274 //
275 // Invalid triangles are *not* removed if / when points are removed; they're
276 // handled by skipping them as necessary in the main loop by checking the
277 // corresponding entry in adjacent for (0, 0) values
278
279 // Compute the initial triangles
280 let mut pq = orig
281 .triangles()
282 .enumerate()
283 .map(|(i, triangle)| VScore {
284 area: triangle.area().abs(),
285 current: i + 1,
286 left: i,
287 right: i + 2,
288 intersector: false,
289 })
290 .collect::<BinaryHeap<VScore<T, bool>>>();
291
292 // While there are still points for which the associated triangle
293 // has an area below the epsilon
294 while let Some(mut smallest) = pq.pop() {
295 if smallest.area > *epsilon {
296 continue;
297 }
298 if counter <= geomtype.initial_min {
299 // we can't remove any more points no matter what
300 break;
301 }
302 let (left, right) = adjacent[smallest.current];
303 // A point in this triangle has been removed since this VScore
304 // was created, so skip it
305 if left as i32 != smallest.left as i32 || right as i32 != smallest.right as i32 {
306 continue;
307 }
308 // if removal of this point causes a self-intersection, we also remove the previous point
309 // that removal alters the geometry, removing the self-intersection
310 // HOWEVER if we're within 2 points of the absolute minimum, we can't remove this point or the next
311 // because we could then no longer form a valid geometry if removal of next also caused an intersection.
312 // The simplification process is thus over.
313 smallest.intersector = tree_intersect(tree, &smallest, &orig.0);
314 if smallest.intersector && counter <= geomtype.min_points {
315 break;
316 }
317 // We've got a valid triangle, and its area is smaller than epsilon, so
318 // remove it from the simulated "linked list"
319 adjacent[smallest.current as usize] = (0, 0);
320 counter -= 1;
321 // Remove stale segments from R* tree
322 let left_point = Point(orig.0[left as usize]);
323 let middle_point = Point(orig.0[smallest.current]);
324 let right_point = Point(orig.0[right as usize]);
325
326 let line_1 = Line::new(left_point, middle_point);
327 let line_2 = Line::new(middle_point, right_point);
328 assert!(tree.remove(&line_1).is_some());
329 assert!(tree.remove(&line_2).is_some());
330
331 // Restore continous line segment
332 tree.insert(Line::new(left_point, right_point));
333
334 // Now recompute the adjacent triangle(s), using left and right adjacent points
335 let (ll, _) = adjacent[left as usize];
336 let (_, rr) = adjacent[right as usize];
337 adjacent[left as usize] = (ll, right);
338 adjacent[right as usize] = (left, rr);
339 let choices = [(ll, left, right), (left, right, rr)];
340 for &(ai, current_point, bi) in &choices {
341 if ai as usize >= max || bi as usize >= max {
342 // Out of bounds, i.e. we're on one edge
343 continue;
344 }
345 let new = Triangle(
346 orig.0[ai as usize],
347 orig.0[current_point as usize],
348 orig.0[bi as usize],
349 );
350 // The current point causes a self-intersection, and this point precedes it
351 // we ensure it gets removed next by demoting its area to negative epsilon
352 let temp_area = if smallest.intersector && (current_point as usize) < smallest.current {
353 -*epsilon
354 } else {
355 new.area().abs()
356 };
357 let new_triangle = VScore {
358 area: temp_area,
359 current: current_point as usize,
360 left: ai as usize,
361 right: bi as usize,
362 intersector: false,
363 };
364
365 // push re-computed triangle onto heap
366 pq.push(new_triangle);
367 }
368 }
369 // Filter out the points that have been deleted, returning remaining points
370 orig.0
371 .iter()
372 .zip(adjacent.iter())
373 .filter_map(|(tup, adj)| if *adj != (0, 0) { Some(*tup) } else { None })
374 .collect()
375 }
376
377 /// is p1 -> p2 -> p3 wound counterclockwise?
ccw<T>(p1: Point<T>, p2: Point<T>, p3: Point<T>) -> bool where T: Float,378 fn ccw<T>(p1: Point<T>, p2: Point<T>, p3: Point<T>) -> bool
379 where
380 T: Float,
381 {
382 (p3.y() - p1.y()) * (p2.x() - p1.x()) > (p2.y() - p1.y()) * (p3.x() - p1.x())
383 }
384
385 /// checks whether line segments with p1-p4 as their start and endpoints touch or cross
cartesian_intersect<T>(p1: Point<T>, p2: Point<T>, p3: Point<T>, p4: Point<T>) -> bool where T: Float,386 fn cartesian_intersect<T>(p1: Point<T>, p2: Point<T>, p3: Point<T>, p4: Point<T>) -> bool
387 where
388 T: Float,
389 {
390 (ccw(p1, p3, p4) ^ ccw(p2, p3, p4)) & (ccw(p1, p2, p3) ^ ccw(p1, p2, p4))
391 }
392
393 /// check whether a triangle's edges intersect with any other edges of the LineString
tree_intersect<T>( tree: &RTree<Line<T>>, triangle: &VScore<T, bool>, orig: &[Coordinate<T>], ) -> bool where T: Float + RTreeNum,394 fn tree_intersect<T>(
395 tree: &RTree<Line<T>>,
396 triangle: &VScore<T, bool>,
397 orig: &[Coordinate<T>],
398 ) -> bool
399 where
400 T: Float + RTreeNum,
401 {
402 let point_a = orig[triangle.left];
403 let point_c = orig[triangle.right];
404 let bounding_rect = Triangle(
405 orig[triangle.left],
406 orig[triangle.current],
407 orig[triangle.right],
408 )
409 .bounding_rect();
410 let br = Point::new(bounding_rect.min().x, bounding_rect.min().y);
411 let tl = Point::new(bounding_rect.max().x, bounding_rect.max().y);
412 tree.locate_in_envelope_intersecting(&rstar::AABB::from_corners(br, tl))
413 .any(|c| {
414 // triangle start point, end point
415 let (ca, cb) = c.points();
416 ca.0 != point_a
417 && ca.0 != point_c
418 && cb.0 != point_a
419 && cb.0 != point_c
420 && cartesian_intersect(ca, cb, Point(point_a), Point(point_c))
421 })
422 }
423
424 /// Simplifies a geometry.
425 ///
426 /// Polygons are simplified by running the algorithm on all their constituent rings. This may
427 /// result in invalid Polygons, and has no guarantee of preserving topology. Multi* objects are
428 /// simplified by simplifying all their constituent geometries individually.
429 pub trait SimplifyVW<T, Epsilon = T> {
430 /// Returns the simplified representation of a geometry, using the [Visvalingam-Whyatt](http://www.tandfonline.com/doi/abs/10.1179/000870493786962263) algorithm
431 ///
432 /// See [here](https://bost.ocks.org/mike/simplify/) for a graphical explanation
433 ///
434 /// # Examples
435 ///
436 /// ```
437 /// use geo::algorithm::simplifyvw::SimplifyVW;
438 /// use geo::{LineString, Point};
439 ///
440 /// let mut vec = Vec::new();
441 /// vec.push(Point::new(5.0, 2.0));
442 /// vec.push(Point::new(3.0, 8.0));
443 /// vec.push(Point::new(6.0, 20.0));
444 /// vec.push(Point::new(7.0, 25.0));
445 /// vec.push(Point::new(10.0, 10.0));
446 /// let linestring = LineString::from(vec);
447 /// let mut compare = Vec::new();
448 /// compare.push(Point::new(5.0, 2.0));
449 /// compare.push(Point::new(7.0, 25.0));
450 /// compare.push(Point::new(10.0, 10.0));
451 /// let ls_compare = LineString::from(compare);
452 /// let simplified = linestring.simplifyvw(&30.0);
453 /// assert_eq!(simplified, ls_compare)
454 /// ```
simplifyvw(&self, epsilon: &T) -> Self where T: Float455 fn simplifyvw(&self, epsilon: &T) -> Self
456 where
457 T: Float;
458 }
459
460 /// Simplifies a geometry, returning the retained _indices_ of the output
461 ///
462 /// This operation uses the Visvalingam-Whyatt algorithm,
463 /// and does **not** guarantee that the returned geometry is valid.
464 pub trait SimplifyVwIdx<T, Epsilon = T> {
465 /// Returns the simplified representation of a geometry, using the [Visvalingam-Whyatt](http://www.tandfonline.com/doi/abs/10.1179/000870493786962263) algorithm
466 ///
467 /// See [here](https://bost.ocks.org/mike/simplify/) for a graphical explanation
468 ///
469 /// # Examples
470 ///
471 /// ```
472 /// use geo::algorithm::simplifyvw::SimplifyVwIdx;
473 /// use geo::{LineString, Point};
474 ///
475 /// let mut vec = Vec::new();
476 /// vec.push(Point::new(5.0, 2.0));
477 /// vec.push(Point::new(3.0, 8.0));
478 /// vec.push(Point::new(6.0, 20.0));
479 /// vec.push(Point::new(7.0, 25.0));
480 /// vec.push(Point::new(10.0, 10.0));
481 /// let linestring = LineString::from(vec);
482 /// let mut compare = Vec::new();
483 /// compare.push(0_usize);
484 /// compare.push(3_usize);
485 /// compare.push(4_usize);
486 /// let simplified = linestring.simplifyvw_idx(&30.0);
487 /// assert_eq!(simplified, compare)
488 /// ```
simplifyvw_idx(&self, epsilon: &T) -> Vec<usize> where T: Float489 fn simplifyvw_idx(&self, epsilon: &T) -> Vec<usize>
490 where
491 T: Float;
492 }
493
494 /// Simplifies a geometry, preserving its topology by removing self-intersections
495 pub trait SimplifyVWPreserve<T, Epsilon = T> {
496 /// Returns the simplified representation of a geometry, using a topology-preserving variant of the
497 /// [Visvalingam-Whyatt](http://www.tandfonline.com/doi/abs/10.1179/000870493786962263) algorithm.
498 ///
499 /// See [here](https://www.jasondavies.com/simplify/) for a graphical explanation.
500 ///
501 /// The topology-preserving algorithm uses an [R* tree](../../../rstar/struct.RTree.html) to
502 /// efficiently find candidate line segments which are tested for intersection with a given triangle.
503 /// If intersections are found, the previous point (i.e. the left component of the current triangle)
504 /// is also removed, altering the geometry and removing the intersection.
505 ///
506 /// In the example below, `(135.0, 68.0)` would be retained by the standard algorithm,
507 /// forming triangle `(0, 1, 3),` which intersects with the segments `(280.0, 19.0),
508 /// (117.0, 48.0)` and `(117.0, 48.0), (300,0, 40.0)`. By removing it,
509 /// a new triangle with indices `(0, 3, 4)` is formed, which does not cause a self-intersection.
510 ///
511 /// **Note**: it is possible for the simplification algorithm to displace a Polygon's interior ring outside its shell.
512 ///
513 /// **Note**: if removal of a point causes a self-intersection, but the geometry only has `n + 2`
514 /// points remaining (4 for a `LineString`, 6 for a `Polygon`), the point is retained and the
515 /// simplification process ends. This is because there is no guarantee that removal of two points will remove
516 /// the intersection, but removal of further points would leave too few points to form a valid geometry.
517 ///
518 /// # Examples
519 ///
520 /// ```
521 /// use geo::algorithm::simplifyvw::SimplifyVWPreserve;
522 /// use geo::{LineString, Point};
523 ///
524 /// let mut vec = Vec::new();
525 /// vec.push(Point::new(10., 60.));
526 /// vec.push(Point::new(135., 68.));
527 /// vec.push(Point::new(94., 48.));
528 /// vec.push(Point::new(126., 31.));
529 /// vec.push(Point::new(280., 19.));
530 /// vec.push(Point::new(117., 48.));
531 /// vec.push(Point::new(300., 40.));
532 /// vec.push(Point::new(301., 10.));
533 /// let linestring = LineString::from(vec);
534 /// let mut compare = Vec::new();
535 /// compare.push(Point::new(10., 60.));
536 /// compare.push(Point::new(126., 31.));
537 /// compare.push(Point::new(280., 19.));
538 /// compare.push(Point::new(117., 48.));
539 /// compare.push(Point::new(300., 40.));
540 /// compare.push(Point::new(301., 10.));
541 /// let ls_compare = LineString::from(compare);
542 /// let simplified = linestring.simplifyvw_preserve(&668.6);
543 /// assert_eq!(simplified, ls_compare)
544 /// ```
simplifyvw_preserve(&self, epsilon: &T) -> Self where T: Float + RTreeNum545 fn simplifyvw_preserve(&self, epsilon: &T) -> Self
546 where
547 T: Float + RTreeNum;
548 }
549
550 impl<T> SimplifyVWPreserve<T> for LineString<T>
551 where
552 T: Float + RTreeNum,
553 {
simplifyvw_preserve(&self, epsilon: &T) -> LineString<T>554 fn simplifyvw_preserve(&self, epsilon: &T) -> LineString<T> {
555 let gt = GeomSettings {
556 initial_min: 2,
557 min_points: 4,
558 geomtype: GeomType::Line,
559 };
560 let mut simplified = vwp_wrapper(>, self, None, epsilon);
561 LineString::from(simplified.pop().unwrap())
562 }
563 }
564
565 impl<T> SimplifyVWPreserve<T> for MultiLineString<T>
566 where
567 T: Float + RTreeNum,
568 {
simplifyvw_preserve(&self, epsilon: &T) -> MultiLineString<T>569 fn simplifyvw_preserve(&self, epsilon: &T) -> MultiLineString<T> {
570 MultiLineString(
571 self.0
572 .iter()
573 .map(|l| l.simplifyvw_preserve(epsilon))
574 .collect(),
575 )
576 }
577 }
578
579 impl<T> SimplifyVWPreserve<T> for Polygon<T>
580 where
581 T: Float + RTreeNum,
582 {
simplifyvw_preserve(&self, epsilon: &T) -> Polygon<T>583 fn simplifyvw_preserve(&self, epsilon: &T) -> Polygon<T> {
584 let gt = GeomSettings {
585 initial_min: 4,
586 min_points: 6,
587 geomtype: GeomType::Ring,
588 };
589 let mut simplified = vwp_wrapper(>, self.exterior(), Some(self.interiors()), epsilon);
590 let exterior = LineString::from(simplified.remove(0));
591 let interiors = simplified.into_iter().map(LineString::from).collect();
592 Polygon::new(exterior, interiors)
593 }
594 }
595
596 impl<T> SimplifyVWPreserve<T> for MultiPolygon<T>
597 where
598 T: Float + RTreeNum,
599 {
simplifyvw_preserve(&self, epsilon: &T) -> MultiPolygon<T>600 fn simplifyvw_preserve(&self, epsilon: &T) -> MultiPolygon<T> {
601 MultiPolygon(
602 self.0
603 .iter()
604 .map(|p| p.simplifyvw_preserve(epsilon))
605 .collect(),
606 )
607 }
608 }
609
610 impl<T> SimplifyVW<T> for LineString<T>
611 where
612 T: Float,
613 {
simplifyvw(&self, epsilon: &T) -> LineString<T>614 fn simplifyvw(&self, epsilon: &T) -> LineString<T> {
615 LineString::from(visvalingam(self, epsilon))
616 }
617 }
618
619 impl<T> SimplifyVwIdx<T> for LineString<T>
620 where
621 T: Float,
622 {
simplifyvw_idx(&self, epsilon: &T) -> Vec<usize>623 fn simplifyvw_idx(&self, epsilon: &T) -> Vec<usize> {
624 visvalingam_indices(self, epsilon)
625 }
626 }
627
628 impl<T> SimplifyVW<T> for MultiLineString<T>
629 where
630 T: Float,
631 {
simplifyvw(&self, epsilon: &T) -> MultiLineString<T>632 fn simplifyvw(&self, epsilon: &T) -> MultiLineString<T> {
633 MultiLineString(self.0.iter().map(|l| l.simplifyvw(epsilon)).collect())
634 }
635 }
636
637 impl<T> SimplifyVW<T> for Polygon<T>
638 where
639 T: Float,
640 {
simplifyvw(&self, epsilon: &T) -> Polygon<T>641 fn simplifyvw(&self, epsilon: &T) -> Polygon<T> {
642 Polygon::new(
643 self.exterior().simplifyvw(epsilon),
644 self.interiors()
645 .iter()
646 .map(|l| l.simplifyvw(epsilon))
647 .collect(),
648 )
649 }
650 }
651
652 impl<T> SimplifyVW<T> for MultiPolygon<T>
653 where
654 T: Float,
655 {
simplifyvw(&self, epsilon: &T) -> MultiPolygon<T>656 fn simplifyvw(&self, epsilon: &T) -> MultiPolygon<T> {
657 MultiPolygon(self.0.iter().map(|p| p.simplifyvw(epsilon)).collect())
658 }
659 }
660
661 #[cfg(test)]
662 mod test {
663 use super::{
664 cartesian_intersect, visvalingam, vwp_wrapper, GeomSettings, GeomType, SimplifyVW,
665 SimplifyVWPreserve,
666 };
667 use crate::{
668 line_string, point, polygon, Coordinate, LineString, MultiLineString, MultiPolygon, Point,
669 Polygon,
670 };
671
672 #[test]
visvalingam_test()673 fn visvalingam_test() {
674 // this is the PostGIS example
675 let ls = line_string![
676 (x: 5.0, y: 2.0),
677 (x: 3.0, y: 8.0),
678 (x: 6.0, y: 20.0),
679 (x: 7.0, y: 25.0),
680 (x: 10.0, y: 10.0)
681 ];
682
683 let correct = vec![(5.0, 2.0), (7.0, 25.0), (10.0, 10.0)];
684 let correct_ls: Vec<_> = correct
685 .iter()
686 .map(|e| Coordinate::from((e.0, e.1)))
687 .collect();
688
689 let simplified = visvalingam(&ls, &30.);
690 assert_eq!(simplified, correct_ls);
691 }
692 #[test]
vwp_intersection_test()693 fn vwp_intersection_test() {
694 // does the intersection check always work
695 let a = point!(x: 1., y: 3.);
696 let b = point!(x: 3., y: 1.);
697 let c = point!(x: 3., y: 3.);
698 let d = point!(x: 1., y: 1.);
699 // cw + ccw
700 assert_eq!(cartesian_intersect(a, b, c, d), true);
701 // ccw + ccw
702 assert_eq!(cartesian_intersect(b, a, c, d), true);
703 // cw + cw
704 assert_eq!(cartesian_intersect(a, b, d, c), true);
705 // ccw + cw
706 assert_eq!(cartesian_intersect(b, a, d, c), true);
707 }
708 #[test]
simple_vwp_test()709 fn simple_vwp_test() {
710 // this LineString will have a self-intersection if the point with the
711 // smallest associated area is removed
712 // the associated triangle is (1, 2, 3), and has an area of 668.5
713 // the new triangle (0, 1, 3) self-intersects with triangle (3, 4, 5)
714 // Point 1 must also be removed giving a final, valid
715 // LineString of (0, 3, 4, 5, 6, 7)
716 let ls = line_string![
717 (x: 10., y:60.),
718 (x: 135., y: 68.),
719 (x: 94., y: 48.),
720 (x: 126., y: 31.),
721 (x: 280., y: 19.),
722 (x: 117., y: 48.),
723 (x: 300., y: 40.),
724 (x: 301., y: 10.)
725 ];
726 let gt = &GeomSettings {
727 initial_min: 2,
728 min_points: 4,
729 geomtype: GeomType::Line,
730 };
731 let simplified = vwp_wrapper(>, &ls, None, &668.6);
732 // this is the correct, non-intersecting LineString
733 let correct = vec![
734 (10., 60.),
735 (126., 31.),
736 (280., 19.),
737 (117., 48.),
738 (300., 40.),
739 (301., 10.),
740 ];
741 let correct_ls: Vec<_> = correct
742 .iter()
743 .map(|e| Coordinate::from((e.0, e.1)))
744 .collect();
745 assert_eq!(simplified[0], correct_ls);
746 }
747 #[test]
retained_vwp_test()748 fn retained_vwp_test() {
749 // we would expect outer[2] to be removed, as its associated area
750 // is below epsilon. However, this causes a self-intersection
751 // with the inner ring, which would also trigger removal of outer[1],
752 // leaving the geometry below min_points. It is thus retained.
753 // Inner should also be reduced, but has points == initial_min for the Polygon type
754 let outer = line_string![
755 (x: -54.4921875, y: 21.289374355860424),
756 (x: -33.5, y: 56.9449741808516),
757 (x: -22.5, y: 44.08758502824516),
758 (x: -19.5, y: 23.241346102386135),
759 (x: -54.4921875, y: 21.289374355860424)
760 ];
761 let inner = line_string![
762 (x: -24.451171875, y: 35.266685523707665),
763 (x: -29.513671875, y: 47.32027765985069),
764 (x: -22.869140625, y: 43.80817468459856),
765 (x: -24.451171875, y: 35.266685523707665)
766 ];
767 let poly = Polygon::new(outer.clone(), vec![inner]);
768 let simplified = poly.simplifyvw_preserve(&95.4);
769 assert_eq!(simplified.exterior(), &outer);
770 }
771 #[test]
remove_inner_point_vwp_test()772 fn remove_inner_point_vwp_test() {
773 // we would expect outer[2] to be removed, as its associated area
774 // is below epsilon. However, this causes a self-intersection
775 // with the inner ring, which would also trigger removal of outer[1],
776 // leaving the geometry below min_points. It is thus retained.
777 // Inner should be reduced to four points by removing inner[2]
778 let outer = line_string![
779 (x: -54.4921875, y: 21.289374355860424),
780 (x: -33.5, y: 56.9449741808516),
781 (x: -22.5, y: 44.08758502824516),
782 (x: -19.5, y: 23.241346102386135),
783 (x: -54.4921875, y: 21.289374355860424)
784 ];
785 let inner = line_string![
786 (x: -24.451171875, y: 35.266685523707665),
787 (x: -40.0, y: 45.),
788 (x: -29.513671875, y: 47.32027765985069),
789 (x: -22.869140625, y: 43.80817468459856),
790 (x: -24.451171875, y: 35.266685523707665)
791 ];
792 let correct_inner = line_string![
793 (x: -24.451171875, y: 35.266685523707665),
794 (x: -40.0, y: 45.0),
795 (x: -22.869140625, y: 43.80817468459856),
796 (x: -24.451171875, y: 35.266685523707665)
797 ];
798 let poly = Polygon::new(outer.clone(), vec![inner]);
799 let simplified = poly.simplifyvw_preserve(&95.4);
800 assert_eq!(simplified.exterior(), &outer);
801 assert_eq!(simplified.interiors()[0], correct_inner);
802 }
803 #[test]
very_long_vwp_test()804 fn very_long_vwp_test() {
805 // simplify an 8k-point LineString, eliminating self-intersections
806 let points = include!("test_fixtures/norway_main.rs");
807 let points_ls: Vec<_> = points.iter().map(|e| Point::new(e[0], e[1])).collect();
808 let gt = &GeomSettings {
809 initial_min: 2,
810 min_points: 4,
811 geomtype: GeomType::Line,
812 };
813 let simplified = vwp_wrapper(>, &points_ls.into(), None, &0.0005);
814 assert_eq!(simplified[0].len(), 3278);
815 }
816
817 #[test]
visvalingam_test_long()818 fn visvalingam_test_long() {
819 // simplify a longer LineString
820 let points = include!("test_fixtures/vw_orig.rs");
821 let points_ls: LineString<_> = points.iter().map(|e| Point::new(e[0], e[1])).collect();
822 let correct = include!("test_fixtures/vw_simplified.rs");
823 let correct_ls: Vec<_> = correct
824 .iter()
825 .map(|e| Coordinate::from((e[0], e[1])))
826 .collect();
827 let simplified = visvalingam(&points_ls, &0.0005);
828 assert_eq!(simplified, correct_ls);
829 }
830 #[test]
visvalingam_preserve_test_long()831 fn visvalingam_preserve_test_long() {
832 // simplify a longer LineString using the preserve variant
833 let points = include!("test_fixtures/vw_orig.rs");
834 let points_ls: LineString<_> = points.iter().map(|e| Point::new(e[0], e[1])).collect();
835 let correct = include!("test_fixtures/vw_simplified.rs");
836 let correct_ls: Vec<_> = correct.iter().map(|e| Point::new(e[0], e[1])).collect();
837 let simplified = LineString::from(points_ls).simplifyvw_preserve(&0.0005);
838 assert_eq!(simplified, LineString::from(correct_ls));
839 }
840 #[test]
visvalingam_test_empty_linestring()841 fn visvalingam_test_empty_linestring() {
842 let vec: Vec<[f32; 2]> = Vec::new();
843 let compare = Vec::new();
844 let simplified = visvalingam(&LineString::from(vec), &1.0);
845 assert_eq!(simplified, compare);
846 }
847 #[test]
visvalingam_test_two_point_linestring()848 fn visvalingam_test_two_point_linestring() {
849 let mut vec = Vec::new();
850 vec.push(Point::new(0.0, 0.0));
851 vec.push(Point::new(27.8, 0.1));
852 let mut compare = Vec::new();
853 compare.push(Coordinate::from((0.0, 0.0)));
854 compare.push(Coordinate::from((27.8, 0.1)));
855 let simplified = visvalingam(&LineString::from(vec), &1.0);
856 assert_eq!(simplified, compare);
857 }
858
859 #[test]
multilinestring()860 fn multilinestring() {
861 // this is the PostGIS example
862 let points = vec![
863 (5.0, 2.0),
864 (3.0, 8.0),
865 (6.0, 20.0),
866 (7.0, 25.0),
867 (10.0, 10.0),
868 ];
869 let points_ls: Vec<_> = points.iter().map(|e| Point::new(e.0, e.1)).collect();
870
871 let correct = vec![(5.0, 2.0), (7.0, 25.0), (10.0, 10.0)];
872 let correct_ls: Vec<_> = correct.iter().map(|e| Point::new(e.0, e.1)).collect();
873
874 let mline = MultiLineString(vec![LineString::from(points_ls)]);
875 assert_eq!(
876 mline.simplifyvw(&30.),
877 MultiLineString(vec![LineString::from(correct_ls)])
878 );
879 }
880
881 #[test]
polygon()882 fn polygon() {
883 let poly = polygon![
884 (x: 0., y: 0.),
885 (x: 0., y: 10.),
886 (x: 5., y: 11.),
887 (x: 10., y: 10.),
888 (x: 10., y: 0.),
889 (x: 0., y: 0.),
890 ];
891
892 let poly2 = poly.simplifyvw(&10.);
893
894 assert_eq!(
895 poly2,
896 polygon![
897 (x: 0., y: 0.),
898 (x: 0., y: 10.),
899 (x: 10., y: 10.),
900 (x: 10., y: 0.),
901 (x: 0., y: 0.),
902 ],
903 );
904 }
905
906 #[test]
multipolygon()907 fn multipolygon() {
908 let mpoly = MultiPolygon(vec![Polygon::new(
909 LineString::from(vec![
910 (0., 0.),
911 (0., 10.),
912 (5., 11.),
913 (10., 10.),
914 (10., 0.),
915 (0., 0.),
916 ]),
917 vec![],
918 )]);
919
920 let mpoly2 = mpoly.simplifyvw(&10.);
921
922 assert_eq!(
923 mpoly2,
924 MultiPolygon(vec![Polygon::new(
925 LineString::from(vec![(0., 0.), (0., 10.), (10., 10.), (10., 0.), (0., 0.)]),
926 vec![],
927 )])
928 );
929 }
930 }
931