Fix numerical bug in point query

Author: radevgit

CHANGELOG.md

@@ -1,5 +1,10 @@
 # Changelog
 
+
+## [0.6.10] - 2025-11-01
+- Fixed bug query_circle_points() and query_nearest_k_points() for very small r
+- Added get_point() 
+
 ## [0.6.9] - 2025-11-01
 - New queries query_circle_points() and query_nearest_k_points() add_point() 
 

Cargo.lock

@@ -1,7 +1,7 @@
 [package]
 
 name = "aabb"
-version = "0.6.9"
+version = "0.6.10"
 description = "Static AABB spatial index for 2D queries"
 rust-version = "1.88"
 edition = "2024"

Cargo.toml

@@ -92,9 +92,10 @@ The Hilbert R-tree stores bounding boxes in a flat array and sorts them by their
 - `HilbertRTreeI32::new()` or `AABBI32::new()` - Create a new empty tree
 - `HilbertRTreeI32::with_capacity(capacity)` or `AABBI32::with_capacity(capacity)` - Create a new tree with preallocated capacity
 - `add(min_x, min_y, max_x, max_y)` - `(f64, i32)` Add a bounding box
-- `add_point(x, y)` - `(f64)` Add a point (convenience method - internally stores as (x, x, y, y))
+- `add_point(x, y)` - `(f64)` Add a point (convenience method - internally stores as (x, y, x, y))
 - `build()` - `(f64, i32)` Build the spatial index (required before querying)
 - `get(item_id)` - `(f64, i32)` Retrieve the bounding box for an item by its ID
+- `get_point(item_id)` - `(f64)` Retrieve a point as (x, y) for items added with `add_point()`
 - `save(path)` - `(f64, i32)` Save the built tree to a file for fast loading later
 - `load(path)` - `(f64, i32)` Load a previously saved tree from a file
 

README.md

@@ -51,6 +51,10 @@ pub struct HilbertRTree {
     total_nodes: usize,
     /// Pre-allocated capacity in bytes (0 if not pre-allocated)
     allocated_capacity: usize,
+    /// Inverse mapping: original item ID -> position in sorted tree
+    /// After build(), sorted_order[original_id] = position in memory
+    /// This allows get(item_id) to find the correct box despite Hilbert sorting
+    pub(crate) sorted_order: Vec<usize>,
 }
 
 const MAX_HILBERT: u32 = u16::MAX as u32;
@@ -119,6 +123,7 @@ impl HilbertRTree {
             position: 0,
             bounds: Box::new(f64::INFINITY, f64::INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY),
             total_nodes: 0,
+            sorted_order: Vec::new(),
         }
     }
 
@@ -181,7 +186,31 @@ impl HilbertRTree {
     /// // Results contain points 0 and 1 within distance 1.5
     /// ```
     pub fn add_point(&mut self, x: f64, y: f64) {
-        self.add(x, x, y, y);
+        self.add(x, y, x, y);
+    }
+
+    /// Retrieves a point that was added with `add_point()`.
+    ///
+    /// This is a convenience method for retrieving point data. Returns `Some((x, y))` if the
+    /// point exists, or `None` if the index is out of bounds.
+    ///
+    /// # Arguments
+    /// * `item_id` - The index of the point (0-based, in order added)
+    ///
+    /// # Example
+    /// ```
+    /// use aabb::prelude::*;
+    /// let mut tree = AABB::with_capacity(2);
+    /// tree.add_point(1.5, 2.5);
+    /// tree.add_point(3.0, 4.0);
+    /// tree.build();
+    ///
+    /// assert_eq!(tree.get_point(0), Some((1.5, 2.5)));
+    /// assert_eq!(tree.get_point(1), Some((3.0, 4.0)));
+    /// assert_eq!(tree.get_point(2), None);
+    /// ```
+    pub fn get_point(&self, item_id: usize) -> Option<(f64, f64)> {
+        self.get(item_id).map(|(min_x, min_y, _max_x, _max_y)| (min_x, min_y))
     }
 
     /// Builds the Hilbert R-tree index
@@ -260,6 +289,9 @@ impl HilbertRTree {
             unsafe {
                 std::ptr::write_unaligned(root_idx_ptr, 0_u32 << 2_u32); // First child at position 0
             }
+            
+            // For single-node case, no sorting happens, so sorted_order is identity mapping
+            self.sorted_order = (0..num_items).collect();
             return;
         }
 
@@ -315,6 +347,16 @@ impl HilbertRTree {
 
         // Initialize leaf indices AFTER sorting - map new position to original box ID
         let indices_start = HEADER_SIZE + total_nodes * size_of::<Box>();
+        
+        // Build inverse mapping: original_id -> sorted_position
+        // This allows get(item_id) to find the correct box despite Hilbert sorting
+        let mut sorted_order = vec![0usize; num_items];
+        for sorted_pos in 0..num_items {
+            let original_id = sort_indices[sorted_pos];
+            sorted_order[original_id] = sorted_pos;
+        }
+        self.sorted_order = sorted_order;
+        
         for i in 0..num_items {
             let idx_ptr = &mut self.data[indices_start + i * size_of::<u32>()] as *mut u8 as *mut u32;
             unsafe {
@@ -1739,7 +1781,16 @@ impl HilbertRTree {
             return None;
         }
 
-        let box_data = self.get_box(item_id);
+        // After build(), boxes are sorted by Hilbert curve order
+        // sorted_order[item_id] gives us the position of this item in the sorted tree
+        let pos = if self.sorted_order.is_empty() {
+            // If tree hasn't been built yet, items are in original order
+            item_id
+        } else {
+            self.sorted_order[item_id]
+        };
+        
+        let box_data = self.get_box(pos);
         Some((box_data.min_x, box_data.min_y, box_data.max_x, box_data.max_y))
     }
 
@@ -1948,6 +1999,7 @@ impl HilbertRTree {
             bounds,
             total_nodes,
             allocated_capacity: data_len,
+            sorted_order: Vec::new(),
         })
     }
 }

src/hilbert_rtree.rs

@@ -322,4 +322,127 @@ mod test_point_optimizations {
         assert!(results[0] != 999);
         assert!(results[1] != 888);
     }
+
+    /// Test get_point() convenience method
+    #[test]
+    fn test_get_point_method() {
+        let mut tree = AABB::with_capacity(4);
+        tree.add_point(0.0, 0.0);
+        tree.add_point(1.5, 2.5);
+        tree.add_point(3.0, 4.0);
+        tree.add_point(-1.0, -1.0);
+        tree.build();
+
+        // Verify get_point returns correct coordinates
+        assert_eq!(tree.get_point(0), Some((0.0, 0.0)));
+        assert_eq!(tree.get_point(1), Some((1.5, 2.5)));
+        assert_eq!(tree.get_point(2), Some((3.0, 4.0)));
+        assert_eq!(tree.get_point(3), Some((-1.0, -1.0)));
+        
+        // Out of bounds should return None
+        assert_eq!(tree.get_point(4), None);
+    }
+
+    /// Test that get_point works before build
+    #[test]
+    fn test_get_point_before_build() {
+        let mut tree = AABB::with_capacity(2);
+        tree.add_point(1.0, 2.0);
+        tree.add_point(3.0, 4.0);
+        
+        // Should work even before build()
+        assert_eq!(tree.get_point(0), Some((1.0, 2.0)));
+        assert_eq!(tree.get_point(1), Some((3.0, 4.0)));
+    }
+
+    /// Test query with very small radius (1e-9) - critical for exact point matching
+    #[test]
+    fn test_query_circle_points_very_small_radius() {
+        let mut tree = AABB::with_capacity(5);
+        tree.add_point(0.0, 0.0);
+        tree.add_point(1.0, 0.0);
+        tree.add_point(0.0, 1.0);
+        tree.add_point(1.0, 1.0);
+        tree.add_point(0.5, 0.5);
+        tree.build();
+
+        // Query with radius 1e-9 should only return exact point
+        let mut results = Vec::new();
+        tree.query_circle_points(0.0, 0.0, 1e-9, &mut results);
+        assert_eq!(results, vec![0], "Query with radius 1e-9 should only return exact point");
+        
+        // Query with radius 1e-6 should still only return exact point
+        tree.query_circle_points(0.0, 0.0, 1e-6, &mut results);
+        assert_eq!(results, vec![0], "Query with radius 1e-6 should only return exact point");
+    }
+
+    /// Test point storage correctness after build with Hilbert sorting
+    #[test]
+    fn test_point_storage_after_build() {
+        let mut tree = AABB::with_capacity(5);
+        
+        // Add points in specific order
+        tree.add_point(0.0, 0.0);
+        tree.add_point(1.0, 0.0);
+        tree.add_point(0.0, 1.0);
+        tree.add_point(1.0, 1.0);
+        tree.add_point(0.5, 0.5);
+        
+        tree.build();
+        
+        // Verify all points are stored correctly
+        assert_eq!(tree.get_point(0), Some((0.0, 0.0)));
+        assert_eq!(tree.get_point(1), Some((1.0, 0.0)));
+        assert_eq!(tree.get_point(2), Some((0.0, 1.0)));
+        assert_eq!(tree.get_point(3), Some((1.0, 1.0)));
+        assert_eq!(tree.get_point(4), Some((0.5, 0.5)));
+    }
+
+    /// Test that add_point stores correct coordinates (not swapped)
+    #[test]
+    fn test_add_point_coordinate_order() {
+        let mut tree = AABB::with_capacity(3);
+        
+        tree.add_point(1.5, 2.5);
+        tree.add_point(3.0, 4.0);
+        tree.add_point(-1.0, -2.0);
+        
+        tree.build();
+        
+        // Verify X and Y are not swapped
+        assert_eq!(tree.get_point(0), Some((1.5, 2.5)));
+        assert_eq!(tree.get_point(1), Some((3.0, 4.0)));
+        assert_eq!(tree.get_point(2), Some((-1.0, -2.0)));
+        
+        // Verify via get() as well
+        assert_eq!(tree.get(0).unwrap(), (1.5, 2.5, 1.5, 2.5));
+        assert_eq!(tree.get(1).unwrap(), (3.0, 4.0, 3.0, 4.0));
+        assert_eq!(tree.get(2).unwrap(), (-1.0, -2.0, -1.0, -2.0));
+    }
+
+    /// Test consistency between query_circle_points and get_point
+    #[test]
+    fn test_query_circle_points_consistency_with_get_point() {
+        let mut tree = AABB::with_capacity(10);
+        
+        // Create grid of points
+        for i in 0..3 {
+            for j in 0..3 {
+                tree.add_point(i as f64, j as f64);
+            }
+        }
+        tree.build();
+
+        // Query for points near (1.0, 1.0) with radius 1.5
+        let mut results = Vec::new();
+        tree.query_circle_points(1.0, 1.0, 1.5, &mut results);
+        
+        // For each result, verify get_point returns the correct coordinates
+        for &idx in &results {
+            if let Some((x, y)) = tree.get_point(idx) {
+                let distance = ((x - 1.0).powi(2) + (y - 1.0).powi(2)).sqrt();
+                assert!(distance <= 1.5, "Point {} at ({}, {}) has distance {}, exceeds radius 1.5", idx, x, y, distance);
+            }
+        }
+    }
 }