Fix bug in CoordPos logic for MonoPoly

+ inc PR suggestions; thanks @frewsxcv
thehappycheese · Jul 19, 2023 · 570c9bf · 570c9bf
1 parent 6532ccc
commit 570c9bf
Show file tree

Hide file tree

Showing 4 changed files with 51 additions and 20 deletions.
diff --git a/geo/benches/utils/random.rs b/geo/benches/utils/random.rs
@@ -1,7 +1,7 @@
 #![allow(unused)]
 use std::f64::consts::PI;
 
-use geo::algorithm::{ConcaveHull, ConvexHull, MapCoords, Rotate};
+use geo::algorithm::{BoundingRect, ConcaveHull, ConvexHull, MapCoords, Rotate};
 use geo::geometry::*;
 
 use rand::{thread_rng, Rng};
@@ -76,29 +76,50 @@ pub fn circular_polygon<R: Rng>(mut rng: R, steps: usize) -> Polygon<f64> {
 pub fn steppy_polygon<R: Rng>(mut rng: R, steps: usize) -> Polygon<f64> {
     let mut ring = Vec::with_capacity(2 * steps);
 
-    let y_step = 1.0;
+    let y_step = 10.0;
     let nudge_std = y_step / 1000.0;
     let mut y = 0.0;
     let normal = Normal::new(0.0, nudge_std * nudge_std).unwrap();
-    let x_shift = 50.0;
+    let x_shift = 100.0;
 
     ring.push((0.0, 0.0).into());
     (0..steps).for_each(|_| {
-        let x: f64 = rng.sample::<f64, _>(Standard) * x_shift / 2.;
+        let x: f64 = rng.sample::<f64, _>(Standard);
         y += y_step;
         ring.push((x, y).into());
     });
-    ring.push((5. * x_shift, y).into());
+    ring.push((x_shift, y).into());
     (0..steps).for_each(|_| {
-        let x: f64 = rng.sample::<f64, _>(Standard) * x_shift;
+        let x: f64 = rng.sample::<f64, _>(Standard);
         y -= y_step;
         // y += normal.sample(&mut rng);
         ring.push((x_shift + x, y).into());
     });
 
-    Polygon::new(LineString(ring), vec![])
+    normalize_polygon(Polygon::new(LineString(ring), vec![]))
+}
+
+/// Normalizes polygon to fit and fill `[-1, 1] X [-1, 1]` square.
+///
+/// Uses `MapCoord` and `BoundingRect`
+pub fn normalize_polygon(poly: Polygon<f64>) -> Polygon<f64> {
+    let bounds = poly.bounding_rect().unwrap();
+    let dims = bounds.max() - bounds.min();
+    let x_scale = 2. / dims.x;
+    let y_scale = 2. / dims.y;
+
+    let x_shift = -bounds.min().x * x_scale - 1.;
+    let y_shift = -bounds.min().y * y_scale - 1.;
+    poly.map_coords(|mut c| {
+        c.x *= x_scale;
+        c.x += x_shift;
+        c.y *= y_scale;
+        c.y += y_shift;
+        c
+    })
 }
 
+#[derive(Debug, Clone)]
 pub struct Samples<T>(Vec<T>);
 impl<T> Samples<T> {
     pub fn sampler<'a>(&'a self) -> impl FnMut() -> &'a T {
@@ -113,4 +134,8 @@ impl<T> Samples<T> {
     pub fn from_fn<F: FnMut() -> T>(size: usize, mut proc: F) -> Self {
         Self((0..size).map(|_| proc()).collect())
     }
+
+    pub fn map<U, F: FnMut(T) -> U>(self, mut proc: F) -> Samples<U> {
+        Samples(self.0.into_iter().map(proc).collect())
+    }
 }
diff --git a/geo/src/algorithm/monotone/builder.rs b/geo/src/algorithm/monotone/builder.rs
@@ -358,11 +358,11 @@ impl<T: GeoNum> Chain<T> {
     }
 
     pub fn finish_with(self, other: Self) -> MonoPoly<T> {
-        // assert!(
-        //     self.0 .0[0] == other.0 .0[0]
-        //         && self.0 .0.last().unwrap() == other.0 .0.last().unwrap(),
-        //     "chains must finish with same start/end points"
-        // );
+        assert!(
+            self.0 .0[0] == other.0 .0[0]
+                && self.0 .0.last().unwrap() == other.0 .0.last().unwrap(),
+            "chains must finish with same start/end points"
+        );
         debug!("finishing {self:?} with {other:?}");
         MonoPoly::new(other.0, self.0)
     }

diff --git a/geo/src/algorithm/monotone/mod.rs b/geo/src/algorithm/monotone/mod.rs
@@ -25,7 +25,7 @@ pub use builder::monotone_subdivision;
 /// # Example
 ///
 /// Construct a `MonotonicPolygons` from a `Polygon`, or a `MultiPolygon` using
-/// `MontonicPolygons::from`, and` query point intersection via the
+/// `MontonicPolygons::from`, and query point intersection via the
 /// `Intersects<Coord>` trait.
 ///
 /// ```rust

diff --git a/geo/src/algorithm/monotone/mono_poly.rs b/geo/src/algorithm/monotone/mono_poly.rs
@@ -50,10 +50,9 @@ impl<T: GeoNum> MonoPoly<T> {
     /// strictly above the bottom chain except at the end-points.  Not all these
     /// conditions are checked, and the algorithm may panic if they are not met.
     pub(super) fn new(top: LineString<T>, bot: LineString<T>) -> Self {
-        // TODO: move these to debug-only asserts
-        assert_eq!(top.0.first(), bot.0.first());
-        assert_eq!(top.0.last(), bot.0.last());
-        assert_ne!(top.0.first(), top.0.last());
+        debug_assert_eq!(top.0.first(), bot.0.first());
+        debug_assert_eq!(top.0.last(), bot.0.last());
+        debug_assert_ne!(top.0.first(), top.0.last());
         let bounds = get_bounding_rect(top.0.iter().chain(bot.0.iter()).cloned()).unwrap();
         Self { top, bot, bounds }
     }
@@ -77,7 +76,8 @@ impl<T: GeoNum> MonoPoly<T> {
 
     /// Get the pair of segments in the chain that intersects the line parallel
     /// to the Y-axis at the given x-coordinate.  Ties are broken by picking the
-    /// segment with smaller x coordinates.
+    /// segment with lower index, i.e. the segment closer to the start of the
+    /// chains.
     pub fn bounding_segment(&self, x: T) -> Option<(Line<T>, Line<T>)> {
         // binary search for the segment that contains the x coordinate.
         let tl_idx = self.top.0.partition_point(|c| c.x < x);
@@ -88,11 +88,17 @@ impl<T: GeoNum> MonoPoly<T> {
         if bl_idx == 0 {
             debug_assert_eq!(tl_idx, 0);
             debug_assert_eq!(self.bot.0[0].x, x);
+            return Some((
+                Line::new(self.top.0[0], self.top.0[1]),
+                Line::new(self.bot.0[0], self.bot.0[1]),
+            ));
+        } else {
+            debug_assert_ne!(tl_idx, 0);
         }
 
         Some((
-            Line::new(self.top.0[tl_idx], self.top.0[tl_idx + 1]),
-            Line::new(self.bot.0[bl_idx], self.bot.0[bl_idx + 1]),
+            Line::new(self.top.0[tl_idx - 1], self.top.0[tl_idx]),
+            Line::new(self.bot.0[bl_idx - 1], self.bot.0[bl_idx]),
         ))
     }