better deduplication

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<f64>,
 }
 
 #[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<County>,
@@ -76,8 +88,10 @@ pub fn find_crossings(
             .bounding_rect()
             .context("Track segment has no bounding box.")?;
 
-        let envelope =
+        let envelope = AABB::from_corners(
-            AABB::from_corners([bbox.min().x, bbox.min().y], [bbox.max().x, bbox.max().y]);
+            [bbox.min().x, bbox.min().y],
+            [bbox.max().x, bbox.max().y],
+        );
 
         let candidates: Vec<County> = 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 =
+        let transitions = reconstruct_transitions(
-            reconstruct_transitions(start_point, end_point, &candidates, &hits);
+            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.
+            // This is the actual snapping step and guarantees that the
-            // If neither has elevation, fall back to zero because Garmin
+            // waypoint lies on the original track segment.
-            // requires an elevation value for waypoint distance handling.
+            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<f64>,
     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<BoundaryHit>,
 ) {
     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
+        // The sampling distance is relative to the GPX segment, so it
-        // administrative area rather than merely touching a boundary.
+        // 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<Crossing>,
 ) -> Vec<Crossing> {
+    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<Crossing> = Vec::new();
 
     for crossing in crossings {
-        let duplicate = result.iter().any(|existing: &Crossing| {
+        let duplicate = result.iter().any(|existing| {
-            existing.segment_index == crossing.segment_index
+            if existing.segment_index != crossing.segment_index {
-                && existing.from.id == crossing.from.id
+                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<f64>, b: Point<f64>) -> bool {
-    let dx = a.x() - b.x();
-    let dy = a.y() - b.y();
-
-    dx * dx + dy * dy < 1e-20
-}

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<f64>, b: Coord<f64>) -> f64 {
     a.x * b.y - a.y * b.x
 }
 
-pub fn interpolate(start: Point<f64>, end: Point<f64>, t: f64) -> Point<f64> {
+/// Interpolate a point on a segment.
+///
+/// `position = 0.0` returns `start`.
+/// `position = 1.0` returns `end`.
+pub fn interpolate(
+    start: Point<f64>,
+    end: Point<f64>,
+    position: f64,
+) -> Point<f64> {
     Point::new(
-        start.x() + t * (end.x() - start.x()),
+        start.x() + position * (end.x() - start.x()),
-        start.y() + t * (end.y() - start.y()),
+        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<f64>, b: Point<f64>) -> 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<f64>,
     p2: Coord<f64>,
@@ -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<BoundaryHit>) -> Vec<BoundaryHit> {
+    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<BoundaryHit>) -> Vec<BoundaryHit> {
 
     let mut result = Vec::new();
 
-    for hit in hits {
+    let mut cluster: Vec<BoundaryHit> = Vec::new();
-        let duplicate = result
-            .iter()
-            .any(|existing: &BoundaryHit| (existing.position - hit.position).abs() < 1e-9);
 
-        if !duplicate {
+    for hit in hits {
-            result.push(hit);
+        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::<f64>()
+            / cluster.len() as f64;
+
+    let x =
+        cluster.iter().map(|hit| hit.point.x()).sum::<f64>()
+            / cluster.len() as f64;
+
+    let y =
+        cluster.iter().map(|hit| hit.point.y()).sum::<f64>()
+            / 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
+/// This guarantees that the returned point lies on the original segment
-/// relying on the coordinates calculated during the intersection operation.
+/// rather than relying on the coordinates calculated during intersection.
-pub fn snap_to_segment(start: Point<f64>, end: Point<f64>, position: f64) -> Point<f64> {
+pub fn snap_to_segment(
-    interpolate(start, end, position.clamp(0.0, 1.0))
+    start: Point<f64>,
+    end: Point<f64>,
+    position: f64,
+) -> Point<f64> {
+    interpolate(
+        start,
+        end,
+        position.clamp(0.0, 1.0),
+    )
 }