blob: af5bbc312939ba79f6a077c742493149fb86fc15 [file] [log] [blame]
#pragma once
#include <mapbox/geojsonvt/convert.hpp>
#include <mapbox/geojsonvt/tile.hpp>
#include <mapbox/geojsonvt/types.hpp>
#include <mapbox/geojsonvt/wrap.hpp>
#include <chrono>
#include <cmath>
#include <map>
#include <unordered_map>
namespace mapbox {
namespace geojsonvt {
using geometry = mapbox::geometry::geometry<double>;
using feature = mapbox::geometry::feature<double>;
using feature_collection = mapbox::geometry::feature_collection<double>;
using geometry_collection = mapbox::geometry::geometry_collection<double>;
using geojson = mapbox::util::variant<geometry, feature, feature_collection>;
struct ToFeatureCollection {
feature_collection operator()(const feature_collection& value) const {
return value;
}
feature_collection operator()(const feature& value) const {
return { value };
}
feature_collection operator()(const geometry& value) const {
return { { value } };
}
};
struct TileOptions {
// simplification tolerance (higher means simpler)
double tolerance = 3;
// tile extent
uint16_t extent = 4096;
// tile buffer on each side
uint16_t buffer = 64;
};
struct Options : TileOptions {
// max zoom to preserve detail on
uint8_t maxZoom = 18;
// max zoom in the tile index
uint8_t indexMaxZoom = 5;
// max number of points per tile in the tile index
uint32_t indexMaxPoints = 100000;
};
const Tile empty_tile{};
inline uint64_t toID(uint8_t z, uint32_t x, uint32_t y) {
return (((1ull << z) * y + x) * 32) + z;
}
inline const Tile geoJSONToTile(const geojson& geojson_,
uint8_t z,
uint32_t x,
uint32_t y,
const TileOptions& options = TileOptions(),
bool wrap = false,
bool clip = false) {
const auto features_ = geojson::visit(geojson_, ToFeatureCollection{});
auto z2 = 1u << z;
auto tolerance = (options.tolerance / options.extent) / z2;
auto features = detail::convert(features_, tolerance);
if (wrap) {
features = detail::wrap(features, double(options.buffer) / options.extent);
}
if (clip) {
const double p = options.buffer / options.extent;
const auto left = detail::clip<0>(features, (x - p) / z2, (x + 1 + p) / z2, -1, 2);
features = detail::clip<1>(left, (y - p) / z2, (y + 1 + p) / z2, -1, 2);
}
return detail::InternalTile({ features, z, x, y, options.extent, tolerance }).tile;
}
class GeoJSONVT {
public:
const Options options;
GeoJSONVT(const mapbox::geometry::feature_collection<double>& features_,
const Options& options_ = Options())
: options(options_) {
const uint32_t z2 = 1u << options.maxZoom;
auto converted = detail::convert(features_, (options.tolerance / options.extent) / z2);
auto features = detail::wrap(converted, double(options.buffer) / options.extent);
splitTile(features, 0, 0, 0);
}
GeoJSONVT(const geojson& geojson_, const Options& options_ = Options())
: GeoJSONVT(geojson::visit(geojson_, ToFeatureCollection{}), options_) {
}
std::map<uint8_t, uint32_t> stats;
uint32_t total = 0;
const Tile& getTile(const uint8_t z, const uint32_t x_, const uint32_t y) {
if (z > options.maxZoom)
throw std::runtime_error("Requested zoom higher than maxZoom: " + std::to_string(z));
const uint32_t z2 = 1u << z;
const uint32_t x = ((x_ % z2) + z2) % z2; // wrap tile x coordinate
const uint64_t id = toID(z, x, y);
auto it = tiles.find(id);
if (it != tiles.end())
return it->second.tile;
it = findParent(z, x, y);
if (it == tiles.end())
throw std::runtime_error("Parent tile not found");
// if we found a parent tile containing the original geometry, we can drill down from it
const auto& parent = it->second;
// drill down parent tile up to the requested one
splitTile(parent.source_features, parent.z, parent.x, parent.y, z, x, y);
it = tiles.find(id);
if (it != tiles.end())
return it->second.tile;
it = findParent(z, x, y);
if (it == tiles.end())
throw std::runtime_error("Parent tile not found");
return empty_tile;
}
const std::unordered_map<uint64_t, detail::InternalTile>& getInternalTiles() const {
return tiles;
}
private:
std::unordered_map<uint64_t, detail::InternalTile> tiles;
std::unordered_map<uint64_t, detail::InternalTile>::iterator
findParent(const uint8_t z, const uint32_t x, const uint32_t y) {
uint8_t z0 = z;
uint32_t x0 = x;
uint32_t y0 = y;
const auto end = tiles.end();
auto parent = end;
while ((parent == end) && (z0 != 0)) {
z0--;
x0 = x0 / 2;
y0 = y0 / 2;
parent = tiles.find(toID(z0, x0, y0));
}
return parent;
}
void splitTile(const detail::vt_features& features,
const uint8_t z,
const uint32_t x,
const uint32_t y,
const uint8_t cz = 0,
const uint32_t cx = 0,
const uint32_t cy = 0) {
const double z2 = 1u << z;
const uint64_t id = toID(z, x, y);
auto it = tiles.find(id);
if (it == tiles.end()) {
const double tolerance =
(z == options.maxZoom ? 0 : options.tolerance / (z2 * options.extent));
it = tiles
.emplace(id,
detail::InternalTile{ features, z, x, y, options.extent, tolerance })
.first;
stats[z] = (stats.count(z) ? stats[z] + 1 : 1);
total++;
// printf("tile z%i-%i-%i\n", z, x, y);
}
auto& tile = it->second;
if (features.empty())
return;
// if it's the first-pass tiling
if (cz == 0u) {
// stop tiling if we reached max zoom, or if the tile is too simple
if (z == options.indexMaxZoom || tile.tile.num_points <= options.indexMaxPoints) {
tile.source_features = features;
return;
}
} else { // drilldown to a specific tile;
// stop tiling if we reached base zoom
if (z == options.maxZoom)
return;
// stop tiling if it's our target tile zoom
if (z == cz) {
tile.source_features = features;
return;
}
// stop tiling if it's not an ancestor of the target tile
const double m = 1u << (cz - z);
if (x != static_cast<uint32_t>(std::floor(cx / m)) ||
y != static_cast<uint32_t>(std::floor(cy / m))) {
tile.source_features = features;
return;
}
}
const double p = 0.5 * options.buffer / options.extent;
const auto& min = tile.bbox.min;
const auto& max = tile.bbox.max;
const auto left = detail::clip<0>(features, (x - p) / z2, (x + 0.5 + p) / z2, min.x, max.x);
splitTile(detail::clip<1>(left, (y - p) / z2, (y + 0.5 + p) / z2, min.y, max.y), z + 1,
x * 2, y * 2, cz, cx, cy);
splitTile(detail::clip<1>(left, (y + 0.5 - p) / z2, (y + 1 + p) / z2, min.y, max.y), z + 1,
x * 2, y * 2 + 1, cz, cx, cy);
const auto right =
detail::clip<0>(features, (x + 0.5 - p) / z2, (x + 1 + p) / z2, min.x, max.x);
splitTile(detail::clip<1>(right, (y - p) / z2, (y + 0.5 + p) / z2, min.y, max.y), z + 1,
x * 2 + 1, y * 2, cz, cx, cy);
splitTile(detail::clip<1>(right, (y + 0.5 - p) / z2, (y + 1 + p) / z2, min.y, max.y), z + 1,
x * 2 + 1, y * 2 + 1, cz, cx, cy);
// if we sliced further down, no need to keep source geometry
tile.source_features = {};
}
};
} // namespace geojsonvt
} // namespace mapbox