diff --git a/apps/gipy/gpconv/src/main.rs b/apps/gipy/gpconv/src/main.rs index a19f0530d..4f2d4a079 100644 --- a/apps/gipy/gpconv/src/main.rs +++ b/apps/gipy/gpconv/src/main.rs @@ -1,4 +1,5 @@ use itertools::Itertools; +use std::collections::HashMap; use std::fs::File; use std::io::{BufReader, Write}; use std::path::Path; @@ -6,6 +7,14 @@ use std::path::Path; use gpx::read; use gpx::{Gpx, Track, TrackSegment}; +fn squared_distance_between(p1: &(f64, f64), p2: &(f64, f64)) -> f64 { + let (x1, y1) = *p1; + let (x2, y2) = *p2; + let dx = x2 - x1; + let dy = y2 - y1; + dx * dx + dy * dy +} + fn points(filename: &str) -> impl Iterator { // This XML file actually exists — try it for yourself! let file = File::open(filename).unwrap(); @@ -22,6 +31,7 @@ fn points(filename: &str) -> impl Iterator { .into_iter() .flat_map(|segment| segment.linestring().points().collect::>()) .map(|point| (point.x(), point.y())) + .dedup() } // returns distance from point p to line passing through points p1 and p2 @@ -32,33 +42,193 @@ fn distance_to_line(p: &(f64, f64), p1: &(f64, f64), p2: &(f64, f64)) -> f64 { let (x2, y2) = *p2; let dx = x2 - x1; let dy = y2 - y1; + //TODO: remove this division by computing fake distances (dx * (y1 - y0) - dy * (x1 - x0)).abs() / (dx * dx + dy * dy).sqrt() } +fn acceptable_angles(p1: &(f64, f64), p2: &(f64, f64), epsilon: f64) -> (f64, f64) { + // first, convert p2's coordinates for p1 as origin + let (x1, y1) = *p1; + let (x2, y2) = *p2; + let (x, y) = (x2 - x1, y2 - y1); + // rotate so that (p1, p2) ends on x axis + let theta = y.atan2(x); + let rx = x * theta.cos() - y * theta.sin(); + let ry = x * theta.sin() + y * theta.cos(); + assert!(ry.abs() <= std::f64::EPSILON); + + // now imagine a line at an angle alpha. + // we want the distance d from (rx, 0) to our line + // we have sin(alpha) = d / rx + // limiting d to epsilon, we solve + // sin(alpha) = e / rx + // and get + // alpha = arcsin(e/rx) + let alpha = (epsilon / rx).asin(); + + // now we just need to rotate back + let a1 = theta + alpha.abs(); + let a2 = theta - alpha.abs(); + assert!(a1 >= a2); + (a1, a2) +} + +// this is like ramer douglas peucker algorithm +// except that we advance from the start without knowing the end. +// each point we meet constrains the chosen segment's angle +// a bit more. +// +fn simplify(mut points: &[(f64, f64)]) -> Vec<(f64, f64)> { + let mut remaining_points = Vec::new(); + while !points.is_empty() { + let (sx, sy) = points.first().unwrap(); + let i = match points + .iter() + .enumerate() + .map(|(i, (x, y))| todo!("compute angles")) + .try_fold( + (0.0f64, std::f64::consts::FRAC_2_PI), + |(amin, amax), (i, (amin2, amax2))| -> Result<(f64, f64), usize> { + let new_amax = amax.min(amax2); + let new_amin = amin.max(amin2); + if new_amin >= new_amax { + Err(i) + } else { + Ok((new_amin, new_amax)) + } + }, + ) { + Err(i) => i, + Ok(_) => points.len(), + }; + remaining_points.push(points.first().cloned().unwrap()); + points = &points[i..]; + } + remaining_points +} + +fn extract_prog_dyn_solution( + points: &[(f64, f64)], + start: usize, + end: usize, + cache: &HashMap<(usize, usize), (Option, usize)>, +) -> Vec<(f64, f64)> { + if let Some(choice) = cache.get(&(start, end)).unwrap().0 { + let mut v1 = extract_prog_dyn_solution(points, start, choice + 1, cache); + let mut v2 = extract_prog_dyn_solution(points, choice, end, cache); + v1.pop(); + v1.append(&mut v2); + v1 + } else { + vec![points[start], points[end - 1]] + } +} + +fn simplify_prog_dyn( + points: &[(f64, f64)], + start: usize, + end: usize, + epsilon: f64, + cache: &mut HashMap<(usize, usize), (Option, usize)>, +) -> usize { + if let Some(val) = cache.get(&(start, end)) { + val.1 + } else { + let res = if end - start <= 2 { + assert_eq!(end - start, 2); + (None, end - start) + } else { + let first_point = &points[start]; + let last_point = &points[end - 1]; + if points[(start + 1)..end] + .iter() + .map(|p| distance_to_line(p, first_point, last_point)) + .all(|d| d <= epsilon) + { + (None, 2) + } else { + // now we test all possible cutting points + ((start + 1)..(end - 1)) //TODO: take middle min + .map(|i| { + let v1 = simplify_prog_dyn(points, start, i + 1, epsilon, cache); + let v2 = simplify_prog_dyn(points, i, end, epsilon, cache); + (Some(i), v1 + v2 - 1) + }) + .min_by_key(|(_, v)| *v) + .unwrap() + } + }; + cache.insert((start, end), res); + res.1 + } +} + fn rdp(points: &[(f64, f64)], epsilon: f64) -> Vec<(f64, f64)> { if points.len() <= 2 { points.iter().copied().collect() } else { - let (index_farthest, farthest_distance) = points - .iter() - .map(|p| distance_to_line(p, points.first().unwrap(), points.last().unwrap())) - .enumerate() - .max_by(|(_, d1), (_, d2)| d1.partial_cmp(d2).unwrap()) - .unwrap(); - if farthest_distance <= epsilon { - vec![ - points.first().copied().unwrap(), - points.last().copied().unwrap(), - ] - } else { - let (start, end) = points.split_at(index_farthest); + if points.first().unwrap() == points.last().unwrap() { + let first = points.first().unwrap(); + let index_farthest = points + .iter() + .enumerate() + .skip(1) + .max_by(|(_, p1), (_, p2)| { + squared_distance_between(first, p1) + .partial_cmp(&squared_distance_between(first, p2)) + .unwrap() + }) + .map(|(i, _)| i) + .unwrap(); + + let start = &points[..(index_farthest + 1)]; + let end = &points[index_farthest..]; let mut res = rdp(start, epsilon); + res.pop(); res.append(&mut rdp(end, epsilon)); res + } else { + let (index_farthest, farthest_distance) = points + .iter() + .map(|p| distance_to_line(p, points.first().unwrap(), points.last().unwrap())) + .enumerate() + .max_by(|(_, d1), (_, d2)| { + if d1.is_nan() { + std::cmp::Ordering::Greater + } else { + if d2.is_nan() { + std::cmp::Ordering::Less + } else { + d1.partial_cmp(d2).unwrap() + } + } + }) + .unwrap(); + if farthest_distance <= epsilon { + vec![ + points.first().copied().unwrap(), + points.last().copied().unwrap(), + ] + } else { + let start = &points[..(index_farthest + 1)]; + let end = &points[index_farthest..]; + let mut res = rdp(start, epsilon); + res.pop(); + res.append(&mut rdp(end, epsilon)); + res + } } } } +fn simplify_path(points: &[(f64, f64)], epsilon: f64) -> Vec<(f64, f64)> { + if points.len() <= 600 { + optimal_simplification(points, epsilon) + } else { + hybrid_simplification(points, epsilon) + } +} + fn convert_coordinates(points: &[(f64, f64)]) -> (f64, f64, Vec<(i32, i32)>) { let xmin = points .iter() @@ -157,53 +327,116 @@ fn save_json>(path: P, points: &[(f64, f64)]) -> std::io::Result< Ok(()) } +fn optimal_simplification(points: &[(f64, f64)], epsilon: f64) -> Vec<(f64, f64)> { + let mut cache = HashMap::new(); + simplify_prog_dyn(&points, 0, points.len(), epsilon, &mut cache); + extract_prog_dyn_solution(&points, 0, points.len(), &cache) +} + +fn hybrid_simplification(points: &[(f64, f64)], epsilon: f64) -> Vec<(f64, f64)> { + if points.len() <= 300 { + optimal_simplification(points, epsilon) + } else { + if points.first().unwrap() == points.last().unwrap() { + let first = points.first().unwrap(); + let index_farthest = points + .iter() + .enumerate() + .skip(1) + .max_by(|(_, p1), (_, p2)| { + squared_distance_between(first, p1) + .partial_cmp(&squared_distance_between(first, p2)) + .unwrap() + }) + .map(|(i, _)| i) + .unwrap(); + + let start = &points[..(index_farthest + 1)]; + let end = &points[index_farthest..]; + let mut res = hybrid_simplification(start, epsilon); + res.pop(); + res.append(&mut hybrid_simplification(end, epsilon)); + res + } else { + let (index_farthest, farthest_distance) = points + .iter() + .map(|p| distance_to_line(p, points.first().unwrap(), points.last().unwrap())) + .enumerate() + .max_by(|(_, d1), (_, d2)| { + if d1.is_nan() { + std::cmp::Ordering::Greater + } else { + if d2.is_nan() { + std::cmp::Ordering::Less + } else { + d1.partial_cmp(d2).unwrap() + } + } + }) + .unwrap(); + if farthest_distance <= epsilon { + vec![ + points.first().copied().unwrap(), + points.last().copied().unwrap(), + ] + } else { + let start = &points[..(index_farthest + 1)]; + let end = &points[index_farthest..]; + let mut res = hybrid_simplification(start, epsilon); + res.pop(); + res.append(&mut hybrid_simplification(end, epsilon)); + res + } + } + } +} + fn main() { let input_file = std::env::args().nth(1).unwrap_or("m.gpx".to_string()); eprintln!("input is {}", input_file); let p = points(&input_file).collect::>(); - let rp = rdp(&p, 0.001); - // let rp = rdp(&p, 0.0001); + eprintln!("initialy we have {} points", p.len()); + //eprintln!("opt is {}", optimal_simplification(&p, 0.0005).len()); + let start = std::time::Instant::now(); + let rp = hybrid_simplification(&p, 0.0005); + eprintln!("hybrid took {:?}", start.elapsed()); + eprintln!("we now have {} points", rp.len()); + let start = std::time::Instant::now(); + eprintln!("rdp would have had {}", rdp(&p, 0.0005).len()); + eprintln!("rdp took {:?}", start.elapsed()); + // let rp = rdp(&p, 0.001); save_coordinates("test.gpc", &rp).unwrap(); - return; - eprintln!("we go from {} to {}", p.len(), rp.len()); - //TODO: assert we don't wrap around the globe - let (xmin, ymin, p) = convert_coordinates(&rp); - // let diffs = compress_coordinates(&p); + let (xmin, xmax) = rp + .iter() + .map(|&(x, _)| x) + .minmax_by(|a, b| a.partial_cmp(b).unwrap()) + .into_option() + .unwrap(); - // save_coordinates("test.gpc", xmin, ymin, &p).unwrap(); + let (ymin, ymax) = rp + .iter() + .map(|&(_, y)| y) + .minmax_by(|a, b| a.partial_cmp(b).unwrap()) + .into_option() + .unwrap(); - // // compress_coordinates(&p); - // let (xmin, xmax) = p - // .iter() - // .map(|&(x, _)| x) - // .minmax_by(|a, b| a.partial_cmp(b).unwrap()) - // .into_option() - // .unwrap(); - - // let (ymin, ymax) = p - // .iter() - // .map(|&(_, y)| y) - // .minmax_by(|a, b| a.partial_cmp(b).unwrap()) - // .into_option() - // .unwrap(); - - // println!( - // "", - // xmin, - // ymin, - // xmax - xmin, - // ymax - ymin - // ); - // print!( - // "", - // xmin, - // ymin, - // xmax - xmin, - // ymax - ymin - // ); - // print!(""); - // println!(""); + println!( + "", + xmin, + ymin, + xmax - xmin, + ymax - ymin + ); + print!( + "", + xmin, + ymin, + xmax - xmin, + ymax - ymin + ); + print!(""); + println!(""); }