gistools/data_structures/
point_cluster.rs

1use super::{IndexPoint, PointIndex, Tile};
2use crate::{
3    data_store::vector::{Vector, VectorStore},
4    geometry::{K_AVG_ANGLE_SPAN, LonLat, S1ChordAngle, S2CellId, S2Point, convert},
5    parsers::FeatureReader,
6};
7use alloc::{
8    collections::{BTreeMap, BTreeSet},
9    string::String,
10};
11use core::cell::RefCell;
12use libm::{floor, log2};
13use s2json::{
14    Face, GetXY, JSONCollection, MValue, PrimitiveValue, Projection, Properties, ValueType,
15    VectorFeature, VectorFeatureType, VectorGeometry, VectorGeometryType, VectorPoint,
16    VectorPointGeometry,
17};
18use serde::{Deserialize, Serialize, de::DeserializeOwned};
19
20/// Options for point clustering
21#[derive(Debug, Default, Clone, Serialize, Deserialize, PartialEq)]
22pub struct ClusterOptions {
23    /// projection to use
24    pub projection: Option<Projection>,
25    /// Name of the layer to build when requesting a tile
26    pub layer_name: Option<String>,
27    /// min zoom to generate clusters on
28    pub minzoom: Option<u8>,
29    /// max zoom level to cluster the points on
30    pub maxzoom: Option<u8>,
31    /// cluster radius in pixels. Assuming a tile size of 512x512, radius is a multiplier
32    /// of a singular pixel "cell"'s average length
33    pub radius: Option<f64>,
34}
35
36/// A cluster is a storage device to maintain groups of information in a cluster
37#[derive(Debug, Default, Clone, Serialize, Deserialize)]
38pub struct Cluster<M: Clone + Default = MValue> {
39    /// The data of the cluster
40    pub data: Option<M>,
41    /// The number of points in the cluster
42    pub count: usize,
43    /// Whether the cluster has been visited. 0 for not visited, 1 for visited
44    visited: u8,
45}
46impl<M: Clone + Default> Cluster<M> {
47    /// Create a cluster given the value (cluster size) and data
48    pub fn new(data: Option<M>, count: usize) -> Self {
49        Cluster { data, count, visited: 0 }
50    }
51}
52
53/// Used for storing features.
54#[derive(Debug, Default, Clone, Serialize, Deserialize, PartialEq, MValue)]
55pub struct ClusterData {
56    /// The number of points that were added into the cluster
57    pub count: usize,
58}
59
60/// Trait for creating a new cluster struct
61pub trait NewCluster: Clone + Default {
62    /// Create the struct adding the count
63    fn new_count(count: usize) -> Self;
64}
65impl NewCluster for ClusterData {
66    fn new_count(count: usize) -> Self {
67        ClusterData { count }
68    }
69}
70impl NewCluster for MValue {
71    fn new_count(count: usize) -> Self {
72        MValue::from([("count".into(), ValueType::Primitive(PrimitiveValue::U64(count as u64)))])
73    }
74}
75
76/// Compare two data items, return true to merge data
77pub type ClusterDataComparitor<M> = fn(a: &Option<M>, b: &Option<M>) -> bool;
78
79/// A cluster point
80pub type ClusterPoint<M> = IndexPoint<Cluster<M>>;
81
82/// A locally allocated Point Cluster Store
83pub type LocalClusterStore<M = MValue> = PointCluster<M, Vector<S2CellId, ClusterPoint<M>>>;
84
85/// # Point Cluster
86///
87/// ## Description
88/// A cluster store to index points at each zoom level
89#[derive(Debug)]
90pub struct PointCluster<
91    M: Clone + Default + Serialize + DeserializeOwned = Properties,
92    S: VectorStore<S2CellId, ClusterPoint<M>> = Vector<S2CellId, ClusterPoint<M>>,
93> {
94    /// The projection
95    pub projection: Projection,
96    /// The layer name
97    pub layer_name: String,
98    /// store which faces are active. 0 face could be entire WM projection
99    pub faces: BTreeSet<Face>,
100    /// min zoom to preserve detail on
101    pub minzoom: u8,
102    /// max zoom to preserve detail on
103    pub maxzoom: u8,
104    /// cluster radius
105    pub radius: f64,
106    /// grid size. Defaults to 512x512 pixel tile
107    pub grid_size: u32,
108    /// Index store for each zoom
109    pub indexes: BTreeMap<u8, RefCell<PointIndex<Cluster<M>, S>>>, // zoom => index
110}
111impl<M: Clone + Default + Serialize + DeserializeOwned, S: VectorStore<S2CellId, ClusterPoint<M>>>
112    PointCluster<M, S>
113{
114    /// Create a new point cluster
115    pub fn new(
116        options: Option<ClusterOptions>,
117        max_zoom_store: Option<PointIndex<Cluster<M>, S>>,
118    ) -> PointCluster<M, S> {
119        let options = options.unwrap_or_default();
120        let mut cluster_store = PointCluster {
121            projection: options.projection.unwrap_or(Projection::S2),
122            layer_name: options.layer_name.unwrap_or(String::from("default")),
123            faces: BTreeSet::new(),
124            minzoom: options.minzoom.unwrap_or(0),
125            maxzoom: options.maxzoom.unwrap_or(16),
126            radius: options.radius.unwrap_or(40.0),
127            grid_size: 512,
128            indexes: BTreeMap::new(),
129        };
130        // one extra zoom incase its a cell search system (bottom zoom isn't clustered to a cell)
131        for zoom in cluster_store.minzoom..=cluster_store.maxzoom + 1 {
132            cluster_store
133                .indexes
134                .insert(zoom, PointIndex::new(None, Some(cluster_store.projection)).into());
135        }
136        // inject a store if provided
137        if let Some(store) = max_zoom_store {
138            cluster_store.indexes.insert(cluster_store.maxzoom + 1, store.into());
139        }
140
141        cluster_store
142    }
143
144    /// Add a point to the maxzoom index. The point is a Point3D
145    pub fn insert(&mut self, id: S2CellId, point: S2Point, data: Option<M>) {
146        let mut maxzoom_index = self.indexes.get(&(self.maxzoom + 1)).unwrap().borrow_mut();
147        let cluster = Cluster::new(data, 1);
148        maxzoom_index.insert(id, point, Some(cluster));
149    }
150
151    /// Add a lon-lat pair as with any shape
152    pub fn insert_point<P: GetXY>(&mut self, point: P, data: Option<M>) {
153        self.insert_vector_point(VectorPoint::new(point.x(), point.y(), None, data));
154    }
155
156    /// Add a lon-lat pair to the cluster
157    pub fn insert_lon_lat(&mut self, mut ll: LonLat<M>) {
158        self.insert_vector_point(ll.take());
159    }
160
161    /// Add a lon-lat pair from a vectorPoint to the cluster
162    pub fn insert_vector_point(&mut self, ll: VectorPoint<M>) {
163        self.insert_feature(JSONCollection::<(), M, M>::VectorFeature(VectorFeature {
164            _type: VectorFeatureType::VectorFeature,
165            geometry: VectorGeometry::Point(VectorPointGeometry {
166                _type: VectorGeometryType::Point,
167                coordinates: ll,
168                is_3d: false,
169                ..Default::default()
170            }),
171            ..Default::default()
172        }));
173    }
174
175    /// Insert an STPoint to the index
176    pub fn insert_face_st(&mut self, face: Face, s: f64, t: f64, data: M) {
177        self.insert_feature(JSONCollection::<(), M, M>::VectorFeature(VectorFeature {
178            _type: VectorFeatureType::S2Feature,
179            face,
180            geometry: VectorGeometry::Point(VectorPointGeometry {
181                _type: VectorGeometryType::Point,
182                coordinates: VectorPoint::new(s, t, None, Some(data)),
183                is_3d: false,
184                ..Default::default()
185            }),
186            ..Default::default()
187        }));
188    }
189
190    /// Add all points from a reader. It will try to use the M-value first, but if it doesn't exist
191    /// it will use the feature properties data
192    pub fn insert_reader<T: Clone, F: FeatureReader<T, M, M>>(&mut self, reader: &F) {
193        for feature in reader.iter() {
194            self.insert_feature(JSONCollection::<T, M, M>::VectorFeature(feature))
195        }
196    }
197
198    /// Add a vector feature. It will try to use the M-value first, but if it doesn't exist
199    /// it will use the feature properties data
200    pub fn insert_feature<T: Clone>(&mut self, data: JSONCollection<T, M, M>) {
201        let features = convert(self.projection, &data, Some(true), Some(true));
202        for feature in features {
203            match feature.geometry {
204                VectorGeometry::Point(geometry) => {
205                    let coordinates = geometry.coordinates;
206                    self.faces.insert(feature.face);
207                    self._insert_face_st(
208                        feature.face.into(),
209                        coordinates.x,
210                        coordinates.y,
211                        coordinates.m.or(Some(feature.properties)),
212                    );
213                }
214                VectorGeometry::MultiPoint(geometry) => {
215                    for point in geometry.coordinates {
216                        self.faces.insert(feature.face);
217                        self._insert_face_st(
218                            feature.face.into(),
219                            point.x,
220                            point.y,
221                            point.m.or(Some(feature.properties.clone())),
222                        );
223                    }
224                }
225                _ => {}
226            }
227        }
228    }
229
230    /// Build the clusters when done adding points
231    pub fn build_clusters(&mut self, cmp_: Option<ClusterDataComparitor<M>>) {
232        let cmp = cmp_.unwrap_or(|_a: &Option<M>, _b: &Option<M>| true);
233        for zoom in (self.minzoom..=self.maxzoom).rev() {
234            self.cluster_radius(zoom, &cmp);
235        }
236        // ensure minzoom is sorted:
237        let mut minzoom_index = self.indexes.get(&self.minzoom).unwrap().borrow_mut();
238        minzoom_index.sort();
239    }
240
241    /// Radial clustering
242    fn cluster_radius(&mut self, zoom: u8, cmp: &ClusterDataComparitor<M>) {
243        let radius = self.get_level_radius(zoom);
244        let mut curr_index = self.indexes.get(&zoom).unwrap().borrow_mut();
245        let mut query_index = self.indexes.get(&(zoom + 1)).unwrap().borrow_mut();
246        query_index.sort();
247        let len = query_index.len();
248        for i in 0..len {
249            let (_, point) = query_index.get_index_mut(i).unwrap();
250            let cluster_data = &mut point.data;
251            if cluster_data.visited == 1 {
252                continue;
253            }
254            cluster_data.visited = 1;
255            let point = point.clone();
256            let cluster_data = &point.data;
257            // setup a new weighted cluster point
258            let mut new_num_points = cluster_data.count;
259            let mut new_cluster_point: S2Point = point.s2point;
260            new_cluster_point *= new_num_points as f64;
261            // joining all points found within radius
262            for (_, found_point) in query_index.s2_search_radius_mut(point.s2point, radius, None) {
263                let found_data = &mut found_point.data;
264                let count = found_data.count;
265                // only add points that match or have not been visited already
266                if !cmp(&cluster_data.data, &found_data.data) || found_data.visited == 1 {
267                    continue;
268                }
269                found_data.visited = 1;
270                // weighted add to new_cluster_point position
271                let fp = found_point.s2point;
272                new_cluster_point += fp * count as f64;
273                new_num_points += count;
274            }
275            // finish position average
276            new_cluster_point /= new_num_points as f64;
277            new_cluster_point.normalize();
278            // store the new cluster point
279            curr_index.insert(
280                (&new_cluster_point).into(),
281                new_cluster_point,
282                Some(Cluster::new(cluster_data.data.clone(), new_num_points)),
283            );
284        }
285    }
286
287    /// Given a tile ID, return the vector tile of all the clustered points
288    pub fn get_tile<C: NewCluster>(&mut self, id: S2CellId) -> Tile<(), M, C> {
289        let mut tile: Tile<(), M, C> = Tile::new(id);
290        let zoom = id.level();
291        if zoom < self.minzoom {
292            return tile;
293        }
294        let (min, max) = id.range();
295        // get the index
296        let index = self.indexes.get(&u8::min(zoom, self.maxzoom));
297        if index.is_none() {
298            return tile;
299        }
300
301        for (_, point) in index.unwrap().borrow().search_range(min, Some(max), None) {
302            let (face, s, t) = point.s2point.to_face_st();
303            let cluster = &point.data;
304            tile.add_feature(
305                VectorFeature {
306                    _type: VectorFeatureType::VectorFeature,
307                    face: face.into(),
308                    properties: cluster.data.clone().unwrap_or_default(),
309                    geometry: VectorGeometry::new_point(
310                        VectorPoint {
311                            x: s,
312                            y: t,
313                            z: None,
314                            m: Some(C::new_count(cluster.count)),
315                            t: None,
316                        },
317                        None,
318                    ),
319                    ..Default::default()
320                },
321                Some(self.layer_name.clone()),
322            );
323        }
324        // transform the geometry to be relative to the tile
325        tile.transform(0., Some(self.maxzoom));
326
327        tile
328    }
329
330    /// Insert an STPoint to the index
331    fn _insert_face_st(&mut self, face: u8, s: f64, t: f64, data: Option<M>) {
332        let point = S2Point::from_face_st(face, s, t);
333        self.insert((&point).into(), point, data);
334    }
335
336    /// Get a S1ChordAngle relative to a tile zoom level
337    fn get_level_radius(&self, zoom: u8) -> S1ChordAngle {
338        let zoom: i32 = zoom as i32;
339        let grid_size = self.grid_size as f64;
340
341        // if the grid is 512 x 512, log2 is 9, meaning the quadtree must split 9 times to analyze
342        // each individual pixel. Make sure we don't dive past 30 levels as that's the limit of the spec.
343        let zoom_grid_cell_level = i32::min(zoom + (floor(log2(grid_size)) as i32), 30);
344        // get the average angle span for all cells at the zoom's cell level
345        // divide by 2 to get half length to act more like a radius instead of a diameter
346        let angle_span = K_AVG_ANGLE_SPAN.get_value(zoom_grid_cell_level) / 2.;
347
348        S1ChordAngle::from_length2(angle_span * self.radius)
349    }
350}
gistools/data_structures/point_cluster.rs

gistools/data_structures/
point_cluster.rs
