jagua_rs/geometry/fail_fast/
poi.rs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
use std::collections::VecDeque;

use ordered_float::NotNan;

use crate::fsize;
use crate::geometry::geo_traits::{CollidesWith, Distance, SeparationDistance, Shape};
use crate::geometry::primitives::aa_rectangle::AARectangle;
use crate::geometry::primitives::circle::Circle;
use crate::geometry::primitives::simple_polygon::SimplePolygon;

/// Generates the Pole of Inaccessibility (PoI). The PoI is the point in the interior of the shape that is farthest from the boundary.
/// The interior is defined as the interior of the `shape` minus the interior of the `poles`.
pub fn generate_next_pole(shape: &SimplePolygon, poles: &[Circle]) -> Circle {
    //Based on Mapbox's "Polylabel" algorithm: <https://github.com/mapbox/polylabel>
    let square_bbox = shape.bbox().inflate_to_square();
    let root = POINode::new(square_bbox, MAX_POI_TREE_DEPTH, shape, poles);
    let mut queue = VecDeque::from([root]);
    let mut best: Option<Circle> = None;
    let distance = |circle: &Option<Circle>| circle.as_ref().map_or(0.0, |c| c.radius);

    while let Some(node) = queue.pop_front() {
        //check if better than current best
        if node.distance > distance(&best) {
            best = Some(Circle::new(node.bbox.centroid(), node.distance));
        }

        //see if worth it to split
        if node.distance_upperbound() > distance(&best) {
            if let Some(children) = node.split(shape, poles) {
                queue.extend(children);
            }
        }
    }
    best.expect("no pole present")
}

///Generates additional poles for a shape alongside the PoI
pub fn generate_additional_surrogate_poles(
    shape: &SimplePolygon,
    n_pole_limits: &[(usize, fsize)],
) -> Vec<Circle> {
    //generate the additional poles
    let additional_poles = {
        let mut all_poles = vec![shape.poi.clone()];
        let mut total_pole_area = shape.poi.area();

        //Generate the poles
        loop {
            let next = generate_next_pole(shape, &all_poles);

            total_pole_area += next.area();
            all_poles.push(next);

            let current_coverage = total_pole_area / shape.area();

            //check if any limit in the number of poles is reached at this coverage
            let active_pole_limit = n_pole_limits
                .iter()
                .filter(|(_, coverage_threshold)| current_coverage > *coverage_threshold)
                .min_by_key(|(n_poles, _)| *n_poles)
                .map(|(n_poles, _)| n_poles);

            if let Some(active_pole_limit) = active_pole_limit {
                if all_poles.len() >= *active_pole_limit {
                    //stop generating if we are above the limit
                    break;
                }
            }
            assert!(
                all_poles.len() < 1000,
                "More than 1000 poles were generated, please check the SPSurrogateConfig"
            )
        }
        all_poles[1..].to_vec()
    };

    //Sort the poles to maximize chance of early fail fast
    let mut unsorted_poles = additional_poles.clone();
    let mut sorted_poles = vec![];

    while !unsorted_poles.is_empty() {
        //find the pole that maximizes the distance to the closest prior pole and is also large
        let (next_pole_index, _) = unsorted_poles
            .iter()
            .enumerate()
            .map(|(i, p)| {
                let prior_poles = sorted_poles.iter().chain([&shape.poi]);

                let min_distance_prior_poles = prior_poles
                    .map(|prior| prior.separation_distance(&p.centroid()).1)
                    .min_by(|a, b| a.partial_cmp(b).unwrap())
                    .unwrap();
                (i, p.radius.powi(2) * min_distance_prior_poles)
            })
            .max_by_key(|(_i, obj)| NotNan::new(*obj).unwrap())
            .unwrap();
        sorted_poles.push(unsorted_poles.remove(next_pole_index));
    }

    sorted_poles
}

const MAX_POI_TREE_DEPTH: usize = 10;

struct POINode {
    pub level: usize,
    pub bbox: AARectangle,
    pub radius: fsize,
    pub distance: fsize,
}

impl POINode {
    pub fn new(bbox: AARectangle, level: usize, poly: &SimplePolygon, poles: &[Circle]) -> Self {
        let radius = bbox.diameter() / 2.0;

        let centroid_inside = poly.collides_with(&bbox.centroid())
            && poles.iter().all(|c| !c.collides_with(&bbox.centroid()));

        let distance = {
            let distance_to_edges = poly.edge_iter().map(|e| e.distance(&bbox.centroid()));

            let distance_to_poles = poles
                .iter()
                .map(|c| c.separation_distance(&bbox.centroid()).1);

            let distance_to_border = distance_to_edges
                .chain(distance_to_poles)
                .fold(fsize::MAX, |acc, d| acc.min(d));

            //if the centroid is outside, distance is counted negative
            match centroid_inside {
                true => distance_to_border,
                false => -distance_to_border,
            }
        };

        Self {
            bbox,
            level,
            radius,
            distance,
        }
    }

    pub fn split(&self, poly: &SimplePolygon, poles: &[Circle]) -> Option<[POINode; 4]> {
        match self.level {
            0 => None,
            _ => Some(
                self.bbox
                    .quadrants()
                    .map(|qd| POINode::new(qd, self.level - 1, poly, poles)),
            ),
        }
    }

    pub fn distance_upperbound(&self) -> fsize {
        self.radius + self.distance
    }
}