gistools/tools/
polylabel.rs

1use crate::data_structures::PriorityQueue;
2use alloc::vec::Vec;
3use core::f64::consts::SQRT_2;
4use libm::{fmax, fmin, sqrt};
5use s2json::{GetXY, MValueCompatible, NewXYM};
6
7/// The metadata inserted into the Vector Feature
8#[derive(Debug, Default, Copy, MValueCompatible, Clone, PartialEq)]
9pub struct PolyLabelMetadata {
10    /// The distance to the label
11    pub distance: f64,
12}
13impl PolyLabelMetadata {
14    /// Create a new PolyLabelMetadata
15    pub fn new(distance: f64) -> Self {
16        PolyLabelMetadata { distance }
17    }
18}
19
20/// A cell in the polygon label algorithm
21#[derive(Debug, Default, Clone)]
22pub struct PolyLabelCell {
23    /// cell center x
24    pub x: f64,
25    /// cell center y
26    pub y: f64,
27    /// half the cell size
28    pub h: f64,
29    /// distance from cell center to polygon
30    pub d: f64,
31    /// max distance to polygon within a cell
32    pub max: f64,
33}
34
35/// # Polylabels
36///
37/// ## Description
38/// Find the labels for a collection of vector polygons
39///
40/// ## Usage
41///
42/// ```rust
43/// use gistools::tools::{PolyLabelMetadata, polylabels};
44/// use s2json::VectorPoint;
45///
46/// let data: Vec<Vec<Vec<VectorPoint<()>>>> = vec![vec![vec![]]];
47/// let emp: Vec<VectorPoint<PolyLabelMetadata>> = polylabels(&data, None);
48///
49/// assert_eq!(emp, vec![VectorPoint::new_xy(0., 0., Some(PolyLabelMetadata::new(0.)))]);
50/// ```
51///
52/// ## Links
53/// - <https://sites.google.com/site/polesofinaccessibility/>
54///
55/// ## Parameters
56/// - `polygons`: the vector multi-polygon to find the label for
57/// - `precision`: the precision of the label
58///
59/// ## Returns
60/// The collection of label positions and the distances to the labels
61pub fn polylabels<P: GetXY, R: NewXYM<PolyLabelMetadata>>(
62    polygons: &[Vec<Vec<P>>],
63    precision: Option<f64>,
64) -> Vec<R> {
65    polygons.iter().map(|polygon| polylabel(polygon, precision)).collect()
66}
67
68/// # Polylabel
69///
70/// ## Description
71/// Find the label for a vector polygon
72///
73/// ## Usage
74///
75/// ```rust
76/// use gistools::tools::{PolyLabelMetadata, polylabel};
77/// use s2json::VectorPoint;
78///
79/// let data: Vec<Vec<VectorPoint<()>>> = vec![vec![]];
80/// let emp: VectorPoint<PolyLabelMetadata> = polylabel(&data, None);
81///
82/// assert_eq!(emp, VectorPoint::new_xy(0., 0., Some(PolyLabelMetadata::new(0.))));
83/// ```
84///
85/// ## Links
86/// - <https://sites.google.com/site/polesofinaccessibility/>
87///
88/// ## Parameters
89/// - `polygon`: the vector polygon to find the label for
90/// - `precision`: the precision of the label
91///
92/// ## Returns
93/// The label position and the distance to the label
94pub fn polylabel<P: GetXY, R: NewXYM<PolyLabelMetadata>>(
95    polygon: &Vec<Vec<P>>,
96    precision: Option<f64>,
97) -> R {
98    let precision = precision.unwrap_or(1.0);
99    // find the bounding box of the outer ring
100    let mut min_x = f64::MAX;
101    let mut min_y = f64::MAX;
102    let mut max_x = f64::MIN;
103    let mut max_y = f64::MIN;
104
105    if polygon.is_empty() || polygon[0].is_empty() {
106        return R::new_xym(0.0, 0.0, PolyLabelMetadata::default());
107    }
108
109    for point in &polygon[0] {
110        let (x, y) = point.xy();
111        if x < min_x {
112            min_x = x;
113        }
114        if y < min_y {
115            min_y = y;
116        }
117        if x > max_x {
118            max_x = x;
119        }
120        if y > max_y {
121            max_y = y;
122        }
123    }
124
125    let width = max_x - min_x;
126    let height = max_y - min_y;
127    let cell_size = fmax(precision, fmin(width, height));
128
129    if cell_size == precision {
130        return R::new_xym(min_x, min_y, PolyLabelMetadata::default());
131    }
132
133    // a priority queue of cells in order of their "potential" (max distance to polygon)
134    let mut cell_queue =
135        PriorityQueue::<PolyLabelCell>::new(|a: &PolyLabelCell, b: &PolyLabelCell| {
136            b.max.partial_cmp(&a.max).unwrap_or(core::cmp::Ordering::Equal)
137        });
138
139    // take centroid as the first best guess
140    let mut best_cell = get_centroid_cell(polygon);
141
142    // second guess: bounding box centroid
143    let bbox_cell = build_cell(min_x + width / 2., min_y + height / 2., 0., polygon);
144    if bbox_cell.d > best_cell.d {
145        best_cell = bbox_cell;
146    }
147
148    // add a cell to the queue
149    let potentially_queue =
150        |x: f64,
151         y: f64,
152         h: f64,
153         best_cell: &mut PolyLabelCell,
154         cell_queue: &mut PriorityQueue<PolyLabelCell>| {
155            let cell = build_cell(x, y, h, polygon);
156            if cell.max > best_cell.d + precision {
157                cell_queue.push(cell.clone());
158            }
159            // update the best cell if we found a better one
160            if cell.d > best_cell.d {
161                *best_cell = cell;
162            }
163        };
164
165    // cover polygon with initial cells
166    let mut h = cell_size / 2.;
167    let mut x = min_x;
168    while x < max_x {
169        let mut y = min_y;
170        while y < max_y {
171            potentially_queue(x + h, y + h, h, &mut best_cell, &mut cell_queue);
172            y += cell_size;
173        }
174        x += cell_size;
175    }
176
177    loop {
178        // pick the most promising cell from the queue
179        let cell = cell_queue.pop();
180        if cell.is_none() {
181            break;
182        }
183        let PolyLabelCell { max, x, y, h: ch, .. } = &cell.unwrap();
184
185        // do not drill down further if there's no chance of a better solution
186        if max - best_cell.d <= precision {
187            break;
188        }
189
190        // split the cell into four cells
191        h = ch / 2.;
192        potentially_queue(x - h, y - h, h, &mut best_cell, &mut cell_queue);
193        potentially_queue(x + h, y - h, h, &mut best_cell, &mut cell_queue);
194        potentially_queue(x - h, y + h, h, &mut best_cell, &mut cell_queue);
195        potentially_queue(x + h, y + h, h, &mut best_cell, &mut cell_queue);
196    }
197
198    R::new_xym(best_cell.x, best_cell.y, PolyLabelMetadata::default())
199}
200
201/// build a cell
202///
203/// ## Parameters
204/// - `x`: the cell x coordinate
205/// - `y`: the cell y coordinate
206/// - `h`: half the cell size
207/// - `polygon`: the vector polygon
208///
209/// ## Returns
210/// The cell
211fn build_cell<P: GetXY>(x: f64, y: f64, h: f64, polygon: &Vec<Vec<P>>) -> PolyLabelCell {
212    let d = point_to_polygon_dist(x, y, polygon);
213    PolyLabelCell { x, y, h, d, max: d + h * SQRT_2 }
214}
215
216/// signed distance from point to polygon outline (negative if point is outside)
217///
218/// ## Parameters
219/// - `x`: the point x coordinate
220/// - `y`: the point y coordinate
221/// - `polygon`: the vector polygon to check
222///
223/// ## Returns
224/// The signed distance
225fn point_to_polygon_dist<P: GetXY>(x: f64, y: f64, polygon: &Vec<Vec<P>>) -> f64 {
226    let mut inside = false;
227    let mut min_dist_sq = f64::MAX;
228
229    for ring in polygon {
230        let len = ring.len();
231        let mut j = len - 1;
232        for i in 0..len {
233            let a = &ring[i];
234            let b = &ring[j];
235
236            if (a.y() > y) != (b.y() > y)
237                && x < ((b.x() - a.x()) * (y - a.y())) / (b.y() - a.y()) + a.x()
238            {
239                inside = !inside;
240            }
241
242            min_dist_sq = fmin(min_dist_sq, get_seg_dist_sq(x, y, a, b));
243            j = i; // Update j to the previous i (j = i++)
244        }
245    }
246
247    if min_dist_sq == 0. { 0. } else { (if inside { 1. } else { -1. }) * sqrt(min_dist_sq) }
248}
249
250/// get polygon centroid
251/// return the centroid as a cell
252fn get_centroid_cell<P: GetXY>(polygon: &Vec<Vec<P>>) -> PolyLabelCell {
253    let mut area = 0.;
254    let mut x = 0.;
255    let mut y = 0.;
256    let points = &polygon[0];
257
258    let len = points.len();
259    let mut j = len - 1;
260    for i in 0..len {
261        let a = &points[i];
262        let b = &points[j];
263        let f = a.x() * b.y() - b.x() * a.y();
264        x += (a.x() + b.x()) * f;
265        y += (a.y() + b.y()) * f;
266        area += f * 3.;
267        j = i; // Update j to the previous i (j = i++)
268    }
269    let centroid = build_cell(x / area, y / area, 0., polygon);
270    if area == 0. || centroid.d < 0. {
271        build_cell(points[0].x(), points[0].y(), 0., polygon)
272    } else {
273        centroid
274    }
275}
276
277/// get squared distance from a point to a segment AB
278/// return the squared distance
279fn get_seg_dist_sq<P: GetXY>(px: f64, py: f64, a: &P, b: &P) -> f64 {
280    let mut x = a.x();
281    let mut y = a.y();
282    let mut dx = b.x() - x;
283    let mut dy = b.y() - y;
284
285    if dx != 0. || dy != 0. {
286        let t = ((px - x) * dx + (py - y) * dy) / (dx * dx + dy * dy);
287
288        if t > 1. {
289            x = b.x();
290            y = b.y();
291        } else if t > 0. {
292            x += dx * t;
293            y += dy * t;
294        }
295    }
296
297    dx = px - x;
298    dy = py - y;
299
300    dx * dx + dy * dy
301}
gistools/tools/polylabel.rs

gistools/tools/
polylabel.rs
