geo-index 0.3.4

use std::cmp::Reverse;
use std::collections::{BinaryHeap, VecDeque};
use std::vec;

#[cfg(feature = "use-geo_0_31")]
use geo_0_31::Geometry;
use geo_traits::{CoordTrait, RectTrait};

use crate::error::Result;
use crate::indices::Indices;
use crate::r#type::IndexableNum;
#[cfg(feature = "use-geo_0_31")]
use crate::rtree::distance::DistanceMetric;
use crate::rtree::index::{RTree, RTreeRef};
use crate::rtree::traversal::{IntersectionIterator, Node};
use crate::rtree::util::upper_bound;
use crate::rtree::RTreeMetadata;
use crate::GeoIndexError;

/// A simple distance metric trait that doesn't depend on geo.
///
/// This trait is used for basic distance calculations without geometry support.
pub trait SimpleDistanceMetric<N: IndexableNum> {
    /// Calculate the distance between two points (x1, y1) and (x2, y2).
    fn distance(&self, x1: N, y1: N, x2: N, y2: N) -> N;

    /// Calculate the distance from a point to a bounding box.
    fn distance_to_bbox(&self, x: N, y: N, min_x: N, min_y: N, max_x: N, max_y: N) -> N;

    /// Return the maximum distance value for this metric.
    fn max_distance(&self) -> N {
        N::max_value()
    }
}

/// A trait for accessing geometries by index.
///
/// This trait allows different storage strategies for geometries (direct storage,
/// WKB decoding, caching, etc.) to be used with spatial index queries.
#[cfg(feature = "use-geo_0_31")]
pub trait GeometryAccessor {
    /// Get the geometry at the given index.
    ///
    /// # Arguments
    /// * `item_index` - Index of the item to retrieve
    ///
    /// # Returns
    /// A reference to the geometry at the given index, or None if the index is out of bounds
    fn get_geometry(&self, item_index: usize) -> Option<&Geometry<f64>>;
}

/// A trait for searching and accessing data out of an RTree.
pub trait RTreeIndex<N: IndexableNum>: Sized {
    /// A slice representing all the bounding boxes of all elements contained within this tree,
    /// including the bounding boxes of each internal node.
    fn boxes(&self) -> &[N];

    /// A slice representing the indices within the `boxes` slice, including internal nodes.
    fn indices(&self) -> Indices<'_>;

    /// Access the metadata describing this RTree
    fn metadata(&self) -> &RTreeMetadata<N>;

    /// The total number of items contained in this RTree.
    fn num_items(&self) -> u32 {
        self.metadata().num_items()
    }

    /// The total number of nodes in this RTree, including both leaf and intermediate nodes.
    fn num_nodes(&self) -> usize {
        self.metadata().num_nodes()
    }

    /// The maximum number of elements in each node.
    fn node_size(&self) -> u16 {
        self.metadata().node_size()
    }

    /// The offsets into [RTreeIndex::boxes] where each level's boxes starts and ends. The tree is
    /// laid out bottom-up, and there's an implicit initial 0. So the boxes of the lowest level of
    /// the tree are located from `boxes[0..self.level_bounds()[0]]`.
    fn level_bounds(&self) -> &[usize] {
        self.metadata().level_bounds()
    }

    /// The number of levels (height) of the tree.
    fn num_levels(&self) -> usize {
        self.level_bounds().len()
    }

    /// The tree is laid out from bottom to top. Level 0 is the _base_ of the tree. Each integer
    /// higher is one level higher of the tree.
    fn boxes_at_level(&self, level: usize) -> Result<&[N]> {
        let level_bounds = self.level_bounds();
        if level >= level_bounds.len() {
            return Err(GeoIndexError::General("Level out of bounds".to_string()));
        }
        let result = if level == 0 {
            &self.boxes()[0..level_bounds[0]]
        } else if level == level_bounds.len() {
            &self.boxes()[level_bounds[level]..]
        } else {
            &self.boxes()[level_bounds[level - 1]..level_bounds[level]]
        };
        Ok(result)
    }

    /// Search an RTree given the provided bounding box.
    ///
    /// Results are the indexes of the inserted objects in insertion order.
    fn search(&self, min_x: N, min_y: N, max_x: N, max_y: N) -> Vec<u32> {
        let boxes = self.boxes();
        let indices = self.indices();
        if boxes.is_empty() {
            return vec![];
        }

        let mut outer_node_index = boxes.len().checked_sub(4);

        let mut queue = VecDeque::with_capacity(self.node_size() as usize);
        let mut results = vec![];

        while let Some(node_index) = outer_node_index {
            // find the end index of the node
            let end = (node_index + self.node_size() as usize * 4)
                .min(upper_bound(node_index, self.level_bounds()));

            // search through child nodes
            for pos in (node_index..end).step_by(4) {
                // Safety: pos was checked before to be within bounds
                // Justification: avoiding bounds check improves performance by up to 30%
                let (node_min_x, node_min_y, node_max_x, node_max_y) = unsafe {
                    let node_min_x = *boxes.get_unchecked(pos);
                    let node_min_y = *boxes.get_unchecked(pos + 1);
                    let node_max_x = *boxes.get_unchecked(pos + 2);
                    let node_max_y = *boxes.get_unchecked(pos + 3);
                    (node_min_x, node_min_y, node_max_x, node_max_y)
                };

                // check if the query box disjoint with the node box
                if max_x < node_min_x
                    || max_y < node_min_y
                    || min_x > node_max_x
                    || min_y > node_max_y
                {
                    continue;
                }

                let index = indices.get(pos >> 2);

                if node_index >= self.num_items() as usize * 4 {
                    queue.push_back(index); // node; add it to the search queue
                } else {
                    // Since the max items of the index is u32, we can coerce to u32
                    results.push(index.try_into().unwrap()); // leaf item
                }
            }

            outer_node_index = queue.pop_front();
        }

        results
    }

    /// Search an RTree given the provided bounding box.
    ///
    /// Results are the indexes of the inserted objects in insertion order.
    fn search_rect(&self, rect: &impl RectTrait<T = N>) -> Vec<u32> {
        self.search(
            rect.min().x(),
            rect.min().y(),
            rect.max().x(),
            rect.max().y(),
        )
    }

    /// Search items in order of distance from the given point.
    ///
    /// This method uses Euclidean distance by default. For other distance metrics,
    /// use [`neighbors_with_distance`].
    ///
    /// ```
    /// use geo_index::rtree::{RTreeBuilder, RTreeIndex, RTreeRef};
    /// use geo_index::rtree::sort::HilbertSort;
    ///
    /// // Create an RTree
    /// let mut builder = RTreeBuilder::<f64>::new(3);
    /// builder.add(0., 0., 2., 2.);
    /// builder.add(1., 1., 3., 3.);
    /// builder.add(2., 2., 4., 4.);
    /// let tree = builder.finish::<HilbertSort>();
    ///
    /// let results = tree.neighbors(5., 5., None, None);
    /// assert_eq!(results, vec![2, 1, 0]);
    /// ```
    fn neighbors(
        &self,
        x: N,
        y: N,
        max_results: Option<usize>,
        max_distance: Option<N>,
    ) -> Vec<u32> {
        // Use simple squared distance for backward compatibility
        struct SimpleSquaredDistance;
        impl<N: IndexableNum> SimpleDistanceMetric<N> for SimpleSquaredDistance {
            fn distance(&self, x1: N, y1: N, x2: N, y2: N) -> N {
                let dx = x2 - x1;
                let dy = y2 - y1;
                dx * dx + dy * dy
            }
            fn distance_to_bbox(&self, x: N, y: N, min_x: N, min_y: N, max_x: N, max_y: N) -> N {
                let dx = axis_dist(x, min_x, max_x);
                let dy = axis_dist(y, min_y, max_y);
                dx * dx + dy * dy
            }
        }
        let simple_distance = SimpleSquaredDistance;
        self.neighbors_with_simple_distance(x, y, max_results, max_distance, &simple_distance)
    }

    /// Search items in order of distance from the given point using a simple distance metric.
    ///
    /// This is the base method for distance-based neighbor searches that works without the geo feature.
    fn neighbors_with_simple_distance<M: SimpleDistanceMetric<N> + ?Sized>(
        &self,
        x: N,
        y: N,
        max_results: Option<usize>,
        max_distance: Option<N>,
        distance_metric: &M,
    ) -> Vec<u32> {
        let boxes = self.boxes();
        if boxes.is_empty() {
            return vec![];
        }

        let indices = self.indices();
        let max_distance = max_distance.unwrap_or(distance_metric.max_distance());

        let mut outer_node_index = boxes.len().checked_sub(4);
        let mut queue = BinaryHeap::new();
        let mut results: Vec<u32> = vec![];

        'outer: while let Some(node_index) = outer_node_index {
            // find the end index of the node
            let end = (node_index + self.node_size() as usize * 4)
                .min(upper_bound(node_index, self.level_bounds()));

            // add child nodes to the queue
            for pos in (node_index..end).step_by(4) {
                let index = indices.get(pos >> 2);

                // Use the custom distance metric for bbox distance calculation
                let dist = distance_metric.distance_to_bbox(
                    x,
                    y,
                    boxes[pos],
                    boxes[pos + 1],
                    boxes[pos + 2],
                    boxes[pos + 3],
                );

                if dist > max_distance {
                    continue;
                }

                if node_index >= self.num_items() as usize * 4 {
                    // node (use even id)
                    queue.push(Reverse(NeighborNode {
                        id: index << 1,
                        dist,
                    }));
                } else {
                    // leaf item (use odd id)
                    // Use consistent distance calculation for both nodes and leaf items
                    queue.push(Reverse(NeighborNode {
                        id: (index << 1) + 1,
                        dist,
                    }));
                }
            }

            // pop items from the queue
            while !queue.is_empty() && queue.peek().is_some_and(|val| (val.0.id & 1) != 0) {
                let dist = queue.peek().unwrap().0.dist;
                if dist > max_distance {
                    break 'outer;
                }
                let item = queue.pop().unwrap();
                results.push((item.0.id >> 1).try_into().unwrap());
                if max_results.is_some_and(|max_results| results.len() == max_results) {
                    break 'outer;
                }
            }

            if let Some(item) = queue.pop() {
                outer_node_index = Some(item.0.id >> 1);
            } else {
                outer_node_index = None;
            }
        }

        results
    }

    /// Search items in order of distance from the given point using a custom distance metric.
    ///
    /// This method allows you to specify a custom distance calculation method, such as
    /// Euclidean, Haversine, or Spheroid distance.
    ///
    /// ```
    /// use geo_index::rtree::{RTreeBuilder, RTreeIndex};
    /// use geo_index::rtree::distance::{EuclideanDistance, HaversineDistance};
    /// use geo_index::rtree::sort::HilbertSort;
    ///
    /// // Create an RTree with geographic coordinates (longitude, latitude)
    /// let mut builder = RTreeBuilder::<f64>::new(3);
    /// builder.add(-74.0, 40.7, -74.0, 40.7); // New York
    /// builder.add(-0.1, 51.5, -0.1, 51.5);   // London
    /// builder.add(139.7, 35.7, 139.7, 35.7); // Tokyo
    /// let tree = builder.finish::<HilbertSort>();
    ///
    /// // Find nearest neighbors using Haversine distance (great-circle distance)
    /// let haversine = HaversineDistance::default();
    /// let results = tree.neighbors_with_distance(-74.0, 40.7, Some(2), None, &haversine);
    /// ```
    #[cfg(feature = "use-geo_0_31")]
    fn neighbors_with_distance<M: DistanceMetric<N> + ?Sized>(
        &self,
        x: N,
        y: N,
        max_results: Option<usize>,
        max_distance: Option<N>,
        distance_metric: &M,
    ) -> Vec<u32> {
        self.neighbors_with_simple_distance(x, y, max_results, max_distance, distance_metric)
    }

    /// Search items in order of distance from the given coordinate.
    fn neighbors_coord(
        &self,
        coord: &impl CoordTrait<T = N>,
        max_results: Option<usize>,
        max_distance: Option<N>,
    ) -> Vec<u32> {
        self.neighbors(coord.x(), coord.y(), max_results, max_distance)
    }

    /// Search items in order of distance from the given coordinate using a custom distance metric.
    #[cfg(feature = "use-geo_0_31")]
    fn neighbors_coord_with_distance<M: DistanceMetric<N> + ?Sized>(
        &self,
        coord: &impl CoordTrait<T = N>,
        max_results: Option<usize>,
        max_distance: Option<N>,
        distance_metric: &M,
    ) -> Vec<u32> {
        self.neighbors_with_distance(
            coord.x(),
            coord.y(),
            max_results,
            max_distance,
            distance_metric,
        )
    }

    /// Search items in order of distance from a query geometry using a distance metric and geometry accessor.
    ///
    /// This method allows searching with geometry-to-geometry distance calculations.
    /// The distance metric defines how distances are computed, and the geometry accessor
    /// provides access to the actual geometries by index.
    ///
    /// ```
    /// use geo_index::rtree::{RTreeBuilder, RTreeIndex};
    /// use geo_index::rtree::distance::{EuclideanDistance, SliceGeometryAccessor};
    /// use geo_index::rtree::sort::HilbertSort;
    /// use geo_0_31::{Point, Geometry};
    ///
    /// // Create an RTree
    /// let mut builder = RTreeBuilder::<f64>::new(3);
    /// builder.add(0., 0., 2., 2.);
    /// builder.add(5., 5., 7., 7.);
    /// builder.add(10., 10., 12., 12.);
    /// let tree = builder.finish::<HilbertSort>();
    ///
    /// // Example geometries
    /// let geometries: Vec<Geometry<f64>> = vec![
    ///     Geometry::Point(Point::new(1.0, 1.0)),
    ///     Geometry::Point(Point::new(6.0, 6.0)),
    ///     Geometry::Point(Point::new(11.0, 11.0)),
    /// ];
    ///
    /// let metric = EuclideanDistance;
    /// let accessor = SliceGeometryAccessor::new(&geometries);
    /// let query_geom = Geometry::Point(Point::new(3.0, 3.0));
    /// let results = tree.neighbors_geometry(&query_geom, None, None, &metric, &accessor);
    /// ```
    #[cfg(feature = "use-geo_0_31")]
    fn neighbors_geometry<M: DistanceMetric<N> + ?Sized, A: GeometryAccessor + ?Sized>(
        &self,
        query_geometry: &Geometry<f64>,
        max_results: Option<usize>,
        max_distance: Option<N>,
        distance_metric: &M,
        accessor: &A,
    ) -> Vec<u32> {
        let boxes = self.boxes();
        if boxes.is_empty() {
            return vec![];
        }

        let indices = self.indices();
        let max_distance = max_distance.unwrap_or(distance_metric.max_distance());

        let mut outer_node_index = boxes.len().checked_sub(4);
        let mut queue = BinaryHeap::new();
        let mut results: Vec<u32> = vec![];

        'outer: while let Some(node_index) = outer_node_index {
            // find the end index of the node
            let end = (node_index + self.node_size() as usize * 4)
                .min(upper_bound(node_index, self.level_bounds()));

            // add child nodes to the queue
            for pos in (node_index..end).step_by(4) {
                let index = indices.get(pos >> 2);

                let dist = if node_index >= self.num_items() as usize * 4 {
                    // For internal nodes, use geometry-to-bbox distance
                    distance_metric.distance_geometry_to_bbox(
                        query_geometry,
                        boxes[pos],
                        boxes[pos + 1],
                        boxes[pos + 2],
                        boxes[pos + 3],
                    )
                } else {
                    // For leaf items, use geometry-to-geometry distance
                    if let Some(item_geom) = accessor.get_geometry(index) {
                        distance_metric.distance_to_geometry(query_geometry, item_geom)
                    } else {
                        distance_metric.max_distance()
                    }
                };

                if dist > max_distance {
                    continue;
                }

                if node_index >= self.num_items() as usize * 4 {
                    // node (use even id)
                    queue.push(Reverse(NeighborNode {
                        id: index << 1,
                        dist,
                    }));
                } else {
                    // leaf item (use odd id)
                    queue.push(Reverse(NeighborNode {
                        id: (index << 1) + 1,
                        dist,
                    }));
                }
            }

            // pop items from the queue
            while !queue.is_empty() && queue.peek().is_some_and(|val| (val.0.id & 1) != 0) {
                let dist = queue.peek().unwrap().0.dist;
                if dist > max_distance {
                    break 'outer;
                }
                let item = queue.pop().unwrap();
                results.push((item.0.id >> 1).try_into().unwrap());
                if max_results.is_some_and(|max_results| results.len() == max_results) {
                    break 'outer;
                }
            }

            if let Some(item) = queue.pop() {
                outer_node_index = Some(item.0.id >> 1);
            } else {
                outer_node_index = None;
            }
        }

        results
    }

    /// Returns an iterator over the indexes of objects in this and another tree that intersect.
    ///
    /// Each returned object is of the form `(u32, u32)`, where the first is the positional
    /// index of the "left" tree and the second is the index of the "right" tree.
    fn intersection_candidates_with_other_tree<'a>(
        &'a self,
        other: &'a impl RTreeIndex<N>,
    ) -> impl Iterator<Item = (u32, u32)> + 'a {
        IntersectionIterator::from_trees(self, other)
    }

    /// Access the root node of the RTree for manual traversal.
    fn root(&self) -> Node<'_, N, Self> {
        Node::from_root(self)
    }
}

/// A wrapper around a node and its distance for use in the priority queue.
#[derive(Debug, Clone, Copy, PartialEq)]
struct NeighborNode<N: IndexableNum> {
    id: usize,
    dist: N,
}

impl<N: IndexableNum> Eq for NeighborNode<N> {}

impl<N: IndexableNum> Ord for NeighborNode<N> {
    fn cmp(&self, other: &Self) -> std::cmp::Ordering {
        // We don't allow NaN. This should only panic on NaN
        self.dist.partial_cmp(&other.dist).unwrap()
    }
}

impl<N: IndexableNum> PartialOrd for NeighborNode<N> {
    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
        Some(self.cmp(other))
    }
}

impl<N: IndexableNum> RTreeIndex<N> for RTree<N> {
    fn boxes(&self) -> &[N] {
        self.metadata.boxes_slice(&self.buffer)
    }

    fn indices(&self) -> Indices<'_> {
        self.metadata.indices_slice(&self.buffer)
    }

    fn metadata(&self) -> &RTreeMetadata<N> {
        &self.metadata
    }
}

impl<N: IndexableNum> RTreeIndex<N> for RTreeRef<'_, N> {
    fn boxes(&self) -> &[N] {
        self.boxes
    }

    fn indices(&self) -> Indices<'_> {
        self.indices
    }

    fn metadata(&self) -> &RTreeMetadata<N> {
        &self.metadata
    }
}

/// 1D distance from a value to a range.
#[inline]
pub(crate) fn axis_dist<N: IndexableNum>(k: N, min: N, max: N) -> N {
    if k < min {
        min - k
    } else if k <= max {
        N::zero()
    } else {
        k - max
    }
}

#[cfg(test)]
mod test {
    // Replication of tests from flatbush js
    mod js {
        use crate::rtree::RTreeIndex;
        use crate::test::{flatbush_js_test_data, flatbush_js_test_index};

        #[test]
        fn performs_bbox_search() {
            let data = flatbush_js_test_data();
            let index = flatbush_js_test_index();
            let ids = index.search(40., 40., 60., 60.);

            let mut results: Vec<usize> = vec![];
            for id in ids {
                results.push(data[4 * id as usize] as usize);
                results.push(data[4 * id as usize + 1] as usize);
                results.push(data[4 * id as usize + 2] as usize);
                results.push(data[4 * id as usize + 3] as usize);
            }

            results.sort();

            let mut expected = vec![
                57, 59, 58, 59, 48, 53, 52, 56, 40, 42, 43, 43, 43, 41, 47, 43,
            ];
            expected.sort();

            assert_eq!(results, expected);
        }
    }
    #[cfg(feature = "use-geo_0_31")]
    mod distance_metrics {
        use crate::rtree::distance::{EuclideanDistance, HaversineDistance};
        use crate::rtree::r#trait::SimpleDistanceMetric;
        use crate::rtree::sort::HilbertSort;
        use crate::rtree::{RTreeBuilder, RTreeIndex};

        #[test]
        fn test_euclidean_distance_neighbors() {
            let mut builder = RTreeBuilder::<f64>::new(3);
            builder.add(0., 0., 1., 1.);
            builder.add(2., 2., 3., 3.);
            builder.add(4., 4., 5., 5.);
            let tree = builder.finish::<HilbertSort>();

            let euclidean = EuclideanDistance;
            let results = tree.neighbors_with_distance(0., 0., None, None, &euclidean);

            // Should return items in order of distance from (0,0)
            assert_eq!(results, vec![0, 1, 2]);
        }

        #[test]
        fn test_haversine_distance_neighbors() {
            let mut builder = RTreeBuilder::<f64>::new(3);
            // Add some geographic points (longitude, latitude)
            builder.add(-74.0, 40.7, -74.0, 40.7); // New York
            builder.add(-0.1, 51.5, -0.1, 51.5); // London
            builder.add(139.7, 35.7, 139.7, 35.7); // Tokyo
            let tree = builder.finish::<HilbertSort>();

            let haversine = HaversineDistance::default();
            let results = tree.neighbors_with_distance(-74.0, 40.7, None, None, &haversine);

            // From New York, should find New York first, then London, then Tokyo
            assert_eq!(results, vec![0, 1, 2]);
        }

        #[test]
        fn test_backward_compatibility() {
            let mut builder = RTreeBuilder::<f64>::new(3);
            builder.add(0., 0., 1., 1.);
            builder.add(2., 2., 3., 3.);
            builder.add(4., 4., 5., 5.);
            let tree = builder.finish::<HilbertSort>();

            // Test that original neighbors method still works
            let results_original = tree.neighbors(0., 0., None, None);

            // Test that new method with Euclidean distance gives same results
            let euclidean = EuclideanDistance;
            let results_new = tree.neighbors_with_distance(0., 0., None, None, &euclidean);

            assert_eq!(results_original, results_new);
        }

        #[test]
        fn test_max_distance_filtering() {
            let mut builder = RTreeBuilder::<f64>::new(3);
            builder.add(0., 0., 1., 1.);
            builder.add(2., 2., 3., 3.);
            builder.add(10., 10., 11., 11.);
            let tree = builder.finish::<HilbertSort>();

            let euclidean = EuclideanDistance;
            // Only find neighbors within distance 5
            let results = tree.neighbors_with_distance(0., 0., None, Some(5.0), &euclidean);

            // Should only find first two items, not the distant third one
            assert_eq!(results.len(), 2);
            assert_eq!(results, vec![0, 1]);
        }

        #[test]
        #[cfg(feature = "use-geo_0_31")]
        fn test_geometry_neighbors_euclidean() {
            use crate::r#type::IndexableNum;
            use crate::rtree::distance::{DistanceMetric, SliceGeometryAccessor};
            use geo_0_31::algorithm::{Distance, Euclidean};
            use geo_0_31::{Geometry, Point};

            let mut builder = RTreeBuilder::<f64>::new(3);
            builder.add(0., 0., 2., 2.); // Item 0
            builder.add(5., 5., 7., 7.); // Item 1
            builder.add(10., 10., 12., 12.); // Item 2
            let tree = builder.finish::<HilbertSort>();

            // Geometries corresponding to the bboxes
            let geometries: Vec<Geometry<f64>> = vec![
                Geometry::Point(Point::new(1.0, 1.0)),   // Item 0
                Geometry::Point(Point::new(6.0, 6.0)),   // Item 1
                Geometry::Point(Point::new(11.0, 11.0)), // Item 2
            ];

            struct SimpleMetric;
            impl<N: IndexableNum> SimpleDistanceMetric<N> for SimpleMetric {
                fn distance(&self, x1: N, y1: N, x2: N, y2: N) -> N {
                    let dx = x2 - x1;
                    let dy = y2 - y1;
                    (dx * dx + dy * dy).sqrt().unwrap_or(N::max_value())
                }
                fn distance_to_bbox(
                    &self,
                    x: N,
                    y: N,
                    min_x: N,
                    min_y: N,
                    max_x: N,
                    max_y: N,
                ) -> N {
                    let dx = if x < min_x {
                        min_x - x
                    } else if x > max_x {
                        x - max_x
                    } else {
                        N::zero()
                    };
                    let dy = if y < min_y {
                        min_y - y
                    } else if y > max_y {
                        y - max_y
                    } else {
                        N::zero()
                    };
                    (dx * dx + dy * dy).sqrt().unwrap_or(N::max_value())
                }
            }
            impl<N: IndexableNum> DistanceMetric<N> for SimpleMetric {
                fn distance_to_geometry(&self, geom1: &Geometry<f64>, geom2: &Geometry<f64>) -> N {
                    N::from_f64(Euclidean.distance(geom1, geom2)).unwrap_or(N::max_value())
                }
            }

            let query_geom = Geometry::Point(Point::new(3.0, 3.0));
            let metric = SimpleMetric;
            let accessor = SliceGeometryAccessor::new(&geometries);
            let results = tree.neighbors_geometry(&query_geom, None, None, &metric, &accessor);

            // Item 0 should be closest to query point (3,3)
            assert_eq!(results[0], 0);
            assert_eq!(results[1], 1);
            assert_eq!(results[2], 2);
        }

        #[test]
        #[cfg(feature = "use-geo_0_31")]
        fn test_geometry_neighbors_linestring() {
            use crate::r#type::IndexableNum;
            use crate::rtree::distance::{DistanceMetric, SliceGeometryAccessor};
            use geo_0_31::algorithm::{Distance, Euclidean};
            use geo_0_31::{coord, Geometry, LineString, Point};

            let mut builder = RTreeBuilder::<f64>::new(3);
            builder.add(0., 0., 10., 0.); // Item 0 - horizontal line
            builder.add(5., 5., 15., 5.); // Item 1 - horizontal line higher up
            builder.add(0., 10., 10., 10.); // Item 2 - horizontal line at top
            let tree = builder.finish::<HilbertSort>();

            // Geometries corresponding to the bboxes
            let geometries: Vec<Geometry<f64>> = vec![
                Geometry::LineString(LineString::new(vec![
                    coord! { x: 0.0, y: 0.0 },
                    coord! { x: 10.0, y: 0.0 },
                ])),
                Geometry::LineString(LineString::new(vec![
                    coord! { x: 5.0, y: 5.0 },
                    coord! { x: 15.0, y: 5.0 },
                ])),
                Geometry::LineString(LineString::new(vec![
                    coord! { x: 0.0, y: 10.0 },
                    coord! { x: 10.0, y: 10.0 },
                ])),
            ];

            struct SimpleMetric;
            impl<N: IndexableNum> SimpleDistanceMetric<N> for SimpleMetric {
                fn distance(&self, x1: N, y1: N, x2: N, y2: N) -> N {
                    let dx = x2 - x1;
                    let dy = y2 - y1;
                    (dx * dx + dy * dy).sqrt().unwrap_or(N::max_value())
                }
                fn distance_to_bbox(
                    &self,
                    x: N,
                    y: N,
                    min_x: N,
                    min_y: N,
                    max_x: N,
                    max_y: N,
                ) -> N {
                    let dx = if x < min_x {
                        min_x - x
                    } else if x > max_x {
                        x - max_x
                    } else {
                        N::zero()
                    };
                    let dy = if y < min_y {
                        min_y - y
                    } else if y > max_y {
                        y - max_y
                    } else {
                        N::zero()
                    };
                    (dx * dx + dy * dy).sqrt().unwrap_or(N::max_value())
                }
            }
            impl<N: IndexableNum> DistanceMetric<N> for SimpleMetric {
                fn distance_to_geometry(&self, geom1: &Geometry<f64>, geom2: &Geometry<f64>) -> N {
                    N::from_f64(Euclidean.distance(geom1, geom2)).unwrap_or(N::max_value())
                }
            }

            let query_geom = Geometry::Point(Point::new(5.0, 2.0));
            let metric = SimpleMetric;
            let accessor = SliceGeometryAccessor::new(&geometries);
            let results = tree.neighbors_geometry(&query_geom, None, None, &metric, &accessor);

            // Item 0 (bottom line) should be closest to point (5, 2)
            assert_eq!(results[0], 0);
        }

        #[test]
        #[cfg(feature = "use-geo_0_31")]
        fn test_geometry_neighbors_with_max_results() {
            use crate::r#type::IndexableNum;
            use crate::rtree::distance::{DistanceMetric, SliceGeometryAccessor};
            use geo_0_31::algorithm::{Distance, Euclidean};
            use geo_0_31::{Geometry, Point};

            let mut builder = RTreeBuilder::<f64>::new(5);
            for i in 0..5 {
                let x = (i * 3) as f64;
                builder.add(x, x, x + 1., x + 1.);
            }
            let tree = builder.finish::<HilbertSort>();

            // Create geometries for each bbox
            let geometries: Vec<Geometry<f64>> = (0..5)
                .map(|i| {
                    let x = (i * 3) as f64;
                    Geometry::Point(Point::new(x + 0.5, x + 0.5))
                })
                .collect();

            struct SimpleMetric;
            impl<N: IndexableNum> SimpleDistanceMetric<N> for SimpleMetric {
                fn distance(&self, x1: N, y1: N, x2: N, y2: N) -> N {
                    let dx = x2 - x1;
                    let dy = y2 - y1;
                    (dx * dx + dy * dy).sqrt().unwrap_or(N::max_value())
                }
                fn distance_to_bbox(
                    &self,
                    x: N,
                    y: N,
                    min_x: N,
                    min_y: N,
                    max_x: N,
                    max_y: N,
                ) -> N {
                    let dx = if x < min_x {
                        min_x - x
                    } else if x > max_x {
                        x - max_x
                    } else {
                        N::zero()
                    };
                    let dy = if y < min_y {
                        min_y - y
                    } else if y > max_y {
                        y - max_y
                    } else {
                        N::zero()
                    };
                    (dx * dx + dy * dy).sqrt().unwrap_or(N::max_value())
                }
            }
            impl<N: IndexableNum> DistanceMetric<N> for SimpleMetric {
                fn distance_to_geometry(&self, geom1: &Geometry<f64>, geom2: &Geometry<f64>) -> N {
                    N::from_f64(Euclidean.distance(geom1, geom2)).unwrap_or(N::max_value())
                }
            }

            let query_geom = Geometry::Point(Point::new(5.0, 5.0));
            let metric = SimpleMetric;
            let accessor = SliceGeometryAccessor::new(&geometries);
            let results = tree.neighbors_geometry(&query_geom, Some(3), None, &metric, &accessor);

            assert_eq!(results.len(), 3);
            // Should get the 3 closest items
        }

        #[test]
        #[cfg(feature = "use-geo_0_31")]
        fn test_geometry_neighbors_haversine() {
            use crate::r#type::IndexableNum;
            use crate::rtree::distance::{DistanceMetric, SliceGeometryAccessor};
            use geo_0_31::algorithm::{Centroid, Distance, Haversine};
            use geo_0_31::{Geometry, Point};

            let mut builder = RTreeBuilder::<f64>::new(3);
            // Geographic bounding boxes (lon, lat)
            builder.add(-74.1, 40.6, -74.0, 40.7); // New York area
            builder.add(-0.2, 51.4, -0.1, 51.5); // London area
            builder.add(139.6, 35.6, 139.7, 35.7); // Tokyo area
            let tree = builder.finish::<HilbertSort>();

            let geometries: Vec<Geometry<f64>> = vec![
                Geometry::Point(Point::new(-74.0, 40.7)), // New York
                Geometry::Point(Point::new(-0.1, 51.5)),  // London
                Geometry::Point(Point::new(139.7, 35.7)), // Tokyo
            ];

            struct HaversineMetric;
            impl<N: IndexableNum> SimpleDistanceMetric<N> for HaversineMetric {
                fn distance(&self, lon1: N, lat1: N, lon2: N, lat2: N) -> N {
                    let p1 = Point::new(lon1.to_f64().unwrap_or(0.0), lat1.to_f64().unwrap_or(0.0));
                    let p2 = Point::new(lon2.to_f64().unwrap_or(0.0), lat2.to_f64().unwrap_or(0.0));
                    N::from_f64(Haversine.distance(p1, p2)).unwrap_or(N::max_value())
                }
                fn distance_to_bbox(
                    &self,
                    lon: N,
                    lat: N,
                    min_lon: N,
                    min_lat: N,
                    max_lon: N,
                    max_lat: N,
                ) -> N {
                    let lon_f = lon.to_f64().unwrap_or(0.0);
                    let lat_f = lat.to_f64().unwrap_or(0.0);
                    let min_lon_f = min_lon.to_f64().unwrap_or(0.0);
                    let min_lat_f = min_lat.to_f64().unwrap_or(0.0);
                    let max_lon_f = max_lon.to_f64().unwrap_or(0.0);
                    let max_lat_f = max_lat.to_f64().unwrap_or(0.0);
                    let closest_lon = lon_f.clamp(min_lon_f, max_lon_f);
                    let closest_lat = lat_f.clamp(min_lat_f, max_lat_f);
                    let point = Point::new(lon_f, lat_f);
                    let closest_point = Point::new(closest_lon, closest_lat);
                    N::from_f64(Haversine.distance(point, closest_point)).unwrap_or(N::max_value())
                }
            }
            impl<N: IndexableNum> DistanceMetric<N> for HaversineMetric {
                fn distance_to_geometry(&self, geom1: &Geometry<f64>, geom2: &Geometry<f64>) -> N {
                    let c1 = geom1.centroid().unwrap_or(Point::new(0.0, 0.0));
                    let c2 = geom2.centroid().unwrap_or(Point::new(0.0, 0.0));
                    N::from_f64(Haversine.distance(c1, c2)).unwrap_or(N::max_value())
                }
            }

            let query_geom = Geometry::Point(Point::new(-74.0, 40.7)); // New York
            let metric = HaversineMetric;
            let accessor = SliceGeometryAccessor::new(&geometries);
            let results = tree.neighbors_geometry(&query_geom, None, None, &metric, &accessor);

            // New York should be closest (distance 0)
            assert_eq!(results[0], 0);
        }
    }
}