From f44754bd1e64c4367fb1fcb3e979738949971a6d Mon Sep 17 00:00:00 2001 From: Austin Gill Date: Mon, 28 Nov 2022 17:07:46 -0600 Subject: [PATCH] Convert Delaunay triangulation to graph --- generative/triangulation.rs | 161 +++++++++++++++++++++++++++++++++++- generative/wkio.rs | 2 +- 2 files changed, 160 insertions(+), 3 deletions(-) diff --git a/generative/triangulation.rs b/generative/triangulation.rs index 8651447..30162dd 100644 --- a/generative/triangulation.rs +++ b/generative/triangulation.rs @@ -1,3 +1,12 @@ +use geo::Point; +use petgraph::{Directed, Undirected}; + +type NodeData = Point; +type EdgeWeight = (); +type NodeIndex = usize; +pub type DiGraph = petgraph::Graph; +pub type Graph = petgraph::Graph; + /// Calculate the Delaunay triangulation of the given point cloud pub fn triangulate(points: impl Iterator) -> Triangulation { let points: Vec = points @@ -67,14 +76,86 @@ impl Triangulation { geo::Polygon::new(hull, vec![]) } - // TODO: Create the graph directly instead of relying on geom2graph? + pub fn digraph(&self) -> DiGraph { + let nodes = self.points.len(); + let edges = self.triangulation.halfedges.len(); + let mut graph = DiGraph::with_capacity(nodes, edges); + + // Add all the nodes + for (_i, point) in self.points.iter().enumerate() { + let point = Point::new(point.x, point.y); + // NOTE: It's important that the _node_index is the same as the index into the + // self.points array! + let _node_index = graph.add_node(point); + // debug_assert_eq!(_node_index.index(), _i); + } + + // Add the hull edges + for window in self.triangulation.hull.windows(2) { + let curr = window[0]; + let next = window[1]; + + graph.add_edge(curr.into(), next.into(), ()); + } + // NOTE: The hull is open and needs to be closed in order to capture the last edge! + if let (Some(first), Some(last)) = ( + self.triangulation.hull.first(), + self.triangulation.hull.last(), + ) { + graph.add_edge( + petgraph::graph::NodeIndex::new(*last), + petgraph::graph::NodeIndex::new(*first), + (), + ); + } + + // Add the interior half-edges + for (dst_i, src_i) in self.triangulation.halfedges.iter().enumerate() { + // This is a hull edge. The half-edges array doesn't contain enough information to + // build the graph edge, which is why we loop over the hull above. + if *src_i == delaunator::EMPTY { + continue; + } + + // dst_i and src_i are indices into the triangles array, which itself contains indices + // into the nodes array. + let src = self.triangulation.triangles[*src_i]; + let dst = self.triangulation.triangles[dst_i]; + + graph.add_edge(src.into(), dst.into(), ()); + } + + graph + } + + pub fn graph(&self) -> Graph { + let digraph = self.digraph(); + let nodes = self.points.len(); + let directed_edges = self.triangulation.halfedges.len(); + let mut graph = Graph::with_capacity(nodes, directed_edges - self.triangulation.hull.len()); + + // Add the nodes + for (_i, node) in digraph.raw_nodes().iter().enumerate() { + let _node_index = graph.add_node(node.weight); + // debug_assert_eq!(_i, _node_index.index()); + } + + // Add the edges. Use update_edge() to avoid duplicates + for edge in digraph.raw_edges() { + #[allow(clippy::unit_arg)] + let _edge_index = graph.update_edge(edge.source(), edge.target(), edge.weight); + } + + graph + } } #[cfg(test)] mod tests { use super::*; use crate::flatten::flatten_geometries_into_points_ref; - use crate::wkio::read_wkt_geometries; + use crate::wkio::{read_wkt_geometries, write_wkt_geometries}; + use delaunator::EMPTY; #[test] fn test_triangulate_unit_square() { @@ -110,4 +191,80 @@ mod tests { let hull = triangulation.hull(); assert_eq!(geo::Geometry::Polygon(hull), geometries[0]); } + + #[test] + fn test_graph() { + let wkt = b"POINT (65.85186826230156 -39.36525618186133)\n\ + POINT (61.35756898870892 -34.85194590696902)\n\ + POINT (38.25507241174806 -24.06029365638358)\n\ + POINT (9.25896849065506 -63.356505724778266)\n\ + POINT (-5.692741678486288 -17.181741298068346)\n\ + POINT (-32.93567551272198 -61.274655097506745)\n"; + let geometries: Vec<_> = read_wkt_geometries(&wkt[..]).collect(); + assert_eq!(geometries.len(), 6); + let points = flatten_geometries_into_points_ref(geometries.iter()); + let triangulation = triangulate(points); + + // NOTE: This all makes much more sense if you draw a picture! Pipe the following through + // render.py: + // cargo test test_graph -- --nocapture | ./tools/render.py + let lines = triangulation.lines().map(geo::Geometry::Line); + write_wkt_geometries(std::io::stdout(), geometries); + write_wkt_geometries(std::io::stdout(), lines); + + let triangles = [3, 5, 4, 4, 2, 3, 2, 1, 3, 1, 0, 3]; + assert_eq!(triangulation.triangulation.triangles, triangles); + let halfedges = [EMPTY, EMPTY, 5, EMPTY, 8, 2, EMPTY, 11, 4, EMPTY, EMPTY, 7]; + assert_eq!(triangulation.triangulation.halfedges, halfedges); + + // NOTE: None of the indices from 'triangles' or 'halfedges' index into 'hull'! You have to + // create the hull half-edges by looping over the hull. However, be aware that the hull is + // an open LINESTRING, and needs to be implicitly closed in order to capture the last edge! + + let hull = [1, 0, 3, 5, 4, 2]; + assert_eq!(triangulation.triangulation.hull, hull); + + let graph = triangulation.digraph(); + // It's not necessary that the edges be compared in order, but that's easiest to implement + // here. + let edges: Vec<_> = graph + .raw_edges() + .iter() + .map(|e| (e.source().index(), e.target().index())) + .collect(); + let expected = vec![ + (1, 0), + (0, 3), + (3, 5), + (5, 4), + (4, 2), + (2, 1), + (3, 4), + (3, 2), + (4, 3), + (3, 1), + (2, 3), + (1, 3), + ]; + assert_eq!(edges, expected); + + let graph = triangulation.graph(); + let edges: Vec<_> = graph + .raw_edges() + .iter() + .map(|e| (e.source().index(), e.target().index())) + .collect(); + let expected = vec![ + (1, 0), + (0, 3), + (3, 5), + (5, 4), + (4, 2), + (2, 1), + (3, 4), + (3, 2), + (3, 1), + ]; + assert_eq!(edges, expected); + } } diff --git a/generative/wkio.rs b/generative/wkio.rs index 1cbebb0..62718bc 100644 --- a/generative/wkio.rs +++ b/generative/wkio.rs @@ -206,7 +206,7 @@ where /// Write the given geometries with the given Writer in WKT format /// /// Each geometry will be written on its own line. -fn write_wkt_geometries(mut writer: W, geometries: G) +pub(crate) fn write_wkt_geometries(mut writer: W, geometries: G) where W: Write, G: IntoIterator>,