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///Generates a set of 'poles' for a shape according to specified coverage limits.
11///See [`compute_pole`] for details on what a 'pole' is.
12pub fn generate_surrogate_poles(
13    shape: &SPolygon,
14    n_pole_limits: &[(usize, f32)],
15) -> Result<Vec<Circle>> {
16    let mut all_poles = vec![shape.poi];
17    let mut total_pole_area = shape.poi.area();
18
19    //Generate the poles until one of the pole number / coverage limits is reached
20    loop {
21        let next = compute_pole(shape, &all_poles)?;
22
23        total_pole_area += next.area();
24        all_poles.push(next);
25
26        let current_coverage = total_pole_area / shape.area;
27
28        //check if any limit in the number of poles is reached at this coverage
29        let active_pole_limit = n_pole_limits
30            .iter()
31            .filter(|(_, coverage_threshold)| current_coverage > *coverage_threshold)
32            .min_by_key(|(n_poles, _)| *n_poles)
33            .map(|(n_poles, _)| n_poles);
34
35        if let Some(active_pole_limit) = active_pole_limit
36            && all_poles.len() >= *active_pole_limit
37        {
38            //stop generating if we are above the limit
39            break;
40        }
41        assert!(
42            all_poles.len() < 1000,
43            "More than 1000 poles were generated, please check the SPSurrogateConfig"
44        )
45    }
46    Ok(all_poles)
47}
48
49/// Computes the *pole* - the largest circle which is both inside of `shape` while being outside all other `poles`.
50/// Closely related to [Pole of Inaccessibility (PoI)](https://en.wikipedia.org/wiki/Pole_of_inaccessibility),
51/// and inspired by Mapbox's [`polylabel`](https://github.com/mapbox/polylabel) algorithm.
52pub fn compute_pole(shape: &SPolygon, poles: &[Circle]) -> Result<Circle> {
53    let square_bbox = shape.bbox.inflate_to_square();
54    let root = POINode::new(square_bbox, MAX_POI_TREE_DEPTH, shape, poles);
55    let mut queue = VecDeque::from([root]);
56    let mut best: Option<Circle> = None;
57    let distance = |circle: &Option<Circle>| circle.as_ref().map_or(0.0, |c| c.radius);
58
59    while let Some(node) = queue.pop_front() {
60        //check if better than current best
61        if node.distance > distance(&best) {
62            best = Some(Circle::try_new(node.bbox.centroid(), node.distance).unwrap());
63        }
64
65        //see if worth it to split
66        if node.distance_upperbound() > distance(&best)
67            && let Some(children) = node.split(shape, poles)
68        {
69            queue.extend(children);
70        }
71    }
72    best.ok_or(anyhow!(
73        "no pole found with {} levels of recursion. Please check the input shape: {:?}",
74        MAX_POI_TREE_DEPTH,
75        &shape.vertices
76    ))
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
