Robust Winding Order

Compute Winding Order using robust predicates.

+ derive `PartialOrd` for `Coordinate`
+ add `is_closed` utility to `LineString`
+ add `robust-0.2.2` crate
+ add `RobustWindingOrder` trait
+ impl `RobustWindingOrder` for `LineString<T: Float>`
+ add tests (copied from `WindingOrder`)

Author: rmanoka

geo-types/src/coordinate.rs

@@ -9,7 +9,7 @@ use approx::{AbsDiffEq, RelativeEq, UlpsEq};
 /// as an envelope, a precision model, and spatial reference system
 /// information), a `Coordinate` only contains ordinate values and accessor
 /// methods.
-#[derive(Eq, PartialEq, Clone, Copy, Debug, Hash)]
+#[derive(Eq, PartialEq, PartialOrd, Clone, Copy, Debug, Hash)]
 #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
 pub struct Coordinate<T>
 where

geo-types/src/line_string.rs

@@ -160,10 +160,8 @@ impl<T: CoordinateType> LineString<T> {
     /// and the value of the first coordinate does not equal the value of the last coordinate, then
     /// a new coordinate is added to the end with the value of the first coordinate.
     pub(crate) fn close(&mut self) {
-        if let (Some(first), Some(last)) = (self.0.first().copied(), self.0.last().copied()) {
-            if first != last {
-                self.0.push(first);
-            }
+        if !self.is_closed() && !self.0.is_empty() {
+            self.0.push(self.0[0]);
         }
     }
 
@@ -181,6 +179,29 @@ impl<T: CoordinateType> LineString<T> {
     pub fn num_coords(&self) -> usize {
         self.0.len()
     }
+
+    /// Checks if the linestring is closed; i.e. the first
+    /// and last points have the same coords.
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// use geo_types::LineString;
+    ///
+    /// let mut coords = vec![(0., 0.), (5., 0.), (0., 0.)];
+    /// let line_string: LineString<f32> = coords.into_iter().collect();
+    /// assert!(line_string.is_closed());
+    /// ```
+    pub fn is_closed(&self) -> bool {
+        // LineString with less than 2 elements can't be closed
+        if self.0.len() < 2 {
+            false
+        } else if self.0.first().unwrap() != self.0.last().unwrap() {
+            false
+        } else {
+            true
+        }
+    }
 }
 
 /// Turn a `Vec` of `Point`-like objects into a `LineString`.

geo/Cargo.toml

@@ -24,6 +24,8 @@ proj = { version = "0.20.3", optional = true }
 
 geo-types = { version = "0.6.0", path = "../geo-types", features = ["rstar"] }
 
+robust = { version = "0.2.2" }
+
 [features]
 default = []
 use-proj = ["proj"]

geo/src/algorithm/mod.rs

@@ -59,6 +59,8 @@ pub mod vincenty_distance;
 pub mod vincenty_length;
 /// Calculate and work with the winding order of `Linestring`s.
 pub mod winding_order;
+/// Calculate the winding order of `LineString`s using robust predicates.
+pub mod robust_winding_order;
 /// Locate a point along a `Line` or `LineString`.
 pub mod line_locate_point;
 /// Interpolate a point along a `Line` or `LineString`.

geo/src/algorithm/robust_winding_order.rs

@@ -0,0 +1,159 @@
+use super::winding_order::*;
+use crate::{Coordinate, CoordinateType, LineString};
+use num_traits::Float;
+
+pub trait RobustWinding<T>: Winding<T>
+where
+    T: CoordinateType,
+{
+    /// Return the winding order of this object
+    fn robust_winding_order(&self) -> Option<WindingOrder>;
+}
+
+/// Compute index of the lexicographically least point.
+/// Should only be called on a non-empty slice.
+fn lexicographically_least_index<T: Copy + PartialOrd>(pts: &[T]) -> usize {
+    assert!(pts.len() > 0);
+
+    let mut min: Option<(usize, T)> = None;
+    for (i, pt) in pts.iter().enumerate() {
+        if let Some((_, min_pt)) = min {
+            if pt < &min_pt {
+                min = Some( (i, *pt) )
+            }
+        } else {
+            min = Some( (i, *pt) )
+        }
+    }
+
+    min.unwrap().0
+}
+
+impl<T: CoordinateType + Float> RobustWinding<T> for LineString<T> {
+    fn robust_winding_order(&self) -> Option<WindingOrder> {
+        // If linestring has at most 2 points, it is either
+        // not closed, or is the same point. Either way, the
+        // WindingOrder is unspecified.
+        if self.num_coords() < 3 { return None; }
+
+        // Open linestrings do not have a winding order.
+        if !self.is_closed() { return None; }
+
+        let i = lexicographically_least_index(&self.0);
+        let mut next = i;
+        while self.0[next] == self.0[i] {
+            next += 1;
+            if next >= self.num_coords() { next = 0; }
+        }
+
+        let mut prev = i;
+        while self.0[prev] == self.0[i] {
+            if prev == 0 {
+                prev = self.num_coords() - 1;
+            } else {
+                prev -= 1;
+            }
+        }
+
+        use robust::{Coord, orient2d};
+        use num_traits::NumCast;
+        let orientation = orient2d(
+            Coord {
+                x: <f64 as NumCast>::from( self.0[prev].x ).unwrap(),
+                y: <f64 as NumCast>::from( self.0[prev].y ).unwrap(),
+            },
+            Coord {
+                x: <f64 as NumCast>::from( self.0[i].x ).unwrap(),
+                y: <f64 as NumCast>::from( self.0[i].y ).unwrap(),
+            },
+            Coord {
+                x: <f64 as NumCast>::from( self.0[next].x ).unwrap(),
+                y: <f64 as NumCast>::from( self.0[next].y ).unwrap(),
+            },
+        );
+
+        if orientation < 0. {
+            Some(WindingOrder::Clockwise)
+        } else if orientation > 0. {
+            Some(WindingOrder::CounterClockwise)
+        } else {
+            None
+        }
+
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use crate::Point;
+
+    #[test]
+    fn robust_winding_order() {
+        // 3 points forming a triangle
+        let a = Point::new(0., 0.);
+        let b = Point::new(2., 0.);
+        let c = Point::new(1., 2.);
+
+        // That triangle, but in clockwise ordering
+        let cw_line = LineString::from(vec![a.0, c.0, b.0, a.0]);
+        // That triangle, but in counterclockwise ordering
+        let ccw_line = LineString::from(vec![a.0, b.0, c.0, a.0]);
+
+        // Verify open linestrings return None
+        assert!(LineString::from(vec![a.0, b.0, c.0])
+            .robust_winding_order()
+            .is_none());
+
+        assert_eq!(cw_line.robust_winding_order(), Some(WindingOrder::Clockwise));
+        assert_eq!(cw_line.is_cw(), true);
+        assert_eq!(cw_line.is_ccw(), false);
+        assert_eq!(
+            ccw_line.robust_winding_order(),
+            Some(WindingOrder::CounterClockwise)
+        );
+        assert_eq!(ccw_line.is_cw(), false);
+        assert_eq!(ccw_line.is_ccw(), true);
+
+        let cw_points1: Vec<_> = cw_line.points_cw().collect();
+        assert_eq!(cw_points1.len(), 4);
+        assert_eq!(cw_points1[0], a);
+        assert_eq!(cw_points1[1], c);
+        assert_eq!(cw_points1[2], b);
+        assert_eq!(cw_points1[3], a);
+
+        let ccw_points1: Vec<_> = cw_line.points_ccw().collect();
+        assert_eq!(ccw_points1.len(), 4);
+        assert_eq!(ccw_points1[0], a);
+        assert_eq!(ccw_points1[1], b);
+        assert_eq!(ccw_points1[2], c);
+        assert_eq!(ccw_points1[3], a);
+
+        assert_ne!(cw_points1, ccw_points1);
+
+        let cw_points2: Vec<_> = ccw_line.points_cw().collect();
+        let ccw_points2: Vec<_> = ccw_line.points_ccw().collect();
+
+        // cw_line and ccw_line are wound differently, but the ordered winding iterator should have
+        // make them similar
+        assert_eq!(cw_points2, cw_points2);
+        assert_eq!(ccw_points2, ccw_points2);
+
+        // test make_clockwise_winding
+        let mut new_line1 = ccw_line.clone();
+        new_line1.make_cw_winding();
+        assert_eq!(new_line1.robust_winding_order(), Some(WindingOrder::Clockwise));
+        assert_eq!(new_line1, cw_line);
+        assert_ne!(new_line1, ccw_line);
+
+        // test make_counterclockwise_winding
+        let mut new_line2 = cw_line.clone();
+        new_line2.make_ccw_winding();
+        assert_eq!(
+            new_line2.robust_winding_order(),
+            Some(WindingOrder::CounterClockwise)
+        );
+        assert_ne!(new_line2, cw_line);
+        assert_eq!(new_line2, ccw_line);
+    }
+}