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

jagua_rs/geometry/fail_fast/
pole.rs
