WIP: Add point-snapping utilities for geometries

Cargo.toml

@@ -62,6 +62,10 @@ path = "tools/pack.rs"
 name = "streamline"
 path = "tools/streamline.rs"
 
+[[bin]]
+name = "snap"
+path = "tools/snap.rs"
+
 [dependencies]
 clap = {version="4.0", features=["derive"]}
 delaunator = "1.0"
@@ -87,5 +91,8 @@ cmake = "0.1"
 fs_extra = "1.3"
 glob = "0.3"
 
+[dev-dependencies]
+float-cmp = "0.9"
+
 [features]
 test-io = []

generative/lib.rs

@@ -3,6 +3,7 @@ pub mod flatten;
 mod geometry_mut_map;
 pub mod graph;
 pub mod io;
+pub mod snap;
 pub mod triangulation;
 
 pub use geometry_mut_map::MapCoordsInPlaceMut;

generative/snap.rs

@@ -0,0 +1,322 @@
+use geo::{Coord, CoordsIter, Geometry, Line, LineString, Polygon, Triangle};
+use kdtree::distance::squared_euclidean;
+use kdtree::KdTree;
+
+// use petgraph::graph::NodeIndex;
+use crate::flatten::{flatten_geometries_into_points_ref, flatten_nested_geometries};
+// use crate::graph::GeometryGraph;
+use crate::MapCoordsInPlaceMut;
+
+pub type GeomKdTree = KdTree<f64, Coord, [f64; 2]>;
+// pub type GraphKdTree = KdTree<f64, NodeIndex, [f64; 2]>;
+
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub enum SnappingStrategy {
+    /// Snap to the closest input point
+    ///
+    /// Not sensitive to the snapping order, but may give less dramatic results
+    ClosestInputPoint(f64),
+    /// Snap to the closest point, modifying the point cloud as the algorithm progresses
+    ///
+    /// Sensitive to the order in which snapping is performed
+    ClosestOutputPoint(f64),
+    /// Snap points to a regular grid, instead of themselves
+    RegularGrid(f64),
+}
+
+pub fn snap_geoms(
+    geoms: impl Iterator<Item = Geometry>,
+    strategy: SnappingStrategy,
+) -> Box<dyn Iterator<Item = Geometry>> {
+    // Build a k-d tree from the given geometries. flatten_geometries_into_points would require
+    // cloning all of the given geometries, so flatten first, and then use the _ref() variant.
+    let geoms = flatten_nested_geometries(geoms);
+    let geoms: Vec<_> = geoms.collect();
+
+    // Short circuit the creation of the k-d tree
+    if let SnappingStrategy::RegularGrid(tolerance) = strategy {
+        let snapped = geoms.into_iter().map(move |g| snap_geom_grid(g, tolerance));
+        return Box::new(snapped);
+    }
+
+    let points = flatten_geometries_into_points_ref(geoms.iter());
+    let mut index = GeomKdTree::new(2);
+    for point in points {
+        let coord: Coord = point.into();
+        let coords = [point.x(), point.y()];
+        let closest = index.nearest(&coords, 1, &squared_euclidean).unwrap();
+        if let Some(closest) = closest.first() {
+            let (distance, _) = closest;
+            if *distance == 0.0 {
+                // This coordinate has already been added
+                continue;
+            }
+        }
+        // Since you don't get the coordinates back when doing nearest neighbor lookups, we
+        // need to store the coordinates again in the point data.
+        index.add(coords, coord).unwrap();
+    }
+
+    // TODO: Should this filter out duplicates introduced by snapping?
+    let snapped = geoms
+        .into_iter()
+        .map(move |g| snap_geom(g, &mut index, &strategy));
+    Box::new(snapped)
+}
+
+pub fn snap_geom(geom: Geometry, index: &mut GeomKdTree, strategy: &SnappingStrategy) -> Geometry {
+    match strategy {
+        SnappingStrategy::ClosestInputPoint(tolerance) => {
+            // TODO: This isn't the best algorithm, because it biases points away from themselves,
+            // meaning two points in range will always _SWITCH PLACES_ instead of _SNAP TOGETHER_.
+            // I think I'd need to maintain an index of snapped vertices, so that if I find a
+            // vertex there, I know I don't need to snap it again.
+            snap_geom_impl(geom, index, *tolerance, false)
+        }
+        SnappingStrategy::ClosestOutputPoint(tolerance) => {
+            snap_geom_impl(geom, index, *tolerance, true)
+        }
+        SnappingStrategy::RegularGrid(tolerance) => snap_geom_grid(geom, *tolerance),
+    }
+}
+
+fn snap_geom_impl(
+    mut geom: Geometry,
+    index: &mut GeomKdTree,
+    tolerance: f64,
+    remove_snapped_point: bool,
+) -> Geometry {
+    geom.map_coords_in_place_mut(|c| snap_coord(c, index, tolerance, remove_snapped_point));
+    filter_duplicate_vertices(geom)
+}
+
+fn snap_coord(
+    coord: Coord,
+    index: &mut GeomKdTree,
+    tolerance: f64,
+    remove_snapped_point: bool,
+) -> Coord {
+    // Find the closest two points in the index, because the first closest should always be ourself.
+    let coords = [coord.x, coord.y];
+    let neighbors = index
+        .within(&coords, tolerance, &squared_euclidean)
+        .unwrap();
+    // We should always find ourselves, or, if move_snapped_point is true, at least find where
+    // ourselves have already been snapped to (because one point in the kd-tree could be multiple
+    // vertices from multiple geometries).
+    debug_assert!(!neighbors.is_empty());
+
+    if !neighbors.is_empty() {
+        let (mut _distance, mut found_coords) = neighbors[0];
+        // We found ourselves. Now look for a neighbor in range
+        if found_coords == &coord && neighbors.len() > 1 {
+            // The next closest point
+            (_distance, found_coords) = neighbors[1];
+        }
+
+        let snapped_coord = *found_coords;
+        if remove_snapped_point {
+            index.remove(&coords, &coord).unwrap();
+        }
+
+        return snapped_coord;
+    }
+
+    coord
+}
+
+fn snap_coord_grid(coord: Coord, tolerance: f64) -> Coord {
+    Coord {
+        x: snap_f64_grid(coord.x, tolerance),
+        y: snap_f64_grid(coord.y, tolerance),
+    }
+}
+
+fn snap_f64_grid(value: f64, tolerance: f64) -> f64 {
+    // This was waaaay harder than I expected :(
+    let rem = value.rem_euclid(tolerance);
+    let floor = value - rem;
+
+    // Need to round away from zero, which takes special handling for negative and positive values
+    let distance = value - floor;
+    let half_tol = 0.5 * tolerance;
+
+    if (value < 0.0 && distance <= half_tol) || distance < half_tol {
+        return floor;
+    }
+
+    floor + tolerance
+}
+
+fn snap_geom_grid(mut geom: Geometry, tolerance: f64) -> Geometry {
+    geom.map_coords_in_place_mut(|c| snap_coord_grid(c, tolerance));
+    filter_duplicate_vertices(geom)
+}
+
+fn filter_duplicate_vertices(geom: Geometry) -> Geometry {
+    match geom {
+        Geometry::Point(_) => geom,
+        Geometry::Line(l) => {
+            if l.start == l.end {
+                Geometry::Point(l.start.into())
+            } else {
+                geom
+            }
+        }
+        Geometry::LineString(ls) => filter_duplicate_ls(ls, false),
+        Geometry::Polygon(p) => {
+            let (exterior, interiors) = p.into_inner();
+            let exterior = filter_duplicate_ls(exterior, true);
+            match exterior {
+                Geometry::Point(_) => exterior,
+                Geometry::Line(_) => exterior,
+                Geometry::LineString(exterior) => {
+                    let mut filtered_interiors = Vec::new();
+                    for interior in interiors {
+                        if let Geometry::LineString(interior) = filter_duplicate_ls(interior, true)
+                        {
+                            filtered_interiors.push(interior);
+                        }
+                    }
+                    Geometry::Polygon(Polygon::new(exterior, filtered_interiors))
+                }
+                _ => unreachable!(),
+            }
+        }
+        Geometry::Rect(r) => filter_duplicate_vertices(Geometry::Polygon(r.to_polygon())),
+        Geometry::Triangle(t) => {
+            let mut unique = Vec::new();
+            for c in t.to_array() {
+                if !unique.contains(&c) {
+                    unique.push(c);
+                }
+            }
+            if unique.len() == 3 {
+                geom
+            } else if unique.len() == 2 {
+                Geometry::Line(Line::new(unique[0], unique[1]))
+            } else {
+                Geometry::Point(unique[0].into())
+            }
+        }
+        _ => unreachable!("flatten_nested_geometries in the call graph prevents MULTI-geometries"),
+    }
+}
+
+fn filter_duplicate_ls(mut ls: LineString, closed: bool) -> Geometry {
+    ls.0.dedup();
+
+    if ls.coords_count() == 1 {
+        Geometry::Point(ls.0[0].into())
+    } else if ls.coords_count() == 2 {
+        Geometry::Line(Line::new(ls.0[0], ls.0[1]))
+    } else if closed && ls.coords_count() == 3 {
+        Geometry::Triangle(Triangle::new(ls.0[0], ls.0[1], ls.0[2]))
+    } else {
+        Geometry::LineString(ls)
+    }
+}
+
+// pub fn snap_graph<D>(graph: GeometryGraph<D>, strategy: SnappingStrategy) -> GeometryGraph<D> {
+//     graph
+// }
+
+#[cfg(test)]
+mod tests {
+    use float_cmp::assert_approx_eq;
+    use geo::{LineString, Point};
+
+    use super::*;
+
+    #[test]
+    fn test_f64_snapping() {
+        let values = [
+            // value, tolerance, expected
+            (0.0, 1.0, 0.0),
+            (0.4, 1.0, 0.0),
+            (0.5, 1.0, 1.0),
+            (1.0, 1.0, 1.0),
+            (-0.0, 1.0, 0.0),
+            (-0.4, 1.0, 0.0),
+            (-0.5, 1.0, -1.0),
+            (-1.5, 1.0, -2.0), // round away from 0.0
+            (-1.0, 1.0, -1.0),
+            (0.0, 0.5, 0.0),
+            (0.24, 0.5, 0.0),
+            (0.25, 0.5, 0.5),
+            (-0.25, 0.5, -0.5), // round away from 0.0
+            (1.4, 0.5, 1.5),
+            (1.2, 0.5, 1.0),
+            (-1.4, 0.5, -1.5),
+            (-1.2, 0.5, -1.0),
+        ];
+        for (value, tolerance, expected) in values {
+            let actual = snap_f64_grid(value, tolerance);
+            assert_approx_eq!(f64, actual, expected);
+        }
+    }
+
+    #[test]
+    fn test_coord_grid_snapping() {
+        let coord = Coord { x: 1.5, y: 0.5 };
+        let expected = Coord { x: 2.0, y: 1.0 };
+        let actual = snap_coord_grid(coord, 1.0);
+        assert_eq!(actual, expected);
+
+        let coord = Coord { x: 1.4, y: -0.6 };
+        let expected = Coord { x: 1.0, y: -1.0 };
+        let actual = snap_coord_grid(coord, 1.0);
+        assert_eq!(actual, expected);
+
+        let coord = Coord { x: 1.4, y: 0.4 };
+        let expected = Coord { x: 1.5, y: 0.5 };
+        let actual = snap_coord_grid(coord, 0.5);
+        assert_eq!(actual, expected);
+    }
+
+    #[test]
+    fn test_snap_two_points() {
+        let geoms = [
+            Geometry::Point(Point::new(0.0, 0.0)),
+            Geometry::Point(Point::new(0.5, 0.0)),
+        ];
+        let expected = [
+            Geometry::Point(Point::new(0.5, 0.0)),
+            Geometry::Point(Point::new(0.0, 0.0)), // Order reversed because they both snapped to
+                                                   // each other
+        ];
+
+        let actual: Vec<_> =
+            snap_geoms(geoms.into_iter(), SnappingStrategy::ClosestInputPoint(0.6)).collect();
+        assert_eq!(actual, expected);
+
+        let geoms = [
+            Geometry::Point(Point::new(0.0, 0.0)),
+            Geometry::Point(Point::new(0.5, 0.0)),
+        ];
+        let expected = [
+            Geometry::Point(Point::new(0.5, 0.0)),
+            Geometry::Point(Point::new(0.5, 0.0)),
+        ];
+
+        let actual: Vec<_> =
+            snap_geoms(geoms.into_iter(), SnappingStrategy::ClosestOutputPoint(0.6)).collect();
+        assert_eq!(actual, expected);
+    }
+
+    #[test]
+    fn test_points_on_linestring_snap_together() {
+        let geoms = [Geometry::LineString(LineString::new(vec![
+            Coord { x: 0.0, y: 0.0 },
+            Coord { x: 0.1, y: 0.0 },
+            Coord { x: 1.0, y: 0.0 },
+        ]))];
+        let expected = [Geometry::Line(Line::new(
+            Coord { x: 0.1, y: 0.0 },
+            Coord { x: 1.0, y: 0.0 },
+        ))];
+        let actual: Vec<_> =
+            snap_geoms(geoms.into_iter(), SnappingStrategy::ClosestOutputPoint(0.5)).collect();
+        assert_eq!(actual, expected);
+    }
+}