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

gistools/data_structures/
point_cluster.rs
