Add a function to use the index to find the nearest lrms of a coord

Signed-off-by: Tristram Gräbener <tristram+git@tristramg.eu>

Author: Tristramg

python/src/lib.rs

@@ -288,9 +288,27 @@ impl Anchor {
     }
 }
 
+impl From<&liblrs::lrm_scale::Anchor> for Anchor {
+    fn from(value: &liblrs::lrm_scale::Anchor) -> Self {
+        let name = match value {
+            liblrs::lrm_scale::Anchor::Named(value) => value.name.clone(),
+            liblrs::lrm_scale::Anchor::Unnamed(_) => "-".to_owned(),
+        };
+        Self {
+            name,
+            position: value.point().map(|p| p.into()),
+            curve_position: value.curve_position(),
+            scale_position: value.scale_position(),
+        }
+    }
+}
+
 #[pyclass]
 /// The result of a projection onto an [`LrmScale`].
 pub struct LrmProjection {
+    #[pyo3(get, set)]
+    /// Handle of the [`Lrm`] where this measurement applies to
+    lrm_handle: usize,
     /// Contains `measure` ([`LrmScaleMeasure`]) and `lrm` ([`LrmHandle`]).
     #[pyo3(get, set)]
     pub measure: LrmScaleMeasure,
@@ -299,18 +317,15 @@ pub struct LrmProjection {
     pub orthogonal_offset: f64,
 }
 
-impl From<&liblrs::lrm_scale::Anchor> for Anchor {
-    fn from(value: &liblrs::lrm_scale::Anchor) -> Self {
-        let name = match value {
-            liblrs::lrm_scale::Anchor::Named(value) => value.name.clone(),
-            liblrs::lrm_scale::Anchor::Unnamed(_) => "-".to_owned(),
-        };
-
+impl From<&liblrs::lrs::LrmProjection> for LrmProjection {
+    fn from(value: &liblrs::lrs::LrmProjection) -> Self {
         Self {
-            name,
-            position: value.point().map(|p| p.into()),
-            curve_position: value.curve_position(),
-            scale_position: value.scale_position(),
+            lrm_handle: value.measure.lrm.0,
+            measure: LrmScaleMeasure {
+                anchor_name: value.measure.measure.anchor_name.to_owned(),
+                scale_offset: value.measure.measure.scale_offset,
+            },
+            orthogonal_offset: value.orthogonal_offset,
         }
     }
 }
@@ -395,13 +410,19 @@ impl Lrs {
             .lrs
             .lookup(point.into(), LrmHandle(lrm_handle))
             .iter()
-            .map(|p| LrmProjection {
-                measure: LrmScaleMeasure {
-                    anchor_name: p.measure.measure.anchor_name.to_owned(),
-                    scale_offset: p.measure.measure.scale_offset,
-                },
-                orthogonal_offset: p.orthogonal_offset,
-            })
+            .map(LrmProjection::from)
+            .collect()
+    }
+
+    /// Projects a [`Point`] on all applicable [`Traversal`]s nearby.
+    /// The [`Point`] must be in the bounding box of the [`Curve`] of the [`Traversal`].
+    /// The result is sorted by `orthogonal_offset`: the nearest [`Lrm`] to the [`Point`] is the first item.
+    fn lookup_lrms(&self, point: Point) -> Vec<LrmProjection> {
+        self.lrs
+            .lrs
+            .lookup_lrms(point.into())
+            .iter()
+            .map(LrmProjection::from)
             .collect()
     }
 

src/builder.rs

@@ -65,6 +65,7 @@ struct TempTraversal {
     id: String,
     curve: SphericalLineStringCurve,
     segments: Vec<SegmentOfTraversal>,
+    lrms: Vec<usize>,
 }
 
 impl TempTraversal {
@@ -248,6 +249,7 @@ impl<'fbb> Builder<'fbb> {
             id: traversal_id.to_owned(),
             curve: SphericalLineStringCurve::new(geo::LineString::new(coords), 100.),
             segments: segments.to_vec(),
+            lrms: vec![],
         });
         self.nodes_of_traversal.push(nodes_of_traversal);
 
@@ -279,6 +281,9 @@ impl<'fbb> Builder<'fbb> {
             projected_anchors: Some(self.project_anchors(&anchors, traversal_index)),
             ..Default::default()
         };
+        self.temp_traversal[traversal_index]
+            .lrms
+            .push(self.lrms.len());
         self.lrms
             .push(LinearReferencingMethod::create(&mut self.fbb, &args));
     }
@@ -544,14 +549,14 @@ impl<'fbb> Builder<'fbb> {
 }
 
 #[cfg(test)]
-mod tests {
+pub(crate) mod tests {
     use crate::lrm_scale::LrmScaleMeasure;
 
     use super::*;
     use approx::assert_relative_eq;
     use geo::{coord, point};
 
-    fn build_traversal(builder: &mut Builder) -> usize {
+    pub fn build_traversal(builder: &mut Builder) -> usize {
         let s1 = builder.add_segment("s1", &[coord! {x: 0., y: 0.}, coord! {x:1., y: 0.}], 0, 1);
         let s2 = builder.add_segment("s2", &[coord! {x: 1., y: 0.}, coord! {x:2., y: 0.}], 1, 2);
         let sot1 = super::SegmentOfTraversal {

src/curves.rs

@@ -16,6 +16,11 @@ use thiserror::Error;
 /// A [`Curve`] can be part of a larger [`Curve`] (e.g. for optimisation purposes and to have better bounding boxes).
 /// The [`Curve`] can be implemented.
 pub trait Curve {
+    /// Is the curve represented in spherical coordinates (like longitude;latitude)
+    ///
+    /// This has an influence on how to compute distances and bearings
+    const IS_SPHERICAL: bool;
+
     /// Builds a new [`Curve`] from a [`LineString`].
     /// `max_extent` is the maximum distance that is considered to be “on the curve”.
     /// `max_extent` plays a role in the bounding box.
@@ -101,6 +106,8 @@ pub struct PlanarLineStringCurve {
 }
 
 impl Curve for PlanarLineStringCurve {
+    const IS_SPHERICAL: bool = false;
+
     fn new(geom: LineString, max_extent: f64) -> Self {
         let length = geom.length(&Euclidean);
         Self {
@@ -357,6 +364,8 @@ impl SphericalLineStringCurve {
 }
 
 impl Curve for SphericalLineStringCurve {
+    const IS_SPHERICAL: bool = true;
+
     fn new(geom: LineString, max_extent: f64) -> Self {
         let length = geom.length(&Geodesic);
         Self {

src/lrs.rs

@@ -10,14 +10,15 @@ use std::cmp::Ordering;
 
 use flatbuffers::{ForwardsUOffset, Vector};
 use geo::orient::Direction;
+use geo_index::rtree::{RTreeIndex, RTreeRef};
 use thiserror::Error;
 
 use crate::curves::{Curve, CurveError};
 use crate::lrm_scale::{
     Anchor, CurvePosition, LrmScale, LrmScaleError, LrmScaleMeasure, ScalePosition,
 };
 use crate::lrs_generated;
-use geo::{LineString, Point, coord, point};
+use geo::{Contains, LineString, Point, coord, point};
 
 /// Used as handle to identify a [`LrmScale`] within a specific [`Lrs`].
 #[derive(Copy, Clone, Debug, Eq, Hash, PartialEq)]
@@ -66,7 +67,7 @@ pub struct Lrs<CurveImpl: Curve> {
     /// All the [`Segment`] of this Lrs
     pub segments: Vec<Segment>,
     /// An RTree spatial index of the LRM extents
-    rtree_data: Option<Vec<u8>>,
+    pub rtree_data: Option<Vec<u8>>,
 }
 
 /// A Node is a topological element of the [`Lrs`] that represents a intersection (or an extremity) of an [`Lrm`]
@@ -400,6 +401,10 @@ pub trait LrsBase {
     /// Projects a [`Point`] on all [`Lrm`] where the [`Point`] is in the bounding box.
     /// The result is sorted by `orthogonal_offset`: the nearest [`Lrm`] to the [`Point`] is the first item.
     fn lookup_lrms(&self, point: Point) -> Vec<LrmProjection>;
+    /// Returns all the traversals whose bounding box include the given point
+    ///
+    /// The function will use the spatial index if it is defined
+    fn traversals_containing(&self, point: Point) -> Vec<TraversalHandle>;
 
     /// Given a [`TraversalPosition`], returns it geographical position ([`Point`]).
     fn locate_traversal(&self, position: TraversalPosition) -> Result<Point, LrsError>;
@@ -482,12 +487,37 @@ impl<CurveImpl: Curve> LrsBase for Lrs<CurveImpl> {
         })
     }
 
+    fn traversals_containing(&self, point: Point) -> Vec<TraversalHandle> {
+        let rtree = self
+            .rtree_data
+            .as_ref()
+            .and_then(|buf| RTreeRef::try_new(buf).ok());
+
+        rtree
+            .map(|tree| {
+                // The rectangles inserted in the rtree are already expanded by the margin
+                // That is why the min and max values are the same.
+                tree.search(point.x(), point.y(), point.x(), point.y())
+                    .iter()
+                    .map(|idx| TraversalHandle(*idx as usize))
+                    .collect()
+            })
+            .unwrap_or(
+                self.traversals
+                    .iter()
+                    .enumerate()
+                    .filter(|(_idx, traversal)| traversal.curve.bbox().contains(&point))
+                    .map(|(idx, _traversal)| TraversalHandle(idx))
+                    .collect(),
+            )
+    }
+
     fn lookup_lrms(&self, point: Point) -> Vec<LrmProjection> {
         let mut result: Vec<_> = self
-            .lrms
+            .traversals_containing(point)
             .iter()
-            .enumerate()
-            .flat_map(|(lrm_idx, _lrm)| self.lookup(point, LrmHandle(lrm_idx)))
+            .flat_map(|traversal_handle| &self.traversals[traversal_handle.0].lrms)
+            .flat_map(|&lrm_handle| self.lookup(point, lrm_handle))
             .collect();
         result.sort_by(|a, b| {
             a.orthogonal_offset

src/lrs_ext.rs

@@ -4,11 +4,8 @@
 use geo::{Coord, Point};
 
 use crate::curves::{Curve, SphericalLineStringCurve};
-use crate::lrm_scale::Anchor;
-use crate::lrm_scale::LrmScaleMeasure;
-use crate::lrs::Properties;
-use crate::lrs::{self, TraversalPosition};
-use crate::lrs::{LrsBase, LrsError};
+use crate::lrm_scale::{Anchor, LrmScaleMeasure};
+use crate::lrs::{self, LrmProjection, LrsBase, LrsError, Properties, TraversalPosition};
 
 type Lrs = lrs::Lrs<SphericalLineStringCurve>;
 
@@ -101,4 +98,72 @@ impl ExtLrs {
     pub fn anchor_properties(&self, lrm_index: usize, anchor_index: usize) -> &Properties {
         self.lrs.lrms[lrm_index].scale.anchors[anchor_index].properties()
     }
+
+    /// Projects a [`Point`] on all [`Lrm`] where the [`Point`] is in the bounding box.
+    /// The result is sorted by `orthogonal_offset`: the nearest [`Lrm`] to the [`Point`] is the first item.
+    pub fn lookup_lrms(&self, point: Point) -> Vec<LrmProjection> {
+        self.lrs.lookup_lrms(point)
+    }
+}
+
+#[cfg(test)]
+pub(crate) mod tests {
+    use geo::{Coord, coord, point};
+
+    use crate::builder::{AnchorOnLrm, Builder, SegmentOfTraversal};
+    use crate::{lrs, properties};
+
+    fn build_lrm(builder: &mut Builder, name: &str, coords: &[Coord]) {
+        let segment_index = builder.add_segment("name", coords, 0, 1);
+        let sot = SegmentOfTraversal {
+            segment_index,
+            reversed: false,
+        };
+        let traversal_index = builder.add_traversal(name, &[sot]);
+        let start_anchor = AnchorOnLrm {
+            anchor_index: builder.add_anchor("start", Some("start"), coords[0], properties!()),
+            distance_along_lrm: 0.0,
+        };
+        let end_anchor = AnchorOnLrm {
+            anchor_index: builder.add_anchor(
+                "end",
+                Some("end"),
+                *coords.last().unwrap(),
+                properties!(),
+            ),
+            distance_along_lrm: 1.0,
+        };
+        builder.add_lrm(
+            name,
+            traversal_index,
+            &[start_anchor, end_anchor],
+            properties!(),
+        );
+    }
+
+    #[test]
+    fn test() {
+        let mut b = Builder::new();
+        build_lrm(&mut b, "lrm1", &[coord! {x:0., y:0.}, coord! {x:2., y:0.}]);
+        build_lrm(&mut b, "lrm2", &[coord! {x:0., y:1.}, coord! {x:2., y:1.}]);
+        build_lrm(&mut b, "lrm3", &[coord! {x:1., y:1.}, coord! {x:1., y:-1.}]);
+        let lrs = b.build_lrs(properties!()).unwrap();
+
+        let nearest0 = lrs.lookup_lrms(point! {x: 1.0001, y:0.0});
+        assert_eq!(nearest0.len(), 2);
+        assert_eq!(nearest0[0].measure.lrm, lrs::LrmHandle(0));
+        assert_eq!(nearest0[1].measure.lrm, lrs::LrmHandle(2));
+
+        let nearest1 = lrs.lookup_lrms(point! {x: 0.5, y:1.0});
+        assert_eq!(nearest1.len(), 1);
+        assert_eq!(nearest1[0].measure.lrm, lrs::LrmHandle(1));
+
+        let nearest2 = lrs.lookup_lrms(point! {x:1., y:-0.0001});
+        assert_eq!(nearest2.len(), 2);
+        assert_eq!(nearest2[0].measure.lrm, lrs::LrmHandle(2));
+        assert_eq!(nearest2[1].measure.lrm, lrs::LrmHandle(0));
+
+        let nearest3 = lrs.lookup_lrms(point! {x:2.35, y:48.98});
+        assert!(nearest3.is_empty());
+    }
 }