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.

Author: rmanoka

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"

geo/src/algorithm/kernels/mod.rs

@@ -1,39 +1,44 @@
 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>,
         r: Coordinate<Self::Scalar>,
     ) -> 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;

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>);
 

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>);
 

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, f64::INFINITY);
+            data.push(<T as HasKernel>::Ker::orient2d(p1, p2, p3));
+        }
+        yd2 = nextafter(yd2, 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);
+}