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
17namespace mapbox {
18namespace supercluster {
19
20using namespace mapbox::geometry;
21
22struct Cluster {
23 point<double> pos;
24 std::uint32_t num_points;
25 std::uint32_t id = 0;
26 bool visited = false;
27
28 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
41namespace kdbush {
42
43using Cluster = mapbox::supercluster::Cluster;
44
45template <>
46struct nth<0, Cluster> {
47 inline static double get(const Cluster &c) {
48 return c.pos.x;
49 };
50};
51template <>
52struct nth<1, Cluster> {
53 inline static double get(const Cluster &c) {
54 return c.pos.y;
55 };
56};
57
58} // namespace kdbush
59
60namespace mapbox {
61namespace supercluster {
62
63#ifdef DEBUG_TIMER
64class Timer {
65public:
66 std::chrono::high_resolution_clock::time_point started;
67 Timer() {
68 started = std::chrono::high_resolution_clock::now();
69 }
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
79struct 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
86class 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
94public:
95 const GeoJSONFeatures features;
96 const Options options;
97
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(args: options.maxZoom + 1, args: 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(x: 2, y: z));
112 zooms.emplace(args&: z, args: Zoom(zooms[z + 1], r));
113#ifdef DEBUG_TIMER
114 timer(std::to_string(zooms[z].clusters.size()) + " clusters");
115#endif
116 }
117 }
118
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(x: 2, y: 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(x: this->options.extent * (c.pos.x * z2 - x)),
131 ::round(x: 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(x: feature);
142 };
143
144 double const top = (y - r) / z2;
145 double const bottom = (y + 1 + r) / z2;
146
147 zoom.tree.range(minX: (x - r) / z2, minY: top, maxX: (x + 1 + r) / z2, maxY: bottom, visitor);
148
149 if (x_ == 0) {
150 x = z2;
151 zoom.tree.range(minX: 1 - r / z2, minY: top, maxX: 1, maxY: bottom, visitor);
152 }
153 if (x_ == z2 - 1) {
154 x = -1;
155 zoom.tree.range(minX: 0, minY: top, maxX: r / z2, maxY: bottom, visitor);
156 }
157
158 return result;
159 }
160
161private:
162 struct Zoom {
163 kdbush::KDBush<Cluster, std::uint32_t> tree;
164 std::vector<Cluster> clusters;
165
166 Zoom() = default;
167
168 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(x: { project(p: f.geometry.get<GeoJSONPoint>()), 1, i++ });
174 }
175
176 tree.fill(points_: clusters);
177 }
178
179 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(qx: p.pos.x, qy: p.pos.y, r, visitor: [&](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(x: { weight / double(num_points), num_points, p.id });
203 }
204
205 tree.fill(points_: clusters);
206 }
207 };
208
209 std::unordered_map<std::uint8_t, Zoom> zooms;
210
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
219 static point<double> project(const GeoJSONPoint &p) {
220 auto lngX = p.x / 360 + 0.5;
221 const double sine = std::sin(x: p.y * M_PI / 180);
222 const double y = 0.5 - 0.25 * std::log(x: (1 + sine) / (1 - sine)) / M_PI;
223 auto latY = std::min(a: std::max(a: y, b: 0.0), b: 1.0);
224 return { lngX, latY };
225 }
226};
227
228} // namespace supercluster
229} // namespace mapbox
230

source code of qtlocation/src/3rdparty/mapbox-gl-native/deps/supercluster/0.2.2/include/supercluster.hpp