diff --git a/src/crossings.rs b/src/crossings.rs index 9dd38a1..55ca991 100644 --- a/src/crossings.rs +++ b/src/crossings.rs @@ -1,6 +1,11 @@ use anyhow::{Context, Result}; use geo::{ - algorithm::bounding_rect::BoundingRect, Contains, Coord, Geometry, LineString, Point, + algorithm::bounding_rect::BoundingRect, + Contains, + Coord, + Geometry, + LineString, + Point, }; use rstar::{AABB, RTree}; use std::cmp::Ordering; @@ -8,13 +13,17 @@ use std::cmp::Ordering; use crate::{ county::County, geometry::{ - deduplicate_hits, interpolate, segment_intersection, snap_to_segment, + deduplicate_hits, + interpolate, + segment_intersection, + snap_to_segment, }, }; #[derive(Debug, Clone)] pub struct BoundaryHit { pub position: f64, + pub point: Point, } #[derive(Debug, Clone)] @@ -47,7 +56,10 @@ pub struct TrackPoint { /// from the original track segment and the intersection parameter, ensuring /// that the resulting waypoint lies exactly on the original segment. /// -/// Elevation is linearly interpolated between the two original track points. +/// Boundary intersections are clustered before transitions are reconstructed. +/// This is important for polygon datasets such as the German VG250 data, +/// where the same shared Landkreis boundary is represented by both adjacent +/// polygons. pub fn find_crossings( track: &[TrackPoint], tree: &RTree, @@ -76,8 +88,10 @@ pub fn find_crossings( .bounding_rect() .context("Track segment has no bounding box.")?; - let envelope = - AABB::from_corners([bbox.min().x, bbox.min().y], [bbox.max().x, bbox.max().y]); + let envelope = AABB::from_corners( + [bbox.min().x, bbox.min().y], + [bbox.max().x, bbox.max().y], + ); let candidates: Vec = tree .locate_in_envelope_intersecting(envelope) @@ -88,26 +102,31 @@ pub fn find_crossings( continue; } - let hits = collect_segment_hits(start_point, end_point, &candidates); + let hits = + collect_segment_hits(start_point, end_point, &candidates); if hits.is_empty() { continue; } - let transitions = - reconstruct_transitions(start_point, end_point, &candidates, &hits); + let transitions = reconstruct_transitions( + start_point, + end_point, + &candidates, + &hits, + ); for transition in transitions { // Reconstruct the point from the original GPX segment. - // This is the actual snapping step. - let snapped_point = - snap_to_segment(start_point, end_point, transition.position); - - // Interpolate the elevation at the exact crossing position. // - // If one track point has no elevation, use the other one. - // If neither has elevation, fall back to zero because Garmin - // requires an elevation value for waypoint distance handling. + // This is the actual snapping step and guarantees that the + // waypoint lies on the original track segment. + let snapped_point = snap_to_segment( + start_point, + end_point, + transition.position, + ); + let elevation = interpolate_elevation( start.elevation, end.elevation, @@ -133,11 +152,18 @@ fn interpolate_elevation( end: Option, position: f64, ) -> f64 { + let position = position.clamp(0.0, 1.0); + match (start, end) { (Some(start), Some(end)) => { - start + (end - start) * position.clamp(0.0, 1.0) + start + (end - start) * position } - (Some(elevation), None) | (None, Some(elevation)) => elevation, + + (Some(elevation), None) | (None, Some(elevation)) => { + elevation + } + + // Garmin expects an elevation value for the waypoint. (None, None) => 0.0, } } @@ -227,10 +253,13 @@ fn collect_ring_intersections( hits: &mut Vec, ) { for edge in ring.lines() { - if let Some((position, _point)) = + if let Some((position, point)) = segment_intersection(start, end, edge.start, edge.end) { - hits.push(BoundaryHit { position }); + hits.push(BoundaryHit { + position, + point, + }); } } } @@ -246,8 +275,8 @@ fn reconstruct_transitions( for hit in hits { // Sample slightly before and after the intersection. // - // This lets us determine whether the track actually changes - // administrative area rather than merely touching a boundary. + // The sampling distance is relative to the GPX segment, so it + // remains independent of the absolute coordinate values. let before_t = (hit.position - 1e-8).max(0.0); let after_t = (hit.position + 1e-8).min(1.0); @@ -261,7 +290,8 @@ fn reconstruct_transitions( continue; }; - // Ignore boundary touches where the track remains in the same county. + // Ignore boundary touches where the track remains in the same + // administrative area. if from.id == to.id { continue; } @@ -289,6 +319,10 @@ fn county_at_point( fn deduplicate_crossings( mut crossings: Vec, ) -> Vec { + if crossings.len() <= 1 { + return crossings; + } + crossings.sort_by(|a, b| { a.segment_index .cmp(&b.segment_index) @@ -299,14 +333,41 @@ fn deduplicate_crossings( }) }); - let mut result = Vec::new(); + // Same tolerance used for boundary-hit clustering. + const POSITION_EPSILON: f64 = 1e-7; + + let mut result: Vec = Vec::new(); for crossing in crossings { - let duplicate = result.iter().any(|existing: &Crossing| { - existing.segment_index == crossing.segment_index - && existing.from.id == crossing.from.id + let duplicate = result.iter().any(|existing| { + if existing.segment_index != crossing.segment_index { + return false; + } + + if (existing.position - crossing.position).abs() + > POSITION_EPSILON + { + return false; + } + + // Same transition reported twice. + if existing.from.id == crossing.from.id && existing.to.id == crossing.to.id - && same_point(existing.point, crossing.point) + { + return true; + } + + // The most important VG250 case: + // + // A -> B + // B -> A + // + // at effectively the same physical location. + // + // This can arise when both sides of a shared boundary produce + // an intersection independently. + existing.from.id == crossing.to.id + && existing.to.id == crossing.from.id }); if !duplicate { @@ -316,10 +377,3 @@ fn deduplicate_crossings( result } - -fn same_point(a: Point, b: Point) -> bool { - let dx = a.x() - b.x(); - let dy = a.y() - b.y(); - - dx * dx + dy * dy < 1e-20 -} diff --git a/src/geometry.rs b/src/geometry.rs index 4c21143..8e28c7e 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -2,17 +2,41 @@ use geo::{Coord, Point}; use crate::crossings::BoundaryHit; +/// Calculate the 2D cross product of two vectors. pub fn cross(a: Coord, b: Coord) -> f64 { a.x * b.y - a.y * b.x } -pub fn interpolate(start: Point, end: Point, t: f64) -> Point { +/// Interpolate a point on a segment. +/// +/// `position = 0.0` returns `start`. +/// `position = 1.0` returns `end`. +pub fn interpolate( + start: Point, + end: Point, + position: f64, +) -> Point { Point::new( - start.x() + t * (end.x() - start.x()), - start.y() + t * (end.y() - start.y()), + start.x() + position * (end.x() - start.x()), + start.y() + position * (end.y() - start.y()), ) } +/// Calculate the squared coordinate distance between two points. +/// +/// This intentionally operates in the coordinate system of the input data. +/// For the usual WGS84 GPX/GeoJSON workflow this is a squared degree +/// distance and is only used for a very small numerical tolerance. +pub fn squared_distance(a: Point, b: Point) -> f64 { + let dx = a.x() - b.x(); + let dy = a.y() - b.y(); + + dx * dx + dy * dy +} + +/// Intersect two line segments. +/// +/// Returns the position on segment `p -> p2` and the intersection point. pub fn segment_intersection( p: Coord, p2: Coord, @@ -45,18 +69,57 @@ pub fn segment_intersection( const EPSILON: f64 = 1e-10; - if !(-EPSILON..=1.0 + EPSILON).contains(&t) || !(-EPSILON..=1.0 + EPSILON).contains(&u) { + if !(-EPSILON..=1.0 + EPSILON).contains(&t) + || !(-EPSILON..=1.0 + EPSILON).contains(&u) + { return None; } let t = t.clamp(0.0, 1.0); - let point = Point::new(p.x + t * r.x, p.y + t * r.y); + let point = Point::new( + p.x + t * r.x, + p.y + t * r.y, + ); Some((t, point)) } +/// Deduplicate and cluster boundary intersections. +/// +/// A shared Landkreis boundary is normally present in both adjacent +/// polygons. Consequently, the same physical crossing can be reported +/// several times with slightly different floating-point positions. +/// +/// We cluster hits using both: +/// +/// 1. their position along the GPX segment, and +/// 2. their actual coordinate distance. +/// +/// The tolerances are deliberately small enough not to merge ordinary +/// separate crossings, while being large enough for the numerical noise +/// introduced by polygon conversion and floating-point intersection. +/// +/// The function preserves the first hit in each cluster and uses the +/// average position for the cluster. The actual crossing point is later +/// reconstructed from the original GPX segment. pub fn deduplicate_hits(mut hits: Vec) -> Vec { + if hits.len() <= 1 { + return hits; + } + + // About one metre at German latitudes when coordinates are WGS84. + // + // This is deliberately conservative. We mainly want to collapse + // duplicate representations of the same shared polygon boundary. + const COORD_EPSILON: f64 = 1.5e-5; + const COORD_EPSILON_SQUARED: f64 = + COORD_EPSILON * COORD_EPSILON; + + // Position epsilon protects against tiny differences in the + // intersection calculation itself. + const POSITION_EPSILON: f64 = 1e-7; + hits.sort_by(|a, b| { a.position .partial_cmp(&b.position) @@ -65,23 +128,76 @@ pub fn deduplicate_hits(mut hits: Vec) -> Vec { let mut result = Vec::new(); - for hit in hits { - let duplicate = result - .iter() - .any(|existing: &BoundaryHit| (existing.position - hit.position).abs() < 1e-9); + let mut cluster: Vec = Vec::new(); - if !duplicate { - result.push(hit); + for hit in hits { + if cluster.is_empty() { + cluster.push(hit); + continue; } + + let representative = cluster + .last() + .expect("cluster cannot be empty"); + + let position_close = + (representative.position - hit.position).abs() + <= POSITION_EPSILON; + + let point_close = + squared_distance(representative.point, hit.point) + <= COORD_EPSILON_SQUARED; + + if position_close && point_close { + cluster.push(hit); + } else { + result.push(merge_hit_cluster(&cluster)); + cluster.clear(); + cluster.push(hit); + } + } + + if !cluster.is_empty() { + result.push(merge_hit_cluster(&cluster)); } result } -/// Recalculate the point directly from the original track segment. +/// Merge several numerical representations of the same boundary hit. +fn merge_hit_cluster(cluster: &[BoundaryHit]) -> BoundaryHit { + debug_assert!(!cluster.is_empty()); + + let position = + cluster.iter().map(|hit| hit.position).sum::() + / cluster.len() as f64; + + let x = + cluster.iter().map(|hit| hit.point.x()).sum::() + / cluster.len() as f64; + + let y = + cluster.iter().map(|hit| hit.point.y()).sum::() + / cluster.len() as f64; + + BoundaryHit { + position, + point: Point::new(x, y), + } +} + +/// Recalculate the point directly from the original GPX segment. /// -/// This guarantees that the returned point lies on the segment rather than -/// relying on the coordinates calculated during the intersection operation. -pub fn snap_to_segment(start: Point, end: Point, position: f64) -> Point { - interpolate(start, end, position.clamp(0.0, 1.0)) +/// This guarantees that the returned point lies on the original segment +/// rather than relying on the coordinates calculated during intersection. +pub fn snap_to_segment( + start: Point, + end: Point, + position: f64, +) -> Point { + interpolate( + start, + end, + position.clamp(0.0, 1.0), + ) }