@@ -41,6 +41,7 @@ pub fn snap_geoms(
4141 }
4242
4343 let points = flatten_geometries_into_points_ref ( geoms. iter ( ) ) ;
44+ // Build the k-d tree, filtering out duplicate points as we go
4445 let mut index = GeomKdTree :: new ( 2 ) ;
4546 for point in points {
4647 let coord: Coord = point. into ( ) ;
@@ -79,32 +80,26 @@ fn snap_geom_impl(mut geom: Geometry, index: &mut GeomKdTree, tolerance: f64) ->
7980 filter_duplicate_vertices ( geom)
8081}
8182
82- fn snap_coord ( coord : Coord , index : & mut GeomKdTree , tolerance : f64 ) -> Coord {
83- // Find the closest two points in the index, because the first closest should always be ourself.
84- let coords = [ coord. x , coord. y ] ;
85- let neighbors = index
86- . within ( & coords, tolerance, & squared_euclidean)
87- . unwrap ( ) ;
88- // We should always find ourselves, or, if move_snapped_point is true, at least find where
89- // ourselves have already been snapped to (because one point in the kd-tree could be multiple
90- // vertices from multiple geometries).
91- debug_assert ! ( !neighbors. is_empty( ) ) ;
83+ fn snap_coord ( to_snap : Coord , index : & mut GeomKdTree , tolerance : f64 ) -> Coord {
84+ let query = [ to_snap. x , to_snap. y ] ;
85+ let neighbors = index. within ( & query, tolerance, & squared_euclidean) . unwrap ( ) ;
9286
9387 if !neighbors. is_empty ( ) {
9488 let ( mut _distance, mut found_coords) = neighbors[ 0 ] ;
95- // We found ourselves. Now look for a neighbor in range
96- if found_coords == & coord && neighbors. len ( ) > 1 {
89+ // If we found ourselves, snap to the next closest point
90+ if found_coords == & to_snap && neighbors. len ( ) > 1 {
9791 // The next closest point
9892 ( _distance, found_coords) = neighbors[ 1 ] ;
9993 }
10094
95+ // Remove the point that we snapped to, so that future snaps don't find it again
10196 let snapped_coord = * found_coords;
102- index. remove ( & coords , & coord ) . unwrap ( ) ;
97+ index. remove ( & query , & to_snap ) . unwrap ( ) ;
10398
104- return snapped_coord;
99+ snapped_coord
100+ } else {
101+ to_snap
105102 }
106-
107- coord
108103}
109104
110105fn snap_coord_grid ( coord : Coord , tolerance : f64 ) -> Coord {
@@ -210,6 +205,15 @@ where
210205 for node_idx in graph. node_indices ( ) {
211206 let node = graph[ node_idx] ;
212207 let coords = [ node. 0 . x , node. 0 . y ] ;
208+
209+ // Don't add duplicate vertices to the index
210+ let closest = index. nearest ( & coords, 1 , & squared_euclidean) . unwrap ( ) ;
211+ if let Some ( closest) = closest. first ( ) {
212+ let ( distance, _) = closest;
213+ if * distance == 0.0 {
214+ continue ;
215+ }
216+ }
213217 index. add ( coords, node_idx) . unwrap ( ) ;
214218 }
215219
@@ -231,8 +235,8 @@ where
231235{
232236 let mut nodes_to_remove = Vec :: new ( ) ;
233237 for node in graph. node_indices ( ) {
234- if let Some ( snapped ) = snap_graph_node ( & mut graph, node, index, tolerance) {
235- nodes_to_remove. push ( snapped ) ;
238+ if let Some ( _snapped_to ) = snap_graph_node ( & mut graph, node, index, tolerance) {
239+ nodes_to_remove. push ( node ) ;
236240 }
237241 }
238242
@@ -246,6 +250,9 @@ where
246250 graph
247251}
248252
253+ /// Snap the given node to the closest other node within the given tolerance
254+ ///
255+ /// Returns the NodeIndex that was snapped to, if the given node was snapped.
249256fn snap_graph_node < D > (
250257 graph : & mut GeometryGraph < D > ,
251258 node_idx : NodeIndex < usize > ,
@@ -259,39 +266,48 @@ where
259266 let nearest_coords = index
260267 . within ( & coords, tolerance, & squared_euclidean)
261268 . unwrap ( ) ;
262- debug_assert ! (
263- !nearest_coords. is_empty( ) ,
264- "We'll always look up at least ourselves"
265- ) ;
266269
267270 // There's no node close enough to snap to
268- if nearest_coords. len ( ) <= 1 {
271+ if nearest_coords. is_empty ( ) {
269272 return None ;
270273 }
271274
272- let ( mut _distance, mut found_idx) = nearest_coords[ 0 ] ;
273- let found_coord = graph[ * found_idx] . 0 ;
274- if found_coord == graph[ node_idx] . 0 && nearest_coords. len ( ) > 1 {
275- ( _distance, found_idx) = nearest_coords[ 1 ] ;
275+ // Find the closest node that isn't the query node itself
276+ let mut snap_to = None ;
277+ for ( _distance, found_idx) in nearest_coords {
278+ if * found_idx != node_idx {
279+ snap_to = Some ( * found_idx) ;
280+ break ;
281+ }
276282 }
277- let found_idx = * found_idx;
278283
279- index. remove ( & coords, & node_idx) . unwrap ( ) ;
284+ // We found a node to snap to
285+ if let Some ( found_idx) = snap_to {
286+ // Remove the snapped from node from the index, but we have to be careful to only remove it
287+ // if we know the coordinates are actually in the index (duplicate coordinates are filtered
288+ // out ahead of time because they would otherwise cause infinite loops here)
289+ if graph[ found_idx] != graph[ node_idx] {
290+ // Remove the node we're snapping from
291+ index. remove ( & coords, & node_idx) . unwrap ( ) ;
292+ }
280293
281- snap_graph_nodes ( graph, node_idx, found_idx) ;
282- Some ( node_idx)
294+ // Snap the two nodes together, updating the adjacencies
295+ snap_graph_nodes ( graph, node_idx, found_idx) ;
296+ Some ( found_idx)
297+ } else {
298+ None
299+ }
283300}
284301
302+ /// Snap `snap_from` to `snap_to`, and update all of `snap_from`s adjacencies
285303fn snap_graph_nodes < D > (
286304 graph : & mut GeometryGraph < D > ,
287305 snap_from : NodeIndex < usize > ,
288306 snap_to : NodeIndex < usize > ,
289307) where
290308 D : EdgeType ,
291309{
292- if snap_from == snap_to || graph[ snap_from] == graph[ snap_to] {
293- return ;
294- }
310+ debug_assert_ne ! ( snap_from, snap_to) ;
295311
296312 let neighbors: Vec < _ > = graph. neighbors ( snap_from) . collect ( ) ;
297313 let mut neighbors_to_snap = Vec :: new ( ) ;
@@ -358,6 +374,7 @@ mod tests {
358374 use float_cmp:: assert_approx_eq;
359375 use geo:: { LineString , Point } ;
360376 use petgraph:: Undirected ;
377+ use pretty_assertions:: assert_eq;
361378
362379 use super :: * ;
363380 use crate :: io:: { read_tgf_graph, write_tgf_graph} ;
@@ -471,32 +488,102 @@ mod tests {
471488 }
472489
473490 #[ test]
474- fn test_snap_graph_closest_simple ( ) {
475- let tgf = b"0 POINT(0 0)\n 1 POINT(1 0)\n 2 POINT(1.1 0)\n 3 POINT(2 0)\n #\n 0 1\n 1 2\n 2 3" ;
491+ fn test_snap_graph_closest_duplicate ( ) {
492+ let tgf = b"\
493+ 0 POINT(0 0)\n \
494+ 1 POINT(0 0)\n \
495+ #\n \
496+ 0 1\n \
497+ ";
476498 let graph = read_tgf_graph :: < Undirected , _ > ( & tgf[ ..] ) ;
477- assert_eq ! ( graph. node_count( ) , 4 ) ;
478- assert_eq ! ( graph. edge_count( ) , 3 ) ;
499+ assert_eq ! ( graph. node_count( ) , 2 ) ;
500+ assert_eq ! ( graph. edge_count( ) , 1 ) ;
479501
480- let tgf = b"0\t POINT(0 0)\n 1\t POINT(2 0)\n 2\t POINT(1.1 0)\n #\n 2\t 1\n 2\t 0\n " ;
481- let expected_tgf = String :: from_utf8_lossy ( tgf) ;
502+ let expected_tgf = "\
503+ 0\t POINT(0 0)\n \
504+ #\n \
505+ ";
482506
483507 let actual = snap_graph ( graph, SnappingStrategy :: ClosestPoint ( 0.2 ) ) ;
484508
485509 let actual_tgf = get_tgf ( & actual) ;
486- assert_eq ! ( actual_tgf, expected_tgf) ;
510+ assert_eq ! ( expected_tgf, actual_tgf) ;
511+ }
512+
513+ #[ test]
514+ fn test_snap_graph_closest_simple ( ) {
515+ let tgf = b"\
516+ 0 POINT(0 0)\n \
517+ 1 POINT(1 0)\n \
518+ 2 POINT(1.1 0)\n \
519+ 3 POINT(2 0)\n \
520+ 4 POINT(1.1 0)\n \
521+ #\n \
522+ 0 1\n \
523+ 1 2\n \
524+ 2 3\n \
525+ 2 4\n \
526+ ";
527+ let graph = read_tgf_graph :: < Undirected , _ > ( & tgf[ ..] ) ;
528+ assert_eq ! ( graph. node_count( ) , 5 ) ;
529+ assert_eq ! ( graph. edge_count( ) , 4 ) ;
530+
531+ let expected_tgf = "\
532+ 0\t POINT(0 0)\n \
533+ 1\t POINT(2 0)\n \
534+ 2\t POINT(1.1 0)\n \
535+ #\n \
536+ 2\t 1\n \
537+ 2\t 0\n \
538+ ";
539+
540+ let actual = snap_graph ( graph, SnappingStrategy :: ClosestPoint ( 0.2 ) ) ;
541+
542+ let actual_tgf = get_tgf ( & actual) ;
543+ assert_eq ! ( expected_tgf, actual_tgf) ;
487544 }
488545
489546 #[ test]
490547 fn test_snap_graph_closest_complex ( ) {
491- let tgf = b"0\t POINT(-0.1 0)\n 1\t POINT(0 0)\n 2\t POINT(0 0.1)\n 3\t POINT(0 -0.1)\n 4\t POINT(2 0)\n #\n 0\t 1\n 2\t 1\n 3\t 1\n 1\t 4\n " ;
548+ let tgf = b"\
549+ 0\t POINT(-0.1 0)\n \
550+ 1\t POINT(0 0)\n \
551+ 2\t POINT(0 0.1)\n \
552+ 3\t POINT(0 -0.1)\n \
553+ 4\t POINT(2 0)\n \
554+ #\n \
555+ 0\t 1\n \
556+ 2\t 1\n \
557+ 3\t 1\n \
558+ 1\t 4\n \
559+ ";
492560 let graph = read_tgf_graph :: < Undirected , _ > ( & tgf[ ..] ) ;
493561
494- let tgf = b"0\t POINT(2 0)\n 1\t POINT(0 -0.1)\n #\n 1\t 0\n " ;
495- let expected_tgf = String :: from_utf8_lossy ( tgf) ;
562+ let expected_tgf = "\
563+ 0\t POINT(2 0)\n \
564+ 1\t POINT(0 -0.1)\n \
565+ #\n \
566+ 1\t 0\n \
567+ ";
496568
497569 let actual = snap_graph ( graph, SnappingStrategy :: ClosestPoint ( 0.11 ) ) ;
498570
499571 let actual_tgf = get_tgf ( & actual) ;
500- assert_eq ! ( actual_tgf, expected_tgf) ;
572+ assert_eq ! ( expected_tgf, actual_tgf) ;
573+ }
574+
575+ #[ test]
576+ fn test_snap_duplicate_vertices_crash ( ) {
577+ let points = [
578+ Geometry :: Point ( Point :: new ( -0.4999999999999998 , 0.8660254037844387 ) ) ,
579+ Geometry :: Point ( Point :: new ( -0.4999999999999998 , 0.8660254037844387 ) ) ,
580+ ] ;
581+
582+ let snapped: Vec < _ > =
583+ snap_geoms ( points. into_iter ( ) , SnappingStrategy :: ClosestPoint ( 0.0 ) ) . collect ( ) ;
584+
585+ // Snapping can't remove duplicate vertices; it can only move the coordinates of a vertex
586+ // to the coordinates of another vertex.
587+ assert_eq ! ( snapped. len( ) , 2 ) ;
501588 }
502589}
0 commit comments