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

gistools/data_structures/
point_cluster.rs
