From 4ceddbcfb9790cf1cbd97f916d990f4dce15a0b6 Mon Sep 17 00:00:00 2001 From: Richard Fairhurst Date: Wed, 16 Oct 2024 09:03:08 +0100 Subject: [PATCH 1/2] Visvalingam-Whyatt simplification --- CMakeLists.txt | 1 + Makefile | 1 + README.md | 1 + docs/CONFIGURATION.md | 1 + include/shared_data.h | 6 +- include/visvalingam.h | 10 ++ src/shared_data.cpp | 8 +- src/tile_worker.cpp | 23 +++- src/visvalingam.cpp | 265 ++++++++++++++++++++++++++++++++++++++++++ 9 files changed, 307 insertions(+), 9 deletions(-) create mode 100644 include/visvalingam.h create mode 100644 src/visvalingam.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 74694e04..c5c2309d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -113,6 +113,7 @@ file(GLOB tilemaker_src_files src/tile_data.cpp src/tilemaker.cpp src/tile_worker.cpp + src/visvalingam.cpp src/way_stores.cpp ) add_executable(tilemaker ${tilemaker_src_files}) diff --git a/Makefile b/Makefile index 8f6f86de..ccbd1b7e 100644 --- a/Makefile +++ b/Makefile @@ -131,6 +131,7 @@ tilemaker: \ src/tile_data.o \ src/tilemaker.o \ src/tile_worker.o \ + src/visvalingam.o \ src/way_stores.o $(CXX) $(CXXFLAGS) -o tilemaker $^ $(INC) $(LIB) $(LDFLAGS) diff --git a/README.md b/README.md index f264d84e..caa1a98c 100644 --- a/README.md +++ b/README.md @@ -126,4 +126,5 @@ Licenses of third-party libraries: - [Simple-Web-Server](https://gitlab.com/eidheim/Simple-Web-Server) is licensed under MIT - [sqlite_modern_cpp](https://github.com/SqliteModernCpp/sqlite_modern_cpp) is licensed under MIT - [streamvbyte](https://github.com/lemire/streamvbyte) is licensed under Apache 2 +- [visvalingam.cpp](https://github.com/felt/tippecanoe/blob/main/visvalingam.cpp) is licensed under MIT - [vtzero](https://github.com/mapbox/vtzero) is licensed under BSD 2-clause diff --git a/docs/CONFIGURATION.md b/docs/CONFIGURATION.md index aae02478..5919f317 100644 --- a/docs/CONFIGURATION.md +++ b/docs/CONFIGURATION.md @@ -77,6 +77,7 @@ You can add optional parameters to layers: * `simplify_level` - how much to simplify features (in degrees of longitude) on the zoom level `simplify_below-1` * `simplify_length` - how much to simplify features (in kilometers) on the zoom level `simplify_below-1`, preceding `simplify_level` * `simplify_ratio` - (optional: the default value is 2.0) the actual simplify level will be `simplify_level * pow(simplify_ratio, (simplify_below-1) - )` +* `simplify_algorithm` - which simplification algorithm to use (defaults to Douglas-Peucker; you can also specify `"visvalingam"`, which can be better for landuse and similar polygons) * `filter_below` - filter areas by minimum size below this zoom level * `filter_area` - minimum size (in square degrees of longitude) for the zoom level `filter_below-1` * `feature_limit` - restrict the number of features written to each tile diff --git a/include/shared_data.h b/include/shared_data.h index d51cebba..7ddb3292 100644 --- a/include/shared_data.h +++ b/include/shared_data.h @@ -23,6 +23,7 @@ struct LayerDef { double simplifyLevel; double simplifyLength; double simplifyRatio; + uint simplifyAlgo; uint filterBelow; double filterArea; uint combinePolygonsBelow; @@ -41,6 +42,9 @@ struct LayerDef { const bool useColumn(std::string &col) { return allSourceColumns || (std::find(sourceColumns.begin(), sourceColumns.end(), col) != sourceColumns.end()); } + + static const uint DOUGLAS_PEUCKER = 0; + static const uint VISVALINGAM = 1; }; ///\brief Defines layers used in map rendering @@ -53,7 +57,7 @@ class LayerDefinition { // Define a layer (as read from the .json file) uint addLayer(std::string name, uint minzoom, uint maxzoom, - uint simplifyBelow, double simplifyLevel, double simplifyLength, double simplifyRatio, + uint simplifyBelow, double simplifyLevel, double simplifyLength, double simplifyRatio, uint simplifyAlgo, uint filterBelow, double filterArea, uint combinePolygonsBelow, bool sortZOrderAscending, uint featureLimit, uint featureLimitBelow, bool combinePoints, const std::string &source, diff --git a/include/visvalingam.h b/include/visvalingam.h new file mode 100644 index 00000000..7bc29da8 --- /dev/null +++ b/include/visvalingam.h @@ -0,0 +1,10 @@ +/*! \file */ +#ifndef _VISVALINGAM_H +#define _VISVALINGAM_H + +// Visvalingam simplify +Linestring simplifyVis(const Linestring &ls, double max_distance); +Polygon simplifyVis(const Polygon &p, double max_distance); +MultiPolygon simplifyVis(const MultiPolygon &mp, double max_distance); + +#endif //_VISVALINGAM_H diff --git a/src/shared_data.cpp b/src/shared_data.cpp index 1692b187..7432e769 100644 --- a/src/shared_data.cpp +++ b/src/shared_data.cpp @@ -135,7 +135,7 @@ void SharedData::writePMTilesBounds() { // Define a layer (as read from the .json file) uint LayerDefinition::addLayer(string name, uint minzoom, uint maxzoom, - uint simplifyBelow, double simplifyLevel, double simplifyLength, double simplifyRatio, + uint simplifyBelow, double simplifyLevel, double simplifyLength, double simplifyRatio, uint simplifyAlgo, uint filterBelow, double filterArea, uint combinePolygonsBelow, bool sortZOrderAscending, uint featureLimit, uint featureLimitBelow, bool combinePoints, const std::string &source, @@ -146,7 +146,7 @@ uint LayerDefinition::addLayer(string name, uint minzoom, uint maxzoom, const std::string &writeTo) { bool isWriteTo = !writeTo.empty(); - LayerDef layer = { name, minzoom, maxzoom, simplifyBelow, simplifyLevel, simplifyLength, simplifyRatio, + LayerDef layer = { name, minzoom, maxzoom, simplifyBelow, simplifyLevel, simplifyLength, simplifyRatio, simplifyAlgo, filterBelow, filterArea, combinePolygonsBelow, sortZOrderAscending, featureLimit, featureLimitBelow, combinePoints, source, sourceColumns, allSourceColumns, indexed, indexName, std::map(), isWriteTo }; @@ -319,6 +319,8 @@ void Config::readConfig(rapidjson::Document &jsonConfig, bool &hasClippingBox, B int featureLimitBelow= it->value.HasMember("feature_limit_below") ? it->value["feature_limit_below"].GetInt() : (maxZoom+1); bool combinePoints = it->value.HasMember("combine_points" ) ? it->value["combine_points" ].GetBool() : true; bool sortZOrderAscending = it->value.HasMember("z_order_ascending") ? it->value["z_order_ascending"].GetBool() : (featureLimit==0); + string algo = it->value.HasMember("simplify_algorithm") ? it->value["simplify_algorithm"].GetString() : ""; + uint simplifyAlgo = algo=="visvalingam" ? LayerDef::VISVALINGAM : LayerDef::DOUGLAS_PEUCKER; string source = it->value.HasMember("source") ? it->value["source"].GetString() : ""; vector sourceColumns; bool allSourceColumns = false; @@ -337,7 +339,7 @@ void Config::readConfig(rapidjson::Document &jsonConfig, bool &hasClippingBox, B string indexName = it->value.HasMember("index_column") ? it->value["index_column"].GetString() : ""; layers.addLayer(layerName, minZoom, maxZoom, - simplifyBelow, simplifyLevel, simplifyLength, simplifyRatio, + simplifyBelow, simplifyLevel, simplifyLength, simplifyRatio, simplifyAlgo, filterBelow, filterArea, combinePolyBelow, sortZOrderAscending, featureLimit, featureLimitBelow, combinePoints, source, sourceColumns, allSourceColumns, indexed, indexName, writeTo); diff --git a/src/tile_worker.cpp b/src/tile_worker.cpp index 958ca923..bd0eb877 100644 --- a/src/tile_worker.cpp +++ b/src/tile_worker.cpp @@ -5,6 +5,7 @@ #include #include #include "helpers.h" +#include "visvalingam.h" using namespace std; extern bool verbose; @@ -100,6 +101,7 @@ void writeMultiLinestring( const OutputObjectID& oo, unsigned zoom, double simplifyLevel, + unsigned simplifyAlgo, const MultiLinestring& mls ) { vtzero::linestring_feature_builder fbuilder{vtLayer}; @@ -114,7 +116,11 @@ void writeMultiLinestring( if (simplifyLevel>0) { for(auto const &ls: mls) { - tmp.push_back(simplify(ls, simplifyLevel)); + if (simplifyAlgo==LayerDef::VISVALINGAM) { + tmp.push_back(simplifyVis(ls, simplifyLevel)); + } else { + tmp.push_back(simplify(ls, simplifyLevel)); + } } toWrite = &tmp; } else { @@ -208,11 +214,16 @@ void writeMultiPolygon( const OutputObjectID& oo, unsigned zoom, double simplifyLevel, + unsigned simplifyAlgo, const MultiPolygon& mp ) { MultiPolygon current = bbox.scaleGeometry(mp); if (simplifyLevel>0) { - current = simplify(current, simplifyLevel/bbox.xscale); + if (simplifyAlgo == LayerDef::VISVALINGAM) { + current = simplifyVis(current, simplifyLevel/bbox.xscale); + } else { + current = simplify(current, simplifyLevel/bbox.xscale); + } geom::remove_spikes(current); } if (geom::is_empty(current)) @@ -264,6 +275,7 @@ void ProcessObjects( OutputObjectsConstIt ooSameLayerEnd, class SharedData& sharedData, double simplifyLevel, + unsigned simplifyAlgo, double filterArea, bool combinePolygons, bool combinePoints, @@ -350,9 +362,9 @@ void ProcessObjects( } if (oo.oo.geomType == LINESTRING_ || oo.oo.geomType == MULTILINESTRING_) - writeMultiLinestring(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, boost::get(g)); + writeMultiLinestring(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, simplifyAlgo, boost::get(g)); else if (oo.oo.geomType == POLYGON_) - writeMultiPolygon(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, boost::get(g)); + writeMultiPolygon(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, simplifyAlgo, boost::get(g)); } } } @@ -436,7 +448,8 @@ void ProcessLayer( if (ld.featureLimit>0 && end-ooListSameLayer.first>ld.featureLimit && zoom3) { diff --git a/src/visvalingam.cpp b/src/visvalingam.cpp new file mode 100644 index 00000000..7ae43e8f --- /dev/null +++ b/src/visvalingam.cpp @@ -0,0 +1,265 @@ +#include "geom.h" +#include + +// Adapted from +// https://github.com/paulmach/orb/blob/dcade4901baea0727377ccf7c4aab2addd92d152/simplify/visvalingam.go +// and from +// https://github.com/felt/tippecanoe/blob/main/visvalingam.cpp + +// The MIT License (MIT) +// +// Copyright (c) 2017 Paul Mach +// +// Permission is hereby granted, free of charge, to any person obtaining +// a copy of this software and associated documentation files (the +// "Software"), to deal in the Software without restriction, including +// without limitation the rights to use, copy, modify, merge, publish, +// distribute, sublicense, and/or sell copies of the Software, and to +// permit persons to whom the Software is furnished to do so, subject +// to the following conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR +// ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF +// CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION +// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +#include +#include + +// Stuff to create the priority queue, or min heap. +// Rewriting it here, vs using the std lib, resulted in a 50% performance bump! + +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 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; + } + } +}; + +template +static double doubleTriangleArea(GeometryType const &ls, int start, int i1, int i2, int i3) { + Point a = ls[i1 + start]; + Point b = ls[i2 + start]; + Point c = ls[i3 + start]; + + return std::abs((b.x() - a.x()) * (c.y() - a.y()) - (b.y() - a.y()) * (c.x() - a.x())); +} + +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; + + // build the initial minheap linked list. + minHeap heap; + + 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; + visItem *item = &linkedListStart; + while (item != NULL) { + output.emplace_back(ls[item->pointIndex + start]); + item = item->next; + } + return output; +} + +Linestring simplifyVis(const Linestring &ls, double max_distance) { + if (ls.size()<3) return ls; + bool isClosed = ls[0].x()==ls[ls.size()-1].x() && ls[0].y()==ls[ls.size()-1].y(); + double threshold = max_distance * max_distance * 4; + return visvalingam(ls, threshold, isClosed ? 3 : 2); +} +Polygon simplifyVis(const Polygon &p, double max_distance) { + Polygon output; + double threshold = max_distance * max_distance * 4; + output.outer() = visvalingam(p.outer(), threshold, 4); + for (const auto &ring : p.inners()) { + output.inners().emplace_back(visvalingam(ring, threshold, 4)); + } + return output; +} +MultiPolygon simplifyVis(const MultiPolygon &mp, double max_distance) { + MultiPolygon output; + for (const auto &p : mp) { + output.emplace_back(simplifyVis(p, max_distance)); + } + make_valid(output); + return output; +} From 5e5210f9d78bffc642fa26fbb67eec0f91926728 Mon Sep 17 00:00:00 2001 From: systemed Date: Fri, 18 Oct 2024 09:15:02 +0100 Subject: [PATCH 2/2] Use visvalingam for ocean and admin boundaries --- resources/config-openmaptiles.json | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/resources/config-openmaptiles.json b/resources/config-openmaptiles.json index 01d23bfd..111e2458 100644 --- a/resources/config-openmaptiles.json +++ b/resources/config-openmaptiles.json @@ -1,7 +1,7 @@ { "layers": { "place": { "minzoom": 0, "maxzoom": 14 }, - "boundary": { "minzoom": 0, "maxzoom": 14, "simplify_below": 12, "simplify_level": 0.0003, "simplify_ratio": 2 }, + "boundary": { "minzoom": 0, "maxzoom": 14, "simplify_below": 12, "simplify_level": 0.0003, "simplify_ratio": 2, "simplify_algorithm": "visvalingam" }, "poi": { "minzoom": 12, "maxzoom": 14 }, "poi_detail": { "minzoom": 14, "maxzoom": 14, "write_to": "poi"}, @@ -17,7 +17,7 @@ "building": { "minzoom": 13, "maxzoom": 14 }, "water": { "minzoom": 6, "maxzoom": 14, "simplify_below": 12, "simplify_level": 0.0003, "simplify_ratio": 2}, - "ocean": { "minzoom": 0, "maxzoom": 14, "source": "coastline/water_polygons.shp", "filter_below": 12, "filter_area": 0.5, "simplify_below": 13, "simplify_level": 0.0001, "simplify_ratio": 2, "write_to": "water" }, + "ocean": { "minzoom": 0, "maxzoom": 14, "source": "coastline/water_polygons.shp", "filter_below": 12, "filter_area": 0.5, "simplify_below": 13, "simplify_level": 0.0001, "simplify_ratio": 2, "simplify_algorithm": "visvalingam", "write_to": "water" }, "water_name": { "minzoom": 14, "maxzoom": 14 }, "water_name_detail": { "minzoom": 14, "maxzoom": 14, "write_to": "water_name" },