opentracker

opentracker/crates/otserver/src/polyline.rs
polyline.rs
⎇
Raw
1//! Trail decimation and encoding.
2//!
3//! A day of walking at one fix per 25 s is ~3500 points. Sent as JSON floats that
4//! is ~100 kB; as a Ramer–Douglas–Peucker-simplified encoded polyline it is a few
5//! kilobytes, and the browser has proportionally less to draw. Both halves matter:
6//! simplification removes points that would render on top of each other anyway,
7//! and the encoding removes the per-number punctuation.
8//!
9//! Coordinates arrive as degrees × 1e7 (the protocol and database representation)
10//! and are encoded at 1e5, which is the polyline format's fixed precision — about
11//! 1.1 m, well under the accuracy of any consumer GPS fix.
12
13/// Points are `(lat_e7, lon_e7)`.
14type Pt = (i64, i64);
15
16/// Simplify to at most `max` points using Ramer–Douglas–Peucker.
17///
18/// The tolerance is searched rather than fixed, because the useful question is
19/// "how much detail fits in the budget?" and the answer depends on how far the
20/// person actually travelled. A fixed ε either keeps too much on a road trip or
21/// destroys a walk around a park.
22pub fn simplify(points: &[Pt], max: usize) -> Vec<Pt> {
23 if points.len() <= max {
24 return points.to_vec();
25 }
26
27 // Ramer–Douglas–Peucker is O(n log n) on a typical track but O(n²) in the
28 // worst case — a stationary phone reporting a heartbeat all week is exactly
29 // that case — and the tolerance search below runs it several times. Uniformly
30 // subsampling first bounds the work at a size where that no longer matters.
31 // At 8× the budget the subsampling itself removes no visible detail: anything
32 // it drops would have been RDP's next victim anyway.
33 let owned;
34 let points = if points.len() > max.saturating_mul(PRE_DECIMATE_FACTOR) {
35 owned = uniform(points, max * PRE_DECIMATE_FACTOR);
36 owned.as_slice()
37 } else {
38 points
39 };
40
41 // ~15 m expressed in units of 1e-7 degrees.
42 let base: f64 = 15.0 * 1e7 / 111_320.0;
43
44 // Grow until something fits, to bracket the answer. Plain doubling would stop
45 // at the *first* tolerance that fits, and on a track with real large-scale
46 // shape that first hit overshoots: a meander can collapse to its two
47 // endpoints, which is a straight line where there was a walk. So this phase
48 // only establishes an upper bound.
49 let mut lo = 0.0f64; // known to yield more than `max`
50 let mut hi = base;
51 let mut best: Option<Vec<Pt>> = None;
52 for _ in 0..24 {
53 let out = rdp(points, hi);
54 if out.len() <= max {
55 best = Some(out);
56 break;
57 }
58 lo = hi;
59 hi *= 4.0;
60 }
61
62 // Then bisect, keeping the largest result that still fits: the most detail the
63 // budget allows, rather than the first thing under it.
64 if let Some(mut candidate) = best {
65 for _ in 0..8 {
66 let mid = (lo + hi) / 2.0;
67 let out = rdp(points, mid);
68 if out.len() <= max {
69 if out.len() >= candidate.len() {
70 candidate = out;
71 }
72 hi = mid;
73 } else {
74 lo = mid;
75 }
76 }
77 return candidate;
78 }
79
80 // No tolerance in 24 quadruplings fit the budget. Unreachable for real data;
81 // uniform sampling keeps the shape far better than an infinite tolerance.
82 uniform(points, max)
83}
84
85/// How much larger than the budget the RDP input may be.
86const PRE_DECIMATE_FACTOR: usize = 8;
87
88/// Evenly spaced subsample of at most `max` points, always keeping the last one so
89/// a trail still ends where the person is.
90fn uniform(points: &[Pt], max: usize) -> Vec<Pt> {
91 if points.len() <= max {
92 return points.to_vec();
93 }
94 let step = points.len().div_ceil(max);
95 let mut out: Vec<Pt> = points.iter().copied().step_by(step).collect();
96 if out.last() != points.last() {
97 out.push(*points.last().expect("non-empty"));
98 }
99 out
100}
101
102fn rdp(points: &[Pt], epsilon: f64) -> Vec<Pt> {
103 if points.len() < 3 {
104 return points.to_vec();
105 }
106 let mut keep = vec![false; points.len()];
107 keep[0] = true;
108 keep[points.len() - 1] = true;
109 // Iterative rather than recursive: a 100k-point input would otherwise be able
110 // to blow the stack, and this runs on attacker-influenced data volumes.
111 let mut stack = vec![(0usize, points.len() - 1)];
112 while let Some((start, end)) = stack.pop() {
113 if end <= start + 1 {
114 continue;
115 }
116 let mut worst = 0.0;
117 let mut worst_i = start;
118 for i in (start + 1)..end {
119 let d = perpendicular_distance(points[i], points[start], points[end]);
120 if d > worst {
121 worst = d;
122 worst_i = i;
123 }
124 }
125 if worst > epsilon {
126 keep[worst_i] = true;
127 stack.push((start, worst_i));
128 stack.push((worst_i, end));
129 }
130 }
131 points
132 .iter()
133 .zip(keep)
134 .filter_map(|(p, k)| k.then_some(*p))
135 .collect()
136}
137
138fn perpendicular_distance(p: Pt, a: Pt, b: Pt) -> f64 {
139 let (px, py) = (p.1 as f64, p.0 as f64);
140 let (ax, ay) = (a.1 as f64, a.0 as f64);
141 let (bx, by) = (b.1 as f64, b.0 as f64);
142 let dx = bx - ax;
143 let dy = by - ay;
144 let len_sq = dx * dx + dy * dy;
145 if len_sq == 0.0 {
146 return ((px - ax).powi(2) + (py - ay).powi(2)).sqrt();
147 }
148 ((dx * (ay - py) - (ax - px) * dy).abs()) / len_sq.sqrt()
149}
150
151/// Google's encoded polyline algorithm at 1e5 precision.
152pub fn encode(points: &[Pt]) -> String {
153 let mut out = String::with_capacity(points.len() * 6);
154 let mut prev_lat = 0i64;
155 let mut prev_lon = 0i64;
156 for &(lat_e7, lon_e7) in points {
157 // 1e7 -> 1e5, rounding rather than truncating so error stays centred.
158 let lat = div_round(lat_e7, 100);
159 let lon = div_round(lon_e7, 100);
160 encode_value(lat - prev_lat, &mut out);
161 encode_value(lon - prev_lon, &mut out);
162 prev_lat = lat;
163 prev_lon = lon;
164 }
165 out
166}
167
168fn div_round(v: i64, d: i64) -> i64 {
169 if v >= 0 {
170 (v + d / 2) / d
171 } else {
172 -((-v + d / 2) / d)
173 }
174}
175
176fn encode_value(value: i64, out: &mut String) {
177 let mut v = if value < 0 { !(value << 1) } else { value << 1 };
178 while v >= 0x20 {
179 out.push(char::from(((0x20 | (v & 0x1f)) + 63) as u8));
180 v >>= 5;
181 }
182 out.push(char::from((v + 63) as u8));
183}
184
185/// Decode. Used by the tests to prove the encoder is reversible, and kept public
186/// because a Rust consumer of `/api/users/:id/track` needs it.
187#[cfg_attr(not(test), allow(dead_code))]
188pub fn decode(s: &str) -> Vec<Pt> {
189 let bytes = s.as_bytes();
190 let mut out = Vec::new();
191 let mut i = 0;
192 let mut lat = 0i64;
193 let mut lon = 0i64;
194 while i < bytes.len() {
195 let Some(dlat) = decode_value(bytes, &mut i) else {
196 break;
197 };
198 let Some(dlon) = decode_value(bytes, &mut i) else {
199 break;
200 };
201 lat += dlat;
202 lon += dlon;
203 out.push((lat * 100, lon * 100));
204 }
205 out
206}
207
208#[cfg_attr(not(test), allow(dead_code))]
209fn decode_value(bytes: &[u8], i: &mut usize) -> Option<i64> {
210 let mut shift = 0;
211 let mut result = 0i64;
212 loop {
213 let b = *bytes.get(*i)? as i64 - 63;
214 *i += 1;
215 result |= (b & 0x1f) << shift;
216 shift += 5;
217 if b < 0x20 {
218 break;
219 }
220 if shift > 60 {
221 return None; // malformed; refuse to shift forever
222 }
223 }
224 Some(if result & 1 != 0 {
225 !(result >> 1)
226 } else {
227 result >> 1
228 })
229}
230
231#[cfg(test)]
232mod tests {
233 use super::*;
234
235 #[test]
236 fn a_known_vector_matches_the_reference_implementation() {
237 // The example from Google's polyline documentation, in 1e7 units.
238 // (38.5, -120.2), (40.7, -120.95), (43.252, -126.453)
239 let points = vec![
240 (385_000_000, -1_202_000_000),
241 (407_000_000, -1_209_500_000),
242 (432_520_000, -1_264_530_000),
243 ];
244 assert_eq!(encode(&points), "_p~iF~ps|U_ulLnnqC_mqNvxq`@");
245 }
246
247 #[test]
248 fn encode_decode_round_trips_within_the_formats_precision() {
249 let points = vec![
250 (525_200_080, 134_050_000),
251 (525_210_000, 134_060_000),
252 (-338_688_000, -1_754_500_000),
253 ];
254 let back = decode(&encode(&points));
255 assert_eq!(back.len(), points.len());
256 for (a, b) in points.iter().zip(&back) {
257 // 1e5 precision means the last two 1e7 digits are lost: ≤ 50 units,
258 // about 0.5 cm. Well inside any GPS accuracy.
259 assert!((a.0 - b.0).abs() <= 50, "lat drifted: {a:?} vs {b:?}");
260 assert!((a.1 - b.1).abs() <= 50, "lon drifted: {a:?} vs {b:?}");
261 }
262 }
263
264 #[test]
265 fn an_empty_track_encodes_to_an_empty_string() {
266 assert_eq!(encode(&[]), "");
267 assert_eq!(decode(""), Vec::<Pt>::new());
268 }
269
270 #[test]
271 fn decoding_garbage_does_not_panic_or_hang() {
272 for s in ["~", "?????", "\u{1}\u{2}\u{3}", &"~".repeat(1000)] {
273 let _ = decode(s);
274 }
275 }
276
277 #[test]
278 fn simplify_keeps_short_tracks_untouched() {
279 let points: Vec<Pt> = (0..10).map(|i| (i * 1000, i * 1000)).collect();
280 assert_eq!(simplify(&points, 2000), points);
281 }
282
283 #[test]
284 fn simplify_respects_the_budget() {
285 // A 20k-point meander: a long sinusoidal path with per-fix jitter on top.
286 // Large-scale shape matters here — a straight line with noise legitimately
287 // simplifies to two points, so it would not test anything.
288 let points: Vec<Pt> = (0..20_000i64)
289 .map(|i| {
290 let jitter = if i % 2 == 0 { 300 } else { -300 };
291 let wave = (5_000_000.0 * ((i as f64) / 300.0).sin()) as i64;
292 (525_200_000 + i * 40 + jitter, 134_050_000 + wave - jitter)
293 })
294 .collect();
295 let out = simplify(&points, 500);
296 assert!(out.len() <= 500, "budget exceeded: {} points", out.len());
297 // And it must actually *use* the budget. Returning two endpoints would
298 // satisfy the limit while turning a wander into a straight line.
299 assert!(
300 out.len() > 250,
301 "budget under-used: only {} points",
302 out.len()
303 );
304 }
305
306 #[test]
307 fn simplify_keeps_the_endpoints() {
308 let points: Vec<Pt> = (0..5_000)
309 .map(|i| (525_200_000 + i * 100, 134_050_000))
310 .collect();
311 let out = simplify(&points, 100);
312 assert_eq!(
313 out.first(),
314 points.first(),
315 "the start of a trail must survive"
316 );
317 assert_eq!(out.last(), points.last(), "the end of a trail must survive");
318 }
319
320 #[test]
321 fn simplify_drops_collinear_interior_points() {
322 // A dead-straight line: everything between the ends is redundant.
323 let points: Vec<Pt> = (0..3_000)
324 .map(|i| (525_200_000 + i * 1_000, 134_050_000))
325 .collect();
326 let out = simplify(&points, 2_000);
327 assert!(
328 out.len() < 50,
329 "a straight line should collapse, got {}",
330 out.len()
331 );
332 }
333
334 #[test]
335 fn simplify_survives_a_track_that_never_moves() {
336 // Every point identical: len_sq == 0 in the distance function.
337 let points: Vec<Pt> = std::iter::repeat_n((525_200_000, 134_050_000), 5_000).collect();
338 let out = simplify(&points, 100);
339 assert!(out.len() <= 100);
340 }
341
342 #[test]
343 fn negative_coordinates_encode_correctly() {
344 // The zig-zag encoding of negatives is the classic place to get an
345 // off-by-one, and half the planet is at a negative longitude.
346 let points = vec![(-338_688_000, -1_754_500_000)];
347 let back = decode(&encode(&points));
348 assert_eq!(back.len(), 1);
349 assert!(back[0].0 < 0 && back[0].1 < 0);
350 }
351}
352