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, f32::min);
105
106            //if the centroid is outside, distance is counted negative
107            if centroid_inside {
108                distance_to_border
109            } else {
110                -distance_to_border
111            }
112        };
113
114        Self {
115            level,
116            bbox,
117            radius,
118            distance,
119        }
120    }
121
122    fn split(&self, poly: &SPolygon, poles: &[Circle]) -> Option<[POINode; 4]> {
123        match self.level {
124            0 => None,
125            _ => Some(
126                self.bbox
127                    .quadrants()
128                    .map(|qd| POINode::new(qd, self.level - 1, poly, poles)),
129            ),
130        }
131    }
132
133    fn distance_upperbound(&self) -> f32 {
134        self.radius + self.distance
135    }
136}
jagua_rs/geometry/fail_fast/pole.rs

jagua_rs/geometry/fail_fast/
pole.rs
