1 #pragma once
2 
3 #include <kdbush.hpp>
4 #include <mapbox/geometry/feature.hpp>
5 #include <mapbox/geometry/point_arithmetic.hpp>
6 
7 #include <algorithm>
8 #include <cmath>
9 #include <cstdint>
10 #include <vector>
11 
12 #ifdef DEBUG_TIMER
13 #include <chrono>
14 #include <iostream>
15 #endif
16 
17 namespace mapbox {
18 namespace supercluster {
19 
20 using namespace mapbox::geometry;
21 
22 struct Cluster {
23     point<double> pos;
24     std::uint32_t num_points;
25     std::uint32_t id = 0;
26     bool visited = false;
27 
Clustermapbox::supercluster::Cluster28     Cluster(point<double> pos_,
29             std::uint32_t num_points_,
30             std::uint32_t id_ = 0,
31             bool visited_ = false)
32         : pos(std::move(pos_)),
33           num_points(num_points_),
34           id(id_),
35           visited(visited_) {}
36 };
37 
38 } // namespace supercluster
39 } // namespace mapbox
40 
41 namespace kdbush {
42 
43 using Cluster = mapbox::supercluster::Cluster;
44 
45 template <>
46 struct nth<0, Cluster> {
getkdbush::nth47     inline static double get(const Cluster &c) {
48         return c.pos.x;
49     };
50 };
51 template <>
52 struct nth<1, Cluster> {
getkdbush::nth53     inline static double get(const Cluster &c) {
54         return c.pos.y;
55     };
56 };
57 
58 } // namespace kdbush
59 
60 namespace mapbox {
61 namespace supercluster {
62 
63 #ifdef DEBUG_TIMER
64 class Timer {
65 public:
66     std::chrono::high_resolution_clock::time_point started;
Timer()67     Timer() {
68         started = std::chrono::high_resolution_clock::now();
69     }
operator ()(std::string msg)70     void operator()(std::string msg) {
71         const auto now = std::chrono::high_resolution_clock::now();
72         const auto ms = std::chrono::duration_cast<std::chrono::microseconds>(now - started);
73         std::cerr << msg << ": " << double(ms.count()) / 1000 << "ms\n";
74         started = now;
75     }
76 };
77 #endif
78 
79 struct Options {
80     std::uint8_t minZoom = 0;   // min zoom to generate clusters on
81     std::uint8_t maxZoom = 16;  // max zoom level to cluster the points on
82     std::uint16_t radius = 40;  // cluster radius in pixels
83     std::uint16_t extent = 512; // tile extent (radius is calculated relative to it)
84 };
85 
86 class Supercluster {
87     using GeoJSONPoint = point<double>;
88     using GeoJSONFeatures = feature_collection<double>;
89 
90     using TilePoint = point<std::int16_t>;
91     using TileFeature = feature<std::int16_t>;
92     using TileFeatures = feature_collection<std::int16_t>;
93 
94 public:
95     const GeoJSONFeatures features;
96     const Options options;
97 
Supercluster(const GeoJSONFeatures & features_,const Options options_=Options ())98     Supercluster(const GeoJSONFeatures &features_, const Options options_ = Options())
99         : features(features_), options(options_) {
100 
101 #ifdef DEBUG_TIMER
102         Timer timer;
103 #endif
104         // convert and index initial points
105         zooms.emplace(options.maxZoom + 1, features);
106 #ifdef DEBUG_TIMER
107         timer(std::to_string(features.size()) + " initial points");
108 #endif
109         for (int z = options.maxZoom; z >= options.minZoom; z--) {
110             // cluster points from the previous zoom level
111             const double r = options.radius / (options.extent * std::pow(2, z));
112             zooms.emplace(z, Zoom(zooms[z + 1], r));
113 #ifdef DEBUG_TIMER
114             timer(std::to_string(zooms[z].clusters.size()) + " clusters");
115 #endif
116         }
117     }
118 
getTile(std::uint8_t z,std::uint32_t x_,std::uint32_t y)119     TileFeatures getTile(std::uint8_t z, std::uint32_t x_, std::uint32_t y) {
120         TileFeatures result;
121         auto &zoom = zooms[limitZoom(z)];
122 
123         std::uint32_t z2 = std::pow(2, z);
124         double const r = static_cast<double>(options.radius) / options.extent;
125         std::int32_t x = static_cast<std::int32_t>(x_);
126 
127         auto visitor = [&, this](const auto &id) {
128             auto const &c = zoom.clusters[id];
129 
130             TilePoint point(::round(this->options.extent * (c.pos.x * z2 - x)),
131                             ::round(this->options.extent * (c.pos.y * z2 - y)));
132             TileFeature feature{ point };
133 
134             if (c.num_points == 1) {
135                 feature.properties = this->features[c.id].properties;
136             } else {
137                 feature.properties["cluster"] = true;
138                 feature.properties["point_count"] = static_cast<std::uint64_t>(c.num_points);
139             }
140 
141             result.push_back(feature);
142         };
143 
144         double const top = (y - r) / z2;
145         double const bottom = (y + 1 + r) / z2;
146 
147         zoom.tree.range((x - r) / z2, top, (x + 1 + r) / z2, bottom, visitor);
148 
149         if (x_ == 0) {
150             x = z2;
151             zoom.tree.range(1 - r / z2, top, 1, bottom, visitor);
152         }
153         if (x_ == z2 - 1) {
154             x = -1;
155             zoom.tree.range(0, top, r / z2, bottom, visitor);
156         }
157 
158         return result;
159     }
160 
161 private:
162     struct Zoom {
163         kdbush::KDBush<Cluster, std::uint32_t> tree;
164         std::vector<Cluster> clusters;
165 
166         Zoom() = default;
167 
Zoommapbox::supercluster::Supercluster::Zoom168         Zoom(const GeoJSONFeatures &features_) {
169             // generate a cluster object for each point
170             std::uint32_t i = 0;
171 
172             for (const auto &f : features_) {
173                 clusters.push_back({ project(f.geometry.get<GeoJSONPoint>()), 1, i++ });
174             }
175 
176             tree.fill(clusters);
177         }
178 
Zoommapbox::supercluster::Supercluster::Zoom179         Zoom(Zoom &previous, double r) {
180             for (auto &p : previous.clusters) {
181                 if (p.visited)
182                     continue;
183                 p.visited = true;
184 
185                 auto num_points = p.num_points;
186                 point<double> weight = p.pos * double(num_points);
187 
188                 // find all nearby points
189                 previous.tree.within(p.pos.x, p.pos.y, r, [&](const auto &id) {
190                     auto &b = previous.clusters[id];
191 
192                     // filter out neighbors that are already processed
193                     if (b.visited)
194                         return;
195                     b.visited = true;
196 
197                     // accumulate coordinates for calculating weighted center
198                     weight += b.pos * double(b.num_points);
199                     num_points += b.num_points;
200                 });
201 
202                 clusters.push_back({ weight / double(num_points), num_points, p.id });
203             }
204 
205             tree.fill(clusters);
206         }
207     };
208 
209     std::unordered_map<std::uint8_t, Zoom> zooms;
210 
limitZoom(std::uint8_t z)211     std::uint8_t limitZoom(std::uint8_t z) {
212         if (z < options.minZoom)
213             return options.minZoom;
214         if (z > options.maxZoom + 1)
215             return options.maxZoom + 1;
216         return z;
217     }
218 
project(const GeoJSONPoint & p)219     static point<double> project(const GeoJSONPoint &p) {
220         auto lngX = p.x / 360 + 0.5;
221         const double sine = std::sin(p.y * M_PI / 180);
222         const double y = 0.5 - 0.25 * std::log((1 + sine) / (1 - sine)) / M_PI;
223         auto latY = std::min(std::max(y, 0.0), 1.0);
224         return { lngX, latY };
225     }
226 };
227 
228 } // namespace supercluster
229 } // namespace mapbox
230