Convert Delaunay triangulation to graph

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<NodeData, EdgeWeight, Directed, NodeIndex>;
+pub type Graph = petgraph::Graph<NodeData, EdgeWeight, Undirected, NodeIndex>;
+
 /// Calculate the Delaunay triangulation of the given point cloud
 pub fn triangulate(points: impl Iterator<Item = geo::Point>) -> Triangulation {
     let points: Vec<delaunator::Point> = 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);
+    }
 }

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<W, G>(mut writer: W, geometries: G)
+pub(crate) fn write_wkt_geometries<W, G>(mut writer: W, geometries: G)
 where
     W: Write,
     G: IntoIterator<Item = Geometry<f64>>,