geo/algorithm/
simplify.rs

1use crate::{Coord, GeoFloat, Line, LineString, MultiLineString, MultiPolygon, Polygon};
2use crate::{CoordsIter, EuclideanDistance};
3
4const LINE_STRING_INITIAL_MIN: usize = 2;
5const POLYGON_INITIAL_MIN: usize = 4;
6
7// Because the RDP algorithm is recursive, we can't assign an index to a point inside the loop
8// instead, we wrap a simple struct around index and point in a wrapper function,
9// passing that around instead, extracting either points or indices on the way back out
10#[derive(Copy, Clone)]
11struct RdpIndex<T>
12where
13    T: GeoFloat,
14{
15    index: usize,
16    coord: Coord<T>,
17}
18
19// Wrapper for the RDP algorithm, returning simplified points
20fn rdp<T, I: Iterator<Item = Coord<T>>, const INITIAL_MIN: usize>(
21    coords: I,
22    epsilon: &T,
23) -> Vec<Coord<T>>
24where
25    T: GeoFloat,
26{
27    // Epsilon must be greater than zero for any meaningful simplification to happen
28    if *epsilon <= T::zero() {
29        return coords.collect::<Vec<Coord<T>>>();
30    }
31    let rdp_indices = &coords
32        .enumerate()
33        .map(|(idx, coord)| RdpIndex { index: idx, coord })
34        .collect::<Vec<RdpIndex<T>>>();
35    let mut simplified_len = rdp_indices.len();
36    let simplified_coords: Vec<_> =
37        compute_rdp::<T, INITIAL_MIN>(rdp_indices, &mut simplified_len, epsilon)
38            .into_iter()
39            .map(|rdpindex| rdpindex.coord)
40            .collect();
41    debug_assert_eq!(simplified_coords.len(), simplified_len);
42    simplified_coords
43}
44
45// Wrapper for the RDP algorithm, returning simplified point indices
46fn calculate_rdp_indices<T, const INITIAL_MIN: usize>(
47    rdp_indices: &[RdpIndex<T>],
48    epsilon: &T,
49) -> Vec<usize>
50where
51    T: GeoFloat,
52{
53    if *epsilon <= T::zero() {
54        return rdp_indices
55            .iter()
56            .map(|rdp_index| rdp_index.index)
57            .collect();
58    }
59
60    let mut simplified_len = rdp_indices.len();
61    let simplified_coords =
62        compute_rdp::<T, INITIAL_MIN>(rdp_indices, &mut simplified_len, epsilon)
63            .into_iter()
64            .map(|rdpindex| rdpindex.index)
65            .collect::<Vec<usize>>();
66    debug_assert_eq!(simplified_len, simplified_coords.len());
67    simplified_coords
68}
69
70// Ramer–Douglas-Peucker line simplification algorithm
71// This function returns both the retained points, and their indices in the original geometry,
72// for more flexible use by FFI implementers
73fn compute_rdp<T, const INITIAL_MIN: usize>(
74    rdp_indices: &[RdpIndex<T>],
75    simplified_len: &mut usize,
76    epsilon: &T,
77) -> Vec<RdpIndex<T>>
78where
79    T: GeoFloat,
80{
81    if rdp_indices.is_empty() {
82        return vec![];
83    }
84
85    let first = rdp_indices[0];
86    let last = rdp_indices[rdp_indices.len() - 1];
87    if rdp_indices.len() == 2 {
88        return vec![first, last];
89    }
90
91    let first_last_line = Line::new(first.coord, last.coord);
92
93    // Find the farthest `RdpIndex` from `first_last_line`
94    let (farthest_index, farthest_distance) = rdp_indices
95        .iter()
96        .enumerate()
97        .take(rdp_indices.len() - 1) // Don't include the last index
98        .skip(1) // Don't include the first index
99        .map(|(index, rdp_index)| (index, rdp_index.coord.euclidean_distance(&first_last_line)))
100        .fold(
101            (0usize, T::zero()),
102            |(farthest_index, farthest_distance), (index, distance)| {
103                if distance >= farthest_distance {
104                    (index, distance)
105                } else {
106                    (farthest_index, farthest_distance)
107                }
108            },
109        );
110    debug_assert_ne!(farthest_index, 0);
111
112    if farthest_distance > *epsilon {
113        // The farthest index was larger than epsilon, so we will recursively simplify subsegments
114        // split by the farthest index.
115        let mut intermediate =
116            compute_rdp::<T, INITIAL_MIN>(&rdp_indices[..=farthest_index], simplified_len, epsilon);
117
118        intermediate.pop(); // Don't include the farthest index twice
119
120        intermediate.extend_from_slice(&compute_rdp::<T, INITIAL_MIN>(
121            &rdp_indices[farthest_index..],
122            simplified_len,
123            epsilon,
124        ));
125        return intermediate;
126    }
127
128    // The farthest index was less than or equal to epsilon, so we will retain only the first
129    // and last indices, resulting in the indices inbetween getting culled.
130
131    // Update `simplified_len` to reflect the new number of indices by subtracting the number
132    // of indices we're culling.
133    let number_culled = rdp_indices.len() - 2;
134    let new_length = *simplified_len - number_culled;
135
136    // If `simplified_len` is now lower than the minimum number of indices needed, then don't
137    // perform the culling and return the original input.
138    if new_length < INITIAL_MIN {
139        return rdp_indices.to_owned();
140    }
141    *simplified_len = new_length;
142
143    // Cull indices between `first` and `last`.
144    vec![first, last]
145}
146
147/// Simplifies a geometry.
148///
149/// The [Ramer–Douglas–Peucker
150/// algorithm](https://en.wikipedia.org/wiki/Ramer–Douglas–Peucker_algorithm) simplifies a
151/// linestring. Polygons are simplified by running the RDP algorithm on all their constituent
152/// rings. This may result in invalid Polygons, and has no guarantee of preserving topology.
153///
154/// Multi* objects are simplified by simplifying all their constituent geometries individually.
155///
156/// An epsilon less than or equal to zero will return an unaltered version of the geometry.
157pub trait Simplify<T, Epsilon = T> {
158    /// Returns the simplified representation of a geometry, using the [Ramer–Douglas–Peucker](https://en.wikipedia.org/wiki/Ramer–Douglas–Peucker_algorithm) algorithm
159    ///
160    /// # Examples
161    ///
162    /// ```
163    /// use geo::Simplify;
164    /// use geo::line_string;
165    ///
166    /// let line_string = line_string![
167    ///     (x: 0.0, y: 0.0),
168    ///     (x: 5.0, y: 4.0),
169    ///     (x: 11.0, y: 5.5),
170    ///     (x: 17.3, y: 3.2),
171    ///     (x: 27.8, y: 0.1),
172    /// ];
173    ///
174    /// let simplified = line_string.simplify(&1.0);
175    ///
176    /// let expected = line_string![
177    ///     (x: 0.0, y: 0.0),
178    ///     (x: 5.0, y: 4.0),
179    ///     (x: 11.0, y: 5.5),
180    ///     (x: 27.8, y: 0.1),
181    /// ];
182    ///
183    /// assert_eq!(expected, simplified)
184    /// ```
185    fn simplify(&self, epsilon: &T) -> Self
186    where
187        T: GeoFloat;
188}
189
190/// Simplifies a geometry, returning the retained _indices_ of the input.
191///
192/// This operation uses the [Ramer–Douglas–Peucker algorithm](https://en.wikipedia.org/wiki/Ramer–Douglas–Peucker_algorithm)
193/// and does not guarantee that the returned geometry is valid.
194///
195/// An epsilon less than or equal to zero will return an unaltered version of the geometry.
196pub trait SimplifyIdx<T, Epsilon = T> {
197    /// Returns the simplified indices of a geometry, using the [Ramer–Douglas–Peucker](https://en.wikipedia.org/wiki/Ramer–Douglas–Peucker_algorithm) algorithm
198    ///
199    /// # Examples
200    ///
201    /// ```
202    /// use geo::SimplifyIdx;
203    /// use geo::line_string;
204    ///
205    /// let line_string = line_string![
206    ///     (x: 0.0, y: 0.0),
207    ///     (x: 5.0, y: 4.0),
208    ///     (x: 11.0, y: 5.5),
209    ///     (x: 17.3, y: 3.2),
210    ///     (x: 27.8, y: 0.1),
211    /// ];
212    ///
213    /// let simplified = line_string.simplify_idx(&1.0);
214    ///
215    /// let expected = vec![
216    ///     0_usize,
217    ///     1_usize,
218    ///     2_usize,
219    ///     4_usize,
220    /// ];
221    ///
222    /// assert_eq!(expected, simplified);
223    /// ```
224    fn simplify_idx(&self, epsilon: &T) -> Vec<usize>
225    where
226        T: GeoFloat;
227}
228
229impl<T> Simplify<T> for LineString<T>
230where
231    T: GeoFloat,
232{
233    fn simplify(&self, epsilon: &T) -> Self {
234        LineString::from(rdp::<_, _, LINE_STRING_INITIAL_MIN>(
235            self.coords_iter(),
236            epsilon,
237        ))
238    }
239}
240
241impl<T> SimplifyIdx<T> for LineString<T>
242where
243    T: GeoFloat,
244{
245    fn simplify_idx(&self, epsilon: &T) -> Vec<usize> {
246        calculate_rdp_indices::<_, LINE_STRING_INITIAL_MIN>(
247            &self
248                .0
249                .iter()
250                .enumerate()
251                .map(|(idx, coord)| RdpIndex {
252                    index: idx,
253                    coord: *coord,
254                })
255                .collect::<Vec<RdpIndex<T>>>(),
256            epsilon,
257        )
258    }
259}
260
261impl<T> Simplify<T> for MultiLineString<T>
262where
263    T: GeoFloat,
264{
265    fn simplify(&self, epsilon: &T) -> Self {
266        MultiLineString::new(self.iter().map(|l| l.simplify(epsilon)).collect())
267    }
268}
269
270impl<T> Simplify<T> for Polygon<T>
271where
272    T: GeoFloat,
273{
274    fn simplify(&self, epsilon: &T) -> Self {
275        Polygon::new(
276            LineString::from(rdp::<_, _, POLYGON_INITIAL_MIN>(
277                self.exterior().coords_iter(),
278                epsilon,
279            )),
280            self.interiors()
281                .iter()
282                .map(|l| {
283                    LineString::from(rdp::<_, _, POLYGON_INITIAL_MIN>(l.coords_iter(), epsilon))
284                })
285                .collect(),
286        )
287    }
288}
289
290impl<T> Simplify<T> for MultiPolygon<T>
291where
292    T: GeoFloat,
293{
294    fn simplify(&self, epsilon: &T) -> Self {
295        MultiPolygon::new(self.iter().map(|p| p.simplify(epsilon)).collect())
296    }
297}
298
299#[cfg(test)]
300mod test {
301    use super::*;
302    use crate::geo_types::coord;
303    use crate::{line_string, polygon};
304
305    #[test]
306    fn recursion_test() {
307        let input = [
308            coord! { x: 8.0, y: 100.0 },
309            coord! { x: 9.0, y: 100.0 },
310            coord! { x: 12.0, y: 100.0 },
311        ];
312        let actual = rdp::<_, _, 2>(input.into_iter(), &1.0);
313        let expected = [coord! { x: 8.0, y: 100.0 }, coord! { x: 12.0, y: 100.0 }];
314        assert_eq!(actual, expected);
315    }
316
317    #[test]
318    fn rdp_test() {
319        let vec = vec![
320            coord! { x: 0.0, y: 0.0 },
321            coord! { x: 5.0, y: 4.0 },
322            coord! { x: 11.0, y: 5.5 },
323            coord! { x: 17.3, y: 3.2 },
324            coord! { x: 27.8, y: 0.1 },
325        ];
326        let compare = vec![
327            coord! { x: 0.0, y: 0.0 },
328            coord! { x: 5.0, y: 4.0 },
329            coord! { x: 11.0, y: 5.5 },
330            coord! { x: 27.8, y: 0.1 },
331        ];
332        let simplified = rdp::<_, _, 2>(vec.into_iter(), &1.0);
333        assert_eq!(simplified, compare);
334    }
335    #[test]
336    fn rdp_test_empty_linestring() {
337        let vec = Vec::new();
338        let compare = Vec::new();
339        let simplified = rdp::<_, _, 2>(vec.into_iter(), &1.0);
340        assert_eq!(simplified, compare);
341    }
342    #[test]
343    fn rdp_test_two_point_linestring() {
344        let vec = vec![coord! { x: 0.0, y: 0.0 }, coord! { x: 27.8, y: 0.1 }];
345        let compare = vec![coord! { x: 0.0, y: 0.0 }, coord! { x: 27.8, y: 0.1 }];
346        let simplified = rdp::<_, _, 2>(vec.into_iter(), &1.0);
347        assert_eq!(simplified, compare);
348    }
349
350    #[test]
351    fn multilinestring() {
352        let mline = MultiLineString::new(vec![LineString::from(vec![
353            (0.0, 0.0),
354            (5.0, 4.0),
355            (11.0, 5.5),
356            (17.3, 3.2),
357            (27.8, 0.1),
358        ])]);
359
360        let mline2 = mline.simplify(&1.0);
361
362        assert_eq!(
363            mline2,
364            MultiLineString::new(vec![LineString::from(vec![
365                (0.0, 0.0),
366                (5.0, 4.0),
367                (11.0, 5.5),
368                (27.8, 0.1),
369            ])])
370        );
371    }
372
373    #[test]
374    fn polygon() {
375        let poly = polygon![
376            (x: 0., y: 0.),
377            (x: 0., y: 10.),
378            (x: 5., y: 11.),
379            (x: 10., y: 10.),
380            (x: 10., y: 0.),
381            (x: 0., y: 0.),
382        ];
383
384        let poly2 = poly.simplify(&2.);
385
386        assert_eq!(
387            poly2,
388            polygon![
389                (x: 0., y: 0.),
390                (x: 0., y: 10.),
391                (x: 10., y: 10.),
392                (x: 10., y: 0.),
393                (x: 0., y: 0.),
394            ],
395        );
396    }
397
398    #[test]
399    fn multipolygon() {
400        let mpoly = MultiPolygon::new(vec![polygon![
401            (x: 0., y: 0.),
402            (x: 0., y: 10.),
403            (x: 5., y: 11.),
404            (x: 10., y: 10.),
405            (x: 10., y: 0.),
406            (x: 0., y: 0.),
407        ]]);
408
409        let mpoly2 = mpoly.simplify(&2.);
410
411        assert_eq!(
412            mpoly2,
413            MultiPolygon::new(vec![polygon![
414                (x: 0., y: 0.),
415                (x: 0., y: 10.),
416                (x: 10., y: 10.),
417                (x: 10., y: 0.),
418                (x: 0., y: 0.)
419            ]]),
420        );
421    }
422
423    #[test]
424    fn simplify_negative_epsilon() {
425        let ls = line_string![
426            (x: 0., y: 0.),
427            (x: 0., y: 10.),
428            (x: 5., y: 11.),
429            (x: 10., y: 10.),
430            (x: 10., y: 0.),
431        ];
432        let simplified = ls.simplify(&-1.0);
433        assert_eq!(ls, simplified);
434    }
435
436    #[test]
437    fn simplify_idx_negative_epsilon() {
438        let ls = line_string![
439            (x: 0., y: 0.),
440            (x: 0., y: 10.),
441            (x: 5., y: 11.),
442            (x: 10., y: 10.),
443            (x: 10., y: 0.),
444        ];
445        let indices = ls.simplify_idx(&-1.0);
446        assert_eq!(vec![0usize, 1, 2, 3, 4], indices);
447    }
448
449    // https://github.com/georust/geo/issues/142
450    #[test]
451    fn simplify_line_string_polygon_initial_min() {
452        let ls = line_string![
453            ( x: 1.4324054e-16, y: 1.4324054e-16 ),
454            ( x: 1.4324054e-16, y: 1.4324054e-16 ),
455            ( x: -5.9730447e26, y: 1.5590374e-27 ),
456            ( x: 1.4324054e-16, y: 1.4324054e-16 ),
457        ];
458        let epsilon: f64 = 3.46e-43;
459
460        // LineString result should be three coordinates
461        let result = ls.simplify(&epsilon);
462        assert_eq!(
463            line_string![
464                ( x: 1.4324054e-16, y: 1.4324054e-16 ),
465                ( x: -5.9730447e26, y: 1.5590374e-27 ),
466                ( x: 1.4324054e-16, y: 1.4324054e-16 ),
467            ],
468            result
469        );
470
471        // Polygon result should be five coordinates
472        let result = Polygon::new(ls, vec![]).simplify(&epsilon);
473        assert_eq!(
474            polygon![
475                ( x: 1.4324054e-16, y: 1.4324054e-16 ),
476                ( x: 1.4324054e-16, y: 1.4324054e-16 ),
477                ( x: -5.9730447e26, y: 1.5590374e-27 ),
478                ( x: 1.4324054e-16, y: 1.4324054e-16 ),
479            ],
480            result,
481        );
482    }
483
484    // https://github.com/georust/geo/issues/995
485    #[test]
486    fn dont_oversimplify() {
487        let unsimplified = line_string![
488            (x: 0.0, y: 0.0),
489            (x: 5.0, y: 4.0),
490            (x: 11.0, y: 5.5),
491            (x: 17.3, y: 3.2),
492            (x: 27.8, y: 0.1)
493        ];
494        let actual = unsimplified.simplify(&30.0);
495        let expected = line_string![
496            (x: 0.0, y: 0.0),
497            (x: 27.8, y: 0.1)
498        ];
499        assert_eq!(actual, expected);
500    }
501}
geo/algorithm/simplify.rs

geo/algorithm/
simplify.rs
