oxigdal_algorithms/vector/
spatial_join.rs

1//! Spatial join operations with spatial indexing
2//!
3//! Efficient spatial joins using R-tree and other spatial indices.
4
5use crate::error::Result;
6use oxigdal_core::vector::Point;
7use rstar::{AABB, PointDistance, RTree, RTreeObject};
8
9/// Spatial join predicate
10#[derive(Debug, Clone, Copy, PartialEq, Eq)]
11pub enum SpatialJoinPredicate {
12    /// Features intersect
13    Intersects,
14    /// First feature contains second
15    Contains,
16    /// First feature is within second
17    Within,
18    /// Features touch (share boundary)
19    Touches,
20    /// Features are within distance
21    WithinDistance,
22}
23
24/// Options for spatial join
25#[derive(Debug, Clone)]
26pub struct SpatialJoinOptions {
27    /// Spatial predicate
28    pub predicate: SpatialJoinPredicate,
29    /// Distance threshold (for WithinDistance predicate)
30    pub distance: f64,
31    /// Whether to build spatial index
32    pub use_index: bool,
33}
34
35impl Default for SpatialJoinOptions {
36    fn default() -> Self {
37        Self {
38            predicate: SpatialJoinPredicate::Intersects,
39            distance: 0.0,
40            use_index: true,
41        }
42    }
43}
44
45/// Result of spatial join
46#[derive(Debug, Clone)]
47pub struct SpatialJoinResult {
48    /// Pairs of matching indices (left_idx, right_idx)
49    pub matches: Vec<(usize, usize)>,
50    /// Number of matches
51    pub num_matches: usize,
52}
53
54/// Indexed point for R-tree
55#[derive(Debug, Clone)]
56struct IndexedPoint {
57    point: Point,
58    index: usize,
59}
60
61impl RTreeObject for IndexedPoint {
62    type Envelope = AABB<[f64; 2]>;
63
64    fn envelope(&self) -> Self::Envelope {
65        AABB::from_point([self.point.coord.x, self.point.coord.y])
66    }
67}
68
69impl PointDistance for IndexedPoint {
70    fn distance_2(&self, point: &[f64; 2]) -> f64 {
71        let dx = self.point.coord.x - point[0];
72        let dy = self.point.coord.y - point[1];
73        dx * dx + dy * dy
74    }
75}
76
77/// Perform spatial join between two point sets
78///
79/// # Arguments
80///
81/// * `left_points` - First set of points
82/// * `right_points` - Second set of points
83/// * `options` - Spatial join options
84///
85/// # Returns
86///
87/// Join result with matching pairs
88///
89/// # Examples
90///
91/// ```
92/// use oxigdal_algorithms::vector::spatial_join::{spatial_join_points, SpatialJoinOptions, SpatialJoinPredicate};
93/// use oxigdal_algorithms::Point;
94/// # use oxigdal_algorithms::error::Result;
95///
96/// # fn main() -> Result<()> {
97/// let left = vec![
98///     Point::new(0.0, 0.0),
99///     Point::new(1.0, 1.0),
100/// ];
101///
102/// let right = vec![
103///     Point::new(0.1, 0.1),
104///     Point::new(10.0, 10.0),
105/// ];
106///
107/// let options = SpatialJoinOptions {
108///     predicate: SpatialJoinPredicate::WithinDistance,
109///     distance: 0.5,
110///     use_index: true,
111/// };
112///
113/// let result = spatial_join_points(&left, &right, &options)?;
114/// assert!(result.num_matches >= 1);
115/// # Ok(())
116/// # }
117/// ```
118pub fn spatial_join_points(
119    left_points: &[Point],
120    right_points: &[Point],
121    options: &SpatialJoinOptions,
122) -> Result<SpatialJoinResult> {
123    if left_points.is_empty() || right_points.is_empty() {
124        return Ok(SpatialJoinResult {
125            matches: Vec::new(),
126            num_matches: 0,
127        });
128    }
129
130    let matches = if options.use_index {
131        // Build R-tree for right points
132        let indexed_points: Vec<IndexedPoint> = right_points
133            .iter()
134            .enumerate()
135            .map(|(idx, point)| IndexedPoint {
136                point: point.clone(),
137                index: idx,
138            })
139            .collect();
140
141        let rtree = RTree::bulk_load(indexed_points);
142
143        // Query R-tree for each left point
144        let mut all_matches = Vec::new();
145
146        for (left_idx, left_point) in left_points.iter().enumerate() {
147            let nearby = match options.predicate {
148                SpatialJoinPredicate::WithinDistance => {
149                    // Query points within distance
150                    let envelope = AABB::from_corners(
151                        [
152                            left_point.coord.x - options.distance,
153                            left_point.coord.y - options.distance,
154                        ],
155                        [
156                            left_point.coord.x + options.distance,
157                            left_point.coord.y + options.distance,
158                        ],
159                    );
160
161                    rtree
162                        .locate_in_envelope(envelope)
163                        .filter(|indexed| {
164                            point_distance(left_point, &indexed.point) <= options.distance
165                        })
166                        .map(|indexed| indexed.index)
167                        .collect::<Vec<_>>()
168                }
169                SpatialJoinPredicate::Intersects => {
170                    // For points, intersects means exactly coincident
171                    let mut matches = Vec::new();
172                    for indexed in rtree.locate_at_point([left_point.coord.x, left_point.coord.y]) {
173                        matches.push(indexed.index);
174                    }
175                    matches
176                }
177                _ => {
178                    // Other predicates not applicable to points
179                    Vec::new()
180                }
181            };
182
183            for right_idx in nearby {
184                all_matches.push((left_idx, right_idx));
185            }
186        }
187
188        all_matches
189    } else {
190        // Brute force comparison
191        let mut all_matches = Vec::new();
192
193        for (left_idx, left_point) in left_points.iter().enumerate() {
194            for (right_idx, right_point) in right_points.iter().enumerate() {
195                if matches_predicate(left_point, right_point, options) {
196                    all_matches.push((left_idx, right_idx));
197                }
198            }
199        }
200
201        all_matches
202    };
203
204    Ok(SpatialJoinResult {
205        num_matches: matches.len(),
206        matches,
207    })
208}
209
210/// Check if two points match the join predicate
211fn matches_predicate(left: &Point, right: &Point, options: &SpatialJoinOptions) -> bool {
212    match options.predicate {
213        SpatialJoinPredicate::Intersects => {
214            (left.coord.x - right.coord.x).abs() < 1e-10
215                && (left.coord.y - right.coord.y).abs() < 1e-10
216        }
217        SpatialJoinPredicate::WithinDistance => point_distance(left, right) <= options.distance,
218        _ => false,
219    }
220}
221
222/// Calculate Euclidean distance between points
223fn point_distance(p1: &Point, p2: &Point) -> f64 {
224    let dx = p1.coord.x - p2.coord.x;
225    let dy = p1.coord.y - p2.coord.y;
226    (dx * dx + dy * dy).sqrt()
227}
228
229/// Nearest neighbor search
230pub fn nearest_neighbor(query: &Point, points: &[Point]) -> Option<(usize, f64)> {
231    let indexed_points: Vec<IndexedPoint> = points
232        .iter()
233        .enumerate()
234        .map(|(idx, point)| IndexedPoint {
235            point: point.clone(),
236            index: idx,
237        })
238        .collect();
239
240    if indexed_points.is_empty() {
241        return None;
242    }
243
244    let rtree = RTree::bulk_load(indexed_points);
245    let nearest = rtree.nearest_neighbor([query.coord.x, query.coord.y])?;
246
247    let distance = point_distance(query, &nearest.point);
248
249    Some((nearest.index, distance))
250}
251
252/// K-nearest neighbors search
253pub fn k_nearest_neighbors(query: &Point, points: &[Point], k: usize) -> Vec<(usize, f64)> {
254    let indexed_points: Vec<IndexedPoint> = points
255        .iter()
256        .enumerate()
257        .map(|(idx, point)| IndexedPoint {
258            point: point.clone(),
259            index: idx,
260        })
261        .collect();
262
263    if indexed_points.is_empty() {
264        return Vec::new();
265    }
266
267    let rtree = RTree::bulk_load(indexed_points);
268
269    rtree
270        .nearest_neighbor_iter([query.coord.x, query.coord.y])
271        .take(k)
272        .map(|indexed| {
273            let dist = point_distance(query, &indexed.point);
274            (indexed.index, dist)
275        })
276        .collect()
277}
278
279/// Range query (all points within distance)
280pub fn range_query(query: &Point, points: &[Point], distance: f64) -> Vec<usize> {
281    let options = SpatialJoinOptions {
282        predicate: SpatialJoinPredicate::WithinDistance,
283        distance,
284        use_index: true,
285    };
286
287    let result = spatial_join_points(std::slice::from_ref(query), points, &options);
288
289    result
290        .map(|r| {
291            r.matches
292                .into_iter()
293                .map(|(_, right_idx)| right_idx)
294                .collect()
295        })
296        .unwrap_or_else(|_| Vec::new())
297}
298
299#[cfg(test)]
300mod tests {
301    use super::*;
302
303    #[test]
304    fn test_spatial_join_within_distance() {
305        let left = vec![Point::new(0.0, 0.0), Point::new(10.0, 10.0)];
306
307        let right = vec![Point::new(0.1, 0.1), Point::new(5.0, 5.0)];
308
309        let options = SpatialJoinOptions {
310            predicate: SpatialJoinPredicate::WithinDistance,
311            distance: 0.5,
312            use_index: true,
313        };
314
315        let result = spatial_join_points(&left, &right, &options);
316        assert!(result.is_ok());
317
318        let join_result = result.expect("Join failed");
319        assert!(join_result.num_matches >= 1);
320    }
321
322    #[test]
323    fn test_nearest_neighbor() {
324        let points = vec![
325            Point::new(0.0, 0.0),
326            Point::new(5.0, 5.0),
327            Point::new(10.0, 10.0),
328        ];
329
330        let query = Point::new(0.1, 0.1);
331        let result = nearest_neighbor(&query, &points);
332
333        assert!(result.is_some());
334
335        let (idx, dist) = result.expect("Nearest neighbor failed");
336        assert_eq!(idx, 0);
337        assert!(dist < 0.2);
338    }
339
340    #[test]
341    fn test_k_nearest_neighbors() {
342        let points = vec![
343            Point::new(0.0, 0.0),
344            Point::new(1.0, 1.0),
345            Point::new(2.0, 2.0),
346            Point::new(10.0, 10.0),
347        ];
348
349        let query = Point::new(0.0, 0.0);
350        let result = k_nearest_neighbors(&query, &points, 2);
351
352        assert_eq!(result.len(), 2);
353        assert_eq!(result[0].0, 0); // First is the query point itself
354    }
355
356    #[test]
357    fn test_range_query() {
358        let points = vec![
359            Point::new(0.0, 0.0),
360            Point::new(0.5, 0.5),
361            Point::new(10.0, 10.0),
362        ];
363
364        let query = Point::new(0.0, 0.0);
365        let result = range_query(&query, &points, 1.0);
366
367        assert!(result.len() >= 2); // Should find points at 0.0 and 0.5
368    }
369
370    #[test]
371    fn test_point_distance() {
372        let p1 = Point::new(0.0, 0.0);
373        let p2 = Point::new(3.0, 4.0);
374
375        let dist = point_distance(&p1, &p2);
376        assert!((dist - 5.0).abs() < 1e-6);
377    }
378}
oxigdal_algorithms/vector/spatial_join.rs

oxigdal_algorithms/vector/
spatial_join.rs
