From d87e746c65963f1bd1e01836de4d269d90df25ea Mon Sep 17 00:00:00 2001 From: Colin Dellow Date: Wed, 3 Jan 2024 18:26:58 -0500 Subject: [PATCH] translate mapbox polylabel to boost This is just the translation of the algorithm and some smoke testing that it seems to do the right thing. It's not integrated into Tilemaker yet, that's blocked on https://github.com/systemed/tilemaker/pull/632 --- Makefile | 5 + include/polylabel.h | 230 ++++++++++++++++++++++++++++++++++++++++ test/polylabel.test.cpp | 44 ++++++++ 3 files changed, 279 insertions(+) create mode 100644 include/polylabel.h create mode 100644 test/polylabel.test.cpp diff --git a/Makefile b/Makefile index 04eb3c42..16b4c12a 100644 --- a/Makefile +++ b/Makefile @@ -164,6 +164,11 @@ test_options_parser: \ test/options_parser.test.o $(CXX) $(CXXFLAGS) -o test.options_parser $^ $(INC) $(LIB) $(LDFLAGS) && ./test.options_parser +test_polylabel: \ + test/polylabel.test.o \ + src/geom.o + $(CXX) $(CXXFLAGS) -o test.polylabel $^ $(INC) $(LIB) $(LDFLAGS) && ./test.polylabel + test_pooled_string: \ src/mmap_allocator.o \ src/pooled_string.o \ diff --git a/include/polylabel.h b/include/polylabel.h new file mode 100644 index 00000000..96fc7ac4 --- /dev/null +++ b/include/polylabel.h @@ -0,0 +1,230 @@ +#pragma once + +// Original source from https://github.com/mapbox/polylabel, licensed +// under ISC. +// +// Adapted to use Boost Geometry instead of MapBox's geometry library. +// + +#include "geom.h" + +#include +#include +#include +#include + +namespace mapbox { + +namespace detail { + +// get squared distance from a point to a segment +double getSegDistSq(const Point& p, + const Point& a, + const Point& b) { + //auto x = a.x; + auto x = a.get<0>(); + //auto y = a.y; + auto y = a.get<1>(); + //auto dx = b.x - x; + auto dx = b.get<0>() - x; + //auto dy = b.y - y; + auto dy = b.get<1>() - y; + + if (dx != 0 || dy != 0) { + + //auto t = ((p.x - x) * dx + (p.y - y) * dy) / (dx * dx + dy * dy); + auto t = ((p.get<0>() - x) * dx + (p.get<1>() - y) * dy) / (dx * dx + dy * dy); + + if (t > 1) { + //x = b.x; + //y = b.y; + x = b.get<0>(); + y = b.get<1>(); + + } else if (t > 0) { + x += dx * t; + y += dy * t; + } + } + + //dx = p.x - x; + //dy = p.y - y; + dx = p.get<0>() - x; + dy = p.get<1>() - y; + + return dx * dx + dy * dy; +} + +// signed distance from point to polygon outline (negative if point is outside) +auto pointToPolygonDist(const Point& point, const Polygon& polygon) { + bool inside = false; + auto minDistSq = std::numeric_limits::infinity(); + + { + 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.y > point.y) != (b.y > point.y) && +// (point.x < (b.x - a.x) * (point.y - a.y) / (b.y - a.y) + a.x)) inside = !inside; + 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; + + minDistSq = std::min(minDistSq, getSegDistSq(point, a, b)); + } + } + + 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.y > point.y) != (b.y > point.y) && +// (point.x < (b.x - a.x) * (point.y - a.y) / (b.y - a.y) + a.x)) inside = !inside; + 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; + + minDistSq = std::min(minDistSq, getSegDistSq(point, a, b)); + } + } + + return (inside ? 1 : -1) * std::sqrt(minDistSq); +} + +struct Cell { + Cell(const Point& c_, double h_, const Polygon& polygon) + : c(c_), + h(h_), + d(pointToPolygonDist(c, polygon)), + max(d + h * std::sqrt(2)) + {} + + Point c; // cell center + double h; // half the cell size + double d; // distance from cell center to polygon + double max; // max distance to polygon within a cell +}; + +// get polygon centroid +Cell getCentroidCell(const Polygon& polygon) { + double area = 0; + //Point c { 0, 0 }; + double cx = 0, cy = 0; + //const auto& ring = polygon.at(0); + const auto& ring = polygon.outer(); + + for (std::size_t i = 0, len = ring.size(), j = len - 1; i < len; j = i++) { + const Point& a = ring[i]; + const Point& b = ring[j]; + auto f = a.get<0>() * b.get<1>() - b.get<0>() * a.get<1>(); + cx += (a.get<0>() + b.get<0>()) * f; + cy += (a.get<1>() + b.get<1>()) * f; + area += f * 3; + } + + Point c { cx, cy }; + // return Cell(area == 0 ? ring.at(0) : c / area, 0, polygon); + // TODO: mapbox's library seems to define division for points? + // https://github.com/mapbox/geometry.hpp/blob/master/include/mapbox/geometry/point.hpp + // doesn't have a division operator, though... + // ah, it's in https://github.com/mapbox/geometry.hpp/blob/12ac5412bf85571852ad1cd7c30456faef8d6464/include/mapbox/geometry/point_arithmetic.hpp#L48-L52 + return Cell(area == 0 ? ring.at(0) : Point { cx / area, cy / area }, 0, polygon); +} + +} // namespace detail + +Point polylabel(const Polygon& polygon, double precision = 1, bool debug = false) { + using namespace detail; + + // find the bounding box of the outer ring + //const geometry::box envelope = geometry::envelope(polygon.at(0)); + Box envelope; + geom::envelope(polygon.outer(), envelope); + +// const Point size { +// envelope.max.x - envelope.min.x, +// envelope.max.y - envelope.min.y +// }; + + 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.x, size.y); + const double cellSize = std::min(size.get<0>(), size.get<1>()); + double h = cellSize / 2; + + // a priority queue of cells in order of their "potential" (max distance to polygon) + auto compareMax = [] (const Cell& a, const Cell& b) { + return a.max < b.max; + }; + using Queue = std::priority_queue, decltype(compareMax)>; + Queue cellQueue(compareMax); + + if (cellSize == 0) { + //return envelope.min; + return envelope.min_corner(); + } + + // cover polygon with initial cells +// for (double x = envelope.min.x; x < envelope.max.x; x += cellSize) { +// for (double y = envelope.min.y; y < envelope.max.y; y += cellSize) { +// cellQueue.push(Cell({x + h, y + h}, h, polygon)); +// } +// } + + 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)); + } + } + + // take centroid as the first best guess + auto bestCell = getCentroidCell(polygon); + + // second guess: bounding box centroid + //Cell bboxCell(envelope.min + size / 2.0, 0, 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; + } + + auto numProbes = cellQueue.size(); + while (!cellQueue.empty()) { + // pick the most promising cell from the queue + auto cell = cellQueue.top(); + cellQueue.pop(); + + // update the best cell if we found a better one + if (cell.d > bestCell.d) { + bestCell = cell; + if (debug) std::cout << "found best " << ::round(1e4 * cell.d) / 1e4 << " after " << numProbes << " probes" << std::endl; + } + + // do not drill down further if there's no chance of a better solution + if (cell.max - bestCell.d <= precision) continue; + + // split the cell into four cells + 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)); + numProbes += 4; + } + + if (debug) { + std::cout << "num probes: " << numProbes << std::endl; + std::cout << "best distance: " << bestCell.d << std::endl; + } + + return bestCell.c; +} + +} // namespace mapbox diff --git a/test/polylabel.test.cpp b/test/polylabel.test.cpp new file mode 100644 index 00000000..75a39ffe --- /dev/null +++ b/test/polylabel.test.cpp @@ -0,0 +1,44 @@ +#include +#include "geom.h" +#include "external/minunit.h" +#include "geojson.h" +#include "./simplify.test.h" +#include "polylabel.h" + +void save(std::string filename, const MultiPolygon& mp) { + GeoJSON json; + json.addGeometry(mp); + json.finalise(); + json.toFile(filename); +} + +MU_TEST(test_simplify) { + //MultiPolygon mp0 = mthope(); + MultiPolygon mp0 = gsmnp(); + for (const auto& p : mp0) { + std::cout << "mp0: poly.outer().size() = " << p.outer().size() << std::endl; + } + save("poly-s0.txt", mp0); + + timespec start, end; + clock_gettime(CLOCK_MONOTONIC, &start); + auto pt = mapbox::polylabel(mp0.at(0)); + std::cout << " point is at " << pt.get<0>() << ", " << pt.get<1>() << std::endl; + + clock_gettime(CLOCK_MONOTONIC, &end); + uint64_t elapsedNs = 1e9 * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec; + std::cout << "took " << std::to_string((uint32_t)(elapsedNs / 1e6)) << " ms" << std::endl; + + +} + +MU_TEST_SUITE(test_suite_simplify) { + MU_RUN_TEST(test_simplify); +} + +int main() { + MU_RUN_SUITE(test_suite_simplify); + MU_REPORT(); + return MU_EXIT_CODE; +} +