geo-polygonize-core 0.51.2

use geo::Coord;
#[cfg(not(target_arch = "wasm32"))]
use multiversion::multiversion;
use wide::f64x4;
use wide::CmpGt;

#[cfg(all(target_arch = "aarch64", target_os = "linux"))]
const SCALAR_MIN_COORDS: usize = 1;
#[cfg(not(all(target_arch = "aarch64", target_os = "linux")))]
const SCALAR_MIN_COORDS: usize = 257;

pub struct SimdRing {
    pub x: Vec<f64>,
    pub y: Vec<f64>,
    len: usize,
}

impl SimdRing {
    pub fn new(coords: &[Coord<f64>]) -> Self {
        let len = coords.len();
        let capacity = (len + 3) & !3;

        let mut x: Vec<f64> = Vec::with_capacity(capacity);
        let mut y: Vec<f64> = Vec::with_capacity(capacity);

        x.extend(coords.iter().map(|c| c.x));
        y.extend(coords.iter().map(|c| c.y));

        let last_x = x.last().copied().unwrap_or(0.0);
        let last_y = y.last().copied().unwrap_or(0.0);

        x.resize(capacity, last_x);
        y.resize(capacity, last_y);

        Self { x, y, len }
    }

    pub fn new_3d(coords: &[crate::types::Coord3D]) -> Self {
        let len = coords.len();
        let capacity = (len + 3) & !3;

        let mut x: Vec<f64> = Vec::with_capacity(capacity);
        let mut y: Vec<f64> = Vec::with_capacity(capacity);

        x.extend(coords.iter().map(|c| c.x));
        y.extend(coords.iter().map(|c| c.y));

        let last_x = x.last().copied().unwrap_or(0.0);
        let last_y = y.last().copied().unwrap_or(0.0);

        x.resize(capacity, last_x);
        y.resize(capacity, last_y);

        Self { x, y, len }
    }

    pub fn contains(&self, point: Coord<f64>) -> bool {
        if self.len == 0 {
            return false;
        }

        if self.len >= SCALAR_MIN_COORDS {
            contains_scalar(&self.x, &self.y, self.len, point)
        } else {
            contains_simd(&self.x, &self.y, self.len, point)
        }
    }
}

fn contains_scalar(x: &[f64], y: &[f64], len: usize, point: Coord<f64>) -> bool {
    let mut crossings = 0;
    for i in 0..len - 1 {
        if ((y[i] > point.y) != (y[i + 1] > point.y))
            && point.x < (x[i + 1] - x[i]) * (point.y - y[i]) / (y[i + 1] - y[i]) + x[i]
        {
            crossings += 1;
        }
    }
    crossings % 2 != 0
}

#[cfg_attr(
    not(target_arch = "wasm32"),
    multiversion(targets(
        "x86_64+avx512f+avx512dq",
        "x86_64+avx2",
        "x86+avx2",
        "x86_64+avx",
        "x86+avx",
        "x86_64+sse2",
        "x86+sse2",
    ))
)]
fn contains_simd(x: &[f64], y: &[f64], len: usize, point: Coord<f64>) -> bool {
    let px = f64x4::splat(point.x);
    let py = f64x4::splat(point.y);

    let n = len - 1; // Number of segments

    let mut i = 0;
    let mut crossings = 0;

    while i < n {
        let remaining = n - i;
        if remaining >= 4 {
            let xi = f64x4::from(&x[i..i + 4]);
            let yi = f64x4::from(&y[i..i + 4]);

            let xj = f64x4::from(&x[i + 1..i + 5]);
            let yj = f64x4::from(&y[i + 1..i + 5]);

            let yi_gt_py = yi.cmp_gt(py);
            let yj_gt_py = yj.cmp_gt(py);
            let in_range = yi_gt_py ^ yj_gt_py;

            let num = (xj - xi) * (py - yi);
            let den = yj - yi;

            let intersect_x = (num / den) + xi;
            let x_cond = intersect_x.cmp_gt(px);

            let is_crossing = in_range & x_cond;

            crossings += is_crossing.move_mask().count_ones();

            i += 4;
        } else {
            let p1x = x[i];
            let p1y = y[i];
            let p2x = x[i + 1];
            let p2y = y[i + 1];

            if ((p1y > point.y) != (p2y > point.y))
                && (point.x < (p2x - p1x) * (point.y - p1y) / (p2y - p1y) + p1x)
            {
                crossings += 1;
            }
            i += 1;
        }
    }

    crossings % 2 != 0
}

#[cfg(test)]
mod tests {
    use super::*;
    use geo::Contains;
    use geo::{Coord, LineString, Polygon};
    use rand::rngs::StdRng;
    use rand::{Rng, SeedableRng};

    #[test]
    fn test_simd_ring_square() {
        let coords = vec![
            Coord { x: 0.0, y: 0.0 },
            Coord { x: 10.0, y: 0.0 },
            Coord { x: 10.0, y: 10.0 },
            Coord { x: 0.0, y: 10.0 },
            Coord { x: 0.0, y: 0.0 },
        ];
        let ring = SimdRing::new(&coords);

        // Inside
        assert!(ring.contains(Coord { x: 5.0, y: 5.0 }));

        // Outside
        assert!(!ring.contains(Coord { x: 15.0, y: 5.0 }));
        assert!(!ring.contains(Coord { x: 5.0, y: 15.0 }));
        assert!(!ring.contains(Coord { x: -5.0, y: 5.0 }));
        assert!(!ring.contains(Coord { x: 5.0, y: -5.0 }));
    }

    #[test]
    fn test_simd_ring_triangle() {
        let coords = vec![
            Coord { x: 0.0, y: 0.0 },
            Coord { x: 10.0, y: 0.0 },
            Coord { x: 0.0, y: 10.0 },
            Coord { x: 0.0, y: 0.0 },
        ];
        let ring = SimdRing::new(&coords);

        // Inside
        assert!(ring.contains(Coord { x: 2.0, y: 2.0 }));
        // Outside
        assert!(!ring.contains(Coord { x: 10.0, y: 10.0 })); // Outside bounding box corner
        assert!(!ring.contains(Coord { x: 5.0, y: 6.0 })); // Outside hypotenuse
    }

    #[test]
    fn test_simd_ring_complex() {
        // A "U" shape polygon
        // 0,0 -> 10,0 -> 10,10 -> 8,10 -> 8,2 -> 2,2 -> 2,10 -> 0,10 -> 0,0
        let coords = vec![
            Coord { x: 0.0, y: 0.0 },
            Coord { x: 10.0, y: 0.0 },
            Coord { x: 10.0, y: 10.0 },
            Coord { x: 8.0, y: 10.0 },
            Coord { x: 8.0, y: 2.0 },
            Coord { x: 2.0, y: 2.0 },
            Coord { x: 2.0, y: 10.0 },
            Coord { x: 0.0, y: 10.0 },
            Coord { x: 0.0, y: 0.0 },
        ];
        let ring = SimdRing::new(&coords);

        // Inside the U base
        assert!(ring.contains(Coord { x: 5.0, y: 1.0 }));
        // Inside the U arms
        assert!(ring.contains(Coord { x: 1.0, y: 5.0 }));
        assert!(ring.contains(Coord { x: 9.0, y: 5.0 }));

        // Inside the "hole" of the U (concave part)
        // x=5, y=5 is in the empty space between arms
        assert!(!ring.contains(Coord { x: 5.0, y: 5.0 }));

        // Outside bounding box
        assert!(!ring.contains(Coord { x: 12.0, y: 5.0 }));
    }

    #[test]
    fn test_simd_ring_empty() {
        let coords: Vec<Coord<f64>> = vec![];
        let ring = SimdRing::new(&coords);
        assert!(!ring.contains(Coord { x: 0.0, y: 0.0 }));

        let coords_single = vec![Coord { x: 0.0, y: 0.0 }];
        let ring_single = SimdRing::new(&coords_single);
        assert!(!ring_single.contains(Coord { x: 0.0, y: 0.0 }));
    }

    #[test]
    fn test_simd_ring_alignment_fallback() {
        // Test varying sizes to exercise SIMD + scalar tail logic
        // A square 0,0 to 10,10.
        // 5 points: 0,0 -> 10,0 -> 10,10 -> 0,10 -> 0,0.
        // len=5. n=4. 1 iteration of SIMD (4). No scalar tail.

        // Let's create a polygon with more points on the edges to force scalar tail.
        // 0,0 -> 5,0 -> 10,0 -> 10,5 -> 10,10 -> 5,10 -> 0,10 -> 0,5 -> 0,0.
        // 9 points. len=9. n=8.
        // 2 iterations of SIMD. No scalar tail.

        // Add one more point: 0,0 -> ... -> 0,2 -> 0,0.
        // 10 points. len=10. n=9.
        // 2 iterations of SIMD (8). 1 scalar tail (9th segment).

        let coords = vec![
            Coord { x: 0.0, y: 0.0 },
            Coord { x: 5.0, y: 0.0 },
            Coord { x: 10.0, y: 0.0 },
            Coord { x: 10.0, y: 5.0 },
            Coord { x: 10.0, y: 10.0 },
            Coord { x: 5.0, y: 10.0 },
            Coord { x: 0.0, y: 10.0 },
            Coord { x: 0.0, y: 5.0 },
            Coord { x: 0.0, y: 2.0 },
            Coord { x: 0.0, y: 0.0 },
        ];
        let ring = SimdRing::new(&coords);

        // Inside
        assert!(ring.contains(Coord { x: 5.0, y: 5.0 }));
        // Outside
        assert!(!ring.contains(Coord { x: 15.0, y: 5.0 }));
    }

    #[test]
    fn test_simd_ring_large_coords() {
        // Test with large coordinates to check numerical stability
        let offset = 1_000_000.0;
        let coords = vec![
            Coord {
                x: offset,
                y: offset,
            },
            Coord {
                x: offset + 10.0,
                y: offset,
            },
            Coord {
                x: offset + 10.0,
                y: offset + 10.0,
            },
            Coord {
                x: offset,
                y: offset + 10.0,
            },
            Coord {
                x: offset,
                y: offset,
            },
        ];
        let ring = SimdRing::new(&coords);

        assert!(ring.contains(Coord {
            x: offset + 5.0,
            y: offset + 5.0
        }));
        assert!(!ring.contains(Coord {
            x: offset + 15.0,
            y: offset + 5.0
        }));
    }

    #[test]
    fn test_simd_ring_vs_geo_random() {
        let mut rng = StdRng::seed_from_u64(42);

        for _ in 0..100 {
            // Generate a random triangle
            let p1 = Coord {
                x: rng.gen_range(0.0..100.0),
                y: rng.gen_range(0.0..100.0),
            };
            let p2 = Coord {
                x: rng.gen_range(0.0..100.0),
                y: rng.gen_range(0.0..100.0),
            };
            let p3 = Coord {
                x: rng.gen_range(0.0..100.0),
                y: rng.gen_range(0.0..100.0),
            };

            let coords = vec![p1, p2, p3, p1];
            let ring = SimdRing::new(&coords);
            let poly = Polygon::new(LineString::new(coords), vec![]);

            // Test random points
            for _ in 0..10 {
                let test_point = Coord {
                    x: rng.gen_range(0.0..100.0),
                    y: rng.gen_range(0.0..100.0),
                };

                // Geo's contains is strict (false for boundary).
                // SimdRing is strict-ish (false for most boundary, true for some vertices).
                // We skip points too close to boundary to avoid flaky tests on undefined behavior.
                // Convert test_point to Point for distance check
                // let point_geo = geo::Point(test_point);
                // Check distance to boundary
                let boundary = poly.exterior();

                let mut min_dist = f64::MAX;
                for i in 0..boundary.0.len() - 1 {
                    let a = boundary.0[i];
                    let b = boundary.0[i + 1];
                    let p = test_point;

                    let l2 = (b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y);
                    let d = if l2 == 0.0 {
                        (p.x - a.x) * (p.x - a.x) + (p.y - a.y) * (p.y - a.y)
                    } else {
                        let t = ((p.x - a.x) * (b.x - a.x) + (p.y - a.y) * (b.y - a.y)) / l2;
                        let t = t.clamp(0.0, 1.0);
                        let proj_x = a.x + t * (b.x - a.x);
                        let proj_y = a.y + t * (b.y - a.y);
                        (p.x - proj_x) * (p.x - proj_x) + (p.y - proj_y) * (p.y - proj_y)
                    };

                    if d < min_dist {
                        min_dist = d;
                    }
                }

                // Using Euclidean distance squared avoiding expensive `.sqrt()` in the hot loop
                if min_dist < 1e-10 {
                    // 1e-5 squared
                    continue;
                }

                let simd_contains = ring.contains(test_point);
                let geo_contains = poly.contains(&test_point);

                assert_eq!(
                    simd_contains, geo_contains,
                    "Mismatch for point {:?} in triangle {:?}, {:?}, {:?}",
                    test_point, p1, p2, p3
                );
            }
        }
    }

    #[test]
    fn adaptive_scalar_matches_simd() {
        let mut coords: Vec<_> = (0..256)
            .map(|i| {
                let angle = std::f64::consts::TAU * i as f64 / 256.0;
                Coord {
                    x: angle.cos() * 10.0,
                    y: angle.sin() * 10.0,
                }
            })
            .collect();
        coords.push(coords[0]);
        let ring = SimdRing::new(&coords);

        for point in [
            Coord { x: 0.0, y: 0.0 },
            Coord { x: 9.0, y: 0.0 },
            Coord { x: 11.0, y: 0.0 },
        ] {
            assert_eq!(
                contains_scalar(&ring.x, &ring.y, ring.len, point),
                contains_simd(&ring.x, &ring.y, ring.len, point)
            );
            assert_eq!(
                ring.contains(point),
                contains_scalar(&ring.x, &ring.y, ring.len, point)
            );
        }
    }

    #[test]
    fn test_simd_ring_boundary_documented_behavior() {
        // Documenting the specific behavior of the current implementation for boundary points.
        // Logic: Ray casting to +infinity X.
        // Vertical edges: Ignored if x != px. If x == px, ray runs along edge (degenerate).
        // If ray passes through a vertex:
        //  - If edges go Up-Up or Down-Down: 2 crossings (even) -> Outside.
        //  - If edges go Up-Down or Down-Up: 1 crossing (odd) -> Inside.
        //  - Vertices with y == py are treated as "below" the ray (strict > check).

        let coords = vec![
            Coord { x: 0.0, y: 0.0 },
            Coord { x: 10.0, y: 0.0 },
            Coord { x: 10.0, y: 10.0 },
            Coord { x: 0.0, y: 10.0 },
            Coord { x: 0.0, y: 0.0 },
        ];
        let ring = SimdRing::new(&coords);

        // (0,0): Vertex.
        // Edge (0,10)->(0,0) comes down to y=0. (10 > 0 is True).
        // Edge (0,0)->(10,0) goes right at y=0. (0 > 0 is False).
        // Crossing check: `true != false` -> 1 crossing. Inside.
        assert!(ring.contains(Coord { x: 0.0, y: 0.0 }));

        // (10,0): Vertex.
        // Edge (0,0)->(10,0) comes to y=0.
        // Edge (10,0)->(10,10) goes up from y=0.
        // Crossing check: `false != true`. 1 crossing.
        // BUT point.x < intersect_x.
        // Ray starts at x=10. Intersect is at x=10.
        // 10 < 10 is False. No crossing counted. Outside.
        assert!(!ring.contains(Coord { x: 10.0, y: 0.0 }));

        // (10,10): Vertex.
        // Edge (10,0)->(10,10) comes up to y=10. (10 > 10 is False).
        // Edge (10,10)->(0,10) goes left at y=10. (10 > 10 is False).
        // Crossing check: `false != false`. 0 crossings. Outside.
        assert!(!ring.contains(Coord { x: 10.0, y: 10.0 }));

        // (0,10): Vertex.
        // Edge (10,10)->(0,10) comes to y=10.
        // Edge (0,10)->(0,0) goes down from y=10.
        // Crossing check: `false != true`. 1 crossing.
        // intersect_x = 0. 0 < 0 False. Outside.
        assert!(!ring.contains(Coord { x: 0.0, y: 10.0 }));

        // Midpoint of bottom edge (5,0).
        // Ray hits (10,0)-(10,10) at x=10. 5 < 10. 1 crossing. Inside.
        assert!(ring.contains(Coord { x: 5.0, y: 0.0 }));

        // Midpoint of top edge (5,10).
        // Ray hits nothing to the right (except parallel edge (10,10)-(0,10) which is ignored).
        // Outside.
        assert!(!ring.contains(Coord { x: 5.0, y: 10.0 }));

        // Midpoint of left edge (0,5).
        // Ray hits right edge at x=10. 0 < 10. Inside.
        assert!(ring.contains(Coord { x: 0.0, y: 5.0 }));

        // Midpoint of right edge (10,5).
        // Ray starts at x=10. Hits right edge at x=10. 10 < 10 False. Outside.
        assert!(!ring.contains(Coord { x: 10.0, y: 5.0 }));
    }
}