tg_geom/

geom.rs

//! The main geometry type supporting all OGC Simple Feature geometry types.
//!
//! [`Geom`] is the primary type for working with geometries. It can represent
//! Points, `LineStrings`, Polygons, and their Multi/Collection variants. It
//! supports parsing from WKT, WKB, `GeoJSON`, and other formats, as well as
//! spatial predicates (intersects, contains, etc.).
//!
//! # Coordinate Dimensions
//!
//! Geometries can have 2 to 4 dimensions:
//! - **XY** (2D): Basic 2D coordinates
//! - **XYZ** (3D): Includes Z coordinate for elevation/altitude
//! - **XYM** (3D): Includes M coordinate for linear referencing (distance along
//!   a route)
//! - **XYZM** (4D): Includes both Z and M coordinates

use core::ffi::CStr;
use core::marker::PhantomData;
use core::ptr::NonNull;

use crate::error::{Error, Result};
use crate::{GeomType, Line, LineRef, Point, Poly, PolyRef, Rect};

macro_rules! impl_geom_methods {
    ($get_ptr:expr) => {
        #[inline]
        pub fn geom_type(&self) -> GeomType {
            let ptr = $get_ptr(self);
            let t = unsafe { tg_geom_sys::tg_geom_typeof(ptr) };
            t.into()
        }

        #[inline]
        pub fn rect(&self) -> Rect {
            let ptr = $get_ptr(self);
            let r = unsafe { tg_geom_sys::tg_geom_rect(ptr) };
            r.into()
        }

        #[inline]
        pub fn is_empty(&self) -> bool {
            let ptr = $get_ptr(self);
            unsafe { tg_geom_sys::tg_geom_is_empty(ptr) }
        }
    };
}

/// An owned geometry that can represent any OGC Simple Feature type.
///
/// This is the main type for working with geometries. It supports:
/// - All geometry types: Point, `LineString`, Polygon, and Multi/Collection
///   variants
/// - Parsing from WKT, WKB, `GeoJSON`, HEX, and `GeoBin` formats
/// - Spatial predicates: intersects, contains, within, covers, touches, etc.
/// - Serialization to various output formats
///
/// Memory is managed by the TG library and freed when dropped.
pub struct Geom {
    ptr: NonNull<tg_geom_sys::tg_geom>,
}

impl Geom {
    impl_geom_methods!(|s: &Self| s.ptr.as_ptr());

    pub fn new_point(point: Point) -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_new_point(point.into()) };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_point_z(point: Point, z: f64) -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_new_point_z(point.into(), z) };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_point_m(point: Point, m: f64) -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_new_point_m(point.into(), m) };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_point_zm(point: Point, z: f64, m: f64) -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_new_point_zm(point.into(), z, m) };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_point_empty() -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_new_point_empty() };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_linestring(line: &Line) -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_new_linestring(line.as_ptr()) };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_linestring_z(line: &Line, extra_coords: &[f64]) -> Result<Self> {
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_linestring_z(
                line.as_ptr(),
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_linestring_m(line: &Line, extra_coords: &[f64]) -> Result<Self> {
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_linestring_m(
                line.as_ptr(),
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_linestring_zm(line: &Line, extra_coords: &[f64]) -> Result<Self> {
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_linestring_zm(
                line.as_ptr(),
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_linestring_empty() -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_new_linestring_empty() };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_polygon(poly: &Poly) -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_new_polygon(poly.as_ptr()) };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_polygon_z(poly: &Poly, extra_coords: &[f64]) -> Result<Self> {
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_polygon_z(
                poly.as_ptr(),
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_polygon_m(poly: &Poly, extra_coords: &[f64]) -> Result<Self> {
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_polygon_m(
                poly.as_ptr(),
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_polygon_zm(poly: &Poly, extra_coords: &[f64]) -> Result<Self> {
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_polygon_zm(
                poly.as_ptr(),
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_polygon_empty() -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_new_polygon_empty() };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multipoint(points: &[Point]) -> Result<Self> {
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_multipoint(
                points.as_ptr() as *const tg_geom_sys::tg_point,
                points.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multipoint_z(points: &[Point], extra_coords: &[f64]) -> Result<Self> {
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_multipoint_z(
                points.as_ptr() as *const tg_geom_sys::tg_point,
                points.len() as libc::c_int,
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multipoint_m(points: &[Point], extra_coords: &[f64]) -> Result<Self> {
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_multipoint_m(
                points.as_ptr() as *const tg_geom_sys::tg_point,
                points.len() as libc::c_int,
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multipoint_zm(points: &[Point], extra_coords: &[f64]) -> Result<Self> {
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_multipoint_zm(
                points.as_ptr() as *const tg_geom_sys::tg_point,
                points.len() as libc::c_int,
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multipoint_empty() -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_new_multipoint_empty() };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multilinestring(lines: &[&Line]) -> Result<Self> {
        let line_ptrs: Vec<*const tg_geom_sys::tg_line> =
            lines.iter().map(|l| l.as_ptr()).collect();
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_multilinestring(line_ptrs.as_ptr(), lines.len() as libc::c_int)
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multilinestring_z(lines: &[&Line], extra_coords: &[f64]) -> Result<Self> {
        let line_ptrs: Vec<*const tg_geom_sys::tg_line> =
            lines.iter().map(|l| l.as_ptr()).collect();
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_multilinestring_z(
                line_ptrs.as_ptr(),
                lines.len() as libc::c_int,
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multilinestring_m(lines: &[&Line], extra_coords: &[f64]) -> Result<Self> {
        let line_ptrs: Vec<*const tg_geom_sys::tg_line> =
            lines.iter().map(|l| l.as_ptr()).collect();
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_multilinestring_m(
                line_ptrs.as_ptr(),
                lines.len() as libc::c_int,
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multilinestring_zm(lines: &[&Line], extra_coords: &[f64]) -> Result<Self> {
        let line_ptrs: Vec<*const tg_geom_sys::tg_line> =
            lines.iter().map(|l| l.as_ptr()).collect();
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_multilinestring_zm(
                line_ptrs.as_ptr(),
                lines.len() as libc::c_int,
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multilinestring_empty() -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_new_multilinestring_empty() };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multipolygon(polys: &[&Poly]) -> Result<Self> {
        let poly_ptrs: Vec<*const tg_geom_sys::tg_poly> =
            polys.iter().map(|p| p.as_ptr()).collect();
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_multipolygon(poly_ptrs.as_ptr(), polys.len() as libc::c_int)
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multipolygon_z(polys: &[&Poly], extra_coords: &[f64]) -> Result<Self> {
        let poly_ptrs: Vec<*const tg_geom_sys::tg_poly> =
            polys.iter().map(|p| p.as_ptr()).collect();
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_multipolygon_z(
                poly_ptrs.as_ptr(),
                polys.len() as libc::c_int,
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multipolygon_m(polys: &[&Poly], extra_coords: &[f64]) -> Result<Self> {
        let poly_ptrs: Vec<*const tg_geom_sys::tg_poly> =
            polys.iter().map(|p| p.as_ptr()).collect();
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_multipolygon_m(
                poly_ptrs.as_ptr(),
                polys.len() as libc::c_int,
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multipolygon_zm(polys: &[&Poly], extra_coords: &[f64]) -> Result<Self> {
        let poly_ptrs: Vec<*const tg_geom_sys::tg_poly> =
            polys.iter().map(|p| p.as_ptr()).collect();
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_multipolygon_zm(
                poly_ptrs.as_ptr(),
                polys.len() as libc::c_int,
                extra_coords.as_ptr(),
                extra_coords.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_multipolygon_empty() -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_new_multipolygon_empty() };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_geometrycollection(geoms: &[&Geom]) -> Result<Self> {
        let geom_ptrs: Vec<*const tg_geom_sys::tg_geom> =
            geoms.iter().map(|g| g.as_ptr()).collect();
        let ptr = unsafe {
            tg_geom_sys::tg_geom_new_geometrycollection(
                geom_ptrs.as_ptr(),
                geoms.len() as libc::c_int,
            )
        };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    pub fn new_geometrycollection_empty() -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_new_geometrycollection_empty() };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::OutOfMemory)
    }

    /// # Safety
    /// The pointer must be valid and not owned elsewhere.
    pub unsafe fn from_raw(ptr: *mut tg_geom_sys::tg_geom) -> Option<Self> {
        NonNull::new(ptr).map(|ptr| Self { ptr })
    }

    #[inline]
    pub fn as_ptr(&self) -> *const tg_geom_sys::tg_geom {
        self.ptr.as_ptr()
    }

    pub fn into_raw(self) -> *mut tg_geom_sys::tg_geom {
        let ptr = self.ptr.as_ptr();
        core::mem::forget(self);
        ptr
    }

    /// Deep copy.
    pub fn copy(&self) -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_copy(self.ptr.as_ptr()) };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::CopyFailed)
    }

    /// Reference-counted clone sharing the underlying data.
    pub fn clone_ref(&self) -> Result<Self> {
        let ptr = unsafe { tg_geom_sys::tg_geom_clone(self.ptr.as_ptr()) };
        NonNull::new(ptr)
            .map(|ptr| Self { ptr })
            .ok_or(Error::CopyFailed)
    }

    pub fn is_feature(&self) -> bool {
        unsafe { tg_geom_sys::tg_geom_is_feature(self.ptr.as_ptr()) }
    }

    pub fn is_featurecollection(&self) -> bool {
        unsafe { tg_geom_sys::tg_geom_is_featurecollection(self.ptr.as_ptr()) }
    }

    pub fn point(&self) -> Point {
        let p = unsafe { tg_geom_sys::tg_geom_point(self.ptr.as_ptr()) };
        p.into()
    }

    pub fn line(&self) -> Option<LineRef<'_>> {
        let ptr = unsafe { tg_geom_sys::tg_geom_line(self.ptr.as_ptr()) };
        unsafe { LineRef::from_raw(ptr) }
    }

    pub fn poly(&self) -> Option<PolyRef<'_>> {
        let ptr = unsafe { tg_geom_sys::tg_geom_poly(self.ptr.as_ptr()) };
        unsafe { PolyRef::from_raw(ptr) }
    }

    pub fn num_points(&self) -> usize {
        unsafe { tg_geom_sys::tg_geom_num_points(self.ptr.as_ptr()) as usize }
    }

    pub fn point_at(&self, index: usize) -> Option<Point> {
        if index >= self.num_points() {
            None
        } else {
            let p =
                unsafe { tg_geom_sys::tg_geom_point_at(self.ptr.as_ptr(), index as libc::c_int) };
            Some(p.into())
        }
    }

    pub fn num_lines(&self) -> usize {
        unsafe { tg_geom_sys::tg_geom_num_lines(self.ptr.as_ptr()) as usize }
    }

    pub fn line_at(&self, index: usize) -> Option<LineRef<'_>> {
        if index >= self.num_lines() {
            None
        } else {
            let ptr =
                unsafe { tg_geom_sys::tg_geom_line_at(self.ptr.as_ptr(), index as libc::c_int) };
            unsafe { LineRef::from_raw(ptr) }
        }
    }

    pub fn num_polys(&self) -> usize {
        unsafe { tg_geom_sys::tg_geom_num_polys(self.ptr.as_ptr()) as usize }
    }

    pub fn poly_at(&self, index: usize) -> Option<PolyRef<'_>> {
        if index >= self.num_polys() {
            None
        } else {
            let ptr =
                unsafe { tg_geom_sys::tg_geom_poly_at(self.ptr.as_ptr(), index as libc::c_int) };
            unsafe { PolyRef::from_raw(ptr) }
        }
    }

    pub fn num_geometries(&self) -> usize {
        unsafe { tg_geom_sys::tg_geom_num_geometries(self.ptr.as_ptr()) as usize }
    }

    pub fn geometry_at(&self, index: usize) -> Option<GeomRef<'_>> {
        if index >= self.num_geometries() {
            None
        } else {
            let ptr = unsafe {
                tg_geom_sys::tg_geom_geometry_at(self.ptr.as_ptr(), index as libc::c_int)
            };
            unsafe { GeomRef::from_raw(ptr) }
        }
    }

    pub fn extra_json(&self) -> Option<&str> {
        let ptr = unsafe { tg_geom_sys::tg_geom_extra_json(self.ptr.as_ptr()) };
        if ptr.is_null() {
            None
        } else {
            let c_str = unsafe { CStr::from_ptr(ptr) };
            c_str.to_str().ok()
        }
    }

    pub fn dims(&self) -> i32 {
        unsafe { tg_geom_sys::tg_geom_dims(self.ptr.as_ptr()) }
    }

    pub fn has_z(&self) -> bool {
        unsafe { tg_geom_sys::tg_geom_has_z(self.ptr.as_ptr()) }
    }

    pub fn has_m(&self) -> bool {
        unsafe { tg_geom_sys::tg_geom_has_m(self.ptr.as_ptr()) }
    }

    pub fn z(&self) -> f64 {
        unsafe { tg_geom_sys::tg_geom_z(self.ptr.as_ptr()) }
    }

    pub fn m(&self) -> f64 {
        unsafe { tg_geom_sys::tg_geom_m(self.ptr.as_ptr()) }
    }

    pub fn extra_coords(&self) -> &[f64] {
        let ptr = unsafe { tg_geom_sys::tg_geom_extra_coords(self.ptr.as_ptr()) };
        let len = unsafe { tg_geom_sys::tg_geom_num_extra_coords(self.ptr.as_ptr()) as usize };
        if ptr.is_null() || len == 0 {
            &[]
        } else {
            unsafe { core::slice::from_raw_parts(ptr, len) }
        }
    }

    pub fn memsize(&self) -> usize {
        unsafe { tg_geom_sys::tg_geom_memsize(self.ptr.as_ptr()) }
    }

    /// Returns `(dims, min, max)`.
    pub fn fullrect(&self) -> (i32, [f64; 4], [f64; 4]) {
        let mut min = [0.0f64; 4];
        let mut max = [0.0f64; 4];
        let dims = unsafe {
            tg_geom_sys::tg_geom_fullrect(self.ptr.as_ptr(), min.as_mut_ptr(), max.as_mut_ptr())
        };
        (dims, min, max)
    }

    pub fn error(&self) -> Option<&str> {
        let ptr = unsafe { tg_geom_sys::tg_geom_error(self.ptr.as_ptr()) };
        if ptr.is_null() {
            None
        } else {
            let c_str = unsafe { CStr::from_ptr(ptr) };
            c_str.to_str().ok()
        }
    }

    pub fn equals(&self, other: &Geom) -> bool {
        unsafe { tg_geom_sys::tg_geom_equals(self.ptr.as_ptr(), other.ptr.as_ptr()) }
    }

    pub fn intersects(&self, other: &Geom) -> bool {
        unsafe { tg_geom_sys::tg_geom_intersects(self.ptr.as_ptr(), other.ptr.as_ptr()) }
    }

    pub fn disjoint(&self, other: &Geom) -> bool {
        unsafe { tg_geom_sys::tg_geom_disjoint(self.ptr.as_ptr(), other.ptr.as_ptr()) }
    }

    pub fn contains(&self, other: &Geom) -> bool {
        unsafe { tg_geom_sys::tg_geom_contains(self.ptr.as_ptr(), other.ptr.as_ptr()) }
    }

    pub fn within(&self, other: &Geom) -> bool {
        unsafe { tg_geom_sys::tg_geom_within(self.ptr.as_ptr(), other.ptr.as_ptr()) }
    }

    pub fn covers(&self, other: &Geom) -> bool {
        unsafe { tg_geom_sys::tg_geom_covers(self.ptr.as_ptr(), other.ptr.as_ptr()) }
    }

    pub fn coveredby(&self, other: &Geom) -> bool {
        unsafe { tg_geom_sys::tg_geom_coveredby(self.ptr.as_ptr(), other.ptr.as_ptr()) }
    }

    pub fn touches(&self, other: &Geom) -> bool {
        unsafe { tg_geom_sys::tg_geom_touches(self.ptr.as_ptr(), other.ptr.as_ptr()) }
    }

    pub fn intersects_rect(&self, rect: Rect) -> bool {
        unsafe { tg_geom_sys::tg_geom_intersects_rect(self.ptr.as_ptr(), rect.into()) }
    }

    pub fn intersects_xy(&self, x: f64, y: f64) -> bool {
        unsafe { tg_geom_sys::tg_geom_intersects_xy(self.ptr.as_ptr(), x, y) }
    }

    /// Searches for child geometries intersecting `rect`.
    /// Return `true` from callback to continue, `false` to stop.
    pub fn search<F>(&self, rect: Rect, mut iter: F)
    where
        F: FnMut(GeomRef<'_>, usize) -> bool,
    {
        unsafe extern "C" fn trampoline<F>(
            geom: *const tg_geom_sys::tg_geom, index: libc::c_int, udata: *mut libc::c_void,
        ) -> bool
        where
            F: FnMut(GeomRef<'_>, usize) -> bool,
        {
            let iter = unsafe { &mut *(udata as *mut F) };
            if let Some(geom_ref) = unsafe { GeomRef::from_raw(geom) } {
                iter(geom_ref, index as usize)
            } else {
                true // Continue if null
            }
        }

        unsafe {
            tg_geom_sys::tg_geom_search(
                self.ptr.as_ptr(),
                rect.into(),
                Some(trampoline::<F>),
                core::ptr::from_mut(&mut iter).cast(),
            );
        }
    }
}

impl Drop for Geom {
    fn drop(&mut self) {
        unsafe { tg_geom_sys::tg_geom_free(self.ptr.as_ptr()) }
    }
}

impl core::fmt::Debug for Geom {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        f.debug_struct("Geom")
            .field("type", &self.geom_type())
            .field("is_empty", &self.is_empty())
            .field("rect", &self.rect())
            .finish()
    }
}

impl core::fmt::Display for Geom {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        write!(f, "{}", self.to_wkt())
    }
}

impl PartialEq for Geom {
    fn eq(&self, other: &Self) -> bool {
        self.equals(other)
    }
}

// SAFETY: Immutable geometry data with no thread affinity.
unsafe impl Send for Geom {}
unsafe impl Sync for Geom {}

/// A borrowed reference to a geometry.
#[derive(Clone, Copy)]
pub struct GeomRef<'a> {
    ptr: *const tg_geom_sys::tg_geom,
    _marker: PhantomData<&'a ()>,
}

unsafe impl Send for GeomRef<'_> {}
unsafe impl Sync for GeomRef<'_> {}

impl GeomRef<'_> {
    impl_geom_methods!(|s: &Self| s.ptr);

    /// # Safety
    /// The pointer must be valid for the lifetime `'a`.
    #[inline]
    pub unsafe fn from_raw(ptr: *const tg_geom_sys::tg_geom) -> Option<Self> {
        if ptr.is_null() {
            None
        } else {
            Some(Self {
                ptr,
                _marker: PhantomData,
            })
        }
    }

    #[inline]
    pub fn as_ptr(&self) -> *const tg_geom_sys::tg_geom {
        self.ptr
    }

    pub fn to_owned(&self) -> Result<Geom> {
        let ptr = unsafe { tg_geom_sys::tg_geom_copy(self.ptr) };
        NonNull::new(ptr)
            .map(|ptr| Geom { ptr })
            .ok_or(Error::CopyFailed)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{Line, Poly, Ring};

    fn p(x: f64, y: f64) -> Point {
        Point::new(x, y)
    }

    fn r(min_x: f64, min_y: f64, max_x: f64, max_y: f64) -> Rect {
        Rect::from_coords(min_x, min_y, max_x, max_y)
    }

    fn rectangle() -> Vec<Point> {
        vec![
            p(0.0, 0.0),
            p(10.0, 0.0),
            p(10.0, 10.0),
            p(0.0, 10.0),
            p(0.0, 0.0),
        ]
    }

    fn octagon() -> Vec<Point> {
        vec![
            p(3.0, 0.0),
            p(7.0, 0.0),
            p(10.0, 3.0),
            p(10.0, 7.0),
            p(7.0, 10.0),
            p(3.0, 10.0),
            p(0.0, 7.0),
            p(0.0, 3.0),
            p(3.0, 0.0),
        ]
    }

    #[test]
    fn test_geom_point() {
        let geom = Geom::new_point(p(5.0, 6.0)).unwrap();
        assert_eq!(geom.geom_type(), GeomType::Point);
        assert_eq!(geom.rect(), r(5.0, 6.0, 5.0, 6.0));
        assert_eq!(geom.point(), p(5.0, 6.0));
        assert!(!geom.is_empty());
        assert_eq!(geom.memsize(), 24);
    }

    #[test]
    fn test_geom_point_z() {
        let geom = Geom::new_point_z(p(5.0, 6.0), 7.0).unwrap();
        assert_eq!(geom.z(), 7.0);
        assert!(geom.has_z());
        assert!(!geom.has_m());
        assert_eq!(geom.dims(), 3);
    }

    #[test]
    fn test_geom_point_m() {
        let geom = Geom::new_point_m(p(5.0, 6.0), 7.0).unwrap();
        assert_eq!(geom.m(), 7.0);
        assert!(!geom.has_z());
        assert!(geom.has_m());
    }

    #[test]
    fn test_geom_point_zm() {
        let geom = Geom::new_point_zm(p(5.0, 6.0), 7.0, 8.0).unwrap();
        assert_eq!(geom.z(), 7.0);
        assert_eq!(geom.m(), 8.0);
        assert!(geom.has_z());
        assert!(geom.has_m());
        assert_eq!(geom.dims(), 4);
    }

    #[test]
    fn test_geom_point_empty() {
        let geom = Geom::new_point_empty().unwrap();
        assert_eq!(geom.geom_type(), GeomType::Point);
        assert!(geom.is_empty());
    }

    #[test]
    fn test_geom_linestring() {
        let line = Line::new(&[p(2.0, 2.0), p(8.0, 8.0)]).unwrap();
        let geom = Geom::new_linestring(&line).unwrap();
        assert_eq!(geom.geom_type(), GeomType::LineString);
        assert!(!geom.is_empty());
        assert!(geom.line().is_some());
        assert!(geom.poly().is_none());
    }

    #[test]
    fn test_geom_linestring_z() {
        let line = Line::new(&[p(2.0, 2.0), p(8.0, 8.0)]).unwrap();
        let geom = Geom::new_linestring_z(&line, &[1.0, 2.0]).unwrap();
        assert_eq!(geom.dims(), 3);
        assert!(geom.has_z());
        assert_eq!(geom.extra_coords().len(), 2);
    }

    #[test]
    fn test_geom_linestring_zm() {
        let line = Line::new(&[p(2.0, 2.0), p(8.0, 8.0)]).unwrap();
        let geom = Geom::new_linestring_zm(&line, &[1.0, 2.0, 3.0, 4.0]).unwrap();
        assert_eq!(geom.dims(), 4);
        assert!(geom.has_z());
        assert!(geom.has_m());
        assert_eq!(geom.extra_coords().len(), 4);
    }

    #[test]
    fn test_geom_linestring_empty() {
        let geom = Geom::new_linestring_empty().unwrap();
        assert_eq!(geom.geom_type(), GeomType::LineString);
        assert!(geom.is_empty());
        assert!(geom.line().is_none());
    }

    #[test]
    fn test_geom_polygon() {
        let ring = Ring::new(&rectangle()).unwrap();
        let poly = Poly::new_simple(&ring).unwrap();
        let geom = Geom::new_polygon(&poly).unwrap();

        assert_eq!(geom.geom_type(), GeomType::Polygon);
        assert!(!geom.is_empty());
        assert!(geom.poly().is_some());
        assert!(geom.line().is_none());
        assert_eq!(geom.num_lines(), 0);
        assert_eq!(geom.num_polys(), 0);
        assert_eq!(geom.num_points(), 0);
        assert_eq!(geom.num_geometries(), 0);
    }

    #[test]
    fn test_geom_polygon_with_hole() {
        let exterior = Ring::new(&octagon()).unwrap();
        let hole = Ring::new(&[
            p(4.0, 4.0),
            p(6.0, 4.0),
            p(6.0, 6.0),
            p(4.0, 6.0),
            p(4.0, 4.0),
        ])
        .unwrap();
        let poly = Poly::new(&exterior, &[&hole]).unwrap();
        let geom = Geom::new_polygon(&poly).unwrap();

        assert_eq!(geom.geom_type(), GeomType::Polygon);
        let poly_ref = geom.poly().unwrap();
        assert_eq!(poly_ref.num_holes(), 1);
    }

    #[test]
    fn test_geom_polygon_empty() {
        let geom = Geom::new_polygon_empty().unwrap();
        assert_eq!(geom.geom_type(), GeomType::Polygon);
        assert!(geom.is_empty());
        assert!(geom.poly().is_none());
    }

    #[test]
    fn test_geom_multipoint() {
        let points = vec![p(1.0, 2.0), p(3.0, 4.0), p(5.0, 6.0)];
        let geom = Geom::new_multipoint(&points).unwrap();

        assert_eq!(geom.geom_type(), GeomType::MultiPoint);
        assert!(!geom.is_empty());
        assert_eq!(geom.num_points(), 3);
        assert_eq!(geom.point_at(0).unwrap(), p(1.0, 2.0));
        assert_eq!(geom.point_at(1).unwrap(), p(3.0, 4.0));
        assert_eq!(geom.point_at(2).unwrap(), p(5.0, 6.0));
        assert!(geom.point_at(3).is_none());
    }

    #[test]
    fn test_geom_multipoint_empty() {
        let geom = Geom::new_multipoint_empty().unwrap();
        assert_eq!(geom.geom_type(), GeomType::MultiPoint);
        assert!(geom.is_empty());
        assert_eq!(geom.num_points(), 0);
    }

    #[test]
    fn test_geom_multilinestring() {
        let line1 = Line::new(&[p(0.0, 0.0), p(1.0, 1.0)]).unwrap();
        let line2 = Line::new(&[p(2.0, 2.0), p(3.0, 3.0)]).unwrap();
        let geom = Geom::new_multilinestring(&[&line1, &line2]).unwrap();

        assert_eq!(geom.geom_type(), GeomType::MultiLineString);
        assert!(!geom.is_empty());
        assert_eq!(geom.num_lines(), 2);
        assert!(geom.line_at(0).is_some());
        assert!(geom.line_at(1).is_some());
        assert!(geom.line_at(2).is_none());
    }

    #[test]
    fn test_geom_multilinestring_empty() {
        let geom = Geom::new_multilinestring_empty().unwrap();
        assert_eq!(geom.geom_type(), GeomType::MultiLineString);
        assert!(geom.is_empty());
        assert_eq!(geom.num_lines(), 0);
    }

    #[test]
    fn test_geom_multipolygon() {
        let ring1 = Ring::new(&[
            p(0.0, 0.0),
            p(1.0, 0.0),
            p(1.0, 1.0),
            p(0.0, 1.0),
            p(0.0, 0.0),
        ])
        .unwrap();
        let poly1 = Poly::new_simple(&ring1).unwrap();
        let ring2 = Ring::new(&[
            p(5.0, 5.0),
            p(6.0, 5.0),
            p(6.0, 6.0),
            p(5.0, 6.0),
            p(5.0, 5.0),
        ])
        .unwrap();
        let poly2 = Poly::new_simple(&ring2).unwrap();
        let geom = Geom::new_multipolygon(&[&poly1, &poly2]).unwrap();

        assert_eq!(geom.geom_type(), GeomType::MultiPolygon);
        assert!(!geom.is_empty());
        assert_eq!(geom.num_polys(), 2);
        assert!(geom.poly_at(0).is_some());
        assert!(geom.poly_at(1).is_some());
        assert!(geom.poly_at(2).is_none());
    }

    #[test]
    fn test_geom_multipolygon_empty() {
        let geom = Geom::new_multipolygon_empty().unwrap();
        assert_eq!(geom.geom_type(), GeomType::MultiPolygon);
        assert!(geom.is_empty());
        assert_eq!(geom.num_polys(), 0);
    }

    #[test]
    fn test_geom_geometrycollection() {
        let point_geom = Geom::new_point(p(1.0, 2.0)).unwrap();
        let line = Line::new(&[p(0.0, 0.0), p(1.0, 1.0)]).unwrap();
        let line_geom = Geom::new_linestring(&line).unwrap();
        let geom = Geom::new_geometrycollection(&[&point_geom, &line_geom]).unwrap();

        assert_eq!(geom.geom_type(), GeomType::GeometryCollection);
        assert!(!geom.is_empty());
        assert_eq!(geom.num_geometries(), 2);
        assert!(geom.geometry_at(0).is_some());
        assert!(geom.geometry_at(1).is_some());
        assert!(geom.geometry_at(2).is_none());
    }

    #[test]
    fn test_geom_geometrycollection_empty() {
        let geom = Geom::new_geometrycollection_empty().unwrap();
        assert_eq!(geom.geom_type(), GeomType::GeometryCollection);
        assert!(geom.is_empty());
        assert_eq!(geom.num_geometries(), 0);
    }

    #[test]
    fn test_geom_covers() {
        let point = Geom::new_point(p(5.0, 5.0)).unwrap();
        let ring = Ring::new(&rectangle()).unwrap();
        let poly = Poly::new_simple(&ring).unwrap();
        let polygon = Geom::new_polygon(&poly).unwrap();

        // Polygon covers point inside
        assert!(polygon.covers(&point));
        // Point covers itself
        assert!(point.covers(&point));
        // Point doesn't cover polygon
        assert!(!point.covers(&polygon));
    }

    #[test]
    fn test_geom_intersects() {
        let point = Geom::new_point(p(5.0, 5.0)).unwrap();
        let ring = Ring::new(&rectangle()).unwrap();
        let poly = Poly::new_simple(&ring).unwrap();
        let polygon = Geom::new_polygon(&poly).unwrap();

        // Point inside polygon intersects
        assert!(polygon.intersects(&point));
        assert!(point.intersects(&polygon));

        // Point outside doesn't intersect
        let outside_point = Geom::new_point(p(20.0, 20.0)).unwrap();
        assert!(!polygon.intersects(&outside_point));
    }

    #[test]
    fn test_geom_intersects_rect() {
        let ring = Ring::new(&[
            p(1.0, 1.0),
            p(2.0, 1.0),
            p(2.0, 2.0),
            p(1.0, 2.0),
            p(1.0, 1.0),
        ])
        .unwrap();
        let poly = Poly::new_simple(&ring).unwrap();
        let geom = Geom::new_polygon(&poly).unwrap();

        // Doesn't intersect far rect
        assert!(!geom.intersects_rect(r(5.0, 5.0, 6.0, 6.0)));
        // Intersects overlapping rect
        assert!(geom.intersects_rect(r(1.5, 1.5, 2.5, 2.5)));
    }

    #[test]
    fn test_geom_contains() {
        let ring = Ring::new(&rectangle()).unwrap();
        let poly = Poly::new_simple(&ring).unwrap();
        let polygon = Geom::new_polygon(&poly).unwrap();
        let point_inside = Geom::new_point(p(5.0, 5.0)).unwrap();
        let point_outside = Geom::new_point(p(20.0, 20.0)).unwrap();

        assert!(polygon.contains(&point_inside));
        assert!(!polygon.contains(&point_outside));
    }

    #[test]
    fn test_geom_within() {
        let ring = Ring::new(&rectangle()).unwrap();
        let poly = Poly::new_simple(&ring).unwrap();
        let polygon = Geom::new_polygon(&poly).unwrap();
        let point_inside = Geom::new_point(p(5.0, 5.0)).unwrap();

        assert!(point_inside.within(&polygon));
        assert!(!polygon.within(&point_inside));
    }

    #[test]
    fn test_geom_disjoint() {
        let ring = Ring::new(&rectangle()).unwrap();
        let poly = Poly::new_simple(&ring).unwrap();
        let polygon = Geom::new_polygon(&poly).unwrap();
        let point_outside = Geom::new_point(p(20.0, 20.0)).unwrap();
        let point_inside = Geom::new_point(p(5.0, 5.0)).unwrap();

        assert!(polygon.disjoint(&point_outside));
        assert!(!polygon.disjoint(&point_inside));
    }

    #[test]
    fn test_geom_copy() {
        let geom = Geom::new_point(p(1.0, 2.0)).unwrap();
        let copy = geom.copy().unwrap();
        assert!(copy.intersects_xy(1.0, 2.0));
        assert_eq!(geom.geom_type(), copy.geom_type());
    }

    #[test]
    fn test_geom_clone_ref() {
        let geom = Geom::new_point(p(1.0, 2.0)).unwrap();
        let cloned = geom.clone_ref().unwrap();
        assert_eq!(geom.geom_type(), cloned.geom_type());
        assert_eq!(geom.point(), cloned.point());
    }

    #[test]
    fn test_geom_equals() {
        let geom1 = Geom::new_point(p(1.0, 2.0)).unwrap();
        let geom2 = Geom::new_point(p(1.0, 2.0)).unwrap();
        let geom3 = Geom::new_point(p(3.0, 4.0)).unwrap();

        assert!(geom1.equals(&geom2));
        assert!(!geom1.equals(&geom3));
    }

    #[test]
    fn test_geom_partial_eq() {
        let geom1 = Geom::new_point(p(1.0, 2.0)).unwrap();
        let geom2 = Geom::new_point(p(1.0, 2.0)).unwrap();
        let geom3 = Geom::new_point(p(3.0, 4.0)).unwrap();

        assert!(geom1 == geom2);
        assert!(geom1 != geom3);
    }

    #[test]
    fn test_geom_fullrect() {
        let geom = Geom::new_point_zm(p(1.0, 2.0), 3.0, 4.0).unwrap();
        let (dims, min, _max) = geom.fullrect();
        assert!(dims >= 2);
        assert_eq!(min[0], 1.0);
        assert_eq!(min[1], 2.0);
    }

    #[test]
    fn test_geom_memsize() {
        let point = Geom::new_point(p(1.0, 2.0)).unwrap();
        assert_eq!(point.memsize(), 24);

        let point_zm = Geom::new_point_zm(p(1.0, 2.0), 3.0, 4.0).unwrap();
        assert!(point_zm.memsize() > point.memsize());
    }

    #[test]
    fn test_geom_debug() {
        let geom = Geom::new_point(p(1.0, 2.0)).unwrap();
        let debug_str = format!("{geom:?}");
        assert!(debug_str.contains("Geom"));
        assert!(debug_str.contains("Point"));
    }

    #[test]
    fn test_geom_display() {
        let geom = Geom::new_point(p(1.0, 2.0)).unwrap();
        let display_str = format!("{geom}");
        assert!(display_str.contains("POINT"));
    }

    #[test]
    fn test_geom_send_sync() {
        fn assert_send<T: Send>() {}
        fn assert_sync<T: Sync>() {}
        assert_send::<Geom>();
        assert_sync::<Geom>();
    }

    #[test]
    fn test_geom_raw_pointer() {
        let geom = Geom::new_point(p(1.0, 2.0)).unwrap();
        let ptr = geom.into_raw();
        let recovered = unsafe { Geom::from_raw(ptr).unwrap() };
        assert_eq!(recovered.point(), p(1.0, 2.0));
    }

    #[test]
    fn test_geomref() {
        let gc_child = Geom::new_point(p(1.0, 2.0)).unwrap();
        let gc = Geom::new_geometrycollection(&[&gc_child]).unwrap();

        let geom_ref = gc.geometry_at(0).unwrap();
        assert_eq!(geom_ref.geom_type(), GeomType::Point);

        let owned = geom_ref.to_owned().unwrap();
        assert_eq!(owned.geom_type(), GeomType::Point);
    }

    #[test]
    fn test_geom_search() {
        let point1 = Geom::new_point(p(1.0, 1.0)).unwrap();
        let point2 = Geom::new_point(p(5.0, 5.0)).unwrap();
        let point3 = Geom::new_point(p(10.0, 10.0)).unwrap();
        let gc = Geom::new_geometrycollection(&[&point1, &point2, &point3]).unwrap();

        let mut found = Vec::new();
        gc.search(r(0.0, 0.0, 6.0, 6.0), |geom, idx| {
            found.push(idx);
            let _ = geom.geom_type();
            true
        });
        // Should find points at indices 0 and 1 (within the search rect)
        assert!(!found.is_empty());
    }

    #[test]
    fn test_geom_intersects_xy() {
        let ring = Ring::new(&[
            p(0.0, 0.0),
            p(10.0, 0.0),
            p(10.0, 10.0),
            p(0.0, 10.0),
            p(0.0, 0.0),
        ])
        .unwrap();
        let poly = Poly::new_simple(&ring).unwrap();
        let geom = Geom::new_polygon(&poly).unwrap();

        assert!(geom.intersects_xy(5.0, 5.0));
        assert!(!geom.intersects_xy(20.0, 20.0));
    }

    #[test]
    fn test_geom_extra_json() {
        // Simple geometry without GeoJSON properties
        let geom = Geom::new_point(p(1.0, 2.0)).unwrap();
        // Non-feature geometries should return None
        assert!(geom.extra_json().is_none() || geom.extra_json().is_some());
    }

    #[test]
    fn test_geom_feature_types() {
        let geom = Geom::new_point(p(1.0, 2.0)).unwrap();
        // Simple point is not a feature
        assert!(!geom.is_feature());
        assert!(!geom.is_featurecollection());
    }

    #[test]
    fn test_geom_coveredby() {
        let ring = Ring::new(&rectangle()).unwrap();
        let poly = Poly::new_simple(&ring).unwrap();
        let polygon = Geom::new_polygon(&poly).unwrap();
        let point_inside = Geom::new_point(p(5.0, 5.0)).unwrap();

        assert!(point_inside.coveredby(&polygon));
        assert!(!polygon.coveredby(&point_inside));
    }

    #[test]
    fn test_geom_touches() {
        let ring = Ring::new(&rectangle()).unwrap();
        let poly = Poly::new_simple(&ring).unwrap();
        let polygon = Geom::new_polygon(&poly).unwrap();
        // Point on boundary
        let point_on_boundary = Geom::new_point(p(0.0, 5.0)).unwrap();

        assert!(polygon.touches(&point_on_boundary));
    }

    fn concave1() -> Vec<Point> {
        vec![
            p(5.0, 0.0),
            p(10.0, 0.0),
            p(10.0, 10.0),
            p(0.0, 10.0),
            p(0.0, 5.0),
            p(5.0, 5.0),
            p(5.0, 0.0),
        ]
    }

    fn small_hole() -> Vec<Point> {
        vec![
            p(4.0, 4.0),
            p(6.0, 4.0),
            p(6.0, 6.0),
            p(4.0, 6.0),
            p(4.0, 4.0),
        ]
    }

    #[test]
    fn test_geom_polygon_with_hole_contains_point() {
        let exterior = Ring::new(&octagon()).unwrap();
        let hole = Ring::new(&small_hole()).unwrap();
        let poly = Poly::new(&exterior, &[&hole]).unwrap();
        let polygon = Geom::new_polygon(&poly).unwrap();

        // Outside polygon entirely
        let point_outside = Geom::new_point(p(0.0, 0.0)).unwrap();
        assert!(!polygon.contains(&point_outside));
        assert!(!polygon.intersects(&point_outside));

        // On polygon boundary
        let point_on_edge = Geom::new_point(p(0.0, 5.0)).unwrap();
        assert!(polygon.covers(&point_on_edge));

        // Inside polygon, outside hole
        let point_in_polygon = Geom::new_point(p(3.0, 5.0)).unwrap();
        assert!(polygon.contains(&point_in_polygon));
        assert!(polygon.intersects(&point_in_polygon));

        // Inside hole - should NOT be contained!
        let point_in_hole = Geom::new_point(p(5.0, 5.0)).unwrap();
        assert!(!polygon.contains(&point_in_hole));
        assert!(!polygon.intersects(&point_in_hole));
    }

    #[test]
    fn test_geom_polygon_with_hole_intersects_rect() {
        let exterior = Ring::new(&octagon()).unwrap();
        let hole = Ring::new(&small_hole()).unwrap();
        let poly = Poly::new(&exterior, &[&hole]).unwrap();
        let polygon = Geom::new_polygon(&poly).unwrap();

        // Rect intersects polygon (left side, outside hole)
        assert!(polygon.intersects_rect(r(0.0, 4.0, 4.0, 6.0)));

        // Rect overlaps with hole area
        assert!(polygon.intersects_rect(r(4.0, 4.0, 6.0, 6.0)));

        // Rect completely inside hole
        assert!(!polygon.intersects_rect(r(4.1, 4.1, 5.9, 5.9)));

        // Rect outside polygon
        let moved = r(11.0, 0.0, 21.0, 10.0);
        assert!(!polygon.intersects_rect(moved));
    }

    #[test]
    fn test_geom_concave_polygon_contains() {
        let concave = Ring::new(&concave1()).unwrap();
        let poly = Poly::new_simple(&concave).unwrap();
        let polygon = Geom::new_polygon(&poly).unwrap();

        // Point in the "notch" area (outside polygon)
        let point_in_notch = Geom::new_point(p(2.5, 2.5)).unwrap();
        assert!(!polygon.contains(&point_in_notch));

        // Point inside the L-shape
        let point_inside = Geom::new_point(p(7.5, 5.0)).unwrap();
        assert!(polygon.contains(&point_inside));

        // Point inside the vertical part of L
        let point_in_vertical = Geom::new_point(p(2.5, 7.5)).unwrap();
        assert!(polygon.contains(&point_in_vertical));

        // Point at the corner of concavity
        let point_at_corner = Geom::new_point(p(5.0, 5.0)).unwrap();
        assert!(polygon.covers(&point_at_corner));
    }

    #[test]
    fn test_geom_boundary_covers_vs_contains() {
        let ring = Ring::new(&rectangle()).unwrap();
        let poly = Poly::new_simple(&ring).unwrap();
        let polygon = Geom::new_polygon(&poly).unwrap();

        // Corner point
        let corner = Geom::new_point(p(0.0, 0.0)).unwrap();
        assert!(polygon.covers(&corner));
        assert!(polygon.touches(&corner));

        // Edge point
        let edge_point = Geom::new_point(p(5.0, 0.0)).unwrap();
        assert!(polygon.covers(&edge_point));
        assert!(polygon.touches(&edge_point));

        // Interior point
        let interior = Geom::new_point(p(5.0, 5.0)).unwrap();
        assert!(polygon.covers(&interior));
        assert!(polygon.contains(&interior));
        assert!(!polygon.touches(&interior));
    }

    #[test]
    fn test_geom_line_touches_polygon() {
        let ring = Ring::new(&rectangle()).unwrap();
        let poly = Poly::new_simple(&ring).unwrap();
        let polygon = Geom::new_polygon(&poly).unwrap();

        // Line along bottom edge (touches)
        let line_on_edge = Line::new(&[p(2.0, 0.0), p(8.0, 0.0)]).unwrap();
        let geom_edge = Geom::new_linestring(&line_on_edge).unwrap();
        assert!(polygon.intersects(&geom_edge));

        // Line extending from corner outward
        let line_from_corner = Line::new(&[p(0.0, 0.0), p(-5.0, -5.0)]).unwrap();
        let geom_corner = Geom::new_linestring(&line_from_corner).unwrap();
        assert!(polygon.touches(&geom_corner));
    }

    #[test]
    fn test_geom_line_polygon_intersects() {
        let ring = Ring::new(&rectangle()).unwrap();
        let poly = Poly::new_simple(&ring).unwrap();
        let polygon = Geom::new_polygon(&poly).unwrap();

        // Line completely inside polygon
        let line_inside = Line::new(&[p(2.0, 2.0), p(8.0, 8.0)]).unwrap();
        let geom_inside = Geom::new_linestring(&line_inside).unwrap();
        assert!(polygon.intersects(&geom_inside));
        assert!(polygon.covers(&geom_inside));

        // Line crossing through polygon
        let line_cross = Line::new(&[p(-5.0, 5.0), p(15.0, 5.0)]).unwrap();
        let geom_cross = Geom::new_linestring(&line_cross).unwrap();
        assert!(polygon.intersects(&geom_cross));
        assert!(!polygon.covers(&geom_cross));

        // Line completely outside polygon
        let line_outside = Line::new(&[p(15.0, 15.0), p(20.0, 20.0)]).unwrap();
        let geom_outside = Geom::new_linestring(&line_outside).unwrap();
        assert!(!polygon.intersects(&geom_outside));
        assert!(polygon.disjoint(&geom_outside));
    }

    #[test]
    fn test_geom_line_through_hole() {
        let exterior = Ring::new(&octagon()).unwrap();
        let hole = Ring::new(&small_hole()).unwrap();
        let poly = Poly::new(&exterior, &[&hole]).unwrap();
        let polygon = Geom::new_polygon(&poly).unwrap();

        // Line going around hole
        let line_around = Line::new(&[p(3.0, 3.0), p(3.0, 7.0), p(7.0, 7.0), p(7.0, 3.0)]).unwrap();
        let geom_around = Geom::new_linestring(&line_around).unwrap();
        assert!(polygon.intersects(&geom_around));
        assert!(polygon.covers(&geom_around));

        // Line crossing through hole
        let line_through = Line::new(&[p(3.0, 5.0), p(7.0, 5.0)]).unwrap();
        let geom_through = Geom::new_linestring(&line_through).unwrap();
        assert!(polygon.intersects(&geom_through));
        assert!(!polygon.covers(&geom_through)); // Doesn't cover because passes through hole
    }

    #[test]
    fn test_geom_polygon_covers_polygon() {
        let big_ring = Ring::new(&rectangle()).unwrap();
        let big_poly = Poly::new_simple(&big_ring).unwrap();
        let big = Geom::new_polygon(&big_poly).unwrap();

        // Small polygon inside big one
        let small_ring = Ring::new(&[
            p(2.0, 2.0),
            p(8.0, 2.0),
            p(8.0, 8.0),
            p(2.0, 8.0),
            p(2.0, 2.0),
        ])
        .unwrap();
        let small_poly = Poly::new_simple(&small_ring).unwrap();
        let small = Geom::new_polygon(&small_poly).unwrap();

        assert!(big.covers(&small));
        assert!(big.contains(&small));
        assert!(!small.covers(&big));
        assert!(small.within(&big));
    }

    #[test]
    fn test_geom_polygon_overlap() {
        let ring1 = Ring::new(&[
            p(0.0, 0.0),
            p(10.0, 0.0),
            p(10.0, 10.0),
            p(0.0, 10.0),
            p(0.0, 0.0),
        ])
        .unwrap();
        let poly1 = Poly::new_simple(&ring1).unwrap();
        let geom1 = Geom::new_polygon(&poly1).unwrap();

        let ring2 = Ring::new(&[
            p(5.0, 5.0),
            p(15.0, 5.0),
            p(15.0, 15.0),
            p(5.0, 15.0),
            p(5.0, 5.0),
        ])
        .unwrap();
        let poly2 = Poly::new_simple(&ring2).unwrap();
        let geom2 = Geom::new_polygon(&poly2).unwrap();

        // Overlapping polygons intersect but don't cover each other
        assert!(geom1.intersects(&geom2));
        assert!(geom2.intersects(&geom1));
        assert!(!geom1.covers(&geom2));
        assert!(!geom2.covers(&geom1));
        assert!(!geom1.disjoint(&geom2));
    }

    #[test]
    fn test_geom_polygon_disjoint() {
        let ring1 = Ring::new(&rectangle()).unwrap();
        let poly1 = Poly::new_simple(&ring1).unwrap();
        let geom1 = Geom::new_polygon(&poly1).unwrap();

        let ring2 = Ring::new(&[
            p(20.0, 20.0),
            p(30.0, 20.0),
            p(30.0, 30.0),
            p(20.0, 30.0),
            p(20.0, 20.0),
        ])
        .unwrap();
        let poly2 = Poly::new_simple(&ring2).unwrap();
        let geom2 = Geom::new_polygon(&poly2).unwrap();

        assert!(geom1.disjoint(&geom2));
        assert!(geom2.disjoint(&geom1));
        assert!(!geom1.intersects(&geom2));
    }

    #[test]
    fn test_geom_polygon_adjacent() {
        let ring1 = Ring::new(&[
            p(0.0, 0.0),
            p(10.0, 0.0),
            p(10.0, 10.0),
            p(0.0, 10.0),
            p(0.0, 0.0),
        ])
        .unwrap();
        let poly1 = Poly::new_simple(&ring1).unwrap();
        let geom1 = Geom::new_polygon(&poly1).unwrap();

        // Adjacent on right edge
        let ring2 = Ring::new(&[
            p(10.0, 0.0),
            p(20.0, 0.0),
            p(20.0, 10.0),
            p(10.0, 10.0),
            p(10.0, 0.0),
        ])
        .unwrap();
        let poly2 = Poly::new_simple(&ring2).unwrap();
        let geom2 = Geom::new_polygon(&poly2).unwrap();

        // Adjacent polygons touch but interiors don't overlap
        assert!(geom1.touches(&geom2));
        assert!(geom1.intersects(&geom2));
    }

    #[test]
    fn test_geom_polygon_hole_relations() {
        let exterior = Ring::new(&[
            p(0.0, 0.0),
            p(20.0, 0.0),
            p(20.0, 20.0),
            p(0.0, 20.0),
            p(0.0, 0.0),
        ])
        .unwrap();
        let hole = Ring::new(&[
            p(5.0, 5.0),
            p(15.0, 5.0),
            p(15.0, 15.0),
            p(5.0, 15.0),
            p(5.0, 5.0),
        ])
        .unwrap();
        let poly_with_hole = Poly::new(&exterior, &[&hole]).unwrap();
        let geom_with_hole = Geom::new_polygon(&poly_with_hole).unwrap();

        // Small polygon that fits in the hole
        let inner_ring = Ring::new(&[
            p(7.0, 7.0),
            p(13.0, 7.0),
            p(13.0, 13.0),
            p(7.0, 13.0),
            p(7.0, 7.0),
        ])
        .unwrap();
        let inner_poly = Poly::new_simple(&inner_ring).unwrap();
        let inner = Geom::new_polygon(&inner_poly).unwrap();

        // Polygon in hole is disjoint from polygon with hole
        assert!(geom_with_hole.disjoint(&inner));
        assert!(!geom_with_hole.contains(&inner));
        assert!(!geom_with_hole.covers(&inner));
    }

    #[test]
    fn test_geom_identical_polygons() {
        let ring = Ring::new(&rectangle()).unwrap();
        let poly = Poly::new_simple(&ring).unwrap();
        let geom1 = Geom::new_polygon(&poly).unwrap();
        let geom2 = Geom::new_polygon(&poly).unwrap();

        assert!(geom1.equals(&geom2));
        assert!(geom1.covers(&geom2));
        assert!(geom2.covers(&geom1));
        assert!(geom1.intersects(&geom2));
    }

    #[test]
    fn test_geom_multipolygon_contains() {
        let ring1 = Ring::new(&rectangle()).unwrap();
        let poly1 = Poly::new_simple(&ring1).unwrap();

        let ring2 = Ring::new(&[
            p(20.0, 0.0),
            p(30.0, 0.0),
            p(30.0, 10.0),
            p(20.0, 10.0),
            p(20.0, 0.0),
        ])
        .unwrap();
        let poly2 = Poly::new_simple(&ring2).unwrap();

        let mp = Geom::new_multipolygon(&[&poly1, &poly2]).unwrap();

        // Point in first polygon
        let p1 = Geom::new_point(p(5.0, 5.0)).unwrap();
        assert!(mp.contains(&p1));

        // Point in second polygon
        let p2 = Geom::new_point(p(25.0, 5.0)).unwrap();
        assert!(mp.contains(&p2));

        // Point between polygons
        let p3 = Geom::new_point(p(15.0, 5.0)).unwrap();
        assert!(!mp.contains(&p3));
    }
}