jagua_rs/geometry/fail_fast/
pole.rs

1use std::collections::VecDeque;
2
3use crate::geometry::geo_traits::{CollidesWith, DistanceTo, SeparationDistance};
4use crate::geometry::primitives::Circle;
5use crate::geometry::primitives::Rect;
6use crate::geometry::primitives::SPolygon;
7
8/// Computes the *pole* - the largest circle which is both inside of `shape` while being outside all other `poles`.
9/// Closely related to [Pole of Inaccessibility (PoI)](https://en.wikipedia.org/wiki/Pole_of_inaccessibility),
10/// and inspired by Mapbox's [`polylabel`](https://github.com/mapbox/polylabel) algorithm.
11pub fn compute_pole(shape: &SPolygon, poles: &[Circle]) -> Circle {
12    let square_bbox = shape.bbox.inflate_to_square();
13    let root = POINode::new(square_bbox, MAX_POI_TREE_DEPTH, shape, poles);
14    let mut queue = VecDeque::from([root]);
15    let mut best: Option<Circle> = None;
16    let distance = |circle: &Option<Circle>| circle.as_ref().map_or(0.0, |c| c.radius);
17
18    while let Some(node) = queue.pop_front() {
19        //check if better than current best
20        if node.distance > distance(&best) {
21            best = Some(Circle::try_new(node.bbox.centroid(), node.distance).unwrap());
22        }
23
24        //see if worth it to split
25        if node.distance_upperbound() > distance(&best) {
26            if let Some(children) = node.split(shape, poles) {
27                queue.extend(children);
28            }
29        }
30    }
31    best.expect("no pole present")
32}
33
34///Generates a set of 'poles' for a shape according to specified coverage limits.
35///See [`compute_pole`] for details on what a 'pole' is.
36pub fn generate_surrogate_poles(shape: &SPolygon, n_pole_limits: &[(usize, f32)]) -> Vec<Circle> {
37    let mut all_poles = vec![shape.poi];
38    let mut total_pole_area = shape.poi.area();
39
40    //Generate the poles until one of the pole number / coverage limits is reached
41    loop {
42        let next = compute_pole(shape, &all_poles);
43
44        total_pole_area += next.area();
45        all_poles.push(next);
46
47        let current_coverage = total_pole_area / shape.area;
48
49        //check if any limit in the number of poles is reached at this coverage
50        let active_pole_limit = n_pole_limits
51            .iter()
52            .filter(|(_, coverage_threshold)| current_coverage > *coverage_threshold)
53            .min_by_key(|(n_poles, _)| *n_poles)
54            .map(|(n_poles, _)| n_poles);
55
56        if let Some(active_pole_limit) = active_pole_limit {
57            if all_poles.len() >= *active_pole_limit {
58                //stop generating if we are above the limit
59                break;
60            }
61        }
62        assert!(
63            all_poles.len() < 1000,
64            "More than 1000 poles were generated, please check the SPSurrogateConfig"
65        )
66    }
67    all_poles
68}
69
70const MAX_POI_TREE_DEPTH: usize = 10;
71
72struct POINode {
73    pub level: usize,
74    pub bbox: Rect,
75    pub radius: f32,
76    pub distance: f32,
77}
78
79impl POINode {
80    fn new(bbox: Rect, level: usize, poly: &SPolygon, poles: &[Circle]) -> Self {
81        let radius = bbox.diameter() / 2.0;
82
83        let centroid_inside = poly.collides_with(&bbox.centroid())
84            && poles.iter().all(|c| !c.collides_with(&bbox.centroid()));
85
86        let distance = {
87            let distance_to_edges = poly.edge_iter().map(|e| e.distance_to(&bbox.centroid()));
88
89            let distance_to_poles = poles
90                .iter()
91                .map(|c| c.separation_distance(&bbox.centroid()).1);
92
93            let distance_to_border = distance_to_edges
94                .chain(distance_to_poles)
95                .fold(f32::MAX, |acc, d| acc.min(d));
96
97            //if the centroid is outside, distance is counted negative
98            match centroid_inside {
99                true => distance_to_border,
100                false => -distance_to_border,
101            }
102        };
103
104        Self {
105            bbox,
106            level,
107            radius,
108            distance,
109        }
110    }
111
112    fn split(&self, poly: &SPolygon, poles: &[Circle]) -> Option<[POINode; 4]> {
113        match self.level {
114            0 => None,
115            _ => Some(
116                self.bbox
117                    .quadrants()
118                    .map(|qd| POINode::new(qd, self.level - 1, poly, poles)),
119            ),
120        }
121    }
122
123    fn distance_upperbound(&self) -> f32 {
124        self.radius + self.distance
125    }
126}
jagua_rs/geometry/fail_fast/pole.rs

jagua_rs/geometry/fail_fast/
pole.rs
