polylabel/
lib.rs

Help
1#![doc(
2    html_logo_url = "https://cdn.rawgit.com/urschrei/polylabel-rs/7a07336e85572eb5faaf0657c2383d7de5620cd8/ell.svg",
3    html_root_url = "https://docs.rs/polylabel-rs/"
4)]
5//! This crate provides a Rust implementation of the [Polylabel](https://github.com/mapbox/polylabel) algorithm
6//! for finding the optimum position of a polygon label.
7//!
8//! ffi bindings are provided: enable the `ffi` and `headers` features when building the crate.
9use geo::{prelude::*, Coord, Rect};
10use geo::{GeoFloat, Point, Polygon};
11use num_traits::FromPrimitive;
12use std::cmp::Ordering;
13use std::collections::BinaryHeap;
14use std::iter::Sum;
15use std::ops::{Deref, DerefMut};
16
17pub mod errors;
18use errors::PolylabelError;
19
20#[cfg(not(target_arch = "wasm32"))]
21#[cfg(feature = "ffi")]
22mod ffi;
23#[cfg(not(target_arch = "wasm32"))]
24#[cfg(feature = "ffi")]
25pub use crate::ffi::{polylabel_ffi, Array, Position, WrapperArray};
26
27/// Represention of a Quadtree node's cells. A node contains four Qcells.
28#[derive(Debug)]
29struct Qcell<T>
30where
31    T: GeoFloat,
32{
33    // The cell's centroid
34    centroid: Point<T>,
35    // Half of the parent node's extent
36    half_extent: T,
37    // Distance from centroid to polygon
38    distance: T,
39    // Maximum distance to polygon within a cell
40    max_distance: T,
41}
42
43impl<T> Qcell<T>
44where
45    T: GeoFloat,
46{
47    fn new(centroid: Point<T>, half_extent: T, polygon: &Polygon<T>) -> Qcell<T> {
48        let two = T::one() + T::one();
49        let distance = signed_distance(centroid, polygon);
50        let max_distance = distance + half_extent * two.sqrt();
51        Qcell {
52            centroid,
53            half_extent,
54            distance,
55            max_distance,
56        }
57    }
58}
59
60impl<T> Ord for Qcell<T>
61where
62    T: GeoFloat,
63{
64    fn cmp(&self, other: &Qcell<T>) -> std::cmp::Ordering {
65        self.max_distance.partial_cmp(&other.max_distance).unwrap()
66    }
67}
68impl<T> PartialOrd for Qcell<T>
69where
70    T: GeoFloat,
71{
72    fn partial_cmp(&self, other: &Qcell<T>) -> Option<Ordering> {
73        Some(self.cmp(other))
74    }
75}
76impl<T> Eq for Qcell<T> where T: GeoFloat {}
77impl<T> PartialEq for Qcell<T>
78where
79    T: GeoFloat,
80{
81    fn eq(&self, other: &Qcell<T>) -> bool
82    where
83        T: GeoFloat,
84    {
85        self.max_distance == other.max_distance
86    }
87}
88
89/// Signed distance from a Qcell's centroid to a Polygon's outline
90/// Returned value is negative if the point is outside the polygon's exterior ring
91fn signed_distance<T>(point: Point<T>, polygon: &Polygon<T>) -> T
92where
93    T: GeoFloat,
94{
95    let inside = polygon.contains(&point);
96    // Use LineString distance, because Polygon distance returns 0.0 for inside
97    let distance = point.euclidean_distance(polygon.exterior());
98    if inside {
99        distance
100    } else {
101        -distance
102    }
103}
104
105struct QuadTree<T>(pub BinaryHeap<Qcell<T>>)
106where
107    T: GeoFloat;
108
109impl<T> Deref for QuadTree<T>
110where
111    T: GeoFloat,
112{
113    type Target = BinaryHeap<Qcell<T>>;
114    fn deref(&self) -> &Self::Target {
115        &self.0
116    }
117}
118impl<T> DerefMut for QuadTree<T>
119where
120    T: GeoFloat,
121{
122    fn deref_mut(&mut self) -> &mut Self::Target {
123        &mut self.0
124    }
125}
126
127impl<T> QuadTree<T>
128where
129    T: GeoFloat,
130{
131    pub fn new(bbox: Rect<T>, half_extent: T, polygon: &Polygon<T>) -> Self {
132        let mut cell_queue: BinaryHeap<Qcell<T>> = BinaryHeap::new();
133
134        let two = T::one() + T::one();
135        let cell_size = half_extent * two;
136
137        let nx = (bbox.width() / cell_size).ceil().to_usize();
138        let ny = (bbox.height() / cell_size).ceil().to_usize();
139
140        match (nx, ny) {
141            (Some(nx), Some(ny)) => {
142                let one = T::one();
143                let delta_mid = Coord { x: one, y: one } * half_extent;
144                let origin = bbox.min();
145                let inital_points = (0..nx)
146                    .flat_map(|x| (0..ny).map(move |y| (x, y)))
147                    .filter_map(|(x, y)| Some((T::from(x)?, T::from(y)?)))
148                    .map(|(x, y)| Coord { x, y } * cell_size)
149                    .map(|delta_cell| origin + delta_cell + delta_mid)
150                    .map(Point::from)
151                    .map(|centroid| Qcell::new(centroid, half_extent, polygon));
152                cell_queue.extend(inital_points);
153            }
154            _ => {
155                // Do nothing, maybe error instead?
156            }
157        }
158
159        Self(cell_queue)
160    }
161
162    pub fn add_quad(&mut self, cell: &Qcell<T>, half_extent: T, polygon: &Polygon<T>) {
163        let new_cells = [
164            (-T::one(), -T::one()),
165            (T::one(), -T::one()),
166            (-T::one(), T::one()),
167            (T::one(), T::one()),
168        ]
169        .map(|(sign_x, sign_y)| (sign_x * half_extent, sign_y * half_extent))
170        .map(|(dx, dy)| Point::new(dx, dy))
171        .map(|delta| cell.centroid + delta)
172        .map(|centroid| Qcell::new(centroid, half_extent, polygon));
173        self.extend(new_cells);
174    }
175}
176
177/// Calculate a Polygon's ideal label position by calculating its ✨pole of inaccessibility✨
178///
179/// The calculation uses an [iterative grid-based algorithm](https://github.com/mapbox/polylabel#how-the-algorithm-works).
180///
181/// # Examples
182///
183/// ```
184/// use polylabel::polylabel;
185/// extern crate geo;
186/// use geo::{Point, LineString, Polygon};
187/// use geo::prelude::*;
188///
189/// // An approximate `L` shape
190/// let coords = vec![
191///    (0.0, 0.0),
192///    (4.0, 0.0),
193///    (4.0, 1.0),
194///    (1.0, 1.0),
195///    (1.0, 4.0),
196///    (0.0, 4.0),
197///    (0.0, 0.0)];
198///
199/// let poly = Polygon::new(coords.into(), vec![]);
200///
201/// // Its centroid lies outside the polygon
202/// assert_eq!(poly.centroid().unwrap(), Point::new(1.3571428571428572, 1.3571428571428572));
203///
204/// let label_position = polylabel(&poly, &0.1).unwrap();
205/// // Optimum label position is inside the polygon
206/// assert_eq!(label_position, Point::new(0.5625, 0.5625));
207/// ```
208///
209pub fn polylabel<T>(polygon: &Polygon<T>, tolerance: &T) -> Result<Point<T>, PolylabelError>
210where
211    T: GeoFloat + FromPrimitive + Sum,
212{
213    // special case for degenerate polygons
214    if polygon.signed_area() == T::zero() {
215        return Ok(Point::new(T::zero(), T::zero()));
216    }
217
218    let bbox = polygon
219        .bounding_rect()
220        .ok_or(PolylabelError::RectCalculation)?;
221    let cell_size = bbox.width().min(bbox.height());
222    // Special case for degenerate polygons
223    if cell_size == T::zero() {
224        return Ok(Point::from(bbox.min()));
225    }
226
227    let two = T::one() + T::one();
228    let mut half_extent = cell_size / two;
229
230    // initial best guess using centroid
231    let centroid = polygon
232        .centroid()
233        .ok_or(PolylabelError::CentroidCalculation)?;
234    let centroid_cell = Qcell::new(centroid, T::zero(), polygon);
235
236    // special case guess for rectangular polygons
237    let bbox_cell = Qcell::new(bbox.centroid(), T::zero(), polygon);
238
239    // deciding which initial guess was better
240    let mut best_cell = if bbox_cell.distance < centroid_cell.distance {
241        bbox_cell
242    } else {
243        centroid_cell
244    };
245
246    // setup priority queue
247    let mut cell_queue = QuadTree::<T>::new(bbox, half_extent, polygon);
248
249    // Now try to find better solutions
250    while let Some(cell) = cell_queue.pop() {
251        // Update the best cell if we find a cell with greater distance
252        if cell.distance > best_cell.distance {
253            best_cell = Qcell { ..cell };
254        }
255
256        // Bail out of this iteration if we can't find a better solution
257        if cell.max_distance - best_cell.distance <= *tolerance {
258            continue;
259        }
260
261        // Otherwise, add a new quadtree node and start again
262        half_extent = cell.half_extent / two;
263        cell_queue.add_quad(&cell, half_extent, polygon);
264    }
265
266    // We've exhausted the queue, so return the best solution we've found
267    Ok(best_cell.centroid)
268}
269
270#[cfg(test)]
271mod tests {
272    use super::{polylabel, Qcell};
273    use geo::prelude::*;
274    use geo::{Point, Polygon};
275    use std::collections::BinaryHeap;
276    #[test]
277    // polygons are those used in Shapely's tests
278    fn test_polylabel() {
279        let coords = include!("../tests/fixtures/poly1.rs");
280        let poly = Polygon::new(coords.into(), vec![]);
281        let res = polylabel(&poly, &10.000).unwrap();
282        assert_eq!(res, Point::new(59.35615556364569, 121.83919629746435));
283    }
284    #[test]
285    // does a concave polygon contain the calculated point?
286    // relevant because the centroid lies outside the polygon in this case
287    fn test_concave() {
288        let coords = include!("../tests/fixtures/poly2.rs");
289        let poly = Polygon::new(coords.into(), vec![]);
290        let res = polylabel(&poly, &1.0).unwrap();
291        assert!(poly.contains(&res));
292    }
293    #[test]
294    fn test_london() {
295        let coords = include!("../tests/fixtures/poly3.rs");
296        let poly = Polygon::new(coords.into(), vec![]);
297        let res = polylabel(&poly, &0.001).unwrap();
298        assert_eq!(res, Point::new(-0.45556816445920356, 51.54848888202887));
299    }
300    #[test]
301    fn polygon_l_test() {
302        // an L shape
303        let coords = vec![
304            (0.0, 0.0),
305            (4.0, 0.0),
306            (4.0, 1.0),
307            (1.0, 1.0),
308            (1.0, 4.0),
309            (0.0, 4.0),
310            (0.0, 0.0),
311        ];
312        let poly = Polygon::new(coords.into(), vec![]);
313        let res = polylabel(&poly, &0.10).unwrap();
314        assert_eq!(res, Point::new(0.5625, 0.5625));
315    }
316    #[test]
317    fn degenerate_polygon_test() {
318        let a_coords = vec![(0.0, 0.0), (1.0, 0.0), (2.0, 0.0), (0.0, 0.0)];
319        let a_poly = Polygon::new(a_coords.into(), vec![]);
320        let a_res = polylabel(&a_poly, &1.0).unwrap();
321        assert_eq!(a_res, Point::new(0.0, 0.0));
322    }
323    #[test]
324    fn degenerate_polygon_test_2() {
325        let b_coords = vec![(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)];
326        let b_poly = Polygon::new(b_coords.into(), vec![]);
327        let b_res = polylabel(&b_poly, &1.0).unwrap();
328        assert_eq!(b_res, Point::new(0.0, 0.0));
329    }
330    #[test]
331    // Is our priority queue behaving as it should?
332    fn test_queue() {
333        let a = Qcell {
334            centroid: Point::new(1.0, 2.0),
335            half_extent: 3.0,
336            distance: 4.0,
337            max_distance: 8.0,
338        };
339        let b = Qcell {
340            centroid: Point::new(1.0, 2.0),
341            half_extent: 3.0,
342            distance: 4.0,
343            max_distance: 7.0,
344        };
345        let c = Qcell {
346            centroid: Point::new(1.0, 2.0),
347            half_extent: 3.0,
348            distance: 4.0,
349            max_distance: 9.0,
350        };
351        let v = vec![a, b, c];
352        let mut q = BinaryHeap::from(v);
353        assert_eq!(q.pop().unwrap().max_distance, 9.0);
354        assert_eq!(q.pop().unwrap().max_distance, 8.0);
355        assert_eq!(q.pop().unwrap().max_distance, 7.0);
356    }
357}
polylabel/lib.rs

polylabel/
lib.rs
