@@ -13,6 +13,8 @@ use crate::extensions::segments::{FromSegments, HullSegments};
13
13
use std:: cmp:: Ordering ;
14
14
use std:: collections:: BinaryHeap ;
15
15
16
+ /// Store triangle information. Score is used for ranking the priority queue which determines
17
+ /// removal order.
16
18
#[ derive( Debug ) ]
17
19
struct VWScore < T : CoordFloat > {
18
20
score : T ,
@@ -42,6 +44,10 @@ impl<T: CoordFloat> PartialEq for VWScore<T> {
42
44
}
43
45
}
44
46
47
+ /// Area and topology preserving Visvalingam-Whyatt algorithm
48
+ /// adapted from the [geo implementation](https://github.com/georust/geo/blob/e8419735b5986f120ddf1de65ac68c1779c3df30/geo/src/algorithm/simplify_vw.rs)
49
+ ///
50
+ ///
45
51
fn visvalingam_preserve < T > ( orig : & LineString < T > , eps : T , min_len : usize ) -> Vec < Coord < T > >
46
52
where
47
53
T : GeoFloat + RTreeNum ,
56
62
let tree: RTree < CachedEnvelope < _ > > =
57
63
RTree :: bulk_load ( orig. lines ( ) . map ( CachedEnvelope :: new) . collect :: < Vec < _ > > ( ) ) ;
58
64
65
+ // Point adjacency. Tuple at index contains indices into `orig`. Negative values or values
66
+ // greater than or equal to max means no next element. (0, 0) sentinel means deleted element.
59
67
let mut adjacent: Vec < _ > = ( 0 ..orig. 0 . len ( ) )
60
68
. map ( |i| {
61
69
if i == 0 {
66
74
} )
67
75
. collect ( ) ;
68
76
77
+ // Store all triangles in a minimum priority queue, based on signed area.
78
+ //
79
+ // Only triangles of positive score are pushed to the heap.
80
+ // Invalid triangles are *not* removed when the corresponding point is removed; they are
81
+ // invalidated using (0, 0) values in `adjacent` and skipped as necessary.
69
82
let mut pq = orig
70
83
. ord_triangles ( )
71
84
. enumerate ( )
@@ -78,16 +91,25 @@ where
78
91
. filter ( |point| point. score >= T :: zero ( ) )
79
92
. collect :: < BinaryHeap < VWScore < T > > > ( ) ;
80
93
94
+ // Iterate over points while there is an associated triangle with area between 0 and epsilon
81
95
while let Some ( smallest) = pq. pop ( ) {
82
96
if smallest. score > eps {
97
+ // Min-heap guarantees all future points have areas greater than epsilon
98
+ break ;
99
+ }
100
+
101
+ if len <= min_len {
102
+ // Further removal would send us below the minimum length
83
103
break ;
84
104
}
85
105
86
106
let ( left, right) = adjacent[ smallest. current ] ;
107
+ // A point in this triangle has been removed since this `VScore` was created, so skip it
87
108
if left != smallest. left as i32 || right != smallest. right as i32 {
88
109
continue ;
89
110
}
90
111
112
+ // Removal of this point would cause self-intersection, so skip it
91
113
if tree_intersect ( & tree, & smallest, & orig. 0 ) {
92
114
continue ;
93
115
}
@@ -96,23 +118,31 @@ where
96
118
let ( _, rr) = adjacent[ right as usize ] ;
97
119
adjacent[ left as usize ] = ( ll, right) ;
98
120
adjacent[ right as usize ] = ( left, rr) ;
121
+ // Remove the point from the adjacency list
99
122
adjacent[ smallest. current ] = ( 0 , 0 ) ;
123
+ // Update the length of the linestring
100
124
len -= 1 ;
125
+ // The rtree is never updated as self-intersection can never occur with stale segments and
126
+ // if a segment were to intersect with a new segment, then it also intersects with a stale
127
+ // segment
101
128
102
- if len == min_len {
103
- break ;
104
- }
105
-
129
+ // Recompute the areas of adjacent triangles(s) using left and right adjacent points,
130
+ // this may add new triangles to the heap
106
131
recompute_triangles ( orig, & mut pq, ll, left, right, rr, max) ;
107
132
}
108
133
134
+ // Filter out deleted points, returning remaining points
109
135
orig. 0
110
136
. iter ( )
111
137
. zip ( adjacent. iter ( ) )
112
138
. filter_map ( |( tup, adj) | if * adj != ( 0 , 0 ) { Some ( * tup) } else { None } )
113
139
. collect ( )
114
140
}
115
141
142
+ /// Check whether the removal of a candidate point would cause a self-intersection.
143
+ ///
144
+ /// To do this efficiently, and rtree is queried for any existing line segments which fall within
145
+ /// the bounding box of the new line segment created.
116
146
fn tree_intersect < T > (
117
147
tree : & RTree < CachedEnvelope < Line < T > > > ,
118
148
triangle : & VWScore < T > ,
@@ -142,6 +172,7 @@ where
142
172
} )
143
173
}
144
174
175
+ /// Recompute adjacent triangle(s) using left and right adjacent points, pushing to the heap
145
176
fn recompute_triangles < T : CoordFloat > (
146
177
orig : & LineString < T > ,
147
178
pq : & mut BinaryHeap < VWScore < T > > ,
@@ -154,6 +185,7 @@ fn recompute_triangles<T: CoordFloat>(
154
185
let choices = [ ( ll, left, right) , ( left, right, rr) ] ;
155
186
for & ( ai, current_point, bi) in & choices {
156
187
if ai as usize >= max || bi as usize >= max {
188
+ // Out of bounds, i.e. we're at an end point
157
189
continue ;
158
190
}
159
191
@@ -164,6 +196,7 @@ fn recompute_triangles<T: CoordFloat>(
164
196
)
165
197
. signed_area ( ) ;
166
198
199
+ // If removal of a point would cause a reduction in signed area, skip it
167
200
if area < T :: zero ( ) {
168
201
continue ;
169
202
}
@@ -178,7 +211,10 @@ fn recompute_triangles<T: CoordFloat>(
178
211
}
179
212
}
180
213
214
+ /// Simplifies a geometry while preserving its topology and area.
181
215
pub trait SimplifyVW < T , Epsilon = T > {
216
+ /// Returns the simplified geometry using a topology and area preserving variant of the
217
+ /// [Visvalingam-Whyatt](https://doi.org/10.1179/000870493786962263) algorithm.
182
218
fn simplify_vw ( & self , eps : Epsilon , len : usize ) -> Self ;
183
219
}
184
220
@@ -196,6 +232,7 @@ where
196
232
T : GeoFloat + RTreeNum + Send + Sync ,
197
233
{
198
234
fn simplify_vw ( & self , eps : T , len : usize ) -> Self {
235
+ // Get convex hull segments, as their endpoints are invariant under reduction
199
236
let segments = self . hull_segments ( ) ;
200
237
201
238
if len > segments. len ( ) {
0 commit comments