From e83bc29f4d983721991fbcf7647f61977ca2e47b 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 1/3] 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; } From 9f48efdfa484f967ddabcf17de51db4c120e434a Mon Sep 17 00:00:00 2001 From: dktue <2878623+dktue@users.noreply.github.com> Date: Tue, 27 Oct 2026 13:56:48 +0100 Subject: [PATCH 2/3] chore: remove TMOPT_VERIFY cross-check code --- include/clip_cache.h | 39 ------ include/geom.h | 24 ---- include/lua_pool.h | 27 ---- include/osm_lua_processing.h | 16 +-- include/polylabel.h | 108 +--------------- src/attribute_store.cpp | 23 ---- src/geom.cpp | 53 -------- src/osm_lua_processing.cpp | 51 -------- src/osm_mem_tiles.cpp | 19 --- src/output_object.cpp | 24 ---- src/pbf_processor.cpp | 11 -- src/shp_processor.cpp | 3 - src/sorted_node_store.cpp | 56 --------- src/tile_coordinates_set.cpp | 26 ---- src/tile_data.cpp | 53 -------- src/tile_worker.cpp | 41 ------- src/visvalingam.cpp | 231 ----------------------------------- 17 files changed, 9 insertions(+), 796 deletions(-) diff --git a/include/clip_cache.h b/include/clip_cache.h index f4c16cfb..66a400ab 100644 --- a/include/clip_cache.h +++ b/include/clip_cache.h @@ -6,10 +6,6 @@ #include "geom.h" #include #include -#ifdef TMOPT_VERIFY -#include -#include -#endif class TileBbox; @@ -57,27 +53,6 @@ class ClipCache { break; } } -#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; } @@ -100,11 +75,6 @@ class ClipCache { // 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 Level& level = s.levels[bbox.zoom]; if ((s.usedMask & (1u << bbox.zoom)) && !(level.index == bbox.index)) { @@ -159,9 +129,6 @@ class ClipCache { 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) { @@ -178,9 +145,6 @@ class ClipCache { 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; @@ -188,9 +152,6 @@ class ClipCache { 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; } diff --git a/include/geom.h b/include/geom.h index 31902fec..a599da42 100644 --- a/include/geom.h +++ b/include/geom.h @@ -8,12 +8,6 @@ using uint = unsigned int; #include #include -#ifdef TMOPT_VERIFY -#include -#include -#include -#include -#endif // boost::geometry #define BOOST_GEOMETRY_NO_ROBUSTNESS @@ -105,25 +99,7 @@ bool fast_clip(MultiPolygon &mp, Box const &bbox); 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/lua_pool.h b/include/lua_pool.h index 063d58a2..7d536225 100644 --- a/include/lua_pool.h +++ b/include/lua_pool.h @@ -21,10 +21,6 @@ #include #include #include -#ifdef TMOPT_VERIFY -#include -#include -#endif class LuaPool { public: @@ -56,10 +52,6 @@ class LuaPool { 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); @@ -100,36 +92,17 @@ class LuaPool { 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; diff --git a/include/osm_lua_processing.h b/include/osm_lua_processing.h index cfc31659..ca22a5f0 100644 --- a/include/osm_lua_processing.h +++ b/include/osm_lua_processing.h @@ -195,16 +195,8 @@ class OsmLuaProcessing { // 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; - } + static bool PrecheckValid(const GeometryT &) { return false; } + static bool PrecheckValid(const MultiPolygon &geom) { return multiPolygonDefinitelyValid(geom); } // ---- Requests from Lua to write this way/node to a vector tile's Layer template @@ -219,7 +211,7 @@ class OsmLuaProcessing { { geom::validity_failure_type failure = geom::validity_failure_type::no_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); + const bool pre = (isRelation || isWay) && PrecheckValid(geom); if (isRelation && !pre && !geom::is_valid(geom,failure)) { if (verbose) std::cout << "Relation " << originalOsmID << " has " << boost_validity_error(failure) << std::endl; } else if (isWay && !pre && !geom::is_valid(geom,failure)) { @@ -235,7 +227,7 @@ class OsmLuaProcessing { // 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)) + if (PrecheckValid(geom) || geom::is_valid(geom)) return CorrectGeometryResult::Reversed; ReverseRings(geom); } diff --git a/include/polylabel.h b/include/polylabel.h index d56356c1..6ce44f16 100644 --- a/include/polylabel.h +++ b/include/polylabel.h @@ -34,11 +34,6 @@ #include #define POLYLABEL_SSE2 1 #endif -#ifdef TMOPT_VERIFY -#include -#include -#include -#endif namespace mapbox { @@ -108,7 +103,7 @@ auto pointToPolygonDist(const Point& point, const Polygon& polygon) { // 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) { +double pointToPolygonDistKScalar(const Point& point, const Polygon& polygon, int& kr, std::size_t& ki) { bool inside = false; auto minDistSq = std::numeric_limits::infinity(); kr = -1; @@ -149,7 +144,7 @@ double pointToPolygonDistKOrig(const Point& point, const Polygon& polygon, int& #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. +// pointToPolygonDistKScalar 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 { @@ -205,7 +200,7 @@ struct PolygonEdges { // 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. +// pointToPolygonDistKScalar 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 @@ -280,37 +275,16 @@ inline double pointToPolygonDistKSoA(const Point& point, const PolygonEdges& e, } #endif -// pointToPolygonDistKOrig, via the edge arrays when polylabel() has built them for this polygon +// pointToPolygonDistKScalar, 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); + return pointToPolygonDistKScalar(point, polygon, kr, ki); } struct Cell { @@ -365,51 +339,6 @@ inline double edgeDistBound(const Point& p, const Cell& near, const Polygon& pol 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 @@ -481,20 +410,7 @@ Point polylabel(const Polygon& polygon, double precision = 0.00001, bool debug = 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; }; @@ -550,20 +466,6 @@ Point polylabel(const Polygon& polygon, double precision = 0.00001, bool debug = // 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/src/attribute_store.cpp b/src/attribute_store.cpp index c81c5027..6e661c8d 100644 --- a/src/attribute_store.cpp +++ b/src/attribute_store.cpp @@ -216,29 +216,6 @@ uint32_t AttributePairStore::addPair(AttributePair& pair, bool isHot) { 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; } } diff --git a/src/geom.cpp b/src/geom.cpp index 406905e9..50f54705 100644 --- a/src/geom.cpp +++ b/src/geom.cpp @@ -11,11 +11,6 @@ #include #include #include -#ifdef TMOPT_VERIFY -#include -#include -#include -#endif #include #include @@ -431,31 +426,6 @@ bool multiPolygonDefinitelyValid(const MultiPolygon &mp) { 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. @@ -472,26 +442,9 @@ static bool repair_one_polygon(const Polygon &p, double minArea, MultiPolygon &o 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; } @@ -535,9 +488,6 @@ bool repair_multi_polygon(MultiPolygon &mp, bool knownInvalid) { 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; } @@ -559,9 +509,6 @@ bool repair_multi_polygon(MultiPolygon &mp, bool knownInvalid) // 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; diff --git a/src/osm_lua_processing.cpp b/src/osm_lua_processing.cpp index 159a7bc7..c719eaac 100644 --- a/src/osm_lua_processing.cpp +++ b/src/osm_lua_processing.cpp @@ -783,15 +783,6 @@ 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) { @@ -801,12 +792,6 @@ const Linestring &OsmLuaProcessing::linestringCached() { 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; } @@ -896,19 +881,6 @@ void OsmLuaProcessing::Layer(const string &layerName, bool area) { 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, /*allowReverse=*/isWay); if(correctionResult == CorrectGeometryResult::Invalid) return; @@ -956,16 +928,6 @@ void OsmLuaProcessing::Layer(const string &layerName, bool area) { 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) { @@ -1121,13 +1083,6 @@ Point OsmLuaProcessing::calculateCentroid(CentroidAlgorithm algorithm) { // (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); @@ -1384,9 +1339,6 @@ AttrKeyCacheEntry& lookupAttrKey(AttributeKeyStore& keyStore, const char* key, s 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 @@ -1421,9 +1373,6 @@ void OsmLuaProcessing::attributeFast(const char *key, size_t keyLen, T val, char 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); diff --git a/src/osm_mem_tiles.cpp b/src/osm_mem_tiles.cpp index 37f465e6..a20bfca8 100644 --- a/src/osm_mem_tiles.cpp +++ b/src/osm_mem_tiles.cpp @@ -156,25 +156,6 @@ void OsmMemTiles::populateMultiPolygon(MultiPolygon& dst, NodeID objectID) { 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 628696ed..3d5703ca 100644 --- a/src/output_object.cpp +++ b/src/output_object.cpp @@ -215,30 +215,6 @@ void OutputObject::writeAttributes( 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 diff --git a/src/pbf_processor.cpp b/src/pbf_processor.cpp index eac83b41..a602f6e0 100644 --- a/src/pbf_processor.cpp +++ b/src/pbf_processor.cpp @@ -429,17 +429,6 @@ bool PbfProcessor::ReadBlock( 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); diff --git a/src/shp_processor.cpp b/src/shp_processor.cpp index f74e0704..f3b3c4f2 100644 --- a/src/shp_processor.cpp +++ b/src/shp_processor.cpp @@ -367,9 +367,6 @@ void ShpProcessor::processShapeGeometry(SHPObject* shape, AttributeIndex attrIdx // Make valid if needs be string 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); diff --git a/src/sorted_node_store.cpp b/src/sorted_node_store.cpp index 1930697e..c6d282e8 100644 --- a/src/sorted_node_store.cpp +++ b/src/sorted_node_store.cpp @@ -263,50 +263,6 @@ namespace SortedNodeStoreTypes { } } -#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), @@ -424,7 +380,6 @@ LatpLon SortedNodeStore::at(const NodeID id) const { return uncompressedNode(basePtr, id); CachedChunk scratch; decodeChunk((const CompressedChunkInfo*)basePtr, scratch); - TMOPT_VERIFY_CHUNK((const CompressedChunkInfo*)basePtr, scratch); return cachedNode(scratch, id); } @@ -441,7 +396,6 @@ LatpLon SortedNodeStore::at(const NodeID id) const { 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); } @@ -476,7 +430,6 @@ void SortedNodeStore::at(const NodeID* ids, size_t n, LatpLon* out) const { 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]; @@ -484,15 +437,6 @@ void SortedNodeStore::at(const NodeID* ids, size_t n, LatpLon* out) const { out[i] = cachedNode(*entry, id); } -#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(); - } - } -#endif } size_t SortedNodeStore::size() const { diff --git a/src/tile_coordinates_set.cpp b/src/tile_coordinates_set.cpp index 462505aa..e6d6e78e 100644 --- a/src/tile_coordinates_set.cpp +++ b/src/tile_coordinates_set.cpp @@ -4,9 +4,6 @@ #include #include #include -#ifdef TMOPT_VERIFY -#include -#endif #ifdef _MSC_VER #include #endif @@ -106,29 +103,6 @@ void PreciseTileCoordinatesSet::forEachSet(const std::function 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) { diff --git a/src/tile_data.cpp b/src/tile_data.cpp index e34251f8..3fff009a 100644 --- a/src/tile_data.cpp +++ b/src/tile_data.cpp @@ -87,9 +87,6 @@ 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; @@ -105,9 +102,6 @@ void TileDataSource::finalize(size_t threadNum) { // 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; @@ -128,12 +122,6 @@ 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); @@ -410,9 +398,6 @@ Geometry TileDataSource::buildWayGeometry(OutputGeometryType const geomType, geom::correct(mp); geom::validity_failure_type failure = geom::validity_failure_type::no_failure; bool valid = multiPolygonDefinitelyValid(mp); -#ifdef TMOPT_VERIFY - precheckVerify(mp, valid, PRECHECK_TDS); -#endif if (!valid) valid = geom::is_valid(mp,failure); if (!valid) { @@ -535,18 +520,6 @@ 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( @@ -558,11 +531,6 @@ void TileDataSource::addGeometryToIndex( tiles.clear(); try { 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; @@ -595,10 +563,6 @@ void TileDataSource::addGeometryToIndex( // Smaller objects - add to each individual tile index 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); @@ -621,11 +585,6 @@ void TileDataSource::addGeometryToIndex( for (const auto& ls : geom) { 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); @@ -652,18 +611,6 @@ void TileDataSource::addGeometryToIndex( } 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 (const TileCoordinates& index : tiles) { diff --git a/src/tile_worker.cpp b/src/tile_worker.cpp index 614ed30c..4f2f21d4 100644 --- a/src/tile_worker.cpp +++ b/src/tile_worker.cpp @@ -283,28 +283,6 @@ void writeMultiPolygon( 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 (!knownPolygonsValid && (simplifyLevel > 0 || verbose) && !geom::is_valid(current, failure)) { @@ -642,25 +620,6 @@ void outputProc( 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/visvalingam.cpp b/src/visvalingam.cpp index 5e79e051..c278f5b2 100644 --- a/src/visvalingam.cpp +++ b/src/visvalingam.cpp @@ -33,129 +33,7 @@ #include #include #include -#ifdef TMOPT_VERIFY -#include -#include -#endif -#ifdef TMOPT_VERIFY -// The original pointer-based implementation, used to cross-check the indexed one. -namespace visOld { - -struct visItem { - double area = 0; // triangle area - int pointIndex = 0; // index of point in original path - - // to keep a virtual linked list to help rebuild the triangle areas as we remove points. - visItem *next = NULL; - visItem *previous = NULL; - - int index = 0; // interal index in heap, for removal and update -}; - -struct minHeap { - std::vector h; - - void Reserve(size_t size) { - h.reserve(size); - } - - void Push(visItem *item) { - item->index = h.size(); - h.push_back(item); - up(item->index); - } - - visItem *Pop() { - visItem *removed = h[0]; - visItem *lastItem = h[h.size() - 1]; - h.pop_back(); - - if (h.size() > 0) { - lastItem->index = 0; - h[0] = lastItem; - down(0); - } - - return removed; - } - - void Update(visItem *item, double area) { - if (item->area > area) { - // area got smaller - item->area = area; - up(item->index); - } else { - // area got larger - item->area = area; - down(item->index); - } - } - - void up(int i) { - visItem *object = h[i]; - while (i > 0) { - int up = ((i + 1) >> 1) - 1; - visItem *parent = h[up]; - - if (parent->area <= object->area) { - // parent is smaller so we're done fixing up the heap. - break; - } - - // swap nodes - parent->index = i; - h[i] = parent; - - object->index = up; - h[up] = object; - - i = up; - } - } - - void down(int i) { - visItem *object = h[i]; - while (1) { - size_t right = (i + 1) << 1; - size_t left = right - 1; - - int down = i; - visItem *child = h[down]; - - // swap with smallest child - if (left < h.size() && h[left]->area < child->area) { - down = left; - child = h[down]; - } - - if (right < h.size() && h[right]->area < child->area) { - down = right; - child = h[down]; - } - - // non smaller, so quit - if (down == i) { - break; - } - - // swap the nodes - child->index = i; - h[child->index] = child; - - object->index = down; - h[down] = object; - - i = down; - } - } -}; - -} -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! @@ -305,100 +183,6 @@ 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 visvalingamOld(const GeometryType &ls, double threshold, size_t retain) { - size_t start = 0; - size_t end = ls.size(); // - 1; - int removed = 0; - threshold *= 2; - - // build the initial minheap linked list. - minHeap heap; - heap.Reserve(end - start); - - visItem linkedListStart; - linkedListStart.area = INFINITY; - linkedListStart.pointIndex = 0; - heap.Push(&linkedListStart); - - // internal path items - std::vector items; - items.resize((end - start)); - - { - visItem *previous = &linkedListStart; - for (size_t i = 1; i < (end - start) - 1; i++) { - visItem *item = &items[i]; - - item->area = doubleTriangleArea(ls, start, i - 1, i, i + 1); - item->pointIndex = i; - item->previous = previous; - - heap.Push(item); - previous->next = item; - previous = item; - } - - // final item - visItem *endItem = &items[(end - start) - 1]; - endItem->area = INFINITY; - endItem->pointIndex = (end - start) - 1; - endItem->previous = previous; - - previous->next = endItem; - heap.Push(endItem); - } - - // run through the reduction process - while (heap.h.size() > 0) { - visItem *current = heap.Pop(); - if (current->area > threshold) { - break; - } else if ((end - start) - removed <= retain) { - break; - } - - visItem *next = current->next; - visItem *previous = current->previous; - - // remove current element from linked list - previous->next = current->next; - next->previous = current->previous; - removed++; - - // figure out the new areas - if (previous->previous != NULL) { - double area = doubleTriangleArea(ls, start, - previous->previous->pointIndex, - previous->pointIndex, - next->pointIndex); - - area = std::max(area, current->area); - heap.Update(previous, area); - } - - if (next->next != NULL) { - double area = doubleTriangleArea(ls, start, - previous->pointIndex, - next->pointIndex, - next->next->pointIndex); - - area = std::max(area, current->area); - heap.Update(next, area); - } - } - - GeometryType output; - output.reserve(end - start - removed); - visItem *item = &linkedListStart; - while (item != NULL) { - output.emplace_back(ls[item->pointIndex + start]); - item = item->next; - } - return output; -} -#endif template GeometryType visvalingam(const GeometryType &ls, double threshold, size_t retain) { @@ -488,21 +272,6 @@ GeometryType visvalingam(const GeometryType &ls, double threshold, size_t retain 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; } From 80199db49021861e4e2dcbaddd530e0fabb3ffbd Mon Sep 17 00:00:00 2001 From: dktue <2878623+dktue@users.noreply.github.com> Date: Tue, 27 Oct 2026 14:07:51 +0100 Subject: [PATCH 3/3] fix: include for std::size_t in coordinates.h --- include/coordinates.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/coordinates.h b/include/coordinates.h index 9d6f458a..d826cddd 100644 --- a/include/coordinates.h +++ b/include/coordinates.h @@ -7,6 +7,7 @@ // // Things that pull in boost::geometry should go in coordinates_geom.h +#include #include #include #include @@ -130,7 +131,7 @@ 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); +void sortUniqueTiles(std::vector &tiles, std::size_t start = 0); #endif //_COORDINATES_H