scirs2_spatial/rtree/
node.rs

1use crate::error::{SpatialError, SpatialResult};
2use scirs2_core::ndarray::{Array1, ArrayView1};
3use std::cmp::Ordering;
4// No imports from std::collections needed
5
6/// Rectangle represents a minimum bounding rectangle (MBR) in N-dimensional space.
7#[derive(Clone, Debug)]
8pub struct Rectangle {
9    /// Minimum coordinates in each dimension
10    pub(crate) min: Array1<f64>,
11    /// Maximum coordinates in each dimension
12    pub(crate) max: Array1<f64>,
13}
14
15impl Rectangle {
16    /// Create a new rectangle from min and max coordinates
17    ///
18    /// # Arguments
19    ///
20    /// * `min` - Minimum coordinates in each dimension
21    /// * `max` - Maximum coordinates in each dimension
22    ///
23    /// # Returns
24    ///
25    /// A `SpatialResult` containing the rectangle if valid, or an error if the input is invalid
26    pub fn new(min: Array1<f64>, max: Array1<f64>) -> SpatialResult<Self> {
27        if min.len() != max.len() {
28            return Err(SpatialError::DimensionError(format!(
29                "Min and max must have the same dimensions. Got {} and {}",
30                min.len(),
31                max.len()
32            )));
33        }
34
35        // Validate that min ≤ max in all dimensions
36        for i in 0..min.len() {
37            if min[i] > max[i] {
38                return Err(SpatialError::ValueError(
39                    format!("Min must be less than or equal to max in all dimensions. Dimension {}: {} > {}", 
40                           i, min[i], max[i])
41                ));
42            }
43        }
44
45        Ok(Rectangle { min, max })
46    }
47
48    /// Create a rectangle from a point (zero-area rectangle)
49    ///
50    /// # Arguments
51    ///
52    /// * `point` - The point coordinates
53    ///
54    /// # Returns
55    ///
56    /// A rectangle containing only the given point
57    pub fn from_point(point: &ArrayView1<f64>) -> Self {
58        Rectangle {
59            min: point.to_owned(),
60            max: point.to_owned(),
61        }
62    }
63
64    /// Get the dimensionality of the rectangle
65    pub fn ndim(&self) -> usize {
66        self.min.len()
67    }
68
69    /// Calculate the area of the rectangle
70    pub fn area(&self) -> f64 {
71        let mut area = 1.0;
72        for i in 0..self.ndim() {
73            area *= self.max[i] - self.min[i];
74        }
75        area
76    }
77
78    /// Calculate the margin (perimeter) of the rectangle
79    pub fn margin(&self) -> f64 {
80        let mut margin = 0.0;
81        for i in 0..self.ndim() {
82            margin += self.max[i] - self.min[i];
83        }
84        2.0 * margin
85    }
86
87    /// Check if this rectangle contains a point
88    ///
89    /// # Arguments
90    ///
91    /// * `point` - The point to check
92    ///
93    /// # Returns
94    ///
95    /// `true` if the point is inside or on the boundary of the rectangle, `false` otherwise
96    pub fn contains_point(&self, point: &ArrayView1<f64>) -> SpatialResult<bool> {
97        if point.len() != self.ndim() {
98            return Err(SpatialError::DimensionError(format!(
99                "Point dimension {} does not match rectangle dimension {}",
100                point.len(),
101                self.ndim()
102            )));
103        }
104
105        for i in 0..self.ndim() {
106            if point[i] < self.min[i] || point[i] > self.max[i] {
107                return Ok(false);
108            }
109        }
110
111        Ok(true)
112    }
113
114    /// Check if this rectangle contains another rectangle
115    ///
116    /// # Arguments
117    ///
118    /// * `other` - The other rectangle to check
119    ///
120    /// # Returns
121    ///
122    /// `true` if the other rectangle is completely inside this rectangle, `false` otherwise
123    pub fn contains_rectangle(&self, other: &Rectangle) -> SpatialResult<bool> {
124        if other.ndim() != self.ndim() {
125            return Err(SpatialError::DimensionError(format!(
126                "Rectangle dimensions do not match: {} and {}",
127                self.ndim(),
128                other.ndim()
129            )));
130        }
131
132        for i in 0..self.ndim() {
133            if other.min[i] < self.min[i] || other.max[i] > self.max[i] {
134                return Ok(false);
135            }
136        }
137
138        Ok(true)
139    }
140
141    /// Check if this rectangle intersects another rectangle
142    ///
143    /// # Arguments
144    ///
145    /// * `other` - The other rectangle to check
146    ///
147    /// # Returns
148    ///
149    /// `true` if the rectangles intersect, `false` otherwise
150    pub fn intersects(&self, other: &Rectangle) -> SpatialResult<bool> {
151        if other.ndim() != self.ndim() {
152            return Err(SpatialError::DimensionError(format!(
153                "Rectangle dimensions do not match: {} and {}",
154                self.ndim(),
155                other.ndim()
156            )));
157        }
158
159        for i in 0..self.ndim() {
160            if self.max[i] < other.min[i] || self.min[i] > other.max[i] {
161                return Ok(false);
162            }
163        }
164
165        Ok(true)
166    }
167
168    /// Calculate the intersection of this rectangle with another rectangle
169    ///
170    /// # Arguments
171    ///
172    /// * `other` - The other rectangle
173    ///
174    /// # Returns
175    ///
176    /// A `SpatialResult` containing the intersection rectangle if it exists,
177    /// or an error if the rectangles do not intersect or have different dimensions
178    pub fn intersection(&self, other: &Rectangle) -> SpatialResult<Rectangle> {
179        if other.ndim() != self.ndim() {
180            return Err(SpatialError::DimensionError(format!(
181                "Rectangle dimensions do not match: {} and {}",
182                self.ndim(),
183                other.ndim()
184            )));
185        }
186
187        let intersects = self.intersects(other)?;
188        if !intersects {
189            return Err(SpatialError::ValueError(
190                "Rectangles do not intersect".into(),
191            ));
192        }
193
194        let mut min = Array1::zeros(self.ndim());
195        let mut max = Array1::zeros(self.ndim());
196
197        for i in 0..self.ndim() {
198            min[i] = f64::max(self.min[i], other.min[i]);
199            max[i] = f64::min(self.max[i], other.max[i]);
200        }
201
202        Rectangle::new(min, max)
203    }
204
205    /// Calculate the minimum bounding rectangle (MBR) containing this rectangle and another
206    ///
207    /// # Arguments
208    ///
209    /// * `other` - The other rectangle
210    ///
211    /// # Returns
212    ///
213    /// A `SpatialResult` containing the enlarged rectangle containing both rectangles,
214    /// or an error if the rectangles have different dimensions
215    pub fn enlarge(&self, other: &Rectangle) -> SpatialResult<Rectangle> {
216        if other.ndim() != self.ndim() {
217            return Err(SpatialError::DimensionError(format!(
218                "Rectangle dimensions do not match: {} and {}",
219                self.ndim(),
220                other.ndim()
221            )));
222        }
223
224        let mut min = Array1::zeros(self.ndim());
225        let mut max = Array1::zeros(self.ndim());
226
227        for i in 0..self.ndim() {
228            min[i] = f64::min(self.min[i], other.min[i]);
229            max[i] = f64::max(self.max[i], other.max[i]);
230        }
231
232        Rectangle::new(min, max)
233    }
234
235    /// Calculate the enlargement area needed to include another rectangle
236    ///
237    /// # Arguments
238    ///
239    /// * `other` - The other rectangle
240    ///
241    /// # Returns
242    ///
243    /// The difference between the area of the enlarged rectangle and this rectangle,
244    /// or an error if the rectangles have different dimensions
245    pub fn enlargement_area(&self, other: &Rectangle) -> SpatialResult<f64> {
246        let enlarged = self.enlarge(other)?;
247        Ok(enlarged.area() - self.area())
248    }
249
250    /// Calculate the minimum distance from this rectangle to a point
251    ///
252    /// # Arguments
253    ///
254    /// * `point` - The point
255    ///
256    /// # Returns
257    ///
258    /// The minimum distance from the rectangle to the point,
259    /// or an error if the point has a different dimension
260    pub fn min_distance_to_point(&self, point: &ArrayView1<f64>) -> SpatialResult<f64> {
261        if point.len() != self.ndim() {
262            return Err(SpatialError::DimensionError(format!(
263                "Point dimension {} does not match rectangle dimension {}",
264                point.len(),
265                self.ndim()
266            )));
267        }
268
269        let mut distance_sq = 0.0;
270
271        for i in 0..self.ndim() {
272            if point[i] < self.min[i] {
273                distance_sq += (point[i] - self.min[i]).powi(2);
274            } else if point[i] > self.max[i] {
275                distance_sq += (point[i] - self.max[i]).powi(2);
276            }
277        }
278
279        Ok(distance_sq.sqrt())
280    }
281
282    /// Calculate the minimum distance from this rectangle to another rectangle
283    ///
284    /// # Arguments
285    ///
286    /// * `other` - The other rectangle
287    ///
288    /// # Returns
289    ///
290    /// The minimum distance between the rectangles,
291    /// or an error if they have different dimensions
292    pub fn min_distance_to_rectangle(&self, other: &Rectangle) -> SpatialResult<f64> {
293        if other.ndim() != self.ndim() {
294            return Err(SpatialError::DimensionError(format!(
295                "Rectangle dimensions do not match: {} and {}",
296                self.ndim(),
297                other.ndim()
298            )));
299        }
300
301        // If the rectangles intersect, the distance is 0
302        if self.intersects(other)? {
303            return Ok(0.0);
304        }
305
306        let mut distance_sq = 0.0;
307
308        for i in 0..self.ndim() {
309            if self.max[i] < other.min[i] {
310                distance_sq += (other.min[i] - self.max[i]).powi(2);
311            } else if self.min[i] > other.max[i] {
312                distance_sq += (self.min[i] - other.max[i]).powi(2);
313            }
314        }
315
316        Ok(distance_sq.sqrt())
317    }
318}
319
320/// Entry type in R-tree node
321#[derive(Clone, Debug)]
322pub(crate) enum Entry<T>
323where
324    T: Clone,
325{
326    /// Leaf entry containing a data record
327    Leaf {
328        /// Bounding rectangle
329        mbr: Rectangle,
330        /// Data associated with this entry
331        data: T,
332        /// Index of the data in the original array
333        index: usize,
334    },
335
336    /// Non-leaf entry pointing to a child node
337    NonLeaf {
338        /// Bounding rectangle
339        mbr: Rectangle,
340        /// Child node
341        child: Box<Node<T>>,
342    },
343}
344
345impl<T: Clone> Entry<T> {
346    /// Get the minimum bounding rectangle of this entry
347    pub(crate) fn mbr(&self) -> &Rectangle {
348        match self {
349            Entry::Leaf { mbr, .. } => mbr,
350            Entry::NonLeaf { mbr, .. } => mbr,
351        }
352    }
353}
354
355/// Node in the R-tree
356#[derive(Clone, Debug)]
357pub(crate) struct Node<T>
358where
359    T: Clone,
360{
361    /// Entries stored in this node
362    pub entries: Vec<Entry<T>>,
363    /// Whether this is a leaf node
364    pub _isleaf: bool,
365    /// Node level in the tree (0 for leaf nodes, increasing towards the root)
366    pub level: usize,
367}
368
369impl<T: Clone> Default for Node<T> {
370    fn default() -> Self {
371        Self {
372            entries: Vec::new(),
373            _isleaf: true,
374            level: 0,
375        }
376    }
377}
378
379impl<T: Clone> Node<T> {
380    /// Create a new node
381    pub fn new(_isleaf: bool, level: usize) -> Self {
382        Node {
383            entries: Vec::new(),
384            _isleaf,
385            level,
386        }
387    }
388
389    /// Get the number of entries in the node
390    pub fn size(&self) -> usize {
391        self.entries.len()
392    }
393
394    /// Calculate the bounding rectangle for the node
395    pub fn mbr(&self) -> SpatialResult<Option<Rectangle>> {
396        if self.entries.is_empty() {
397            return Ok(None);
398        }
399
400        let mut result = self.entries[0].mbr().clone();
401
402        for i in 1..self.size() {
403            result = result.enlarge(self.entries[i].mbr())?;
404        }
405
406        Ok(Some(result))
407    }
408}
409
410/// Entry with distance for nearest neighbor search priority queue
411#[derive(Clone, Debug)]
412pub(crate) struct EntryWithDistance<T>
413where
414    T: Clone,
415{
416    /// Entry
417    pub entry: Entry<T>,
418    /// Distance to the query point
419    pub distance: f64,
420}
421
422impl<T: Clone> PartialEq for EntryWithDistance<T> {
423    fn eq(&self, other: &Self) -> bool {
424        self.distance == other.distance
425    }
426}
427
428impl<T: Clone> Eq for EntryWithDistance<T> {}
429
430impl<T: Clone> PartialOrd for EntryWithDistance<T> {
431    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
432        Some(self.cmp(other))
433    }
434}
435
436impl<T: Clone> Ord for EntryWithDistance<T> {
437    fn cmp(&self, other: &Self) -> Ordering {
438        // Reverse ordering - smaller distances come first
439        other
440            .distance
441            .partial_cmp(&self.distance)
442            .unwrap_or(Ordering::Equal)
443    }
444}
445
446/// R-tree for spatial indexing and queries
447///
448/// # Examples
449///
450/// ```
451/// use scirs2_spatial::rtree::RTree;
452/// use scirs2_core::ndarray::array;
453///
454/// // Create points
455/// let points = array![
456///     [0.0, 0.0],
457///     [1.0, 1.0],
458///     [2.0, 2.0],
459///     [3.0, 3.0],
460///     [4.0, 4.0],
461/// ];
462///
463/// // Build R-tree with points
464/// let mut rtree = RTree::new(2, 4, 8).unwrap();
465/// for (i, point) in points.rows().into_iter().enumerate() {
466///     rtree.insert(point.to_owned(), i).unwrap();
467/// }
468///
469/// // Search for points within a range
470/// let query_min = array![0.5, 0.5];
471/// let query_max = array![2.5, 2.5];
472///
473/// // The views need to be passed explicitly in the actual code
474/// let query_min_view = query_min.view();
475/// let query_max_view = query_max.view();
476/// let results = rtree.search_range(&query_min_view, &query_max_view).unwrap();
477///
478/// println!("Found {} points in range", results.len());
479///
480/// // Find the nearest neighbors to a point
481/// let query_point = array![2.5, 2.5];
482/// let query_point_view = query_point.view();
483/// let nearest = rtree.nearest(&query_point_view, 2).unwrap();
484///
485/// println!("Nearest points: {:?}", nearest);
486/// ```
487#[derive(Clone, Debug)]
488pub struct RTree<T>
489where
490    T: Clone,
491{
492    /// Root node of the tree
493    pub(crate) root: Node<T>,
494    /// Number of dimensions
495    ndim: usize,
496    /// Minimum number of entries in each node (except the root)
497    pub(crate) min_entries: usize,
498    /// Maximum number of entries in each node
499    pub(crate) maxentries: usize,
500    /// Number of data points in the tree
501    size: usize,
502    /// Height of the tree
503    height: usize,
504}
505
506impl<T: Clone> RTree<T> {
507    /// Create a new R-tree
508    ///
509    /// # Arguments
510    ///
511    /// * `ndim` - Number of dimensions
512    /// * `min_entries` - Minimum number of entries in each node (except the root)
513    /// * `maxentries` - Maximum number of entries in each node
514    ///
515    /// # Returns
516    ///
517    /// A `SpatialResult` containing the new R-tree, or an error if the parameters are invalid
518    pub fn new(ndim: usize, min_entries: usize, maxentries: usize) -> SpatialResult<Self> {
519        if ndim == 0 {
520            return Err(SpatialError::ValueError(
521                "Number of dimensions must be positive".into(),
522            ));
523        }
524
525        if min_entries < 1 || min_entries > maxentries / 2 {
526            return Err(SpatialError::ValueError(format!(
527                "min_entries must be between 1 and maxentries/2, got: {min_entries}"
528            )));
529        }
530
531        if maxentries < 2 {
532            return Err(SpatialError::ValueError(format!(
533                "maxentries must be at least 2, got: {maxentries}"
534            )));
535        }
536
537        Ok(RTree {
538            root: Node::new(true, 0),
539            ndim,
540            min_entries,
541            maxentries,
542            size: 0,
543            height: 1,
544        })
545    }
546
547    /// Get the number of dimensions in the R-tree
548    pub fn ndim(&self) -> usize {
549        self.ndim
550    }
551
552    /// Get the number of data points in the R-tree
553    pub fn size(&self) -> usize {
554        self.size
555    }
556
557    /// Get the height of the R-tree
558    pub fn height(&self) -> usize {
559        self.height
560    }
561
562    /// Check if the R-tree is empty
563    pub fn is_empty(&self) -> bool {
564        self.size == 0
565    }
566
567    /// Clear the R-tree, removing all data points
568    pub fn clear(&mut self) {
569        self.root = Node::new(true, 0);
570        self.size = 0;
571        self.height = 1;
572    }
573
574    /// Increment the size of the tree (internal use only)
575    pub(crate) fn increment_size(&mut self) {
576        self.size += 1;
577    }
578
579    /// Decrement the size of the tree (internal use only)
580    pub(crate) fn decrement_size(&mut self) {
581        if self.size > 0 {
582            self.size -= 1;
583        }
584    }
585
586    /// Increment the height of the tree (internal use only)
587    pub(crate) fn increment_height(&mut self) {
588        self.height += 1;
589    }
590
591    /// Decrement the height of the tree (internal use only)
592    #[allow(dead_code)]
593    pub(crate) fn decrement_height(&mut self) {
594        if self.height > 0 {
595            self.height -= 1;
596        }
597    }
598}
599
600#[cfg(test)]
601mod tests {
602    use super::*;
603    use approx::assert_relative_eq;
604    use scirs2_core::ndarray::array;
605
606    #[test]
607    fn test_rectangle_creation() {
608        let min = array![0.0, 0.0, 0.0];
609        let max = array![1.0, 2.0, 3.0];
610
611        let rect = Rectangle::new(min.clone(), max.clone()).unwrap();
612
613        assert_eq!(rect.ndim(), 3);
614        assert_eq!(rect.min, min);
615        assert_eq!(rect.max, max);
616
617        // Test invalid rectangle (min > max)
618        let min_invalid = array![0.0, 3.0, 0.0];
619        let max_invalid = array![1.0, 2.0, 3.0];
620
621        let result = Rectangle::new(min_invalid, max_invalid);
622        assert!(result.is_err());
623
624        // Test dimension mismatch
625        let min_dim_mismatch = array![0.0, 0.0];
626        let max_dim_mismatch = array![1.0, 2.0, 3.0];
627
628        let result = Rectangle::new(min_dim_mismatch, max_dim_mismatch);
629        assert!(result.is_err());
630    }
631
632    #[test]
633    fn test_rectangle_area_and_margin() {
634        let min = array![0.0, 0.0, 0.0];
635        let max = array![1.0, 2.0, 3.0];
636
637        let rect = Rectangle::new(min, max).unwrap();
638
639        assert_relative_eq!(rect.area(), 6.0);
640        assert_relative_eq!(rect.margin(), 12.0);
641    }
642
643    #[test]
644    fn test_rectangle_contains_point() {
645        let min = array![0.0, 0.0];
646        let max = array![1.0, 1.0];
647
648        let rect = Rectangle::new(min, max).unwrap();
649
650        // Inside
651        let point_inside = array![0.5, 0.5];
652        assert!(rect.contains_point(&point_inside.view()).unwrap());
653
654        // On boundary
655        let point_boundary = array![0.0, 0.5];
656        assert!(rect.contains_point(&point_boundary.view()).unwrap());
657
658        // Outside
659        let point_outside = array![2.0, 0.5];
660        assert!(!rect.contains_point(&point_outside.view()).unwrap());
661
662        // Dimension mismatch
663        let point_dim_mismatch = array![0.5, 0.5, 0.5];
664        assert!(rect.contains_point(&point_dim_mismatch.view()).is_err());
665    }
666
667    #[test]
668    fn test_rectangle_intersects() {
669        let rect1 = Rectangle::new(array![0.0, 0.0], array![2.0, 2.0]).unwrap();
670
671        // Overlap
672        let rect2 = Rectangle::new(array![1.0, 1.0], array![3.0, 3.0]).unwrap();
673        assert!(rect1.intersects(&rect2).unwrap());
674
675        // Touch
676        let rect3 = Rectangle::new(array![2.0, 0.0], array![3.0, 2.0]).unwrap();
677        assert!(rect1.intersects(&rect3).unwrap());
678
679        // No intersection
680        let rect4 = Rectangle::new(array![3.0, 3.0], array![4.0, 4.0]).unwrap();
681        assert!(!rect1.intersects(&rect4).unwrap());
682
683        // Dimension mismatch
684        let rect5 = Rectangle::new(array![0.0, 0.0, 0.0], array![1.0, 1.0, 1.0]).unwrap();
685        assert!(rect1.intersects(&rect5).is_err());
686    }
687
688    #[test]
689    fn test_rectangle_enlarge() {
690        let rect1 = Rectangle::new(array![1.0, 1.0], array![2.0, 2.0]).unwrap();
691        let rect2 = Rectangle::new(array![0.0, 0.0], array![3.0, 1.5]).unwrap();
692
693        let enlarged = rect1.enlarge(&rect2).unwrap();
694
695        assert_eq!(enlarged.min, array![0.0, 0.0]);
696        assert_eq!(enlarged.max, array![3.0, 2.0]);
697
698        // Calculate enlargement area
699        let enlargement = rect1.enlargement_area(&rect2).unwrap();
700        let original_area = rect1.area();
701        let enlarged_area = enlarged.area();
702
703        assert_relative_eq!(enlargement, enlarged_area - original_area);
704    }
705
706    #[test]
707    fn test_rectangle_distance() {
708        let rect = Rectangle::new(array![1.0, 1.0], array![3.0, 3.0]).unwrap();
709
710        // Inside: distance should be 0
711        let point_inside = array![2.0, 2.0];
712        assert_relative_eq!(
713            rect.min_distance_to_point(&point_inside.view()).unwrap(),
714            0.0
715        );
716
717        // Outside
718        let point_outside = array![0.0, 0.0];
719        assert_relative_eq!(
720            rect.min_distance_to_point(&point_outside.view()).unwrap(),
721            2.0_f64.sqrt()
722        );
723
724        // Other rectangle - no intersection
725        let rect2 = Rectangle::new(array![4.0, 4.0], array![5.0, 5.0]).unwrap();
726        assert_relative_eq!(
727            rect.min_distance_to_rectangle(&rect2).unwrap(),
728            2.0_f64.sqrt()
729        );
730
731        // Other rectangle - intersection (distance = 0)
732        let rect3 = Rectangle::new(array![2.0, 2.0], array![4.0, 4.0]).unwrap();
733        assert_relative_eq!(rect.min_distance_to_rectangle(&rect3).unwrap(), 0.0);
734    }
735
736    #[test]
737    fn test_rtree_creation() {
738        // Valid parameters
739        let rtree: RTree<usize> = RTree::new(2, 2, 5).unwrap();
740        assert_eq!(rtree.ndim(), 2);
741        assert_eq!(rtree.size(), 0);
742        assert_eq!(rtree.height(), 1);
743        assert!(rtree.is_empty());
744
745        // Invalid parameters
746        assert!(RTree::<usize>::new(0, 2, 5).is_err()); // ndim = 0
747        assert!(RTree::<usize>::new(2, 0, 5).is_err()); // min_entries = 0
748        assert!(RTree::<usize>::new(2, 3, 5).is_err()); // min_entries > maxentries/2
749        assert!(RTree::<usize>::new(2, 2, 1).is_err()); // maxentries < 2
750    }
751}
scirs2_spatial/rtree/node.rs

scirs2_spatial/rtree/
node.rs
