rustial_engine/

simplify.rs

// ---------------------------------------------------------------------------
//! # Douglas-Peucker line simplification for vector LOD
//!
//! This module provides the [Douglas-Peucker][dp] algorithm for reducing
//! the vertex count of polylines and polygon rings while preserving their
//! visual shape within a caller-specified tolerance.
//!
//! ## Intended use
//!
//! The engine uses this during zoom-dependent level-of-detail (LOD)
//! reduction in [`VectorLayer`](crate::layers::VectorLayer) and
//! exposes it through the public API so that host applications can
//! pre-simplify very large datasets before feeding them into layers.
//!
//! ## Coordinate space
//!
//! All distance calculations are performed in the **degree plane**
//! (lon as X, lat as Y).  This is a deliberate trade-off:
//!
//! - At the equator, 1 degree of longitude ~ 111 km and 1 degree of
//!   latitude ~ 111 km, so the metric is approximately isotropic.
//! - At higher latitudes, longitude degrees shrink by `cos(lat)`, so
//!   east-west features will be simplified more aggressively than
//!   north-south features at the same `epsilon`.  For map display
//!   purposes this is acceptable because Web Mercator stretches
//!   east-west distances proportionally, cancelling the distortion on
//!   screen.
//! - For truly metric simplification, project to Web Mercator meters
//!   first, simplify, then unproject.  The engine may do this
//!   internally in a future release.
//!
//! ## Altitude
//!
//! The altitude component (`GeoCoord::alt`) is **carried through** but
//! does **not** participate in the distance calculation.  A vertex is
//! retained or discarded based solely on its 2D (lat/lon) deviation
//! from the polyline spine.
//!
//! ## Algorithm complexity
//!
//! - Worst case: O(n^2) when every recursive split isolates a single
//!   point (e.g. a zigzag polyline).
//! - Average case on natural geographic data: O(n log n).
//!
//! [dp]: https://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm
// ---------------------------------------------------------------------------

use rustial_math::GeoCoord;

// ---------------------------------------------------------------------------
// Public API
// ---------------------------------------------------------------------------

/// Simplify a polyline using the Douglas-Peucker algorithm.
///
/// `epsilon` is the maximum allowed perpendicular distance **in degrees**
/// (see [module docs](self) for the coordinate-space rationale).  Smaller
/// values retain more detail; zero retains every vertex.
///
/// # Returns
///
/// A new `Vec<GeoCoord>` containing only the vertices that exceed the
/// tolerance threshold, plus the first and last vertices which are
/// always retained.
///
/// # Edge cases
///
/// | Input | Result |
/// |-------|--------|
/// | empty slice | empty `Vec` |
/// | 1 vertex | that vertex (copied) |
/// | 2 vertices | both vertices (copied) |
/// | `epsilon <= 0.0` | all vertices (no simplification) |
/// | `epsilon` is NaN | all vertices (no simplification -- NaN comparisons always false) |
///
/// # Example
///
/// ```
/// use rustial_engine::simplify_douglas_peucker;
/// use rustial_engine::GeoCoord;
///
/// let line = vec![
///     GeoCoord::from_lat_lon(0.0, 0.0),
///     GeoCoord::from_lat_lon(0.0, 0.5),  // collinear -- will be removed
///     GeoCoord::from_lat_lon(0.0, 1.0),
/// ];
/// let simplified = simplify_douglas_peucker(&line, 0.01);
/// assert_eq!(simplified.len(), 2);
/// ```
pub fn simplify_douglas_peucker(coords: &[GeoCoord], epsilon: f64) -> Vec<GeoCoord> {
    if coords.len() <= 2 {
        return coords.to_vec();
    }

    let mut keep = vec![false; coords.len()];
    keep[0] = true;
    keep[coords.len() - 1] = true;

    dp_recurse(coords, 0, coords.len() - 1, epsilon, &mut keep);

    coords
        .iter()
        .zip(keep.iter())
        .filter_map(|(c, &k)| if k { Some(*c) } else { None })
        .collect()
}

/// Simplify a closed polygon ring using Douglas-Peucker.
///
/// Like [`simplify_douglas_peucker`] but guarantees the result remains a
/// valid ring:
///
/// - The first and last vertices are always retained.
/// - If the ring has a closing duplicate (first == last), it is preserved
///   in the output so the ring stays closed.
/// - The result is guaranteed to have at least 3 distinct vertices (plus
///   the optional closing duplicate), because a polygon with fewer than
///   3 vertices is degenerate.  If simplification would reduce the ring
///   below this minimum, the original ring is returned unchanged.
///
/// # Example
///
/// ```
/// use rustial_engine::simplify_polygon_ring;
/// use rustial_engine::GeoCoord;
///
/// let ring = vec![
///     GeoCoord::from_lat_lon(0.0, 0.0),
///     GeoCoord::from_lat_lon(0.0, 1.0),
///     GeoCoord::from_lat_lon(1.0, 1.0),
///     GeoCoord::from_lat_lon(1.0, 0.0),
///     GeoCoord::from_lat_lon(0.0, 0.0), // closing vertex
/// ];
/// let simplified = simplify_polygon_ring(&ring, 0.01);
/// // All corners are significant; result keeps them.
/// assert!(simplified.len() >= 4); // 4 corners or 4 + closing
/// ```
pub fn simplify_polygon_ring(coords: &[GeoCoord], epsilon: f64) -> Vec<GeoCoord> {
    if coords.len() <= 3 {
        return coords.to_vec();
    }

    // Detect whether the ring has a closing duplicate vertex.
    let n = coords.len();
    let has_closing = n > 1
        && (coords[0].lat - coords[n - 1].lat).abs() < 1e-12
        && (coords[0].lon - coords[n - 1].lon).abs() < 1e-12;

    // Work on the open ring (without closing duplicate).
    let open = if has_closing {
        &coords[..n - 1]
    } else {
        coords
    };

    let simplified = simplify_douglas_peucker(open, epsilon);

    // A polygon ring must have at least 3 distinct vertices.
    if simplified.len() < 3 {
        return coords.to_vec();
    }

    if has_closing {
        // Re-close the ring.
        let mut result = simplified;
        result.push(result[0]);
        result
    } else {
        simplified
    }
}

// ---------------------------------------------------------------------------
// Internal: recursive Douglas-Peucker
// ---------------------------------------------------------------------------

/// Recursive core of the Douglas-Peucker algorithm.
///
/// Finds the vertex between `start` and `end` with the greatest
/// perpendicular distance to the line segment `coords[start]`--
/// `coords[end]`.  If that distance exceeds `epsilon`, the vertex is
/// marked for retention and the algorithm recurses on both halves.
fn dp_recurse(coords: &[GeoCoord], start: usize, end: usize, epsilon: f64, keep: &mut [bool]) {
    if end <= start + 1 {
        return;
    }

    let mut max_dist = 0.0_f64;
    let mut max_idx = start;

    let a = coords[start];
    let b = coords[end];

    for (i, coord) in coords.iter().enumerate().take(end).skip(start + 1) {
        let d = perpendicular_distance(coord, &a, &b);
        if d > max_dist {
            max_dist = d;
            max_idx = i;
        }
    }

    if max_dist > epsilon {
        keep[max_idx] = true;
        dp_recurse(coords, start, max_idx, epsilon, keep);
        dp_recurse(coords, max_idx, end, epsilon, keep);
    }
}

// ---------------------------------------------------------------------------
// Internal: geometry helpers
// ---------------------------------------------------------------------------

/// Perpendicular distance from point `p` to the line through `a`--`b`,
/// computed in the degree plane (lon = X, lat = Y).
///
/// ## Formula
///
/// Given the line direction `d = b - a`, the signed area of the
/// parallelogram formed by `d` and `(p - a)` is:
///
/// ```text
/// cross = (p.lon - a.lon) * (b.lat - a.lat)
///       - (p.lat - a.lat) * (b.lon - a.lon)
/// ```
///
/// The perpendicular distance is `|cross| / |d|`.
///
/// ## Degenerate case
///
/// When `a` and `b` coincide (segment length < 1e-12 degrees, i.e.
/// sub-millimeter), the function falls back to the Euclidean distance
/// from `p` to `a`.
fn perpendicular_distance(p: &GeoCoord, a: &GeoCoord, b: &GeoCoord) -> f64 {
    let dx = b.lon - a.lon;
    let dy = b.lat - a.lat;
    let len_sq = dx * dx + dy * dy;

    // Degenerate: a and b are effectively the same point.
    // 1e-24 in degree^2 corresponds to ~1e-12 degrees (~0.1 mm).
    if len_sq < 1e-24 {
        let ex = p.lon - a.lon;
        let ey = p.lat - a.lat;
        return (ex * ex + ey * ey).sqrt();
    }

    let cross = ((p.lon - a.lon) * dy - (p.lat - a.lat) * dx).abs();
    cross / len_sq.sqrt()
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;

    // -- simplify_douglas_peucker -----------------------------------------

    #[test]
    fn empty_input() {
        let result = simplify_douglas_peucker(&[], 0.1);
        assert!(result.is_empty());
    }

    #[test]
    fn single_point() {
        let coords = vec![GeoCoord::from_lat_lon(51.0, 17.0)];
        let result = simplify_douglas_peucker(&coords, 0.1);
        assert_eq!(result.len(), 1);
        assert_eq!(result[0], coords[0]);
    }

    #[test]
    fn no_simplification_for_two_points() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(1.0, 1.0),
        ];
        let result = simplify_douglas_peucker(&coords, 0.1);
        assert_eq!(result.len(), 2);
    }

    #[test]
    fn straight_line_simplified() {
        // Three collinear points: the middle one should be removed.
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 0.5),
            GeoCoord::from_lat_lon(0.0, 1.0),
        ];
        let result = simplify_douglas_peucker(&coords, 0.01);
        assert_eq!(result.len(), 2);
        assert_eq!(result[0], coords[0]);
        assert_eq!(result[1], coords[2]);
    }

    #[test]
    fn bent_line_keeps_vertex() {
        // A 90-degree bend: the apex must be retained at tight tolerance.
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(1.0, 0.5),
            GeoCoord::from_lat_lon(0.0, 1.0),
        ];
        let result = simplify_douglas_peucker(&coords, 0.01);
        assert_eq!(result.len(), 3);
    }

    #[test]
    fn large_epsilon_simplifies_everything() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.1, 0.5),
            GeoCoord::from_lat_lon(0.0, 1.0),
            GeoCoord::from_lat_lon(-0.1, 1.5),
            GeoCoord::from_lat_lon(0.0, 2.0),
        ];
        let result = simplify_douglas_peucker(&coords, 10.0);
        assert_eq!(result.len(), 2);
    }

    #[test]
    fn zero_epsilon_retains_all() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.1, 0.5),
            GeoCoord::from_lat_lon(0.0, 1.0),
        ];
        let result = simplify_douglas_peucker(&coords, 0.0);
        assert_eq!(result.len(), 3);
    }

    #[test]
    fn negative_epsilon_retains_all() {
        // Negative epsilon means the `max_dist > epsilon` check always
        // passes (any positive distance > negative), so all vertices
        // are retained.
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.1, 0.5),
            GeoCoord::from_lat_lon(0.0, 1.0),
        ];
        let result = simplify_douglas_peucker(&coords, -1.0);
        assert_eq!(result.len(), 3);
    }

    #[test]
    fn altitude_is_preserved() {
        let coords = vec![
            GeoCoord::new(0.0, 0.0, 100.0),
            GeoCoord::new(1.0, 0.5, 200.0), // significant bend
            GeoCoord::new(0.0, 1.0, 300.0),
        ];
        let result = simplify_douglas_peucker(&coords, 0.01);
        assert_eq!(result.len(), 3);
        assert_eq!(result[0].alt, 100.0);
        assert_eq!(result[1].alt, 200.0);
        assert_eq!(result[2].alt, 300.0);
    }

    #[test]
    fn many_collinear_points() {
        // 100 points on a straight east-west line at lat=0.
        let coords: Vec<GeoCoord> = (0..100)
            .map(|i| GeoCoord::from_lat_lon(0.0, i as f64 * 0.01))
            .collect();
        let result = simplify_douglas_peucker(&coords, 0.001);
        assert_eq!(result.len(), 2, "all collinear points should be removed");
        assert_eq!(result[0], coords[0]);
        assert_eq!(result[1], coords[99]);
    }

    #[test]
    fn zigzag_retains_peaks() {
        // A zigzag with 5 peaks should retain all of them at tight epsilon.
        let mut coords = Vec::new();
        for i in 0..11 {
            let lon = i as f64 * 0.1;
            let lat = if i % 2 == 0 { 0.0 } else { 1.0 };
            coords.push(GeoCoord::from_lat_lon(lat, lon));
        }
        let result = simplify_douglas_peucker(&coords, 0.01);
        assert_eq!(
            result.len(),
            coords.len(),
            "tight epsilon retains all zigzag vertices"
        );
    }

    // -- simplify_polygon_ring --------------------------------------------

    #[test]
    fn ring_preserves_square() {
        let ring = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 0.0), // closing
        ];
        let result = simplify_polygon_ring(&ring, 0.01);
        // All 4 corners are significant.
        assert_eq!(result.len(), 5, "4 corners + closing vertex");
        // Must still be closed.
        assert_eq!(result[0], result[result.len() - 1]);
    }

    #[test]
    fn ring_simplifies_collinear_edges() {
        // A rectangle with extra collinear midpoints on each edge.
        let ring = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 0.5), // collinear
            GeoCoord::from_lat_lon(0.0, 1.0),
            GeoCoord::from_lat_lon(0.5, 1.0), // collinear
            GeoCoord::from_lat_lon(1.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 0.5), // collinear
            GeoCoord::from_lat_lon(1.0, 0.0),
            GeoCoord::from_lat_lon(0.5, 0.0), // collinear
            GeoCoord::from_lat_lon(0.0, 0.0), // closing
        ];
        let result = simplify_polygon_ring(&ring, 0.01);
        // The result must be shorter than the original (collinear midpoints
        // are removed) and must remain a valid closed ring.
        assert!(
            result.len() < ring.len(),
            "expected fewer vertices, got {} (original {})",
            result.len(),
            ring.len()
        );
        assert!(
            result.len() >= 4,
            "ring must have at least 3 corners + close"
        );
        assert_eq!(result[0], result[result.len() - 1], "ring must be closed");
    }

    #[test]
    fn ring_open_without_closing_vertex() {
        // An open ring (no closing duplicate).
        let ring = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 0.0),
        ];
        let result = simplify_polygon_ring(&ring, 0.01);
        assert_eq!(result.len(), 4);
    }

    #[test]
    fn ring_too_small_is_returned_unchanged() {
        // 3 vertices (triangle) cannot be further simplified.
        let ring = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 0.0),
        ];
        let result = simplify_polygon_ring(&ring, 100.0);
        assert_eq!(result.len(), 3);
    }

    #[test]
    fn ring_degenerate_returns_unchanged() {
        // 2 vertices: degenerate ring, returned as-is.
        let ring = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(1.0, 1.0),
        ];
        let result = simplify_polygon_ring(&ring, 0.01);
        assert_eq!(result.len(), 2);
    }

    // -- perpendicular_distance -------------------------------------------

    #[test]
    fn distance_on_line_is_zero() {
        let a = GeoCoord::from_lat_lon(0.0, 0.0);
        let b = GeoCoord::from_lat_lon(0.0, 1.0);
        let p = GeoCoord::from_lat_lon(0.0, 0.5); // on the line
        let d = perpendicular_distance(&p, &a, &b);
        assert!(d < 1e-12, "point on line should have ~0 distance, got {d}");
    }

    #[test]
    fn distance_off_line_is_positive() {
        let a = GeoCoord::from_lat_lon(0.0, 0.0);
        let b = GeoCoord::from_lat_lon(0.0, 1.0);
        let p = GeoCoord::from_lat_lon(1.0, 0.5); // 1 degree north of line
        let d = perpendicular_distance(&p, &a, &b);
        assert!((d - 1.0).abs() < 1e-10, "expected ~1.0 degree, got {d}");
    }

    #[test]
    fn distance_coincident_endpoints() {
        // When a == b, distance degenerates to Euclidean from p to a.
        let a = GeoCoord::from_lat_lon(0.0, 0.0);
        let b = GeoCoord::from_lat_lon(0.0, 0.0);
        let p = GeoCoord::from_lat_lon(3.0, 4.0);
        let d = perpendicular_distance(&p, &a, &b);
        assert!((d - 5.0).abs() < 1e-10, "expected 5.0, got {d}");
    }
}