gistools/geometry/tools/polys/
path_builder.rs

1use crate::geometry::{RingIntersection, RingIntersectionLookup};
2use alloc::{
3    collections::{BTreeMap, BTreeSet},
4    rc::Rc,
5    vec,
6    vec::Vec,
7};
8use core::{
9    cell::RefCell,
10    cmp::Ordering,
11    f64::consts::{PI, TAU},
12    mem::take,
13};
14use libm::atan2;
15use s2json::{BBox, FullXY, Point};
16
17/// Reconstructing a poly line that interacts with intersections
18#[derive(Debug)]
19pub struct PolyPath<P: FullXY> {
20    /// helps down the road to spot duplicate pulls of this Path
21    pub id: usize,
22    /// True outer ring
23    pub outer: Option<Vec<P>>,
24    /// Outers that have already been consumed
25    pub old_outers: Vec<BBox>,
26    /// Holes
27    pub holes: Vec<Vec<P>>,
28    /// indexes of the polygons in the multipolygon. So we can quickly consume holes.
29    pub polys_consumed: BTreeSet<usize>,
30    /// Bounding box
31    pub bbox: BBox,
32}
33/// A reference to a PolyPath wrapped in an RC & RefCell
34pub type PolyPathRef<P> = Rc<RefCell<PolyPath<P>>>;
35impl<P: FullXY> PolyPath<P> {
36    /// Create a new PolyPath as a PolyPathRef
37    pub fn new_ref(
38        ring: Vec<P>,
39        polys_consumed: BTreeSet<usize>,
40        is_outer: bool,
41        bbox: Option<BBox>,
42    ) -> PolyPathRef<P> {
43        Rc::new(RefCell::new(PolyPath::new(ring, polys_consumed, is_outer, bbox)))
44    }
45
46    /// Create a new PolyPath
47    pub fn new(
48        ring: Vec<P>,
49        polys_consumed: BTreeSet<usize>,
50        is_outer: bool,
51        bbox: Option<BBox>,
52    ) -> Self {
53        let bbox = bbox.unwrap_or(BBox::from_linestring(&ring));
54        let mut outer = None;
55        let mut holes = vec![];
56        if is_outer {
57            outer = Some(ring);
58        } else {
59            holes.push(ring);
60        }
61        PolyPath { id: 0, outer, old_outers: vec![], holes, polys_consumed, bbox }
62    }
63
64    /// Get the path as a vector polygon
65    ///
66    /// ## Returns
67    /// The resultant poly if it exists
68    pub fn get_path(&mut self) -> Option<Vec<Vec<P>>> {
69        self.outer.as_ref()?;
70        let outer = self.outer.as_mut().unwrap();
71        if outer.len() < 4 {
72            return None;
73        }
74        let mut res = vec![take(outer)];
75        for hole in self.holes.iter_mut() {
76            if hole.len() < 4 {
77                continue;
78            }
79            res.push(take(hole));
80        }
81        Some(res)
82    }
83}
84
85/// The next poly chunk
86#[derive(Debug, Clone, PartialEq)]
87pub struct NextRingChunk<P: FullXY> {
88    /// The intersection point
89    pub int_point: Point,
90    /// The next chunk
91    pub chunk: RingChunkRef<P>,
92}
93
94/// A path/piece/chunk from a polygon
95#[derive(Debug, Clone, PartialEq)]
96pub struct RingChunk<P: FullXY> {
97    /// If the chunk has been visited. Used by connection algorithms
98    pub visted: bool,
99    /// The index of the polygon it belongs to
100    pub poly_index: usize,
101    /// The index of the ring
102    pub ring_index: usize,
103    /// The bounding box
104    pub bbox: BBox,
105    /// Always stars with either the beginning of the poly ring OR an intersection point.
106    pub mid: Vec<P>,
107    /// The intersection point this chunk starts at
108    pub from: Point,
109    /// The intersection point this chunk ends at
110    pub to: Point,
111    /// can point to just one or multiple chunks. Many polys can touch the same point.
112    /// If none provided could be a start-end point
113    pub next: Option<NextRingChunk<P>>,
114    /// Precomputed from angle. Useful for when intersections are computed
115    pub from_angle: f64,
116    /// Precomputed to angle. Useful for when intersections are computed
117    pub to_angle: f64,
118}
119/// A reference to a RingChunk wrapped in an RC & RefCell
120pub type RingChunkRef<P> = Rc<RefCell<RingChunk<P>>>;
121impl<P: FullXY> RingChunk<P> {
122    /// Create a new RingChunk
123    #[allow(clippy::too_many_arguments)]
124    pub fn new(
125        poly_index: usize,
126        ring_index: usize,
127        bbox: BBox,
128        mid: Vec<P>,
129        from: Point,
130        to: Point,
131        from_angle: f64,
132        to_angle: f64,
133    ) -> Rc<RefCell<Self>> {
134        Rc::new(RefCell::new(RingChunk {
135            visted: false,
136            poly_index,
137            ring_index,
138            bbox,
139            mid,
140            from,
141            to,
142            next: None,
143            from_angle,
144            to_angle,
145        }))
146    }
147
148    /// Check if two chunks are equal
149    pub fn equal_chunk(&self, other: &RingChunkRef<P>) -> bool {
150        let other = &other.borrow();
151        (self.ring_index > 0) == (other.ring_index > 0)
152            && self.from == other.from
153            && self.to == other.to
154            && self.mid == other.mid
155    }
156}
157
158/// Intersection Point
159#[derive(Debug, Clone, PartialEq)]
160pub struct IntersectionPoint<P: FullXY> {
161    /// The intersection point
162    pub point: Point,
163    /// The chunks whose end point intersect this
164    pub from: Vec<RingChunkRef<P>>,
165    /// The chunks whose start point intersect this
166    pub to: Vec<RingChunkRef<P>>,
167}
168/// A reference to a IntersectionPoint wrapped in an RC & RefCell
169pub type IntersectionPointRef<P> = Rc<RefCell<IntersectionPoint<P>>>;
170impl<P: FullXY> IntersectionPoint<P> {
171    /// Create a new IntersectionPoint wrapped in an RC & RefCell
172    pub fn new(point: Point) -> IntersectionPointRef<P> {
173        Rc::new(RefCell::new(IntersectionPoint { point, from: vec![], to: vec![] }))
174    }
175}
176
177/// Intersection Lookup for chunks
178#[derive(Debug, Clone, PartialEq)]
179pub struct InterPointLookup<P: FullXY> {
180    /// The lookup
181    pub lookup: BTreeMap<Point, IntersectionPointRef<P>>,
182}
183impl<P: FullXY> Default for InterPointLookup<P> {
184    fn default() -> Self {
185        Self { lookup: BTreeMap::new() }
186    }
187}
188impl<P: FullXY> InterPointLookup<P> {
189    /// Get the intersection point
190    pub fn get(&mut self, point: Point) -> IntersectionPointRef<P> {
191        self.lookup.entry(point).or_insert_with(|| IntersectionPoint::new(point)).clone()
192    }
193
194    /// Link two points to eachother
195    #[allow(clippy::too_many_arguments)]
196    pub fn link_ints(
197        &mut self,
198        poly_index: usize,
199        ring_index: usize,
200        from: Point,
201        to: Point,
202        mid: Vec<P>,
203        from_angle: Option<f64>,
204        to_angle: Option<f64>,
205    ) -> RingChunkRef<P> {
206        // first build a chunk
207        let bbox = BBox::from_linestring(&mid).merge(&BBox::from_linestring(&[from, to]));
208        let from_angle =
209            from_angle.unwrap_or(angle(&mid.last().map(Point::from).unwrap_or(from), &to));
210        let to_angle =
211            to_angle.unwrap_or(angle(&mid.first().map(Point::from).unwrap_or(to), &from));
212        let chunk =
213            RingChunk::new(poly_index, ring_index, bbox, mid, from, to, from_angle, to_angle);
214        self.get(from).borrow_mut().to.push(chunk.clone());
215        self.get(to).borrow_mut().from.push(chunk.clone());
216        chunk
217    }
218}
219
220/// Build the PolyPaths and PolyChunks
221///
222/// ## Parameters
223/// - `vector_polygons`: the collection of polygons
224/// - `ring_intersect_lookup`: the ring intersection lookup for all rings in the multipolygon collection
225///
226/// ## Returns
227/// The PolyPaths, their lookups, and PolyChunks
228#[allow(clippy::type_complexity)]
229pub fn build_paths_and_chunks<P: FullXY>(
230    vector_polygons: &[Vec<Vec<P>>],
231    ring_intersect_lookup: &mut RingIntersectionLookup,
232) -> (
233    Vec<PolyPathRef<P>>,
234    BTreeMap<usize, PolyPathRef<P>>,
235    Vec<RingChunkRef<P>>,
236    InterPointLookup<P>,
237    Vec<BBox>,
238) {
239    // Setup result. Paths are the final structure of joined polygons.
240    let mut paths: Vec<PolyPathRef<P>> = vec![];
241    // Lookup is a helper for quickly finding paths in the future
242    let mut path_lookup: BTreeMap<usize, PolyPathRef<P>> = BTreeMap::new();
243    // let mut outer_ring_bboxes: VecBBox> = new Array(vectorPolygons.length);
244    let mut outer_ring_bboxes: Vec<BBox> = vec![Default::default(); vector_polygons.len()];
245
246    // 2) Build Poly Pieces
247    // If no intersections for the poly_index+Ring_index -> push as completed ring (into paths)
248    let mut chunks: Vec<RingChunkRef<P>> = vec![];
249    let mut int_lookup = InterPointLookup::<P>::default();
250    for (p_i, poly) in vector_polygons.iter().enumerate() {
251        for (r_i, ring) in poly.iter().enumerate() {
252            let intersections = ring_intersect_lookup.get_mut(&p_i).and_then(|r| r.get_mut(&r_i));
253            // Case 1: Insert into paths because it's already completed or expand existing path
254            if intersections.is_none() || intersections.as_ref().unwrap().is_empty() {
255                if let Some(existing_path) = path_lookup.get(&p_i) {
256                    let existing_path = &mut existing_path.borrow_mut();
257                    if r_i == 0 {
258                        existing_path.bbox.merge_in_place(&BBox::from_linestring(ring));
259                        existing_path.outer = Some(ring.clone());
260                        outer_ring_bboxes[p_i] = existing_path.bbox;
261                    } else {
262                        existing_path.holes.push(ring.clone());
263                    }
264                } else {
265                    let path =
266                        PolyPath::new_ref(ring.clone(), BTreeSet::from([p_i]), r_i == 0, None);
267                    if r_i == 0 {
268                        outer_ring_bboxes[p_i] = path.borrow().bbox;
269                    }
270                    path_lookup.insert(p_i, path.clone());
271                    paths.push(path);
272                }
273                continue;
274            }
275            // Case 2: Handle the intersections and build RingChunks
276            if r_i == 0 {
277                outer_ring_bboxes[p_i] = BBox::from_linestring(ring);
278            }
279            let intersections: Vec<&RingIntersection> =
280                intersections.as_ref().unwrap().iter().filter(|i| i.t != 0.).collect();
281            let mut curr_index = 0;
282            let mut int_index = 0;
283            while curr_index < ring.len() - 1 {
284                // if we are still working with intersections, build points with them
285                if let Some(cur_inter) = intersections.get(int_index) {
286                    // until we get to the next intersection, we link the points
287                    if curr_index != cur_inter.from {
288                        let start = curr_index;
289                        while curr_index != cur_inter.from {
290                            curr_index += 1;
291                        }
292                        let mid = ring[start + 1..curr_index].to_vec();
293                        let from = &ring[start];
294                        let to = &ring[curr_index];
295                        let chunk =
296                            int_lookup.link_ints(p_i, r_i, from.into(), to.into(), mid, None, None);
297                        chunks.push(chunk);
298                    }
299                    // now build links with the intersections until we get to the next intersection that isn't the same index
300                    let mut from: Point = (&ring[curr_index]).into();
301                    let mut cur_int = Some(cur_inter);
302                    while let Some(c_i) = cur_int
303                        && c_i.from == curr_index
304                    {
305                        // NOTE: For robustness, we have to store the angles we found when studying the intersections.
306                        // We make decisions about the polygons during the analysis of the intersections using
307                        // robust predicates. otherwise we would actually compute slightly different angles
308                        // that could percieve the intersection lines as swapped (non-existent).
309                        if from != *c_i.point.borrow() {
310                            chunks.push(int_lookup.link_ints(
311                                p_i,
312                                r_i,
313                                from,
314                                *c_i.point.borrow(),
315                                vec![],
316                                Some(invert_angle(c_i.t_angle)),
317                                Some(c_i.t_angle),
318                            ));
319                        }
320                        int_index += 1;
321                        from = *c_i.point.borrow();
322                        cur_int = intersections.get(int_index);
323                    }
324                    // if the intersection t is not 1, then we need to link the point to the end of the curr_index
325                    let next: Point = (&ring[curr_index + 1]).into();
326                    if from != next {
327                        let ang = intersections[int_index - 1].t_angle;
328                        chunks.push(int_lookup.link_ints(
329                            p_i,
330                            r_i,
331                            from,
332                            next,
333                            vec![],
334                            Some(invert_angle(ang)),
335                            Some(ang),
336                        ));
337                    }
338                } else {
339                    // no intersection, just build the point
340                    chunks.push(int_lookup.link_ints(
341                        p_i,
342                        r_i,
343                        (&ring[curr_index]).into(),
344                        (&ring[curr_index + 1]).into(),
345                        vec![],
346                        None,
347                        None,
348                    ));
349                }
350                curr_index += 1;
351            }
352        }
353    }
354
355    // sort the chunks by leftmost bboxes then bottom most if there is a tie
356    chunks.sort_by(|a, b| {
357        let BBox { left: a_left, bottom: a_bottom, .. } = a.borrow().bbox;
358        let BBox { left: b_left, bottom: b_bottom, .. } = b.borrow().bbox;
359        a_left
360            .partial_cmp(&b_left)
361            .unwrap_or(Ordering::Equal)
362            .then(a_bottom.partial_cmp(&b_bottom).unwrap_or(Ordering::Equal))
363    });
364
365    (paths, path_lookup, chunks, int_lookup, outer_ring_bboxes)
366}
367
368#[derive(Debug)]
369struct IntPair<P: FullXY> {
370    from: RingChunkRef<P>,
371    to: RingChunkRef<P>,
372    angle: f64,
373}
374
375/// Given an of intersection, find the best way to connect the from->to chunks
376///
377/// ## Parameters
378/// - `intersection`: the intersection to analyze
379pub fn merge_intersection_pairs<P: FullXY>(intersection: &IntersectionPointRef<P>) {
380    let IntersectionPoint { from, to, point, .. } = &mut *intersection.borrow_mut();
381    let int_point = *point;
382    if from.is_empty() || to.is_empty() {
383        return;
384    }
385    if from.len() == 1 && to.len() == 1 {
386        // connect the two chunks and move on
387        from[0].borrow_mut().next = Some(NextRingChunk { chunk: to[0].clone(), int_point });
388        return;
389    }
390
391    // remove "duplicate"/"same" chunks
392    let mut froms: Vec<RingChunkRef<P>> = vec![];
393    for c in from {
394        if c.borrow().visted {
395            continue;
396        }
397        if !froms.iter().any(|r| r.borrow().equal_chunk(c)) {
398            froms.push(c.clone());
399        } else {
400            c.borrow_mut().visted = true;
401        }
402    }
403    let mut tos: Vec<RingChunkRef<P>> = vec![];
404    for c in to {
405        if c.borrow().visted {
406            continue;
407        }
408        if !tos.iter().any(|r| r.borrow().equal_chunk(c)) {
409            tos.push(c.clone());
410        } else {
411            c.borrow_mut().visted = true;
412        }
413    }
414
415    // build pairs
416    let mut pairs: Vec<IntPair<P>> = vec![];
417    for f in froms.iter() {
418        for t in tos.iter() {
419            let mut angle = t.borrow().to_angle - f.borrow().from_angle;
420            angle = if angle < 0. { angle + TAU } else { angle };
421            pairs.push(IntPair { from: f.clone(), to: t.clone(), angle });
422        }
423    }
424    pairs.sort_by(|a, b| a.angle.partial_cmp(&b.angle).unwrap_or(Ordering::Equal));
425
426    for IntPair { from, to, .. } in &pairs {
427        let from = &mut from.borrow_mut();
428        if from.visted || to.borrow().visted {
429            continue;
430        }
431        from.next = Some(NextRingChunk { chunk: to.clone(), int_point });
432        from.visted = true;
433        to.borrow_mut().visted = true;
434    }
435
436    // cleanup visited
437    for f in froms {
438        f.borrow_mut().visted = false;
439    }
440    for t in tos {
441        t.borrow_mut().visted = false;
442    }
443}
444
445/// Returns the absolute angle between points A->B->C
446///
447/// ## Parameters
448/// - `a`: First point
449/// - `b`: Second Point
450///
451/// ## Returns
452/// Angle in degrees [-PI, PI]
453fn angle<P: FullXY, Q: FullXY>(a: &P, b: &Q) -> f64 {
454    atan2(a.y() - b.y(), a.x() - b.x())
455}
456
457/// Returns the absolute angle between points A->B->C
458///
459/// ## Parameters
460/// - `angle`: Angle in degrees [-PI, PI]
461///
462/// ## Returns
463/// Angle in degrees [-PI, PI]
464fn invert_angle(angle: f64) -> f64 {
465    if angle >= 0. { angle - PI } else { angle + PI }
466}
gistools/geometry/tools/polys/path_builder.rs

gistools/geometry/tools/polys/
path_builder.rs
