anime/
lib.rs

1pub mod interpolate;
2mod overlap;
3pub mod structs;
4
5use crate::{
6    overlap::*, overlap_range, solve_no_x_overlap, solve_no_y_overlap, structs::*, x_range,
7    y_range, TarLine,
8};
9use geo::{BoundingRect, EuclideanDistance, EuclideanLength};
10use rstar::primitives::{CachedEnvelope, GeomWithData};
11use std::{cell::OnceCell, collections::BTreeMap, error::Error, fmt::Display};
12
13/// Anime Error Type
14#[derive(Debug, Clone, Copy)]
15pub enum AnimeError {
16    IncorrectLength,
17    MatchesNotFound,
18}
19
20impl Display for AnimeError {
21    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
22        match self {
23            AnimeError::IncorrectLength => write!(f, "Variable to interpolate must have the same number of observations as the `target` lines"),
24            AnimeError::MatchesNotFound => write!(f, "`matches` needs to be instantiated with `self.find_matches()`"),
25        }
26    }
27}
28
29impl Error for AnimeError {}
30
31/// R* Tree for source geometries
32pub type SourceTree = rstar::RTree<GeomWithData<CachedEnvelope<geo_types::Line>, (usize, f64)>>;
33
34/// R* Tree for target geometries
35pub type TargetTree = rstar::RTree<GeomWithData<CachedEnvelope<TarLine>, (usize, f64)>>;
36
37/// Represents a partial source <-> target match
38#[derive(Debug, Clone)]
39pub struct MatchCandidate {
40    /// The index of the target geometry as an i32
41    pub index: i32,
42    /// The amount of shared length between two geometries
43    pub shared_len: f64,
44}
45
46/// Stores match length
47///
48/// The BTreeMap key is the index of the source geometry
49/// whereas the entry contains
50pub type MatchesMap = BTreeMap<i32, Vec<MatchCandidate>>;
51
52/// Approximate Network Matching, Integration, and Enrichment
53///
54/// This struct contains all of the information needed to perform
55/// and store the results of the ANIME algorithm.
56///
57/// The `source_tree` and `target_tree` are used to perform the
58/// partial matching based on the `distance_tolerance` and
59/// `angle_tolerance`. The results of the matching
60/// are stored in the `BTreeMap`.
61///
62/// The lengths, represented as `Vec<f64>` are required for the
63/// integration of attributes.
64#[derive(Clone, Debug)]
65pub struct Anime {
66    pub distance_tolerance: f64,
67    pub angle_tolerance: f64,
68    pub source_tree: SourceTree,
69    pub source_lens: Vec<f64>,
70    pub target_tree: TargetTree,
71    pub target_lens: Vec<f64>,
72    pub matches: OnceCell<MatchesMap>,
73}
74
75impl Anime {
76    /// Load source and target `LineString` geometries
77    ///
78    /// This creates two R* Trees using cached envelopes for each component
79    /// line in a LineString. In addition to the envelope, the slope and
80    /// index of the LineString is stored.
81    pub fn load_geometries(
82        source: impl Iterator<Item = geo_types::LineString>,
83        target: impl Iterator<Item = geo_types::LineString>,
84        distance_tolerance: f64,
85        angle_tolerance: f64,
86    ) -> Self {
87        let mut source_lens = Vec::new();
88        let mut target_lens = Vec::new();
89        let source_tree = create_source_rtree(source, &mut source_lens);
90        let target_tree = create_target_rtree(target, &mut target_lens, distance_tolerance);
91        Self {
92            distance_tolerance,
93            angle_tolerance,
94            source_tree,
95            source_lens,
96            target_tree,
97            target_lens,
98            matches: OnceCell::new(),
99        }
100    }
101
102    /// Find candidate matches between source and target
103    ///
104    /// The matches can only be found once for each source and target pair.
105    pub fn find_matches(&mut self) -> Result<&mut Anime, MatchesMap> {
106        let mut matches: MatchesMap = BTreeMap::new();
107        let candidates = self
108            .source_tree
109            .intersection_candidates_with_other_tree(&self.target_tree);
110
111        candidates.for_each(|(cx, cy)| {
112            let xbb = cx.geom().bounding_rect();
113            let ybb = cy.geom().0.bounding_rect();
114
115            // extract cached slopes and index positions
116            let (i, x_slope) = cx.data;
117            let (j, y_slope) = cy.data;
118
119            // convert calculated slopes to degrees
120            let x_deg = x_slope.atan().to_degrees();
121            let y_deg = y_slope.atan().to_degrees();
122
123            // compare slopes:
124            let is_tolerant = (x_deg - y_deg).abs() < self.angle_tolerance;
125
126            // if the slopes are within tolerance then we check for overlap
127            if is_tolerant {
128                let xx_range = x_range(&xbb);
129                let xy_range = x_range(&ybb);
130                let x_overlap = overlap_range(xx_range, xy_range);
131                let y_overlap = overlap_range(y_range(&xbb), y_range(&ybb));
132
133                // if theres overlap then we do a distance based check
134                // following, check that they're within distance tolerance,
135                // if so, calculate the shared length
136                if x_overlap.is_some() || y_overlap.is_some() {
137                    // calculate the distance from the line segment
138                    // if its within our threshold we include it;
139                    let d = cy.geom().distance(cx.geom());
140
141                    // if distance is less than or equal to tolerance, add the key
142                    if d <= self.distance_tolerance {
143                        let shared_len = if x_slope.atan().to_degrees() <= 45.0 {
144                            if x_overlap.is_some() {
145                                let (p1, p2) =
146                                    solve_no_y_overlap(x_overlap.unwrap(), cx.geom(), &x_slope);
147
148                                p1.euclidean_distance(&p2)
149                            } else {
150                                0.0
151                            }
152                        } else if y_overlap.is_some() {
153                            let (p1, p2) =
154                                solve_no_x_overlap(y_overlap.unwrap(), cx.geom(), &x_slope);
155                            p1.euclidean_distance(&p2)
156                        } else {
157                            0.0
158                        };
159                        // add 1 for R indexing
160                        // ensures that no duplicates are inserted. Creates a new empty vector is needed
161                        let entry = matches.entry(j as i32).or_default();
162
163                        if let Some(tuple) = entry.iter_mut().find(|x| x.index == i as i32) {
164                            tuple.shared_len += shared_len;
165                        } else {
166                            entry.extend(std::iter::once(MatchCandidate {
167                                index: i as i32,
168                                shared_len,
169                            }));
170                        }
171                    }
172                }
173            }
174        });
175        self.matches.set(matches)?;
176        Ok(self)
177    }
178}
179
180fn create_source_rtree(
181    x: impl Iterator<Item = geo_types::LineString>,
182    source_lens: &mut Vec<f64>,
183) -> SourceTree {
184    let to_insert = x
185        .enumerate()
186        .flat_map(|(i, xi)| {
187            let xi_len = xi.euclidean_length();
188            source_lens.push(xi_len);
189            let components = xi
190                .lines()
191                .map(|li| {
192                    let slope = li.slope();
193                    let env = CachedEnvelope::new(li);
194                    GeomWithData::new(env, (i, slope))
195                })
196                .collect::<Vec<GeomWithData<_, _>>>();
197            components
198        })
199        .collect::<Vec<_>>();
200
201    rstar::RTree::bulk_load(to_insert)
202}
203
204fn create_target_rtree(
205    y: impl Iterator<Item = geo_types::LineString>,
206    target_lens: &mut Vec<f64>,
207    dist: f64,
208) -> TargetTree {
209    let to_insert = y
210        .enumerate()
211        .flat_map(|(i, yi)| {
212            let yi_len = yi.euclidean_length();
213            target_lens.push(yi_len);
214            let components = yi
215                .lines()
216                .map(|li| {
217                    let tl = TarLine(li, dist);
218                    let slope = li.slope();
219                    let env = CachedEnvelope::new(tl);
220                    GeomWithData::new(env, (i, slope))
221                })
222                .collect::<Vec<GeomWithData<_, _>>>();
223            components
224        })
225        .collect::<Vec<_>>();
226
227    rstar::RTree::bulk_load(to_insert)
228}
anime/lib.rs

anime/
lib.rs
