open-standards/map/map-rs/crates/map-index/src/polygon.rs

// S2Polygon — polygon on the sphere (outer boundary + holes).
//
// Ported from joule-s2 with two corrected algorithms and an S2Region impl:
//   * containment now uses an orientation-robust spherical winding-angle test
//     (the upstream ray test cast its arc 1e-10 away from the query point, so it
//     crossed nothing and reported every interior point as "outside");
//   * area now uses the van Oosterom–Strackee signed-area fan (the upstream
//     Girard formulation returned the complement, ~4π − A);
//   * `impl S2Region for S2Polygon` lets the region coverer (SPEC §4) cover an
//     arbitrary polygon — a boundary, an airspace footprint, a geofence — not
//     just caps and rectangles.

use crate::cell_id::S2CellId;
use crate::point::S2Point;
use crate::region::{S2Cap, S2Region};
use serde::{Deserialize, Serialize};

/// A polygon on the unit sphere, represented as one or more loops.
/// The first loop is the outer boundary; subsequent loops are holes.
#[derive(Clone, Debug, Serialize, Deserialize)]
pub struct S2Polygon {
    /// Each loop is a sequence of vertices. The last vertex connects back to the first.
    pub loops: Vec<Vec<S2Point>>,
}

impl S2Polygon {
    /// Create a polygon from loops.
    pub fn new(loops: Vec<Vec<S2Point>>) -> Self {
        Self { loops }
    }

    /// Create a simple polygon (no holes) from a single loop.
    pub fn from_loop(vertices: Vec<S2Point>) -> Self {
        Self { loops: vec![vertices] }
    }

    /// Test whether a point is inside the polygon. The point is inside iff it is
    /// inside the outer boundary and inside no hole.
    pub fn contains_point(&self, point: &S2Point) -> bool {
        if self.loops.is_empty() {
            return false;
        }
        if !loop_contains_point(&self.loops[0], point) {
            return false;
        }
        for hole in &self.loops[1..] {
            if loop_contains_point(hole, point) {
                return false;
            }
        }
        true
    }

    /// Area of the polygon in steradians (outer loop minus holes).
    pub fn get_area(&self) -> f64 {
        if self.loops.is_empty() {
            return 0.0;
        }
        let mut total = loop_signed_area(&self.loops[0]).abs();
        for hole in &self.loops[1..] {
            total -= loop_signed_area(hole).abs();
        }
        total.max(0.0)
    }

    /// Area-weighted centroid (planar triangulation approximation).
    pub fn get_centroid(&self) -> S2Point {
        if self.loops.is_empty() {
            return S2Point::new(0.0, 0.0, 0.0);
        }
        let (mut cx, mut cy, mut cz, mut total_weight) = (0.0, 0.0, 0.0, 0.0);
        for lp in &self.loops {
            if lp.len() < 3 {
                continue;
            }
            for i in 1..lp.len() - 1 {
                let (a, b, c) = (lp[0], lp[i], lp[i + 1]);
                let area = b.add(a.neg()).cross(c.add(a.neg())).norm() * 0.5;
                cx += (a.x + b.x + c.x) / 3.0 * area;
                cy += (a.y + b.y + c.y) / 3.0 * area;
                cz += (a.z + b.z + c.z) / 3.0 * area;
                total_weight += area;
            }
        }
        if total_weight < 1e-30 {
            return S2Point::new(0.0, 0.0, 0.0);
        }
        S2Point::new(cx / total_weight, cy / total_weight, cz / total_weight).normalize()
    }

    /// The minimal-ish bounding cap over the outer loop: centred on the mean
    /// direction of the outer vertices, with radius the largest angle to any of
    /// them. Used as a provably-safe superset for `may_intersect_cell` — a cell
    /// the bounding cap misses cannot intersect the polygon.
    pub fn bounding_cap(&self) -> Option<S2Cap> {
        let outer = self.loops.first()?;
        if outer.len() < 3 {
            return None;
        }
        let (mut sx, mut sy, mut sz) = (0.0, 0.0, 0.0);
        for v in outer {
            sx += v.x;
            sy += v.y;
            sz += v.z;
        }
        let center = S2Point::new(sx, sy, sz).normalize();
        let radius = outer.iter().fold(0.0_f64, |r, v| r.max(center.angle(*v)));
        Some(S2Cap::new(center, radius))
    }

    /// Number of vertices across all loops.
    pub fn num_vertices(&self) -> usize {
        self.loops.iter().map(|l| l.len()).sum()
    }

    /// Number of loops.
    pub fn num_loops(&self) -> usize {
        self.loops.len()
    }
}

impl S2Region for S2Polygon {
    fn contains_point(&self, p: &S2Point) -> bool {
        S2Polygon::contains_point(self, p)
    }

    /// Conservative (never a false negative): the polygon can only intersect a
    /// cell its bounding cap intersects.
    fn may_intersect_cell(&self, cell: S2CellId) -> bool {
        match self.bounding_cap() {
            Some(cap) => cap.may_intersect_cell(cell),
            None => false,
        }
    }
}

/// The tangent direction at `p` toward `v` (the component of `v` orthogonal to
/// `p`, normalized). `None` if `v` coincides with `±p`.
fn tangent_toward(p: S2Point, v: S2Point) -> Option<S2Point> {
    let t = v.add(p.scale(-v.dot(p)));
    if t.norm() < 1e-12 {
        None
    } else {
        Some(t.normalize())
    }
}

/// Orientation-robust spherical point-in-loop: sum the signed angles subtended
/// at `point` by consecutive edges. The total is ±2π inside and ≈0 outside,
/// regardless of winding direction — so `|total| > π` ⇒ inside.
fn loop_contains_point(vertices: &[S2Point], point: &S2Point) -> bool {
    let n = vertices.len();
    if n < 3 {
        return false;
    }
    let p = *point;
    let mut total = 0.0;
    for i in 0..n {
        let (ta, tb) = match (tangent_toward(p, vertices[i]), tangent_toward(p, vertices[(i + 1) % n])) {
            (Some(a), Some(b)) => (a, b),
            // The query point coincides with a vertex — treat as on-boundary/inside.
            _ => return true,
        };
        let dot = ta.dot(tb).clamp(-1.0, 1.0);
        let signed = ta.cross(tb).dot(p);
        total += signed.atan2(dot);
    }
    total.abs() > std::f64::consts::PI
}

/// Signed area of a loop in steradians via the van Oosterom–Strackee fan from
/// vertex 0. Sign encodes orientation (positive = counterclockwise); signed
/// triangle contributions cancel correctly for non-convex loops, as the planar
/// shoelace does.
fn loop_signed_area(vertices: &[S2Point]) -> f64 {
    let n = vertices.len();
    if n < 3 {
        return 0.0;
    }
    let a = vertices[0];
    let mut sum = 0.0;
    for i in 1..n - 1 {
        let b = vertices[i];
        let c = vertices[i + 1];
        let num = a.dot(b.cross(c));
        let den = 1.0 + a.dot(b) + b.dot(c) + c.dot(a);
        sum += 2.0 * num.atan2(den);
    }
    sum
}

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

    fn make_square_loop(half_size: f64) -> Vec<S2Point> {
        vec![
            S2Point::from_lat_lng(half_size, half_size),
            S2Point::from_lat_lng(half_size, -half_size),
            S2Point::from_lat_lng(-half_size, -half_size),
            S2Point::from_lat_lng(-half_size, half_size),
        ]
    }

    #[test]
    fn simple_polygon_contains_center() {
        let poly = S2Polygon::from_loop(make_square_loop(0.01));
        assert!(poly.contains_point(&S2Point::from_lat_lng(0.0, 0.0)));
    }

    #[test]
    fn simple_polygon_excludes_far_point() {
        let poly = S2Polygon::from_loop(make_square_loop(0.01));
        assert!(!poly.contains_point(&S2Point::from_lat_lng(1.0, 1.0)));
    }

    #[test]
    fn polygon_area_positive_and_reasonable() {
        let poly = S2Polygon::from_loop(make_square_loop(0.1));
        let area = poly.get_area();
        // A ~0.2 × 0.2 rad square ≈ 0.04 sr.
        assert!(area > 0.0, "area should be positive, got {area}");
        assert!(area < 1.0, "area should be reasonable, got {area}");
        assert!((area - 0.04).abs() < 0.01, "area ≈ 0.04 sr, got {area}");
    }

    #[test]
    fn winding_is_orientation_independent() {
        // Same square, reversed winding — containment must agree.
        let mut rev = make_square_loop(0.05);
        rev.reverse();
        let poly = S2Polygon::from_loop(rev);
        assert!(poly.contains_point(&S2Point::from_lat_lng(0.0, 0.0)));
        assert!(!poly.contains_point(&S2Point::from_lat_lng(0.5, 0.5)));
    }

    #[test]
    fn polygon_with_hole() {
        let outer = make_square_loop(0.1);
        let inner = vec![
            S2Point::from_lat_lng(0.02, 0.02),
            S2Point::from_lat_lng(0.02, -0.02),
            S2Point::from_lat_lng(-0.02, -0.02),
            S2Point::from_lat_lng(-0.02, 0.02),
        ];
        let poly = S2Polygon::new(vec![outer, inner]);
        assert!(!poly.contains_point(&S2Point::from_lat_lng(0.0, 0.0)), "center is in the hole");
        assert!(poly.contains_point(&S2Point::from_lat_lng(0.05, 0.05)), "between boundaries is inside");
    }

    #[test]
    fn bounding_cap_encloses_vertices() {
        let lp = make_square_loop(0.1);
        let poly = S2Polygon::from_loop(lp.clone());
        let cap = poly.bounding_cap().unwrap();
        for v in &lp {
            assert!(cap.contains_point(v), "cap must enclose every outer vertex");
        }
    }

    #[test]
    fn centroid_near_center() {
        let poly = S2Polygon::from_loop(make_square_loop(0.01));
        let (lat, lng) = poly.get_centroid().to_lat_lng();
        assert!(lat.abs() < 0.01 && lng.abs() < 0.01);
    }

    #[test]
    fn num_vertices_and_empty() {
        let poly = S2Polygon::from_loop(make_square_loop(0.01));
        assert_eq!(poly.num_vertices(), 4);
        assert_eq!(poly.num_loops(), 1);
        let empty = S2Polygon::new(vec![]);
        assert!(!empty.contains_point(&S2Point::new(1.0, 0.0, 0.0)));
        assert!((empty.get_area() - 0.0).abs() < 1e-14);
    }
}