From 805aa839194fccf306815cf62fcf8af20fb3581d Mon Sep 17 00:00:00 2001 From: Austin Gill Date: Sun, 24 Dec 2023 17:52:55 -0600 Subject: [PATCH] WIP: Add point-snapping utilities --- Cargo.toml | 7 ++ generative/lib.rs | 1 + generative/snap.rs | 281 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 289 insertions(+) create mode 100644 generative/snap.rs diff --git a/Cargo.toml b/Cargo.toml index 61cf6ad..c9c218b 100644 --- a/Cargo.toml +++ b/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 = [] diff --git a/generative/lib.rs b/generative/lib.rs index 9677b4a..7940198 100644 --- a/generative/lib.rs +++ b/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; diff --git a/generative/snap.rs b/generative/snap.rs new file mode 100644 index 0000000..bb94c21 --- /dev/null +++ b/generative/snap.rs @@ -0,0 +1,281 @@ +use geo::{Coord, Geometry}; +use kdtree::distance::squared_euclidean; +use kdtree::KdTree; + +use crate::flatten::{flatten_geometries_into_points_ref, flatten_nested_geometries}; +use crate::MapCoordsInPlaceMut; + +pub type GeomKdTree = KdTree; + +#[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, + strategy: SnappingStrategy, +) -> Box> { + // 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) +} + +// NOTE: This can change the geometry type! +fn filter_duplicate_vertices(geom: Geometry) -> Geometry { + // todo!() + geom +} + +// TODO: Snap graphs +// use petgraph::graph::NodeIndex; +// use crate::graph::GeometryGraph; +// pub type GraphKdTree = KdTree; +// // TODO: Move or reference? +// pub fn snap_graph(graph: GeometryGraph, strategy: SnappingStrategy) -> GeometryGraph { +// graph +// } + +// pub fn snap_graph_immutable_input( +// graph: GeometryGraph, +// index: &GeomKdTree, +// ) -> GeometryGraph { +// graph +// } + +// pub fn snap_graph_mutable_output( +// graph: GeometryGraph, +// index: &mut GeomKdTree, +// ) -> GeometryGraph { +// 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::LineString(LineString::new(vec![ + Coord { x: 0.1, y: 0.0 }, // TODO: This should be removed + 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); + } +}