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

gistools/geometry/tools/polys/
path_builder.rs
