From d9c195ad80e90f6d34e1e675c0164f16a7f5d326 Mon Sep 17 00:00:00 2001
From: Rajsekar Manokaran <rajsekar@gmail.com>
Date: Sat, 22 Aug 2020 03:30:46 +0530
Subject: [PATCH 1/8] 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`)
---
 geo-types/src/coordinate.rs               |   2 +-
 geo-types/src/line_string.rs              |  29 +++-
 geo/Cargo.toml                            |   2 +
 geo/src/algorithm/mod.rs                  |   2 +
 geo/src/algorithm/robust_winding_order.rs | 182 ++++++++++++++++++++++
 5 files changed, 212 insertions(+), 5 deletions(-)
 create mode 100644 geo/src/algorithm/robust_winding_order.rs

diff --git a/geo-types/src/coordinate.rs b/geo-types/src/coordinate.rs
index d6f2f7b8f..54755a290 100644
--- a/geo-types/src/coordinate.rs
+++ b/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
diff --git a/geo-types/src/line_string.rs b/geo-types/src/line_string.rs
index 406fb3476..e59c1d7fb 100644
--- a/geo-types/src/line_string.rs
+++ b/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`.
diff --git a/geo/Cargo.toml b/geo/Cargo.toml
index cf26c306e..65026f813 100644
--- a/geo/Cargo.toml
+++ b/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"]
diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs
index 50885e147..7e162d44c 100644
--- a/geo/src/algorithm/mod.rs
+++ b/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`.
diff --git a/geo/src/algorithm/robust_winding_order.rs b/geo/src/algorithm/robust_winding_order.rs
new file mode 100644
index 000000000..40bc8db02
--- /dev/null
+++ b/geo/src/algorithm/robust_winding_order.rs
@@ -0,0 +1,182 @@
+use super::winding_order::*;
+use crate::{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 increment = |x: &mut usize| {
+            *x += 1;
+            if *x >= self.num_coords() { *x = 0; }
+        };
+
+        let decrement = |x: &mut usize| {
+            if *x == 0 {
+                *x = self.num_coords() - 1;
+            } else { *x -= 1; }
+        };
+
+        let i = lexicographically_least_index(&self.0);
+
+        let mut next = i;
+        increment(&mut next);
+        while self.0[next] == self.0[i] {
+            if next == i {
+                // We've looped too much. There aren't
+                // enough unique coords to compute orientation.
+                return None;
+            }
+            increment(&mut next);
+        }
+
+        let mut prev = i;
+        decrement(&mut prev);
+        while self.0[prev] == self.0[i] {
+            if prev == i {
+                // We've looped too much. There aren't enough
+                // unique coords to compute orientation.
+                //
+                // Note: we actually don't need this check as
+                // the previous loop succeeded, and so we have
+                // at least two distinct elements in the list
+                return None;
+            }
+            decrement(&mut prev);
+        }
+
+        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);
+    }
+}

From 262de5ff1e95c47f10763f26f88fd169e91ea122 Mon Sep 17 00:00:00 2001
From: Rajsekar Manokaran <rajsekar@gmail.com>
Date: Sat, 22 Aug 2020 20:59:54 +0530
Subject: [PATCH 2/8] Abstract predicates to Kernel trait

+ add `Kernel` trait (for now only requires `orient2d`)
+ provide `RobustKernel`, and `SimpleKernel` abilities

+ implement `RobustWinding<Kernel>` for `LineString`s
+ add tests
---
 geo-types/src/line_string.rs              |   9 +-
 geo/src/algorithm/kernels/mod.rs          |  18 ++++
 geo/src/algorithm/kernels/robust.rs       |  43 +++++++++
 geo/src/algorithm/kernels/simple.rs       |  27 ++++++
 geo/src/algorithm/mod.rs                  |   3 +
 geo/src/algorithm/robust_winding_order.rs | 106 ++++++----------------
 6 files changed, 122 insertions(+), 84 deletions(-)
 create mode 100644 geo/src/algorithm/kernels/mod.rs
 create mode 100644 geo/src/algorithm/kernels/robust.rs
 create mode 100644 geo/src/algorithm/kernels/simple.rs

diff --git a/geo-types/src/line_string.rs b/geo-types/src/line_string.rs
index e59c1d7fb..17c830482 100644
--- a/geo-types/src/line_string.rs
+++ b/geo-types/src/line_string.rs
@@ -193,14 +193,7 @@ impl<T: CoordinateType> LineString<T> {
     /// 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
-        }
+        self.0.len() > 1 && self.0.first() == self.0.last()
     }
 }
 
diff --git a/geo/src/algorithm/kernels/mod.rs b/geo/src/algorithm/kernels/mod.rs
new file mode 100644
index 000000000..68be23ed4
--- /dev/null
+++ b/geo/src/algorithm/kernels/mod.rs
@@ -0,0 +1,18 @@
+use crate::{CoordinateType, Coordinate};
+use super::winding_order::WindingOrder;
+
+pub trait Kernel {
+    type Scalar: CoordinateType;
+
+    fn orient2d(
+        p: Coordinate<Self::Scalar>,
+        q: Coordinate<Self::Scalar>,
+        r: Coordinate<Self::Scalar>,
+    ) -> Option<WindingOrder>;
+}
+
+pub mod robust;
+pub use self::robust::RobustKernel;
+
+pub mod simple;
+pub use self::simple::SimpleKernel;
diff --git a/geo/src/algorithm/kernels/robust.rs b/geo/src/algorithm/kernels/robust.rs
new file mode 100644
index 000000000..a78c0d84f
--- /dev/null
+++ b/geo/src/algorithm/kernels/robust.rs
@@ -0,0 +1,43 @@
+use super::Kernel;
+use crate::Coordinate;
+use crate::algorithm::winding_order::WindingOrder;
+use std::marker::PhantomData;
+
+#[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>,
+    ) -> Option<WindingOrder> {
+        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. {
+            Some(WindingOrder::Clockwise)
+        } else if orientation > 0. {
+            Some(WindingOrder::CounterClockwise)
+        } else {
+            None
+        }
+    }
+}
diff --git a/geo/src/algorithm/kernels/simple.rs b/geo/src/algorithm/kernels/simple.rs
new file mode 100644
index 000000000..a5461c8c8
--- /dev/null
+++ b/geo/src/algorithm/kernels/simple.rs
@@ -0,0 +1,27 @@
+use super::Kernel;
+use crate::{Coordinate, CoordinateType};
+use crate::algorithm::winding_order::WindingOrder;
+use std::marker::PhantomData;
+
+#[derive(Default)]
+pub struct SimpleKernel<T>(PhantomData<T>);
+
+impl<T: CoordinateType> Kernel for SimpleKernel<T> {
+    type Scalar = T;
+
+    fn orient2d(
+        p: Coordinate<Self::Scalar>,
+        q: Coordinate<Self::Scalar>,
+        r: Coordinate<Self::Scalar>,
+    ) -> Option<WindingOrder> {
+
+        let res = (q.x - p.x) * (r.y - q.y) - (q.y - p.y) * (r.x - q.x);
+        if res > T::zero() {
+            Some(WindingOrder::CounterClockwise)
+        } else if res < T::zero() {
+            Some(WindingOrder::Clockwise)
+        } else {
+            None
+        }
+    }
+}
diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs
index 7e162d44c..6836559b8 100644
--- a/geo/src/algorithm/mod.rs
+++ b/geo/src/algorithm/mod.rs
@@ -1,3 +1,6 @@
+/// Kernels to compute various predicates
+pub mod kernels;
+
 /// Calculate the area of the surface of a `Geometry`.
 pub mod area;
 /// Calculate the bearing to another `Point`, in degrees.
diff --git a/geo/src/algorithm/robust_winding_order.rs b/geo/src/algorithm/robust_winding_order.rs
index 40bc8db02..9fcf04ab5 100644
--- a/geo/src/algorithm/robust_winding_order.rs
+++ b/geo/src/algorithm/robust_winding_order.rs
@@ -1,10 +1,9 @@
 use super::winding_order::*;
+use super::kernels::Kernel;
 use crate::{CoordinateType, LineString};
-use num_traits::Float;
 
-pub trait RobustWinding<T>: Winding<T>
-where
-    T: CoordinateType,
+pub trait RobustWinding<K>
+where K: Kernel
 {
     /// Return the winding order of this object
     fn robust_winding_order(&self) -> Option<WindingOrder>;
@@ -29,7 +28,10 @@ fn lexicographically_least_index<T: Copy + PartialOrd>(pts: &[T]) -> usize {
     min.unwrap().0
 }
 
-impl<T: CoordinateType + Float> RobustWinding<T> for LineString<T> {
+impl<T, K> RobustWinding<K> for LineString<T>
+where T: CoordinateType,
+      K: Kernel<Scalar = 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
@@ -78,30 +80,11 @@ impl<T: CoordinateType + Float> RobustWinding<T> for LineString<T> {
             decrement(&mut prev);
         }
 
-        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
-        }
+        K::orient2d(
+            self.0[prev],
+            self.0[i],
+            self.0[next]
+        )
 
     }
 }
@@ -110,6 +93,7 @@ impl<T: CoordinateType + Float> RobustWinding<T> for LineString<T> {
 mod test {
     use super::*;
     use crate::Point;
+    use crate::algorithm::kernels::*;
 
     #[test]
     fn robust_winding_order() {
@@ -124,59 +108,29 @@ mod test {
         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());
+        let open_line = LineString::from(vec![a.0, b.0, c.0]);
+        assert!(
+            RobustWinding::<RobustKernel<_>>::robust_winding_order(&open_line)
+                .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(),
+            RobustWinding::<RobustKernel<_>>::robust_winding_order(&cw_line),
+            Some(WindingOrder::Clockwise)
+        );
+        assert_eq!(
+            RobustWinding::<RobustKernel<_>>::robust_winding_order(&ccw_line),
             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();
+
+        // This is simple enough that SimpleKernel works
+        assert_eq!(
+            RobustWinding::<SimpleKernel<_>>::robust_winding_order(&cw_line),
+            Some(WindingOrder::Clockwise)
+        );
         assert_eq!(
-            new_line2.robust_winding_order(),
+            RobustWinding::<SimpleKernel<_>>::robust_winding_order(&ccw_line),
             Some(WindingOrder::CounterClockwise)
         );
-        assert_ne!(new_line2, cw_line);
-        assert_eq!(new_line2, ccw_line);
     }
 }

From 7523d683358e17513461baef9fe5c9c823937074 Mon Sep 17 00:00:00 2001
From: Rajsekar Manokaran <rajsekar@gmail.com>
Date: Sat, 22 Aug 2020 21:35:08 +0530
Subject: [PATCH 3/8] Introduce HasKernel

Allows assigning Kernel for each of the basic coordinate
types.
---
 geo/src/algorithm/kernels/mod.rs          | 21 +++++++++++
 geo/src/algorithm/robust_winding_order.rs | 45 ++++++++++++++++-------
 2 files changed, 53 insertions(+), 13 deletions(-)

diff --git a/geo/src/algorithm/kernels/mod.rs b/geo/src/algorithm/kernels/mod.rs
index 68be23ed4..dff9334a0 100644
--- a/geo/src/algorithm/kernels/mod.rs
+++ b/geo/src/algorithm/kernels/mod.rs
@@ -11,8 +11,29 @@ pub trait Kernel {
     ) -> Option<WindingOrder>;
 }
 
+/// Marker trait
+pub trait HasKernel: CoordinateType {
+    type Ker: Kernel<Scalar = Self>;
+}
+
 pub mod robust;
 pub use self::robust::RobustKernel;
 
 pub mod simple;
 pub use self::simple::SimpleKernel;
+
+impl HasKernel for f64 {
+    type Ker = RobustKernel<f64>;
+}
+
+impl HasKernel for f32 {
+    type Ker = RobustKernel<f32>;
+}
+
+impl HasKernel for i64 {
+    type Ker = SimpleKernel<i64>;
+}
+
+impl HasKernel for i32 {
+    type Ker = SimpleKernel<i32>;
+}
diff --git a/geo/src/algorithm/robust_winding_order.rs b/geo/src/algorithm/robust_winding_order.rs
index 9fcf04ab5..fbb836a4e 100644
--- a/geo/src/algorithm/robust_winding_order.rs
+++ b/geo/src/algorithm/robust_winding_order.rs
@@ -1,10 +1,8 @@
 use super::winding_order::*;
-use super::kernels::Kernel;
+use super::kernels::*;
 use crate::{CoordinateType, LineString};
 
-pub trait RobustWinding<K>
-where K: Kernel
-{
+pub trait RobustWinding {
     /// Return the winding order of this object
     fn robust_winding_order(&self) -> Option<WindingOrder>;
 }
@@ -28,8 +26,8 @@ fn lexicographically_least_index<T: Copy + PartialOrd>(pts: &[T]) -> usize {
     min.unwrap().0
 }
 
-impl<T, K> RobustWinding<K> for LineString<T>
-where T: CoordinateType,
+impl<T, K> RobustWinding for LineString<T>
+where T: CoordinateType + HasKernel<Ker = K>,
       K: Kernel<Scalar = T>
 {
     fn robust_winding_order(&self) -> Option<WindingOrder> {
@@ -96,7 +94,7 @@ mod test {
     use crate::algorithm::kernels::*;
 
     #[test]
-    fn robust_winding_order() {
+    fn robust_winding_float() {
         // 3 points forming a triangle
         let a = Point::new(0., 0.);
         let b = Point::new(2., 0.);
@@ -110,27 +108,48 @@ mod test {
         // Verify open linestrings return None
         let open_line = LineString::from(vec![a.0, b.0, c.0]);
         assert!(
-            RobustWinding::<RobustKernel<_>>::robust_winding_order(&open_line)
+            open_line.robust_winding_order()
                 .is_none()
         );
 
         assert_eq!(
-            RobustWinding::<RobustKernel<_>>::robust_winding_order(&cw_line),
+            cw_line.robust_winding_order(),
             Some(WindingOrder::Clockwise)
         );
         assert_eq!(
-            RobustWinding::<RobustKernel<_>>::robust_winding_order(&ccw_line),
+            ccw_line.robust_winding_order(),
             Some(WindingOrder::CounterClockwise)
         );
 
-        // This is simple enough that SimpleKernel works
+    }
+
+    #[test]
+    fn robust_winding_integers() {
+        // 3 points forming a triangle
+        let a = Point::new(0i64, 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
+        let open_line = LineString::from(vec![a.0, b.0, c.0]);
+        assert!(
+            open_line.robust_winding_order()
+                .is_none()
+        );
+
         assert_eq!(
-            RobustWinding::<SimpleKernel<_>>::robust_winding_order(&cw_line),
+            cw_line.robust_winding_order(),
             Some(WindingOrder::Clockwise)
         );
         assert_eq!(
-            RobustWinding::<SimpleKernel<_>>::robust_winding_order(&ccw_line),
+            ccw_line.robust_winding_order(),
             Some(WindingOrder::CounterClockwise)
         );
+
     }
 }

From c65b03a3d09df9a94282d63d2f74bd57d50ed518 Mon Sep 17 00:00:00 2001
From: Rajsekar Manokaran <rajsekar@gmail.com>
Date: Sun, 23 Aug 2020 13:16:35 +0530
Subject: [PATCH 4/8] Add test to verify robustness

Tests to produce png output of output from `orient2d` over a
grid of close-by floats.  These are ignored by default, and
should be explicitly run.

Edits: fix some stable-rust compilation error in CI
---
 geo/Cargo.toml                            |   2 +
 geo/src/algorithm/kernels/mod.rs          |  37 ++---
 geo/src/algorithm/kernels/robust.rs       |   5 +
 geo/src/algorithm/kernels/simple.rs       |   3 +
 geo/src/algorithm/kernels/test.rs         | 167 ++++++++++++++++++++++
 geo/src/algorithm/robust_winding_order.rs | 157 ++++++++++++++------
 geo/src/algorithm/winding_order.rs        |   2 +-
 7 files changed, 314 insertions(+), 59 deletions(-)
 create mode 100644 geo/src/algorithm/kernels/test.rs

diff --git a/geo/Cargo.toml b/geo/Cargo.toml
index 65026f813..cf812bdfc 100644
--- a/geo/Cargo.toml
+++ b/geo/Cargo.toml
@@ -34,6 +34,8 @@ use-serde = ["serde", "geo-types/serde"]
 [dev-dependencies]
 approx = "0.3.0"
 criterion = { version = "0.3" }
+png = "0.16.7"
+float_extras = "0.1.6"
 
 [[bench]]
 name = "area"
diff --git a/geo/src/algorithm/kernels/mod.rs b/geo/src/algorithm/kernels/mod.rs
index dff9334a0..a5b8ed011 100644
--- a/geo/src/algorithm/kernels/mod.rs
+++ b/geo/src/algorithm/kernels/mod.rs
@@ -1,9 +1,13 @@
 use crate::{CoordinateType, Coordinate};
 use super::winding_order::WindingOrder;
 
+/// 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>,
@@ -11,29 +15,30 @@ pub trait Kernel {
     ) -> Option<WindingOrder>;
 }
 
-/// Marker trait
+/// Marker trait to assign Kernel for scalars
 pub trait HasKernel: CoordinateType {
     type Ker: Kernel<Scalar = Self>;
 }
 
+#[macro_export]
+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);
 
-impl HasKernel for f64 {
-    type Ker = RobustKernel<f64>;
-}
-
-impl HasKernel for f32 {
-    type Ker = RobustKernel<f32>;
-}
-
-impl HasKernel for i64 {
-    type Ker = SimpleKernel<i64>;
-}
-
-impl HasKernel for i32 {
-    type Ker = SimpleKernel<i32>;
-}
+#[cfg(test)]
+mod test;
diff --git a/geo/src/algorithm/kernels/robust.rs b/geo/src/algorithm/kernels/robust.rs
index a78c0d84f..2f482b3b5 100644
--- a/geo/src/algorithm/kernels/robust.rs
+++ b/geo/src/algorithm/kernels/robust.rs
@@ -3,6 +3,11 @@ use crate::Coordinate;
 use crate::algorithm::winding_order::WindingOrder;
 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>);
 
diff --git a/geo/src/algorithm/kernels/simple.rs b/geo/src/algorithm/kernels/simple.rs
index a5461c8c8..a505dfcf5 100644
--- a/geo/src/algorithm/kernels/simple.rs
+++ b/geo/src/algorithm/kernels/simple.rs
@@ -3,6 +3,9 @@ use crate::{Coordinate, CoordinateType};
 use crate::algorithm::winding_order::WindingOrder;
 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>);
 
diff --git a/geo/src/algorithm/kernels/test.rs b/geo/src/algorithm/kernels/test.rs
new file mode 100644
index 000000000..673cd71c6
--- /dev/null
+++ b/geo/src/algorithm/kernels/test.rs
@@ -0,0 +1,167 @@
+use super::*;
+use num_traits::*;
+
+// Define a SimpleFloat to test non-robustness of naive
+// implementation.
+#[derive(Copy, Clone, PartialEq, PartialOrd)]
+struct SimpleFloat(pub f64);
+
+use std::ops::*;
+use crate::has_kernel;
+
+macro_rules! impl_ops {
+	  ($t:ident, $m:ident) => {
+		    impl $t for SimpleFloat {
+            type Output = Self;
+            fn $m(self, rhs: Self) -> Self {
+                SimpleFloat((self.0).$m(rhs.0))
+            }
+        }
+	  };
+}
+
+impl_ops!(Rem, rem);
+impl_ops!(Div, div);
+impl_ops!(Mul, mul);
+impl_ops!(Sub, sub);
+impl_ops!(Add, add);
+
+impl One for SimpleFloat {
+    fn one() -> Self {
+        SimpleFloat(One::one())
+    }
+}
+impl Zero for SimpleFloat {
+    fn zero() -> Self {
+        SimpleFloat(Zero::zero())
+    }
+
+    fn is_zero(&self) -> bool {
+        Zero::is_zero(&self.0)
+    }
+}
+
+
+impl Num for SimpleFloat {
+    type FromStrRadixErr = ParseFloatError;
+
+    fn from_str_radix(str: &str, radix: u32) -> Result<Self, Self::FromStrRadixErr> {
+        Num::from_str_radix(str, radix).map(SimpleFloat)
+    }
+}
+
+macro_rules! tp_method {
+	  ($t:ident, $m:ident) => {
+        fn $m(&self) -> Option<$t> {
+            Some(self.0 as $t)
+        }
+	  };
+}
+
+impl ToPrimitive for SimpleFloat{
+    tp_method!(i64, to_i64);
+    tp_method!(u64, to_u64);
+    tp_method!(f64, to_f64);
+}
+
+impl NumCast for SimpleFloat {
+    fn from<T: ToPrimitive>(n: T) -> Option<Self> {
+        NumCast::from(n).map(SimpleFloat)
+    }
+}
+
+impl From<f64> for SimpleFloat {
+    fn from(f: f64) -> Self {
+        SimpleFloat(f)
+    }
+}
+
+has_kernel!(SimpleFloat, SimpleKernel);
+
+fn orient2d_tests<T: From<f64> + CoordinateType + HasKernel>(
+    x1: f64, y1: f64,
+    x2: f64, y2: f64,
+    x3: f64, y3: f64,
+    width: usize, height: usize,
+) -> Vec<Option<WindingOrder>> {
+    let p1 = Coordinate{
+        x: <T as From<_>>::from(x1),
+        y: <T as From<_>>::from(y1),
+    };
+    let p3 = Coordinate{
+        x: <T as From<_>>::from(x3),
+        y: <T as From<_>>::from(y3),
+    };
+
+    use float_extras::f64::nextafter;
+    let mut yd2 = y2;
+    let mut data = Vec::with_capacity(width * height);
+
+    for _ in 0..height {
+        let mut xd2 = x2;
+        for _ in 0..width {
+            let p2 = Coordinate{
+                x: <T as From<_>>::from(xd2),
+                y: <T as From<_>>::from(yd2),
+            };
+            xd2 = nextafter(xd2, std::f64::INFINITY);
+            data.push(<T as HasKernel>::Ker::orient2d(p1, p2, p3));
+        }
+        yd2 = nextafter(yd2, std::f64::INFINITY);
+    }
+
+    data
+}
+
+use std::path::Path;
+fn write_png(
+    data: &[Option<WindingOrder>],
+    path: &Path,
+    width: usize, height: usize,
+) {
+    assert_eq!(data.len(), width * height);
+
+    use std::fs::File;
+    use std::io::BufWriter;
+
+    let file = File::create(path).unwrap();
+    let ref mut w = BufWriter::new(file);
+
+    let mut encoder = png::Encoder::new(w, width as u32, height as u32);
+    encoder.set_color(png::ColorType::Grayscale);
+    encoder.set_depth(png::BitDepth::Eight);
+
+    let mut writer = encoder.write_header().unwrap();
+    let data = data.iter().map(|w| {
+        match w {
+            Some(WindingOrder::Clockwise) => 0u8,
+            None => 127,
+            Some(WindingOrder::CounterClockwise) => 255,
+        }
+    }).collect::<Vec<_>>();
+    writer.write_image_data(&data).unwrap();
+}
+
+#[test]
+#[ignore]
+fn test_naive() {
+    let data = orient2d_tests::<SimpleFloat>(
+        12.0, 12.0,
+        0.5, 0.5,
+        24.0, 24.0,
+        256, 256
+    );
+    write_png(&data, Path::new("naive-orientation-map.png"), 256, 256);
+}
+
+#[test]
+#[ignore]
+fn test_robust() {
+    let data = orient2d_tests::<f64>(
+        12.0, 12.0,
+        0.5, 0.5,
+        24.0, 24.0,
+        256, 256
+    );
+    write_png(&data, Path::new("robust-orientation-map.png"), 256, 256);
+}
diff --git a/geo/src/algorithm/robust_winding_order.rs b/geo/src/algorithm/robust_winding_order.rs
index fbb836a4e..448a42724 100644
--- a/geo/src/algorithm/robust_winding_order.rs
+++ b/geo/src/algorithm/robust_winding_order.rs
@@ -1,10 +1,59 @@
-use super::winding_order::*;
+use super::winding_order::{WindingOrder, Points};
 use super::kernels::*;
 use crate::{CoordinateType, LineString};
+use crate::utils::EitherIter;
 
 pub trait RobustWinding {
+    type Scalar: CoordinateType;
+
     /// Return the winding order of this object
-    fn robust_winding_order(&self) -> Option<WindingOrder>;
+    fn winding_order(&self) -> Option<WindingOrder>;
+
+    /// True iff this clockwise
+    fn is_cw(&self) -> bool {
+        self.winding_order() == Some(WindingOrder::Clockwise)
+    }
+
+    /// True iff this is wound counterclockwise
+    fn is_ccw(&self) -> bool {
+        self.winding_order() == Some(WindingOrder::CounterClockwise)
+    }
+
+    /// Iterate over the points in a clockwise order
+    ///
+    /// The object isn't changed, and the points are returned either in order, or in reverse
+    /// order, so that the resultant order makes it appear clockwise
+    fn points_cw(&self) -> Points<Self::Scalar>;
+
+    /// Iterate over the points in a counter-clockwise order
+    ///
+    /// The object isn't changed, and the points are returned either in order, or in reverse
+    /// order, so that the resultant order makes it appear counter-clockwise
+    fn points_ccw(&self) -> Points<Self::Scalar>;
+
+    /// Change this objects's points so they are in clockwise winding order
+    fn make_cw_winding(&mut self);
+
+    /// Change this line's points so they are in counterclockwise winding order
+    fn make_ccw_winding(&mut self);
+
+    /// Return a clone of this object, but in the specified winding order
+    fn clone_to_winding_order(&self, winding_order: WindingOrder) -> Self
+    where
+        Self: Sized + Clone,
+    {
+        let mut new: Self = self.clone();
+        new.make_winding_order(winding_order);
+        new
+    }
+
+    /// Change the winding order so that it is in this winding order
+    fn make_winding_order(&mut self, winding_order: WindingOrder) {
+        match winding_order {
+            WindingOrder::Clockwise => self.make_cw_winding(),
+            WindingOrder::CounterClockwise => self.make_ccw_winding(),
+        }
+    }
 }
 
 /// Compute index of the lexicographically least point.
@@ -30,7 +79,9 @@ impl<T, K> RobustWinding for LineString<T>
 where T: CoordinateType + HasKernel<Ker = K>,
       K: Kernel<Scalar = T>
 {
-    fn robust_winding_order(&self) -> Option<WindingOrder> {
+    type Scalar = T;
+
+    fn 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.
@@ -66,15 +117,9 @@ where T: CoordinateType + HasKernel<Ker = K>,
         let mut prev = i;
         decrement(&mut prev);
         while self.0[prev] == self.0[i] {
-            if prev == i {
-                // We've looped too much. There aren't enough
-                // unique coords to compute orientation.
-                //
-                // Note: we actually don't need this check as
-                // the previous loop succeeded, and so we have
-                // at least two distinct elements in the list
-                return None;
-            }
+            // Note: we don't need to check if prev == i as
+            // the previous loop succeeded, and so we have
+            // at least two distinct elements in the list
             decrement(&mut prev);
         }
 
@@ -85,13 +130,48 @@ where T: CoordinateType + HasKernel<Ker = K>,
         )
 
     }
+
+    /// Iterate over the points in a clockwise order
+    ///
+    /// The Linestring isn't changed, and the points are returned either in order, or in reverse
+    /// order, so that the resultant order makes it appear clockwise
+    fn points_cw(&self) -> Points<Self::Scalar> {
+        match self.winding_order() {
+            Some(WindingOrder::CounterClockwise) => Points(EitherIter::B(self.points_iter().rev())),
+            _ => Points(EitherIter::A(self.points_iter())),
+        }
+    }
+
+    /// Iterate over the points in a counter-clockwise order
+    ///
+    /// The Linestring isn't changed, and the points are returned either in order, or in reverse
+    /// order, so that the resultant order makes it appear counter-clockwise
+    fn points_ccw(&self) -> Points<Self::Scalar> {
+        match self.winding_order() {
+            Some(WindingOrder::Clockwise) => Points(EitherIter::B(self.points_iter().rev())),
+            _ => Points(EitherIter::A(self.points_iter())),
+        }
+    }
+
+    /// Change this line's points so they are in clockwise winding order
+    fn make_cw_winding(&mut self) {
+        if let Some(WindingOrder::CounterClockwise) = self.winding_order() {
+            self.0.reverse();
+        }
+    }
+
+    /// Change this line's points so they are in counterclockwise winding order
+    fn make_ccw_winding(&mut self) {
+        if let Some(WindingOrder::Clockwise) = self.winding_order() {
+            self.0.reverse();
+        }
+    }
 }
 
 #[cfg(test)]
 mod test {
     use super::*;
     use crate::Point;
-    use crate::algorithm::kernels::*;
 
     #[test]
     fn robust_winding_float() {
@@ -100,56 +180,49 @@ mod test {
         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
-        let open_line = LineString::from(vec![a.0, b.0, c.0]);
+        let mut ls = LineString::from(vec![a.0, b.0, c.0]);
         assert!(
-            open_line.robust_winding_order()
+            ls.winding_order()
                 .is_none()
         );
 
+        ls.0.push(ls.0[0]);
         assert_eq!(
-            cw_line.robust_winding_order(),
-            Some(WindingOrder::Clockwise)
+            ls.winding_order(),
+            Some(WindingOrder::CounterClockwise)
         );
+
+        ls.make_cw_winding();
         assert_eq!(
-            ccw_line.robust_winding_order(),
-            Some(WindingOrder::CounterClockwise)
+            ls.winding_order(),
+            Some(WindingOrder::Clockwise)
         );
 
     }
 
     #[test]
-    fn robust_winding_integers() {
+    fn robust_winding_integer() {
         // 3 points forming a triangle
         let a = Point::new(0i64, 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
-        let open_line = LineString::from(vec![a.0, b.0, c.0]);
-        assert!(
-            open_line.robust_winding_order()
-                .is_none()
-        );
+        let mut ls = LineString::from(vec![a.0, b.0, c.0]);
+        assert!(ls.winding_order().is_none());
+
+        ls.0.push(ls.0[0]);
+        assert!(ls.is_ccw());
+
+        let ccw_ls: Vec<_> = ls.points_ccw().collect();
+
+        ls.make_cw_winding();
+        assert!(ls.is_cw());
 
         assert_eq!(
-            cw_line.robust_winding_order(),
-            Some(WindingOrder::Clockwise)
-        );
-        assert_eq!(
-            ccw_line.robust_winding_order(),
-            Some(WindingOrder::CounterClockwise)
+            &ls.points_ccw().collect::<Vec<_>>(),
+            &ccw_ls,
         );
-
     }
 }
diff --git a/geo/src/algorithm/winding_order.rs b/geo/src/algorithm/winding_order.rs
index 431d47c9b..1ec7804bc 100644
--- a/geo/src/algorithm/winding_order.rs
+++ b/geo/src/algorithm/winding_order.rs
@@ -42,7 +42,7 @@ where
 }
 
 /// Iterates through a list of `Point`s
-pub struct Points<'a, T>(EitherIter<Point<T>, PointsIter<'a, T>, Rev<PointsIter<'a, T>>>)
+pub struct Points<'a, T>(pub(crate) EitherIter<Point<T>, PointsIter<'a, T>, Rev<PointsIter<'a, T>>>)
 where
     T: CoordinateType + 'a;
 

From 9211c5a33016cd05d0ac2f31a454a8b45d897fd9 Mon Sep 17 00:00:00 2001
From: Rajsekar Manokaran <rajsekar@gmail.com>
Date: Sun, 23 Aug 2020 21:47:09 +0530
Subject: [PATCH 5/8] Add Graham Convex Hull algorithm

Unlike quick hull, this doesn't need `T: Float`, and thus
works for both integer and float scalars.

+ move `lexicographically_least_index` to utils
+ make `LineString::close` public
+ make separate enum `Orientation` for `Kernel::orient2d`
+ add quick hull tests to graham hull
---
 geo-types/src/line_string.rs              |   2 +-
 geo/src/algorithm/graham_hull.rs          | 160 ++++++++++++++++++++++
 geo/src/algorithm/kernels/mod.rs          |  29 +++-
 geo/src/algorithm/kernels/robust.rs       |  11 +-
 geo/src/algorithm/kernels/simple.rs       |  19 +--
 geo/src/algorithm/kernels/test.rs         |  11 +-
 geo/src/algorithm/mod.rs                  |   1 +
 geo/src/algorithm/robust_winding_order.rs |  29 ++--
 geo/src/utils.rs                          |  19 +++
 9 files changed, 228 insertions(+), 53 deletions(-)
 create mode 100644 geo/src/algorithm/graham_hull.rs

diff --git a/geo-types/src/line_string.rs b/geo-types/src/line_string.rs
index 17c830482..80702795c 100644
--- a/geo-types/src/line_string.rs
+++ b/geo-types/src/line_string.rs
@@ -159,7 +159,7 @@ 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) {
+    pub fn close(&mut self) {
         if !self.is_closed() && !self.0.is_empty() {
             self.0.push(self.0[0]);
         }
diff --git a/geo/src/algorithm/graham_hull.rs b/geo/src/algorithm/graham_hull.rs
new file mode 100644
index 000000000..c61b964a6
--- /dev/null
+++ b/geo/src/algorithm/graham_hull.rs
@@ -0,0 +1,160 @@
+use crate::{Coordinate, CoordinateType, 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::robust_winding_order::RobustWinding;
+        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
+}
+
+pub trait ConvexHull {
+    type Scalar: CoordinateType;
+    fn convex_hull(&self) -> LineString<Self::Scalar>;
+}
+
+
+#[cfg(test)]
+mod test {
+    use super::*;
+
+    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);
+    }
+}
diff --git a/geo/src/algorithm/kernels/mod.rs b/geo/src/algorithm/kernels/mod.rs
index a5b8ed011..4b4ec722e 100644
--- a/geo/src/algorithm/kernels/mod.rs
+++ b/geo/src/algorithm/kernels/mod.rs
@@ -1,5 +1,11 @@
 use crate::{CoordinateType, Coordinate};
-use super::winding_order::WindingOrder;
+
+#[derive(Debug, PartialEq, Eq, PartialOrd, Ord)]
+pub enum Orientation {
+    CounterClockwise,
+    Clockwise,
+    Colinear,
+}
 
 /// Kernel trait to provide predicates to operate on
 /// different scalar types.
@@ -12,7 +18,26 @@ pub trait Kernel {
         p: Coordinate<Self::Scalar>,
         q: Coordinate<Self::Scalar>,
         r: Coordinate<Self::Scalar>,
-    ) -> Option<WindingOrder>;
+    ) -> 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
diff --git a/geo/src/algorithm/kernels/robust.rs b/geo/src/algorithm/kernels/robust.rs
index 2f482b3b5..4a580a1bd 100644
--- a/geo/src/algorithm/kernels/robust.rs
+++ b/geo/src/algorithm/kernels/robust.rs
@@ -1,6 +1,5 @@
-use super::Kernel;
+use super::{Kernel, Orientation};
 use crate::Coordinate;
-use crate::algorithm::winding_order::WindingOrder;
 use std::marker::PhantomData;
 
 /// Robust kernel that uses [fast robust
@@ -19,7 +18,7 @@ impl<T: Float> Kernel for RobustKernel<T> {
         p: Coordinate<Self::Scalar>,
         q: Coordinate<Self::Scalar>,
         r: Coordinate<Self::Scalar>,
-    ) -> Option<WindingOrder> {
+    ) -> Orientation {
         use robust::{orient2d, Coord};
 
         let orientation = orient2d(
@@ -38,11 +37,11 @@ impl<T: Float> Kernel for RobustKernel<T> {
         );
 
         if orientation < 0. {
-            Some(WindingOrder::Clockwise)
+            Orientation::Clockwise
         } else if orientation > 0. {
-            Some(WindingOrder::CounterClockwise)
+            Orientation::CounterClockwise
         } else {
-            None
+            Orientation::Colinear
         }
     }
 }
diff --git a/geo/src/algorithm/kernels/simple.rs b/geo/src/algorithm/kernels/simple.rs
index a505dfcf5..a50c08a59 100644
--- a/geo/src/algorithm/kernels/simple.rs
+++ b/geo/src/algorithm/kernels/simple.rs
@@ -1,6 +1,5 @@
 use super::Kernel;
-use crate::{Coordinate, CoordinateType};
-use crate::algorithm::winding_order::WindingOrder;
+use crate::CoordinateType;
 use std::marker::PhantomData;
 
 /// Simple kernel provides the direct implementation of the
@@ -11,20 +10,4 @@ pub struct SimpleKernel<T>(PhantomData<T>);
 
 impl<T: CoordinateType> Kernel for SimpleKernel<T> {
     type Scalar = T;
-
-    fn orient2d(
-        p: Coordinate<Self::Scalar>,
-        q: Coordinate<Self::Scalar>,
-        r: Coordinate<Self::Scalar>,
-    ) -> Option<WindingOrder> {
-
-        let res = (q.x - p.x) * (r.y - q.y) - (q.y - p.y) * (r.x - q.x);
-        if res > T::zero() {
-            Some(WindingOrder::CounterClockwise)
-        } else if res < T::zero() {
-            Some(WindingOrder::Clockwise)
-        } else {
-            None
-        }
-    }
 }
diff --git a/geo/src/algorithm/kernels/test.rs b/geo/src/algorithm/kernels/test.rs
index 673cd71c6..a78388103 100644
--- a/geo/src/algorithm/kernels/test.rs
+++ b/geo/src/algorithm/kernels/test.rs
@@ -83,7 +83,7 @@ fn orient2d_tests<T: From<f64> + CoordinateType + HasKernel>(
     x2: f64, y2: f64,
     x3: f64, y3: f64,
     width: usize, height: usize,
-) -> Vec<Option<WindingOrder>> {
+) -> Vec<Orientation> {
     let p1 = Coordinate{
         x: <T as From<_>>::from(x1),
         y: <T as From<_>>::from(y1),
@@ -115,7 +115,7 @@ fn orient2d_tests<T: From<f64> + CoordinateType + HasKernel>(
 
 use std::path::Path;
 fn write_png(
-    data: &[Option<WindingOrder>],
+    data: &[Orientation],
     path: &Path,
     width: usize, height: usize,
 ) {
@@ -133,10 +133,11 @@ fn write_png(
 
     let mut writer = encoder.write_header().unwrap();
     let data = data.iter().map(|w| {
+        use Orientation::*;
         match w {
-            Some(WindingOrder::Clockwise) => 0u8,
-            None => 127,
-            Some(WindingOrder::CounterClockwise) => 255,
+            Clockwise => 0u8,
+            Colinear => 127,
+            CounterClockwise => 255,
         }
     }).collect::<Vec<_>>();
     writer.write_image_data(&data).unwrap();
diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs
index 6836559b8..817ae6fe2 100644
--- a/geo/src/algorithm/mod.rs
+++ b/geo/src/algorithm/mod.rs
@@ -17,6 +17,7 @@ pub mod closest_point;
 pub mod contains;
 /// Calculate the convex hull of a `Geometry`.
 pub mod convexhull;
+pub mod graham_hull;
 /// Calculate the minimum Euclidean distance between two `Geometries`.
 pub mod euclidean_distance;
 /// Calculate the length of a planar line between two `Geometries`.
diff --git a/geo/src/algorithm/robust_winding_order.rs b/geo/src/algorithm/robust_winding_order.rs
index 448a42724..229daab77 100644
--- a/geo/src/algorithm/robust_winding_order.rs
+++ b/geo/src/algorithm/robust_winding_order.rs
@@ -1,4 +1,4 @@
-use super::winding_order::{WindingOrder, Points};
+use super::winding_order::{Points, WindingOrder};
 use super::kernels::*;
 use crate::{CoordinateType, LineString};
 use crate::utils::EitherIter;
@@ -56,24 +56,6 @@ pub trait RobustWinding {
     }
 }
 
-/// 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, K> RobustWinding for LineString<T>
 where T: CoordinateType + HasKernel<Ker = K>,
@@ -101,6 +83,7 @@ where T: CoordinateType + HasKernel<Ker = K>,
             } else { *x -= 1; }
         };
 
+        use crate::utils::lexicographically_least_index;
         let i = lexicographically_least_index(&self.0);
 
         let mut next = i;
@@ -123,11 +106,15 @@ where T: CoordinateType + HasKernel<Ker = K>,
             decrement(&mut prev);
         }
 
-        K::orient2d(
+        match K::orient2d(
             self.0[prev],
             self.0[i],
             self.0[next]
-        )
+        ) {
+            Orientation::CounterClockwise => Some(WindingOrder::CounterClockwise),
+            Orientation::Clockwise => Some(WindingOrder::Clockwise),
+            _ => None
+        }
 
     }
 
diff --git a/geo/src/utils.rs b/geo/src/utils.rs
index 1da1c5636..b2bb2d170 100644
--- a/geo/src/utils.rs
+++ b/geo/src/utils.rs
@@ -127,6 +127,25 @@ where
     }
 }
 
+/// Compute index of the lexicographically least point.
+/// Should only be called on a non-empty slice.
+pub 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
+}
+
 #[cfg(test)]
 mod test {
     use super::{partial_max, partial_min};

From 27b9af9074a789dd52b54fdce80100e189a98e05 Mon Sep 17 00:00:00 2001
From: Rajsekar Manokaran <rajsekar@gmail.com>
Date: Mon, 24 Aug 2020 02:53:17 +0530
Subject: [PATCH 6/8] Minor code move around

+ make `kernels` module `pub(crate)`
+ remove `kernels` test in favor of adding it in the
  `robust` crate
+ replace existing `Winding` trait with robust variant
+ fix `Orient` trait to use new `Winding` trait
---
 geo-types/src/line_string.rs              |   5 +-
 geo/src/algorithm/area.rs                 |  38 +++-
 geo/src/algorithm/graham_hull.rs          |   2 +-
 geo/src/algorithm/kernels/mod.rs          |   3 -
 geo/src/algorithm/kernels/test.rs         | 168 ----------------
 geo/src/algorithm/mod.rs                  |   4 +-
 geo/src/algorithm/orient.rs               |  13 +-
 geo/src/algorithm/robust_winding_order.rs | 215 ---------------------
 geo/src/algorithm/winding_order.rs        | 222 ++++++++++------------
 9 files changed, 154 insertions(+), 516 deletions(-)
 delete mode 100644 geo/src/algorithm/kernels/test.rs
 delete mode 100644 geo/src/algorithm/robust_winding_order.rs

diff --git a/geo-types/src/line_string.rs b/geo-types/src/line_string.rs
index 80702795c..4e967e3c5 100644
--- a/geo-types/src/line_string.rs
+++ b/geo-types/src/line_string.rs
@@ -181,7 +181,8 @@ impl<T: CoordinateType> LineString<T> {
     }
 
     /// Checks if the linestring is closed; i.e. the first
-    /// and last points have the same coords.
+    /// and last points have the same coords. Note that a
+    /// single point is considered closed.
     ///
     /// # Examples
     ///
@@ -193,7 +194,7 @@ impl<T: CoordinateType> LineString<T> {
     /// assert!(line_string.is_closed());
     /// ```
     pub fn is_closed(&self) -> bool {
-        self.0.len() > 1 && self.0.first() == self.0.last()
+        !self.0.is_empty() && self.0.first() == self.0.last()
     }
 }
 
diff --git a/geo/src/algorithm/area.rs b/geo/src/algorithm/area.rs
index 43fa0a7d4..502bc5183 100644
--- a/geo/src/algorithm/area.rs
+++ b/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.
 ///
diff --git a/geo/src/algorithm/graham_hull.rs b/geo/src/algorithm/graham_hull.rs
index c61b964a6..0ff694a54 100644
--- a/geo/src/algorithm/graham_hull.rs
+++ b/geo/src/algorithm/graham_hull.rs
@@ -20,7 +20,7 @@ where
         // Nothing to build with fewer than four points. We
         // remove repeated points if any, and ensure ccw
         // invariant.
-        use super::robust_winding_order::RobustWinding;
+        use super::winding_order::Winding;
         let mut ls: LineString<T> = points
             .iter()
             .enumerate()
diff --git a/geo/src/algorithm/kernels/mod.rs b/geo/src/algorithm/kernels/mod.rs
index 4b4ec722e..418c2d88a 100644
--- a/geo/src/algorithm/kernels/mod.rs
+++ b/geo/src/algorithm/kernels/mod.rs
@@ -64,6 +64,3 @@ pub mod simple;
 pub use self::simple::SimpleKernel;
 has_kernel!(i64, SimpleKernel);
 has_kernel!(i32, SimpleKernel);
-
-#[cfg(test)]
-mod test;
diff --git a/geo/src/algorithm/kernels/test.rs b/geo/src/algorithm/kernels/test.rs
deleted file mode 100644
index a78388103..000000000
--- a/geo/src/algorithm/kernels/test.rs
+++ /dev/null
@@ -1,168 +0,0 @@
-use super::*;
-use num_traits::*;
-
-// Define a SimpleFloat to test non-robustness of naive
-// implementation.
-#[derive(Copy, Clone, PartialEq, PartialOrd)]
-struct SimpleFloat(pub f64);
-
-use std::ops::*;
-use crate::has_kernel;
-
-macro_rules! impl_ops {
-	  ($t:ident, $m:ident) => {
-		    impl $t for SimpleFloat {
-            type Output = Self;
-            fn $m(self, rhs: Self) -> Self {
-                SimpleFloat((self.0).$m(rhs.0))
-            }
-        }
-	  };
-}
-
-impl_ops!(Rem, rem);
-impl_ops!(Div, div);
-impl_ops!(Mul, mul);
-impl_ops!(Sub, sub);
-impl_ops!(Add, add);
-
-impl One for SimpleFloat {
-    fn one() -> Self {
-        SimpleFloat(One::one())
-    }
-}
-impl Zero for SimpleFloat {
-    fn zero() -> Self {
-        SimpleFloat(Zero::zero())
-    }
-
-    fn is_zero(&self) -> bool {
-        Zero::is_zero(&self.0)
-    }
-}
-
-
-impl Num for SimpleFloat {
-    type FromStrRadixErr = ParseFloatError;
-
-    fn from_str_radix(str: &str, radix: u32) -> Result<Self, Self::FromStrRadixErr> {
-        Num::from_str_radix(str, radix).map(SimpleFloat)
-    }
-}
-
-macro_rules! tp_method {
-	  ($t:ident, $m:ident) => {
-        fn $m(&self) -> Option<$t> {
-            Some(self.0 as $t)
-        }
-	  };
-}
-
-impl ToPrimitive for SimpleFloat{
-    tp_method!(i64, to_i64);
-    tp_method!(u64, to_u64);
-    tp_method!(f64, to_f64);
-}
-
-impl NumCast for SimpleFloat {
-    fn from<T: ToPrimitive>(n: T) -> Option<Self> {
-        NumCast::from(n).map(SimpleFloat)
-    }
-}
-
-impl From<f64> for SimpleFloat {
-    fn from(f: f64) -> Self {
-        SimpleFloat(f)
-    }
-}
-
-has_kernel!(SimpleFloat, SimpleKernel);
-
-fn orient2d_tests<T: From<f64> + CoordinateType + HasKernel>(
-    x1: f64, y1: f64,
-    x2: f64, y2: f64,
-    x3: f64, y3: f64,
-    width: usize, height: usize,
-) -> Vec<Orientation> {
-    let p1 = Coordinate{
-        x: <T as From<_>>::from(x1),
-        y: <T as From<_>>::from(y1),
-    };
-    let p3 = Coordinate{
-        x: <T as From<_>>::from(x3),
-        y: <T as From<_>>::from(y3),
-    };
-
-    use float_extras::f64::nextafter;
-    let mut yd2 = y2;
-    let mut data = Vec::with_capacity(width * height);
-
-    for _ in 0..height {
-        let mut xd2 = x2;
-        for _ in 0..width {
-            let p2 = Coordinate{
-                x: <T as From<_>>::from(xd2),
-                y: <T as From<_>>::from(yd2),
-            };
-            xd2 = nextafter(xd2, std::f64::INFINITY);
-            data.push(<T as HasKernel>::Ker::orient2d(p1, p2, p3));
-        }
-        yd2 = nextafter(yd2, std::f64::INFINITY);
-    }
-
-    data
-}
-
-use std::path::Path;
-fn write_png(
-    data: &[Orientation],
-    path: &Path,
-    width: usize, height: usize,
-) {
-    assert_eq!(data.len(), width * height);
-
-    use std::fs::File;
-    use std::io::BufWriter;
-
-    let file = File::create(path).unwrap();
-    let ref mut w = BufWriter::new(file);
-
-    let mut encoder = png::Encoder::new(w, width as u32, height as u32);
-    encoder.set_color(png::ColorType::Grayscale);
-    encoder.set_depth(png::BitDepth::Eight);
-
-    let mut writer = encoder.write_header().unwrap();
-    let data = data.iter().map(|w| {
-        use Orientation::*;
-        match w {
-            Clockwise => 0u8,
-            Colinear => 127,
-            CounterClockwise => 255,
-        }
-    }).collect::<Vec<_>>();
-    writer.write_image_data(&data).unwrap();
-}
-
-#[test]
-#[ignore]
-fn test_naive() {
-    let data = orient2d_tests::<SimpleFloat>(
-        12.0, 12.0,
-        0.5, 0.5,
-        24.0, 24.0,
-        256, 256
-    );
-    write_png(&data, Path::new("naive-orientation-map.png"), 256, 256);
-}
-
-#[test]
-#[ignore]
-fn test_robust() {
-    let data = orient2d_tests::<f64>(
-        12.0, 12.0,
-        0.5, 0.5,
-        24.0, 24.0,
-        256, 256
-    );
-    write_png(&data, Path::new("robust-orientation-map.png"), 256, 256);
-}
diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs
index 817ae6fe2..ad5d72b35 100644
--- a/geo/src/algorithm/mod.rs
+++ b/geo/src/algorithm/mod.rs
@@ -1,5 +1,5 @@
 /// Kernels to compute various predicates
-pub mod kernels;
+pub(crate) mod kernels;
 
 /// Calculate the area of the surface of a `Geometry`.
 pub mod area;
@@ -63,8 +63,6 @@ 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`.
diff --git a/geo/src/algorithm/orient.rs b/geo/src/algorithm/orient.rs
index 9c5c04169..41cd5226b 100644
--- a/geo/src/algorithm/orient.rs
+++ b/geo/src/algorithm/orient.rs
@@ -1,8 +1,9 @@
 use crate::{CoordinateType, MultiPolygon, Polygon};
+use super::kernels::*;
 
 use crate::algorithm::winding_order::{Winding, WindingOrder};
 
-pub trait Orient<T> {
+pub trait Orient {
     /// Orients a Polygon's exterior and interior rings according to convention
     ///
     /// By default, the exterior ring of a Polygon is oriented counter-clockwise, and any interior
@@ -65,18 +66,18 @@ pub trait Orient<T> {
     fn orient(&self, orientation: Direction) -> Self;
 }
 
-impl<T> Orient<T> for Polygon<T>
+impl<T> Orient for Polygon<T>
 where
-    T: CoordinateType,
+    T: CoordinateType + HasKernel ,
 {
     fn orient(&self, direction: Direction) -> Polygon<T> {
         orient(self, direction)
     }
 }
 
-impl<T> Orient<T> for MultiPolygon<T>
+impl<T> Orient for MultiPolygon<T>
 where
-    T: CoordinateType,
+    T: CoordinateType + HasKernel,
 {
     fn orient(&self, direction: Direction) -> MultiPolygon<T> {
         MultiPolygon(self.0.iter().map(|poly| poly.orient(direction)).collect())
@@ -99,7 +100,7 @@ pub enum Direction {
 // and the interior ring(s) will be oriented clockwise
 fn orient<T>(poly: &Polygon<T>, direction: Direction) -> Polygon<T>
 where
-    T: CoordinateType,
+    T: CoordinateType + HasKernel,
 {
     let interiors = poly
         .interiors()
diff --git a/geo/src/algorithm/robust_winding_order.rs b/geo/src/algorithm/robust_winding_order.rs
deleted file mode 100644
index 229daab77..000000000
--- a/geo/src/algorithm/robust_winding_order.rs
+++ /dev/null
@@ -1,215 +0,0 @@
-use super::winding_order::{Points, WindingOrder};
-use super::kernels::*;
-use crate::{CoordinateType, LineString};
-use crate::utils::EitherIter;
-
-pub trait RobustWinding {
-    type Scalar: CoordinateType;
-
-    /// Return the winding order of this object
-    fn winding_order(&self) -> Option<WindingOrder>;
-
-    /// True iff this clockwise
-    fn is_cw(&self) -> bool {
-        self.winding_order() == Some(WindingOrder::Clockwise)
-    }
-
-    /// True iff this is wound counterclockwise
-    fn is_ccw(&self) -> bool {
-        self.winding_order() == Some(WindingOrder::CounterClockwise)
-    }
-
-    /// Iterate over the points in a clockwise order
-    ///
-    /// The object isn't changed, and the points are returned either in order, or in reverse
-    /// order, so that the resultant order makes it appear clockwise
-    fn points_cw(&self) -> Points<Self::Scalar>;
-
-    /// Iterate over the points in a counter-clockwise order
-    ///
-    /// The object isn't changed, and the points are returned either in order, or in reverse
-    /// order, so that the resultant order makes it appear counter-clockwise
-    fn points_ccw(&self) -> Points<Self::Scalar>;
-
-    /// Change this objects's points so they are in clockwise winding order
-    fn make_cw_winding(&mut self);
-
-    /// Change this line's points so they are in counterclockwise winding order
-    fn make_ccw_winding(&mut self);
-
-    /// Return a clone of this object, but in the specified winding order
-    fn clone_to_winding_order(&self, winding_order: WindingOrder) -> Self
-    where
-        Self: Sized + Clone,
-    {
-        let mut new: Self = self.clone();
-        new.make_winding_order(winding_order);
-        new
-    }
-
-    /// Change the winding order so that it is in this winding order
-    fn make_winding_order(&mut self, winding_order: WindingOrder) {
-        match winding_order {
-            WindingOrder::Clockwise => self.make_cw_winding(),
-            WindingOrder::CounterClockwise => self.make_ccw_winding(),
-        }
-    }
-}
-
-
-impl<T, K> RobustWinding for LineString<T>
-where T: CoordinateType + HasKernel<Ker = K>,
-      K: Kernel<Scalar = T>
-{
-    type Scalar = T;
-
-    fn 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 increment = |x: &mut usize| {
-            *x += 1;
-            if *x >= self.num_coords() { *x = 0; }
-        };
-
-        let decrement = |x: &mut usize| {
-            if *x == 0 {
-                *x = self.num_coords() - 1;
-            } else { *x -= 1; }
-        };
-
-        use crate::utils::lexicographically_least_index;
-        let i = lexicographically_least_index(&self.0);
-
-        let mut next = i;
-        increment(&mut next);
-        while self.0[next] == self.0[i] {
-            if next == i {
-                // We've looped too much. There aren't
-                // enough unique coords to compute orientation.
-                return None;
-            }
-            increment(&mut next);
-        }
-
-        let mut prev = i;
-        decrement(&mut prev);
-        while self.0[prev] == self.0[i] {
-            // Note: we don't need to check if prev == i as
-            // the previous loop succeeded, and so we have
-            // at least two distinct elements in the list
-            decrement(&mut prev);
-        }
-
-        match K::orient2d(
-            self.0[prev],
-            self.0[i],
-            self.0[next]
-        ) {
-            Orientation::CounterClockwise => Some(WindingOrder::CounterClockwise),
-            Orientation::Clockwise => Some(WindingOrder::Clockwise),
-            _ => None
-        }
-
-    }
-
-    /// Iterate over the points in a clockwise order
-    ///
-    /// The Linestring isn't changed, and the points are returned either in order, or in reverse
-    /// order, so that the resultant order makes it appear clockwise
-    fn points_cw(&self) -> Points<Self::Scalar> {
-        match self.winding_order() {
-            Some(WindingOrder::CounterClockwise) => Points(EitherIter::B(self.points_iter().rev())),
-            _ => Points(EitherIter::A(self.points_iter())),
-        }
-    }
-
-    /// Iterate over the points in a counter-clockwise order
-    ///
-    /// The Linestring isn't changed, and the points are returned either in order, or in reverse
-    /// order, so that the resultant order makes it appear counter-clockwise
-    fn points_ccw(&self) -> Points<Self::Scalar> {
-        match self.winding_order() {
-            Some(WindingOrder::Clockwise) => Points(EitherIter::B(self.points_iter().rev())),
-            _ => Points(EitherIter::A(self.points_iter())),
-        }
-    }
-
-    /// Change this line's points so they are in clockwise winding order
-    fn make_cw_winding(&mut self) {
-        if let Some(WindingOrder::CounterClockwise) = self.winding_order() {
-            self.0.reverse();
-        }
-    }
-
-    /// Change this line's points so they are in counterclockwise winding order
-    fn make_ccw_winding(&mut self) {
-        if let Some(WindingOrder::Clockwise) = self.winding_order() {
-            self.0.reverse();
-        }
-    }
-}
-
-#[cfg(test)]
-mod test {
-    use super::*;
-    use crate::Point;
-
-    #[test]
-    fn robust_winding_float() {
-        // 3 points forming a triangle
-        let a = Point::new(0., 0.);
-        let b = Point::new(2., 0.);
-        let c = Point::new(1., 2.);
-
-        // Verify open linestrings return None
-        let mut ls = LineString::from(vec![a.0, b.0, c.0]);
-        assert!(
-            ls.winding_order()
-                .is_none()
-        );
-
-        ls.0.push(ls.0[0]);
-        assert_eq!(
-            ls.winding_order(),
-            Some(WindingOrder::CounterClockwise)
-        );
-
-        ls.make_cw_winding();
-        assert_eq!(
-            ls.winding_order(),
-            Some(WindingOrder::Clockwise)
-        );
-
-    }
-
-    #[test]
-    fn robust_winding_integer() {
-        // 3 points forming a triangle
-        let a = Point::new(0i64, 0);
-        let b = Point::new(2, 0);
-        let c = Point::new(1, 2);
-
-        // Verify open linestrings return None
-        let mut ls = LineString::from(vec![a.0, b.0, c.0]);
-        assert!(ls.winding_order().is_none());
-
-        ls.0.push(ls.0[0]);
-        assert!(ls.is_ccw());
-
-        let ccw_ls: Vec<_> = ls.points_ccw().collect();
-
-        ls.make_cw_winding();
-        assert!(ls.is_cw());
-
-        assert_eq!(
-            &ls.points_ccw().collect::<Vec<_>>(),
-            &ccw_ls,
-        );
-    }
-}
diff --git a/geo/src/algorithm/winding_order.rs b/geo/src/algorithm/winding_order.rs
index 1ec7804bc..a1c8919a0 100644
--- a/geo/src/algorithm/winding_order.rs
+++ b/geo/src/algorithm/winding_order.rs
@@ -1,46 +1,9 @@
+use super::kernels::*;
+use crate::{Point, CoordinateType, LineString};
 use crate::utils::EitherIter;
-use crate::{CoordinateType, LineString, Point};
 use geo_types::PointsIter;
 use std::iter::Rev;
 
-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
-}
-
 /// Iterates through a list of `Point`s
 pub struct Points<'a, T>(pub(crate) EitherIter<Point<T>, PointsIter<'a, T>, Rev<PointsIter<'a, T>>>)
 where
@@ -64,11 +27,9 @@ pub enum WindingOrder {
     CounterClockwise,
 }
 
-/// Calculate, and work with, the winding order
-pub trait Winding<T>
-where
-    T: CoordinateType,
-{
+pub trait Winding {
+    type Scalar: CoordinateType;
+
     /// Return the winding order of this object
     fn winding_order(&self) -> Option<WindingOrder>;
 
@@ -86,13 +47,13 @@ where
     ///
     /// The object isn't changed, and the points are returned either in order, or in reverse
     /// order, so that the resultant order makes it appear clockwise
-    fn points_cw(&self) -> Points<T>;
+    fn points_cw(&self) -> Points<Self::Scalar>;
 
     /// Iterate over the points in a counter-clockwise order
     ///
     /// The object isn't changed, and the points are returned either in order, or in reverse
     /// order, so that the resultant order makes it appear counter-clockwise
-    fn points_ccw(&self) -> Points<T>;
+    fn points_ccw(&self) -> Points<Self::Scalar>;
 
     /// Change this objects's points so they are in clockwise winding order
     fn make_cw_winding(&mut self);
@@ -119,31 +80,73 @@ where
     }
 }
 
-impl<T> Winding<T> for LineString<T>
-where
-    T: CoordinateType,
+
+impl<T, K> Winding for LineString<T>
+where T: CoordinateType + HasKernel<Ker = K>,
+      K: Kernel<Scalar = T>
 {
-    /// Returns the winding order of this line
-    /// None if the winding order is undefined.
+    type Scalar = T;
+
     fn winding_order(&self) -> Option<WindingOrder> {
-        let shoelace = twice_signed_ring_area(self);
-        if shoelace < T::zero() {
-            Some(WindingOrder::Clockwise)
-        } else if shoelace > T::zero() {
-            Some(WindingOrder::CounterClockwise)
-        } else if shoelace == T::zero() {
-            None
-        } else {
-            // make compiler stop complaining
-            unreachable!()
+        // 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 increment = |x: &mut usize| {
+            *x += 1;
+            if *x >= self.num_coords() { *x = 0; }
+        };
+
+        let decrement = |x: &mut usize| {
+            if *x == 0 {
+                *x = self.num_coords() - 1;
+            } else { *x -= 1; }
+        };
+
+        use crate::utils::lexicographically_least_index;
+        let i = lexicographically_least_index(&self.0);
+
+        let mut next = i;
+        increment(&mut next);
+        while self.0[next] == self.0[i] {
+            if next == i {
+                // We've looped too much. There aren't
+                // enough unique coords to compute orientation.
+                return None;
+            }
+            increment(&mut next);
+        }
+
+        let mut prev = i;
+        decrement(&mut prev);
+        while self.0[prev] == self.0[i] {
+            // Note: we don't need to check if prev == i as
+            // the previous loop succeeded, and so we have
+            // at least two distinct elements in the list
+            decrement(&mut prev);
+        }
+
+        match K::orient2d(
+            self.0[prev],
+            self.0[i],
+            self.0[next]
+        ) {
+            Orientation::CounterClockwise => Some(WindingOrder::CounterClockwise),
+            Orientation::Clockwise => Some(WindingOrder::Clockwise),
+            _ => None
         }
+
     }
 
     /// Iterate over the points in a clockwise order
     ///
     /// The Linestring isn't changed, and the points are returned either in order, or in reverse
     /// order, so that the resultant order makes it appear clockwise
-    fn points_cw(&self) -> Points<T> {
+    fn points_cw(&self) -> Points<Self::Scalar> {
         match self.winding_order() {
             Some(WindingOrder::CounterClockwise) => Points(EitherIter::B(self.points_iter().rev())),
             _ => Points(EitherIter::A(self.points_iter())),
@@ -154,7 +157,7 @@ where
     ///
     /// The Linestring isn't changed, and the points are returned either in order, or in reverse
     /// order, so that the resultant order makes it appear counter-clockwise
-    fn points_ccw(&self) -> Points<T> {
+    fn points_ccw(&self) -> Points<Self::Scalar> {
         match self.winding_order() {
             Some(WindingOrder::Clockwise) => Points(EitherIter::B(self.points_iter().rev())),
             _ => Points(EitherIter::A(self.points_iter())),
@@ -179,73 +182,58 @@ where
 #[cfg(test)]
 mod test {
     use super::*;
+    use crate::Point;
 
     #[test]
-    fn winding_order() {
+    fn robust_winding_float() {
         // 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])
-            .winding_order()
-            .is_none());
+        let mut ls = LineString::from(vec![a.0, b.0, c.0]);
+        assert!(
+            ls.winding_order()
+                .is_none()
+        );
 
-        assert_eq!(cw_line.winding_order(), Some(WindingOrder::Clockwise));
-        assert_eq!(cw_line.is_cw(), true);
-        assert_eq!(cw_line.is_ccw(), false);
+        ls.0.push(ls.0[0]);
         assert_eq!(
-            ccw_line.winding_order(),
+            ls.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.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();
+
+        ls.make_cw_winding();
         assert_eq!(
-            new_line2.winding_order(),
-            Some(WindingOrder::CounterClockwise)
+            ls.winding_order(),
+            Some(WindingOrder::Clockwise)
+        );
+
+    }
+
+    #[test]
+    fn robust_winding_integer() {
+        // 3 points forming a triangle
+        let a = Point::new(0i64, 0);
+        let b = Point::new(2, 0);
+        let c = Point::new(1, 2);
+
+        // Verify open linestrings return None
+        let mut ls = LineString::from(vec![a.0, b.0, c.0]);
+        assert!(ls.winding_order().is_none());
+
+        ls.0.push(ls.0[0]);
+        assert!(ls.is_ccw());
+
+        let ccw_ls: Vec<_> = ls.points_ccw().collect();
+
+        ls.make_cw_winding();
+        assert!(ls.is_cw());
+
+        assert_eq!(
+            &ls.points_ccw().collect::<Vec<_>>(),
+            &ccw_ls,
         );
-        assert_ne!(new_line2, cw_line);
-        assert_eq!(new_line2, ccw_line);
     }
 }

From 9f1af5412f3161c1268a0f632bea5c8fd8954e3e Mon Sep 17 00:00:00 2001
From: Rajsekar Manokaran <rajsekar@gmail.com>
Date: Tue, 25 Aug 2020 21:02:56 +0530
Subject: [PATCH 7/8] Minor fixes

+ add more docs
+ add more types to Kernel (i16, i128, isize)
---
 geo/Cargo.toml                   |  2 --
 geo/src/algorithm/kernels/mod.rs | 13 +++++++++++--
 geo/src/utils.rs                 |  6 ++++--
 3 files changed, 15 insertions(+), 6 deletions(-)

diff --git a/geo/Cargo.toml b/geo/Cargo.toml
index cf812bdfc..65026f813 100644
--- a/geo/Cargo.toml
+++ b/geo/Cargo.toml
@@ -34,8 +34,6 @@ use-serde = ["serde", "geo-types/serde"]
 [dev-dependencies]
 approx = "0.3.0"
 criterion = { version = "0.3" }
-png = "0.16.7"
-float_extras = "0.1.6"
 
 [[bench]]
 name = "area"
diff --git a/geo/src/algorithm/kernels/mod.rs b/geo/src/algorithm/kernels/mod.rs
index 418c2d88a..1febdc857 100644
--- a/geo/src/algorithm/kernels/mod.rs
+++ b/geo/src/algorithm/kernels/mod.rs
@@ -45,7 +45,11 @@ pub trait HasKernel: CoordinateType {
     type Ker: Kernel<Scalar = Self>;
 }
 
-#[macro_export]
+// 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 {
@@ -59,8 +63,13 @@ 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);
diff --git a/geo/src/utils.rs b/geo/src/utils.rs
index b2bb2d170..6d8126517 100644
--- a/geo/src/utils.rs
+++ b/geo/src/utils.rs
@@ -127,8 +127,10 @@ where
     }
 }
 
-/// Compute index of the lexicographically least point.
-/// Should only be called on a non-empty slice.
+/// Compute index of the lexicographically least point. In
+/// other words, point with minimum `x` coord, and breaking
+/// ties with minimum `y` coord. Should only be called on a
+/// non-empty slice with no `nan` coordinates.
 pub fn lexicographically_least_index<T: Copy + PartialOrd>(pts: &[T]) -> usize {
     assert!(pts.len() > 0);
 

From 51e8a113aa14a28f15bf4f5262c966e5ea7655c7 Mon Sep 17 00:00:00 2001
From: Rajsekar Manokaran <rajsekar@gmail.com>
Date: Tue, 25 Aug 2020 21:49:58 +0530
Subject: [PATCH 8/8] Remove unnecessary additions from PR

+ revert back to not deriving `PartialOrd` for `Coordinate`
+ remove unused `ConvexHull` trait in graham_hull impl
---
 geo-types/src/coordinate.rs      |  2 +-
 geo/src/algorithm/graham_hull.rs |  9 ++-------
 geo/src/utils.rs                 | 22 ++++++++--------------
 3 files changed, 11 insertions(+), 22 deletions(-)

diff --git a/geo-types/src/coordinate.rs b/geo-types/src/coordinate.rs
index 54755a290..d6f2f7b8f 100644
--- a/geo-types/src/coordinate.rs
+++ b/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, PartialOrd, Clone, Copy, Debug, Hash)]
+#[derive(Eq, PartialEq, Clone, Copy, Debug, Hash)]
 #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
 pub struct Coordinate<T>
 where
diff --git a/geo/src/algorithm/graham_hull.rs b/geo/src/algorithm/graham_hull.rs
index 0ff694a54..df055351a 100644
--- a/geo/src/algorithm/graham_hull.rs
+++ b/geo/src/algorithm/graham_hull.rs
@@ -1,4 +1,4 @@
-use crate::{Coordinate, CoordinateType, LineString};
+use crate::{Coordinate, LineString};
 use super::kernels::*;
 
 // Utility function: swap idx to head(0th position), remove
@@ -87,15 +87,10 @@ where
     output
 }
 
-pub trait ConvexHull {
-    type Scalar: CoordinateType;
-    fn convex_hull(&self) -> LineString<Self::Scalar>;
-}
-
-
 #[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] {
diff --git a/geo/src/utils.rs b/geo/src/utils.rs
index 6d8126517..5caacc3e8 100644
--- a/geo/src/utils.rs
+++ b/geo/src/utils.rs
@@ -1,6 +1,7 @@
 //! Internal utility functions, types, and data structures.
 
 use crate::contains::Contains;
+use geo_types::{CoordinateType, Coordinate};
 
 /// Partition a mutable slice in-place so that it contains all elements for
 /// which `predicate(e)` is `true`, followed by all elements for which
@@ -131,21 +132,14 @@ where
 /// other words, point with minimum `x` coord, and breaking
 /// ties with minimum `y` coord. Should only be called on a
 /// non-empty slice with no `nan` coordinates.
-pub 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) )
-        }
-    }
+pub fn lexicographically_least_index<T: CoordinateType>(pts: &[Coordinate<T>]) -> usize {
+
+    pts.iter().enumerate().min_by(
+        |(_, p), (_, q)| p.x.partial_cmp(&q.x).unwrap().then(
+            p.y.partial_cmp(&q.y).unwrap()
+        )
+    ).unwrap().0
 
-    min.unwrap().0
 }
 
 #[cfg(test)]