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; +} +