From ce774d2099f5c8bd971038b389d4558e13b606af Mon Sep 17 00:00:00 2001 From: dktue <2878623+dktue@users.noreply.github.com> Date: Mon, 12 Oct 2026 07:03:23 +0200 Subject: [PATCH 1/2] feat: add MapLibre Tile (MLT) output Add an optional MLT writer, built with -DTILEMAKER_WITH_MLT=ON and selected with "tile_format": "mlt". The tile writing code is templated over an MVT (vtzero) and an MLT backend; MLT tiles are encoded with the maplibre-tile-spec C++ encoder, fetched at a pinned commit. - PMTiles tile type 0x06, MBTiles format=mlt, .mlt directory output - feature order (ZOrder) preserved; mixed int/float columns kept numeric - fast encoder settings for compressed output - CMake build moves to C++20, with fixes for Boost 1.74 Asio and the bundled sqlite_modern_cpp; Dockerfile and Linux CI build with MLT Co-Authored-By: Claude Opus 5.5 (1M context) --- .github/workflows/ci.yml | 6 +- CMakeLists.txt | 28 ++++- Dockerfile | 8 +- docs/CONFIGURATION.md | 9 ++ docs/INSTALL.md | 6 ++ include/external/sqlite_modern_cpp.h | 8 +- include/mlt_writer.h | 121 +++++++++++++++++++++ include/output_object.h | 9 +- include/pmtiles.h | 3 + include/shared_data.h | 5 +- resources/config-schema.json | 1 + src/mlt_writer.cpp | 63 +++++++++++ src/output_object.cpp | 29 ++++- src/pmtiles.cpp | 2 +- src/shared_data.cpp | 16 ++- src/tile_worker.cpp | 151 ++++++++++++++++++--------- src/tilemaker.cpp | 5 + test/config_validator.test.cpp | 34 ++++++ 18 files changed, 439 insertions(+), 65 deletions(-) create mode 100644 include/mlt_writer.h create mode 100644 src/mlt_writer.cpp diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index cb04206a..9a41708a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -419,7 +419,7 @@ jobs: run: | mkdir build cd build - cmake -DTILEMAKER_BUILD_STATIC=ON -DCMAKE_BUILD_TYPE=Release -DVCPKG_TARGET_TRIPLET="${{ matrix.triplet }}" -DCMAKE_TOOLCHAIN_FILE="$VCPKG_INSTALLATION_ROOT/scripts/buildsystems/vcpkg.cmake" -DCMAKE_CXX_COMPILER=g++ .. + cmake -DTILEMAKER_BUILD_STATIC=ON -DTILEMAKER_WITH_MLT=ON -DCMAKE_BUILD_TYPE=Release -DVCPKG_TARGET_TRIPLET="${{ matrix.triplet }}" -DCMAKE_TOOLCHAIN_FILE="$VCPKG_INSTALLATION_ROOT/scripts/buildsystems/vcpkg.cmake" -DCMAKE_CXX_COMPILER=g++ .. cmake --build . strip tilemaker @@ -449,6 +449,9 @@ jobs: run_tilemaker "${{ env.AREA }}.mbtiles" "${tilemaker}" "${{ env.AREA }}.osm.pbf" --config=resources/config-openmaptiles.json --process=resources/process-openmaptiles.lua --output="${{ env.AREA }}.mbtiles" --verbose --store "${{ runner.temp }}/store" run_tilemaker "${{ env.AREA }}-repeat.mbtiles" "${tilemaker}" "${{ env.AREA }}.osm.pbf" --config=resources/config-openmaptiles.json --process=resources/process-openmaptiles.lua --output="${{ env.AREA }}-repeat.mbtiles" --verbose --store "${{ runner.temp }}/store-repeat" + jq '.settings.tile_format = "mlt"' resources/config-openmaptiles.json > "${{ runner.temp }}/config-openmaptiles-mlt.json" + run_tilemaker "${{ env.AREA }}-mlt.pmtiles" "${tilemaker}" "${{ env.AREA }}.osm.pbf" --config="${{ runner.temp }}/config-openmaptiles-mlt.json" --process=resources/process-openmaptiles.lua --output="${{ env.AREA }}-mlt.pmtiles" --verbose + - name: Upload generated tiles uses: actions/upload-artifact@v7 with: @@ -459,6 +462,7 @@ jobs: ${{ env.AREA }}.pmtiles ${{ env.AREA }}-repeat.mbtiles ${{ env.AREA }}-repeat.pmtiles + ${{ env.AREA }}-mlt.pmtiles - name: 'Upload compiled executable' uses: actions/upload-artifact@v7 diff --git a/CMakeLists.txt b/CMakeLists.txt index 46ffd341..3207fb67 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -10,6 +10,7 @@ string(STRIP "${tilemaker_version}" tilemaker_version) project(tilemaker VERSION ${tilemaker_version}) OPTION(TILEMAKER_BUILD_STATIC "Attempt to link dependencies static" OFF) +OPTION(TILEMAKER_WITH_MLT "Support writing MapLibre Tiles (MLT); needs CMake 3.25+" OFF) list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake") @@ -69,7 +70,25 @@ else() include_directories(${LUAJIT_INCLUDE_DIR}) endif() -set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD 20) +# tilemaker doesn't use coroutines, and Boost.Asio's C++20 coroutine support +# doesn't compile with some Boost versions (e.g. 1.74 in Debian bookworm) +add_compile_definitions(BOOST_ASIO_DISABLE_CO_AWAIT) + +if(TILEMAKER_WITH_MLT) + # The C++ MLT encoder from the MapLibre Tile spec repository. Only the submodules + # the encoder needs are fetched; set FETCHCONTENT_SOURCE_DIR_MLT to use a local checkout. + include(FetchContent) + set(MLT_WITH_JSON OFF CACHE BOOL "" FORCE) + set(MLT_WITH_TESTS OFF CACHE BOOL "" FORCE) + set(MLT_WITH_TOOLS OFF CACHE BOOL "" FORCE) + FetchContent_Declare(mlt + GIT_REPOSITORY https://github.com/maplibre/maplibre-tile-spec.git + GIT_TAG a9fa488f54ba5e7edf0151cf61f7b297d76dffe2 + GIT_SUBMODULES cpp/vendor/fsst cpp/vendor/earcut + SOURCE_SUBDIR cpp) + FetchContent_MakeAvailable(mlt) +endif() if(NOT TM_VERSION) set(TM_VERSION v${PROJECT_VERSION}) @@ -185,8 +204,15 @@ file(GLOB tilemaker_src_files src/visvalingam.cpp src/way_stores.cpp ) +if(TILEMAKER_WITH_MLT) + list(APPEND tilemaker_src_files src/mlt_writer.cpp) +endif() add_executable(tilemaker ${tilemaker_src_files}) target_include_directories(tilemaker PRIVATE include) +if(TILEMAKER_WITH_MLT) + target_compile_definitions(tilemaker PRIVATE TILEMAKER_WITH_MLT) + target_link_libraries(tilemaker mlt-cpp-encoder) +endif() target_include_directories(tilemaker PRIVATE ${CMAKE_BINARY_DIR}) # for generated files target_link_libraries(tilemaker ${THREAD_LIB} ${CMAKE_DL_LIBS} diff --git a/Dockerfile b/Dockerfile index 93850bc4..6fb8d748 100644 --- a/Dockerfile +++ b/Dockerfile @@ -15,7 +15,9 @@ RUN apt-get update && DEBIAN_FRONTEND=noninteractive apt-get install -y --no-ins libboost-system-dev \ luarocks \ rapidjson-dev \ - cmake && \ + cmake \ + git \ + ca-certificates && \ rm -rf /var/lib/apt/lists/* && \ luarocks install luaflock @@ -36,9 +38,9 @@ RUN mkdir build && \ tm_version_arg="-DTM_VERSION=$TILEMAKER_BUILD_VERSION"; \ fi && \ if [ -z "$BUILD_DEBUG" ]; then \ - cmake -DCMAKE_BUILD_TYPE=Release -DBoost_USE_DEBUG_RUNTIME=OFF $tm_version_arg ..; \ + cmake -DCMAKE_BUILD_TYPE=Release -DBoost_USE_DEBUG_RUNTIME=OFF -DTILEMAKER_WITH_MLT=ON $tm_version_arg ..; \ else \ - cmake -DCMAKE_BUILD_TYPE=Debug $tm_version_arg ..; \ + cmake -DCMAKE_BUILD_TYPE=Debug -DTILEMAKER_WITH_MLT=ON $tm_version_arg ..; \ fi; \ cmake --build . --parallel $(nproc) && \ if [ -z "$BUILD_DEBUG" ]; then \ diff --git a/docs/CONFIGURATION.md b/docs/CONFIGURATION.md index 82fad41b..6c191059 100644 --- a/docs/CONFIGURATION.md +++ b/docs/CONFIGURATION.md @@ -38,6 +38,7 @@ It also includes these global settings: * `bounding_box` (optional) - the bounding box to output, in [minlon, minlat, maxlon, maxlat] order * `default_view` (optional) - the default location for the client to view, in [lon, lat, zoom] order (MBTiles only) * `mvt_version` (optional) - the version of the [Mapbox Vector Tile](https://github.com/mapbox/vector-tile-spec) spec to use; defaults to 2 +* `tile_format` (optional) - `"mvt"` (default) for Mapbox Vector Tiles, or `"mlt"` for [MapLibre Tiles](https://maplibre.org/maplibre-tile-spec/). MLT output needs tilemaker to be built with MLT support (see below and INSTALL.md) * `combine_below` (deprecated) - whether to merge all linestrings in the tile with the same attributes: will be done at zoom levels below that specified here (e.g. `"combine_below": 14` to merge at z1-13). Please use the layer-specific parameter `combine_lines_below` instead. A typical config file would look like this: @@ -68,6 +69,14 @@ All options are compulsory unless stated otherwise. If tilemaker baulks at the J By default tilemaker expects to find this file at config.json, but you can specify another filename with the `--config` command-line option. +### MapLibre Tiles (MLT) + +With `"tile_format": "mlt"`, tilemaker writes [MapLibre Tiles](https://maplibre.org/maplibre-tile-spec/) instead of Mapbox Vector Tiles. Your Lua profile and layer configuration stay the same. MLT tiles are typically smaller than MVT tiles, and are supported by MapLibre GL JS and MapLibre Native. In your style, add `"encoding": "mlt"` to the vector source. + +* .pmtiles output is marked with the MapLibre Tile type, .mbtiles output has `format` set to `mlt`, and tiles written to a directory get the `.mlt` extension. +* MLT stores each attribute as a typed column per layer. If the same attribute is written as an integer for some features and as a number with decimals for others in one layer, the integers are stored as doubles. Other mixed types (e.g. a number for some features and a string for others) are stored as strings. +* `--merge` isn't supported with MLT output. + ### Advanced layer configuration You can add optional parameters to layers: diff --git a/docs/INSTALL.md b/docs/INSTALL.md index d19f3bcb..bc021829 100644 --- a/docs/INSTALL.md +++ b/docs/INSTALL.md @@ -52,6 +52,12 @@ You can optionally use cmake to build: make sudo make install +To be able to write [MapLibre Tiles](https://maplibre.org/maplibre-tile-spec/) (MLT), configure with `-DTILEMAKER_WITH_MLT=ON`. This needs cmake 3.25 or later and a C++20 compiler, and downloads the MLT encoder from GitHub during configuration (set `FETCHCONTENT_SOURCE_DIR_MLT` to use a local checkout of maplibre-tile-spec instead): + + cmake -DTILEMAKER_WITH_MLT=ON .. + +The Makefile build doesn't include MLT support. + ### Docker **The Dockerfile is not formally supported by project maintainers and you are encouraged to send pull requests to fix any issues you encounter.** diff --git a/include/external/sqlite_modern_cpp.h b/include/external/sqlite_modern_cpp.h index 40fd24d4..20c49798 100644 --- a/include/external/sqlite_modern_cpp.h +++ b/include/external/sqlite_modern_cpp.h @@ -259,7 +259,13 @@ namespace sqlite { ~database_binder() noexcept(false) { /* Will be executed if no >>op is found, but not if an exception is in mid flight */ - if(!execution_started && !std::uncaught_exception() && _stmt) { +#if __cplusplus >= 201703L || (defined(_MSVC_LANG) && _MSVC_LANG >= 201703L) + // std::uncaught_exception() was removed in C++20 + const bool unwinding = std::uncaught_exceptions() > 0; +#else + const bool unwinding = std::uncaught_exception(); +#endif + if(!execution_started && !unwinding && _stmt) { execute(); } } diff --git a/include/mlt_writer.h b/include/mlt_writer.h new file mode 100644 index 00000000..bc02a589 --- /dev/null +++ b/include/mlt_writer.h @@ -0,0 +1,121 @@ +/*! \file */ +#ifndef _MLT_WRITER_H +#define _MLT_WRITER_H + +#include +#include +#include +#include + +// Builders that collect the features of one tile for the MapLibre Tile (MLT) encoder. +// They mirror the subset of the vtzero builder API that tile_worker.cpp uses, so the +// same templated code writes either MVT or MLT. + +class MltLayerBuilder { +public: + MltLayerBuilder(const std::string& name, uint32_t extent) { + layer.name = name; + layer.extent = extent; + } + + mlt::Encoder::Layer layer; +}; + +class MltFeatureBuilder { +public: + explicit MltFeatureBuilder(MltLayerBuilder& layerBuilder) : layerBuilder(layerBuilder) { + feature.id = std::nullopt; + } + + void set_id(uint64_t id) { feature.id = id; } + + void add_property(const std::string& key, const std::string& value) { feature.properties.emplace(key, value); } + void add_property(const std::string& key, bool value) { feature.properties.emplace(key, value); } + void add_property(const std::string& key, int value) { feature.properties.emplace(key, static_cast(value)); } + void add_property(const std::string& key, float value) { feature.properties.emplace(key, value); } + +protected: + void commitFeature() { layerBuilder.layer.features.push_back(std::move(feature)); } + + MltLayerBuilder& layerBuilder; + mlt::Encoder::Feature feature; +}; + +class MltPointFeatureBuilder : public MltFeatureBuilder { +public: + using MltFeatureBuilder::MltFeatureBuilder; + + void add_points(uint32_t count) { feature.geometry.coordinates.reserve(count); } + void set_point(int x, int y) { feature.geometry.coordinates.push_back({x, y}); } + + void commit() { + auto& g = feature.geometry; + g.type = g.coordinates.size() == 1 ? mlt::Encoder::GeometryType::POINT : mlt::Encoder::GeometryType::MULTIPOINT; + commitFeature(); + } +}; + +class MltLinestringFeatureBuilder : public MltFeatureBuilder { +public: + using MltFeatureBuilder::MltFeatureBuilder; + + void add_linestring(uint32_t count) { + feature.geometry.parts.emplace_back(); + feature.geometry.parts.back().reserve(count); + } + void set_point(int x, int y) { feature.geometry.parts.back().push_back({x, y}); } + + void commit() { + auto& g = feature.geometry; + if (g.parts.size() == 1) { + g.type = mlt::Encoder::GeometryType::LINESTRING; + g.coordinates = std::move(g.parts.front()); + g.parts.clear(); + } else { + g.type = mlt::Encoder::GeometryType::MULTILINESTRING; + } + commitFeature(); + } +}; + +class MltPolygonFeatureBuilder : public MltFeatureBuilder { +public: + using MltFeatureBuilder::MltFeatureBuilder; + + // count includes the closing point, as with vtzero; MLT rings are stored without it + void add_ring(uint32_t count, bool outer) { + auto& g = feature.geometry; + if (outer) { + g.parts.emplace_back(); + g.partRingSizes.emplace_back(); + } + g.partRingSizes.back().push_back(count - 1); + ringPointsLeft = count; + } + void set_point(int x, int y) { + if (--ringPointsLeft > 0) feature.geometry.parts.back().push_back({x, y}); + } + + void commit() { + auto& g = feature.geometry; + if (g.parts.size() == 1) { + g.type = mlt::Encoder::GeometryType::POLYGON; + g.coordinates = std::move(g.parts.front()); + g.ringSizes = std::move(g.partRingSizes.front()); + g.parts.clear(); + g.partRingSizes.clear(); + } else { + g.type = mlt::Encoder::GeometryType::MULTIPOLYGON; + } + commitFeature(); + } + +private: + uint32_t ringPointsLeft = 0; +}; + +/// Encode the collected layers as an MLT tile. Empty layers are skipped. +/// Set compressed when the tile will be gzip/deflate compressed afterwards. +void encodeMltTile(std::vector& layers, bool includeIds, bool compressed, std::string& output); + +#endif //_MLT_WRITER_H diff --git a/include/output_object.h b/include/output_object.h index 4598d2d0..e0f4ac22 100644 --- a/include/output_object.h +++ b/include/output_object.h @@ -74,7 +74,14 @@ class OutputObject { vtzero::feature_builder& fbuilder, char zoom ) const; - +#ifdef TILEMAKER_WITH_MLT + void writeAttributes( + const AttributeStore& attributeStore, + class MltFeatureBuilder& fbuilder, + char zoom + ) const; +#endif + bool compatible(const OutputObject &other); }; diff --git a/include/pmtiles.h b/include/pmtiles.h index ea9c2491..7e8040c8 100644 --- a/include/pmtiles.h +++ b/include/pmtiles.h @@ -25,6 +25,8 @@ struct TileOffset { #define TINY_LENGTH 100 // Expire the tiny cache when it reaches this size #define TINY_MAX_SIZE 10000 +// PMTiles v3 tile type for MapLibre Tiles (not yet in external/pmtiles.hpp) +#define TILETYPE_MLT 0x6 class PMTiles { @@ -34,6 +36,7 @@ class PMTiles { pmtiles::headerv3 header; bool isSparse = true; + uint8_t tileType = pmtiles::TILETYPE_MVT; void open(std::string &filename); void saveTile(int zoom, int x, int y, std::string &data); diff --git a/include/shared_data.h b/include/shared_data.h index 5cfe13af..d4b1b8fb 100644 --- a/include/shared_data.h +++ b/include/shared_data.h @@ -77,12 +77,15 @@ class LayerDefinition { std::string serialiseToJSON() const; }; +enum class TileFormat { MVT, MLT }; + ///\brief Config read from JSON to control behavior of program class Config { - + public: class LayerDefinition layers; uint baseZoom, startZoom, endZoom; + TileFormat tileFormat; uint mvtVersion, combineBelow; bool includeID, compress, gzip, highResolution; std::string compressOpt; diff --git a/resources/config-schema.json b/resources/config-schema.json index ca8841c7..df45a2c3 100644 --- a/resources/config-schema.json +++ b/resources/config-schema.json @@ -17,6 +17,7 @@ "high_resolution": { "type": "boolean" }, "combine_below": { "type": "integer", "minimum": 0 }, "mvt_version": { "type": "integer", "minimum": 1 }, + "tile_format": { "enum": ["mvt", "mlt"] }, "bounding_box": { "type": "array", "minItems": 4, diff --git a/src/mlt_writer.cpp b/src/mlt_writer.cpp new file mode 100644 index 00000000..03cd2681 --- /dev/null +++ b/src/mlt_writer.cpp @@ -0,0 +1,63 @@ +/*! \file */ +#include "mlt_writer.h" +#include +#include +#include + +namespace { + +// MLT stores each property as a typed column per layer, and the encoder falls back to a +// string column when a key holds both integers and floats. In MVT every value carries its +// own type, so profiles legitimately mix them (e.g. AttributeInteger and AttributeNumeric +// for the same key). Promote the integers to double so such columns stay numeric. +void widenMixedNumericColumns(mlt::Encoder::Layer& layer) { + constexpr uint8_t SEEN_INT = 1, SEEN_FLOAT = 2; + std::unordered_map seen; + for (const auto& feature : layer.features) { + for (const auto& [key, value] : feature.properties) { + if (std::holds_alternative(value)) seen[key] |= SEEN_INT; + else if (std::holds_alternative(value)) seen[key] |= SEEN_FLOAT; + } + } + + std::unordered_set mixed; + for (const auto& [key, kinds] : seen) { + if (kinds == (SEEN_INT | SEEN_FLOAT)) mixed.emplace(key); + } + if (mixed.empty()) return; + + for (auto& feature : layer.features) { + for (const auto& key : mixed) { + auto it = feature.properties.find(key); + if (it != feature.properties.end() && std::holds_alternative(it->second)) { + it->second = static_cast(std::get(it->second)); + } + } + } +} + +} // namespace + +void encodeMltTile(std::vector& layers, bool includeIds, bool compressed, std::string& output) { + thread_local const mlt::Encoder encoder; + + const auto config = mlt::EncoderConfig::with([&](mlt::EncoderConfig& c) { + c.includeIds = includeIds; + // Feature order within a layer carries the profile's ZOrder, so it must be preserved + c.sortFeatures = false; + if (compressed) { + // The encoder's size-driven choices (FSST strings, Morton/Hilbert vertex dictionaries, + // per-stream integer encodings for geometry) make up most of its time, but once the tile + // is gzip/deflate compressed they save little or nothing. Encoding geometry with plain + // varints skips the dictionary candidates entirely. + c.useFsst = false; + c.enableMortonEncoding = false; + c.geometryEncodingOption = mlt::IntegerEncodingOption::PLAIN; + } + }); + + for (auto& layer : layers) widenMixedNumericColumns(layer); + + std::vector encoded = encoder.encode(layers, config); + output.assign(reinterpret_cast(encoded.data()), encoded.size()); +} diff --git a/src/output_object.cpp b/src/output_object.cpp index 93eef495..e94168b0 100644 --- a/src/output_object.cpp +++ b/src/output_object.cpp @@ -6,6 +6,9 @@ #include "helpers.h" #include "coordinates_geom.h" #include +#ifdef TILEMAKER_WITH_MLT +#include "mlt_writer.h" +#endif using namespace std; namespace geom = boost::geometry; @@ -31,11 +34,13 @@ std::ostream& operator<<(std::ostream& os, OutputGeometryType geomType) return os; } -void OutputObject::writeAttributes( +template +static void writeAttributesTo( const AttributeStore& attributeStore, - vtzero::feature_builder& fbuilder, + AttributeIndex attributes, + FeatureBuilder& fbuilder, char zoom -) const { +) { auto attr = attributeStore.getUnsafe(attributes); for(auto const &it: attr) { @@ -58,6 +63,24 @@ void OutputObject::writeAttributes( } } +void OutputObject::writeAttributes( + const AttributeStore& attributeStore, + vtzero::feature_builder& fbuilder, + char zoom +) const { + writeAttributesTo(attributeStore, attributes, fbuilder, zoom); +} + +#ifdef TILEMAKER_WITH_MLT +void OutputObject::writeAttributes( + const AttributeStore& attributeStore, + MltFeatureBuilder& fbuilder, + char zoom +) const { + writeAttributesTo(attributeStore, attributes, fbuilder, zoom); +} +#endif + bool OutputObject::compatible(const OutputObject &other) { return geomType == other.geomType && z_order == other.z_order && diff --git a/src/pmtiles.cpp b/src/pmtiles.cpp index 8cafa4be..c9561b59 100644 --- a/src/pmtiles.cpp +++ b/src/pmtiles.cpp @@ -75,7 +75,7 @@ void PMTiles::close(std::string &metadata) { header.clustered = false; header.internal_compression = pmtiles::COMPRESSION_GZIP; header.tile_compression = pmtiles::COMPRESSION_GZIP; - header.tile_type = pmtiles::TILETYPE_MVT; + header.tile_type = tileType; // write header std::string headerString = header.serialize(); diff --git a/src/shared_data.cpp b/src/shared_data.cpp index d3843562..d15a0bfd 100644 --- a/src/shared_data.cpp +++ b/src/shared_data.cpp @@ -22,7 +22,7 @@ void SharedData::writeMBTilesProjectData() { mbtiles.writeMetadata("type", "baselayer"); mbtiles.writeMetadata("version", config.projectVersion); mbtiles.writeMetadata("description", config.projectDesc); - mbtiles.writeMetadata("format", "pbf"); + mbtiles.writeMetadata("format", config.tileFormat == TileFormat::MLT ? "mlt" : "pbf"); mbtiles.writeMetadata("minzoom", to_string(config.startZoom)); mbtiles.writeMetadata("maxzoom", to_string(config.endZoom)); @@ -230,6 +230,7 @@ Config::Config() { clippingBoxFromJSON = false; baseZoom = 0; combineBelow = 0; + tileFormat = TileFormat::MVT; } Config::~Config() { } @@ -266,6 +267,19 @@ void Config::readConfig(rapidjson::Document &jsonConfig, bool &hasClippingBox, B compressOpt = jsonConfig["settings"]["compress"].GetString(); combineBelow = jsonConfig["settings"].HasMember("combine_below") ? jsonConfig["settings"]["combine_below"].GetUint() : 0; mvtVersion = jsonConfig["settings"].HasMember("mvt_version") ? jsonConfig["settings"]["mvt_version"].GetUint() : 2; + std::string tileFormatOpt = jsonConfig["settings"].HasMember("tile_format") ? jsonConfig["settings"]["tile_format"].GetString() : "mvt"; + if (tileFormatOpt == "mvt") { tileFormat = TileFormat::MVT; } + else if (tileFormatOpt == "mlt") { +#ifdef TILEMAKER_WITH_MLT + tileFormat = TileFormat::MLT; +#else + cerr << "\"tile_format\": \"mlt\" requires tilemaker to be built with MLT support (cmake -DTILEMAKER_WITH_MLT=ON)." << endl; + exit (EXIT_FAILURE); +#endif + } else { + cerr << "\"tile_format\" should be any of \"mvt\",\"mlt\" in JSON file." << endl; + exit (EXIT_FAILURE); + } projectName = jsonConfig["settings"]["name"].GetString(); projectVersion = jsonConfig["settings"]["version"].GetString(); projectDesc = jsonConfig["settings"]["description"].GetString(); diff --git a/src/tile_worker.cpp b/src/tile_worker.cpp index 30d09ebe..891dd00f 100644 --- a/src/tile_worker.cpp +++ b/src/tile_worker.cpp @@ -7,9 +7,32 @@ #include "helpers.h" #include "visvalingam.h" #include "simplify_buildings.h" +#ifdef TILEMAKER_WITH_MLT +#include "mlt_writer.h" +#endif using namespace std; extern bool verbose; +// The tile is written through one of these backends. Both expose the same builder API, +// so the geometry and attribute code below is shared between the output formats. +struct MvtBackend { + using LayerBuilder = vtzero::layer_builder; + using PointFeature = vtzero::point_feature_builder; + using LinestringFeature = vtzero::linestring_feature_builder; + using PolygonFeature = vtzero::polygon_feature_builder; + static void addRing(PolygonFeature& fbuilder, uint32_t points, bool /*outer*/) { fbuilder.add_ring(points); } +}; + +#ifdef TILEMAKER_WITH_MLT +struct MltBackend { + using LayerBuilder = MltLayerBuilder; + using PointFeature = MltPointFeatureBuilder; + using LinestringFeature = MltLinestringFeatureBuilder; + using PolygonFeature = MltPolygonFeatureBuilder; + static void addRing(PolygonFeature& fbuilder, uint32_t points, bool outer) { fbuilder.add_ring(points, outer); } +}; +#endif + thread_local bool enabledUserSignal = false; thread_local MultiPolygon scaledMultiPolygon; typedef std::vector::const_iterator OutputObjectsConstIt; @@ -95,10 +118,11 @@ void RemovePartsBelowSize(MultiPolygon &g, double filterArea) { } } +template void writeMultiLinestring( const AttributeStore& attributeStore, const SharedData& sharedData, - vtzero::layer_builder& vtLayer, + typename Backend::LayerBuilder& vtLayer, const TileBbox& bbox, const OutputObjectID& oo, unsigned zoom, @@ -106,7 +130,7 @@ void writeMultiLinestring( unsigned simplifyAlgo, const MultiLinestring& mls ) { - vtzero::linestring_feature_builder fbuilder{vtLayer}; + typename Backend::LinestringFeature fbuilder{vtLayer}; if (sharedData.config.includeID && oo.id) fbuilder.set_id(oo.id); @@ -174,9 +198,11 @@ void writeMultiLinestring( } } +template bool writeRing( - vtzero::polygon_feature_builder& fbuilder, - const Ring& ring + typename Backend::PolygonFeature& fbuilder, + const Ring& ring, + bool outer ) { // vtzero doesn't like zero-length segments in rings pair lastXy = std::make_pair(0, 0); @@ -194,7 +220,7 @@ bool writeRing( return false; bool firstPoint = true; - fbuilder.add_ring(points); + Backend::addRing(fbuilder, points, outer); for (auto it = ring.rbegin(); it != ring.rend(); ++it) { const Point& point = *it; @@ -212,10 +238,11 @@ bool writeRing( } +template void writeMultiPolygon( const AttributeStore& attributeStore, const SharedData& sharedData, - vtzero::layer_builder& vtLayer, + typename Backend::LayerBuilder& vtLayer, const TileBbox& bbox, const OutputObjectID& oo, unsigned zoom, @@ -271,7 +298,7 @@ void writeMultiPolygon( } } - vtzero::polygon_feature_builder fbuilder{vtLayer}; + typename Backend::PolygonFeature fbuilder{vtLayer}; if (sharedData.config.includeID && oo.id) fbuilder.set_id(oo.id); @@ -282,14 +309,14 @@ void writeMultiPolygon( const Ring& ring = geom::exterior_ring(p); // If we failed to write the outer, no need to write the inners. - if (!writeRing(fbuilder, ring)) + if (!writeRing(fbuilder, ring, true)) continue; hadPoly = true; const InteriorRing& interiors = geom::interior_rings(p); for (const Ring& ring : interiors) - writeRing(fbuilder, ring); + writeRing(fbuilder, ring, false); } if (hadPoly) { @@ -300,6 +327,7 @@ void writeMultiPolygon( } } +template void ProcessObjects( TileDataSource* source, const AttributeStore& attributeStore, @@ -314,7 +342,7 @@ void ProcessObjects( bool combinePolygons, unsigned zoom, const TileBbox &bbox, - vtzero::layer_builder& vtLayer + typename Backend::LayerBuilder& vtLayer ) { for (auto jt = ooSameLayerBegin; jt != ooSameLayerEnd; ++jt) { OutputObjectID oo = *jt; @@ -353,7 +381,7 @@ void ProcessObjects( continue; } - vtzero::point_feature_builder fbuilder{vtLayer}; + typename Backend::PointFeature fbuilder{vtLayer}; if (sharedData.config.includeID && oo.id) fbuilder.set_id(oo.id); fbuilder.add_points(multipoint.size()); @@ -410,9 +438,9 @@ void ProcessObjects( } if (oo.oo.geomType == LINESTRING_ || oo.oo.geomType == MULTILINESTRING_) - writeMultiLinestring(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, simplifyAlgo, boost::get(g)); + writeMultiLinestring(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, simplifyAlgo, boost::get(g)); else if (oo.oo.geomType == POLYGON_) - writeMultiPolygon(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, simplifyAlgo, boost::get(g)); + writeMultiPolygon(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, simplifyAlgo, boost::get(g)); } } } @@ -433,23 +461,12 @@ OutputObjectsConstItPair getObjectsAtSubLayer( return equal_range(data.begin(), data.end(), layerNum, layerComp()); } - - -void ProcessLayer( - const SourceList& sources, - const AttributeStore& attributeStore, - TileCoordinates index, - uint zoom, - const std::vector>& data, - vtzero::vector_tile existingTile, - vtzero::tile_builder& tile, - const TileBbox& bbox, - const std::vector& ltx, - SharedData& sharedData +// Copy a layer of a tile we're merging into (MVT only) +void CopyExistingLayer( + vtzero::vector_tile& existingTile, + const std::string& layerName, + vtzero::layer_builder& vtLayer ) { - std::string layerName = sharedData.layers.layers[ltx.at(0)].name; - vtzero::layer_builder vtLayer{tile, layerName, sharedData.config.mvtVersion, bbox.hires ? 8192u : 4096u}; - vtzero::layer existingLayer = existingTile.get_layer_by_name(layerName); if (existingLayer) { while (auto feature = existingLayer.next_feature()) { @@ -463,6 +480,21 @@ void ProcessLayer( fb.commit(); } } +} + +template +void ProcessLayer( + const SourceList& sources, + const AttributeStore& attributeStore, + TileCoordinates index, + uint zoom, + const std::vector>& data, + typename Backend::LayerBuilder& vtLayer, + const TileBbox& bbox, + const std::vector& ltx, + SharedData& sharedData +) { + const std::string& layerName = sharedData.layers.layers[ltx.at(0)].name; //TileCoordinate tileX = index.x; TileCoordinate tileY = index.y; @@ -493,7 +525,7 @@ void ProcessLayer( auto ooListSameLayer = getObjectsAtSubLayer(data[i], layerNum); auto end = ooListSameLayer.second; if (ld.featureLimit>0 && end-ooListSameLayer.first>ld.featureLimit && zoom(sources[i], attributeStore, ooListSameLayer.first, end, sharedData, simplifyLevel, ld.simplifyAlgo, filterArea, ld.combinePoints, zoom < ld.combineLinesBelow, zoom < ld.combinePolygonsBelow, zoom, bbox, vtLayer); @@ -518,9 +550,6 @@ void outputProc( TileCoordinates coordinates, uint zoom ) { - // Create tile - vtzero::tile_builder tile; - TileBbox bbox(coordinates, zoom, sharedData.config.highResolution && zoom==sharedData.config.endZoom, zoom==sharedData.config.endZoom); if (sharedData.config.clippingBoxFromJSON && ( sharedData.config.maxLon <= bbox.minLon || @@ -529,13 +558,6 @@ void outputProc( sharedData.config.minLat >= bbox.maxLat)) return; - // Read existing tile if merging - std::string rawExistingTile; - if (sharedData.mergeSqlite) { - sharedData.mbtiles.readTileAndUncompress(rawExistingTile, zoom, bbox.index.x, bbox.index.y, sharedData.config.compress, sharedData.config.gzip); - } - vtzero::vector_tile existingTile{rawExistingTile}; - // Loop through layers #ifndef _WIN32 if (!enabledUserSignal) { @@ -545,39 +567,64 @@ void outputProc( #endif signalStop=false; - for (auto lt = sharedData.layers.layerOrder.begin(); lt != sharedData.layers.layerOrder.end(); ++lt) { - if (signalStop) break; - ProcessLayer(sources, attributeStore, coordinates, zoom, data, existingTile, tile, bbox, *lt, sharedData); + const uint32_t extent = bbox.hires ? 8192u : 4096u; + string outputdata, compressed; + +#ifdef TILEMAKER_WITH_MLT + if (sharedData.config.tileFormat == TileFormat::MLT) { + std::vector layers; + for (auto lt = sharedData.layers.layerOrder.begin(); lt != sharedData.layers.layerOrder.end(); ++lt) { + if (signalStop) break; + MltLayerBuilder layer{sharedData.layers.layers[lt->at(0)].name, extent}; + ProcessLayer(sources, attributeStore, coordinates, zoom, data, layer, bbox, *lt, sharedData); + if (!layer.layer.features.empty()) layers.push_back(std::move(layer.layer)); + } + // .pmtiles output is always gzipped + const bool compressed = sharedData.config.compress || sharedData.outputMode == OptionsParser::OutputMode::PMTiles; + encodeMltTile(layers, sharedData.config.includeID, compressed, outputdata); + } else +#endif + { + vtzero::tile_builder tile; + + // Read existing tile if merging + std::string rawExistingTile; + if (sharedData.mergeSqlite) { + sharedData.mbtiles.readTileAndUncompress(rawExistingTile, zoom, bbox.index.x, bbox.index.y, sharedData.config.compress, sharedData.config.gzip); + } + vtzero::vector_tile existingTile{rawExistingTile}; + + for (auto lt = sharedData.layers.layerOrder.begin(); lt != sharedData.layers.layerOrder.end(); ++lt) { + if (signalStop) break; + const std::string& layerName = sharedData.layers.layers[lt->at(0)].name; + vtzero::layer_builder vtLayer{tile, layerName, sharedData.config.mvtVersion, extent}; + CopyExistingLayer(existingTile, layerName, vtLayer); + ProcessLayer(sources, attributeStore, coordinates, zoom, data, vtLayer, bbox, *lt, sharedData); + } + tile.serialize(outputdata); } // Write to file or sqlite - string outputdata, compressed; if (sharedData.outputMode == OptionsParser::OutputMode::MBTiles) { // Write to sqlite - //tile.SerializeToString(&outputdata); - tile.serialize(outputdata); - if (sharedData.config.compress) { compressed = compress_string(outputdata, Z_DEFAULT_COMPRESSION, sharedData.config.gzip); } sharedData.mbtiles.saveTile(zoom, bbox.index.x, bbox.index.y, sharedData.config.compress ? &compressed : &outputdata, sharedData.mergeSqlite); } else if (sharedData.outputMode == OptionsParser::OutputMode::PMTiles) { // Write to pmtiles - tile.serialize(outputdata); sharedData.pmtiles.saveTile(zoom, bbox.index.x, bbox.index.y, outputdata); } else { // Write to file stringstream dirname, filename; dirname << sharedData.outputFile << "/" << zoom << "/" << bbox.index.x; - filename << sharedData.outputFile << "/" << zoom << "/" << bbox.index.x << "/" << bbox.index.y << ".pbf"; + filename << sharedData.outputFile << "/" << zoom << "/" << bbox.index.x << "/" << bbox.index.y + << (sharedData.config.tileFormat == TileFormat::MLT ? ".mlt" : ".pbf"); boost::filesystem::create_directories(dirname.str()); fstream outfile(filename.str(), ios::out | ios::trunc | ios::binary); if (sharedData.config.compress) { - //tile.SerializeToString(&outputdata); - tile.serialize(outputdata); outfile << compress_string(outputdata, Z_DEFAULT_COMPRESSION, sharedData.config.gzip); } else { - tile.serialize(outputdata); outfile << outputdata; } outfile.close(); diff --git a/src/tilemaker.cpp b/src/tilemaker.cpp index f4891425..be7dca11 100644 --- a/src/tilemaker.cpp +++ b/src/tilemaker.cpp @@ -360,6 +360,10 @@ int main(const int argc, const char* argv[]) { sharedData.outputFile = options.outputFile; sharedData.outputMode = options.outputMode; sharedData.mergeSqlite = options.mergeSqlite; + if (sharedData.mergeSqlite && config.tileFormat == TileFormat::MLT) { + cerr << "--merge is not supported with \"tile_format\": \"mlt\"" << endl; + return 0; + } // ---- Initialise mbtiles/pmtiles if required @@ -368,6 +372,7 @@ int main(const int argc, const char* argv[]) { sharedData.writeMBTilesProjectData(); } else if (sharedData.outputMode == OptionsParser::OutputMode::PMTiles) { sharedData.pmtiles.open(sharedData.outputFile); + if (config.tileFormat == TileFormat::MLT) sharedData.pmtiles.tileType = TILETYPE_MLT; } // ---- Write out data diff --git a/test/config_validator.test.cpp b/test/config_validator.test.cpp index 6bf2116f..0d8e36aa 100644 --- a/test/config_validator.test.cpp +++ b/test/config_validator.test.cpp @@ -153,6 +153,39 @@ MU_TEST(test_invalid_source_columns) { mu_check(error.find("oneOf") != std::string::npos); } +MU_TEST(test_tile_format) { + const char* config = R"({ + "settings": { + "basezoom": 14, + "minzoom": 0, + "maxzoom": 14, + "include_ids": false, + "compress": "gzip", + "name": "Test", + "version": "1", + "description": "Test config", + "tile_format": "%s" + }, + "layers": { + "water": { + "minzoom": 0, + "maxzoom": 14 + } + } + })"; + char json[1024]; + + std::string error; + snprintf(json, sizeof(json), config, "mlt"); + mu_check(validateString(json, error)); + mu_check(error.empty()); + + error.clear(); + snprintf(json, sizeof(json), config, "geojson"); + mu_check(!validateString(json, error)); + mu_check(error.find("#/settings/tile_format") != std::string::npos); +} + MU_TEST(test_bundled_configs) { std::vector configs = { "resources/config-coastline.json", @@ -173,6 +206,7 @@ MU_TEST_SUITE(test_suite_config_validator) { MU_RUN_TEST(test_missing_multiple_required_settings); MU_RUN_TEST(test_invalid_layer_type); MU_RUN_TEST(test_invalid_source_columns); + MU_RUN_TEST(test_tile_format); MU_RUN_TEST(test_bundled_configs); } From 29f01329e617e8ccdb8109b9bbda0846609cbbd1 Mon Sep 17 00:00:00 2001 From: dktue <2878623+dktue@users.noreply.github.com> Date: Wed, 14 Oct 2026 10:57:04 +0200 Subject: [PATCH 2/2] perf: deterministic tile batching, lower memory, faster ingest and tile output --- include/attribute_store.h | 46 +++- include/clip_cache.h | 209 +++++++++++--- include/coordinates.h | 4 + include/coordinates_geom.h | 3 + include/geom.h | 40 ++- include/geometry_cache.h | 34 ++- include/lua_pool.h | 139 ++++++++++ include/node_store.h | 6 + include/osm_lua_processing.h | 103 ++++++- include/osm_mem_tiles.h | 10 + include/osm_store.h | 88 ++++-- include/output_object.h | 44 ++- include/pbf_processor.h | 19 +- include/pbf_reader.h | 1 + include/polylabel.h | 388 ++++++++++++++++++++++++-- include/sorted_node_store.h | 5 + include/tag_map.h | 26 +- include/tile_coordinates_set.h | 15 +- include/tile_data.h | 3 + include/visvalingam.h | 1 + src/attribute_store.cpp | 209 +++++++++++--- src/coordinates.cpp | 29 ++ src/coordinates_geom.cpp | 96 +++++++ src/geom.cpp | 398 ++++++++++++++++++++++++++- src/mmap_allocator.cpp | 53 +++- src/osm_lua_processing.cpp | 480 +++++++++++++++++++++++++++++---- src/osm_mem_tiles.cpp | 44 ++- src/output_object.cpp | 187 ++++++++++++- src/pbf_processor.cpp | 106 +++++++- src/pbf_reader.cpp | 22 +- src/shp_processor.cpp | 98 ++++++- src/sorted_node_store.cpp | 461 +++++++++++++++++++++++-------- src/sorted_way_store.cpp | 6 +- src/tag_map.cpp | 136 +++++----- src/tile_coordinates_set.cpp | 123 ++++++++- src/tile_data.cpp | 143 ++++++++-- src/tile_sorting.cpp | 188 +++++++++++-- src/tile_worker.cpp | 101 +++++-- src/tilemaker.cpp | 61 +++-- src/visvalingam.cpp | 274 ++++++++++++++++++- 40 files changed, 3846 insertions(+), 553 deletions(-) create mode 100644 include/lua_pool.h diff --git a/include/attribute_store.h b/include/attribute_store.h index c6020a41..295c495c 100644 --- a/include/attribute_store.h +++ b/include/attribute_store.h @@ -196,11 +196,14 @@ class AttributePairStore { lookupsUncached(0) { // The "hot" shard has a capacity of 64K, the others are unbounded. - pairs.push_back(DequeMap(1 << 16)); + // reserve + emplace: DequeMap's move constructor is not noexcept, so + // vector growth would copy every shard. + pairs.reserve(ATTRIBUTE_SHARDS); + pairs.emplace_back(1 << 16); // Reserve offset 0 as a sentinel pairs[0].add(AttributePair(0, false, 0)); for (size_t i = 1; i < ATTRIBUTE_SHARDS; i++) - pairs.push_back(DequeMap()); + pairs.emplace_back(); } void finalize() { finalized = true; } @@ -345,6 +348,20 @@ struct AttributeSet { return getValueAtIndex(actualIndex); } + // Calls fn(pairIndex) for every pair, in the same order as getPair(0..numPairs()-1) + template + void forEachPair(F&& fn) const { + if (useVector) { + for (const uint32_t v : intValues) + fn(v); + return; + } + + for (size_t i = 0; i < 8; i++) + if (isSet(i)) + fn(getValueAtIndex(i)); + } + AttributeSet(): useVector(false) { for (int i = 0; i < sizeof(shortValues)/sizeof(shortValues[0]); i++) shortValues[i] = 0; @@ -402,6 +419,26 @@ struct AttributeSet { struct AttributeStore { AttributeIndex add(AttributeSet &attributes); std::vector getUnsafe(AttributeIndex index) const; + // As getUnsafe, but calls fn(const AttributePair&) for each pair in order, without building a vector + template + void forEachPairUnsafe(AttributeIndex index, F&& fn) const { + forEachPairIndexUnsafe(index, [&](uint32_t pairIndex) { + fn(pairStore.getPairUnsafe(pairIndex)); + }); + } + // As forEachPairUnsafe, but calls fn(pairIndex) (for pairStore.getPairUnsafe) + template + void forEachPairIndexUnsafe(AttributeIndex index, F&& fn) const { + const AttributeSet* attrSet; + try { + uint32_t shard = index >> (32 - SHARD_BITS); + uint32_t offset = index & (~(~0u << (32 - SHARD_BITS))); + attrSet = &sets[shard].at(offset); + } catch (std::out_of_range &err) { + throw std::runtime_error("Failed to fetch attributes at index "+std::to_string(index)); + } + attrSet->forEachPair(fn); + } void reset(); // used for testing size_t size() const; void reportSize() const; @@ -411,6 +448,11 @@ struct AttributeStore { void addAttribute(AttributeSet& attributeSet, std::string const &key, double v, char minzoom); void addAttribute(AttributeSet& attributeSet, std::string const &key, int v, char minzoom); void addAttribute(AttributeSet& attributeSet, std::string const &key, bool v, char minzoom); + // As above, with the key already resolved by keyStore.key2index(key) (not overloads: kaguya::LuaRef keys convert to anything) + void addAttributeWithKeyIndex(AttributeSet& attributeSet, uint16_t keyIndex, std::string const &key, const protozero::data_view v, char minzoom); + void addAttributeWithKeyIndex(AttributeSet& attributeSet, uint16_t keyIndex, double v, char minzoom); + void addAttributeWithKeyIndex(AttributeSet& attributeSet, uint16_t keyIndex, int v, char minzoom); + void addAttributeWithKeyIndex(AttributeSet& attributeSet, uint16_t keyIndex, bool v, char minzoom); AttributeStore(): finalized(false), diff --git a/include/clip_cache.h b/include/clip_cache.h index 353356dd..f4c16cfb 100644 --- a/include/clip_cache.h +++ b/include/clip_cache.h @@ -5,75 +5,196 @@ #include "coordinates_geom.h" #include "geom.h" #include +#include +#ifdef TMOPT_VERIFY +#include +#include +#endif class TileBbox; +// Identifies the tile batch the current thread is working on. The clip cache +// is per-thread and is reset whenever the batch changes, so a tile can only +// reuse clips of its ancestors from the same batch. This makes output depend +// only on the batch partition, not on thread count or timing. +extern thread_local uint64_t clipCacheBatchId; +// Set when the current batch can't reuse clips (a single tile), so add() is skipped. +extern thread_local bool clipCacheNoStore; + template class ClipCache { public: ClipCache(size_t threadNum, unsigned int baseZoom): - baseZoom(baseZoom), - clipCache(threadNum * 16), - clipCacheMutex(threadNum * 16), - clipCacheSize(threadNum * 16) { + baseZoom(baseZoom) { } - const std::shared_ptr get(uint zoom, TileCoordinate x, TileCoordinate y, NodeID objectID) const{ - // Look for a previously clipped version at z-1, z-2, ... - - std::lock_guard lock(clipCacheMutex[objectID % clipCacheMutex.size()]); - while (zoom > 0) { - zoom--; - x /= 2; - y /= 2; - const auto& cache = clipCache[objectID % clipCache.size()]; - const auto& rv = cache.find(std::make_tuple(zoom, TileCoordinates(x, y), objectID)); - if (rv != cache.end()) { - return rv->second; + // An entry holds the clipped geometry and whether it is known to be valid + // (passed geom::is_valid). Callers that don't track validity get false. + struct Entry { + std::shared_ptr geom; + bool valid = false; + }; + + const std::shared_ptr get(uint zoom, TileCoordinate x, TileCoordinate y, NodeID objectID) const { + return getEntry(zoom, x, y, objectID).geom; + } + + Entry getEntry(uint zoom, TileCoordinate x, TileCoordinate y, NodeID objectID) const { + // Look for a previously clipped version at z-1, z-2, ... (nearest ancestor first) + LocalStore& s = local(zoom, TileCoordinates(x, y)); + Entry rv; + const Key key{(const void*)this, objectID}; + for (uint z = zoom < MAX_LEVELS ? zoom : MAX_LEVELS; z-- > 0; ) { + if (!(s.usedMask & (1u << z))) + continue; + const Level& level = s.levels[z]; + const uint shift = zoom - z; + if (level.index.x != ((uint64_t)x >> shift) || level.index.y != ((uint64_t)y >> shift)) + continue; + const auto it = level.m.find(key); + if (it != level.m.end()) { + rv = it->second; + break; } } - - return nullptr; +#ifdef TMOPT_VERIFY + { + Entry old; + uint vz = zoom; + TileCoordinate vx = x, vy = y; + while (vz > 0) { + vz--; + vx /= 2; + vy /= 2; + const auto& it = s.shadow.find(std::make_tuple((uint16_t)vz, (const void*)this, TileCoordinates(vx, vy), objectID)); + if (it != s.shadow.end()) { + old = it->second; + break; + } + } + if (old.geom != rv.geom || old.valid != rv.valid) { + std::cerr << "TMOPT_VERIFY: ClipCache mismatch z" << zoom << " " << x << "/" << y << " id " << objectID << std::endl; + abort(); + } + } +#endif + return rv; } - void add(const TileBbox& bbox, const NodeID objectID, const T& output) { + void add(const TileBbox& bbox, const NodeID objectID, const T& output, bool valid = false) { // The point of caching is to reuse the clip, so caching at the terminal zoom is // pointless. - if (bbox.zoom == baseZoom) + if (bbox.zoom == baseZoom || clipCacheNoStore) + return; + + // Levels beyond MAX_LEVELS could never be looked up by a child tile anyway + // (tile zooms are far below this). + if (bbox.zoom >= MAX_LEVELS) return; std::shared_ptr copy = std::make_shared(); boost::geometry::assign(*copy, output); - size_t index = objectID % clipCacheMutex.size(); - std::vector> objects; - std::lock_guard lock(clipCacheMutex[index]); - auto& cache = clipCache[index]; - // Reset the cache periodically so it doesn't grow without bound. - // - // I also tried boost's lru_cache -- but it seemed to perform worse, maybe - // due to the bookkeeping? We could try authoring a bounded map that - // evicts in FIFO order, which will have less bookkeeping. - clipCacheSize[index]++; - if (clipCacheSize[index] > 1024) { - clipCacheSize[index] = 0; - // Copy the map's contents to a vector so that calling .clear() - // and releasing the lock can happen separately from running the - // destructors of all of the objects. - objects.reserve(1025); - for (const auto& x : cache) - objects.push_back(x.second); - cache.clear(); - } + LocalStore& s = local(bbox.zoom, bbox.index); + // Safety cap (pruning normally keeps the store small). Count-based, so + // the result stays deterministic. + if (s.total > (1 << 18)) + s.clearFrom(0); +#ifdef TMOPT_VERIFY + if (s.shadow.size() > (1 << 18)) + s.shadow.clear(); + s.shadow[std::make_tuple((uint16_t)bbox.zoom, (const void*)this, bbox.index, objectID)] = Entry{copy, valid}; +#endif - cache[std::make_tuple(bbox.zoom, bbox.index, objectID)] = copy; + Level& level = s.levels[bbox.zoom]; + if ((s.usedMask & (1u << bbox.zoom)) && !(level.index == bbox.index)) { + // Can't happen: local() just dropped every level >= bbox.zoom unless + // bbox is the current tile. Kept so a level only ever holds one tile. + s.total -= level.m.size(); + level.m.clear(); + } + s.usedMask |= (1u << bbox.zoom); + level.index = bbox.index; + auto ins = level.m.emplace(Key{(const void*)this, objectID}, Entry{copy, valid}); + if (ins.second) + s.total++; + else + ins.first->second = Entry{copy, valid}; } private: + static const uint MAX_LEVELS = 32; + + // Entries of one ClipCache instance (`owner`) for one object. + struct Key { + const void* owner; + NodeID id; + bool operator==(const Key& o) const { return owner == o.owner && id == o.id; } + }; + struct KeyHash { + size_t operator()(const Key& k) const { + uint64_t h = (uint64_t)k.id * 0x9E3779B97F4A7C15ull ^ (uint64_t)(uintptr_t)k.owner; + h ^= h >> 29; + h *= 0xBF58476D1CE4E5B9ull; + h ^= h >> 32; + return (size_t)h; + } + }; + // All entries stored at one zoom. Only one tile per zoom can be live at a + // time (see below), so each level holds a single tile index. + struct Level { + TileCoordinates index = TileCoordinates(0, 0); + std::unordered_map m; + }; + + // Per-thread store, shared by all ClipCache instances on the thread (hence + // keyed on `this`), cleared when the thread starts a new batch. + // Tiles within a batch come in depth-first order, so when a new tile at + // zoom z starts, no entry at zoom >= z can be an ancestor of it or of any + // later tile: those are dropped, leaving only the current tile's ancestors. + struct LocalStore { + uint64_t gen = ~0ull; + uint zoom = 0; + TileCoordinates index = TileCoordinates(0, 0); + uint32_t usedMask = 0; + size_t total = 0; + Level levels[MAX_LEVELS]; +#ifdef TMOPT_VERIFY + std::map, Entry> shadow; +#endif + + // Drop every level >= zoom. + void clearFrom(uint zoom) { + for (uint z = zoom; z < MAX_LEVELS; z++) { + if (usedMask & (1u << z)) { + total -= levels[z].m.size(); + levels[z].m.clear(); + } + } + usedMask &= zoom >= MAX_LEVELS ? ~0u : ((1u << zoom) - 1); + } + }; + static LocalStore& local(uint zoom, TileCoordinates index) { + static thread_local LocalStore s; + if (s.gen != clipCacheBatchId) { + s.clearFrom(0); +#ifdef TMOPT_VERIFY + s.shadow.clear(); +#endif + s.gen = clipCacheBatchId; + } else if (s.zoom == zoom && s.index == index) { + return s; + } + s.zoom = zoom; + s.index = index; + s.clearFrom(zoom); +#ifdef TMOPT_VERIFY + s.shadow.erase(s.shadow.lower_bound(std::make_tuple((uint16_t)zoom, (const void*)nullptr, TileCoordinates(0, 0), (NodeID)0)), s.shadow.end()); +#endif + return s; + } + unsigned int baseZoom; - std::vector, std::shared_ptr>> clipCache; - mutable std::vector clipCacheMutex; - std::vector clipCacheSize; }; #endif diff --git a/include/coordinates.h b/include/coordinates.h index 38390125..9d6f458a 100644 --- a/include/coordinates.h +++ b/include/coordinates.h @@ -127,6 +127,10 @@ double meter2degp(double meter, double latp); // the range between smallest y and largest y is filled, for each x void fillCoveredTiles(std::unordered_set &tileSet); +// Same for a vector sorted by TileCoordinatesCompare without duplicates; the result stays sorted and unique. +void fillCoveredTiles(std::vector &tiles); +// Sort (TileCoordinatesCompare) and dedupe tiles[start, end). +void sortUniqueTiles(std::vector &tiles, size_t start = 0); #endif //_COORDINATES_H diff --git a/include/coordinates_geom.h b/include/coordinates_geom.h index 1c035d9c..f233ad3e 100644 --- a/include/coordinates_geom.h +++ b/include/coordinates_geom.h @@ -6,6 +6,9 @@ void insertIntermediateTiles(Linestring const &points, uint baseZoom, std::unordered_set &tileSet); void insertIntermediateTiles(Ring const &points, uint baseZoom, std::unordered_set &tileSet); +// Vector variants: append the tiles to `tiles`, then sort and dedupe the appended range [old size, end). +void insertIntermediateTiles(Linestring const &points, uint baseZoom, std::vector &tiles); +void insertIntermediateTiles(Ring const &points, uint baseZoom, std::vector &tiles); // ------------------------------------------------------ // Helper class for dealing with spherical Mercator tiles diff --git a/include/geom.h b/include/geom.h index 44f00be8..31902fec 100644 --- a/include/geom.h +++ b/include/geom.h @@ -8,6 +8,12 @@ using uint = unsigned int; #include #include +#ifdef TMOPT_VERIFY +#include +#include +#include +#include +#endif // boost::geometry #define BOOST_GEOMETRY_NO_ROBUSTNESS @@ -82,14 +88,42 @@ void make_valid(MultiPolygon &mp); // make_valid first (preserves area), then a zero-width buffer as a last // resort. Returns true if mp is valid afterwards; on failure mp is left as the // best-effort input so callers never regress. -bool repair_multi_polygon(MultiPolygon &mp); +bool repair_multi_polygon(MultiPolygon &mp, bool knownInvalid = false); void union_many(std::vector &mps); Point intersect_edge(Point const &a, Point const &b, char edge, Box const &bbox); char bit_code(Point const &p, Box const &bbox); -void fast_clip(Ring &points, Box const &bbox); -void fast_clip(MultiPolygon &mp, Box const &bbox); +// fast_clip returns true if the geometry was modified (false means identity). +bool fast_clip(Ring &points, Box const &bbox); +bool fast_clip(MultiPolygon &mp, Box const &bbox); + +// Cheap conservative validity tests, run before geom::is_valid. They return +// true only when the geometry is certainly valid (is_valid would return true +// with no_failure); false means "unsure" and the caller must run is_valid. +// Tolerances are relative to the coordinate magnitude. +bool polygonDefinitelyValid(const Polygon &p); +bool multiPolygonDefinitelyValid(const MultiPolygon &mp); + +// Call sites of the precheck (TMOPT_VERIFY hit statistics). +enum PrecheckSite { PRECHECK_TDS = 0, PRECHECK_REPAIR_ONE, PRECHECK_REPAIR_MP, PRECHECK_REPAIR_MP_PART, + PRECHECK_CHECK_RELATION, PRECHECK_CHECK_WAY, PRECHECK_REVERSED, PRECHECK_SHP, PRECHECK_SITES }; + +#ifdef TMOPT_VERIFY +// Per-call-site hit counters, and a cross-check of the precheck against Boost. +void precheckVerifyCount(int site, bool hit); +template +void precheckVerify(const GeometryT &g, bool hit, int site) { + precheckVerifyCount(site, hit); + if (hit && !boost::geometry::is_valid(g)) { + std::string reason; + boost::geometry::is_valid(g, reason); + std::cerr << "TMOPT_VERIFY: precheck site " << site << " disagrees with is_valid (" << reason << "): " + << std::setprecision(17) << boost::geometry::wkt(g) << std::endl; + std::abort(); + } +} +#endif #endif //_GEOM_TYPES_H diff --git a/include/geometry_cache.h b/include/geometry_cache.h index 607ff5f3..e8c39d9c 100644 --- a/include/geometry_cache.h +++ b/include/geometry_cache.h @@ -9,16 +9,19 @@ #include "coordinates.h" #include -#define NUM_BUCKETS 256 -// Keep bucket size small so linear search is fast -#define BUCKET_SIZE 32 - +#include template class GeometryCache { public: + // 64 x 32 entries: the cache only lives for one tile batch (see clear()), so a + // smaller table loses few hits and pins much less heap in every tile thread. + // Keep bucket size small so linear search is fast. + enum : size_t { NUM_BUCKETS = 64, BUCKET_SIZE = 32 }; + GeometryCache(): offsets(NUM_BUCKETS), + filled(NUM_BUCKETS), ids(NUM_BUCKETS * BUCKET_SIZE), geometries(NUM_BUCKETS * BUCKET_SIZE) {} @@ -35,16 +38,37 @@ class GeometryCache { void add(NodeID objectID, std::shared_ptr geom) { const size_t bucket = objectID % NUM_BUCKETS; size_t offset = bucket * BUCKET_SIZE + offsets[bucket]; - geometries[offset] = geom; + geometries[offset] = std::move(geom); ids[offset] = objectID; + if (filled[bucket] == 0) + dirtyBuckets.push_back(bucket); + if (filled[bucket] < BUCKET_SIZE) + filled[bucket]++; offsets[bucket]++; if (offsets[bucket] == BUCKET_SIZE) offsets[bucket] = 0; } + // Drop all entries (and release their geometries). Only touches the + // buckets used since the last clear, so it is cheap on a lightly used cache. + void clear() { + for (size_t bucket : dirtyBuckets) { + const size_t start = bucket * BUCKET_SIZE; + for (size_t i = start; i < start + filled[bucket]; i++) { + ids[i] = 0; + geometries[i].reset(); + } + filled[bucket] = 0; + offsets[bucket] = 0; + } + dirtyBuckets.clear(); + } + private: std::vector offsets; + std::vector filled; // used slots per bucket since clear() + std::vector dirtyBuckets; // buckets with filled > 0 std::vector ids; // CONSIDER: we could use a raw pointer - getOrBuildLinestring controls the lifetime std::vector> geometries; diff --git a/include/lua_pool.h b/include/lua_pool.h new file mode 100644 index 00000000..063d58a2 --- /dev/null +++ b/include/lua_pool.h @@ -0,0 +1,139 @@ +/*! \file */ +#ifndef _LUA_POOL_H +#define _LUA_POOL_H + +// Size-class pool allocator for a lua_State (lua_Alloc). +// +// Lua makes a very large number of small, short-lived allocations (strings, +// tables, closures) while running way_function etc. Each lua_State belongs to +// one OsmLuaProcessing object; it is never used by two threads at once +// (per-thread states are thread_local; the main-thread state is only used +// from the main thread or under the shapefile/GeoJSON attributeMutex), so the +// pool needs no locks. +// +// Blocks of up to MaxSmall bytes are carved from 64 KB slabs in 16-byte size +// classes and recycled through per-class free lists. Larger blocks go to +// malloc/realloc/free. Slabs are only released when the pool is destroyed, +// which must happen after lua_close(). + +#include +#include +#include +#include +#include +#ifdef TMOPT_VERIFY +#include +#include +#endif + +class LuaPool { +public: + static const size_t Granule = 16; + static const size_t MaxSmall = 512; + static const size_t NumClasses = MaxSmall / Granule; + static const size_t SlabSize = 64 * 1024; + + LuaPool() : bump(nullptr), bumpEnd(nullptr) { + for (size_t i = 0; i < NumClasses; i++) freeLists[i] = nullptr; + } + ~LuaPool() { + for (void *slab : slabs) std::free(slab); + } + + // lua_Alloc. When ptr==NULL, osize encodes the object type, not a size. + static void *alloc(void *ud, void *ptr, size_t osize, size_t nsize) { + LuaPool &pool = *static_cast(ud); + if (ptr && osize == 0) { + // Lua always passes the block size; only a misbehaving C module (via lua_getallocf) gets here + std::fputs("LuaPool: allocator called with osize 0 (C module using lua_getallocf incorrectly)\n", stderr); + std::abort(); + } + if (nsize == 0) { + if (ptr) pool.release(ptr, osize); + return nullptr; + } + if (!ptr) return pool.allocate(nsize); + + bool oSmall = osize <= MaxSmall, nSmall = nsize <= MaxSmall; + if (oSmall && nSmall && classOf(osize) == classOf(nsize)) { +#ifdef TMOPT_VERIFY + pool.verifyFree(ptr, osize); + pool.live[ptr] = nsize; +#endif + return ptr; + } + if (!oSmall && !nSmall) return std::realloc(ptr, nsize); + + void *n = pool.allocate(nsize); + if (!n) return nsize <= osize ? ptr : nullptr; // shrinking must not fail (Lua 5.1/5.2 assume it); else Lua keeps the old block + std::memcpy(n, ptr, osize < nsize ? osize : nsize); + pool.release(ptr, osize); + return n; + } + +private: + struct FreeBlock { FreeBlock *next; }; + + LuaPool(const LuaPool&); + LuaPool& operator=(const LuaPool&); + + static size_t classOf(size_t n) { return (n - 1) / Granule; } // n in [1, MaxSmall] + + void *allocate(size_t n) { + if (n > MaxSmall) return std::malloc(n); + size_t c = classOf(n); + void *p; + if (freeLists[c]) { + FreeBlock *b = freeLists[c]; + freeLists[c] = b->next; + p = b; + } else { + size_t bytes = (c + 1) * Granule; + if (bump == nullptr || (size_t)(bumpEnd - bump) < bytes) { + char *slab = static_cast(std::malloc(SlabSize)); + if (!slab) return nullptr; + try { slabs.push_back(slab); } + catch (...) { std::free(slab); return nullptr; } + bump = slab; + bumpEnd = slab + SlabSize; + } + p = bump; + bump += bytes; + } +#ifdef TMOPT_VERIFY + live[p] = n; +#endif + return p; + } + + void release(void *p, size_t osize) { + if (osize > MaxSmall) { std::free(p); return; } +#ifdef TMOPT_VERIFY + verifyFree(p, osize); + live.erase(p); +#endif + size_t c = classOf(osize); + FreeBlock *b = static_cast(p); + b->next = freeLists[c]; + freeLists[c] = b; + } + +#ifdef TMOPT_VERIFY + // Lua must hand back every small block with the size it was last given. + void verifyFree(void *p, size_t osize) { + auto it = live.find(p); + if (it == live.end() || it->second != osize) { + std::cerr << "TMOPT_VERIFY LuaPool: bad free " << p << " osize " << osize + << " recorded " << (it == live.end() ? 0 : it->second) << std::endl; + std::abort(); + } + } + std::unordered_map live; +#endif + + FreeBlock *freeLists[NumClasses]; + char *bump, *bumpEnd; + std::vector slabs; +}; + +#endif //_LUA_POOL_H diff --git a/include/node_store.h b/include/node_store.h index 76fe18b3..fb5e6d74 100644 --- a/include/node_store.h +++ b/include/node_store.h @@ -23,6 +23,12 @@ class NodeStore // Accessors virtual size_t size() const = 0; virtual LatpLon at(NodeID i) const = 0; + // Batched lookup: out[i] = at(ids[i]) for i < n. Throws exactly what + // at() throws for the first missing node. + virtual void at(const NodeID* ids, size_t n, LatpLon* out) const { + for (size_t i = 0; i < n; i++) + out[i] = at(ids[i]); + } virtual bool contains(size_t shard, NodeID id) const = 0; virtual NodeStore& shard(size_t shard) = 0; diff --git a/include/osm_lua_processing.h b/include/osm_lua_processing.h index 912652d8..cfc31659 100644 --- a/include/osm_lua_processing.h +++ b/include/osm_lua_processing.h @@ -31,6 +31,7 @@ extern "C" { } #include "external/kaguya.hpp" +#include "lua_pool.h" // FIXME: why is this global ? extern bool verbose; @@ -166,18 +167,78 @@ class OsmLuaProcessing { CentroidAlgorithm defaultCentroidAlgorithm() const { return CentroidAlgorithm::Polylabel; } CentroidAlgorithm parseCentroidAlgorithm(const std::string& algorithm) const; Point calculateCentroid(CentroidAlgorithm algorithm); + Point calculateCentroidUncached(CentroidAlgorithm algorithm); + + // Reversed: the only defect was ring orientation; geom now has every ring reversed + // and is valid (make_valid would produce the same rings). + enum class CorrectGeometryResult: char { Invalid = 0, Valid = 1, Corrected = 2, Reversed = 3 }; + + // Reverse every ring of a multipolygon in place. Only multipolygons are supported. + // Returns false (without reversing) if any ring is tiny: make_valid drops rings with + // |area| <= 1E-12, so those must take the make_valid path to get the same result. + template + static bool ReverseRings(GeometryT &geom) { return false; } + static bool ReverseRings(MultiPolygon &geom) { + for (const auto &p : geom) { + if (std::abs(boost::geometry::area(p.outer())) <= 1E-12) return false; + for (const auto &inner : p.inners()) + if (std::abs(boost::geometry::area(inner)) <= 1E-12) return false; + } + for (auto &p : geom) { + std::reverse(p.outer().begin(), p.outer().end()); + for (auto &inner : p.inners()) + std::reverse(inner.begin(), inner.end()); + } + return true; + } + + // Cheap conservative "certainly valid" test before geom::is_valid; only areal + // geometries have one (false = unsure, run is_valid). + template + static bool PrecheckValid(const GeometryT &, int) { return false; } + static bool PrecheckValid(const MultiPolygon &geom, int site) { + const bool hit = multiPolygonDefinitelyValid(geom); +#ifdef TMOPT_VERIFY + precheckVerify(geom, hit, site); +#else + (void)site; +#endif + return hit; + } - enum class CorrectGeometryResult: char { Invalid = 0, Valid = 1, Corrected = 2 }; // ---- Requests from Lua to write this way/node to a vector tile's Layer template - CorrectGeometryResult CorrectGeometry(GeometryT &geom) + CorrectGeometryResult CorrectGeometry(GeometryT &geom, bool allowReverse = false) + { + return CorrectGeometryAfterCheck(geom, CheckGeometry(geom), allowReverse); + } + + // Validity check part of CorrectGeometry (reports the failure if verbose) + template + geom::validity_failure_type CheckGeometry(const GeometryT &geom) { geom::validity_failure_type failure = geom::validity_failure_type::no_failure; - if (isRelation && !geom::is_valid(geom,failure)) { + // When the precheck proves validity, failure stays no_failure: Boost's result. + const bool pre = (isRelation || isWay) && PrecheckValid(geom, isRelation ? PRECHECK_CHECK_RELATION : PRECHECK_CHECK_WAY); + if (isRelation && !pre && !geom::is_valid(geom,failure)) { if (verbose) std::cout << "Relation " << originalOsmID << " has " << boost_validity_error(failure) << std::endl; - } else if (isWay && !geom::is_valid(geom,failure)) { + } else if (isWay && !pre && !geom::is_valid(geom,failure)) { if (verbose && failure!=22) std::cout << "Way " << originalOsmID << " has " << boost_validity_error(failure) << std::endl; } + return failure; + } + + // Correction part of CorrectGeometry, given the result of CheckGeometry(geom) + template + CorrectGeometryResult CorrectGeometryAfterCheck(GeometryT &geom, geom::validity_failure_type failure, bool allowReverse = false) + { + // Fast path: orientation-only defect. Reversing the rings is what make_valid + // produces in this case; fall back to make_valid if it isn't enough. + if (allowReverse && failure == boost::geometry::failure_wrong_orientation && ReverseRings(geom)) { + if (PrecheckValid(geom, PRECHECK_REVERSED) || geom::is_valid(geom)) + return CorrectGeometryResult::Reversed; + ReverseRings(geom); + } if (failure==boost::geometry::failure_spikes) geom::remove_spikes(geom); @@ -204,6 +265,8 @@ class OsmLuaProcessing { void AttributeNumeric(const std::string &key, const double val, const char minzoom); void AttributeBoolean(const std::string &key, const bool val, const char minzoom); void AttributeInteger(const std::string &key, const int val, const char minzoom); + // Fast path for the native Lua bindings: key given as a Lua string buffer + template void attributeFast(const char *key, size_t keyLen, T val, char minzoom, uint type); void MinZoom(const double z); void ZOrder(const double z); void Score(const double score); @@ -265,6 +328,7 @@ class OsmLuaProcessing { multiLinestringInited = false; polygonInited = false; multiPolygonInited = false; + centroidMemoValid[0] = centroidMemoValid[1] = false; relationAccepted = false; relationList.clear(); relationSubscript = -1; @@ -275,7 +339,10 @@ class OsmLuaProcessing { isPostScanRelation = false; } - void removeAttributeIfNeeded(const std::string& key); + void removeAttributeIfNeeded(uint16_t keyIndex); + + // Call a global Lua function with no arguments; false if it isn't a function + bool callLuaGlobal(const char *name, bool keepGlobals); void noteIfDecluttered(LatpLon point); @@ -287,6 +354,19 @@ class OsmLuaProcessing { OSMStore &osmStore; // global OSM store + // Owns the lua_State that luaState wraps (non-owning). Declared in this order so the + // pool is constructed first and freed last, after lua_close(). + struct LuaStateHolder { + lua_State *L; + explicit LuaStateHolder(lua_State *L) : L(L) { } + ~LuaStateHolder() { if (L) lua_close(L); } + private: + LuaStateHolder(const LuaStateHolder&); + LuaStateHolder& operator=(const LuaStateHolder&); + }; + static lua_State *newLuaState(LuaPool &pool); + LuaPool luaPool; + LuaStateHolder luaStateHolder; kaguya::State luaState; bool supportsRemappingShapefiles; bool supportsReadingRelations; @@ -321,6 +401,14 @@ class OsmLuaProcessing { bool multiPolygonInited; geom::model::polygon areaPolygonCache; + // Scratch geometries reused across objects (allocation reuse only) + Polygon centroidScratchPolygon; + MultiPolygon layerScratchMultiPolygon; + + // calculateCentroid() result for the current object, per CentroidAlgorithm + Point centroidMemo[2]; + bool centroidMemoValid[2] = { false, false }; + NodeID lastStoredGeometryId; OutputGeometryType lastStoredGeometryType; @@ -333,12 +421,13 @@ class OsmLuaProcessing { // plus the score from Score(). Sparse - most objects add nothing here. struct PendingDeclutter { uint32_t outputIndex; LatpLon point; int32_t score; }; std::vector declutterOutputs; - std::vector outputKeys; + std::vector outputKeys; const PbfReader::Relation* currentRelation; const boost::container::flat_map* currentPostScanTags; // for postScan only const std::vector* stringTable; - std::vector finalizeOutputs(); + std::vector finalizedOutputs; // reused by finalizeOutputs() + std::vector &finalizeOutputs(); bool materializeGeometries; bool wayEmitted; diff --git a/include/osm_mem_tiles.h b/include/osm_mem_tiles.h index a01b1cd9..e8d7ed9b 100644 --- a/include/osm_mem_tiles.h +++ b/include/osm_mem_tiles.h @@ -14,6 +14,10 @@ #define IS_NODE(x) (((x) >> TILE_DATA_ID_SIZE) == (USE_NODE_STORE >> TILE_DATA_ID_SIZE)) #define USE_WAY_STORE (1ull << TILE_DATA_ID_SIZE) #define IS_WAY(x) (((x) >> TILE_DATA_ID_SIZE) == (USE_WAY_STORE >> TILE_DATA_ID_SIZE)) +// Way-store polygon whose ring must be reversed (only defect was orientation). +// Only created for POLYGON_ outputs. +#define USE_WAY_STORE_REVERSED (3ull << TILE_DATA_ID_SIZE) +#define IS_WAY_REV(x) (((x) >> TILE_DATA_ID_SIZE) == (USE_WAY_STORE_REVERSED >> TILE_DATA_ID_SIZE)) #define OSM_ID(x) ((x) & 0b1111111111111111111111111111111111) class NodeStore; @@ -53,6 +57,12 @@ class OsmMemTiles : public TileDataSource { void populateLinestring(Linestring& ls, NodeID objectID) const; Linestring& getOrBuildLinestring(NodeID objectID) const; void populateMultiPolygon(MultiPolygon& dst, NodeID objectID) override; + // Way-store polygons are only emitted by OsmLuaProcessing::Layer when + // CorrectGeometry found them Valid (or Reversed: valid after reversing the ring) + // on the same lon/latp coordinates. + bool polygonKnownValid(NodeID objectID) const override { + return objectID >= OSM_THRESHOLD && (IS_WAY(objectID) || IS_WAY_REV(objectID)); + } const NodeStore& nodeStore; const WayStore& wayStore; diff --git a/include/osm_store.h b/include/osm_store.h index 56f797fd..4f50cd5f 100644 --- a/include/osm_store.h +++ b/include/osm_store.h @@ -7,6 +7,8 @@ #include "mmap_allocator.h" #include "relation_roles.h" +#include +#include #include #include #include @@ -51,45 +53,88 @@ struct DataViewLessThan { class UsedWays { private: - std::vector usedList; + // Sparse list of used way IDs: only member ways of accepted relations + // are inserted, which is a tiny fraction of all ways. Appended unsorted + // during the relation scan, then sorted/deduplicated by finalize(). + // If the used set is large relative to the ID range (planet-scale inputs), + // finalize() switches to a dense bitmap, which is then smaller and faster. + std::vector ids; + std::vector dense; + bool isDense = false; + size_t dedupAt = 1 << 20; + std::atomic sorted{true}; mutable std::mutex mutex; + void dedupLocked() { + std::sort(ids.begin(), ids.end()); + ids.erase(std::unique(ids.begin(), ids.end()), ids.end()); + } + + void finalizeLocked() { + dedupLocked(); + if (!ids.empty() && ids.size() * 64 > ids.back()) { + dense.resize(ids.back() + 1); + for (WayID id : ids) dense[id] = true; + isDense = true; + ids.clear(); + } + ids.shrink_to_fit(); + sorted.store(true, std::memory_order_release); + } + public: bool inited = false; - // Size the vector to a reasonable estimate, to avoid resizing on the fly - // TODO: it'd be preferable if UsedWays didn't know about compact mode -- - // instead, use an efficient data structure if numNodes < 1B, otherwise - // use a large bitvector void reserve(bool compact, size_t numNodes) { std::lock_guard lock(mutex); if (inited) return; inited = true; - if (compact) { - // If we're running in compact mode, way count is roughly 1/9th of node count... say 1/8 to be safe - usedList.reserve(numNodes/8); - } else { - // Otherwise, we could have anything up to the current max node ID (approaching 2**30 in summer 2021) - // 2**31 is 0.25GB with a vector - usedList.reserve(pow(2,31)); - } + ids.reserve(1 << 16); } - + // Mark a way as used void insert(WayID wayid) { std::lock_guard lock(mutex); - if (wayid>=usedList.size()) usedList.resize(wayid+256); - usedList[wayid] = true; + if (isDense) { + if (wayid >= dense.size()) dense.resize(wayid + 1); + dense[wayid] = true; + return; + } + ids.push_back(wayid); + sorted.store(false, std::memory_order_relaxed); + // Bound memory when ways are members of many relations (e.g. routes) + if (ids.size() >= dedupAt) { + dedupLocked(); + dedupAt = std::max(dedupAt, ids.size() * 2); + } + } + + // Sort and deduplicate; call once all inserts are done, before any at() + void finalize() { + std::lock_guard lock(mutex); + finalizeLocked(); } - + // See if a way is used - bool at(WayID wayid) const { - return (wayid>=usedList.size()) ? false : usedList[wayid]; + bool at(WayID wayid) { + if (!sorted.load(std::memory_order_acquire)) { + std::lock_guard lock(mutex); + if (!sorted.load(std::memory_order_relaxed)) finalizeLocked(); + } + if (isDense) + return wayid < dense.size() && dense[wayid]; + return std::binary_search(ids.begin(), ids.end(), wayid); } - + void clear() { std::lock_guard lock(mutex); - usedList.clear(); + ids.clear(); + ids.shrink_to_fit(); + dense.clear(); + dense.shrink_to_fit(); + isDense = false; + dedupAt = 1 << 20; + sorted.store(true, std::memory_order_release); } }; @@ -305,6 +350,7 @@ class OSMStore bool way_is_used(WayID i) { return used_ways.at(i); } void ensureUsedWaysInited(); + void finalizeUsedWays() { used_ways.finalize(); } using tag_map_t = boost::container::flat_map; diff --git a/include/output_object.h b/include/output_object.h index e0f4ac22..146e39b2 100644 --- a/include/output_object.h +++ b/include/output_object.h @@ -69,16 +69,19 @@ class OutputObject { this->attributes = attributes; } + // memo (MVT only, optional) must belong to the layer_builder that fbuilder writes to void writeAttributes( const AttributeStore& attributeStore, vtzero::feature_builder& fbuilder, - char zoom + char zoom, + class MvtAttributeMemo* memo = nullptr ) const; #ifdef TILEMAKER_WITH_MLT void writeAttributes( const AttributeStore& attributeStore, class MltFeatureBuilder& fbuilder, - char zoom + char zoom, + class MvtAttributeMemo* memo = nullptr // unused ) const; #endif @@ -87,6 +90,43 @@ class OutputObject { }; #pragma pack(pop) +// Remembers, for one vtzero::layer_builder at one zoom, the key/value table +// indices that each attribute set's properties resolve to, so that features +// sharing an attribute set don't re-encode and re-find their keys and values. +// layer_builder::add_key/add_value are find-or-append, so a remembered index +// is exactly what vtzero would return again: the tile bytes are unchanged. +// Call reset() whenever a new layer_builder is started. +class MvtAttributeMemo { +public: + void reset(vtzero::layer_builder& layer, char zoom); + +private: + friend class OutputObject; + + // Open-addressing uint32_t -> (a, b) map. clear() is O(1): slots carry a generation. + struct Table { + struct Slot { uint32_t gen, key, a, b; }; + std::vector slots; + uint32_t gen = 1, used = 0; + + const Slot* find(uint32_t key) const; + void insert(uint32_t key, uint32_t a, uint32_t b); // key must not be present + void clear(); + }; + + // Returns the [key, value] index pairs for the set (count pairs) + const uint32_t* resolve(const AttributeStore& attributeStore, AttributeIndex index, uint32_t& count); + void resolvePair(const AttributeStore& attributeStore, uint32_t pairIndex, uint32_t& key, uint32_t& value); + + vtzero::layer_builder* layer = nullptr; + char zoom = 0; + uint32_t layerGen = 0; // validity of keyGen entries + Table sets; // AttributeIndex -> [begin, end) in entries + Table pairs; // pair index -> key, value index (key NO_INDEX: not written at this zoom) + std::vector entries; + std::vector keyGen, keyIdx; // by AttributeKeyStore index +}; + struct OutputObjectID { OutputObject oo; uint64_t id; diff --git a/include/pbf_processor.h b/include/pbf_processor.h index abd43bc6..2668b79e 100644 --- a/include/pbf_processor.h +++ b/include/pbf_processor.h @@ -8,6 +8,8 @@ #include #include #include +#include +#include #include #include "osm_store.h" #include "significant_tags.h" @@ -121,10 +123,25 @@ class PbfProcessor /// Find a string in the dictionary static int findStringPosition(const PbfReader::PrimitiveBlock& pb, const std::string& str); - + + // Subdivided blocks (chunks > 1): the first chunk task to arrive decompresses + // the block once; the other chunk tasks share the buffer. The entry is erased + // when the last chunk releases it, so each pass (shard) re-reads it. + struct SharedBlock { + std::shared_ptr data; + size_t remaining; + bool ready; + bool failed; + }; + std::shared_ptr acquireSharedBlock(std::istream& infile, const BlockMetadata& blockMetadata); + void releaseSharedBlock(const BlockMetadata& blockMetadata); + OSMStore &osmStore; std::mutex ioMutex; std::atomic compactWarningIssued; + std::mutex sharedBlocksMutex; + std::condition_variable sharedBlocksCv; + std::map sharedBlocks; }; int ReadPbfBoundingBox(const std::string &inputFile, double &minLon, double &maxLon, diff --git a/include/pbf_reader.h b/include/pbf_reader.h index 9af930c5..756c9776 100644 --- a/include/pbf_reader.h +++ b/include/pbf_reader.h @@ -277,6 +277,7 @@ namespace PbfReader { public: BlobHeader readBlobHeader(std::istream& input); protozero::data_view readBlob(int32_t datasize, std::istream& input); + protozero::data_view readBlob(int32_t datasize, std::istream& input, std::string& output, bool copyRaw); HeaderBlock readHeaderBlock(protozero::data_view data); HeaderBBox readHeaderBBox(protozero::data_view data); PrimitiveBlock& readPrimitiveBlock(protozero::data_view data); diff --git a/include/polylabel.h b/include/polylabel.h index de3a8bc7..d56356c1 100644 --- a/include/polylabel.h +++ b/include/polylabel.h @@ -25,6 +25,20 @@ #include #include #include +#include +#include +#include +// 2-wide SSE2 distance kernel. Only where the compiler cannot contract mul+add into FMA (that +// would round differently from the scalar original); elsewhere the original scalar loop is used. +#if (defined(__x86_64__) || defined(_M_X64)) && !defined(__FMA__) && !defined(__FMA4__) && !defined(__AVX2__) && (!defined(__FLT_EVAL_METHOD__) || __FLT_EVAL_METHOD__ == 0) +#include +#define POLYLABEL_SSE2 1 +#endif +#ifdef TMOPT_VERIFY +#include +#include +#include +#endif namespace mapbox { @@ -92,14 +106,224 @@ auto pointToPolygonDist(const Point& point, const Polygon& polygon) { return (inside ? 1 : -1) * std::sqrt(minDistSq); } +// as pointToPolygonDist (bit-identical result), also reports the nearest edge: +// ring kr (-1 = outer, else inner index) and edge (ring[ki], ring[ki == 0 ? len - 1 : ki - 1]) +double pointToPolygonDistKOrig(const Point& point, const Polygon& polygon, int& kr, std::size_t& ki) { + bool inside = false; + auto minDistSq = std::numeric_limits::infinity(); + kr = -1; + ki = 0; + + { + const auto& ring = polygon.outer(); + for (std::size_t i = 0, len = ring.size(), j = len - 1; i < len; j = i++) { + const auto& a = ring[i]; + const auto& b = ring[j]; + + if ((a.get<1>() > point.get<1>()) != (b.get<1>() > point.get<1>()) && + (point.get<0>() < (b.get<0>() - a.get<0>()) * (point.get<1>() - a.get<1>()) / (b.get<1>() - a.get<1>()) + a.get<0>())) inside = !inside; + + // same selection as std::min(minDistSq, dd), including NaN handling + const double dd = getSegDistSq(point, a, b); + if (dd < minDistSq) { minDistSq = dd; ki = i; } + } + } + + int r = 0; + for (const auto& ring : polygon.inners()) { + for (std::size_t i = 0, len = ring.size(), j = len - 1; i < len; j = i++) { + const auto& a = ring[i]; + const auto& b = ring[j]; + + if ((a.get<1>() > point.get<1>()) != (b.get<1>() > point.get<1>()) && + (point.get<0>() < (b.get<0>() - a.get<0>()) * (point.get<1>() - a.get<1>()) / (b.get<1>() - a.get<1>()) + a.get<0>())) inside = !inside; + + const double dd = getSegDistSq(point, a, b); + if (dd < minDistSq) { minDistSq = dd; kr = r; ki = i; } + } + r++; + } + + return (inside ? 1 : -1) * std::sqrt(minDistSq); +} + +#ifdef POLYLABEL_SSE2 +// The edges of the polygon being labelled, one entry per (a = ring[i], b = ring[j]) exactly as +// pointToPolygonDistKOrig iterates them (outer ring first, then inners in order), in blocks of two. +// An odd count is padded with a NaN entry: its distance is NaN (never selected by the min) and its +// crossing test compares false, so it changes neither the distance nor the parity. +struct EdgePair { + double ax[2], ay[2], bx[2], by[2], dx[2], dy[2], den[2]; +}; + +struct PolygonEdges { + const Polygon* polygon = nullptr; // polygon the arrays describe + std::vector pairs; + std::vector ring; // per entry: -1 = outer, else inner index + std::vector edge; // per entry: index i of a within its ring + + void build(const Polygon& p) { + std::size_t n = p.outer().size(); + for (const auto& inner : p.inners()) n += inner.size(); + const double nan = std::numeric_limits::quiet_NaN(); + pairs.assign((n + 1) / 2, EdgePair { {nan, nan}, {nan, nan}, {nan, nan}, {nan, nan}, {nan, nan}, {nan, nan}, {nan, nan} }); + ring.resize(n); + edge.resize(n); + std::size_t k = 0; + auto add = [&](const Ring& rg, int r) { + for (std::size_t i = 0, len = rg.size(), j = len - 1; i < len; j = i++, k++) { + EdgePair& q = pairs[k >> 1]; + const std::size_t l = k & 1; + q.ax[l] = rg[i].get<0>(); + q.ay[l] = rg[i].get<1>(); + q.bx[l] = rg[j].get<0>(); + q.by[l] = rg[j].get<1>(); + q.dx[l] = q.bx[l] - q.ax[l]; + q.dy[l] = q.by[l] - q.ay[l]; + q.den[l] = q.dx[l] * q.dx[l] + q.dy[l] * q.dy[l]; + ring[k] = r; + edge[k] = i; + } + }; + add(p.outer(), -1); + int r = 0; + for (const auto& inner : p.inners()) add(inner, r++); + polygon = &p; + } + + // don't keep huge buffers alive after an exceptionally large polygon + void release() { + polygon = nullptr; + if (pairs.capacity() > (1u << 12)) { + std::vector().swap(pairs); + std::vector().swap(ring); + std::vector().swap(edge); + } + } +}; + +// edges of the polygon currently inside polylabel() on this thread (nullptr otherwise) +static thread_local const PolygonEdges* activeEdges = nullptr; + +// pointToPolygonDistKOrig over the edge arrays, two edges per step and without branches. +// Every lane evaluates the same IEEE operations in the same operand order as getSegDistSq and the +// crossing test: t = ((p.x - a.x) * dx + (p.y - a.y) * dy) / (dx * dx + dy * dy), then +// x = t > 1 ? b.x : (t > 0 ? a.x + dx * t : a.x) (a degenerate edge gives t = NaN, so x = a.x as +// in the original's dx == dy == 0 branch). The minimum of the same set of values is exact in any +// order, MINPD keeps the accumulator for a NaN distance like 'dd < minDistSq', and the crossing +// parity is an XOR (order-independent). +inline double pointToPolygonDistKSoA(const Point& point, const PolygonEdges& e, int& kr, std::size_t& ki) { + const __m128d px = _mm_set1_pd(point.get<0>()); + const __m128d py = _mm_set1_pd(point.get<1>()); + const __m128d one = _mm_set1_pd(1.0); + const __m128d zero = _mm_setzero_pd(); + __m128d mn = _mm_set1_pd(std::numeric_limits::infinity()); + __m128d par = _mm_setzero_pd(); + __m128i idx = _mm_set_epi64x(1, 0); + __m128i best = _mm_setzero_si128(); + const __m128i two = _mm_set_epi64x(2, 2); + + const EdgePair* q = e.pairs.data(); + for (const EdgePair* end = q + e.pairs.size(); q != end; q++) { + const __m128d ax = _mm_loadu_pd(q->ax); + const __m128d ay = _mm_loadu_pd(q->ay); + const __m128d bx = _mm_loadu_pd(q->bx); + const __m128d by = _mm_loadu_pd(q->by); + const __m128d dx = _mm_loadu_pd(q->dx); + const __m128d dy = _mm_loadu_pd(q->dy); + const __m128d pxa = _mm_sub_pd(px, ax); + const __m128d pya = _mm_sub_pd(py, ay); + + // (a.y > p.y) != (b.y > p.y) && p.x < (b.x - a.x) * (p.y - a.y) / (b.y - a.y) + a.x + const __m128d straddle = _mm_xor_pd(_mm_cmpgt_pd(ay, py), _mm_cmpgt_pd(by, py)); + const __m128d xi = _mm_add_pd(_mm_div_pd(_mm_mul_pd(dx, pya), dy), ax); + par = _mm_xor_pd(par, _mm_and_pd(straddle, _mm_cmplt_pd(px, xi))); + + // getSegDistSq + const __m128d t = _mm_div_pd(_mm_add_pd(_mm_mul_pd(pxa, dx), _mm_mul_pd(pya, dy)), _mm_loadu_pd(q->den)); + const __m128d gt1 = _mm_cmpgt_pd(t, one); + const __m128d gt0 = _mm_cmpgt_pd(t, zero); + __m128d x = _mm_or_pd(_mm_and_pd(gt0, _mm_add_pd(ax, _mm_mul_pd(dx, t))), _mm_andnot_pd(gt0, ax)); + __m128d y = _mm_or_pd(_mm_and_pd(gt0, _mm_add_pd(ay, _mm_mul_pd(dy, t))), _mm_andnot_pd(gt0, ay)); + x = _mm_or_pd(_mm_and_pd(gt1, bx), _mm_andnot_pd(gt1, x)); + y = _mm_or_pd(_mm_and_pd(gt1, by), _mm_andnot_pd(gt1, y)); + const __m128d ex = _mm_sub_pd(px, x); + const __m128d ey = _mm_sub_pd(py, y); + const __m128d dd = _mm_add_pd(_mm_mul_pd(ex, ex), _mm_mul_pd(ey, ey)); + + // per lane: if (dd < mn) { mn = dd; best = idx; } + const __m128i lt = _mm_castpd_si128(_mm_cmplt_pd(dd, mn)); + mn = _mm_min_pd(dd, mn); + best = _mm_or_si128(_mm_and_si128(lt, idx), _mm_andnot_si128(lt, best)); + idx = _mm_add_epi64(idx, two); + } + + double m2[2]; + _mm_storeu_pd(m2, mn); + std::int64_t b2[2]; + _mm_storeu_si128(reinterpret_cast<__m128i*>(b2), best); + const int pm = _mm_movemask_pd(par); + const bool inside = ((pm ^ (pm >> 1)) & 1) != 0; + + // both lanes are non-NaN (they start at +inf); either lane achieving the minimum is a nearest edge + const int lane = m2[1] < m2[0] ? 1 : 0; + const double minDistSq = m2[lane]; + if (e.ring.empty()) { + kr = -1; + ki = 0; + } else { + kr = e.ring[b2[lane]]; + ki = e.edge[b2[lane]]; + } + + return (inside ? 1 : -1) * std::sqrt(minDistSq); +} +#endif + +// pointToPolygonDistKOrig, via the edge arrays when polylabel() has built them for this polygon +double pointToPolygonDistK(const Point& point, const Polygon& polygon, int& kr, std::size_t& ki) { +#ifdef POLYLABEL_SSE2 + const PolygonEdges* e = activeEdges; + if (e && e->polygon == &polygon) { + const double d = pointToPolygonDistKSoA(point, *e, kr, ki); +#ifdef TMOPT_VERIFY + int okr; + std::size_t oki; + const double od = pointToPolygonDistKOrig(point, polygon, okr, oki); + bool bad = std::memcmp(&d, &od, sizeof(d)) != 0; + if (!bad && std::isfinite(d)) { + // the reported edge must achieve the minimum (edgeDistBound relies on it being an edge) + const auto& ring = kr < 0 ? polygon.outer() : polygon.inners()[kr]; + const double ed = std::sqrt(getSegDistSq(point, ring[ki], ring[ki == 0 ? ring.size() - 1 : ki - 1])); + bad = ed != std::fabs(d); + } + if (bad) { + std::cerr << "TMOPT_VERIFY polylabel: SSE2 distance differs from original" << std::endl; + std::abort(); + } + static struct EvalCount { + std::atomic n{0}; + ~EvalCount() { std::cerr << "TMOPT_VERIFY polylabel SSE2 distance evals: " << n.load() << std::endl; } + } evalCount; + evalCount.n++; +#endif + return d; + } +#endif + return pointToPolygonDistKOrig(point, polygon, kr, ki); +} + struct Cell { Cell(const Point& c_, double h_, const Polygon& polygon) : c(c_), h(h_), - d(pointToPolygonDist(c, polygon)), + d(pointToPolygonDistK(c, polygon, kr, ki)), max(d + h * std::sqrt(2)) {} + // nearest edge of the polygon outline to c (declared first: written while d is computed) + std::size_t ki; + int kr; Point c; // cell center double h; // half the cell size double d; // distance from cell center to polygon @@ -131,6 +355,62 @@ Cell getCentroidCell(const Polygon& polygon) { return Cell(area == 0 ? ring.at(0) : Point { ox + cx / area, oy + cy / area }, 0, polygon); } +// Upper bound on pointToPolygonDist(p, polygon), from the edge nearest to `near.c`. The full scan +// computes min over all edges of getSegDistSq(p, a, b) with the same arguments, so its result is +// <= this edge's value; the margin covers FP contraction differences between call sites. +inline double edgeDistBound(const Point& p, const Cell& near, const Polygon& polygon) { + const auto& ring = near.kr < 0 ? polygon.outer() : polygon.inners()[near.kr]; + const Point& a = ring[near.ki]; + const Point& b = ring[near.ki == 0 ? ring.size() - 1 : near.ki - 1]; + return std::sqrt(getSegDistSq(p, a, b)) * (1 + 1e-12) + 1e-11; +} + +#ifdef TMOPT_VERIFY +// unmodified original algorithm, for cross-checking +Point polylabelOrig(const Polygon& polygon, double precision) { + Box envelope; + geom::envelope(polygon.outer(), envelope); + const Point size { + envelope.max_corner().get<0>() - envelope.min_corner().get<0>(), + envelope.max_corner().get<1>() - envelope.min_corner().get<1>() + }; + const double cellSize = std::min(size.get<0>(), size.get<1>()); + double h = cellSize / 2; + auto compareMax = [] (const Cell& a, const Cell& b) { + if (a.max != b.max) return a.max < b.max; + if (a.d != b.d) return a.d < b.d; + if (a.h != b.h) return a.h > b.h; + if (a.c.get<0>() != b.c.get<0>()) return a.c.get<0>() > b.c.get<0>(); + return a.c.get<1>() > b.c.get<1>(); + }; + using Queue = std::priority_queue, decltype(compareMax)>; + Queue cellQueue(compareMax); + if (cellSize == 0) return envelope.min_corner(); + for (double x = envelope.min_corner().get<0>(); x < envelope.max_corner().get<0>(); x += cellSize) + for (double y = envelope.min_corner().get<1>(); y < envelope.max_corner().get<1>(); y += cellSize) + cellQueue.push(Cell({x + h, y + h}, h, polygon)); + auto bestCell = getCentroidCell(polygon); + Cell bboxCell( + Point { + envelope.min_corner().get<0>() + size.get<0>() / 2.0, + envelope.min_corner().get<1>() + size.get<1>() / 2.0 + }, 0, polygon); + if (bboxCell.d > bestCell.d) bestCell = bboxCell; + while (!cellQueue.empty()) { + auto cell = cellQueue.top(); + cellQueue.pop(); + if (cell.d > bestCell.d) bestCell = cell; + if (cell.max - bestCell.d <= precision) continue; + h = cell.h / 2; + cellQueue.push(Cell({cell.c.get<0>() - h, cell.c.get<1>() - h}, h, polygon)); + cellQueue.push(Cell({cell.c.get<0>() + h, cell.c.get<1>() - h}, h, polygon)); + cellQueue.push(Cell({cell.c.get<0>() - h, cell.c.get<1>() + h}, h, polygon)); + cellQueue.push(Cell({cell.c.get<0>() + h, cell.c.get<1>() + h}, h, polygon)); + } + return bestCell.c; +} +#endif + } // namespace detail Point polylabel(const Polygon& polygon, double precision = 0.00001, bool debug = false) { @@ -156,23 +436,25 @@ Point polylabel(const Polygon& polygon, double precision = 0.00001, bool debug = if (a.c.get<0>() != b.c.get<0>()) return a.c.get<0>() > b.c.get<0>(); return a.c.get<1>() > b.c.get<1>(); }; - using Queue = std::priority_queue, decltype(compareMax)>; - Queue cellQueue(compareMax); if (cellSize == 0) { return envelope.min_corner(); } - // cover polygon with initial cells - - for (double x = envelope.min_corner().get<0>(); x < envelope.max_corner().get<0>(); x += cellSize) { - for (double y = envelope.min_corner().get<1>(); y < envelope.max_corner().get<1>(); y += cellSize) { - cellQueue.push(Cell({x + h, y + h}, h, polygon)); - } - } +#ifdef POLYLABEL_SSE2 + // flatten the polygon's edges once for all distance evaluations of this call + thread_local PolygonEdges edges; + edges.build(polygon); + activeEdges = &edges; + struct EdgesGuard { + ~EdgesGuard() { activeEdges = nullptr; edges.release(); } + } edgesGuard; +#endif // take centroid as the first best guess - auto bestCell = getCentroidCell(polygon); + // (computed before the initial grid so that grid cells can be pruned; Cell construction is pure) + const Cell centroidCell = getCentroidCell(polygon); + auto bestCell = centroidCell; // second guess: bounding box centroid Cell bboxCell( @@ -184,11 +466,56 @@ Point polylabel(const Polygon& polygon, double precision = 0.00001, bool debug = bestCell = bboxCell; } - auto numProbes = cellQueue.size(); + // max-heap of cells, same algorithms as std::priority_queue (identical pop sequence); + // reused per thread to avoid reallocation (polylabel does not recurse) + thread_local std::vector cellQueue; + cellQueue.clear(); + auto pushCell = [&](const Point& p, double hc) { + cellQueue.push_back(Cell(p, hc, polygon)); + std::push_heap(cellQueue.begin(), cellQueue.end(), compareMax); + }; + + // A cell whose distance is bounded by ub can neither replace bestCell (needs d > bestCell.d) + // nor be subdivided when popped (cell.max - bestCell.d <= precision), as bestCell.d only grows. + // Such a cell is a no-op; dropping it leaves the pop order of the others unchanged. + const double sqrt2 = std::sqrt(2); + auto provablyNoOp = [&](const Point& p, double ub, double hc) { + if (!(ub <= bestCell.d && (ub + hc * sqrt2) - bestCell.d <= precision)) return false; +#ifdef TMOPT_VERIFY + Cell full(p, hc, polygon); + if (full.d > bestCell.d || full.max - bestCell.d > precision || full.d > ub) { + std::cerr << "TMOPT_VERIFY polylabel: pruned cell is not a no-op" << std::endl; + std::abort(); + } + static struct PruneCount { + std::atomic n{0}; + ~PruneCount() { std::cerr << "TMOPT_VERIFY polylabel pruned cells: " << n.load() << std::endl; } + } pruneCount; + pruneCount.n++; +#else + (void)p; +#endif + return true; + }; + + // cover polygon with initial cells + + std::size_t numProbes = 0; + for (double x = envelope.min_corner().get<0>(); x < envelope.max_corner().get<0>(); x += cellSize) { + for (double y = envelope.min_corner().get<1>(); y < envelope.max_corner().get<1>(); y += cellSize) { + numProbes++; + const Point p {x + h, y + h}; + const double ub = std::min(edgeDistBound(p, bboxCell, polygon), edgeDistBound(p, centroidCell, polygon)); + if (provablyNoOp(p, ub, h)) continue; + pushCell(p, h); + } + } + while (!cellQueue.empty()) { // pick the most promising cell from the queue - auto cell = cellQueue.top(); - cellQueue.pop(); + std::pop_heap(cellQueue.begin(), cellQueue.end(), compareMax); + auto cell = cellQueue.back(); + cellQueue.pop_back(); // update the best cell if we found a better one if (cell.d > bestCell.d) { @@ -200,11 +527,18 @@ Point polylabel(const Polygon& polygon, double precision = 0.00001, bool debug = if (cell.max - bestCell.d <= precision) continue; // split the cell into four cells + // (children whose distance is bounded by the parent's nearest edge may be provable no-ops) h = cell.h / 2; - cellQueue.push(Cell({cell.c.get<0>() - h, cell.c.get<1>() - h}, h, polygon)); - cellQueue.push(Cell({cell.c.get<0>() + h, cell.c.get<1>() - h}, h, polygon)); - cellQueue.push(Cell({cell.c.get<0>() - h, cell.c.get<1>() + h}, h, polygon)); - cellQueue.push(Cell({cell.c.get<0>() + h, cell.c.get<1>() + h}, h, polygon)); + const Point children[4] = { + {cell.c.get<0>() - h, cell.c.get<1>() - h}, + {cell.c.get<0>() + h, cell.c.get<1>() - h}, + {cell.c.get<0>() - h, cell.c.get<1>() + h}, + {cell.c.get<0>() + h, cell.c.get<1>() + h} + }; + for (const Point& p : children) { + if (provablyNoOp(p, edgeDistBound(p, cell, polygon), h)) continue; + pushCell(p, h); + } numProbes += 4; } @@ -213,6 +547,24 @@ Point polylabel(const Polygon& polygon, double precision = 0.00001, bool debug = std::cout << "best distance: " << bestCell.d << std::endl; } + // don't keep a huge buffer alive after an exceptionally large polygon + if (cellQueue.capacity() > (1u << 12)) std::vector().swap(cellQueue); + +#ifdef TMOPT_VERIFY + { +#ifdef POLYLABEL_SSE2 + activeEdges = nullptr; // polylabelOrig must use pointToPolygonDistKOrig +#endif + const Point ref = polylabelOrig(polygon, precision); + const double r[2] = { ref.get<0>(), ref.get<1>() }; + const double g[2] = { bestCell.c.get<0>(), bestCell.c.get<1>() }; + if (std::memcmp(r, g, sizeof(r)) != 0) { + std::cerr << "TMOPT_VERIFY polylabel: result differs from original" << std::endl; + std::abort(); + } + } +#endif + return bestCell.c; } diff --git a/include/sorted_node_store.h b/include/sorted_node_store.h index 61fdfad3..a4b5e896 100644 --- a/include/sorted_node_store.h +++ b/include/sorted_node_store.h @@ -63,6 +63,7 @@ class SortedNodeStore : public NodeStore void reopen() override; void finalize(size_t threadNum) override; LatpLon at(NodeID i) const override; + void at(const NodeID* ids, size_t n, LatpLon* out) const override; size_t size() const override; void batchStart() override; void insert(const std::vector& elements) override; @@ -80,6 +81,10 @@ class SortedNodeStore : public NodeStore // chunk is sufficiently large. bool compressNodes; + // Unique across all stores and reopen() calls; tags the thread-local + // decoded-chunk cache so it never serves data from a previous store. + uint64_t generation; + mutable std::mutex orphanageMutex; std::vector groups; std::vector> allocatedMemory; diff --git a/include/tag_map.h b/include/tag_map.h index c9510445..4f7df1d6 100644 --- a/include/tag_map.h +++ b/include/tag_map.h @@ -1,6 +1,7 @@ #ifndef _TAG_MAP_H #define _TAG_MAP_H +#include #include #include #include @@ -61,15 +62,26 @@ class TagMap { Iterator end() const; private: - uint32_t ensureString( - std::vector>& vector, - const protozero::data_view& value - ); + static constexpr size_t SHARDS = 16; + static_assert((SHARDS & (SHARDS - 1)) == 0, "SHARDS must be a power of two"); + + // Keys are stored flat (pointer + size inline) with their valueLoc, so a + // lookup miss can reject a candidate by size without dereferencing it. + struct KeySlot { + const char* data; + uint32_t size; + uint32_t valueLoc; + }; + + uint32_t ensureKey(const protozero::data_view& key); + uint32_t ensureValue(const protozero::data_view& value); + std::array, SHARDS> keys; + std::array, SHARDS> values; - std::vector> keys; - std::vector> key2value; - std::vector> values; + // Bit i set => keys[i] (resp. values[i]) may be non-empty. + uint16_t usedKeyShards = 0; + uint16_t usedValueShards = 0; }; #endif //_TAG_MAP_H diff --git a/include/tile_coordinates_set.h b/include/tile_coordinates_set.h index acc7a45b..3d45c50e 100644 --- a/include/tile_coordinates_set.h +++ b/include/tile_coordinates_set.h @@ -2,6 +2,9 @@ #define TILE_COORDINATES_SET_H #include +#include +#include +#include #include "coordinates.h" // Interface representing a bitmap of tiles of interest at a given zoom. @@ -11,6 +14,8 @@ class TileCoordinatesSet { virtual void set(TileCoordinate x, TileCoordinate y) = 0; virtual size_t size() const = 0; virtual size_t zoom() const = 0; + // Visit every set tile in x-major, y-minor order. + virtual void forEachSet(const std::function& f) const; }; // Read-write implementation for precise sets; maximum zoom is @@ -18,14 +23,22 @@ class TileCoordinatesSet { class PreciseTileCoordinatesSet : public TileCoordinatesSet { public: PreciseTileCoordinatesSet(unsigned int zoom); + ~PreciseTileCoordinatesSet(); + PreciseTileCoordinatesSet(const PreciseTileCoordinatesSet&) = delete; + PreciseTileCoordinatesSet& operator=(const PreciseTileCoordinatesSet&) = delete; bool test(TileCoordinate x, TileCoordinate y) const override; size_t size() const override; size_t zoom() const override; void set(TileCoordinate x, TileCoordinate y) override; + void forEachSet(const std::function& f) const override; private: unsigned int zoom_; - std::vector tiles; + uint64_t numBits; + size_t numWords; + uint64_t* tiles; // calloc'd, numWords words + // 1 bit per block of 64 words of `tiles`; set when the block may hold a set bit. + std::vector summary; }; // Read-only implementation for a lossy set. Used when zoom is diff --git a/include/tile_data.h b/include/tile_data.h index cd85b608..5bc9fa67 100644 --- a/include/tile_data.h +++ b/include/tile_data.h @@ -378,6 +378,9 @@ class TileDataSource { } virtual void populateMultiPolygon(MultiPolygon& dst, NodeID objectID); + // True if the polygon populateMultiPolygon returns for objectID is known to + // have passed geom::is_valid on exactly the same coordinates. + virtual bool polygonKnownValid(NodeID objectID) const { return false; } unsigned int getIndexZoom() const { return indexZoom; } inline size_t getId(NodeID id) const { diff --git a/include/visvalingam.h b/include/visvalingam.h index 7bc29da8..eb43d5bc 100644 --- a/include/visvalingam.h +++ b/include/visvalingam.h @@ -6,5 +6,6 @@ Linestring simplifyVis(const Linestring &ls, double max_distance); Polygon simplifyVis(const Polygon &p, double max_distance); MultiPolygon simplifyVis(const MultiPolygon &mp, double max_distance); +MultiPolygon simplifyVis(const MultiPolygon &mp, double max_distance, bool &allValid); #endif //_VISVALINGAM_H diff --git a/src/attribute_store.cpp b/src/attribute_store.cpp index 1f1f215c..c81c5027 100644 --- a/src/attribute_store.cpp +++ b/src/attribute_store.cpp @@ -2,6 +2,7 @@ #include #include +#include // AttributeKeyStore thread_local std::map tlsKeys2Index; @@ -100,66 +101,175 @@ const AttributePair& AttributePairStore::getPairUnsafe(uint32_t i) const { thread_local uint64_t tlsPairLookups = 0; thread_local uint64_t tlsPairLookupsUncached = 0; -thread_local std::vector cachedAttributePairPointers(256); -thread_local std::vector cachedAttributePairIndexes(256); -uint32_t AttributePairStore::addPair(AttributePair& pair, bool isHot) { - if (isHot) { - { - // First, check our thread-local map. - const auto& index = tlsHotShard.find(pair); - if (index != -1) - return index; +// Thread-local open-addressing index of pairs we have already looked up or +// added: maps (pair, isHot) to the index the global store returned. Entries +// point at elements of pairs[shard], which a deque never moves, and each +// pointer was taken under the shard lock, so it can be dereferenced without +// a lock. The index for a given (pair, isHot) never changes (the hot shard +// only grows, and once full it stays full), so cached results are exact. +namespace { +struct PairIndexEntry { + const AttributePair* pair; + uint32_t tag; + uint32_t index; +}; + +const uint32_t PAIR_INDEX_MIN_SLOTS = 1 << 14; +const uint32_t PAIR_INDEX_MAX_SLOTS = 1 << 16; +const uint32_t PAIR_INDEX_PROBES = 8; + +// Bumped by AttributeStore::reset() to invalidate every thread's table. +std::atomic pairIndexGeneration(0); + +struct PairIndexTable { + std::vector slots; + uint32_t mask = 0; + uint32_t count = 0; + uint32_t victim = 0; + uint64_t generation = 0; + const void* owner = nullptr; + + static uint32_t home(uint32_t tag, uint32_t mask) { + return (tag * 0x9E3779B1u) >> 16 & mask; + } + + void ensure(const void* store) { + const uint64_t gen = pairIndexGeneration.load(std::memory_order_acquire); + if (owner == store && generation == gen && !slots.empty()) + return; + slots.assign(PAIR_INDEX_MIN_SLOTS, PairIndexEntry{nullptr, 0, 0}); + mask = PAIR_INDEX_MIN_SLOTS - 1; + count = 0; + owner = store; + generation = gen; + } + + bool find(const AttributePair& pair, uint32_t tag, uint32_t& index) const { + uint32_t slot = home(tag, mask); + for (uint32_t i = 0; i < PAIR_INDEX_PROBES; i++) { + const PairIndexEntry& e = slots[(slot + i) & mask]; + if (e.pair == nullptr) + return false; + if (e.tag == tag && *e.pair == pair) { + index = e.index; + return true; + } } + return false; + } - // Not found, ensure our local map is up-to-date for future calls, - // and fall through to the main map. - if (!tlsHotShard.full()) { - std::lock_guard lock(pairsMutex[0]); - tlsHotShard = pairs[0]; + void place(const PairIndexEntry& entry) { + uint32_t slot = home(entry.tag, mask); + for (uint32_t i = 0; i < PAIR_INDEX_PROBES; i++) { + PairIndexEntry& e = slots[(slot + i) & mask]; + if (e.pair == nullptr) { + e = entry; + count++; + return; + } + } + // Window full: overwrite an entry (never creates an empty slot, so + // lookups that stop at an empty slot stay correct). + slots[(slot + (victim++ % PAIR_INDEX_PROBES)) & mask] = entry; + } + + void insert(const AttributePair* pair, uint32_t tag, uint32_t index) { + if (slots.size() < PAIR_INDEX_MAX_SLOTS && (count + 1) * 2 > slots.size()) { + std::vector old; + old.swap(slots); + slots.assign(old.size() * 2, PairIndexEntry{nullptr, 0, 0}); + mask = slots.size() - 1; + count = 0; + for (const auto& e : old) + if (e.pair != nullptr) + place(e); + } + place(PairIndexEntry{pair, tag, index}); + } +}; + +thread_local PairIndexTable tlsPairIndex; + +inline size_t coldShardForHash(size_t hash) { + size_t shard = hash % ATTRIBUTE_SHARDS; + // Shard 0 is for hot pairs -- pick another shard if it gets selected. + if (shard == 0) shard = (hash >> 8) % ATTRIBUTE_SHARDS; + if (shard == 0) shard = (hash >> 16) % ATTRIBUTE_SHARDS; + if (shard == 0) shard = (hash >> 24) % ATTRIBUTE_SHARDS; + if (shard == 0) shard = 1; + return shard; +} +} + +uint32_t AttributePairStore::addPair(AttributePair& pair, bool isHot) { + const size_t hash = pair.hash(); + // The low bit of the tag records which path (hot/cold) produced the index. + const uint32_t tag = ((uint32_t)((uint64_t)hash ^ ((uint64_t)hash >> 32)) & ~1u) | (isHot ? 1u : 0u); + + PairIndexTable& table = tlsPairIndex; + table.ensure(this); + { + uint32_t cached; + if (table.find(pair, tag, cached)) { + if (cached >> (32 - SHARD_BITS) != 0) { + tlsPairLookups++; + if (tlsPairLookups % 1024 == 0) + lookups += 1024; + } +#ifdef TMOPT_VERIFY + { + int64_t slow = -1; + bool done = false; + if (isHot) { + std::lock_guard lock(pairsMutex[0]); + const auto index = pairs[0].find(pair); + if (index != -1) { slow = index; done = true; } + else if (!pairs[0].full()) { slow = -2; done = true; } + } + if (!done) { + const size_t shard = coldShardForHash(hash); + std::lock_guard lock(pairsMutex[shard]); + const auto index = pairs[shard].find(pair); + if (index != -1) + slow = (shard << (32 - SHARD_BITS)) + index; + } + if (slow != (int64_t)cached) { + std::cerr << "TMOPT_VERIFY: AttributePairStore::addPair mismatch " << cached << " vs " << slow << std::endl; + std::abort(); + } + } +#endif + return cached; } + } + if (isHot) { // This might be a popular pair, worth re-using. // Have we already assigned it a hot ID? std::lock_guard lock(pairsMutex[0]); const auto& index = pairs[0].find(pair); - if (index != -1) + if (index != -1) { + table.insert(&pairs[0][index], tag, index); return index; + } if (!pairs[0].full()) { pair.ensureStringIsOwned(); uint32_t offset = pairs[0].add(pair); uint32_t rv = (0 << (32 - SHARD_BITS)) + offset; + table.insert(&pairs[0][offset], tag, rv); return rv; } } // This is either not a hot key, or there's no room for in the hot shard. // Throw it on the pile with the rest of the pairs. - size_t hash = pair.hash(); - - const size_t candidateIndex = hash % cachedAttributePairPointers.size(); - // Before taking a lock, see if we've seen this attribute pair recently. - tlsPairLookups++; if (tlsPairLookups % 1024 == 0) { lookups += 1024; } - - { - const AttributePair* candidate = cachedAttributePairPointers[candidateIndex]; - - if (candidate != nullptr && *candidate == pair) - return cachedAttributePairIndexes[candidateIndex]; - } - - - size_t shard = hash % ATTRIBUTE_SHARDS; - // Shard 0 is for hot pairs -- pick another shard if it gets selected. - if (shard == 0) shard = (hash >> 8) % ATTRIBUTE_SHARDS; - if (shard == 0) shard = (hash >> 16) % ATTRIBUTE_SHARDS; - if (shard == 0) shard = (hash >> 24) % ATTRIBUTE_SHARDS; - if (shard == 0) shard = 1; + const size_t shard = coldShardForHash(hash); std::lock_guard lock(pairsMutex[shard]); @@ -170,9 +280,7 @@ uint32_t AttributePairStore::addPair(AttributePair& pair, bool isHot) { const auto& index = pairs[shard].find(pair); if (index != -1) { const uint32_t rv = (shard << (32 - SHARD_BITS)) + index; - cachedAttributePairPointers[candidateIndex] = &pairs[shard][index]; - cachedAttributePairIndexes[candidateIndex] = rv; - + table.insert(&pairs[shard][index], tag, rv); return rv; } @@ -183,6 +291,7 @@ uint32_t AttributePairStore::addPair(AttributePair& pair, bool isHot) { throw std::out_of_range("pair shard overflow"); uint32_t rv = (shard << (32 - SHARD_BITS)) + offset; + table.insert(&pairs[shard][offset], tag, rv); return rv; }; @@ -259,6 +368,27 @@ void AttributeStore::addAttribute(AttributeSet& attributeSet, std::string const bool isHot = v >= 0 && v <= 25; attributeSet.addPair(pairStore.addPair(kv, isHot)); } +void AttributeStore::addAttributeWithKeyIndex(AttributeSet& attributeSet, uint16_t keyIndex, std::string const &key, const protozero::data_view v, char minzoom) { + PooledString ps(&v); + AttributePair kv(keyIndex, ps, minzoom); + bool isHot = AttributePair::isHot(key, v); + attributeSet.addPair(pairStore.addPair(kv, isHot)); +} +void AttributeStore::addAttributeWithKeyIndex(AttributeSet& attributeSet, uint16_t keyIndex, bool v, char minzoom) { + AttributePair kv(keyIndex,v,minzoom); + bool isHot = true; // All bools are eligible to be hot pairs + attributeSet.addPair(pairStore.addPair(kv, isHot)); +} +void AttributeStore::addAttributeWithKeyIndex(AttributeSet& attributeSet, uint16_t keyIndex, double v, char minzoom) { + AttributePair kv(keyIndex,v,minzoom); + bool isHot = v >= 0 && v <= 25 && ceil(v) == v; // Whole numbers in 0..25 are eligible to be hot pairs + attributeSet.addPair(pairStore.addPair(kv, isHot)); +} +void AttributeStore::addAttributeWithKeyIndex(AttributeSet& attributeSet, uint16_t keyIndex, int v, char minzoom) { + AttributePair kv(keyIndex,v,minzoom); + bool isHot = v >= 0 && v <= 25; + attributeSet.addPair(pairStore.addPair(kv, isHot)); +} void AttributeSet::finalize() { // Ensure that values are sorted, giving us a canonical representation, @@ -447,8 +577,7 @@ void AttributeStore::reset() { for (int i = 0; i < cachedAttributeSetPointers.size(); i++) cachedAttributeSetPointers[i] = nullptr; - for (int i = 0; i < cachedAttributePairPointers.size(); i++) - cachedAttributePairPointers[i] = nullptr; + pairIndexGeneration++; } void AttributeStore::finalize() { diff --git a/src/coordinates.cpp b/src/coordinates.cpp index 3e8ea2f2..4a79f83e 100644 --- a/src/coordinates.cpp +++ b/src/coordinates.cpp @@ -65,3 +65,32 @@ void fillCoveredTiles(std::unordered_set& tileSet) { prevX = tileX, prevY = tileY; } } + +void sortUniqueTiles(std::vector &tiles, size_t start) { + if (tiles.size() - start > 1) { + std::sort(tiles.begin() + start, tiles.end(), TileCoordinatesCompare()); + tiles.erase(std::unique(tiles.begin() + start, tiles.end()), tiles.end()); + } +} + +void fillCoveredTiles(std::vector &tiles) { + // Fast exit when no column has a y gap (the common case) + size_t n = tiles.size(), i = 1; + while (i < n && !(tiles[i].x == tiles[i-1].x && tiles[i].y > tiles[i-1].y + 1)) i++; + if (i >= n) return; + + // Rebuild in order: copy the prefix, then each tile preceded by the y-gap tiles of its column + thread_local std::vector filled; + filled.clear(); + filled.insert(filled.end(), tiles.begin(), tiles.begin() + i); + for (; i < n; i++) { + const TileCoordinates prev = tiles[i-1], cur = tiles[i]; + if (cur.x == prev.x) { + for (TileCoordinate fillY = prev.y+1; fillY < cur.y; fillY++) + filled.push_back(TileCoordinates(cur.x, fillY)); + } + filled.push_back(cur); + } + tiles.swap(filled); + if (filled.capacity() > (1u << 16)) std::vector().swap(filled); +} diff --git a/src/coordinates_geom.cpp b/src/coordinates_geom.cpp index 05180d84..16887856 100644 --- a/src/coordinates_geom.cpp +++ b/src/coordinates_geom.cpp @@ -199,3 +199,99 @@ void insertIntermediateTiles(Ring const &points, uint baseZoom, std::unordered_s { impl_insertIntermediateTiles(points, baseZoom, tileSet); } + +// Vector variant: same tiles as above, without hashing. Each vertex tile is computed once +// (the previous vertex tile is carried over), and lon2tilexf/latp2tileyf are inlined: +// multiplying by the exact power of two 2^z gives the same bits as scalbn(v, z). +static inline void pushTile(std::vector &tiles, TileCoordinate x, TileCoordinate y) { + if (tiles.empty() || tiles.back().x != x || tiles.back().y != y) + tiles.push_back(TileCoordinates(x, y)); +} + +template +void impl_insertIntermediateTilesVec(T const &points, uint baseZoom, std::vector &tiles) { + const size_t start = tiles.size(); + const double scale = scalbn(1.0, (int)(uint8_t)baseZoom); + TileCoordinate tileX1 = 0, tileY1 = 0; + for (auto it = points.begin(); it != points.end(); ++it) { + // Calculate p2 tile (same arithmetic as lon2tilexf/latp2tileyf) + double p2y = std::min(MaxLatp, std::max(MinLatp, it->y())); + double tileXf2 = ((it->x()+180.0) * (1/360.0)) * scale, tileYf2 = ((180.0-p2y) * (1/360.0)) * scale; + TileCoordinate tileX2 = static_cast(tileXf2), tileY2 = static_cast(tileYf2); + if (it == points.begin()) { // first point, so no line + pushTile(tiles, tileX2, tileY2); + tileX1 = tileX2; tileY1 = tileY2; + continue; + } + + // Supercover line algorithm from http://eugen.dedu.free.fr/projects/bresenham/ + int i; // loop counter + int ystep, xstep; // the step on y and x axis + int error; // the error accumulated during the increment + int errorprev; // *vision the previous value of the error variable + int y = tileY1, x = tileX1; // the line points + int ddy, ddx; // compulsory variables: the double values of dy and dx + int dx = tileX2 - tileX1; + int dy = tileY2 - tileY1; + + if (dy < 0) { ystep = -1; dy = -dy; } else { ystep = 1; } + if (dx < 0) { xstep = -1; dx = -dx; } else { xstep = 1; } + + ddy = 2 * dy; // work with double values for full precision + ddx = 2 * dx; + if (ddx >= ddy) { // first octant (0 <= slope <= 1) + errorprev = error = dx; // start in the middle of the square + for (i=0 ; i < dx ; i++) { // do not use the first point (already done) + x += xstep; + error += ddy; + if (error > ddx){ // increment y if AFTER the middle ( > ) + y += ystep; + error -= ddx; + if (error + errorprev < ddx) // bottom square also + pushTile(tiles, x, y-ystep); + else if (error + errorprev > ddx) // left square also + pushTile(tiles, x-xstep, y); + else { // corner: bottom and left squares also + pushTile(tiles, x, y-ystep); + pushTile(tiles, x-xstep, y); + } + } + pushTile(tiles, x, y); + errorprev = error; + } + } else { // the same as above + errorprev = error = dy; + for (i=0 ; i < dy ; i++){ + y += ystep; + error += ddx; + if (error > ddy){ + x += xstep; + error -= ddy; + if (error + errorprev < ddy) + pushTile(tiles, x-xstep, y); + else if (error + errorprev > ddy) + pushTile(tiles, x, y-ystep); + else{ + pushTile(tiles, x-xstep, y); + pushTile(tiles, x, y-ystep); + } + } + pushTile(tiles, x, y); + errorprev = error; + } + } + pushTile(tiles, tileX2, tileY2); + tileX1 = tileX2; tileY1 = tileY2; + } + sortUniqueTiles(tiles, start); +} + +void insertIntermediateTiles(Linestring const &points, uint baseZoom, std::vector &tiles) +{ + impl_insertIntermediateTilesVec(points, baseZoom, tiles); +} + +void insertIntermediateTiles(Ring const &points, uint baseZoom, std::vector &tiles) +{ + impl_insertIntermediateTilesVec(points, baseZoom, tiles); +} diff --git a/src/geom.cpp b/src/geom.cpp index 5d135e72..406905e9 100644 --- a/src/geom.cpp +++ b/src/geom.cpp @@ -8,6 +8,15 @@ #include "geometry/correct.hpp" +#include +#include +#include +#ifdef TMOPT_VERIFY +#include +#include +#include +#endif + #include #include @@ -146,11 +155,348 @@ void make_valid(MultiPolygon &mp) mp = result; } +// --------------- +// Conservative validity precheck (see geom.h). +// Every decision is taken only when it is clearly away from the degenerate +// case. An orientation determinant counts as non-zero only if it exceeds +// (a) a relative bound on its own rounding error, (b) the point-to-line +// distance tolerance `tol` (relative to the coordinate magnitude, far above +// Boost's epsilon point equality) and (c) Boost's side_by_triangle zero band +// (epsilon * max(1, |coordinate differences|), whichever point it permutes +// first) with a wide margin. Anything borderline answers "unsure". + +namespace { + +const double kPrecheckEps = 1e-9; +const std::size_t kPrecheckMaxHoles = 256; +const std::size_t kPrecheckRetained = 1 << 16; + +struct PrecheckSeg { + double minX, maxX, minY, maxY; + const Point *a; // segment a[0] -> a[1] + uint32_t ring, idx, count; // ring number, segment index, segments in the ring +}; +struct PrecheckBox { + double minX, maxX, minY, maxY; +}; +struct PrecheckScratch { + std::vector segs; + std::vector boxes; + void trim() { + if (segs.capacity() > kPrecheckRetained) std::vector().swap(segs); + if (boxes.capacity() > kPrecheckRetained) std::vector().swap(boxes); + } +}; +thread_local PrecheckScratch precheckScratch; + +inline double max3(double a, double b, double c) { return std::max(a, std::max(b, c)); } + +// Orientation of c relative to a->b: +1 counter-clockwise, -1 clockwise, 0 when too close to call. +inline int clearOrient(const Point &a, const Point &b, const Point &c, double tol) { + const double dx = b.x() - a.x(), dy = b.y() - a.y(); + const double px = c.x() - a.x(), py = c.y() - a.y(); + const double qx = c.x() - b.x(), qy = c.y() - b.y(); + const double t1 = dx * py, t2 = dy * px; + const double d = t1 - t2; + const double f = std::max(max3(1.0, std::fabs(dx), std::fabs(dy)), + std::max(max3(std::fabs(px), std::fabs(py), std::fabs(qx)), std::fabs(qy))); + const double bound = max3(kPrecheckEps * (std::fabs(t1) + std::fabs(t2)), + tol * (std::fabs(dx) + std::fabs(dy)), + 64 * std::numeric_limits::epsilon() * f * f); + if (std::fabs(d) > bound) return d > 0 ? 1 : -1; + return 0; +} + +// Coordinates finite; accumulates the largest magnitude. +inline bool ringFinite(const Ring &r, double &maxAbs) { + for (const Point &p : r) { + const double x = p.x(), y = p.y(); + if (!std::isfinite(x) || !std::isfinite(y)) return false; + maxAbs = max3(maxAbs, std::fabs(x), std::fabs(y)); + } + return true; +} + +// Single-ring tests: closed, >= 4 points, consecutive points either exactly +// equal (Boost allows duplicate points; such zero-length segments are skipped) +// or clearly distinct, area clearly non-zero with the orientation Boost +// expects, no spikes. Appends the ring's segments and returns its bbox. +bool ringClear(const Ring &r, bool interior, double tol, uint32_t ringIdx, std::vector &segs, PrecheckBox &box) { + const std::size_t n = r.size(); + if (n < 4 || n >= std::numeric_limits::max()) return false; + const std::size_t m = n - 1; // segments: i runs r[i] -> r[i+1] + if (r[0].x() != r[m].x() || r[0].y() != r[m].y()) return false; + + const std::size_t first = segs.size(); + box.minX = box.minY = std::numeric_limits::infinity(); + box.maxX = box.maxY = -std::numeric_limits::infinity(); + for (std::size_t i = 0; i < m; i++) { + const Point &a = r[i], &b = r[i + 1]; + if (a.x() == b.x() && a.y() == b.y()) continue; + if (std::fabs(b.x() - a.x()) <= tol && std::fabs(b.y() - a.y()) <= tol) return false; + PrecheckSeg s; + s.minX = std::min(a.x(), b.x()); s.maxX = std::max(a.x(), b.x()); + s.minY = std::min(a.y(), b.y()); s.maxY = std::max(a.y(), b.y()); + s.a = &a; + s.ring = ringIdx; s.idx = static_cast(segs.size() - first); + segs.push_back(s); + box.minX = std::min(box.minX, s.minX); box.maxX = std::max(box.maxX, s.maxX); + box.minY = std::min(box.minY, s.minY); box.maxY = std::max(box.maxY, s.maxY); + } + const std::size_t count = segs.size() - first; + if (count < 3) return false; + for (std::size_t k = first; k < segs.size(); k++) segs[k].count = static_cast(count); + + // Signed area (doubled) relative to r[0], positive = counter-clockwise. It must + // clearly exceed both our rounding error and that of Boost's (x1+x2)*(y1-y2) sum. + const double x0 = r[0].x(), y0 = r[0].y(); + double area = 0, areaScale = 0, boostScale = 0; + for (std::size_t i = 1; i + 1 < m; i++) { + const double t1 = (r[i].x() - x0) * (r[i + 1].y() - y0); + const double t2 = (r[i].y() - y0) * (r[i + 1].x() - x0); + area += t1 - t2; + areaScale += std::fabs(t1) + std::fabs(t2); + } + for (std::size_t i = 0; i < m; i++) + boostScale += (std::fabs(r[i].x()) + std::fabs(r[i + 1].x())) * std::fabs(r[i].y() - r[i + 1].y()); + const double relEps = std::max(kPrecheckEps, 8.0 * (n + 4) * std::numeric_limits::epsilon()); + if (!(std::fabs(area) > relEps * areaScale && std::fabs(area) > relEps * boostScale)) return false; + const bool clockwise = boost::geometry::point_order::value == boost::geometry::clockwise; + const bool wantPositive = clockwise == interior; // clockwise polygon: CW outer, CCW inners + if (wantPositive ? area <= 0 : area >= 0) return false; + + // Consecutive segments a->b, b->c only share b unless they fold back (spike). + for (std::size_t k = 0; k < count; k++) { + const PrecheckSeg &prev = segs[first + (k == 0 ? count - 1 : k - 1)]; + const Point &a = prev.a[0]; + const Point &b = segs[first + k].a[0]; + const Point &c = segs[first + k].a[1]; + if (clearOrient(a, b, c, tol) != 0) continue; + // (nearly) collinear: must clearly continue forward + const double t1 = (b.x() - a.x()) * (c.x() - b.x()); + const double t2 = (b.y() - a.y()) * (c.y() - b.y()); + if (!(t1 + t2 > kPrecheckEps * (std::fabs(t1) + std::fabs(t2)))) return false; + } + return true; +} + +inline bool boxesApart(const PrecheckBox &a, const PrecheckBox &b, double tol) { + return a.maxX + tol < b.minX || b.maxX + tol < a.minX || a.maxY + tol < b.minY || b.maxY + tol < a.minY; +} + +// Two segments are clearly disjoint (adjacent segments of a ring are handled by the spike test). +inline bool segsClear(const PrecheckSeg &s, const PrecheckSeg &t, double tol) { + if (s.ring == t.ring) { + const uint32_t lo = std::min(s.idx, t.idx), hi = std::max(s.idx, t.idx); + if (hi - lo == 1 || (lo == 0 && hi == s.count - 1)) return true; + } + if (s.maxX + tol < t.minX || t.maxX + tol < s.minX || s.maxY + tol < t.minY || t.maxY + tol < s.minY) + return true; + const Point &p1 = s.a[0], &p2 = s.a[1], &q1 = t.a[0], &q2 = t.a[1]; + const int o1 = clearOrient(p1, p2, q1, tol), o2 = clearOrient(p1, p2, q2, tol); + if (o1 != 0 && o1 == o2) return true; // q entirely on one side of line p + const int o3 = clearOrient(q1, q2, p1, tol), o4 = clearOrient(q1, q2, p2, tol); + return o3 != 0 && o3 == o4; // p entirely on one side of line q; else crossing, touching or unsure +} + +// All pairs of segments clearly disjoint: sweep over segments sorted by minX. +// The answer does not depend on the order (every x-overlapping pair is tested). +bool segmentsClear(std::vector &segs, double tol) { + const std::size_t n = segs.size(); + if (n <= 32) { + for (std::size_t i = 0; i < n; i++) + for (std::size_t j = i + 1; j < n; j++) + if (!segsClear(segs[i], segs[j], tol)) return false; + return true; + } + std::sort(segs.begin(), segs.end(), [](const PrecheckSeg &a, const PrecheckSeg &b) { return a.minX < b.minX; }); + for (std::size_t i = 0; i < n; i++) { + const double limit = segs[i].maxX + tol; + for (std::size_t j = i + 1; j < n && segs[j].minX <= limit; j++) + if (!segsClear(segs[i], segs[j], tol)) return false; + } + return true; +} + +// Crossing-number test: 1 = h clearly inside r, -1 = clearly outside, 0 = unsure. +int pointClearlyInRing(const Point &h, const Ring &r, double tol) { + bool inside = false; + const double hx = h.x(), hy = h.y(); + for (std::size_t i = 0; i + 1 < r.size(); i++) { + const Point &a = r[i], &b = r[i + 1]; + if ((a.y() > hy + tol && b.y() > hy + tol) || (a.y() < hy - tol && b.y() < hy - tol)) continue; + if (std::fabs(a.y() - hy) <= tol || std::fabs(b.y() - hy) <= tol) return 0; + // a and b are clearly on opposite sides of the horizontal line through h + if (std::max(a.x(), b.x()) < hx - tol) continue; + if (std::min(a.x(), b.x()) > hx + tol) { inside = !inside; continue; } + const int o = clearOrient(a, b, h, tol); + if (o == 0) return 0; + if ((b.y() > a.y()) == (o > 0)) inside = !inside; // the ray to +x crosses a->b + } + return inside ? 1 : -1; +} + +// Location of a ring that is known not to touch r: tries a few of its vertices. +int ringClearlyInRing(const Ring &inner, const Ring &r, double tol) { + const std::size_t m = inner.size() - 1; + const std::size_t tries = std::min(m, 4); + for (std::size_t k = 0; k < tries; k++) { + const int loc = pointClearlyInRing(inner[k * m / tries], r, tol); + if (loc != 0) return loc; + } + return 0; +} + +bool polygonClear(const Polygon &p, PrecheckScratch &scratch, PrecheckBox &outerBox) { + const Ring &outer = p.outer(); + const auto &inners = p.inners(); + if (inners.size() > kPrecheckMaxHoles) return false; + double maxAbs = 0; + if (!ringFinite(outer, maxAbs)) return false; + for (const auto &inner : inners) + if (!ringFinite(inner, maxAbs)) return false; + const double tol = kPrecheckEps * std::max(1.0, maxAbs); + + std::vector &segs = scratch.segs; + segs.clear(); + PrecheckBox innerBoxes[kPrecheckMaxHoles]; + if (!ringClear(outer, false, tol, 0, segs, outerBox)) return false; + for (std::size_t i = 0; i < inners.size(); i++) + if (!ringClear(inners[i], true, tol, static_cast(i + 1), segs, innerBoxes[i])) return false; + + // No two rings touch, and no ring touches itself. + if (!segmentsClear(segs, tol)) return false; + + // Holes lie inside the outer ring and are not nested in each other. + for (std::size_t i = 0; i < inners.size(); i++) { + if (boxesApart(innerBoxes[i], outerBox, tol)) return false; + if (ringClearlyInRing(inners[i], outer, tol) != 1) return false; + } + for (std::size_t i = 0; i < inners.size(); i++) + for (std::size_t j = i + 1; j < inners.size(); j++) { + if (boxesApart(innerBoxes[i], innerBoxes[j], tol)) continue; + if (ringClearlyInRing(inners[i], inners[j], tol) != -1 || + ringClearlyInRing(inners[j], inners[i], tol) != -1) return false; + } + return true; +} + +bool multiPolygonClear(const MultiPolygon &mp, PrecheckScratch &scratch) { + if (mp.empty()) return false; + PrecheckBox box; + if (mp.size() == 1) return polygonClear(mp[0], scratch, box); + + // Polygons must have clearly disjoint bboxes (outer rings; holes lie inside them). + std::vector &boxes = scratch.boxes; + boxes.clear(); + double maxAbs = 0; + for (const auto &p : mp) { + if (p.outer().empty()) return false; + box.minX = box.minY = std::numeric_limits::infinity(); + box.maxX = box.maxY = -std::numeric_limits::infinity(); + for (const Point &pt : p.outer()) { + const double x = pt.x(), y = pt.y(); + if (!std::isfinite(x) || !std::isfinite(y)) return false; + box.minX = std::min(box.minX, x); box.maxX = std::max(box.maxX, x); + box.minY = std::min(box.minY, y); box.maxY = std::max(box.maxY, y); + } + maxAbs = std::max(maxAbs, max3(std::fabs(box.minX), std::fabs(box.maxX), std::max(std::fabs(box.minY), std::fabs(box.maxY)))); + boxes.push_back(box); + } + const double tol = kPrecheckEps * std::max(1.0, maxAbs); + std::sort(boxes.begin(), boxes.end(), [](const PrecheckBox &a, const PrecheckBox &b) { return a.minX < b.minX; }); + for (std::size_t i = 0; i < boxes.size(); i++) { + const double limit = boxes[i].maxX + tol; + for (std::size_t j = i + 1; j < boxes.size() && boxes[j].minX <= limit; j++) + if (!boxesApart(boxes[i], boxes[j], tol)) return false; + } + + for (const auto &p : mp) + if (!polygonClear(p, scratch, box)) return false; + return true; +} + +} // namespace + +bool polygonDefinitelyValid(const Polygon &p) { + PrecheckBox box; + const bool ok = polygonClear(p, precheckScratch, box); + precheckScratch.trim(); + return ok; +} + +bool multiPolygonDefinitelyValid(const MultiPolygon &mp) { + const bool ok = multiPolygonClear(mp, precheckScratch); + precheckScratch.trim(); + return ok; +} + +#ifdef TMOPT_VERIFY +static std::atomic precheckHits[PRECHECK_SITES], precheckCalls[PRECHECK_SITES]; +static const char *precheckSiteName(int site) { + static const char *names[PRECHECK_SITES] = { "tds_buildWayGeometry", "repair_one_polygon", "repair_multi_polygon", + "repair_multi_polygon_part", "CheckGeometry_relation", "CheckGeometry_way", "CorrectGeometry_reversed", "shp" }; + return names[site]; +} +static void precheckReport(int site) { + std::cerr << ("TMOPT_VERIFY: precheck " + std::string(precheckSiteName(site)) + " hits " + std::to_string(precheckHits[site].load()) + + " of " + std::to_string(precheckCalls[site].load()) + "\n"); +} +namespace { +struct PrecheckFinalReport { + ~PrecheckFinalReport() { + for (int s = 0; s < PRECHECK_SITES; s++) + if (precheckCalls[s].load()) precheckReport(s); + } +} precheckFinalReport; +} +void precheckVerifyCount(int site, bool hit) { + if (hit) precheckHits[site]++; + const uint64_t c = ++precheckCalls[site]; + if (c >= 1024 && (c & (c - 1)) == 0) precheckReport(site); +} +#endif + // Repair a single (possibly invalid) polygon in an area-preserving way and // append the resulting valid polygon(s) to `out`. Returns true on success. // `minArea` is the lower bound on the repaired area we are willing to accept. static bool repair_one_polygon(const Polygon &p, double minArea, MultiPolygon &out) { + // 0) Fast path: a ring without holes whose only fault is its orientation + // (e.g. after the y-flip into tile coordinates). On such a simple ring the + // dissolve finds no intersections and just reverses it, so do that directly. + if (p.inners().empty()) { + Polygon rev = p; + std::reverse(rev.outer().begin(), rev.outer().end()); + const double revArea = std::abs(geom::area(rev)); + bool revValid = false; + if (revArea > 1E-12 && revArea >= minArea) { + const bool pre = polygonDefinitelyValid(rev); +#ifdef TMOPT_VERIFY + precheckVerify(rev, pre, PRECHECK_REPAIR_ONE); +#endif + revValid = pre || geom::is_valid(rev); + } + if (revValid) { +#ifdef TMOPT_VERIFY + MultiPolygon fixed; + geometry::correct(p, fixed, 1E-12); + bool same = fixed.size() == 1 && fixed[0].inners().empty() && + fixed[0].outer().size() == rev.outer().size(); + for (std::size_t i = 0; same && i < rev.outer().size(); ++i) + same = fixed[0].outer()[i].x() == rev.outer()[i].x() && + fixed[0].outer()[i].y() == rev.outer()[i].y(); + if (!same) { + std::cerr << "TMOPT_VERIFY repair_one_polygon mismatch: fast " << boost::geometry::wkt(rev) + << " dissolve " << boost::geometry::wkt(fixed) << std::endl; + abort(); + } +#endif + out.push_back(std::move(rev)); + return true; + } + } + // 1) Dissolve (resolves self-intersections of this single polygon). try { MultiPolygon fixed; @@ -185,9 +531,15 @@ static bool repair_one_polygon(const Polygon &p, double minArea, MultiPolygon &o return false; } -bool repair_multi_polygon(MultiPolygon &mp) +bool repair_multi_polygon(MultiPolygon &mp, bool knownInvalid) { - if (geom::is_valid(mp)) return true; + if (!knownInvalid) { + const bool pre = multiPolygonDefinitelyValid(mp); +#ifdef TMOPT_VERIFY + precheckVerify(mp, pre, PRECHECK_REPAIR_MP); +#endif + if (pre || geom::is_valid(mp)) return true; + } // Repair PER POLYGON, area-preserving. Running make_valid/buffer on the whole // multipolygon can catastrophically COLLAPSE large/complex inputs: a clipped @@ -203,9 +555,17 @@ bool repair_multi_polygon(MultiPolygon &mp) MultiPolygon out; bool allValid = true; for (const auto &p : mp) { - if (geom::is_valid(p)) { - out.push_back(p); - continue; + // With a single polygon, the whole-multipolygon check above already + // established it is invalid. + if (mp.size() > 1) { + const bool pre = polygonDefinitelyValid(p); +#ifdef TMOPT_VERIFY + precheckVerify(p, pre, PRECHECK_REPAIR_MP_PART); +#endif + if (pre || geom::is_valid(p)) { + out.push_back(p); + continue; + } } // Lenient threshold: resolving a self-intersection legitimately changes a // single polygon's (shoelace) area, so anything down to half the original @@ -319,7 +679,9 @@ char bit_code(Point const &p, Box const &bbox) { } // Sutherland-Hodgeman polygon clipping algorithm -void fast_clip(Ring &points, Box const &bbox, Ring &result) { +// Returns true if any clipping was applied (false: points left untouched). +bool fast_clip(Ring &points, Box const &bbox, Ring &result) { + bool modified = false; // clip against each side of the clip rectangle for (char edge = 1; edge <= 8; edge *= 2) { bool needsClip = false; @@ -331,6 +693,7 @@ void fast_clip(Ring &points, Box const &bbox, Ring &result) { } if (!needsClip) continue; + modified = true; result.clear(); result.reserve(points.size() + 4); Point prev = points[points.size() - 1]; @@ -349,36 +712,43 @@ void fast_clip(Ring &points, Box const &bbox, Ring &result) { points.swap(result); if (points.size()==0) break; } + return modified; } -void fast_clip(Ring &points, Box const &bbox) { +bool fast_clip(Ring &points, Box const &bbox) { Ring result; - fast_clip(points, bbox, result); + return fast_clip(points, bbox, result); } // Wrappers for polygon/multipolygon -void fast_clip(Polygon &polygon, Box const &bbox) { +bool fast_clip(Polygon &polygon, Box const &bbox) { Ring result; - fast_clip(polygon.outer(), bbox, result); + bool modified = fast_clip(polygon.outer(), bbox, result); if (polygon.outer().empty()) { + modified = modified || !polygon.inners().empty(); polygon.inners().resize(0); - return; + return modified; } for (auto &inner: polygon.inners()) { - fast_clip(inner, bbox, result); + if (fast_clip(inner, bbox, result)) modified = true; } + size_t innersBefore = polygon.inners().size(); polygon.inners().erase(std::remove_if( polygon.inners().begin(), polygon.inners().end(), [](const Ring &inner) -> bool { return inner.empty(); }), polygon.inners().end()); + return modified || polygon.inners().size() != innersBefore; } -void fast_clip(MultiPolygon &mp, Box const &bbox) { +bool fast_clip(MultiPolygon &mp, Box const &bbox) { + bool modified = false; for (auto &polygon: mp) { - fast_clip(polygon, bbox); + if (fast_clip(polygon, bbox)) modified = true; } + size_t sizeBefore = mp.size(); mp.erase(std::remove_if( mp.begin(), mp.end(), [](const Polygon &poly) -> bool { return poly.outer().empty(); }), mp.end()); + return modified || mp.size() != sizeBefore; } diff --git a/src/mmap_allocator.cpp b/src/mmap_allocator.cpp index 42d2bf44..a68808a7 100644 --- a/src/mmap_allocator.cpp +++ b/src/mmap_allocator.cpp @@ -12,6 +12,12 @@ #include +#if defined(__unix__) || defined(__APPLE__) +#include +#else +#include +#endif + struct mmap_file { @@ -52,15 +58,30 @@ struct mmap_dir_t thread_local mmap_file_ptr mmap_file_thread_ptr; +// Anonymous memory that is committed lazily by the OS (pages are zero and only +// become resident when touched), unlike std::vector which zero-fills eagerly. +struct lazy_region +{ + uint8_t *data; + size_t size; + + lazy_region(size_t size); + ~lazy_region(); + lazy_region(lazy_region const &) = delete; + lazy_region &operator=(lazy_region const &) = delete; +}; + struct mmap_shm { std::mutex mutex; - std::vector region; + lazy_region region; boost::interprocess::managed_external_buffer buffer; static size_t mmap_file_size; mmap_shm(size_t size); + mmap_shm(mmap_shm const &) = delete; + mmap_shm &operator=(mmap_shm const &) = delete; static void open(size_t add_size); static void close(); @@ -158,9 +179,37 @@ bool mmap_dir_t::is_open() } } +lazy_region::lazy_region(size_t size) + : data(nullptr), size(size) +{ +#if defined(__unix__) || defined(__APPLE__) + int flags = MAP_PRIVATE | MAP_ANONYMOUS; +#ifdef MAP_NORESERVE + flags |= MAP_NORESERVE; +#endif + void *p = ::mmap(nullptr, size, PROT_READ | PROT_WRITE, flags, -1, 0); + if(p == MAP_FAILED) + throw std::bad_alloc(); +#else + void *p = std::calloc(size, 1); + if(p == nullptr) + throw std::bad_alloc(); +#endif + data = static_cast(p); +} + +lazy_region::~lazy_region() +{ +#if defined(__unix__) || defined(__APPLE__) + ::munmap(data, size); +#else + std::free(data); +#endif +} + mmap_shm::mmap_shm(size_t size) : region(size) - , buffer(boost::interprocess::create_only, region.data(), region.size()) + , buffer(boost::interprocess::create_only, region.data, region.size) { } void mmap_shm::open(size_t add_size) diff --git a/src/osm_lua_processing.cpp b/src/osm_lua_processing.cpp index b95b56d2..159a7bc7 100644 --- a/src/osm_lua_processing.cpp +++ b/src/osm_lua_processing.cpp @@ -1,5 +1,6 @@ #include #include +#include #include "osm_lua_processing.h" #include "attribute_store.h" @@ -160,26 +161,124 @@ kaguya::LuaTable rawAllKeys() { return getAllTags(*g_luaState, &tags); } -bool rawHolds(const KnownTagKey& key) { +bool rawHasTags() { return osmLuaProcessing->HasTags(); } +void rawSetTag(const std::string &key, const std::string &value) { return osmLuaProcessing->SetTag(key, value); } +// Native lua_CFunction versions of Holds/Find, bypassing kaguya's invoker. +// Semantics match the former kaguya-bound versions (KnownTagKey): kaguya does +// no type check for a single (non-overloaded) function, so the argument is +// read with lua_tolstring (numbers are coerced; other types give NULL/0 and +// simply miss). +static int luaHolds(lua_State* L) { + size_t size = 0; + const char* buffer = lua_tolstring(L, 1, &size); + if (osmLuaProcessing->isPostScanRelation) { - return osmLuaProcessing->Holds(key.stringValue); + bool threw = false, rv = false; + try { + rv = osmLuaProcessing->Holds(std::string(buffer ? buffer : "", size)); + } catch (std::exception& e) { + lua_pushstring(L, e.what()); + threw = true; + } + if (threw) return lua_error(L); + lua_pushboolean(L, rv); + return 1; } - return key.found; + lua_pushboolean(L, osmLuaProcessing->currentTags->getKey(buffer, size) >= 0); + return 1; +} + +// Upvalue 1: the empty string, returned on a miss. +static int luaFind(lua_State* L) { + size_t size = 0; + const char* buffer = lua_tolstring(L, 1, &size); + + if (osmLuaProcessing->isPostScanRelation) { + bool threw = false; + try { + const std::string rv = osmLuaProcessing->Find(std::string(buffer ? buffer : "", size)); + lua_pushlstring(L, rv.data(), rv.size()); + } catch (std::exception& e) { + lua_pushstring(L, e.what()); + threw = true; + } + if (threw) return lua_error(L); + return 1; + } + + const TagMap* tags = osmLuaProcessing->currentTags; + int64_t tagLoc = tags->getKey(buffer, size); + if (tagLoc >= 0) { + const protozero::data_view* value = tags->getValueFromKey(tagLoc); + lua_pushlstring(L, value->data(), value->size()); + } else { + lua_pushvalue(L, lua_upvalueindex(1)); + } + return 1; +} + +// Native lua_CFunction versions of Attribute*/MinZoom/ZOrder. Each wraps the +// kaguya binding registered under the same name (see wrapKaguyaBinding): +// upvalue 1 is kaguya's function tuple, upvalue 2 kaguya's closure. Calls whose +// arguments exactly match one overload (kaguya's strictCheckType, right arg +// count) are handled here; anything else goes to kaguya's own invoker on the +// same stack, so overload resolution, conversions and errors are unchanged. +static int kaguyaFallback(lua_State* L) { + return lua_tocfunction(L, lua_upvalueindex(2))(L); } -bool rawHasTags() { return osmLuaProcessing->HasTags(); } -void rawSetTag(const std::string &key, const std::string &value) { return osmLuaProcessing->SetTag(key, value); } -const std::string rawFind(const KnownTagKey& key) { - if (osmLuaProcessing->isPostScanRelation) - return osmLuaProcessing->Find(key.stringValue); - if (key.found) { - auto value = *(osmLuaProcessing->currentTags->getValueFromKey(key.index)); - return std::string(value.data(), value.size()); +// Exceptions are reported as kaguya's invoker does +template static int luaGuarded(lua_State* L, F&& fn) { + try { + fn(); + return 0; + } catch (std::exception& e) { + kaguya::util::traceBack(L, e.what()); + } catch (...) { + kaguya::util::traceBack(L, "Unknown exception"); } + return lua_error(L); +} + +template static int luaAttributeT(lua_State* L) { + const int n = lua_gettop(L); + if ((n != 2 && n != 3) || + !kaguya::lua_type_traits::strictCheckType(L, 1) || + !kaguya::lua_type_traits::strictCheckType(L, 2) || + (n == 3 && !kaguya::lua_type_traits::strictCheckType(L, 3))) + return kaguyaFallback(L); + + size_t keyLen = 0; + const char* key = lua_tolstring(L, 1, &keyLen); + const T val = kaguya::lua_type_traits::get(L, 2); + const char minzoom = n == 3 ? kaguya::lua_type_traits::get(L, 3) : 0; + return luaGuarded(L, [&]() { osmLuaProcessing->attributeFast(key, keyLen, val, minzoom, Type); }); +} - return EMPTY_STRING; +static int luaMinZoom(lua_State* L) { + if (lua_type(L, 1) != LUA_TNUMBER) return kaguyaFallback(L); + const double z = kaguya::lua_type_traits::get(L, 1); + return luaGuarded(L, [&]() { osmLuaProcessing->MinZoom(z); }); } + +static int luaZOrder(lua_State* L) { + if (lua_type(L, 1) != LUA_TNUMBER) return kaguyaFallback(L); + const double z = kaguya::lua_type_traits::get(L, 1); + return luaGuarded(L, [&]() { osmLuaProcessing->ZOrder(z); }); +} + +// Replace the kaguya-bound global `name` with native function fn, keeping +// kaguya's closure and function tuple as upvalues for kaguyaFallback. +static void wrapKaguyaBinding(lua_State* L, const char* name, lua_CFunction fn) { + lua_getglobal(L, name); + lua_getupvalue(L, -1, 1); + lua_pushvalue(L, -2); + lua_pushcclosure(L, fn, 2); + lua_setglobal(L, name); + lua_pop(L, 1); +} + std::vector rawFindIntersecting(const std::string &layerName) { return osmLuaProcessing->FindIntersecting(layerName); } bool rawIntersects(const std::string& layerName) { return osmLuaProcessing->Intersects(layerName); } std::vector rawFindCovering(const std::string& layerName) { return osmLuaProcessing->FindCovering(layerName); } @@ -223,6 +322,84 @@ int lua_error_handler(int errCode, const char *errMessage) // ---- initialization routines +// Same setup as kaguya::State's owning constructor (and luaL_newstate), but with LuaPool +// as the allocator: kaguya's allocator interface drops osize, which the pool needs. +#if !defined(LUAJIT) && LUA_VERSION_NUM <= 504 +#define TM_LUA_POOL 1 +#endif + +#if defined(TM_LUA_POOL) && LUA_VERSION_NUM == 504 +// Lua 5.4 luaL_newstate() also installs these (static in lauxlib.c): warnings default off +static void tmLuaWarnfOff(void *ud, const char *message, int tocont); +static void tmLuaWarnfOn(void *ud, const char *message, int tocont); +static void tmLuaWarnfCont(void *ud, const char *message, int tocont); +static int tmLuaWarnCheckControl(lua_State *L, const char *message, int tocont) { + if (tocont || *(message++) != '@') return 0; + if (strcmp(message, "off") == 0) lua_setwarnf(L, tmLuaWarnfOff, L); + else if (strcmp(message, "on") == 0) lua_setwarnf(L, tmLuaWarnfOn, L); + return 1; +} +static void tmLuaWarnfOff(void *ud, const char *message, int tocont) { + tmLuaWarnCheckControl((lua_State *)ud, message, tocont); +} +static void tmLuaWarnfCont(void *ud, const char *message, int tocont) { + lua_State *L = (lua_State *)ud; + lua_writestringerror("%s", message); + if (tocont) lua_setwarnf(L, tmLuaWarnfCont, L); + else { + lua_writestringerror("%s", "\n"); + lua_setwarnf(L, tmLuaWarnfOn, L); + } +} +static void tmLuaWarnfOn(void *ud, const char *message, int tocont) { + if (tmLuaWarnCheckControl((lua_State *)ud, message, tocont)) return; + lua_writestringerror("%s", "Lua warning: "); + tmLuaWarnfCont(ud, message, tocont); +} +#endif + +// Copies of kaguya::State's private initializing_panic/default_panic/stderror_out +static int tmLuaInitializingPanic(lua_State *L) { + kaguya::ErrorHandler::throwDefaultError(lua_status(L), lua_tostring(L, -1)); + return 0; +} +static int tmLuaDefaultPanic(lua_State *L) { + if (kaguya::ErrorHandler::handle(lua_status(L), L)) return 0; + fprintf(stderr, "PANIC: unprotected error in call to Lua API (%s)\n", lua_tostring(L, -1)); + fflush(stderr); + return 0; +} +static void tmLuaStderrorOut(int status, const char* message) { + std::cerr << message << std::endl; +} + +lua_State *OsmLuaProcessing::newLuaState(LuaPool &pool) { +#ifdef TM_LUA_POOL + lua_State *L = lua_newstate(&LuaPool::alloc, &pool); +#if LUA_VERSION_NUM == 504 + if (L) lua_setwarnf(L, tmLuaWarnfOff, L); +#endif + if (!L) L = luaL_newstate(); // e.g. LuaJIT not detected as such: behave like base +#else + lua_State *L = luaL_newstate(); +#endif + if (!L) return L; + // kaguya::State::init(AllLoadLibs()) + lua_atpanic(L, &tmLuaInitializingPanic); + try { + if (!kaguya::ErrorHandler::getHandler(L)) { + kaguya::util::ScopedSavedStack save(L); + kaguya::ErrorHandler::registerHandler(L, &tmLuaStderrorOut); + } + luaL_openlibs(L); + lua_atpanic(L, &tmLuaDefaultPanic); + } catch (const kaguya::LuaException&) { + lua_close(L); + L = nullptr; + } + return L; +} + OsmLuaProcessing::OsmLuaProcessing( OSMStore &osmStore, const class Config &configIn, @@ -235,6 +412,8 @@ OsmLuaProcessing::OsmLuaProcessing( bool materializeGeometries, bool isFirst) : osmStore(osmStore), + luaStateHolder(newLuaState(luaPool)), + luaState(luaStateHolder.L), shpMemTiles(shpMemTiles), osmMemTiles(osmMemTiles), attributeStore(attributeStore), @@ -259,8 +438,12 @@ OsmLuaProcessing::OsmLuaProcessing( luaState["OsmType"] = &rawOsmType; luaState["AllKeys"] = &rawAllKeys; luaState["AllTags"] = &rawAllTags; - luaState["Holds"] = &rawHolds; - luaState["Find"] = &rawFind; + lua_State* L = luaState.state(); + lua_pushcclosure(L, luaHolds, 0); + lua_setglobal(L, "Holds"); + lua_pushlstring(L, "", 0); + lua_pushcclosure(L, luaFind, 1); + lua_setglobal(L, "Find"); luaState["HasTags"] = &rawHasTags; luaState["SetTag"] = &rawSetTag; luaState["FindIntersecting"] = &rawFindIntersecting; @@ -295,6 +478,12 @@ OsmLuaProcessing::OsmLuaProcessing( luaState["MinZoom"] = &rawMinZoom; luaState["ZOrder"] = &rawZOrder; + wrapKaguyaBinding(L, "Attribute", luaAttributeT); + wrapKaguyaBinding(L, "AttributeNumeric", luaAttributeT); + wrapKaguyaBinding(L, "AttributeBoolean", luaAttributeT); + wrapKaguyaBinding(L, "AttributeInteger", luaAttributeT); + wrapKaguyaBinding(L, "MinZoom", luaMinZoom); + wrapKaguyaBinding(L, "ZOrder", luaZOrder); luaState["Score"] = &rawScore; luaState["Accept"] = &rawAccept; luaState["NextRelation"] = &rawNextRelation; @@ -594,10 +783,30 @@ double OsmLuaProcessing::Length() { } // Cached geometries creation +#ifdef TMOPT_VERIFY +template +static bool tmVerifySamePoints(const A &a, const B &b) { + if (a.size() != b.size()) return false; + for (size_t i = 0; i < a.size(); i++) + if (!(a[i].x() == b[i].x() && a[i].y() == b[i].y())) return false; + return true; +} +#endif + const Linestring &OsmLuaProcessing::linestringCached() { if (!linestringInited) { linestringInited = true; - linestringCache = osmStore.llListLinestring(llVecPtr->cbegin(),llVecPtr->cend()); + // Same as osmStore.llListLinestring(), refilling the cache in place + linestringCache.clear(); + linestringCache.reserve(llVecPtr->size()); + for (const LatpLon &ll : *llVecPtr) + linestringCache.push_back(geom::make(ll.lon/10000000.0, ll.latp/10000000.0)); +#ifdef TMOPT_VERIFY + if (!tmVerifySamePoints(linestringCache, osmStore.llListLinestring(llVecPtr->cbegin(),llVecPtr->cend()))) { + std::cerr << "TMOPT_VERIFY linestringCached mismatch on way " << originalOsmID << std::endl; + std::abort(); + } +#endif } return linestringCache; } @@ -667,11 +876,9 @@ void OsmLuaProcessing::Layer(const string &layerName, bool area) { else if (geomType==POLYGON_) { // polygon - MultiPolygon mp; - if (isRelation) { try { - mp = multiPolygonCached(); + MultiPolygon mp = multiPolygonCached(); if(CorrectGeometry(mp) == CorrectGeometryResult::Invalid) return; NodeID id = osmMemTiles.storeMultiPolygon(mp); OutputObject oo(geomType, layers.layerMap[layerName], id, 0, layerMinZoom); @@ -684,18 +891,35 @@ void OsmLuaProcessing::Layer(const string &layerName, bool area) { else if (isWay) { //Is there a more efficient way to do this? const Linestring &ls = linestringCached(); - Polygon p; - p.outer().reserve(ls.size()); - geom::assign_points(p, ls); - mp.reserve(1); - mp.push_back(std::move(p)); + // Reused scratch: a single polygon with ls as its outer ring and no inners + MultiPolygon &mp = layerScratchMultiPolygon; + mp.resize(1); + mp[0].outer().assign(ls.begin(), ls.end()); + mp[0].inners().clear(); +#ifdef TMOPT_VERIFY + { + Polygon p; + p.outer().reserve(ls.size()); + geom::assign_points(p, ls); + MultiPolygon old; + old.push_back(std::move(p)); + if (!(old.size() == mp.size() && tmVerifySamePoints(old[0].outer(), mp[0].outer()) && old[0].inners().empty() && mp[0].inners().empty())) { + std::cerr << "TMOPT_VERIFY Layer polygon scratch mismatch on way " << originalOsmID << std::endl; + std::abort(); + } + } +#endif - auto correctionResult = CorrectGeometry(mp); + auto correctionResult = CorrectGeometry(mp, /*allowReverse=*/isWay); if(correctionResult == CorrectGeometryResult::Invalid) return; NodeID id = 0; if (!materializeGeometries && correctionResult == CorrectGeometryResult::Valid) { id = USE_WAY_STORE | originalOsmID; wayEmitted = true; + } else if (!materializeGeometries && correctionResult == CorrectGeometryResult::Reversed) { + // Only defect was ring orientation: rebuild from the way store, reversed. + id = USE_WAY_STORE_REVERSED | originalOsmID; + wayEmitted = true; } else id = osmMemTiles.storeMultiPolygon(mp); OutputObject oo(geomType, layers.layerMap[layerName], id, 0, layerMinZoom); @@ -722,9 +946,26 @@ void OsmLuaProcessing::Layer(const string &layerName, bool area) { } else if (geomType==LINESTRING_) { // linestring - Linestring ls = linestringCached(); - - auto correctionResult = CorrectGeometry(ls); + // Check the cached linestring in place; copy it only if it needs correcting + const Linestring &cachedLs = linestringCached(); + geom::validity_failure_type failure = CheckGeometry(cachedLs); + Linestring correctedLs; + auto correctionResult = CorrectGeometryResult::Valid; + if (failure != geom::validity_failure_type::no_failure) { + correctedLs = cachedLs; + correctionResult = CorrectGeometryAfterCheck(correctedLs, failure); + } + const Linestring &ls = failure != geom::validity_failure_type::no_failure ? correctedLs : cachedLs; +#ifdef TMOPT_VERIFY + { + Linestring old = linestringCached(); + auto oldResult = CorrectGeometry(old); + if (oldResult != correctionResult || (oldResult != CorrectGeometryResult::Invalid && !tmVerifySamePoints(old, ls))) { + std::cerr << "TMOPT_VERIFY Layer linestring mismatch on way " << originalOsmID << std::endl; + std::abort(); + } + } +#endif if(correctionResult == CorrectGeometryResult::Invalid) return; if (isWay && !isRelation) { @@ -876,10 +1117,28 @@ void OsmLuaProcessing::noteIfDecluttered(LatpLon point) { } Point OsmLuaProcessing::calculateCentroid(CentroidAlgorithm algorithm) { + // Memoized per object: profiles often ask for the same centroid more than once + // (e.g. housenumber and POI on the same building). Exceptions are not memoized. + const int a = static_cast(algorithm); + if (centroidMemoValid[a]) { +#ifdef TMOPT_VERIFY + Point check = calculateCentroidUncached(algorithm); + if (memcmp(&check, ¢roidMemo[a], sizeof(Point)) != 0) { + std::cerr << "TMOPT_VERIFY centroid memo mismatch on " << originalOsmID << std::endl; + std::abort(); + } +#endif + return centroidMemo[a]; + } + centroidMemo[a] = calculateCentroidUncached(algorithm); + centroidMemoValid[a] = true; + return centroidMemo[a]; +} + +Point OsmLuaProcessing::calculateCentroidUncached(CentroidAlgorithm algorithm) { Point centroid; if (isRelation) { - MultiPolygon tmp; - tmp = multiPolygonCached(); + const MultiPolygon &tmp = multiPolygonCached(); if (algorithm == CentroidAlgorithm::Polylabel) { int index = 0; @@ -914,9 +1173,10 @@ Point OsmLuaProcessing::calculateCentroid(CentroidAlgorithm algorithm) { } } else { const Linestring &ls = linestringCached(); - Polygon p; - p.outer().reserve(ls.size()); - geom::assign_points(p, ls); + // Reused scratch polygon: ls as the outer ring, no inners + Polygon &p = centroidScratchPolygon; + p.outer().assign(ls.begin(), ls.end()); + p.inners().clear(); if (algorithm == CentroidAlgorithm::Polylabel) { // CONSIDER: pick precision intelligently @@ -965,17 +1225,17 @@ void OsmLuaProcessing::SetTag(const std::string &key, const std::string &value) osmStore.scannedRelations.set_relation_tag(originalOsmID, key, value); } -void OsmLuaProcessing::removeAttributeIfNeeded(const string& key) { +void OsmLuaProcessing::removeAttributeIfNeeded(uint16_t keyIndex) { // Does it exist? for (int i = 0; i < outputKeys.size(); i++) { - if (outputKeys[i] == key) { + if (outputKeys[i] == keyIndex) { AttributeSet& set = outputs.back().second; - set.removePairWithKey(attributeStore.pairStore, attributeStore.keyStore.key2index(key)); + set.removePairWithKey(attributeStore.pairStore, keyIndex); return; } } - outputKeys.push_back(key); + outputKeys.push_back(keyIndex); } // Force a new ID @@ -986,29 +1246,33 @@ void OsmLuaProcessing::ModifyId(const int newId) { // Set attributes in a vector tile's Attributes table void OsmLuaProcessing::Attribute(const string &key, const protozero::data_view val, const char minzoom) { if (outputs.size()==0) { ProcessingError("Can't add Attribute if no Layer set"); return; } - removeAttributeIfNeeded(key); - attributeStore.addAttribute(outputs.back().second, key, val, minzoom); + const uint16_t keyIndex = attributeStore.keyStore.key2index(key); + removeAttributeIfNeeded(keyIndex); + attributeStore.addAttributeWithKeyIndex(outputs.back().second, keyIndex, key, val, minzoom); setVectorLayerMetadata(outputs.back().first.layer, key, 0); } void OsmLuaProcessing::AttributeNumeric(const string &key, const double val, const char minzoom) { if (outputs.size()==0) { ProcessingError("Can't add Attribute if no Layer set"); return; } - removeAttributeIfNeeded(key); - attributeStore.addAttribute(outputs.back().second, key, val, minzoom); + const uint16_t keyIndex = attributeStore.keyStore.key2index(key); + removeAttributeIfNeeded(keyIndex); + attributeStore.addAttributeWithKeyIndex(outputs.back().second, keyIndex, val, minzoom); setVectorLayerMetadata(outputs.back().first.layer, key, 1); } void OsmLuaProcessing::AttributeBoolean(const string &key, const bool val, const char minzoom) { if (outputs.size()==0) { ProcessingError("Can't add Attribute if no Layer set"); return; } - removeAttributeIfNeeded(key); - attributeStore.addAttribute(outputs.back().second, key, val, minzoom); + const uint16_t keyIndex = attributeStore.keyStore.key2index(key); + removeAttributeIfNeeded(keyIndex); + attributeStore.addAttributeWithKeyIndex(outputs.back().second, keyIndex, val, minzoom); setVectorLayerMetadata(outputs.back().first.layer, key, 2); } void OsmLuaProcessing::AttributeInteger(const string &key, const int val, const char minzoom) { if (outputs.size()==0) { ProcessingError("Can't add Attribute if no Layer set"); return; } - removeAttributeIfNeeded(key); - attributeStore.addAttribute(outputs.back().second, key, val, minzoom); + const uint16_t keyIndex = attributeStore.keyStore.key2index(key); + removeAttributeIfNeeded(keyIndex); + attributeStore.addAttributeWithKeyIndex(outputs.back().second, keyIndex, val, minzoom); setVectorLayerMetadata(outputs.back().first.layer, key, 3); } @@ -1094,6 +1358,105 @@ void OsmLuaProcessing::setVectorLayerMetadata(const uint_least8_t layer, const s alreadySeen[layer].insert(key); } +// Per-thread cache of attribute keys, indexed by the address of the Lua +// string buffer. A hit needs the same address, length and bytes, so a reused +// address or a non-interned string simply misses. seenLayers caches, per +// layer, that setVectorLayerMetadata has already run for this key on this +// thread (alreadySeen is thread_local and never cleared, so a set bit stays +// true; an evicted entry just starts again with no bits set). +namespace { +struct AttrKeyCacheEntry { + const char* ptr = nullptr; + size_t len = 0; + std::string key; + uint16_t keyIndex = 0; + std::vector seenLayers; +}; +const size_t ATTR_KEY_CACHE_SIZE = 512; +thread_local std::vector attrKeyCache; +thread_local const AttributeKeyStore* attrKeyCacheStore = nullptr; + +AttrKeyCacheEntry& lookupAttrKey(AttributeKeyStore& keyStore, const char* key, size_t keyLen) { + if (attrKeyCacheStore != &keyStore) { + attrKeyCache.assign(ATTR_KEY_CACHE_SIZE, AttrKeyCacheEntry()); + attrKeyCacheStore = &keyStore; + } + const uintptr_t p = reinterpret_cast(key); + AttrKeyCacheEntry& e = attrKeyCache[((p >> 4) ^ (p >> 13)) & (ATTR_KEY_CACHE_SIZE - 1)]; + if (e.ptr == key && e.len == keyLen && e.key.size() == keyLen && memcmp(e.key.data(), key, keyLen) == 0) { +#ifdef TMOPT_VERIFY + if (keyStore.key2index(e.key) != e.keyIndex) { std::cerr << "TMOPT_VERIFY: keyIndex mismatch for " << e.key << std::endl; abort(); } +#endif + return e; + } + // resolve the index first: key2index may throw (too many keys), and a pcall in Lua could continue + std::string k(key, keyLen); + const uint16_t idx = keyStore.key2index(k); + e.key.swap(k); + e.keyIndex = idx; + e.ptr = key; + e.len = keyLen; + e.seenLayers.clear(); + return e; +} + +inline void addAttributeByIndex(AttributeStore& store, AttributeSet& set, const AttrKeyCacheEntry& e, const protozero::data_view val, char minzoom) { + store.addAttributeWithKeyIndex(set, e.keyIndex, e.key, val, minzoom); +} +template inline void addAttributeByIndex(AttributeStore& store, AttributeSet& set, const AttrKeyCacheEntry& e, T val, char minzoom) { + store.addAttributeWithKeyIndex(set, e.keyIndex, val, minzoom); +} +} + +// Same as Attribute/AttributeNumeric/AttributeBoolean/AttributeInteger, with +// the key index and vector_layers bookkeeping served from the key cache. +template +void OsmLuaProcessing::attributeFast(const char *key, size_t keyLen, T val, char minzoom, uint type) { + if (outputs.size()==0) { ProcessingError("Can't add Attribute if no Layer set"); return; } + AttrKeyCacheEntry& e = lookupAttrKey(attributeStore.keyStore, key, keyLen); + removeAttributeIfNeeded(e.keyIndex); + addAttributeByIndex(attributeStore, outputs.back().second, e, val, minzoom); + + const uint_least8_t layer = outputs.back().first.layer; + const size_t word = layer / 64; + const uint64_t bit = uint64_t(1) << (layer % 64); + if (e.seenLayers.size() > word && (e.seenLayers[word] & bit)) { +#ifdef TMOPT_VERIFY + if (alreadySeen.size() <= layer || !alreadySeen[layer].count(e.key)) { std::cerr << "TMOPT_VERIFY: seen-bit mismatch for " << e.key << std::endl; abort(); } +#endif + return; + } + setVectorLayerMetadata(layer, e.key, type); + if (e.seenLayers.size() <= word) e.seenLayers.resize(word + 1, 0); + e.seenLayers[word] |= bit; +} + +// Call global Lua function `name` with no arguments, as kaguya does for a +// function value (its pcall has no message handler; on error the registered +// error handler gets the status and message). Results are discarded. +// keepGlobals: luaState[name]() keeps the globals table on the stack below the +// call (visible in the error handler's stack dump); calling through a +// kaguya::LuaFunction doesn't. +bool OsmLuaProcessing::callLuaGlobal(const char *name, bool keepGlobals) { + lua_State* L = luaState.state(); + const int top = lua_gettop(L); +#if LUA_VERSION_NUM >= 502 + lua_pushglobaltable(L); +#else + lua_pushvalue(L, LUA_GLOBALSINDEX); +#endif + lua_getfield(L, -1, name); + if (!keepGlobals) lua_remove(L, -2); + if (lua_type(L, -1) != LUA_TFUNCTION) { + lua_settop(L, top); + return false; + } + const int status = lua_pcall(L, 0, 0, 0); + kaguya::except::checkErrorAndThrow(status, L); + lua_settop(L, top); + return true; +} + // Scan relation (but don't write geometry) // return true if we want it, false if we don't bool OsmLuaProcessing::scanRelation(WayID id, const TagMap& tags) { @@ -1143,7 +1506,7 @@ bool OsmLuaProcessing::setNode(NodeID id, LatpLon node, const TagMap& tags) { //Start Lua processing for node try { - luaState["node_function"](); + if (!callLuaGlobal("node_function", true)) luaState["node_function"](); } catch(luaProcessingException &e) { std::cerr << "Lua error on node " << originalOsmID << std::endl; exit(1); @@ -1172,6 +1535,7 @@ bool OsmLuaProcessing::setWay(WayID wayId, LatpLonVec const &llVec, const TagMap outerWayVecPtr = nullptr; innerWayVecPtr = nullptr; linestringInited = polygonInited = multiPolygonInited = false; + centroidMemoValid[0] = centroidMemoValid[1] = false; if (supportsReadingRelations && osmStore.scannedRelations.way_in_any_relations(wayId)) { relationList = osmStore.scannedRelations.relations_for_way(wayId); @@ -1192,9 +1556,11 @@ bool OsmLuaProcessing::setWay(WayID wayId, LatpLonVec const &llVec, const TagMap if (ok) { //Start Lua processing for way try { - kaguya::LuaFunction way_function = luaState["way_function"]; - kaguya::LuaRef ret = way_function(); - assert(!ret); + if (!callLuaGlobal("way_function", false)) { + kaguya::LuaFunction way_function = luaState["way_function"]; + kaguya::LuaRef ret = way_function(); + assert(!ret); + } } catch(luaProcessingException &e) { std::cerr << "Lua error on way " << originalOsmID << std::endl; exit(1); @@ -1202,7 +1568,7 @@ bool OsmLuaProcessing::setWay(WayID wayId, LatpLonVec const &llVec, const TagMap } if (!this->empty()) { - std::vector objects = finalizeOutputs(); + std::vector &objects = finalizeOutputs(); if (!objects.empty()) osmMemTiles.addGeometryToIndex(linestringCached(), objects, originalOsmID); return wayEmitted; } @@ -1240,10 +1606,11 @@ void OsmLuaProcessing::setRelation( // Start Lua processing for relation if (!isNativeMP && !supportsWritingRelations) return; try { - if (isNativeMP && supportsWritingWays) - luaState["way_function"](); - else if (!isNativeMP && supportsWritingRelations) - luaState["relation_function"](); + if (isNativeMP && supportsWritingWays) { + if (!callLuaGlobal("way_function", true)) luaState["way_function"](); + } else if (!isNativeMP && supportsWritingRelations) { + if (!callLuaGlobal("relation_function", true)) luaState["relation_function"](); + } } catch(luaProcessingException &e) { std::cerr << "Lua error on relation " << originalOsmID << std::endl; exit(1); @@ -1251,7 +1618,7 @@ void OsmLuaProcessing::setRelation( if (this->empty()) return; try { - std::vector objects = finalizeOutputs(); + std::vector &objects = finalizeOutputs(); if (objects.empty()) return; if (isClosed) { osmMemTiles.addGeometryToIndex(multiPolygonCached(), objects, originalOsmID); @@ -1282,8 +1649,9 @@ SignificantTags OsmLuaProcessing::GetSignificantWayKeys() { } -std::vector OsmLuaProcessing::finalizeOutputs() { - std::vector list; +std::vector &OsmLuaProcessing::finalizeOutputs() { + std::vector &list = finalizedOutputs; + list.clear(); list.reserve(this->outputs.size()); size_t nextDeclutter = 0; for (size_t i = 0; i < this->outputs.size(); i++) { diff --git a/src/osm_mem_tiles.cpp b/src/osm_mem_tiles.cpp index 0e9a6f07..37f465e6 100644 --- a/src/osm_mem_tiles.cpp +++ b/src/osm_mem_tiles.cpp @@ -49,7 +49,7 @@ Geometry OsmMemTiles::buildWayGeometry( const NodeID objectID, const TileBbox &bbox ) { - if (objectID < OSM_THRESHOLD || (geomType == POLYGON_ && IS_WAY(objectID))) { + if (objectID < OSM_THRESHOLD || (geomType == POLYGON_ && (IS_WAY(objectID) || IS_WAY_REV(objectID)))) { return TileDataSource::buildWayGeometry(geomType, objectID, bbox); } @@ -115,6 +115,14 @@ Linestring& OsmMemTiles::getOrBuildLinestring(NodeID objectID) const { // This is safe, because this function is the only thing that can // create/destroy entries in the cache, and the cache is thread-local. + // The cache is a pure memo of the way store; empty it whenever this thread + // starts a new tile batch so its linestrings don't pin heap between batches. + thread_local uint64_t lastBatch = 0; + if (lastBatch != clipCacheBatchId) { + lastBatch = clipCacheBatchId; + linestringCache.clear(); + } + Linestring* cachedEntry = linestringCache.get(objectID); if (cachedEntry != nullptr) return *cachedEntry; @@ -138,13 +146,35 @@ void OsmMemTiles::populateMultiPolygon(MultiPolygon& dst, NodeID objectID) { // // There are some exceptions, but some rough measurements showed caching // hurt overall throughput. - Linestring ls; - populateLinestring(ls, objectID); - Polygon p; - p.outer().reserve(ls.size()); - geom::assign_points(p, ls); + wayStore.at(OSM_ID(objectID), wayNodes); dst.reserve(dst.size() + 1); - dst.push_back(std::move(p)); + dst.emplace_back(); + Ring& outer = dst.back().outer(); + outer.reserve(wayNodes.size()); + for (const LatpLon& node : wayNodes) + outer.push_back(boost::geometry::make(node.lon/10000000.0, node.latp/10000000.0)); + if (IS_WAY_REV(objectID)) + std::reverse(outer.begin(), outer.end()); + +#ifdef TMOPT_VERIFY + { + Linestring ls; + populateLinestring(ls, objectID); + Polygon p; + p.outer().reserve(ls.size()); + geom::assign_points(p, ls); + if (IS_WAY_REV(objectID)) + std::reverse(p.outer().begin(), p.outer().end()); + const Ring& got = dst.back().outer(); + bool same = got.size() == p.outer().size() && dst.back().inners().empty(); + for (size_t i = 0; same && i < got.size(); i++) + same = got[i].x() == p.outer()[i].x() && got[i].y() == p.outer()[i].y(); + if (!same) { + std::cerr << "TMOPT_VERIFY: populateMultiPolygon mismatch for " << objectID << std::endl; + abort(); + } + } +#endif } void OsmMemTiles::Clear() { diff --git a/src/output_object.cpp b/src/output_object.cpp index e94168b0..628696ed 100644 --- a/src/output_object.cpp +++ b/src/output_object.cpp @@ -6,6 +6,7 @@ #include "helpers.h" #include "coordinates_geom.h" #include +#include #ifdef TILEMAKER_WITH_MLT #include "mlt_writer.h" #endif @@ -41,15 +42,14 @@ static void writeAttributesTo( FeatureBuilder& fbuilder, char zoom ) { - auto attr = attributeStore.getUnsafe(attributes); - - for(auto const &it: attr) { - if (it->minzoom > zoom) continue; + attributeStore.forEachPairUnsafe(attributes, [&](const AttributePair& pair) { + const AttributePair* it = &pair; + if (it->minzoom > zoom) return; // TODO: consider taking a data view that is stable // Look for key const std::string& key = attributeStore.keyStore.getKeyUnsafe(it->keyIndex); - + if (it->hasStringValue()) { fbuilder.add_property(key, it->stringValue()); } else if (it->hasBoolValue()) { @@ -60,22 +60,193 @@ static void writeAttributesTo( } else if (it->hasFloatValue()) { fbuilder.add_property(key, it->floatValue()); } + }); +} + +// MvtAttributeMemo + +static const uint32_t MEMO_MAX_ENTRIES = 1 << 16; // key/value indices, i.e. 32K properties +static const uint32_t MEMO_MAX_SETS = 1 << 13; +static const uint32_t MEMO_MAX_PAIRS = 1 << 14; +static const size_t MEMO_KEYS = 1 << 9; // AttributePair::keyIndex is 9 bits +static const uint32_t NO_INDEX = ~0u; + +static inline uint32_t memoHash(uint32_t x) { + // AttributeIndex/pair index is shard << 18 | offset: mix the high bits into the low ones + x ^= x >> 16; x *= 0x7feb352dU; + x ^= x >> 15; x *= 0x846ca68bU; + x ^= x >> 16; + return x; +} + +const MvtAttributeMemo::Table::Slot* MvtAttributeMemo::Table::find(uint32_t key) const { + if (slots.empty()) return nullptr; + const size_t mask = slots.size() - 1; + size_t i = memoHash(key) & mask; + while (slots[i].gen == gen) { + if (slots[i].key == key) return &slots[i]; + i = (i + 1) & mask; } + return nullptr; +} + +void MvtAttributeMemo::Table::insert(uint32_t key, uint32_t a, uint32_t b) { + if ((used + 1) * 2 > slots.size()) { + std::vector old; + old.swap(slots); + slots.assign(old.empty() ? 256 : old.size() * 2, Slot{0, 0, 0, 0}); + const size_t mask = slots.size() - 1; + for (const Slot& s : old) { + if (s.gen != gen) continue; + size_t i = memoHash(s.key) & mask; + while (slots[i].gen == gen) i = (i + 1) & mask; + slots[i] = s; + } + } + const size_t mask = slots.size() - 1; + size_t i = memoHash(key) & mask; + while (slots[i].gen == gen) i = (i + 1) & mask; + slots[i] = Slot{gen, key, a, b}; + used++; +} + +void MvtAttributeMemo::Table::clear() { + used = 0; + if (++gen == 0) { + for (Slot& s : slots) s.gen = 0; + gen = 1; + } +} + +void MvtAttributeMemo::reset(vtzero::layer_builder& layer, char zoom) { + this->layer = &layer; + this->zoom = zoom; + if (keyGen.empty()) { + keyGen.assign(MEMO_KEYS, 0); + keyIdx.assign(MEMO_KEYS, 0); + } + if (++layerGen == 0) { + std::fill(keyGen.begin(), keyGen.end(), 0); + layerGen = 1; + } + sets.clear(); + pairs.clear(); + entries.clear(); +} + +// Adds the pair's key and value to the layer tables exactly as writeAttributesTo would +void MvtAttributeMemo::resolvePair(const AttributeStore& attributeStore, uint32_t pairIndex, uint32_t& key, uint32_t& value) { + key = NO_INDEX; + value = NO_INDEX; + const AttributePair& it = attributeStore.pairStore.getPairUnsafe(pairIndex); + if (it.minzoom > zoom) return; + + if (it.hasStringValue()) { + const PooledString& ps = it.pooledString(); + value = layer->add_value(vtzero::encoded_property_value{vtzero::data_view{ps.data(), ps.size()}}).value(); + } else if (it.hasBoolValue()) { + value = layer->add_value(vtzero::encoded_property_value{it.boolValue()}).value(); + } else if (it.hasIntValue()) { + value = layer->add_value(vtzero::encoded_property_value{it.intValue()}).value(); + } else if (it.hasFloatValue()) { + value = layer->add_value(vtzero::encoded_property_value{it.floatValue()}).value(); + } else { + return; + } + + const uint16_t keyIndex = it.keyIndex; + if (keyGen[keyIndex] == layerGen) { + key = keyIdx[keyIndex]; + } else { + const std::string& keyString = attributeStore.keyStore.getKeyUnsafe(keyIndex); + key = layer->add_key(vtzero::data_view{keyString.data(), keyString.size()}).value(); + keyGen[keyIndex] = layerGen; + keyIdx[keyIndex] = key; + } +} + +const uint32_t* MvtAttributeMemo::resolve(const AttributeStore& attributeStore, AttributeIndex index, uint32_t& count) { + if (const Table::Slot* s = sets.find(index)) { + count = (s->b - s->a) / 2; + return entries.data() + s->a; + } + + // Miss: resolve each pair in order (the first time it is seen in this layer) + if (entries.size() >= MEMO_MAX_ENTRIES || sets.used >= MEMO_MAX_SETS) { + entries.clear(); + sets.clear(); + } + if (pairs.used >= MEMO_MAX_PAIRS) + pairs.clear(); + const uint32_t begin = entries.size(); + attributeStore.forEachPairIndexUnsafe(index, [&](uint32_t pairIndex) { + uint32_t key, value; + if (const Table::Slot* s = pairs.find(pairIndex)) { + key = s->a; + value = s->b; + } else { + resolvePair(attributeStore, pairIndex, key, value); + pairs.insert(pairIndex, key, value); + } + if (key == NO_INDEX) return; + entries.push_back(key); + entries.push_back(value); + }); + const uint32_t end = entries.size(); + sets.insert(index, begin, end); + + count = (end - begin) / 2; + return entries.data() + begin; } void OutputObject::writeAttributes( const AttributeStore& attributeStore, vtzero::feature_builder& fbuilder, - char zoom + char zoom, + MvtAttributeMemo* memo ) const { - writeAttributesTo(attributeStore, attributes, fbuilder, zoom); + if (!memo || !memo->layer || memo->zoom != zoom) { + writeAttributesTo(attributeStore, attributes, fbuilder, zoom); + return; + } + + uint32_t count; + const uint32_t* kv = memo->resolve(attributeStore, attributes, count); + for (uint32_t i = 0; i < count; i++) + fbuilder.add_property(vtzero::index_value_pair{kv[2*i], kv[2*i+1]}); + +#ifdef TMOPT_VERIFY + // Re-resolve through vtzero's find-or-append: it must return the same + // indices (and therefore add nothing to the layer tables). + uint32_t j = 0; + attributeStore.forEachPairUnsafe(attributes, [&](const AttributePair& it) { + if (it.minzoom > zoom) return; + vtzero::index_value value; + if (it.hasStringValue()) value = memo->layer->add_value(vtzero::encoded_property_value{it.stringValue()}); + else if (it.hasBoolValue()) value = memo->layer->add_value(vtzero::encoded_property_value{it.boolValue()}); + else if (it.hasIntValue()) value = memo->layer->add_value(vtzero::encoded_property_value{it.intValue()}); + else if (it.hasFloatValue()) value = memo->layer->add_value(vtzero::encoded_property_value{it.floatValue()}); + else return; + const vtzero::index_value key = memo->layer->add_key(vtzero::data_view{attributeStore.keyStore.getKeyUnsafe(it.keyIndex)}); + if (j >= count || key.value() != kv[2*j] || value.value() != kv[2*j+1]) { + std::cerr << "TMOPT_VERIFY writeAttributes memo mismatch" << std::endl; + abort(); + } + j++; + }); + if (j != count) { + std::cerr << "TMOPT_VERIFY writeAttributes memo count mismatch" << std::endl; + abort(); + } +#endif } #ifdef TILEMAKER_WITH_MLT void OutputObject::writeAttributes( const AttributeStore& attributeStore, MltFeatureBuilder& fbuilder, - char zoom + char zoom, + MvtAttributeMemo* /*memo*/ ) const { writeAttributesTo(attributeStore, attributes, fbuilder, zoom); } diff --git a/src/pbf_processor.cpp b/src/pbf_processor.cpp index 1ae42e62..eac83b41 100644 --- a/src/pbf_processor.cpp +++ b/src/pbf_processor.cpp @@ -97,7 +97,7 @@ bool PbfProcessor::ReadWays( LatpLonVec llVec; std::vector nodeVec; - for (PbfReader::Way pbfWay : pg.ways()) { + for (const PbfReader::Way& pbfWay : pg.ways()) { tags.reset(); readTags(pbfWay, pb, tags); @@ -124,8 +124,23 @@ bool PbfProcessor::ReadWays( nodeVec.reserve(pbfWay.refs.size()); bool skipToNext = false; + bool resolved = false; - for (int k=0; k 0 && (effectiveShards <= 1 || osmStore.nodes.contains(shard, pbfWay.refs[0]))) { + try { + llVec.resize(numRefs); + osmStore.nodes.at(pbfWay.refs.data(), numRefs, llVec.data()); + nodeVec.assign(pbfWay.refs.begin(), pbfWay.refs.end()); + resolved = true; + } catch (std::out_of_range &) { + llVec.clear(); + } + } + + for (int k=0; !resolved && k 1 && !osmStore.nodes.contains(shard, nodeId)) { @@ -208,7 +223,7 @@ bool PbfProcessor::ScanRelations(OsmLuaProcessing& output, PbfReader::PrimitiveG int typeKey = findStringPosition(pb, "type"); int mpKey = findStringPosition(pb, "multipolygon"); - for (PbfReader::Relation pbfRelation : pg.relations()) { + for (const PbfReader::Relation& pbfRelation : pg.relations()) { bool isMultiPolygon = relationIsType(pbfRelation, typeKey, mpKey); bool isAccepted = false; WayID relid = static_cast(pbfRelation.id); @@ -281,7 +296,7 @@ bool PbfProcessor::ReadRelations( int outerKey= findStringPosition(pb, "outer"); if (typeKey >-1) { int j = -1; - for (PbfReader::Relation pbfRelation : pg.relations()) { + for (const PbfReader::Relation& pbfRelation : pg.relations()) { j++; if (j % blockMetadata.chunks != blockMetadata.chunk) continue; @@ -344,6 +359,50 @@ bool PbfProcessor::ReadRelations( return true; } +std::shared_ptr PbfProcessor::acquireSharedBlock(std::istream& infile, const BlockMetadata& blockMetadata) { + std::unique_lock lock(sharedBlocksMutex); + auto it = sharedBlocks.find(blockMetadata.offset); + if (it == sharedBlocks.end()) { + // First chunk task to arrive: decompress the block outside the lock. + // std::map iterators stay valid while other entries are added/erased, + // and this entry can't be erased before we release it. + it = sharedBlocks.emplace(blockMetadata.offset, SharedBlock{nullptr, blockMetadata.chunks, false, false}).first; + lock.unlock(); + + std::shared_ptr data; + try { + data = std::make_shared(); + infile.seekg(blockMetadata.offset); + reader.readBlob(blockMetadata.length, infile, *data, true); + } catch (...) { + lock.lock(); + it->second.ready = true; + it->second.failed = true; + sharedBlocksCv.notify_all(); + throw; + } + + lock.lock(); + it->second.data = data; + it->second.ready = true; + sharedBlocksCv.notify_all(); + return it->second.data; + } + + SharedBlock& entry = it->second; + sharedBlocksCv.wait(lock, [&entry]() { return entry.ready; }); + if (entry.failed) + throw std::runtime_error("failed to read shared PBF block"); + return entry.data; +} + +void PbfProcessor::releaseSharedBlock(const BlockMetadata& blockMetadata) { + std::lock_guard lock(sharedBlocksMutex); + auto it = sharedBlocks.find(blockMetadata.offset); + if (it != sharedBlocks.end() && --it->second.remaining == 0) + sharedBlocks.erase(it); +} + // Returns true when block was completely handled, thus could be omited by another phases. bool PbfProcessor::ReadBlock( std::istream& infile, @@ -357,11 +416,36 @@ bool PbfProcessor::ReadBlock( uint effectiveShards ) { - infile.seekg(blockMetadata.offset); - - protozero::data_view blob = reader.readBlob(blockMetadata.length, infile); + // Subdivided block: decompress once and share the buffer among its chunk tasks. + std::shared_ptr sharedBlob; + struct SharedBlockReleaser { + PbfProcessor* processor; + const BlockMetadata* blockMetadata; + ~SharedBlockReleaser() { if (processor) processor->releaseSharedBlock(*blockMetadata); } + } sharedBlockReleaser{nullptr, &blockMetadata}; + + protozero::data_view blob; + if (blockMetadata.chunks > 1) { + sharedBlob = acquireSharedBlock(infile, blockMetadata); + sharedBlockReleaser.processor = this; + blob = { sharedBlob->data(), sharedBlob->size() }; +#ifdef TMOPT_VERIFY + { + infile.clear(); + infile.seekg(blockMetadata.offset); + protozero::data_view fresh = reader.readBlob(blockMetadata.length, infile); + if (fresh.size() != blob.size() || memcmp(fresh.data(), blob.data(), blob.size()) != 0) { + std::cerr << "TMOPT_VERIFY: shared block mismatch at offset " << blockMetadata.offset << std::endl; + abort(); + } + } +#endif + } else { + infile.seekg(blockMetadata.offset); + blob = reader.readBlob(blockMetadata.length, infile); + } PbfReader::PrimitiveBlock& pb = reader.readPrimitiveBlock(blob); - if (infile.eof()) { + if (!sharedBlob && infile.eof()) { return true; } @@ -725,6 +809,11 @@ int PbfProcessor::ReadPbfFile( } pool.join(); + { + // drop any shared-block entries left over (e.g. re-subdivided blocks) so none outlive this pass/file + std::lock_guard l(sharedBlocksMutex); + sharedBlocks.clear(); + } #ifdef CLOCK_MONOTONIC clock_gettime(CLOCK_MONOTONIC, &end); @@ -734,6 +823,7 @@ int PbfProcessor::ReadPbfFile( } if(phase == ReadPhase::RelationScan) { + osmStore.finalizeUsedWays(); auto output = generate_output(); output->postScanRelations(); } diff --git a/src/pbf_reader.cpp b/src/pbf_reader.cpp index d12cfb28..75acca09 100644 --- a/src/pbf_reader.cpp +++ b/src/pbf_reader.cpp @@ -63,6 +63,12 @@ PbfReader::BlobHeader PbfReader::PbfReader::readBlobHeader(std::istream& input) } protozero::data_view PbfReader::PbfReader::readBlob(int32_t datasize, std::istream& input) { + return readBlob(datasize, input, blobStorage2, false); +} + +// Like readBlob, but the decompressed data is stored in `output`, which the caller +// owns; with copyRaw, uncompressed blobs are copied there too. +protozero::data_view PbfReader::PbfReader::readBlob(int32_t datasize, std::istream& input, std::string& output, bool copyRaw) { blobStorage.resize(datasize); input.read(&blobStorage[0], datasize); if (input.eof()) @@ -99,15 +105,19 @@ protozero::data_view PbfReader::PbfReader::readBlob(int32_t datasize, std::istre if (dataType == BlobDataType::None) throw std::runtime_error("Blob: missing data"); - if (dataType == BlobDataType::Raw) - return view; + if (dataType == BlobDataType::Raw) { + if (!copyRaw) + return view; + output.assign(view.data(), view.size()); + return { output.data(), output.size() }; + } if (rawSize != -1) - blobStorage2.resize(rawSize); + output.resize(rawSize); else - blobStorage2.clear(); - decompress_string(blobStorage2, view.data(), view.size(), false); - return { blobStorage2.data(), blobStorage2.size() }; + output.clear(); + decompress_string(output, view.data(), view.size(), false); + return { output.data(), output.size() }; } PbfReader::HeaderBBox PbfReader::PbfReader::readHeaderBBox(protozero::data_view data) { diff --git a/src/shp_processor.cpp b/src/shp_processor.cpp index 4cae1179..f74e0704 100644 --- a/src/shp_processor.cpp +++ b/src/shp_processor.cpp @@ -2,6 +2,9 @@ #include #include +#include +#include +#include #include #include #include @@ -26,6 +29,72 @@ namespace { using ShpHandlePtr = std::unique_ptr::type, ShpHandleCloser>; using DbfHandlePtr = std::unique_ptr::type, DbfHandleCloser>; + + struct FileCloser { + void operator()(FILE *f) const { + fclose(f); + } + }; + using FilePtr = std::unique_ptr; + + // Reads record bounding boxes straight from the .shp record header, so that + // records outside the clipping box can be skipped without SHPReadObject. + // Uses its own read-only file handle; shapelib's handle is not touched. + class ShpBoundsReader { + public: + ShpBoundsReader(const string &filename, SHPHandle shp) : shp(shp) { + const uint16_t one = 1; + uint8_t firstByte; + memcpy(&firstByte, &one, 1); + if (firstByte != 1) return; // big-endian host: no prefilter + // Same .shp path resolution as SHPOpen: strip extension, try .shp then .SHP + string base = filename; + size_t slash = base.find_last_of("/\\"); + size_t dot = base.find_last_of('.'); + if (dot != string::npos && dot > 0 && (slash == string::npos || dot > slash)) base = base.substr(0, dot); + file.reset(fopen((base + ".shp").c_str(), "rb")); + if (!file) file.reset(fopen((base + ".SHP").c_str(), "rb")); + if (file) setvbuf(file.get(), nullptr, _IONBF, 0); + } + + // Returns true and fills the bounds if record i has a header bbox that matches + // what SHPReadObject would report as dfXMin/dfYMin/dfXMax/dfYMax. + bool bounds(int i, double &xmin, double &ymin, double &xmax, double &ymax) { + if (!file || i < 0 || i >= shp->nRecords) return false; + uint64_t offset = shp->panRecOffset[i]; + uint64_t size = shp->panRecSize[i]; // content bytes, excluding the 8-byte record header + if (offset < 100 || size < 36) return false; +#ifdef _MSC_VER + if (_fseeki64(file.get(), (__int64)offset, SEEK_SET) != 0) return false; +#else + if (fseeko(file.get(), (off_t)offset, SEEK_SET) != 0) return false; +#endif + uint8_t buf[44]; + if (fread(buf, 1, 44, file.get()) != 44) return false; + // Record header: big-endian content length in 16-bit words + uint64_t contentLength = ((uint64_t)buf[4] << 24) | ((uint64_t)buf[5] << 16) | ((uint64_t)buf[6] << 8) | buf[7]; + if (contentLength * 2 != size) return false; + int32_t shapeType = (int32_t)((uint32_t)buf[8] | ((uint32_t)buf[9] << 8) | ((uint32_t)buf[10] << 16) | ((uint32_t)buf[11] << 24)); + switch (shapeType) { + case SHPT_ARC: case SHPT_POLYGON: case SHPT_MULTIPOINT: + case SHPT_ARCZ: case SHPT_POLYGONZ: case SHPT_MULTIPOINTZ: + case SHPT_ARCM: case SHPT_POLYGONM: case SHPT_MULTIPOINTM: + case SHPT_MULTIPATCH: + break; + default: + return false; + } + memcpy(&xmin, buf + 12, 8); + memcpy(&ymin, buf + 20, 8); + memcpy(&xmax, buf + 28, 8); + memcpy(&ymax, buf + 36, 8); + return true; + } + + private: + SHPHandle shp; + FilePtr file; + }; } /* @@ -165,18 +234,31 @@ void ShpProcessor::read(class LayerDef &layer, uint layerNum) std::exception_ptr error; std::mutex errorMutex; + // Check shape bounds are outside clippingBox + auto outsideClippingBox = [&](double xmin, double ymin, double xmax, double ymax) { + Box shapeBox(Point(xmin, lat2latp(ymin)), Point(xmax, lat2latp(ymax))); + return shapeBox.min_corner().get<0>() > clippingBox.max_corner().get<0>() || + shapeBox.max_corner().get<0>() < clippingBox.min_corner().get<0>() || + shapeBox.min_corner().get<1>() > clippingBox.max_corner().get<1>() || + shapeBox.max_corner().get<1>() < clippingBox.min_corner().get<1>(); + }; + ShpBoundsReader boundsReader(filename, shp.get()); + boost::asio::thread_pool pool(threadNum); for (int i=0; i shape(rawShape, SHPDestroyObject); // Check shape is in clippingBox - Box shapeBox(Point(shape->dfXMin, lat2latp(shape->dfYMin)), Point(shape->dfXMax, lat2latp(shape->dfYMax))); - if (shapeBox.min_corner().get<0>() > clippingBox.max_corner().get<0>() || - shapeBox.max_corner().get<0>() < clippingBox.min_corner().get<0>() || - shapeBox.min_corner().get<1>() > clippingBox.max_corner().get<1>() || - shapeBox.max_corner().get<1>() < clippingBox.min_corner().get<1>()) { + if (outsideClippingBox(shape->dfXMin, shape->dfYMin, shape->dfXMax, shape->dfYMax)) { continue; } @@ -284,7 +366,11 @@ void ShpProcessor::processShapeGeometry(SHPObject* shape, AttributeIndex attrIdx // Make valid if needs be string reason; - if (!geom::is_valid(multi, reason)) { + const bool precheckValid = multiPolygonDefinitelyValid(multi); +#ifdef TMOPT_VERIFY + precheckVerify(multi, precheckValid, PRECHECK_SHP); +#endif + if (!precheckValid && !geom::is_valid(multi, reason)) { if (verbose) cerr << "Shapefile entity " << shape->nShapeId << " type " << shapeType << " is invalid. Parts:" << shape->nParts << ". Reason:" << reason; make_valid(multi); if (verbose) { diff --git a/src/sorted_node_store.cpp b/src/sorted_node_store.cpp index 1b9c452e..1930697e 100644 --- a/src/sorted_node_store.cpp +++ b/src/sorted_node_store.cpp @@ -9,32 +9,311 @@ #include "external/libpopcnt.h" #include "external/streamvbyte.h" #include "external/streamvbyte_zigzag.h" +#ifdef _MSC_VER +#include +#endif namespace SortedNodeStoreTypes { const uint16_t GroupSize = 256; const uint16_t ChunkSize = 256; const uint16_t ChunkAlignment = 16; const uint32_t ChunkCompressed = 1 << 31; - const uint8_t ChunkCacheSize = 4; + + // Direct-mapped, per-thread cache of decoded compressed chunks. + // Decoded values are stored compactly in rank order (the order of the set + // bits of the chunk's node mask), together with the mask and the rank at + // the start of each mask word, so a hit is a bit test, a popcount and a load. + const size_t ChunkCacheBits = 8; + const size_t ChunkCacheSize = 1 << ChunkCacheBits; + const uint64_t ChunkCacheEmpty = ~0ull; struct CachedChunk { - CachedChunk(): id(-1) {} + uint64_t nodeMask[4]; + uint32_t rankBase[4]; + LatpLon nodes[ChunkSize]; + }; + + struct ChunkCache { + ChunkCache() { + for (size_t i = 0; i < ChunkCacheSize; i++) + keys[i] = ChunkCacheEmpty; + } - int64_t id; - std::array lons; - std::array latps; + // Keys are kept apart from the (large) entries so that untouched + // entries are never paged in. + uint64_t keys[ChunkCacheSize]; + CachedChunk entries[ChunkCacheSize]; }; + thread_local ChunkCache* chunkCache = nullptr; + thread_local std::unique_ptr chunkCacheOwner; + + ChunkCache& getChunkCache() { + if (chunkCache == nullptr) { + chunkCacheOwner.reset(new ChunkCache()); + chunkCache = chunkCacheOwner.get(); + } + return *chunkCache; + } + + inline size_t chunkCacheIndex(uint64_t key) { + return (size_t)((key * 0x9E3779B97F4A7C15ull) >> (64 - ChunkCacheBits)); + } + + std::atomic nextGeneration(1); + +#if (defined(__GNUC__) || defined(__clang__)) && defined(__x86_64__) && !defined(__POPCNT__) && \ + defined(HAVE_ASM_POPCNT) && defined(HAVE_CPUID) + // Without -mpopcnt, __builtin_popcountll is a call into libgcc. Use the + // POPCNT instruction (libpopcnt's popcnt64) when the CPU has it. + const bool cpuHasPopcnt = (get_cpuid() & bit_POPCNT) != 0; +#define SORTED_NODE_STORE_RUNTIME_POPCNT +#endif + + inline uint64_t popcount64(uint64_t x) { +#if defined(SORTED_NODE_STORE_RUNTIME_POPCNT) + if (cpuHasPopcnt) + return popcnt64(x); + return ::popcount64(x); +#elif defined(__GNUC__) || defined(__clang__) + return __builtin_popcountll(x); +#elif defined(_MSC_VER) && defined(_M_X64) + return __popcnt64(x); +#else + x = x - ((x >> 1) & 0x5555555555555555ull); + x = (x & 0x3333333333333333ull) + ((x >> 2) & 0x3333333333333333ull); + x = (x + (x >> 4)) & 0x0F0F0F0F0F0F0F0Full; + return (x * 0x0101010101010101ull) >> 56; +#endif + } + + inline unsigned countTrailingZeros64(uint64_t x) { +#if defined(__GNUC__) || defined(__clang__) + return __builtin_ctzll(x); +#elif defined(_MSC_VER) && defined(_M_X64) + unsigned long index; + _BitScanForward64(&index, x); + return index; +#else + return popcount64((x & (0 - x)) - 1); +#endif + } + + // Load 8 mask bytes as a word whose bit j is bit (j % 8) of byte (j / 8), + // matching how the masks are addressed bytewise. Compiles to a plain load + // on little-endian targets. + inline uint64_t loadMaskWord(const uint8_t* p) { + return (uint64_t)p[0] | ((uint64_t)p[1] << 8) | ((uint64_t)p[2] << 16) | ((uint64_t)p[3] << 24) | + ((uint64_t)p[4] << 32) | ((uint64_t)p[5] << 40) | ((uint64_t)p[6] << 48) | ((uint64_t)p[7] << 56); + } + + // Number of set bits in the 256-bit mask below position `bit`. + inline size_t maskRank(const uint8_t* mask, size_t bit) { + const size_t word = bit / 64; + size_t rv = 0; + for (size_t i = 0; i < word; i++) + rv += popcount64(loadMaskWord(mask + i * 8)); + rv += popcount64(loadMaskWord(mask + word * 8) & ((1ull << (bit % 64)) - 1)); + return rv; + } + + inline size_t maskCount(const uint8_t* mask) { + return popcount64(loadMaskWord(mask)) + popcount64(loadMaskWord(mask + 8)) + + popcount64(loadMaskWord(mask + 16)) + popcount64(loadMaskWord(mask + 24)); + } + + inline bool maskTest(const uint8_t* mask, size_t bit) { + return mask[bit / 8] & (1 << (bit % 8)); + } + + // Locate the chunk holding `id`; throws std::out_of_range if its group or + // chunk is missing. + inline const ChunkInfoBase* findChunk(GroupInfo* const* groups, NodeID id) { + const size_t groupIndex = id / (GroupSize * ChunkSize); + const size_t chunk = (id % (GroupSize * ChunkSize)) / ChunkSize; + + const GroupInfo* groupPtr = groups[groupIndex]; + + if (groupPtr == nullptr) { + throw std::out_of_range("SortedNodeStore::at(" + std::to_string(id) + ") uses non-existent group " + std::to_string(groupIndex)); + } + + if (!maskTest(groupPtr->chunkMask, chunk)) + throw std::out_of_range("SortedNodeStore: node " + std::to_string(id) + " missing, no chunk"); + const size_t chunkOffset = maskRank(groupPtr->chunkMask, chunk); + + uint16_t scaledOffset = groupPtr->chunkOffsets[chunkOffset]; + return (const ChunkInfoBase*)(((const char *)(groupPtr->chunkOffsets + maskCount(groupPtr->chunkMask))) + (scaledOffset * ChunkAlignment)); + } + + inline LatpLon uncompressedNode(const ChunkInfoBase* basePtr, NodeID id) { + const UncompressedChunkInfo* ptr = (const UncompressedChunkInfo*)basePtr; + const size_t nodeIndex = id % ChunkSize; + if (!maskTest(ptr->nodeMask, nodeIndex)) + throw std::out_of_range("SortedNodeStore: node " + std::to_string(id) + " missing, no node"); + const size_t nodeOffset = maskRank(ptr->nodeMask, nodeIndex); + + return ptr->nodes[nodeOffset]; + } + + inline LatpLon cachedNode(const CachedChunk& cached, NodeID id) { + const size_t nodeIndex = id % ChunkSize; + const size_t w = nodeIndex / 64; + const uint64_t bit = 1ull << (nodeIndex % 64); + const uint64_t word = cached.nodeMask[w]; + if (!(word & bit)) + throw std::out_of_range("SortedNodeStore: node " + std::to_string(id) + " missing, no node"); + return cached.nodes[cached.rankBase[w] + popcount64(word & (bit - 1))]; + } + + // Same as svb_zigzag_decode_32 in streamvbyte_zigzag.c. + inline int32_t zigzagDecode(uint32_t val) { + return (val >> 1) ^ (0-(val & 1)); + } + + // prev + delta with wrap-around, as zigzag_delta_decode computes it. + inline int32_t addDelta(int32_t prev, uint32_t zigzagged) { + return (int32_t)((uint32_t)prev + (uint32_t)zigzagDecode(zigzagged)); + } + + const uint32_t svbLengthMask[4] = { 0xFFu, 0xFFFFu, 0xFFFFFFu, 0xFFFFFFFFu }; + + // Read a streamvbyte value of (code + 1) bytes at p; p + 4 must not be + // beyond the end of the stream. Assumes little endian, as streamvbyte does. + inline uint32_t svbLoadUnchecked(const uint8_t* p, uint32_t code) { + uint32_t val; + memcpy(&val, p, 4); + return val & svbLengthMask[code]; + } + + // Read a streamvbyte value of (code + 1) bytes at p, never reading at or + // beyond `end`. + inline uint32_t svbLoad(const uint8_t* p, const uint8_t* end, uint32_t code) { + if (end - p >= 4) + return svbLoadUnchecked(p, code); + uint32_t val = 0; + for (uint32_t b = 0; b <= code; b++) + val |= (uint32_t)p[b] << (8 * b); + return val; + } + + // Decode a compressed chunk into `out`, in rank order. + void decodeChunk(const CompressedChunkInfo* ptr, CachedChunk& out) { + const size_t latpSize = (ptr->flags >> 10) & ((1 << 10) - 1); + const size_t lonSize = ptr->flags & ((1 << 10) - 1); + + uint32_t count = 0; + for (size_t w = 0; w < 4; w++) { + const uint64_t word = loadMaskWord(ptr->nodeMask + w * 8); + out.nodeMask[w] = word; + out.rankBase[w] = count; + count += popcount64(word); + } + const size_t n = count - 1; + + const uint8_t* latpData = ptr->data; + const uint8_t* lonData = ptr->data + latpSize; + int32_t latp = ptr->firstLatp; + int32_t lon = ptr->firstLon; + out.nodes[0] = { latp, lon }; + LatpLon* dst = out.nodes + 1; + + if (n < 64) { + // Fused scalar streamvbyte + zigzag delta decode of both streams. + const size_t keyLen = (n + 3) / 4; + const uint8_t* latpEnd = latpData + latpSize; + const uint8_t* lonEnd = lonData + lonSize; + const uint8_t* latpPtr = latpData + keyLen; + const uint8_t* lonPtr = lonData + keyLen; + // Each value takes at least one byte, so a value followed by at + // least 3 more lies at least 4 bytes before the end of its stream. + const size_t unchecked = n > 3 ? n - 3 : 0; + size_t k = 0; + for (; k < unchecked; k++) { + const unsigned shift = (k & 3) * 2; + const uint32_t latpCode = (latpData[k >> 2] >> shift) & 3; + const uint32_t lonCode = (lonData[k >> 2] >> shift) & 3; + latp = addDelta(latp, svbLoadUnchecked(latpPtr, latpCode)); + lon = addDelta(lon, svbLoadUnchecked(lonPtr, lonCode)); + latpPtr += latpCode + 1; + lonPtr += lonCode + 1; + dst[k] = { latp, lon }; + } + for (; k < n; k++) { + const unsigned shift = (k & 3) * 2; + const uint32_t latpCode = (latpData[k >> 2] >> shift) & 3; + const uint32_t lonCode = (lonData[k >> 2] >> shift) & 3; + latp = addDelta(latp, svbLoad(latpPtr, latpEnd, latpCode)); + lon = addDelta(lon, svbLoad(lonPtr, lonEnd, lonCode)); + latpPtr += latpCode + 1; + lonPtr += lonCode + 1; + dst[k] = { latp, lon }; + } + return; + } + + uint32_t latpDeltas[ChunkSize]; + uint32_t lonDeltas[ChunkSize]; + streamvbyte_decode(latpData, latpDeltas, n); + streamvbyte_decode(lonData, lonDeltas, n); + for (size_t k = 0; k < n; k++) { + latp = addDelta(latp, latpDeltas[k]); + lon = addDelta(lon, lonDeltas[k]); + dst[k] = { latp, lon }; + } + } + +#ifdef TMOPT_VERIFY + // Decode the chunk the original way (streamvbyte + zigzag_delta_decode, + // scattered by position) and compare with the compact entry. + void verifyDecodedChunk(const CompressedChunkInfo* ptr, const CachedChunk& got) { + size_t latpSize = (ptr->flags >> 10) & ((1 << 10) - 1); + const size_t count = maskCount(ptr->nodeMask); + const size_t n = count - 1; + uint32_t recovdata[ChunkSize]; + int32_t latps[ChunkSize]; + int32_t lons[ChunkSize]; + streamvbyte_decode(ptr->data, recovdata, n); + latps[0] = ptr->firstLatp; + zigzag_delta_decode(recovdata, &latps[1], n, latps[0]); + streamvbyte_decode(ptr->data + latpSize, recovdata, n); + lons[0] = ptr->firstLon; + zigzag_delta_decode(recovdata, &lons[1], n, lons[0]); + + size_t k = 0; + for (size_t pos = 0; pos < ChunkSize; pos++) { + const bool present = maskTest(ptr->nodeMask, pos); + const bool gotPresent = got.nodeMask[pos / 64] & (1ull << (pos % 64)); + if (present != gotPresent) { + std::cerr << "TMOPT_VERIFY: SortedNodeStore mask mismatch at " << pos << std::endl; + abort(); + } + if (!present) + continue; + const LatpLon v = cachedNode(got, pos); + if (k >= count || v.latp != latps[k] || v.lon != lons[k]) { + std::cerr << "TMOPT_VERIFY: SortedNodeStore decode mismatch at " << pos << std::endl; + abort(); + } + k++; + } + if (k != count) { + std::cerr << "TMOPT_VERIFY: SortedNodeStore count mismatch" << std::endl; + abort(); + } + } +#define TMOPT_VERIFY_CHUNK(ptr, entry) verifyDecodedChunk(ptr, entry) +#else +#define TMOPT_VERIFY_CHUNK(ptr, entry) do {} while (0) +#endif + struct ThreadStorage { ThreadStorage(): collectingOrphans(true), groupStart(-1), localNodes(nullptr), arenaSpace(0), - arenaPtr(nullptr) { - for (uint8_t i = 0; i < ChunkCacheSize; ++i) - cacheOrder[i] = i; - } + arenaPtr(nullptr) {} // When SortedNodeStore first starts, it's not confident that it has seen an // entire segment, so it's in "collecting orphans" mode. Once it crosses a // threshold of 64K elements, it ceases to be in this mode. @@ -45,9 +324,6 @@ namespace SortedNodeStoreTypes { uint64_t groupStart = -1; std::vector* localNodes = nullptr; - std::array cacheChunks; - std::array cacheOrder; - uint32_t arenaSpace = 0; char* arenaPtr = nullptr; }; @@ -64,13 +340,6 @@ namespace SortedNodeStoreTypes { auto& rv = threadStorage.back(); return rv.second; } - - void promoteCacheSlot(ThreadStorage& tls, size_t orderIndex) { - const uint8_t slot = tls.cacheOrder[orderIndex]; - for (size_t i = orderIndex; i > 0; --i) - tls.cacheOrder[i] = tls.cacheOrder[i - 1]; - tls.cacheOrder[0] = slot; - } } using namespace SortedNodeStoreTypes; @@ -83,6 +352,7 @@ SortedNodeStore::SortedNodeStore(bool compressNodes): compressNodes(compressNode void SortedNodeStore::reopen() { allocatedMemory.clear(); + generation = nextGeneration++; totalNodes = 0; totalGroups = 0; @@ -146,102 +416,83 @@ bool SortedNodeStore::contains(size_t shard, NodeID id) const { } LatpLon SortedNodeStore::at(const NodeID id) const { - const size_t groupIndex = id / (GroupSize * ChunkSize); - const size_t chunk = (id % (GroupSize * ChunkSize)) / ChunkSize; - const uint64_t chunkMaskByte = chunk / 8; - const uint64_t chunkMaskBit = chunk % 8; - - const uint64_t nodeMaskByte = (id % ChunkSize) / 8; - const uint64_t nodeMaskBit = id % 8; - - GroupInfo* groupPtr = groups[groupIndex]; - - if (groupPtr == nullptr) { - throw std::out_of_range("SortedNodeStore::at(" + std::to_string(id) + ") uses non-existent group " + std::to_string(groupIndex)); + const uint64_t chunkId = id / ChunkSize; + if ((chunkId >> 32) != 0) { + // Not cacheable (chunk ids >= 2^32): decode into scratch space. + const ChunkInfoBase* basePtr = findChunk(groups.data(), id); + if (!(basePtr->flags & ChunkCompressed)) + return uncompressedNode(basePtr, id); + CachedChunk scratch; + decodeChunk((const CompressedChunkInfo*)basePtr, scratch); + TMOPT_VERIFY_CHUNK((const CompressedChunkInfo*)basePtr, scratch); + return cachedNode(scratch, id); } - size_t chunkOffset = 0; - { - chunkOffset = popcnt(groupPtr->chunkMask, chunkMaskByte); - uint8_t maskByte = groupPtr->chunkMask[chunkMaskByte]; - maskByte = maskByte & ((1 << chunkMaskBit) - 1); - chunkOffset += popcnt(&maskByte, 1); + // Fast path: the chunk was decoded recently on this thread. + const uint64_t key = (generation << 32) | chunkId; + ChunkCache& cache = getChunkCache(); + const size_t cacheIndex = chunkCacheIndex(key); + if (cache.keys[cacheIndex] == key) + return cachedNode(cache.entries[cacheIndex], id); + + const ChunkInfoBase* basePtr = findChunk(groups.data(), id); + if (!(basePtr->flags & ChunkCompressed)) + return uncompressedNode(basePtr, id); + + CachedChunk& entry = cache.entries[cacheIndex]; + decodeChunk((const CompressedChunkInfo*)basePtr, entry); + TMOPT_VERIFY_CHUNK((const CompressedChunkInfo*)basePtr, entry); + cache.keys[cacheIndex] = key; + return cachedNode(entry, id); +} - if (!(groupPtr->chunkMask[chunkMaskByte] & (1 << chunkMaskBit))) - throw std::out_of_range("SortedNodeStore: node " + std::to_string(id) + " missing, no chunk"); - } +void SortedNodeStore::at(const NodeID* ids, size_t n, LatpLon* out) const { + ChunkCache& cache = getChunkCache(); + // Consecutive ids usually share a chunk: reuse its cache entry directly. + uint64_t lastChunkId = ~0ull; + const CachedChunk* entry = nullptr; + + for (size_t i = 0; i < n; i++) { + const NodeID id = ids[i]; + const uint64_t chunkId = id / ChunkSize; + if (chunkId == lastChunkId) { + out[i] = cachedNode(*entry, id); + continue; + } - uint16_t scaledOffset = groupPtr->chunkOffsets[chunkOffset]; - ChunkInfoBase* basePtr = (ChunkInfoBase*)(((char *)(groupPtr->chunkOffsets + popcnt(groupPtr->chunkMask, 32))) + (scaledOffset * ChunkAlignment)); + if ((chunkId >> 32) != 0) { + lastChunkId = ~0ull; + out[i] = at(id); + continue; + } - if (basePtr->flags & ChunkCompressed) { - CompressedChunkInfo* ptr = (CompressedChunkInfo*)basePtr; - size_t latpSize = (ptr->flags >> 10) & ((1 << 10) - 1); - // TODO: we don't actually need the lonSize to decompress the data. - // May as well store it as a sanity check for now. - // size_t lonSize = ptr->flags & ((1 << 10) - 1); - size_t n = popcnt(ptr->nodeMask, 32) - 1; - - const size_t neededChunk = groupIndex * ChunkSize + chunk; - - // Keep a few recently decoded chunks per worker. Way geometry tends to - // revisit adjacent chunks, but not always the immediately preceding one. - ThreadStorage& tls = s(this); - size_t orderIndex = ChunkCacheSize; - for (size_t i = 0; i < ChunkCacheSize; ++i) { - if (tls.cacheChunks[tls.cacheOrder[i]].id == neededChunk) { - orderIndex = i; - break; + const uint64_t key = (generation << 32) | chunkId; + const size_t cacheIndex = chunkCacheIndex(key); + if (cache.keys[cacheIndex] != key) { + const ChunkInfoBase* basePtr = findChunk(groups.data(), id); + if (!(basePtr->flags & ChunkCompressed)) { + lastChunkId = ~0ull; + out[i] = uncompressedNode(basePtr, id); + continue; } + decodeChunk((const CompressedChunkInfo*)basePtr, cache.entries[cacheIndex]); + TMOPT_VERIFY_CHUNK((const CompressedChunkInfo*)basePtr, cache.entries[cacheIndex]); + cache.keys[cacheIndex] = key; } + entry = &cache.entries[cacheIndex]; + lastChunkId = chunkId; + out[i] = cachedNode(*entry, id); + } - size_t cacheSlot; - if (orderIndex == ChunkCacheSize) { - cacheSlot = tls.cacheOrder[ChunkCacheSize - 1]; - CachedChunk& cached = tls.cacheChunks[cacheSlot]; - cached.id = neededChunk; - - uint8_t* latpData = ptr->data; - uint8_t* lonData = ptr->data + latpSize; - uint32_t recovdata[256] = {0}; - - streamvbyte_decode(latpData, recovdata, n); - cached.latps[0] = ptr->firstLatp; - zigzag_delta_decode(recovdata, &cached.latps[1], n, cached.latps[0]); - - streamvbyte_decode(lonData, recovdata, n); - cached.lons[0] = ptr->firstLon; - zigzag_delta_decode(recovdata, &cached.lons[1], n, cached.lons[0]); - - promoteCacheSlot(tls, ChunkCacheSize - 1); - } else { - cacheSlot = tls.cacheOrder[orderIndex]; - if (orderIndex != 0) - promoteCacheSlot(tls, orderIndex); +#ifdef TMOPT_VERIFY + for (size_t i = 0; i < n; i++) { + const LatpLon expected = at(ids[i]); + if (!(expected == out[i])) { + std::cerr << "TMOPT_VERIFY: SortedNodeStore batched at mismatch for node " << ids[i] << std::endl; + abort(); } - - size_t nodeOffset = 0; - nodeOffset = popcnt(ptr->nodeMask, nodeMaskByte); - uint8_t maskByte = ptr->nodeMask[nodeMaskByte]; - maskByte = maskByte & ((1 << nodeMaskBit) - 1); - nodeOffset += popcnt(&maskByte, 1); - if (!(ptr->nodeMask[nodeMaskByte] & (1 << nodeMaskBit))) - throw std::out_of_range("SortedNodeStore: node " + std::to_string(id) + " missing, no node"); - - const CachedChunk& cached = tls.cacheChunks[cacheSlot]; - return { cached.latps[nodeOffset], cached.lons[nodeOffset] }; } - - UncompressedChunkInfo* ptr = (UncompressedChunkInfo*)basePtr; - size_t nodeOffset = 0; - nodeOffset = popcnt(ptr->nodeMask, nodeMaskByte); - uint8_t maskByte = ptr->nodeMask[nodeMaskByte]; - maskByte = maskByte & ((1 << nodeMaskBit) - 1); - nodeOffset += popcnt(&maskByte, 1); - if (!(ptr->nodeMask[nodeMaskByte] & (1 << nodeMaskBit))) - throw std::out_of_range("SortedNodeStore: node " + std::to_string(id) + " missing, no node"); - - return ptr->nodes[nodeOffset]; +#endif } size_t SortedNodeStore::size() const { diff --git a/src/sorted_way_store.cpp b/src/sorted_way_store.cpp index daba19c3..02f66147 100644 --- a/src/sorted_way_store.cpp +++ b/src/sorted_way_store.cpp @@ -202,10 +202,8 @@ void SortedWayStore::at(WayID id, std::vector& rv) const { ThreadStorage& tls = s(this); SortedWayStore::decodeWay(wayPtr->flags, wayPtr->data, tls.decodedWay); const std::vector& nodes = tls.decodedWay; - rv.clear(); - rv.reserve(nodes.size()); - for (const NodeID& node : nodes) - rv.push_back(nodeStore.at(node)); + rv.resize(nodes.size()); + nodeStore.at(nodes.data(), nodes.size(), rv.data()); } void SortedWayStore::insertLatpLons(std::vector &newWays) { diff --git a/src/tag_map.cpp b/src/tag_map.cpp index c3424fb7..cc63958f 100644 --- a/src/tag_map.cpp +++ b/src/tag_map.cpp @@ -1,95 +1,84 @@ #include "tag_map.h" #include #include +#include TagMap::TagMap() { - keys.resize(16); - key2value.resize(16); - values.resize(16); } void TagMap::reset() { - for (int i = 0; i < 16; i++) { - keys[i].clear(); - key2value[i].clear(); - values[i].clear(); - } + // Only clear shards that were used since the last reset. + for (unsigned i = 0, mask = usedKeyShards; mask; i++, mask >>= 1) + if (mask & 1) + keys[i].clear(); + for (unsigned i = 0, mask = usedValueShards; mask; i++, mask >>= 1) + if (mask & 1) + values[i].clear(); + usedKeyShards = 0; + usedValueShards = 0; } bool TagMap::empty() const { - for (int i = 0; i < keys.size(); i++) - if (keys[i].size() > 0) - return false; - - return true; -} -const std::size_t hashString(const std::string& str) { - // This is a pretty crappy hash function in terms of bit - // avalanching and distribution of output values. - // - // But it's very good in terms of speed, which turns out - // to be the important measure. - std::size_t hash = str.size(); - if (hash >= 4) - hash ^= *(uint32_t*)str.data(); - - return hash; + return usedKeyShards == 0; } -const std::size_t hashString(const char* str, size_t size) { +static inline std::size_t hashString(const char* str, size_t size) { // This is a pretty crappy hash function in terms of bit // avalanching and distribution of output values. // // But it's very good in terms of speed, which turns out // to be the important measure. std::size_t hash = size; - if (hash >= 4) - hash ^= *(uint32_t*)str; + if (hash >= 4) { + uint32_t word; + memcpy(&word, str, sizeof(word)); + hash ^= word; + } return hash; } -uint32_t TagMap::ensureString( - std::vector>& vector, - const protozero::data_view& value -) { - std::size_t hash = hashString(value.data(), value.size()); - - const uint16_t shard = hash % vector.size(); - for (int i = 0; i < vector[shard].size(); i++) - if (*(vector[shard][i]) == value) +uint32_t TagMap::ensureKey(const protozero::data_view& key) { + const uint16_t shard = hashString(key.data(), key.size()) & (SHARDS - 1); + std::vector& bucket = keys[shard]; + const uint32_t size = key.size(); + for (size_t i = 0; i < bucket.size(); i++) + if (bucket[i].size == size && memcmp(bucket[i].data, key.data(), size) == 0) return shard << 16 | i; - vector[shard].push_back(&value); - return shard << 16 | (vector[shard].size() - 1); + bucket.push_back(KeySlot{key.data(), size, 0}); + usedKeyShards |= 1u << shard; + return shard << 16 | (bucket.size() - 1); } +uint32_t TagMap::ensureValue(const protozero::data_view& value) { + const uint16_t shard = hashString(value.data(), value.size()) & (SHARDS - 1); + std::vector& bucket = values[shard]; + for (size_t i = 0; i < bucket.size(); i++) + if (bucket[i] == value) + return shard << 16 | i; -void TagMap::addTag(const protozero::data_view& key, const protozero::data_view& value) { - uint32_t valueLoc = ensureString(values, value); - uint32_t keyLoc = ensureString(keys, key); - + bucket.push_back(value); + usedValueShards |= 1u << shard; + return shard << 16 | (bucket.size() - 1); +} - const uint16_t shard = keyLoc >> 16; - const uint16_t pos = keyLoc; - if (key2value[shard].size() <= pos) { - key2value[shard].resize(pos + 1); - } +void TagMap::addTag(const protozero::data_view& key, const protozero::data_view& value) { + uint32_t valueLoc = ensureValue(value); + uint32_t keyLoc = ensureKey(key); - key2value[shard][pos] = valueLoc; + keys[keyLoc >> 16][keyLoc & 0xFFFF].valueLoc = valueLoc; } int64_t TagMap::getKey(const char* key, size_t size) const { // Return -1 if key not found, else return its keyLoc. - std::size_t hash = hashString(key, size); - - const uint16_t shard = hash % keys.size(); - for (int i = 0; i < keys[shard].size(); i++) { - const protozero::data_view& candidate = *keys[shard][i]; - if (candidate.size() != size) + const uint16_t shard = hashString(key, size) & (SHARDS - 1); + const std::vector& bucket = keys[shard]; + for (size_t i = 0; i < bucket.size(); i++) { + if (bucket[i].size != size) continue; - if (memcmp(candidate.data(), key, size) == 0) + if (memcmp(bucket[i].data, key, size) == 0) return shard << 16 | i; } @@ -98,15 +87,13 @@ int64_t TagMap::getKey(const char* key, size_t size) const { int64_t TagMap::getValue(const char* value, size_t size) const { // Return -1 if value not found, else return its valueLoc. - std::size_t hash = hashString(value, size); - - const uint16_t shard = hash % values.size(); - for (int i = 0; i < values[shard].size(); i++) { - const protozero::data_view& candidate = *values[shard][i]; - if (candidate.size() != size) + const uint16_t shard = hashString(value, size) & (SHARDS - 1); + const std::vector& bucket = values[shard]; + for (size_t i = 0; i < bucket.size(); i++) { + if (bucket[i].size() != size) continue; - if (memcmp(candidate.data(), value, size) == 0) + if (memcmp(bucket[i].data(), value, size) == 0) return shard << 16 | i; } @@ -114,23 +101,22 @@ int64_t TagMap::getValue(const char* value, size_t size) const { } const protozero::data_view* TagMap::getValueFromKey(uint32_t keyLoc) const { - const uint32_t valueLoc = key2value[keyLoc >> 16][keyLoc & 0xFFFF]; - return values[valueLoc >> 16][valueLoc & 0xFFFF]; + const uint32_t valueLoc = keys[keyLoc >> 16][keyLoc & 0xFFFF].valueLoc; + return &values[valueLoc >> 16][valueLoc & 0xFFFF]; } const protozero::data_view* TagMap::getValue(uint32_t valueLoc) const { - return values[valueLoc >> 16][valueLoc & 0xFFFF]; + return &values[valueLoc >> 16][valueLoc & 0xFFFF]; } boost::container::flat_map TagMap::exportToBoostMap() const { boost::container::flat_map rv; - for (int i = 0; i < keys.size(); i++) { - for (int j = 0; j < keys[i].size(); j++) { - uint32_t valueLoc = key2value[i][j]; - auto key = *keys[i][j]; - auto value = *values[valueLoc >> 16][valueLoc & 0xFFFF]; - rv[std::string(key.data(), key.size())] = std::string(value.data(), value.size()); + for (size_t i = 0; i < keys.size(); i++) { + for (size_t j = 0; j < keys[i].size(); j++) { + const KeySlot& key = keys[i][j]; + const protozero::data_view& value = *getValue(key.valueLoc); + rv[std::string(key.data, key.size)] = std::string(value.data(), value.size()); } } @@ -165,9 +151,9 @@ void TagMap::Iterator::operator++() { } Tag TagMap::Iterator::operator*() const { - const uint32_t valueLoc = map.key2value[shard][offset]; + const KeySlot& key = map.keys[shard][offset]; return Tag{ - *map.keys[shard][offset], - *map.getValue(valueLoc) + protozero::data_view(key.data, key.size), + *map.getValue(key.valueLoc) }; } diff --git a/src/tile_coordinates_set.cpp b/src/tile_coordinates_set.cpp index 2f63e468..462505aa 100644 --- a/src/tile_coordinates_set.cpp +++ b/src/tile_coordinates_set.cpp @@ -1,37 +1,142 @@ #include "tile_coordinates_set.h" #include #include +#include +#include +#include +#ifdef TMOPT_VERIFY +#include +#endif +#ifdef _MSC_VER +#include +#endif + +static inline unsigned int countTrailingZeros(uint64_t v) { +#if defined(__GNUC__) || defined(__clang__) + return __builtin_ctzll(v); +#elif defined(_MSC_VER) && defined(_M_X64) + unsigned long idx; + _BitScanForward64(&idx, v); + return idx; +#else + unsigned int n = 0; + while (!(v & 1)) { v >>= 1; n++; } + return n; +#endif +} + +static inline unsigned int popCount(uint64_t v) { +#if defined(__GNUC__) || defined(__clang__) + return __builtin_popcountll(v); +#else + unsigned int n = 0; + while (v) { v &= v - 1; n++; } + return n; +#endif +} + +void TileCoordinatesSet::forEachSet(const std::function& f) const { + const TileCoordinate dim = 1u << zoom(); + for (TileCoordinate x = 0; x < dim; x++) + for (TileCoordinate y = 0; y < dim; y++) + if (test(x, y)) + f(x, y); +} PreciseTileCoordinatesSet::PreciseTileCoordinatesSet(unsigned int zoom): zoom_(zoom), - tiles((1 << zoom) * (1 << zoom)) {} + numBits((uint64_t)(1 << zoom) * (1 << zoom)), + numWords((numBits + 63) / 64), + tiles(nullptr), + summary(((numWords + 63) / 64 + 63) / 64) { + // calloc: untouched pages stay unmapped (a z14 set is 32 MB). + tiles = static_cast(std::calloc(numWords, sizeof(uint64_t))); + if (tiles == nullptr) + throw std::bad_alloc(); +} + +PreciseTileCoordinatesSet::~PreciseTileCoordinatesSet() { + std::free(tiles); +} bool PreciseTileCoordinatesSet::test(TileCoordinate x, TileCoordinate y) const { uint64_t loc = x * (1 << zoom_) + y; - if (loc >= tiles.size()) + if (loc >= numBits) return false; - return tiles[loc]; + return (tiles[loc >> 6] >> (loc & 63)) & 1; } size_t PreciseTileCoordinatesSet::zoom() const { return zoom_; } +// Call f(w) for each word index w in a block marked in the summary, ascending. +template +static inline void forEachMarkedWord(const std::vector& summary, size_t numWords, F f) { + for (size_t s = 0; s < summary.size(); s++) { + uint64_t bits = summary[s]; + while (bits) { + const size_t block = s * 64 + countTrailingZeros(bits); + bits &= bits - 1; + const size_t wEnd = std::min((block + 1) * 64, numWords); + for (size_t w = block * 64; w < wEnd; w++) + f(w); + } + } +} + size_t PreciseTileCoordinatesSet::size() const { size_t rv = 0; - for (int i = 0; i < tiles.size(); i++) - if (tiles[i]) - rv++; - + forEachMarkedWord(summary, numWords, [&](size_t w) { rv += popCount(tiles[w]); }); return rv; } void PreciseTileCoordinatesSet::set(TileCoordinate x, TileCoordinate y) { uint64_t loc = x * (1 << zoom_) + y; - if (loc >= tiles.size()) + if (loc >= numBits) return; - tiles[loc] = true; + const uint64_t w = loc >> 6; + tiles[w] |= (uint64_t)1 << (loc & 63); + // 1 summary bit per block of 64 words + summary[w >> 12] |= (uint64_t)1 << ((w >> 6) & 63); +} + +void PreciseTileCoordinatesSet::forEachSet(const std::function& f) const { + // loc = x * 2^zoom + y, so ascending loc order is x-major, y-minor. + // Only blocks marked in the summary can hold set bits; visiting them in + // ascending order gives the same sequence as a full scan. +#ifdef TMOPT_VERIFY + { + std::vector full, fast; + for (size_t w = 0; w < numWords; w++) { + uint64_t word = tiles[w]; + while (word) { + full.push_back((uint64_t)w * 64 + countTrailingZeros(word)); + word &= word - 1; + } + } + forEachMarkedWord(summary, numWords, [&](size_t w) { + uint64_t word = tiles[w]; + while (word) { + fast.push_back((uint64_t)w * 64 + countTrailingZeros(word)); + word &= word - 1; + } + }); + if (full != fast || full.size() != size()) { + std::fprintf(stderr, "TMOPT_VERIFY: PreciseTileCoordinatesSet z%u mismatch %zu vs %zu\n", zoom_, full.size(), fast.size()); + std::abort(); + } + } +#endif + forEachMarkedWord(summary, numWords, [&](size_t w) { + uint64_t word = tiles[w]; + while (word) { + uint64_t loc = (uint64_t)w * 64 + countTrailingZeros(word); + word &= word - 1; + f((TileCoordinate)(loc >> zoom_), (TileCoordinate)(loc & ((1ull << zoom_) - 1))); + } + }); } LossyTileCoordinatesSet::LossyTileCoordinatesSet(unsigned int zoom, const TileCoordinatesSet& underlying) : zoom_(zoom), tiles(underlying), scale(1 << (zoom - underlying.zoom())) { diff --git a/src/tile_data.cpp b/src/tile_data.cpp index 727756ee..e34251f8 100644 --- a/src/tile_data.cpp +++ b/src/tile_data.cpp @@ -1,4 +1,6 @@ #include +#include +#include #include #include #include "tile_data.h" @@ -83,6 +85,11 @@ TileDataSource::TileDataSource(size_t threadNum, unsigned int indexZoom, bool in } thread_local std::vector>* tlsPendingSmallIndexObjects = nullptr; +thread_local uint64_t clipCacheBatchId = 0; +thread_local bool clipCacheNoStore = false; +#ifdef TMOPT_VERIFY +static std::atomic anySourceFinalized(false); +#endif void TileDataSource::finalize(size_t threadNum) { uint64_t finalized = 0; @@ -93,6 +100,15 @@ void TileDataSource::finalize(size_t threadNum) { } } + // Nothing adds to the small index after finalize (shapefiles, GeoJSON, + // Lua and declutter all run before), so release the pending storage. + // TLS pointers of other (finished) threads now dangle; they are never used again. + std::deque>>().swap(pendingSmallIndexObjects); + tlsPendingSmallIndexObjects = nullptr; +#ifdef TMOPT_VERIFY + anySourceFinalized = true; +#endif + std::cout << "indexed " << finalized << " contended objects" << std::endl; finalizeObjects(name(), threadNum, indexZoom, objects.begin(), objects.end(), lowZoomObjects); @@ -112,6 +128,13 @@ void TileDataSource::addObjectToSmallIndex(const TileCoordinates& index, const O const size_t z6index = z6x * CLUSTER_ZOOM_WIDTH + z6y; auto& mutex = objectsMutex[z6index % objectsMutex.size()]; +#ifdef TMOPT_VERIFY + if (anySourceFinalized) { + std::cerr << "TMOPT_VERIFY: addObjectToSmallIndex called after finalize" << std::endl; + abort(); + } +#endif + if (mutex.try_lock()) { addObjectToSmallIndexUnsafe(index, oo, id); mutex.unlock(); @@ -314,7 +337,9 @@ Geometry TileDataSource::buildWayGeometry(OutputGeometryType const geomType, case POLYGON_: { // Look for a previously clipped version at z-1, z-2, ... - std::shared_ptr cachedClip = multiPolygonClipCache.get(bbox.zoom, bbox.index.x, bbox.index.y, objectID); + auto cachedEntry = multiPolygonClipCache.getEntry(bbox.zoom, bbox.index.x, bbox.index.y, objectID); + std::shared_ptr cachedClip = cachedEntry.geom; + bool inputKnownValid = cachedClip == nullptr ? polygonKnownValid(objectID) : cachedEntry.valid; MultiPolygon uncached; @@ -375,10 +400,21 @@ Geometry TileDataSource::buildWayGeometry(OutputGeometryType const geomType, mp = std::move(uncached); else geom::assign(mp, input); - fast_clip(mp, box); + bool modified = fast_clip(mp, box); + if (inputKnownValid && !modified) { + // fast_clip was the identity on a geometry that already passed + // is_valid, so correct() is a no-op and it is still valid. + multiPolygonClipCache.add(bbox, objectID, mp, true); + return mp; + } geom::correct(mp); geom::validity_failure_type failure = geom::validity_failure_type::no_failure; - bool valid = geom::is_valid(mp,failure); + bool valid = multiPolygonDefinitelyValid(mp); +#ifdef TMOPT_VERIFY + precheckVerify(mp, valid, PRECHECK_TDS); +#endif + if (!valid) + valid = geom::is_valid(mp,failure); if (!valid) { if (verbose) { // Build the whole line first and emit it with a single stream write: @@ -422,7 +458,7 @@ Geometry TileDataSource::buildWayGeometry(OutputGeometryType const geomType, } } - multiPolygonClipCache.add(bbox, objectID, mp); + multiPolygonClipCache.add(bbox, objectID, mp, valid); return mp; } @@ -492,19 +528,45 @@ std::vector TileDataSource::getObjectsForTile( // ------------------------------------ // Add geometries to tile/large indices +namespace { +// Per-thread scratch for tile lists: clear() keeps the capacity, very large buffers are released. +thread_local std::vector tlsIndexTiles, tlsIndexTilesTmp; +inline void releaseIfHuge(std::vector& v) { + if (v.capacity() > (1u << 16)) std::vector().swap(v); +} + +#ifdef TMOPT_VERIFY +void verifyIndexTiles(const unordered_set& expected, const std::vector& got, const char* where) { + bool ok = expected.size() == got.size() && + std::is_sorted(got.begin(), got.end(), TileCoordinatesCompare()) && + std::adjacent_find(got.begin(), got.end()) == got.end(); + for (const auto& t : got) if (!ok || !expected.count(t)) { ok = false; break; } + if (!ok) { + std::cerr << "TMOPT_VERIFY: tile list mismatch in " << where << " (expected " << expected.size() << ", got " << got.size() << ")" << std::endl; + abort(); + } +} +#endif +} + void TileDataSource::addGeometryToIndex( const Linestring& geom, const std::vector& outputs, const uint64_t id ) { - unordered_set tileSet; + std::vector& tiles = tlsIndexTiles; + tiles.clear(); try { - insertIntermediateTiles(geom, indexZoom, tileSet); + insertIntermediateTiles(geom, indexZoom, tiles); +#ifdef TMOPT_VERIFY + unordered_set verifySet; + insertIntermediateTiles(geom, indexZoom, verifySet); + verifyIndexTiles(verifySet, tiles, "Linestring"); +#endif bool polygonExists = false; TileCoordinate minTileX = std::numeric_limits::max(), maxTileX = 0, minTileY = std::numeric_limits::max(), maxTileY = 0; - for (auto it = tileSet.begin(); it != tileSet.end(); ++it) { - TileCoordinates index = *it; + for (const TileCoordinates& index : tiles) { minTileX = std::min(index.x, minTileX); minTileY = std::min(index.y, minTileY); maxTileX = std::max(index.x, maxTileX); @@ -531,9 +593,14 @@ void TileDataSource::addGeometryToIndex( addObjectToLargeIndex(box, output, id); } else { // Smaller objects - add to each individual tile index - if (!tilesetFilled) { fillCoveredTiles(tileSet); tilesetFilled = true; } - for (auto it = tileSet.begin(); it != tileSet.end(); ++it) { - TileCoordinates index = *it; + if (!tilesetFilled) { + fillCoveredTiles(tiles); tilesetFilled = true; +#ifdef TMOPT_VERIFY + fillCoveredTiles(verifySet); + verifyIndexTiles(verifySet, tiles, "Linestring fill"); +#endif + } + for (const TileCoordinates& index : tiles) { addObjectToSmallIndex(index, output, id); } } @@ -542,6 +609,7 @@ void TileDataSource::addGeometryToIndex( } catch(std::out_of_range &err) { cerr << "Error calculating intermediate tiles: " << err.what() << endl; } + releaseIfHuge(tiles); } void TileDataSource::addGeometryToIndex( @@ -549,16 +617,22 @@ void TileDataSource::addGeometryToIndex( const std::vector& outputs, const uint64_t id ) { + std::vector& tiles = tlsIndexTiles; for (const auto& ls : geom) { - unordered_set tileSet; - insertIntermediateTiles(ls, indexZoom, tileSet); - for (auto it = tileSet.begin(); it != tileSet.end(); ++it) { - TileCoordinates index = *it; + tiles.clear(); + insertIntermediateTiles(ls, indexZoom, tiles); +#ifdef TMOPT_VERIFY + unordered_set verifySet; + insertIntermediateTiles(ls, indexZoom, verifySet); + verifyIndexTiles(verifySet, tiles, "MultiLinestring"); +#endif + for (const TileCoordinates& index : tiles) { for (const auto& output : outputs) { addObjectToSmallIndex(index, output, id); } } } + releaseIfHuge(tiles); } void TileDataSource::addGeometryToIndex( @@ -566,29 +640,40 @@ void TileDataSource::addGeometryToIndex( std::vector& outputs, const uint64_t id ) { - unordered_set tileSet; - bool singleOuter = geom.size()==1; + // Each outer's tiles are sorted, unique and filled; concatenate them, then sort+unique if there are several outers + std::vector& tiles = tlsIndexTiles; + std::vector& tilesTmp = tlsIndexTilesTmp; + tiles.clear(); for (const auto& poly : geom) { - unordered_set tileSetTmp; - insertIntermediateTiles(poly.outer(), indexZoom, tileSetTmp); - fillCoveredTiles(tileSetTmp); - if (singleOuter) { - tileSet = std::move(tileSetTmp); - } else { - tileSet.insert(tileSetTmp.begin(), tileSetTmp.end()); + tilesTmp.clear(); + insertIntermediateTiles(poly.outer(), indexZoom, tilesTmp); + fillCoveredTiles(tilesTmp); + tiles.insert(tiles.end(), tilesTmp.begin(), tilesTmp.end()); + } + if (geom.size() > 1) sortUniqueTiles(tiles); + releaseIfHuge(tilesTmp); +#ifdef TMOPT_VERIFY + { + unordered_set verifySet; + for (const auto& poly : geom) { + unordered_set tileSetTmp; + insertIntermediateTiles(poly.outer(), indexZoom, tileSetTmp); + fillCoveredTiles(tileSetTmp); + verifySet.insert(tileSetTmp.begin(), tileSetTmp.end()); } + verifyIndexTiles(verifySet, tiles, "MultiPolygon"); } +#endif TileCoordinate minTileX = std::numeric_limits::max(), maxTileX = 0, minTileY = std::numeric_limits::max(), maxTileY = 0; - for (auto it = tileSet.begin(); it != tileSet.end(); ++it) { - TileCoordinates index = *it; + for (const TileCoordinates& index : tiles) { minTileX = std::min(index.x, minTileX); minTileY = std::min(index.y, minTileY); maxTileX = std::max(index.x, maxTileX); maxTileY = std::max(index.y, maxTileY); } - const size_t tileSetSize = tileSet.size(); + const size_t tileSetSize = tiles.size(); for (auto& output : outputs) { if (tileSetSize >= 16) { // Larger objects - add to rtree @@ -599,12 +684,12 @@ void TileDataSource::addGeometryToIndex( addObjectToLargeIndex(box, output, id); } else { // Smaller objects - add to each individual tile index - for (auto it = tileSet.begin(); it != tileSet.end(); ++it) { - TileCoordinates index = *it; + for (const TileCoordinates& index : tiles) { addObjectToSmallIndex(index, output, id); } } } + releaseIfHuge(tiles); } NodeID TileDataSource::storePoint(const Point& input) { diff --git a/src/tile_sorting.cpp b/src/tile_sorting.cpp index 1d9dc58a..0d5cb839 100644 --- a/src/tile_sorting.cpp +++ b/src/tile_sorting.cpp @@ -1,9 +1,98 @@ #include #include #include +#include +#include +#include +#include +#include #include "tile_data_base.h" #include "append_vector.h" #include +#include +#include +#include +#include +#include +#include + +static_assert(sizeof(Z6Offset) <= 2, "interleaveKey assumes Z6Offset has at most 16 bits"); +static const unsigned Z6OffsetBits = 8 * sizeof(Z6Offset); + +// Spread the low 16 bits of v so that bit i moves to bit 2*i. +static inline uint32_t spreadBits16(uint32_t v) { + v &= 0x0000FFFF; + v = (v | (v << 8)) & 0x00FF00FF; + v = (v | (v << 4)) & 0x0F0F0F0F; + v = (v | (v << 2)) & 0x33333333; + v = (v | (v << 1)) & 0x55555555; + return v; +} + +// Bit-interleave the low B bits of x and y, MSB first, with the x bit above +// the y bit at each level. Bits of x and y at or above B are kept as a +// prefix (all of x's high bits, then all of y's high bits), so the key has +// 2 * Z6OffsetBits bits. Requires x, y < 2^Z6OffsetBits and B <= Z6OffsetBits. +static inline uint32_t interleaveKey(uint32_t x, uint32_t y, unsigned B) { + const uint32_t lowMask = (1u << B) - 1; + const unsigned highBits = Z6OffsetBits - B; + const uint32_t low = (spreadBits16(x & lowMask) << 1) | spreadBits16(y & lowMask); + return (x >> B) << (highBits + 2 * B) | + (y >> B) << (2 * B) | + low; +} + +// The original ordering: compare x then y at each zoom from CLUSTER_ZOOM to indexZoom. +static inline bool zoomLoopLess(size_t aX, size_t aY, size_t bX, size_t bY, unsigned int indexZoom) { + for (size_t z = CLUSTER_ZOOM; z <= indexZoom; z++) { + const auto aXz = aX / (1 << (indexZoom - z)); + const auto bXz = bX / (1 << (indexZoom - z)); + if (aXz != bXz) + return aXz < bXz; + + const auto aYz = aY / (1 << (indexZoom - z)); + const auto bYz = bY / (1 << (indexZoom - z)); + + if (aYz != bYz) + return aYz < bYz; + } + + return false; +} + +#ifndef NDEBUG +// Debug-only: verify that comparing interleaved keys matches the zoom-loop comparator. +static void checkInterleaveKey() { + static std::once_flag once; + std::call_once(once, []() { + std::mt19937 rng(12345); + const uint32_t maxVal = (1u << Z6OffsetBits) - 1; + std::uniform_int_distribution dist(0, maxVal); + for (unsigned indexZoom = CLUSTER_ZOOM; indexZoom <= CLUSTER_ZOOM + Z6OffsetBits; indexZoom++) { + const unsigned B = indexZoom - CLUSTER_ZOOM; + for (int i = 0; i < 200000; i++) { + uint32_t aX = dist(rng), aY = dist(rng), bX = dist(rng), bY = dist(rng); + // Also exercise near-equal pairs + if (i & 1) { bX = aX ^ (1u << (rng() % Z6OffsetBits)); bY = (i & 2) ? aY : bY; } + if ((i % 7) == 0) { bX = aX; bY = aY ^ (1u << (rng() % Z6OffsetBits)); } + bX &= maxVal; bY &= maxVal; + const uint32_t ka = interleaveKey(aX, aY, B), kb = interleaveKey(bX, bY, B); + if (zoomLoopLess(aX, aY, bX, bY, indexZoom) != (ka < kb) || + zoomLoopLess(bX, bY, aX, aY, indexZoom) != (kb < ka)) { + std::cerr << "interleaveKey mismatch at indexZoom " << indexZoom << std::endl; + assert(false); + } + } + } + }); +} +#endif + +// The interleaved key has 2 * Z6OffsetBits bits; it is sorted by a single +// counting-sort pass over all possible key values. +static const size_t keyValues = size_t(1) << (2 * Z6OffsetBits); +// Below this, a comparison sort on the key is cheaper than counting. +static const size_t countingSortMinSize = 65536; template void sortOutputObjects( const unsigned int indexZoom, @@ -30,33 +119,73 @@ template void sortOutputObjects( // better to assign chunks of `objects` to each thread. // // That's a future performance improvement, so deferring for now. - boost::sort::block_indirect_sort( - begin, end, - [indexZoom](const OO& a, const OO& b) { - // Cluster by parent zoom, so that a subsequent search - // can find a contiguous range of entries for any tile - // at zoom 6 or higher. - const size_t aX = a.x; - const size_t aY = a.y; - const size_t bX = b.x; - const size_t bY = b.y; - for (size_t z = CLUSTER_ZOOM; z <= indexZoom; z++) { - const auto aXz = aX / (1 << (indexZoom - z)); - const auto bXz = bX / (1 << (indexZoom - z)); - if (aXz != bXz) - return aXz < bXz; - - const auto aYz = aY / (1 << (indexZoom - z)); - const auto bYz = bY / (1 << (indexZoom - z)); - - if (aYz != bYz) - return aYz < bYz; - } - - return false; - }, - threadNum - ); + // + // Cluster by parent zoom, so that a subsequent search can find a + // contiguous range of entries for any tile at zoom 6 or higher. + // i.e. compare x then y at each zoom from CLUSTER_ZOOM to indexZoom, + // which is the same as comparing the bit-interleaved (x, y) key. + if (indexZoom < CLUSTER_ZOOM) + return; // the zoom-loop comparator treats all objects as equal + + // Levels with a shift of Z6OffsetBits or more compare 0 with 0. + const unsigned B = std::min(indexZoom - CLUSTER_ZOOM, Z6OffsetBits); +#ifndef NDEBUG + checkInterleaveKey(); +#endif + auto key = [B](const OO& o) -> uint32_t { return interleaveKey(o.x, o.y, B); }; + const size_t n = end - begin; + + if (n < countingSortMinSize) { + boost::sort::pdqsort(begin, end, [&key](const OO& a, const OO& b) { return key(a) < key(b); }); + return; + } + + // Parallel counting sort on the key: each chunk copies its objects + // into a contiguous buffer and counts its keys, then scatters them + // back into their final positions. + const size_t chunks = std::max(1, std::min(threadNum, n / countingSortMinSize)); + const size_t chunkSize = (n + chunks - 1) / chunks; + // OO has no default constructor, so use raw storage for the copies. + static_assert(std::is_trivially_copyable::value && std::is_trivially_destructible::value, + "the counting sort buffer copies objects as raw storage"); + std::unique_ptr::type[]> storage( + new typename std::aligned_storage::type[n]); + OO* buffer = reinterpret_cast(storage.get()); + std::vector offsets(chunks * keyValues, 0); + + auto runChunks = [&](const std::function& fn) { + boost::asio::thread_pool pool(chunks); + for (size_t c = 0; c < chunks; c++) { + const size_t from = c * chunkSize, to = std::min(n, from + chunkSize); + boost::asio::post(pool, [&fn, c, from, to]() { fn(c, from, to); }); + } + pool.join(); + }; + + runChunks([&](size_t c, size_t from, size_t to) { + size_t* counts = &offsets[c * keyValues]; + auto it = begin + from; + for (size_t i = from; i < to; i++, ++it) { + new (&buffer[i]) OO(*it); + counts[key(buffer[i])]++; + } + }); + + // Turn counts into starting positions: by key, then by chunk. + size_t pos = 0; + for (size_t k = 0; k < keyValues; k++) + for (size_t c = 0; c < chunks; c++) { + size_t& slot = offsets[c * keyValues + k]; + const size_t count = slot; + slot = pos; + pos += count; + } + + runChunks([&](size_t c, size_t from, size_t to) { + size_t* next = &offsets[c * keyValues]; + for (size_t i = from; i < to; i++) + begin[next[key(buffer[i])]++] = buffer[i]; + }); } template void sortOutputObjects( @@ -127,7 +256,10 @@ void sortTileCoordinates( return aY < bY; } - for (size_t z = CLUSTER_ZOOM; z <= baseZoom; z++) { + // Compare beyond baseZoom too, so tiles above baseZoom get a total, + // depth-first order (parents before children) rather than ties. + const size_t maxZoom = std::max(baseZoom, std::max(aZoom, bZoom)); + for (size_t z = CLUSTER_ZOOM; z <= maxZoom; z++) { // Translate both a and b to zoom z, compare. // First, sanity check: can we translate it to this zoom? if (aZoom < z || bZoom < z) { diff --git a/src/tile_worker.cpp b/src/tile_worker.cpp index 891dd00f..614ed30c 100644 --- a/src/tile_worker.cpp +++ b/src/tile_worker.cpp @@ -128,7 +128,8 @@ void writeMultiLinestring( unsigned zoom, double simplifyLevel, unsigned simplifyAlgo, - const MultiLinestring& mls + const MultiLinestring& mls, + MvtAttributeMemo* attrMemo ) { typename Backend::LinestringFeature fbuilder{vtLayer}; @@ -192,7 +193,7 @@ void writeMultiLinestring( if (hadLine) { // add the properties - oo.oo.writeAttributes(attributeStore, fbuilder, zoom); + oo.oo.writeAttributes(attributeStore, fbuilder, zoom, attrMemo); // call commit() when you are done fbuilder.commit(); } @@ -248,30 +249,65 @@ void writeMultiPolygon( unsigned zoom, double simplifyLevel, unsigned simplifyAlgo, - const MultiPolygon& mp + const MultiPolygon& mp, + MvtAttributeMemo* attrMemo ) { bbox.scaleGeometry(scaledMultiPolygon, mp); MultiPolygon ¤t = scaledMultiPolygon; + // Set when every polygon of `current` is individually valid after the + // Visvalingam repair and remove_spikes changed nothing. Then the validity + // check/repair below is a no-op: repair_multi_polygon keeps individually valid + // polygons of a multipolygon as they are, and a single valid polygon is valid. + bool knownPolygonsValid = false; if (simplifyLevel>0) { if (simplifyAlgo == LayerDef::VISVALINGAM) { - current = simplifyVis(current, simplifyLevel/bbox.xscale); - } else if (simplifyAlgo == LayerDef::BUILDINGS) { - if (current.size() > 1 || boost::geometry::num_points(current)>5) { - geom::correct(current); // self-intersections can break simplification - simplifyBuildings(current, simplifyLevel/bbox.xscale); - } + bool allValid = false; + current = simplifyVis(current, simplifyLevel/bbox.xscale, allValid); + const std::size_t numPoints = boost::geometry::num_points(current); + geom::remove_spikes(current); + knownPolygonsValid = allValid && !verbose && boost::geometry::num_points(current) == numPoints; } else { - current = simplify(current, simplifyLevel/bbox.xscale); + if (simplifyAlgo == LayerDef::BUILDINGS) { + if (current.size() > 1 || boost::geometry::num_points(current)>5) { + geom::correct(current); // self-intersections can break simplification + simplifyBuildings(current, simplifyLevel/bbox.xscale); + } + } else { + current = simplify(current, simplifyLevel/bbox.xscale); + } + geom::remove_spikes(current); } - geom::remove_spikes(current); } if (geom::is_empty(current)) return; geom::correct(current); +#ifdef TMOPT_VERIFY + if (knownPolygonsValid) { + MultiPolygon check = current; + if (!geom::is_valid(check)) repair_multi_polygon(check, true); + bool same = check.size() == current.size(); + for (std::size_t i = 0; same && i < check.size(); ++i) { + const Polygon &a = check[i], &b = current[i]; + same = a.inners().size() == b.inners().size() && a.outer().size() == b.outer().size(); + for (std::size_t j = 0; same && j < a.outer().size(); ++j) + same = a.outer()[j].x() == b.outer()[j].x() && a.outer()[j].y() == b.outer()[j].y(); + for (std::size_t r = 0; same && r < a.inners().size(); ++r) { + same = a.inners()[r].size() == b.inners()[r].size(); + for (std::size_t j = 0; same && j < a.inners()[r].size(); ++j) + same = a.inners()[r][j].x() == b.inners()[r][j].x() && a.inners()[r][j].y() == b.inners()[r][j].y(); + } + } + if (!same) { + cerr << "TMOPT_VERIFY writeMultiPolygon skip mismatch" << endl; + abort(); + } + } +#endif + geom::validity_failure_type failure; - if (!geom::is_valid(current, failure)) { + if (!knownPolygonsValid && (simplifyLevel > 0 || verbose) && !geom::is_valid(current, failure)) { if (verbose) { cout << "output multipolygon has " << boost_validity_error(failure) << endl; @@ -285,7 +321,7 @@ void writeMultiPolygon( // Simplification can turn a valid input into a self-intersecting/spiky // one; such polygons are silently dropped by many renderers (missing // features). Repair (dissolve, then zero-width buffer) before writing. - bool repaired = repair_multi_polygon(current); + bool repaired = repair_multi_polygon(current, true); // already known invalid if (geom::is_empty(current)) return; @@ -321,7 +357,7 @@ void writeMultiPolygon( if (hadPoly) { // add the properties - oo.oo.writeAttributes(attributeStore, fbuilder, zoom); + oo.oo.writeAttributes(attributeStore, fbuilder, zoom, attrMemo); // call commit() when you are done fbuilder.commit(); } @@ -342,7 +378,8 @@ void ProcessObjects( bool combinePolygons, unsigned zoom, const TileBbox &bbox, - typename Backend::LayerBuilder& vtLayer + typename Backend::LayerBuilder& vtLayer, + MvtAttributeMemo* attrMemo ) { for (auto jt = ooSameLayerBegin; jt != ooSameLayerEnd; ++jt) { OutputObjectID oo = *jt; @@ -392,7 +429,7 @@ void ProcessObjects( for (const auto &point : multipoint) fbuilder.set_point(point.first, point.second); - oo.oo.writeAttributes(attributeStore, fbuilder, zoom); + oo.oo.writeAttributes(attributeStore, fbuilder, zoom, attrMemo); fbuilder.commit(); oo = *jt; @@ -438,9 +475,9 @@ void ProcessObjects( } if (oo.oo.geomType == LINESTRING_ || oo.oo.geomType == MULTILINESTRING_) - writeMultiLinestring(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, simplifyAlgo, boost::get(g)); + writeMultiLinestring(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, simplifyAlgo, boost::get(g), attrMemo); else if (oo.oo.geomType == POLYGON_) - writeMultiPolygon(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, simplifyAlgo, boost::get(g)); + writeMultiPolygon(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, simplifyAlgo, boost::get(g), attrMemo); } } } @@ -492,7 +529,8 @@ void ProcessLayer( typename Backend::LayerBuilder& vtLayer, const TileBbox& bbox, const std::vector& ltx, - SharedData& sharedData + SharedData& sharedData, + MvtAttributeMemo* attrMemo = nullptr ) { const std::string& layerName = sharedData.layers.layers[ltx.at(0)].name; @@ -528,7 +566,7 @@ void ProcessLayer( ProcessObjects(sources[i], attributeStore, ooListSameLayer.first, end, sharedData, simplifyLevel, ld.simplifyAlgo, - filterArea, ld.combinePoints, zoom < ld.combineLinesBelow, zoom < ld.combinePolygonsBelow, zoom, bbox, vtLayer); + filterArea, ld.combinePoints, zoom < ld.combineLinesBelow, zoom < ld.combinePolygonsBelow, zoom, bbox, vtLayer, attrMemo); } } if (verbose && std::time(0)-start>3) { @@ -594,14 +632,35 @@ void outputProc( } vtzero::vector_tile existingTile{rawExistingTile}; + thread_local MvtAttributeMemo attrMemo; for (auto lt = sharedData.layers.layerOrder.begin(); lt != sharedData.layers.layerOrder.end(); ++lt) { if (signalStop) break; const std::string& layerName = sharedData.layers.layers[lt->at(0)].name; vtzero::layer_builder vtLayer{tile, layerName, sharedData.config.mvtVersion, extent}; + attrMemo.reset(vtLayer, zoom); CopyExistingLayer(existingTile, layerName, vtLayer); - ProcessLayer(sources, attributeStore, coordinates, zoom, data, vtLayer, bbox, *lt, sharedData); + ProcessLayer(sources, attributeStore, coordinates, zoom, data, vtLayer, bbox, *lt, sharedData, &attrMemo); } tile.serialize(outputdata); +#ifdef TMOPT_VERIFY + { + // Encode the tile again without the attribute memo: the bytes must be identical + vtzero::tile_builder verifyTile; + for (auto lt = sharedData.layers.layerOrder.begin(); lt != sharedData.layers.layerOrder.end(); ++lt) { + if (signalStop) break; + const std::string& layerName = sharedData.layers.layers[lt->at(0)].name; + vtzero::layer_builder vtLayer{verifyTile, layerName, sharedData.config.mvtVersion, extent}; + CopyExistingLayer(existingTile, layerName, vtLayer); + ProcessLayer(sources, attributeStore, coordinates, zoom, data, vtLayer, bbox, *lt, sharedData, nullptr); + } + std::string verifyData; + verifyTile.serialize(verifyData); + if (verifyData != outputdata) { + cerr << "TMOPT_VERIFY attribute memo tile mismatch at " << zoom << "/" << bbox.index.x << "/" << bbox.index.y << endl; + abort(); + } + } +#endif } // Write to file or sqlite diff --git a/src/tilemaker.cpp b/src/tilemaker.cpp index be7dca11..7a7ae1e5 100644 --- a/src/tilemaker.cpp +++ b/src/tilemaker.cpp @@ -33,6 +33,10 @@ #ifndef _MSC_VER #include #endif +#ifdef __GLIBC__ +#include +#include +#endif #include "geom.h" #include "node_stores.h" @@ -90,6 +94,17 @@ void sortTileCoordinates( * Worker threads write the output tiles, and start in the outputProc function. */ int main(const int argc, const char* argv[]) { +#ifdef __GLIBC__ + // glibc raises its mmap threshold dynamically (up to 32/64 MB) once a large + // mmapped chunk is freed; multi-MB buffers then land in per-thread arenas + // that are never trimmed, inflating peak RSS. Setting the trim threshold + // explicitly also disables that dynamic adjustment (the mmap threshold + // stays at its 128 KB default), unless the user tunes malloc via the environment. + if (!getenv("MALLOC_TRIM_THRESHOLD_") && !getenv("MALLOC_MMAP_THRESHOLD_") && !getenv("GLIBC_TUNABLES")) { + mallopt(M_TRIM_THRESHOLD, 4 << 20); + } +#endif + // ---- Read command-line options OptionsParser::Options options; try { @@ -454,27 +469,22 @@ int main(const int argc, const char* argv[]) { const auto& zoomResult = zoomResults[zoom]; int numTiles = 0; - for (int x = 0; x < 1 << zoom; x++) { - for (int y = 0; y < 1 << zoom; y++) { - if (!zoomResult->test(x, y)) - continue; - - if (hasClippingBox) { - bool isInAWhollyCoveredZ6Tile = false; - if (zoom >= 6) { - TileCoordinate z6x = x / (1 << (zoom - 6)); - TileCoordinate z6y = y / (1 << (zoom - 6)); - isInAWhollyCoveredZ6Tile = coveredZ6Tiles.test(z6x, z6y); - } - - if(!isInAWhollyCoveredZ6Tile && !boost::geometry::intersects(TileBbox(TileCoordinates(x, y), zoom, false, false).getTileBox(), clippingBox)) - continue; + zoomResult->forEachSet([&](TileCoordinate x, TileCoordinate y) { + if (hasClippingBox) { + bool isInAWhollyCoveredZ6Tile = false; + if (zoom >= 6) { + TileCoordinate z6x = x / (1 << (zoom - 6)); + TileCoordinate z6y = y / (1 << (zoom - 6)); + isInAWhollyCoveredZ6Tile = coveredZ6Tiles.test(z6x, z6y); } - tileCoordinates.push_back(std::make_pair(zoom, TileCoordinates(x, y))); - numTiles++; + if(!isInAWhollyCoveredZ6Tile && !boost::geometry::intersects(TileBbox(TileCoordinates(x, y), zoom, false, false).getTileBox(), clippingBox)) + return; } - } + + tileCoordinates.push_back(std::make_pair(zoom, TileCoordinates(x, y))); + numTiles++; + }); std::cout << " (" << numTiles; #ifdef CLOCK_MONOTONIC @@ -492,14 +502,22 @@ int main(const int argc, const char* argv[]) { // Cluster tiles: breadth-first for z0..z5, depth-first for z6 sortTileCoordinates(config.baseZoom, options.threadNum, tileCoordinates); + std::atomic batchCounter(0); + const uint cutZoom = std::max(12, config.startZoom); std::size_t batchSize = 0; for(std::size_t startIndex = 0; startIndex < tileCoordinates.size(); startIndex += batchSize) { // Compute how many tiles should be assigned to this batch -- // higher-zoom tiles are cheaper to compute, lower-zoom tiles more expensive. batchSize = 0; size_t weight = 0; - while (weight < 1000 && startIndex + batchSize < tileCoordinates.size()) { + // Batches are a fixed partition: each z<=11 tile is its own batch, and + // z12 subtrees are grouped (~2 per batch), so z13/z14 tiles usually share + // a batch with their z12 ancestor (batch-local ClipCache can reuse its clips). + // The partition doesn't depend on thread count, so neither does the output. + while (startIndex + batchSize < tileCoordinates.size()) { const auto& zoom = tileCoordinates[startIndex + batchSize].first; + if (batchSize > 0 && (weight >= 1000 || (weight >= 50 && zoom <= cutZoom))) + break; if (zoom > 12) weight++; else if (zoom > 11) @@ -512,7 +530,10 @@ int main(const int argc, const char* argv[]) { batchSize++; } - boost::asio::post(pool, [=, &tileCoordinates, &pool, &sharedData, &sources, &attributeStore, &io_mutex, &tilesWritten, &lastTilesWritten]() { + boost::asio::post(pool, [=, &tileCoordinates, &pool, &sharedData, &sources, &attributeStore, &io_mutex, &tilesWritten, &lastTilesWritten, &batchCounter]() { + clipCacheBatchId = ++batchCounter; + // A single-tile batch has no descendants to reuse its clips. + clipCacheNoStore = batchSize == 1; std::vector tileTimings; std::size_t endIndex = std::min(tileCoordinates.size(), startIndex + batchSize); for(std::size_t i = startIndex; i < endIndex; ++i) { diff --git a/src/visvalingam.cpp b/src/visvalingam.cpp index f6b8b635..5e79e051 100644 --- a/src/visvalingam.cpp +++ b/src/visvalingam.cpp @@ -1,4 +1,5 @@ #include "geom.h" +#include "visvalingam.h" #include // Adapted from @@ -31,9 +32,15 @@ #include #include +#include +#ifdef TMOPT_VERIFY +#include +#include +#endif -// Stuff to create the priority queue, or min heap. -// Rewriting it here, vs using the std lib, resulted in a 50% performance bump! +#ifdef TMOPT_VERIFY +// The original pointer-based implementation, used to cross-check the indexed one. +namespace visOld { struct visItem { double area = 0; // triangle area @@ -144,6 +151,151 @@ struct minHeap { } }; +} +using visOld::visItem; +using visOld::minHeap; + +#endif + +// Stuff to create the priority queue, or min heap. +// Rewriting it here, vs using the std lib, resulted in a 50% performance bump! +// +// Items are referred to by index; the heap stores (area, item) pairs and +// heapPos[item] is the item's slot in the heap. The heap operations (and so +// the tie-breaking between equal areas) are the same as in the original +// pointer-based version. + +namespace { + +struct visHeapEntry { + double area; // triangle area + int item; // item index +}; + +struct visScratch { + std::vector h; + std::vector heapPos; // item -> index in h + std::vector next; // virtual linked list, -1 = none + std::vector prev; + + void prepare(size_t numItems) { + h.clear(); + h.reserve(numItems); + heapPos.resize(numItems); + next.resize(numItems); + prev.resize(numItems); + } + + // Don't keep huge buffers alive after an unusually large ring. + void trim() { + const size_t cap = 1 << 16; + if (h.capacity() > cap) { + std::vector().swap(h); + std::vector().swap(heapPos); + std::vector().swap(next); + std::vector().swap(prev); + } + } + + void Push(int item, double area) { + int i = h.size(); + heapPos[item] = i; + h.push_back(visHeapEntry{area, item}); + up(i); + } + + visHeapEntry Pop() { + visHeapEntry removed = h[0]; + visHeapEntry lastItem = h[h.size() - 1]; + h.pop_back(); + + if (h.size() > 0) { + heapPos[lastItem.item] = 0; + h[0] = lastItem; + down(0); + } + + return removed; + } + + void Update(int item, double area) { + int i = heapPos[item]; + if (h[i].area > area) { + // area got smaller + h[i].area = area; + up(i); + } else { + // area got larger + h[i].area = area; + down(i); + } + } + + void up(int i) { + visHeapEntry object = h[i]; + while (i > 0) { + int up = ((i + 1) >> 1) - 1; + const visHeapEntry parent = h[up]; + + if (parent.area <= object.area) { + // parent is smaller so we're done fixing up the heap. + break; + } + + // swap nodes + heapPos[parent.item] = i; + h[i] = parent; + + heapPos[object.item] = up; + h[up] = object; + + i = up; + } + } + + void down(int i) { + visHeapEntry object = h[i]; + const size_t size = h.size(); + while (1) { + size_t right = (i + 1) << 1; + size_t left = right - 1; + + int down = i; + double childArea = h[down].area; + + // swap with smallest child + if (left < size && h[left].area < childArea) { + down = left; + childArea = h[down].area; + } + + if (right < size && h[right].area < childArea) { + down = right; + childArea = h[down].area; + } + + // non smaller, so quit + if (down == i) { + break; + } + + // swap the nodes + const visHeapEntry child = h[down]; + heapPos[child.item] = i; + h[i] = child; + + heapPos[object.item] = down; + h[down] = object; + + i = down; + } + } +}; + +thread_local visScratch visTls; + +} + template static double doubleTriangleArea(GeometryType const &ls, int start, int i1, int i2, int i3) { Point const &a = ls[i1 + start]; @@ -153,8 +305,9 @@ static double doubleTriangleArea(GeometryType const &ls, int start, int i1, int return std::abs((b.x() - a.x()) * (c.y() - a.y()) - (b.y() - a.y()) * (c.x() - a.x())); } +#ifdef TMOPT_VERIFY template -GeometryType visvalingam(const GeometryType &ls, double threshold, size_t retain) { +GeometryType visvalingamOld(const GeometryType &ls, double threshold, size_t retain) { size_t start = 0; size_t end = ls.size(); // - 1; int removed = 0; @@ -245,6 +398,113 @@ GeometryType visvalingam(const GeometryType &ls, double threshold, size_t retain } return output; } +#endif + +template +GeometryType visvalingam(const GeometryType &ls, double threshold, size_t retain) { + size_t start = 0; + size_t end = ls.size(); // - 1; + int removed = 0; + threshold *= 2; + if (ls.empty()) return ls; + + const int n = end - start; + // Items 0..n-1 are the path points; item n is the list head, which + // (like the original's separate linkedListStart) refers to point 0. + const int head = n; + auto pointIndex = [head](int item) { return item == head ? 0 : item; }; + + visScratch &s = visTls; + s.prepare(n + 1); + std::vector &next = s.next; + std::vector &prev = s.prev; + + // build the initial minheap linked list. + prev[head] = -1; + s.Push(head, INFINITY); + + { + int previous = head; + for (int i = 1; i < n - 1; i++) { + prev[i] = previous; + s.Push(i, doubleTriangleArea(ls, start, i - 1, i, i + 1)); + next[previous] = i; + previous = i; + } + + // final item + int endItem = n - 1; + prev[endItem] = previous; + next[endItem] = -1; + + next[previous] = endItem; + s.Push(endItem, INFINITY); + } + + // run through the reduction process + while (s.h.size() > 0) { + visHeapEntry current = s.Pop(); + if (current.area > threshold) { + break; + } else if ((end - start) - removed <= retain) { + break; + } + + int nextItem = next[current.item]; + int previous = prev[current.item]; + + // remove current element from linked list + next[previous] = nextItem; + prev[nextItem] = previous; + removed++; + + // figure out the new areas + if (prev[previous] != -1) { + double area = doubleTriangleArea(ls, start, + pointIndex(prev[previous]), + pointIndex(previous), + pointIndex(nextItem)); + + area = std::max(area, current.area); + s.Update(previous, area); + } + + if (next[nextItem] != -1) { + double area = doubleTriangleArea(ls, start, + pointIndex(previous), + pointIndex(nextItem), + pointIndex(next[nextItem])); + + area = std::max(area, current.area); + s.Update(nextItem, area); + } + } + + GeometryType output; + output.reserve(end - start - removed); + int item = head; + while (item != -1) { + output.emplace_back(ls[pointIndex(item) + start]); + item = next[item]; + } + s.trim(); +#ifdef TMOPT_VERIFY + { + GeometryType expected = visvalingamOld(ls, threshold / 2, retain); + bool same = expected.size() == output.size(); + for (size_t k = 0; same && k < output.size(); k++) { + double a[2] = { output[k].x(), output[k].y() }; + double b[2] = { expected[k].x(), expected[k].y() }; + same = std::memcmp(a, b, sizeof(a)) == 0; + } + if (!same) { + std::cerr << "TMOPT_VERIFY: visvalingam mismatch, n=" << ls.size() << " out=" << output.size() << " expected=" << expected.size() << std::endl; + std::abort(); + } + } +#endif + return output; +} Linestring simplifyVis(const Linestring &ls, double max_distance) { if (ls.size()<3) return ls; @@ -263,6 +523,12 @@ Polygon simplifyVis(const Polygon &p, double max_distance) { return output; } MultiPolygon simplifyVis(const MultiPolygon &mp, double max_distance) { + bool allValid; + return simplifyVis(mp, max_distance, allValid); +} +// allValid: set to repair_multi_polygon's result, i.e. true when every polygon +// of the output is individually valid. +MultiPolygon simplifyVis(const MultiPolygon &mp, double max_distance, bool &allValid) { MultiPolygon output; output.reserve(mp.size()); for (const auto &p : mp) { @@ -274,6 +540,6 @@ MultiPolygon simplifyVis(const MultiPolygon &mp, double max_distance) { // almost all of its area (which showed up as missing lake tiles at low // zoom). repair_multi_polygon keeps the simplified geometry untouched if a // repair would not preserve the covered area. - repair_multi_polygon(output); + allValid = repair_multi_polygon(output); return output; }