gistools/geometry/tools/polys/
dekink.rs

1use crate::geometry::{
2    Area, InterPointLookup, NextRingChunk, PolyPath, PolyPathRef, RingChunkRef,
3    RingIntersectionLookup, Segment, build_paths_and_chunks, merge_intersection_pairs,
4    polygons_intersections_lookup, polyline_in_polyline,
5};
6use alloc::{
7    collections::{BTreeMap, BTreeSet},
8    vec,
9    vec::Vec,
10};
11use core::{cmp::Ordering, mem::take};
12use s2json::{BBox, FullXY};
13
14/// Given a collection of polygons, if any of the polygons are kinked, dekink them
15///
16/// ## Parameters
17/// `polygons`: the polygons are from either a VectorFeature, VectorPolygonGeometry, or raw VectorPolygon
18///
19/// ## Returns
20/// The dekinked polygons
21pub fn dekink_polygons<P: FullXY>(polygons: &[Vec<Vec<P>>]) -> Option<(Vec<Vec<Vec<P>>>, BBox)> {
22    polygons.iter().filter_map(|p| dekink_polygon(p)).try_fold(
23        (vec![], BBox::default()),
24        |(mut res, mut bbox), mut p| {
25            res.append(&mut p.0);
26            bbox.merge_in_place(&p.1);
27            Some((res, bbox))
28        },
29    )
30}
31
32/// Given a polygon, if the polygon is kinked, dekink it
33///
34/// ## Parameters
35/// - `polygon`: the polygon as either a VectorFeature, VectorPolygonGeometry, or raw VectorPolygon
36///
37/// ## Returns
38/// The dekinked polygon
39pub fn dekink_polygon<P: FullXY>(polygon: &[Vec<P>]) -> Option<(Vec<Vec<Vec<P>>>, BBox)> {
40    // not enough data, just clone
41    if polygon.is_empty() {
42        return None;
43    }
44    let vector_polygons = vec![polygon.to_vec()];
45
46    // 1) build intersections `[poly_index][ring_index] -> Intersections`. Store where on the ring other rings intersect
47    let mut ring_int_lookup: RingIntersectionLookup = polygons_intersections_lookup(
48        &vector_polygons,
49        Some(|seg1: &Segment, seg2: &Segment| -> bool {
50            // if same id ignore
51            seg2.id != seg1.id &&
52        // only pass forward not backward
53        seg2.id > seg1.id &&
54        // TODO: At some point intersections of inner rings against the outer ring should be considered.
55        // For now the ring_index must be the same, polygonsIntersectionsLookup should return
56        // two problem sets down the road, one for cleaning individual rings, and one for fixing holes that go out of bounds
57        seg2.ring_index == seg1.ring_index
58        }),
59    );
60
61    // 2) Build Poly Pieces
62    // Setup result paths with chunks that are the final structure of joined polygons.
63    // Lookup is a helper for quickly finding the right path in the future as paths can consume multiple polygons
64    // If no intersections for the poly_index+Ring_index -> it's immediately consumed into paths. Otherwise it's a chunk
65    let (mut paths, _, chnks, mut int_lookup, _b) =
66        build_paths_and_chunks(&vector_polygons, &mut ring_int_lookup);
67
68    // 3) Consume chunks into PolyPaths
69    // If no intersections for the poly_index+Ring_index -> push as completed ring
70    build_paths_from_chunks(&mut paths, &mut int_lookup, &chnks);
71
72    // 4) Convert PolyPaths into the resultant MultiPolygon
73    let coordinates: Vec<Vec<Vec<P>>> =
74        paths.iter_mut().filter_map(|p| p.borrow_mut().get_path()).collect();
75    let bbox = paths.iter().map(|p| p.borrow().bbox).fold(BBox::default(), |mut acc, b| {
76        acc.merge_in_place(&b);
77        acc
78    });
79    if coordinates.is_empty() { None } else { Some((coordinates, bbox)) }
80}
81
82/// Simplified ring with guide of how it was rebuild
83#[derive(Debug)]
84struct Ring<P: FullXY> {
85    linestring: Vec<P>,
86    is_ccw: bool,
87    is_hole: bool,
88    bbox: BBox,
89    area: f64, // bbox area
90}
91impl<P: FullXY> Ring<P> {
92    fn new(linestring: Vec<P>, is_ccw: bool, is_hole: bool, bbox: BBox, area: f64) -> Self {
93        Self { linestring, is_ccw, is_hole, bbox, area }
94    }
95}
96/// Collection of rings sorted by poly_index
97/// [poly_index: number]: Collection of Rings for that polygon
98type RingStore<P> = BTreeMap<usize, Vec<Ring<P>>>;
99
100/// Given a set of chunks, build a set of paths
101///
102/// ## Parameters
103/// - `paths`: a set of paths to add to
104/// - `intersections`: all intersections
105/// - `chunks`: a set of chunks
106fn build_paths_from_chunks<P: FullXY>(
107    paths: &mut Vec<PolyPathRef<P>>,
108    intersections: &mut InterPointLookup<P>,
109    chunks: &Vec<RingChunkRef<P>>,
110) {
111    let mut ring_store: RingStore<P> = BTreeMap::new();
112    // store existing paths and reset.
113    for path in paths.iter() {
114        let PolyPath { outer, polys_consumed, bbox, holes, .. } = &mut *path.borrow_mut();
115        let poly_index = polys_consumed.first().unwrap();
116        let poly_store = ring_store.entry(*poly_index).or_insert(vec![]);
117        if let Some(outer) = outer {
118            poly_store.push(Ring::new(take(outer), true, false, *bbox, bbox.area()));
119        }
120        for hole in holes {
121            poly_store.push(Ring::new(take(hole), false, true, *bbox, bbox.area()));
122        }
123    }
124    paths.clear();
125    // for each intersections, connect all the from and to, smallest angle between from->to first slowly work your way through
126    for int in intersections.lookup.values() {
127        merge_intersection_pairs(int);
128    }
129    // run through all chunks, if unvisited, add to paths
130    for chunk in chunks {
131        if chunk.borrow().visted {
132            continue;
133        }
134        // follow along a chunk until we find our start point again
135        let start = chunk.borrow().from;
136        let mut curr_chunk: RingChunkRef<P> = chunk.clone();
137        let mut linestring: Vec<P> = vec![P::new_xy(start.0, start.1)];
138        let mut bbox = curr_chunk.borrow().bbox;
139        loop {
140            if curr_chunk.borrow().visted {
141                break;
142            }
143            curr_chunk.borrow_mut().visted = true;
144            linestring.append(&mut curr_chunk.borrow_mut().mid);
145            bbox.merge_in_place(&curr_chunk.borrow().bbox);
146            if curr_chunk.borrow().next.is_none() {
147                break;
148            }
149            let (next_chunk, int_point) = {
150                let chnk = &mut curr_chunk.borrow_mut();
151                let NextRingChunk { chunk, int_point } = chnk.next.as_ref().unwrap();
152                (chunk.clone(), *int_point)
153            };
154            linestring.push(P::new_xy(int_point.0, int_point.1));
155            curr_chunk = next_chunk;
156            if int_point == start {
157                break;
158            }
159        }
160        let area = (&linestring).area(Some(1.));
161        if area == 0. || linestring.len() < 4 || linestring.first() != linestring.last() {
162            continue;
163        }
164        // now build the path or add to an existing path
165        let is_ccw = area > 0.;
166        let is_hole = chunk.borrow().ring_index != 0;
167        // store in correct location
168        let poly_store = ring_store.entry(chunk.borrow().poly_index).or_insert(vec![]);
169        poly_store.push(Ring::new(linestring, is_ccw, is_hole, bbox, area));
170    }
171    // For each ring_store, build out polys and store them in paths
172    for rings in ring_store.values_mut() {
173        ring_set_to_paths(paths, rings);
174    }
175}
176
177/// Convert a set of rings into a set of paths
178///
179/// ## Parameters
180/// - `paths`: the collection of paths to store the results in
181/// - `ring_set`: the current set of rings re-built from a polygon
182fn ring_set_to_paths<P: FullXY>(paths: &mut Vec<PolyPathRef<P>>, ring_set: &mut [Ring<P>]) {
183    if ring_set.is_empty() {
184        return;
185    }
186    // sort by bbox area desc and prep real outers store
187    ring_set.sort_by(|a, b| b.area.partial_cmp(&a.area).unwrap_or(Ordering::Equal));
188    let mut actual_outers: Vec<PolyPathRef<P>> = vec![];
189    // filter by case type
190
191    // store all true outers
192    {
193        let mut outers: Vec<&mut Ring<P>> =
194            ring_set.iter_mut().filter(|r| !r.is_hole && r.is_ccw).collect();
195        // store all true outers. We know the first outer is the largest original outer ring.
196        // The future outers are holes either kink inside or outside the original. Only store kinks that are outside the original
197        if !outers.is_empty() {
198            let first = outers.remove(0);
199            // store the first one without the ring we still need it for now
200            actual_outers.push(PolyPath::new_ref(vec![], BTreeSet::new(), true, Some(first.bbox)));
201            // check all the others are not in the first
202            for outer in outers {
203                if outer.bbox.inside(&first.bbox)
204                    && polyline_in_polyline(&outer.linestring, &first.linestring)
205                {
206                    continue;
207                }
208                actual_outers.push(PolyPath::new_ref(
209                    take(&mut outer.linestring),
210                    BTreeSet::new(),
211                    true,
212                    Some(outer.bbox),
213                ));
214            }
215            // now store the first one's ring
216            actual_outers[0].borrow_mut().outer = Some(take(&mut first.linestring));
217        }
218    }
219
220    // If outer in `outers_maybe_hole` is inside an actual outer, it's a hole; Otherwise it's another outer
221    {
222        let outers_maybe_hole: Vec<&mut Ring<P>> =
223            ring_set.iter_mut().filter(|r| !r.is_hole && !r.is_ccw).collect();
224        for Ring { linestring, bbox, .. } in outers_maybe_hole {
225            let mut found = false;
226            for actual_outer in actual_outers.iter() {
227                let actual_outer = &mut actual_outer.borrow_mut();
228                if bbox.inside(&actual_outer.bbox)
229                    && polyline_in_polyline(linestring, actual_outer.outer.as_ref().unwrap())
230                {
231                    // store the hole in this outer
232                    actual_outer.holes.push(take(linestring));
233                    found = true;
234                    break;
235                }
236            }
237            if found {
238                continue;
239            }
240            // otherwise, it's a new outer
241            linestring.reverse();
242            actual_outers.push(PolyPath::new_ref(
243                take(linestring),
244                BTreeSet::new(),
245                true,
246                Some(*bbox),
247            ))
248        }
249    }
250
251    // now organize holes
252    {
253        let holes: Vec<&mut Ring<P>> =
254            ring_set.iter_mut().filter(|r| r.is_hole && !r.is_ccw).collect();
255        if !actual_outers.is_empty() {
256            for Ring { linestring, bbox, .. } in holes {
257                if actual_outers.len() == 1 {
258                    actual_outers[0].borrow_mut().holes.push(take(linestring));
259                } else {
260                    // find the outer this hole belongs to
261                    for actual_outer in actual_outers.iter() {
262                        if bbox.inside(&actual_outer.borrow().bbox)
263                            && polyline_in_polyline(
264                                linestring,
265                                actual_outer.borrow().outer.as_ref().unwrap(),
266                            )
267                        {
268                            // store the hole in this outer
269                            actual_outer.borrow_mut().holes.push(take(linestring));
270                            break;
271                        }
272                    }
273                }
274            }
275        }
276    }
277
278    // Now store all actual_outers we built
279    for actual_outer in actual_outers.iter() {
280        paths.push(actual_outer.clone());
281    }
282}
gistools/geometry/tools/polys/dekink.rs

gistools/geometry/tools/polys/
dekink.rs
