oxigdal-algorithms 0.1.4

//! Geometric union operations
//!
//! This module implements robust geometric union algorithms for
//! combining geometric features. Union operations are fundamental
//! for spatial aggregation, dissolve operations, and cartographic
//! generalization.
//!
//! # Algorithms
//!
//! - **Polygon-Polygon Union**: Weiler-Atherton clipping
//! - **Cascaded Union**: Efficient union of multiple polygons
//! - **Boundary Merging**: Combines overlapping polygons
//!
//! # Examples
//!
//! ```
//! use oxigdal_algorithms::vector::{union_polygon, Coordinate, LineString, Polygon};
//! # use oxigdal_algorithms::error::Result;
//!
//! # fn main() -> Result<()> {
//! let coords1 = vec![
//!     Coordinate::new_2d(0.0, 0.0),
//!     Coordinate::new_2d(5.0, 0.0),
//!     Coordinate::new_2d(5.0, 5.0),
//!     Coordinate::new_2d(0.0, 5.0),
//!     Coordinate::new_2d(0.0, 0.0),
//! ];
//! let ext1 = LineString::new(coords1)?;
//! let poly1 = Polygon::new(ext1, vec![])?;
//! # Ok(())
//! # }
//! ```

use crate::error::{AlgorithmError, Result};
use crate::vector::intersection::{intersect_linestrings, point_in_polygon};
use crate::vector::pool::{PoolGuard, get_pooled_polygon};
use oxigdal_core::vector::{Coordinate, Polygon};

#[cfg(feature = "std")]
use std::vec::Vec;

/// Computes union of two polygons
///
/// Returns a simplified union polygon. This implementation handles:
/// - Disjoint polygons (returns both)
/// - One polygon containing the other
/// - Overlapping polygons (simplified boundary merge)
///
/// # Arguments
///
/// * `poly1` - First polygon
/// * `poly2` - Second polygon
///
/// # Returns
///
/// Vector of polygons representing the union
///
/// # Errors
///
/// Returns error if polygons are invalid
pub fn union_polygon(poly1: &Polygon, poly2: &Polygon) -> Result<Vec<Polygon>> {
    // Check validity
    if poly1.exterior.coords.len() < 4 {
        return Err(AlgorithmError::InsufficientData {
            operation: "union_polygon",
            message: "poly1 exterior must have at least 4 coordinates".to_string(),
        });
    }

    if poly2.exterior.coords.len() < 4 {
        return Err(AlgorithmError::InsufficientData {
            operation: "union_polygon",
            message: "poly2 exterior must have at least 4 coordinates".to_string(),
        });
    }

    // Compute bounding boxes for quick checks
    let bounds1 = poly1.bounds();
    let bounds2 = poly2.bounds();

    if let (Some((min_x1, min_y1, max_x1, max_y1)), Some((min_x2, min_y2, max_x2, max_y2))) =
        (bounds1, bounds2)
    {
        // Check if bounding boxes don't overlap
        if max_x1 < min_x2 || max_x2 < min_x1 || max_y1 < min_y2 || max_y2 < min_y1 {
            // Disjoint - return both polygons
            return Ok(vec![poly1.clone(), poly2.clone()]);
        }
    }

    // Find all intersection points between boundaries
    let intersection_points = intersect_linestrings(&poly1.exterior, &poly2.exterior)?;

    if intersection_points.is_empty() {
        // Either one contains the other, or they're disjoint
        // Check containment
        if point_in_polygon(&poly1.exterior.coords[0], poly2)? {
            // poly1 is completely inside poly2
            return Ok(vec![poly2.clone()]);
        }

        if point_in_polygon(&poly2.exterior.coords[0], poly1)? {
            // poly2 is completely inside poly1
            return Ok(vec![poly1.clone()]);
        }

        // Disjoint - return both
        return Ok(vec![poly1.clone(), poly2.clone()]);
    }

    // Delegate to Weiler-Atherton clipping engine for the overlapping case
    crate::vector::clipping::clip_polygons(
        poly1,
        poly2,
        crate::vector::clipping::ClipOperation::Union,
    )
}

/// Computes union of multiple polygons (alias for clarity)
///
/// # Arguments
///
/// * `polygons` - Vector of polygons to union
///
/// # Returns
///
/// Vector of polygons representing the union
///
/// # Errors
///
/// Returns error if any polygon is invalid
pub fn union_polygons(polygons: &[Polygon]) -> Result<Vec<Polygon>> {
    cascaded_union(polygons)
}

/// Computes cascaded union of multiple polygons
///
/// Uses a hierarchical approach to efficiently compute the union
/// of many polygons. More efficient than pairwise union for large
/// collections.
///
/// # Arguments
///
/// * `polygons` - Vector of polygons to union
///
/// # Returns
///
/// Vector of polygons representing the union
///
/// # Errors
///
/// Returns error if any polygon is invalid
pub fn cascaded_union(polygons: &[Polygon]) -> Result<Vec<Polygon>> {
    if polygons.is_empty() {
        return Ok(vec![]);
    }

    if polygons.len() == 1 {
        return Ok(vec![polygons[0].clone()]);
    }

    // Validate all polygons
    for (i, poly) in polygons.iter().enumerate() {
        if poly.exterior.coords.len() < 4 {
            return Err(AlgorithmError::InsufficientData {
                operation: "cascaded_union",
                message: format!("polygon {} exterior must have at least 4 coordinates", i),
            });
        }
    }

    // For cascaded union of potentially disjoint polygons, we use a simpler
    // approach: attempt pairwise unions and collect results.
    // This avoids infinite loops when polygons don't actually merge.
    let mut current = polygons.to_vec();
    let max_iterations = (polygons.len() as f64).log2().ceil() as usize + 1;
    let mut iteration = 0;

    while current.len() > 1 && iteration < max_iterations {
        let mut next_level = Vec::new();
        let initial_count = current.len();

        // Process pairs
        let mut i = 0;
        while i < current.len() {
            if i + 1 < current.len() {
                // Union pair
                let union_result = union_polygon(&current[i], &current[i + 1])?;
                next_level.extend(union_result);
                i += 2;
            } else {
                // Odd one out
                next_level.push(current[i].clone());
                i += 1;
            }
        }

        current = next_level;
        iteration += 1;

        // If no progress was made (all disjoint), stop to avoid infinite loop
        if current.len() >= initial_count {
            break;
        }
    }

    Ok(current)
}

/// Merges touching or overlapping polygons
///
/// Combines polygons that share boundaries or overlap.
/// Useful for dissolve operations.
///
/// # Arguments
///
/// * `polygons` - Vector of polygons to merge
/// * `tolerance` - Distance tolerance for considering polygons as touching
///
/// # Returns
///
/// Vector of merged polygons
///
/// # Errors
///
/// Returns error if any polygon is invalid
pub fn merge_polygons(polygons: &[Polygon], tolerance: f64) -> Result<Vec<Polygon>> {
    if polygons.is_empty() {
        return Ok(vec![]);
    }

    if tolerance < 0.0 {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "tolerance",
            message: "tolerance must be non-negative".to_string(),
        });
    }

    // Simplified implementation: group by proximity and union each group
    let mut merged = Vec::new();
    let mut used = vec![false; polygons.len()];

    for i in 0..polygons.len() {
        if used[i] {
            continue;
        }

        let mut group = vec![polygons[i].clone()];
        used[i] = true;

        // Find all polygons that touch or overlap with current group
        for j in (i + 1)..polygons.len() {
            if used[j] {
                continue;
            }

            // Check if any polygon in group is close to polygons[j]
            let mut is_close = false;
            for poly in &group {
                if polygons_are_close(poly, &polygons[j], tolerance) {
                    is_close = true;
                    break;
                }
            }

            if is_close {
                group.push(polygons[j].clone());
                used[j] = true;
            }
        }

        // Union the group
        let union_result = cascaded_union(&group)?;
        merged.extend(union_result);
    }

    Ok(merged)
}

/// Checks if two polygons are close (within tolerance)
fn polygons_are_close(poly1: &Polygon, poly2: &Polygon, tolerance: f64) -> bool {
    // Check bounding box proximity
    if let (Some(b1), Some(b2)) = (poly1.bounds(), poly2.bounds()) {
        let (min_x1, min_y1, max_x1, max_y1) = b1;
        let (min_x2, min_y2, max_x2, max_y2) = b2;

        // Check if bounding boxes are within tolerance
        let x_gap = if max_x1 < min_x2 {
            min_x2 - max_x1
        } else if max_x2 < min_x1 {
            min_x1 - max_x2
        } else {
            0.0
        };

        let y_gap = if max_y1 < min_y2 {
            min_y2 - max_y1
        } else if max_y2 < min_y1 {
            min_y1 - max_y2
        } else {
            0.0
        };

        x_gap <= tolerance && y_gap <= tolerance
    } else {
        false
    }
}

/// Computes the signed area of a polygon
fn compute_polygon_area(polygon: &Polygon) -> f64 {
    let mut area = ring_area(&polygon.exterior.coords);

    // Subtract hole areas
    for hole in &polygon.interiors {
        area -= ring_area(&hole.coords);
    }

    area.abs()
}

/// Computes the signed area of a ring using the shoelace formula
fn ring_area(coords: &[Coordinate]) -> f64 {
    if coords.len() < 3 {
        return 0.0;
    }

    let mut area = 0.0;
    let n = coords.len();

    for i in 0..n {
        let j = (i + 1) % n;
        area += coords[i].x * coords[j].y;
        area -= coords[j].x * coords[i].y;
    }

    area / 2.0
}

/// Computes the convex hull of a set of points
///
/// Uses Graham scan algorithm.
///
/// # Arguments
///
/// * `points` - Vector of coordinates
///
/// # Returns
///
/// Vector of coordinates forming the convex hull (counter-clockwise)
///
/// # Errors
///
/// Returns error if fewer than 3 points provided
pub fn convex_hull(points: &[Coordinate]) -> Result<Vec<Coordinate>> {
    if points.len() < 3 {
        return Err(AlgorithmError::InsufficientData {
            operation: "convex_hull",
            message: "need at least 3 points for convex hull".to_string(),
        });
    }

    // Find the point with lowest y-coordinate (and leftmost if tied)
    let mut lowest_idx = 0;
    for i in 1..points.len() {
        if points[i].y < points[lowest_idx].y
            || (points[i].y == points[lowest_idx].y && points[i].x < points[lowest_idx].x)
        {
            lowest_idx = i;
        }
    }

    let pivot = points[lowest_idx];
    let mut sorted_points: Vec<(usize, Coordinate, f64)> = points
        .iter()
        .enumerate()
        .filter(|(i, _)| *i != lowest_idx)
        .map(|(i, p)| {
            let angle = (p.y - pivot.y).atan2(p.x - pivot.x);
            (i, *p, angle)
        })
        .collect();

    // Sort by polar angle
    sorted_points.sort_by(|a, b| a.2.partial_cmp(&b.2).unwrap_or(std::cmp::Ordering::Equal));

    let mut hull = vec![pivot];

    for (_, point, _) in sorted_points {
        // Remove points that would create a clockwise turn
        while hull.len() >= 2 {
            let n = hull.len();
            let cross = cross_product(&hull[n - 2], &hull[n - 1], &point);
            if cross <= 0.0 {
                hull.pop();
            } else {
                break;
            }
        }
        hull.push(point);
    }

    Ok(hull)
}

/// Computes cross product for orientation test
fn cross_product(o: &Coordinate, a: &Coordinate, b: &Coordinate) -> f64 {
    (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x)
}

//
// Pooled union operations for reduced allocations
//

/// Computes union of two polygons using object pooling
///
/// This is the pooled version that returns the first result polygon from
/// a thread-local pool. For multiple result polygons (disjoint case),
/// only the first is pooled.
///
/// # Arguments
///
/// * `poly1` - First polygon
/// * `poly2` - Second polygon
///
/// # Returns
///
/// A pooled polygon guard representing the union (first result if multiple)
///
/// # Errors
///
/// Returns error if polygons are invalid
///
/// # Example
///
/// ```
/// use oxigdal_algorithms::vector::{union_polygon_pooled, Coordinate, LineString, Polygon};
///
/// let coords1 = vec![
///     Coordinate::new_2d(0.0, 0.0),
///     Coordinate::new_2d(5.0, 0.0),
///     Coordinate::new_2d(5.0, 5.0),
///     Coordinate::new_2d(0.0, 5.0),
///     Coordinate::new_2d(0.0, 0.0),
/// ];
/// let ext1 = LineString::new(coords1)?;
/// let poly1 = Polygon::new(ext1, vec![])?;
/// # let coords2 = vec![
/// #     Coordinate::new_2d(3.0, 0.0),
/// #     Coordinate::new_2d(8.0, 0.0),
/// #     Coordinate::new_2d(8.0, 5.0),
/// #     Coordinate::new_2d(3.0, 5.0),
/// #     Coordinate::new_2d(3.0, 0.0),
/// # ];
/// # let ext2 = LineString::new(coords2)?;
/// # let poly2 = Polygon::new(ext2, vec![])?;
/// let result = union_polygon_pooled(&poly1, &poly2)?;
/// # Ok::<(), oxigdal_algorithms::error::AlgorithmError>(())
/// ```
pub fn union_polygon_pooled(
    poly1: &Polygon,
    poly2: &Polygon,
) -> Result<PoolGuard<'static, Polygon>> {
    let results = union_polygon(poly1, poly2)?;

    // Get a pooled polygon and copy the first result into it
    if let Some(result) = results.first() {
        let mut poly = get_pooled_polygon();
        poly.exterior = result.exterior.clone();
        poly.interiors = result.interiors.clone();
        Ok(poly)
    } else {
        Err(AlgorithmError::InsufficientData {
            operation: "union_polygon_pooled",
            message: "union resulted in no polygons".to_string(),
        })
    }
}

/// Computes cascaded union of multiple polygons using object pooling
///
/// Efficiently combines multiple polygons into a single union polygon
/// using object pooling to reduce allocations.
///
/// # Arguments
///
/// * `polygons` - Slice of polygons to union
///
/// # Returns
///
/// A pooled polygon guard representing the union
///
/// # Errors
///
/// Returns error if any polygon is invalid
///
/// # Example
///
/// ```
/// use oxigdal_algorithms::vector::{union_polygons_pooled, Coordinate, LineString, Polygon};
///
/// let coords = vec![
///     Coordinate::new_2d(0.0, 0.0),
///     Coordinate::new_2d(5.0, 0.0),
///     Coordinate::new_2d(5.0, 5.0),
///     Coordinate::new_2d(0.0, 5.0),
///     Coordinate::new_2d(0.0, 0.0),
/// ];
/// let ext = LineString::new(coords)?;
/// let poly = Polygon::new(ext, vec![])?;
/// let polygons = vec![poly];
/// let result = union_polygons_pooled(&polygons)?;
/// # Ok::<(), oxigdal_algorithms::error::AlgorithmError>(())
/// ```
pub fn union_polygons_pooled(polygons: &[Polygon]) -> Result<PoolGuard<'static, Polygon>> {
    let results = union_polygons(polygons)?;

    // Get a pooled polygon and copy the first result into it
    if let Some(result) = results.first() {
        let mut poly = get_pooled_polygon();
        poly.exterior = result.exterior.clone();
        poly.interiors = result.interiors.clone();
        Ok(poly)
    } else {
        Err(AlgorithmError::InsufficientData {
            operation: "union_polygons_pooled",
            message: "union resulted in no polygons".to_string(),
        })
    }
}

/// Computes cascaded union using object pooling (alias)
///
/// This is an alias for `union_polygons_pooled` for consistency with
/// the non-pooled API.
pub fn cascaded_union_pooled(polygons: &[Polygon]) -> Result<PoolGuard<'static, Polygon>> {
    union_polygons_pooled(polygons)
}

/// Computes convex hull and returns as pooled polygon
///
/// Computes the convex hull of a set of points and returns it as a
/// pooled polygon to reduce allocations.
///
/// # Arguments
///
/// * `points` - Points to compute convex hull from
///
/// # Returns
///
/// A pooled polygon representing the convex hull
///
/// # Errors
///
/// Returns error if fewer than 3 points provided
///
/// # Example
///
/// ```
/// use oxigdal_algorithms::vector::{convex_hull_pooled, Coordinate};
///
/// let points = vec![
///     Coordinate::new_2d(0.0, 0.0),
///     Coordinate::new_2d(5.0, 0.0),
///     Coordinate::new_2d(2.5, 5.0),
///     Coordinate::new_2d(2.5, 2.5),
/// ];
/// let hull = convex_hull_pooled(&points)?;
/// # Ok::<(), oxigdal_algorithms::error::AlgorithmError>(())
/// ```
pub fn convex_hull_pooled(points: &[Coordinate]) -> Result<PoolGuard<'static, Polygon>> {
    use oxigdal_core::vector::LineString;

    let hull_coords = convex_hull(points)?;

    // Create a polygon from the hull coordinates
    let mut poly = get_pooled_polygon();
    poly.exterior.coords.clear();
    poly.exterior.coords.extend(hull_coords);

    // Close the ring if not already closed
    if let (Some(&first), Some(&last)) = (poly.exterior.coords.first(), poly.exterior.coords.last())
    {
        if first.x != last.x || first.y != last.y {
            poly.exterior.coords.push(first);
        }
    }

    // Ensure minimum 4 points for valid polygon
    while poly.exterior.len() < 4 {
        if let Some(&first) = poly.exterior.coords.first() {
            poly.exterior.coords.push(first);
        }
    }

    Ok(poly)
}

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

    fn create_square(x: f64, y: f64, size: f64) -> Result<Polygon> {
        let coords = vec![
            Coordinate::new_2d(x, y),
            Coordinate::new_2d(x + size, y),
            Coordinate::new_2d(x + size, y + size),
            Coordinate::new_2d(x, y + size),
            Coordinate::new_2d(x, y),
        ];
        let exterior = LineString::new(coords).map_err(|e| AlgorithmError::Core(e))?;
        Polygon::new(exterior, vec![]).map_err(|e| AlgorithmError::Core(e))
    }

    #[test]
    fn test_union_polygon_disjoint() {
        let poly1 = create_square(0.0, 0.0, 5.0);
        let poly2 = create_square(10.0, 10.0, 5.0);

        assert!(poly1.is_ok() && poly2.is_ok());
        if let (Ok(p1), Ok(p2)) = (poly1, poly2) {
            let result = union_polygon(&p1, &p2);
            assert!(result.is_ok());
            if let Ok(polys) = result {
                assert_eq!(polys.len(), 2); // Disjoint - return both
            }
        }
    }

    #[test]
    fn test_union_polygon_contained() {
        let poly1 = create_square(0.0, 0.0, 10.0);
        let poly2 = create_square(2.0, 2.0, 3.0);

        assert!(poly1.is_ok() && poly2.is_ok());
        if let (Ok(p1), Ok(p2)) = (poly1, poly2) {
            let result = union_polygon(&p1, &p2);
            assert!(result.is_ok());
            if let Ok(polys) = result {
                assert_eq!(polys.len(), 1); // One contains the other
            }
        }
    }

    #[test]
    fn test_union_polygon_overlapping() {
        let poly1 = create_square(0.0, 0.0, 5.0);
        let poly2 = create_square(3.0, 0.0, 5.0);

        assert!(poly1.is_ok() && poly2.is_ok());
        if let (Ok(p1), Ok(p2)) = (poly1, poly2) {
            let result = union_polygon(&p1, &p2);
            assert!(result.is_ok());
            // Simplified implementation returns one polygon
            if let Ok(polys) = result {
                assert!(!polys.is_empty());
            }
        }
    }

    #[test]
    fn test_cascaded_union_empty() {
        let result = cascaded_union(&[]);
        assert!(result.is_ok());
        if let Ok(polys) = result {
            assert!(polys.is_empty());
        }
    }

    #[test]
    fn test_cascaded_union_single() {
        let poly = create_square(0.0, 0.0, 5.0);
        assert!(poly.is_ok());

        if let Ok(p) = poly {
            let result = cascaded_union(std::slice::from_ref(&p));
            assert!(result.is_ok());
            if let Ok(polys) = result {
                assert_eq!(polys.len(), 1);
            }
        }
    }

    #[test]
    fn test_cascaded_union_multiple() {
        let poly1 = create_square(0.0, 0.0, 5.0);
        let poly2 = create_square(10.0, 0.0, 5.0);
        let poly3 = create_square(20.0, 0.0, 5.0);

        assert!(poly1.is_ok() && poly2.is_ok() && poly3.is_ok());
        if let (Ok(p1), Ok(p2), Ok(p3)) = (poly1, poly2, poly3) {
            let result = cascaded_union(&[p1, p2, p3]);
            assert!(result.is_ok());
            // All disjoint, should return 3 polygons
            if let Ok(polys) = result {
                assert!(!polys.is_empty());
            }
        }
    }

    #[test]
    fn test_merge_polygons() {
        let poly1 = create_square(0.0, 0.0, 5.0);
        let poly2 = create_square(10.0, 0.0, 5.0);

        assert!(poly1.is_ok() && poly2.is_ok());
        if let (Ok(p1), Ok(p2)) = (poly1, poly2) {
            let result = merge_polygons(&[p1, p2], 1.0);
            assert!(result.is_ok());
        }
    }

    #[test]
    fn test_compute_polygon_area() {
        let poly = create_square(0.0, 0.0, 10.0);
        assert!(poly.is_ok());

        if let Ok(p) = poly {
            let area = compute_polygon_area(&p);
            assert!((area - 100.0).abs() < 1e-10);
        }
    }

    #[test]
    fn test_ring_area() {
        let coords = vec![
            Coordinate::new_2d(0.0, 0.0),
            Coordinate::new_2d(10.0, 0.0),
            Coordinate::new_2d(10.0, 10.0),
            Coordinate::new_2d(0.0, 10.0),
            Coordinate::new_2d(0.0, 0.0),
        ];
        let area = ring_area(&coords);
        assert!((area.abs() - 100.0).abs() < 1e-10);
    }

    #[test]
    fn test_convex_hull_triangle() {
        let points = vec![
            Coordinate::new_2d(0.0, 0.0),
            Coordinate::new_2d(4.0, 0.0),
            Coordinate::new_2d(2.0, 3.0),
        ];
        let result = convex_hull(&points);
        assert!(result.is_ok());
        if let Ok(hull) = result {
            assert_eq!(hull.len(), 3);
        }
    }

    #[test]
    fn test_convex_hull_square_with_interior() {
        let points = vec![
            Coordinate::new_2d(0.0, 0.0),
            Coordinate::new_2d(4.0, 0.0),
            Coordinate::new_2d(4.0, 4.0),
            Coordinate::new_2d(0.0, 4.0),
            Coordinate::new_2d(2.0, 2.0), // Interior point
        ];
        let result = convex_hull(&points);
        assert!(result.is_ok());
        if let Ok(hull) = result {
            assert_eq!(hull.len(), 4); // Interior point should be excluded
        }
    }

    #[test]
    fn test_convex_hull_insufficient_points() {
        let points = vec![Coordinate::new_2d(0.0, 0.0), Coordinate::new_2d(1.0, 1.0)];
        let result = convex_hull(&points);
        assert!(result.is_err());
    }

    #[test]
    fn test_polygons_are_close() {
        let poly1 = create_square(0.0, 0.0, 5.0);
        let poly2 = create_square(5.5, 0.0, 5.0);

        assert!(poly1.is_ok() && poly2.is_ok());
        if let (Ok(p1), Ok(p2)) = (poly1, poly2) {
            assert!(polygons_are_close(&p1, &p2, 1.0)); // Within tolerance
            assert!(!polygons_are_close(&p1, &p2, 0.1)); // Outside tolerance
        }
    }

    #[test]
    fn test_cascaded_union_many_disjoint_no_infinite_loop() {
        // Regression test for infinite loop bug when all polygons are disjoint.
        // This test ensures cascaded_union terminates quickly even with many
        // disjoint polygons that cannot be merged.
        let mut polygons = Vec::new();

        // Create 10 disjoint squares (well separated)
        for i in 0..10 {
            let x = (i as f64) * 20.0; // 20 units apart = well separated
            if let Ok(poly) = create_square(x, 0.0, 5.0) {
                polygons.push(poly);
            }
        }

        let start = std::time::Instant::now();
        let result = cascaded_union(&polygons);
        let elapsed = start.elapsed();

        // Should complete in well under 1 second (was infinite loop before fix)
        assert!(
            elapsed.as_secs() < 1,
            "cascaded_union took too long: {elapsed:?}"
        );
        assert!(result.is_ok());

        if let Ok(results) = result {
            // All disjoint, so should return all 10 polygons
            assert_eq!(results.len(), 10);
        }
    }
}