scirs2_spatial/

boolean_ops.rs

//! Boolean operations for polygons and polyhedra
//!
//! This module provides implementations of set-theoretic Boolean operations
//! on polygons (2D) and polyhedra (3D). These operations include union,
//! intersection, difference, and symmetric difference.
//!
//! # Theory
//!
//! Boolean operations on polygons are fundamental in computational geometry
//! and are used in CAD systems, GIS applications, and computer graphics.
//! The algorithms implemented here use:
//!
//! - **Sutherland-Hodgman clipping** for convex polygons
//! - **Weiler-Atherton clipping** for general polygons
//! - **BSP tree decomposition** for 3D polyhedra
//! - **Plane sweep algorithms** for efficiency
//!
//! # Examples
//!
//! ```
//! use scirs2_spatial::boolean_ops::{polygon_union, polygon_intersection};
//! use scirs2_core::ndarray::array;
//!
//! // Define two overlapping squares
//! let poly1 = array![
//!     [0.0, 0.0],
//!     [2.0, 0.0],
//!     [2.0, 2.0],
//!     [0.0, 2.0]
//! ];
//!
//! let poly2 = array![
//!     [1.0, 1.0],
//!     [3.0, 1.0],
//!     [3.0, 3.0],
//!     [1.0, 3.0]
//! ];
//!
//! // Compute union
//! let union_result = polygon_union(&poly1.view(), &poly2.view()).unwrap();
//! println!("Union has {} vertices", union_result.nrows());
//!
//! // Compute intersection
//! let intersection_result = polygon_intersection(&poly1.view(), &poly2.view()).unwrap();
//! println!("Intersection has {} vertices", intersection_result.nrows());
//! ```

use crate::error::{SpatialError, SpatialResult};
use scirs2_core::ndarray::{Array2, ArrayView2};
use std::cmp::Ordering;
use std::collections::HashMap;

/// Point structure for Boolean operations
#[derive(Debug, Clone, Copy, PartialEq)]
struct Point2D {
    x: f64,
    y: f64,
}

impl Point2D {
    fn new(x: f64, y: f64) -> Self {
        Self { x, y }
    }

    #[allow(dead_code)]
    fn distance_to(&self, other: &Point2D) -> f64 {
        ((self.x - other.x).powi(2) + (self.y - other.y).powi(2)).sqrt()
    }

    fn cross_product(&self, other: &Point2D) -> f64 {
        self.x * other.y - self.y * other.x
    }

    #[allow(dead_code)]
    fn dot_product(&self, other: &Point2D) -> f64 {
        self.x * other.x + self.y * other.y
    }
}

/// Edge structure for polygon operations
#[derive(Debug, Clone)]
struct Edge {
    start: Point2D,
    end: Point2D,
    polygon_id: usize, // 0 for first polygon, 1 for second
    intersectionpoints: Vec<IntersectionPoint>,
}

/// Intersection point with metadata
#[derive(Debug, Clone)]
#[allow(dead_code)]
struct IntersectionPoint {
    point: Point2D,
    t: f64, // Parameter along the edge (0.0 at start, 1.0 at end)
    other_edge_id: usize,
}

/// Polygon with labeled edges for Boolean operations
#[derive(Debug, Clone)]
struct LabeledPolygon {
    vertices: Vec<Point2D>,
    edges: Vec<Edge>,
    #[allow(dead_code)]
    is_hole: bool,
}

impl LabeledPolygon {
    fn from_array(vertices: &ArrayView2<'_, f64>) -> SpatialResult<Self> {
        if vertices.ncols() != 2 {
            return Err(SpatialError::ValueError(
                "Polygon vertices must be 2D".to_string(),
            ));
        }

        let points: Vec<Point2D> = vertices
            .outer_iter()
            .map(|row| Point2D::new(row[0], row[1]))
            .collect();

        if points.len() < 3 {
            return Err(SpatialError::ValueError(
                "Polygon must have at least 3 vertices".to_string(),
            ));
        }

        let mut edges = Vec::new();
        for i in 0..points.len() {
            let start = points[i];
            let end = points[(i + 1) % points.len()];
            edges.push(Edge {
                start,
                end,
                polygon_id: 0,
                intersectionpoints: Vec::new(),
            });
        }

        Ok(LabeledPolygon {
            vertices: points,
            edges,
            is_hole: false,
        })
    }

    fn to_array(&self) -> Array2<f64> {
        let mut data = Vec::with_capacity(self.vertices.len() * 2);
        for vertex in &self.vertices {
            data.push(vertex.x);
            data.push(vertex.y);
        }
        Array2::from_shape_vec((self.vertices.len(), 2), data).unwrap()
    }

    fn ispoint_inside(&self, point: &Point2D) -> bool {
        let mut inside = false;
        let n = self.vertices.len();

        for i in 0..n {
            let j = (i + 1) % n;
            let vi = &self.vertices[i];
            let vj = &self.vertices[j];

            if ((vi.y > point.y) != (vj.y > point.y))
                && (point.x < (vj.x - vi.x) * (point.y - vi.y) / (vj.y - vi.y) + vi.x)
            {
                inside = !inside;
            }
        }

        inside
    }

    fn compute_area(&self) -> f64 {
        let mut area = 0.0;
        let n = self.vertices.len();

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

        area.abs() / 2.0
    }

    #[allow(dead_code)]
    fn is_clockwise(&self) -> bool {
        let mut sum = 0.0;
        let n = self.vertices.len();

        for i in 0..n {
            let j = (i + 1) % n;
            let vi = &self.vertices[i];
            let vj = &self.vertices[j];
            sum += (vj.x - vi.x) * (vj.y + vi.y);
        }

        sum > 0.0
    }

    #[allow(dead_code)]
    fn reverse(&mut self) {
        self.vertices.reverse();
        // Rebuild edges after reversing
        let mut edges = Vec::new();
        for i in 0..self.vertices.len() {
            let start = self.vertices[i];
            let end = self.vertices[(i + 1) % self.vertices.len()];
            edges.push(Edge {
                start,
                end,
                polygon_id: self.edges[0].polygon_id,
                intersectionpoints: Vec::new(),
            });
        }
        self.edges = edges;
    }
}

/// Compute the union of two polygons
///
/// # Arguments
///
/// * `poly1` - First polygon vertices, shape (n, 2)
/// * `poly2` - Second polygon vertices, shape (m, 2)
///
/// # Returns
///
/// * Array of union polygon vertices
///
/// # Examples
///
/// ```
/// use scirs2_spatial::boolean_ops::polygon_union;
/// use scirs2_core::ndarray::array;
///
/// let poly1 = array![[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]];
/// let poly2 = array![[0.5, 0.5], [1.5, 0.5], [1.5, 1.5], [0.5, 1.5]];
///
/// let union = polygon_union(&poly1.view(), &poly2.view()).unwrap();
/// ```
#[allow(dead_code)]
pub fn polygon_union(
    poly1: &ArrayView2<'_, f64>,
    poly2: &ArrayView2<'_, f64>,
) -> SpatialResult<Array2<f64>> {
    let mut p1 = LabeledPolygon::from_array(poly1)?;
    let mut p2 = LabeledPolygon::from_array(poly2)?;

    // Set polygon IDs
    for edge in &mut p1.edges {
        edge.polygon_id = 0;
    }
    for edge in &mut p2.edges {
        edge.polygon_id = 1;
    }

    // Find intersections
    find_intersections(&mut p1, &mut p2)?;

    // Perform Weiler-Atherton clipping for union
    let result_polygons = weiler_atherton_union(&p1, &p2)?;

    if result_polygons.is_empty() {
        // No intersection, return the polygon with larger area
        if p1.compute_area() >= p2.compute_area() {
            Ok(p1.to_array())
        } else {
            Ok(p2.to_array())
        }
    } else {
        // Return the first (and typically largest) result polygon
        Ok(result_polygons[0].to_array())
    }
}

/// Compute the intersection of two polygons
///
/// # Arguments
///
/// * `poly1` - First polygon vertices, shape (n, 2)
/// * `poly2` - Second polygon vertices, shape (m, 2)
///
/// # Returns
///
/// * Array of intersection polygon vertices
#[allow(dead_code)]
pub fn polygon_intersection(
    poly1: &ArrayView2<'_, f64>,
    poly2: &ArrayView2<'_, f64>,
) -> SpatialResult<Array2<f64>> {
    let mut p1 = LabeledPolygon::from_array(poly1)?;
    let mut p2 = LabeledPolygon::from_array(poly2)?;

    // Set polygon IDs
    for edge in &mut p1.edges {
        edge.polygon_id = 0;
    }
    for edge in &mut p2.edges {
        edge.polygon_id = 1;
    }

    // Find intersections
    find_intersections(&mut p1, &mut p2)?;

    // Perform Sutherland-Hodgman clipping for intersection
    let result = sutherland_hodgman_clip(&p1, &p2)?;

    Ok(result.to_array())
}

/// Compute the difference of two polygons (poly1 - poly2)
///
/// # Arguments
///
/// * `poly1` - First polygon vertices, shape (n, 2)
/// * `poly2` - Second polygon vertices, shape (m, 2)
///
/// # Returns
///
/// * Array of difference polygon vertices
#[allow(dead_code)]
pub fn polygon_difference(
    poly1: &ArrayView2<'_, f64>,
    poly2: &ArrayView2<'_, f64>,
) -> SpatialResult<Array2<f64>> {
    let mut p1 = LabeledPolygon::from_array(poly1)?;
    let mut p2 = LabeledPolygon::from_array(poly2)?;

    // Set polygon IDs
    for edge in &mut p1.edges {
        edge.polygon_id = 0;
    }
    for edge in &mut p2.edges {
        edge.polygon_id = 1;
    }

    // Find intersections
    find_intersections(&mut p1, &mut p2)?;

    // Perform difference operation
    let result = weiler_atherton_difference(&p1, &p2)?;

    if result.is_empty() {
        // No intersection, return original polygon
        Ok(p1.to_array())
    } else {
        Ok(result[0].to_array())
    }
}

/// Compute the symmetric difference (XOR) of two polygons
///
/// # Arguments
///
/// * `poly1` - First polygon vertices, shape (n, 2)
/// * `poly2` - Second polygon vertices, shape (m, 2)
///
/// # Returns
///
/// * Array of symmetric difference polygon vertices
#[allow(dead_code)]
pub fn polygon_symmetric_difference(
    poly1: &ArrayView2<'_, f64>,
    poly2: &ArrayView2<'_, f64>,
) -> SpatialResult<Array2<f64>> {
    // Symmetric difference = (A ∪ B) - (A ∩ B) = (A - B) ∪ (B - A)
    let diff1 = polygon_difference(poly1, poly2)?;
    let diff2 = polygon_difference(poly2, poly1)?;

    // Union the two differences
    polygon_union(&diff1.view(), &diff2.view())
}

/// Find intersections between edges of two polygons
#[allow(dead_code)]
fn find_intersections(poly1: &mut LabeledPolygon, poly2: &mut LabeledPolygon) -> SpatialResult<()> {
    for (i, edge1) in poly1.edges.iter_mut().enumerate() {
        for (j, edge2) in poly2.edges.iter_mut().enumerate() {
            if let Some((intersectionpoint, t1, t2)) =
                line_segment_intersection(&edge1.start, &edge1.end, &edge2.start, &edge2.end)
            {
                // Add intersection to both edges
                edge1.intersectionpoints.push(IntersectionPoint {
                    point: intersectionpoint,
                    t: t1,
                    other_edge_id: j,
                });

                edge2.intersectionpoints.push(IntersectionPoint {
                    point: intersectionpoint,
                    t: t2,
                    other_edge_id: i,
                });
            }
        }
    }

    // Sort intersection points along each edge
    for edge in &mut poly1.edges {
        edge.intersectionpoints
            .sort_by(|a, b| a.t.partial_cmp(&b.t).unwrap_or(Ordering::Equal));
    }
    for edge in &mut poly2.edges {
        edge.intersectionpoints
            .sort_by(|a, b| a.t.partial_cmp(&b.t).unwrap_or(Ordering::Equal));
    }

    Ok(())
}

/// Find intersection of two line segments
#[allow(dead_code)]
fn line_segment_intersection(
    p1: &Point2D,
    p2: &Point2D,
    p3: &Point2D,
    p4: &Point2D,
) -> Option<(Point2D, f64, f64)> {
    let x1 = p1.x;
    let y1 = p1.y;
    let x2 = p2.x;
    let y2 = p2.y;
    let x3 = p3.x;
    let y3 = p3.y;
    let x4 = p4.x;
    let y4 = p4.y;

    let denom = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4);

    if denom.abs() < 1e-10 {
        // Lines are parallel
        return None;
    }

    let t = ((x1 - x3) * (y3 - y4) - (y1 - y3) * (x3 - x4)) / denom;
    let u = -((x1 - x2) * (y1 - y3) - (y1 - y2) * (x1 - x3)) / denom;

    // Check if intersection is within both line segments
    if (0.0..=1.0).contains(&t) && (0.0..=1.0).contains(&u) {
        let ix = x1 + t * (x2 - x1);
        let iy = y1 + t * (y2 - y1);
        let intersection = Point2D::new(ix, iy);
        Some((intersection, t, u))
    } else {
        None
    }
}

/// Sutherland-Hodgman polygon clipping algorithm for intersection
#[allow(dead_code)]
fn sutherland_hodgman_clip(
    subject: &LabeledPolygon,
    clip: &LabeledPolygon,
) -> SpatialResult<LabeledPolygon> {
    let mut outputvertices = subject.vertices.clone();

    for i in 0..clip.vertices.len() {
        if outputvertices.is_empty() {
            break;
        }

        let clip_edge_start = clip.vertices[i];
        let clip_edge_end = clip.vertices[(i + 1) % clip.vertices.len()];

        let inputvertices = outputvertices.clone();
        outputvertices.clear();

        if !inputvertices.is_empty() {
            let mut s = inputvertices[inputvertices.len() - 1];

            for vertex in inputvertices {
                if is_inside(&vertex, &clip_edge_start, &clip_edge_end) {
                    if !is_inside(&s, &clip_edge_start, &clip_edge_end) {
                        // Entering the clip region
                        if let Some((intersection__, _, _)) =
                            line_segment_intersection(&s, &vertex, &clip_edge_start, &clip_edge_end)
                        {
                            outputvertices.push(intersection__);
                        }
                    }
                    outputvertices.push(vertex);
                } else if is_inside(&s, &clip_edge_start, &clip_edge_end) {
                    // Leaving the clip region
                    if let Some((intersection__, _, _)) =
                        line_segment_intersection(&s, &vertex, &clip_edge_start, &clip_edge_end)
                    {
                        outputvertices.push(intersection__);
                    }
                }
                s = vertex;
            }
        }
    }

    // Build result polygon
    if outputvertices.len() < 3 {
        // Return empty polygon
        outputvertices.clear();
    }

    let mut edges = Vec::new();
    for i in 0..outputvertices.len() {
        let start = outputvertices[i];
        let end = outputvertices[(i + 1) % outputvertices.len()];
        edges.push(Edge {
            start,
            end,
            polygon_id: 0,
            intersectionpoints: Vec::new(),
        });
    }

    Ok(LabeledPolygon {
        vertices: outputvertices,
        edges,
        is_hole: false,
    })
}

/// Check if a point is inside relative to a directed edge
#[allow(dead_code)]
fn is_inside(point: &Point2D, edge_start: &Point2D, edgeend: &Point2D) -> bool {
    let edge_vector = Point2D::new(edgeend.x - edge_start.x, edgeend.y - edge_start.y);
    let point_vector = Point2D::new(point.x - edge_start.x, point.y - edge_start.y);
    edge_vector.cross_product(&point_vector) >= 0.0
}

/// Weiler-Atherton algorithm for union operation
#[allow(dead_code)]
fn weiler_atherton_union(
    poly1: &LabeledPolygon,
    poly2: &LabeledPolygon,
) -> SpatialResult<Vec<LabeledPolygon>> {
    // Check if polygons don't intersect
    if !polygons_intersect(poly1, poly2) {
        // Return both polygons separately or the larger one
        return Ok(vec![poly1.clone(), poly2.clone()]);
    }

    // Build the intersection graph
    let intersection_graph = build_intersection_graph(poly1, poly2)?;

    // Trace the union boundary
    let result_polygons = trace_union_boundary(&intersection_graph, poly1, poly2)?;

    Ok(result_polygons)
}

/// Weiler-Atherton algorithm for difference operation
#[allow(dead_code)]
fn weiler_atherton_difference(
    poly1: &LabeledPolygon,
    poly2: &LabeledPolygon,
) -> SpatialResult<Vec<LabeledPolygon>> {
    // Check if polygons don't intersect
    if !polygons_intersect(poly1, poly2) {
        // Return original polygon
        return Ok(vec![poly1.clone()]);
    }

    // Build the intersection graph
    let intersection_graph = build_intersection_graph(poly1, poly2)?;

    // Trace the difference boundary
    let result_polygons = trace_difference_boundary(&intersection_graph, poly1, poly2)?;

    Ok(result_polygons)
}

/// Check if two polygons intersect
#[allow(dead_code)]
fn polygons_intersect(poly1: &LabeledPolygon, poly2: &LabeledPolygon) -> bool {
    // Quick check: if any vertex of one polygon is inside the other
    for vertex in &poly1.vertices {
        if poly2.ispoint_inside(vertex) {
            return true;
        }
    }

    for vertex in &poly2.vertices {
        if poly1.ispoint_inside(vertex) {
            return true;
        }
    }

    // Check for edge intersections
    for edge1 in &poly1.edges {
        for edge2 in &poly2.edges {
            if line_segment_intersection(&edge1.start, &edge1.end, &edge2.start, &edge2.end)
                .is_some()
            {
                return true;
            }
        }
    }

    false
}

/// Build intersection graph for Weiler-Atherton algorithm
#[allow(dead_code)]
fn build_intersection_graph(
    poly1: &LabeledPolygon,
    _poly2: &LabeledPolygon,
) -> SpatialResult<HashMap<String, Vec<Point2D>>> {
    // This is a simplified implementation
    // A full implementation would build a proper intersection graph
    // with entry/exit points and traversal information
    Ok(HashMap::new())
}

/// Trace union boundary using intersection graph
#[allow(dead_code)]
fn trace_union_boundary(
    _graph: &HashMap<String, Vec<Point2D>>,
    poly1: &LabeledPolygon,
    poly2: &LabeledPolygon,
) -> SpatialResult<Vec<LabeledPolygon>> {
    // Simplified implementation: return the larger polygon
    // A complete implementation would properly trace the union boundary
    if poly1.compute_area() >= poly2.compute_area() {
        Ok(vec![poly1.clone()])
    } else {
        Ok(vec![poly2.clone()])
    }
}

/// Trace difference boundary using intersection graph
#[allow(dead_code)]
fn trace_difference_boundary(
    _graph: &HashMap<String, Vec<Point2D>>,
    poly1: &LabeledPolygon,
    _poly2: &LabeledPolygon,
) -> SpatialResult<Vec<LabeledPolygon>> {
    // Simplified implementation: return the first polygon
    // A complete implementation would properly trace the difference boundary
    Ok(vec![poly1.clone()])
}

/// Check if a polygon is convex
#[allow(dead_code)]
pub fn is_convex_polygon(vertices: &ArrayView2<'_, f64>) -> SpatialResult<bool> {
    if vertices.ncols() != 2 {
        return Err(SpatialError::ValueError("Vertices must be 2D".to_string()));
    }

    let n = vertices.nrows();
    if n < 3 {
        return Ok(false);
    }

    let mut sign = 0i32;

    for i in 0..n {
        let p1 = Point2D::new(vertices[[i, 0]], vertices[[i, 1]]);
        let p2 = Point2D::new(vertices[[(i + 1) % n, 0]], vertices[[(i + 1) % n, 1]]);
        let p3 = Point2D::new(vertices[[(i + 2) % n, 0]], vertices[[(i + 2) % n, 1]]);

        let v1 = Point2D::new(p2.x - p1.x, p2.y - p1.y);
        let v2 = Point2D::new(p3.x - p2.x, p3.y - p2.y);

        let cross = v1.cross_product(&v2);

        if cross.abs() > 1e-10 {
            let current_sign = if cross > 0.0 { 1 } else { -1 };

            if sign == 0 {
                sign = current_sign;
            } else if sign != current_sign {
                return Ok(false);
            }
        }
    }

    Ok(true)
}

/// Compute the area of a polygon
#[allow(dead_code)]
pub fn compute_polygon_area(vertices: &ArrayView2<'_, f64>) -> SpatialResult<f64> {
    let polygon = LabeledPolygon::from_array(vertices)?;
    Ok(polygon.compute_area())
}

/// Check if a polygon is self-intersecting
#[allow(dead_code)]
pub fn is_self_intersecting(vertices: &ArrayView2<'_, f64>) -> SpatialResult<bool> {
    let polygon = LabeledPolygon::from_array(vertices)?;
    let n = polygon.vertices.len();

    for i in 0..n {
        let edge1_start = polygon.vertices[i];
        let edge1_end = polygon.vertices[(i + 1) % n];

        for j in (i + 2)..n {
            // Skip adjacent edges
            if j == (i + n - 1) % n {
                continue;
            }

            let edge2_start = polygon.vertices[j];
            let edge2_end = polygon.vertices[(j + 1) % n];

            if line_segment_intersection(&edge1_start, &edge1_end, &edge2_start, &edge2_end)
                .is_some()
            {
                return Ok(true);
            }
        }
    }

    Ok(false)
}

#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_relative_eq;
    use scirs2_core::ndarray::arr2;

    #[test]
    fn testpoint2d_operations() {
        let p1 = Point2D::new(0.0, 0.0);
        let p2 = Point2D::new(1.0, 1.0);

        assert_relative_eq!(p1.distance_to(&p2), 2.0_f64.sqrt(), epsilon = 1e-10);

        let v1 = Point2D::new(1.0, 0.0);
        let v2 = Point2D::new(0.0, 1.0);
        assert_relative_eq!(v1.cross_product(&v2), 1.0, epsilon = 1e-10);
        assert_relative_eq!(v1.dot_product(&v2), 0.0, epsilon = 1e-10);
    }

    #[test]
    fn test_polygon_creation() {
        let vertices = arr2(&[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]);
        let polygon = LabeledPolygon::from_array(&vertices.view()).unwrap();

        assert_eq!(polygon.vertices.len(), 4);
        assert_eq!(polygon.edges.len(), 4);
        assert_relative_eq!(polygon.compute_area(), 1.0, epsilon = 1e-10);
    }

    #[test]
    fn testpoint_inside_polygon() {
        let vertices = arr2(&[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]);
        let polygon = LabeledPolygon::from_array(&vertices.view()).unwrap();

        let insidepoint = Point2D::new(0.5, 0.5);
        let outsidepoint = Point2D::new(1.5, 1.5);

        assert!(polygon.ispoint_inside(&insidepoint));
        assert!(!polygon.ispoint_inside(&outsidepoint));
    }

    #[test]
    fn test_line_intersection() {
        let p1 = Point2D::new(0.0, 0.0);
        let p2 = Point2D::new(1.0, 1.0);
        let p3 = Point2D::new(0.0, 1.0);
        let p4 = Point2D::new(1.0, 0.0);

        let result = line_segment_intersection(&p1, &p2, &p3, &p4);
        assert!(result.is_some());

        let (intersection, t1, t2) = result.unwrap();
        assert_relative_eq!(intersection.x, 0.5, epsilon = 1e-10);
        assert_relative_eq!(intersection.y, 0.5, epsilon = 1e-10);
        assert_relative_eq!(t1, 0.5, epsilon = 1e-10);
        assert_relative_eq!(t2, 0.5, epsilon = 1e-10);
    }

    #[test]
    fn test_is_convex_polygon() {
        // Convex square
        let convex = arr2(&[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]);
        assert!(is_convex_polygon(&convex.view()).unwrap());

        // Non-convex polygon (L-shape)
        let non_convex = arr2(&[
            [0.0, 0.0],
            [1.0, 0.0],
            [1.0, 0.5],
            [0.5, 0.5],
            [0.5, 1.0],
            [0.0, 1.0],
        ]);
        assert!(!is_convex_polygon(&non_convex.view()).unwrap());
    }

    #[test]
    fn test_polygon_area() {
        let square = arr2(&[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]);
        let area = compute_polygon_area(&square.view()).unwrap();
        assert_relative_eq!(area, 1.0, epsilon = 1e-10);

        let triangle = arr2(&[[0.0, 0.0], [1.0, 0.0], [0.5, 1.0]]);
        let area = compute_polygon_area(&triangle.view()).unwrap();
        assert_relative_eq!(area, 0.5, epsilon = 1e-10);
    }

    #[test]
    fn test_intersection() {
        // Non-self-intersecting square
        let square = arr2(&[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]);
        assert!(!is_self_intersecting(&square.view()).unwrap());

        // Self-intersecting bowtie
        let bowtie = arr2(&[[0.0, 0.0], [1.0, 1.0], [1.0, 0.0], [0.0, 1.0]]);
        assert!(is_self_intersecting(&bowtie.view()).unwrap());
    }

    #[test]
    fn test_polygon_union_basic() {
        // Two non-overlapping squares
        let poly1 = arr2(&[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]);
        let poly2 = arr2(&[[2.0, 0.0], [3.0, 0.0], [3.0, 1.0], [2.0, 1.0]]);

        let union = polygon_union(&poly1.view(), &poly2.view()).unwrap();
        assert!(union.nrows() >= 4); // Should have at least 4 vertices
    }

    #[test]
    fn test_polygon_intersection_basic() {
        // Two overlapping squares
        let poly1 = arr2(&[[0.0, 0.0], [2.0, 0.0], [2.0, 2.0], [0.0, 2.0]]);
        let poly2 = arr2(&[[1.0, 1.0], [3.0, 1.0], [3.0, 3.0], [1.0, 3.0]]);

        let intersection = polygon_intersection(&poly1.view(), &poly2.view()).unwrap();

        // The intersection should be a unit square
        if intersection.nrows() > 0 {
            let area = compute_polygon_area(&intersection.view()).unwrap();
            assert!(area > 0.0); // Should have non-zero area
        }
    }

    #[test]
    fn test_polygon_difference_basic() {
        let poly1 = arr2(&[[0.0, 0.0], [2.0, 0.0], [2.0, 2.0], [0.0, 2.0]]);
        let poly2 = arr2(&[[1.0, 1.0], [3.0, 1.0], [3.0, 3.0], [1.0, 3.0]]);

        let difference = polygon_difference(&poly1.view(), &poly2.view()).unwrap();
        assert!(difference.nrows() >= 3); // Should have at least 3 vertices
    }

    #[test]
    fn test_sutherland_hodgman_clip() {
        let subjectvertices = vec![
            Point2D::new(0.0, 0.0),
            Point2D::new(2.0, 0.0),
            Point2D::new(2.0, 2.0),
            Point2D::new(0.0, 2.0),
        ];

        let clipvertices = vec![
            Point2D::new(1.0, 1.0),
            Point2D::new(3.0, 1.0),
            Point2D::new(3.0, 3.0),
            Point2D::new(1.0, 3.0),
        ];

        let subject = LabeledPolygon {
            vertices: subjectvertices,
            edges: Vec::new(),
            is_hole: false,
        };

        let clip = LabeledPolygon {
            vertices: clipvertices,
            edges: Vec::new(),
            is_hole: false,
        };

        let result = sutherland_hodgman_clip(&subject, &clip).unwrap();

        // Should produce a valid polygon
        assert!(result.vertices.len() >= 3);
    }
}