Robust Winding Order

geo-types/src/line_string.rs

@@ -159,11 +159,9 @@ impl<T: CoordinateType> LineString<T> {
     /// Close the `LineString`. Specifically, if the `LineString` has is at least one coordinate,
     /// 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);
-            }
+    pub fn close(&mut self) {
+        if !self.is_closed() && !self.0.is_empty() {
+            self.0.push(self.0[0]);
         }
     }
 
@@ -181,6 +179,23 @@ 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. Note that a
+    /// single point is considered closed.
+    ///
+    /// # 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 {
+        !self.0.is_empty() && self.0.first() == self.0.last()
+    }
 }
 
 /// 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/area.rs

@@ -4,7 +4,43 @@ use crate::{
 };
 use num_traits::Float;
 
-use crate::algorithm::winding_order::twice_signed_ring_area;
+pub(crate) fn twice_signed_ring_area<T>(linestring: &LineString<T>) -> T
+where
+    T: CoordinateType,
+{
+    // LineString with less than 3 points is empty, or a
+    // single point, or is not closed.
+    if linestring.0.len() < 3 {
+        return T::zero();
+    }
+
+    // Above test ensures the vector has at least 2 elements.
+    // We check if linestring is closed, and return 0 otherwise.
+    if linestring.0.first().unwrap() != linestring.0.last().unwrap() {
+        return T::zero();
+    }
+
+    // Use a reasonable shift for the line-string coords
+    // to avoid numerical-errors when summing the
+    // determinants.
+    //
+    // Note: we can't use the `Centroid` trait as it
+    // requries `T: Float` and in fact computes area in the
+    // implementation. Another option is to use the average
+    // of the coordinates, but it is not fool-proof to
+    // divide by the length of the linestring (eg. a long
+    // line-string with T = u8)
+    let shift = linestring.0[0];
+
+    let mut tmp = T::zero();
+    for line in linestring.lines() {
+        use crate::algorithm::map_coords::MapCoords;
+        let line = line.map_coords(|&(x, y)| (x - shift.x, y - shift.y));
+        tmp = tmp + line.determinant();
+    }
+
+    tmp
+}
 
 /// Signed and unsigned planar area of a geometry.
 ///

geo/src/algorithm/graham_hull.rs

@@ -0,0 +1,155 @@
+use crate::{Coordinate, LineString};
+use super::kernels::*;
+
+// Utility function: swap idx to head(0th position), remove
+// head (modifies the slice), and return head as a reference
+fn swap_remove_to_first<'a, T>(slice: &mut &'a mut [T], idx: usize) -> &'a mut T {
+    let tmp = std::mem::replace(slice, &mut []);
+    tmp.swap(0, idx);
+    let (h, t) = tmp.split_first_mut().unwrap();
+    *slice = t;
+    h
+}
+
+/// Graham Scan: https://en.wikipedia.org/wiki/Graham_scan
+pub fn graham_hull<T>(mut points: &mut [Coordinate<T>], include_colinear: bool) -> LineString<T>
+where
+    T: HasKernel,
+{
+    if points.len() < 4 {
+        // Nothing to build with fewer than four points. We
+        // remove repeated points if any, and ensure ccw
+        // invariant.
+        use super::winding_order::Winding;
+        let mut ls: LineString<T> = points
+            .iter()
+            .enumerate()
+            .filter_map(|(i, pt)| {
+                if i == 0 || pt != &points[i - 1] {
+                    Some(*pt)
+                } else {
+                    None
+                }
+            })
+            .collect();
+        ls.close();
+        ls.make_ccw_winding();
+        return ls;
+    }
+
+    // Allocate output vector
+    let mut output = Vec::with_capacity(points.len());
+
+    // Find lexicographically least point
+    use crate::utils::lexicographically_least_index;
+    let min_idx = lexicographically_least_index(points);
+    let head = swap_remove_to_first(&mut points, min_idx);
+    output.push(*head);
+
+    // Sort rest of the points by angle it makes with head
+    // point. If two points are colinear with head, we sort
+    // by distance. We use kernel predicates here.
+    let cmp = |q: &Coordinate<T>, r: &Coordinate<T>| {
+        use std::cmp::Ordering;
+        match T::Ker::orient2d(*q, *head, *r) {
+            Orientation::CounterClockwise => Ordering::Greater,
+            Orientation::Clockwise => Ordering::Less,
+            Orientation::Colinear => {
+                let dist1 = T::Ker::square_euclidean_distance(*head, *q);
+                let dist2 = T::Ker::square_euclidean_distance(*head, *r);
+                dist1.partial_cmp(&dist2).unwrap()
+            }
+        }
+    };
+    points.sort_unstable_by(cmp);
+
+    for pt in points.iter() {
+        while output.len() > 1 {
+            let len = output.len();
+            match T::Ker::orient2d(output[len-2], output[len-1], *pt) {
+                Orientation::CounterClockwise => { break; }
+                Orientation::Clockwise => { output.pop(); }
+                Orientation::Colinear => {
+                    if include_colinear {
+                        break;
+                    } else {
+                        output.pop();
+                    }
+                }
+            }
+        }
+        output.push(*pt);
+    }
+
+    // Close and output the line string
+    let mut output = LineString(output);
+    output.close();
+    output
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use geo_types::CoordinateType;
+
+    fn is_ccw_convex<T: CoordinateType + HasKernel>(mut ls: &[Coordinate<T>]) -> bool {
+        if ls.len() > 1 && ls[0] == ls[ls.len() - 1] {
+            ls = &ls[1..];
+        }
+        let n = ls.len();
+        if n < 3 { return true; }
+
+        ls.iter().enumerate().all(|(i, coord)| {
+            let np = ls[( i + 1 ) % n];
+            let nnp = ls[( i + 2 ) % n];
+            T::Ker::orient2d(*coord, np, nnp) == Orientation::CounterClockwise
+        })
+    }
+
+    fn test_convexity<T: CoordinateType + HasKernel>(initial: &[(T, T)]) {
+        let mut v: Vec<_> = initial.iter().map(|e| Coordinate::from((e.0, e.1))).collect();
+        let hull = graham_hull(&mut v, false);
+        assert!(is_ccw_convex(&hull.0));
+    }
+
+    #[test]
+    fn test_graham_hull_ccw() {
+        let initial = vec![
+            (1.0, 0.0), (2.0, 1.0), (1.75, 1.1),
+            (1.0, 2.0), (0.0, 1.0), (1.0, 0.0),
+        ];
+        test_convexity(&initial);
+    }
+
+    #[test]
+    fn graham_hull_test1() {
+        let v: Vec<_> = vec![
+            (0, 0), (4, 0), (4, 1),
+            (1, 1), (1, 4), (0, 4), (0, 0),
+        ];
+        test_convexity(&v);
+    }
+
+    #[test]
+    fn graham_hull_test2() {
+        let v = vec![
+            (0, 10), (1, 1), (10, 0),
+            (1, -1), (0, -10), (-1, -1),
+            (-10, 0), (-1, 1), (0, 10),
+        ];
+        test_convexity(&v);
+    }
+
+    #[test]
+    fn graham_test_complex() {
+        let v = include!("test_fixtures/poly1.rs");
+        test_convexity(&v);
+
+    }
+
+    #[test]
+    fn quick_hull_test_complex_2() {
+        let coords = include!("test_fixtures/poly2.rs");
+        test_convexity(&coords);
+    }
+}

geo/src/algorithm/kernels/mod.rs

@@ -0,0 +1,75 @@
+use crate::{CoordinateType, Coordinate};
+
+#[derive(Debug, PartialEq, Eq, PartialOrd, Ord)]
+pub enum Orientation {
+    CounterClockwise,
+    Clockwise,
+    Colinear,
+}
+
+/// Kernel trait to provide predicates to operate on
+/// different scalar types.
+pub trait Kernel {
+    type Scalar: CoordinateType;
+
+    /// Gives the orientation of 3 2-dimensional points:
+    /// ccw, cw or colinear (None)
+    fn orient2d(
+        p: Coordinate<Self::Scalar>,
+        q: Coordinate<Self::Scalar>,
+        r: Coordinate<Self::Scalar>,
+    ) -> Orientation {
+
+        let res = (q.x - p.x) * (r.y - q.y) - (q.y - p.y) * (r.x - q.x);
+        use num_traits::Zero;
+        if res > Zero::zero() {
+            Orientation::CounterClockwise
+        } else if res < Zero::zero() {
+            Orientation::Clockwise
+        } else {
+            Orientation::Colinear
+        }
+
+    }
+
+    fn square_euclidean_distance(
+        p: Coordinate<Self::Scalar>,
+        q: Coordinate<Self::Scalar>,
+    ) -> Self::Scalar {
+        (p.x - q.x) * (p.x - q.x) + (p.y - q.y) * (p.y - q.y)
+    }
+}
+
+/// Marker trait to assign Kernel for scalars
+pub trait HasKernel: CoordinateType {
+    type Ker: Kernel<Scalar = Self>;
+}
+
+// Helper macro to implement `HasKernel` on a a scalar type
+// `T` (first arg.) by assigning the second arg. It expects
+// the second arg. to be a type that takes one generic
+// parameter that is `T`.
+#[macro_use]
+macro_rules! has_kernel {
+	  ($t:ident, $k:ident) => {
+        impl $crate::algorithm::kernels::HasKernel for $t {
+            type Ker = $k<$t>;
+        }
+	  };
+}
+
+pub mod robust;
+pub use self::robust::RobustKernel;
+has_kernel!(f64, RobustKernel);
+has_kernel!(f32, RobustKernel);
+
+pub mod simple;
+pub use self::simple::SimpleKernel;
+
+has_kernel!(i64, SimpleKernel);
+has_kernel!(i32, SimpleKernel);
+has_kernel!(i16, SimpleKernel);
+has_kernel!(isize, SimpleKernel);
+
+#[cfg(has_i128)]
+has_kernel!(i128, SimpleKernel);

geo/src/algorithm/kernels/robust.rs

@@ -0,0 +1,47 @@
+use super::{Kernel, Orientation};
+use crate::Coordinate;
+use std::marker::PhantomData;
+
+/// Robust kernel that uses [fast robust
+/// predicates](//www.cs.cmu.edu/~quake/robust.html) to
+/// provide robust floating point predicates. Should only be
+/// used with types that can _always_ be casted to `f64`
+/// _without loss in precision_.
+#[derive(Default)]
+pub struct RobustKernel<T>(PhantomData<T>);
+
+use num_traits::{Float, NumCast};
+impl<T: Float> Kernel for RobustKernel<T> {
+    type Scalar = T;
+
+    fn orient2d(
+        p: Coordinate<Self::Scalar>,
+        q: Coordinate<Self::Scalar>,
+        r: Coordinate<Self::Scalar>,
+    ) -> Orientation {
+        use robust::{orient2d, Coord};
+
+        let orientation = orient2d(
+            Coord {
+                x: <f64 as NumCast>::from( p.x ).unwrap(),
+                y: <f64 as NumCast>::from( p.y ).unwrap(),
+            },
+            Coord {
+                x: <f64 as NumCast>::from( q.x ).unwrap(),
+                y: <f64 as NumCast>::from( q.y ).unwrap(),
+            },
+            Coord {
+                x: <f64 as NumCast>::from( r.x ).unwrap(),
+                y: <f64 as NumCast>::from( r.y ).unwrap(),
+            },
+        );
+
+        if orientation < 0. {
+            Orientation::Clockwise
+        } else if orientation > 0. {
+            Orientation::CounterClockwise
+        } else {
+            Orientation::Colinear
+        }
+    }
+}

geo/src/algorithm/kernels/simple.rs

@@ -0,0 +1,13 @@
+use super::Kernel;
+use crate::CoordinateType;
+use std::marker::PhantomData;
+
+/// Simple kernel provides the direct implementation of the
+/// predicates. These are meant to be used with exact
+/// arithmetic signed tpyes (eg. i32, i64).
+#[derive(Default)]
+pub struct SimpleKernel<T>(PhantomData<T>);
+
+impl<T: CoordinateType> Kernel for SimpleKernel<T> {
+    type Scalar = T;
+}