add polygon center utility, based on Mapbox's polylabel

Tests/CMakeLists.txt

@@ -166,3 +166,6 @@ endif()
 
 #---- SunriseSunsetTest
 osmscout_test_project(NAME SunriseSunsetTest SOURCES src/SunriseSunsetTest.cpp)
+
+#---- PolygonCenterTest
+osmscout_test_project(NAME PolygonCenterTest SOURCES src/PolygonCenter.cpp)

Tests/meson.build

@@ -215,6 +215,13 @@ SunriseSunset = executable('SunriseSunset',
              link_with: [osmscout],
              install: false)
 
+PolygonCenter = executable('PolygonCenter',
+             'src/PolygonCenter.cpp',
+             include_directories: [testIncDir, osmscoutIncDir],
+             dependencies: [mathDep, threadDep, openmpDep],
+             link_with: [osmscout],
+             install: false)
+
 TilingTest = executable('TilingTest',
              'src/TilingTest.cpp',
              include_directories: [testIncDir, osmscoutIncDir],
@@ -302,6 +309,7 @@ else
 endif
 
 test('Check correctness of NumberSet class', NumberSet)
+test('Check PolygonCenter utility', PolygonCenter)
 test('Check scan conversion code', ScanConversion)
 test('Check string utils', StringUtils)
 test('Check tiling calculation code', TilingTest)

libosmscout/CMakeLists.txt

@@ -90,6 +90,7 @@ set(HEADER_FILES_UTIL
     include/osmscout/util/NumberSet.h
     include/osmscout/util/ObjectPool.h
     include/osmscout/util/Parsing.h
+    include/osmscout/util/PolygonCenter.h
     include/osmscout/util/ProcessingQueue.h
     include/osmscout/util/Progress.h
     include/osmscout/util/Projection.h
@@ -164,6 +165,7 @@ set(SOURCE_FILES
     src/osmscout/util/Number.cpp
     src/osmscout/util/NumberSet.cpp
     src/osmscout/util/Parsing.cpp
+    src/osmscout/util/PolygonCenter.cpp
     src/osmscout/util/Progress.cpp
     src/osmscout/util/Projection.cpp
     src/osmscout/util/StopClock.cpp

libosmscout/include/meson.build

@@ -33,6 +33,7 @@ osmscoutHeader = [
             'osmscout/util/NumberSet.h',
             'osmscout/util/ObjectPool.h',
             'osmscout/util/Parsing.h',
+            'osmscout/util/PolygonCenter.h',
             'osmscout/util/ProcessingQueue.h',
             'osmscout/util/Progress.h',
             'osmscout/util/Projection.h',

libosmscout/include/osmscout/util/PolygonCenter.h

@@ -0,0 +1,47 @@
+#ifndef OSMSCOUT_POLYGONCENTER_H
+#define OSMSCOUT_POLYGONCENTER_H
+
+/*
+  This source is part of the libosmscout library
+  Copyright (C) 2022 Lukas Karas
+  Copyright (C) 2016 Mapbox
+
+  ISC License (compatible with GNU LGPL)
+
+  Permission to use, copy, modify, and/or distribute this software for any purpose
+  with or without fee is hereby granted, provided that the above copyright notice
+  and this permission notice appear in all copies.
+
+  THE SOFTWARE IS PROVIDED "AS IS" AND ISC DISCLAIMS ALL WARRANTIES WITH REGARD TO
+  THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+  IN NO EVENT SHALL ISC BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR
+  CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA
+  OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
+  ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS
+  SOFTWARE.
+ */
+
+#include <osmscout/GeoCoord.h>
+#include <osmscout/Point.h>
+#include <osmscout/Area.h>
+#include <osmscout/CoreImportExport.h>
+
+namespace osmscout {
+
+/**
+ * \ingroup Util
+ *
+ * A fast algorithm for finding polygon pole of inaccessibility, the most distant internal point from the polygon
+ * outline (not to be confused with centroid).
+ * Useful for optimal placement of a text label on a polygon.
+ *
+ * Based on PolyLabel algorithm from Mapbox
+ *  - https://blog.mapbox.com/a-new-algorithm-for-finding-a-visual-center-of-a-polygon-7c77e6492fbc
+ *  - https://github.com/mapbox/polylabel
+ */
+extern OSMSCOUT_API GeoCoord PolygonCenter(const Area& area, double precision = 1);
+
+extern OSMSCOUT_API GeoCoord PolygonCenter(const std::vector<Point>& polygon, double precision = 1);
+
+}
+#endif //OSMSCOUT_POLYGONCENTER_H

libosmscout/src/meson.build

@@ -24,6 +24,7 @@ osmscoutSrc = [
             'src/osmscout/util/Number.cpp',
             'src/osmscout/util/NumberSet.cpp',
             'src/osmscout/util/Parsing.cpp',
+            'src/osmscout/util/PolygonCenter.cpp',
             'src/osmscout/util/Progress.cpp',
             'src/osmscout/util/Projection.cpp',
             'src/osmscout/util/StopClock.cpp',

libosmscout/src/osmscout/util/PolygonCenter.cpp

@@ -0,0 +1,206 @@
+/*
+  This source is part of the libosmscout library
+  Copyright (C) 2022 Lukas Karas
+  Copyright (C) 2016 Mapbox
+
+  ISC License (compatible with GNU LGPL)
+
+  Permission to use, copy, modify, and/or distribute this software for any purpose
+  with or without fee is hereby granted, provided that the above copyright notice
+  and this permission notice appear in all copies.
+
+  THE SOFTWARE IS PROVIDED "AS IS" AND ISC DISCLAIMS ALL WARRANTIES WITH REGARD TO
+  THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+  IN NO EVENT SHALL ISC BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR
+  CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA
+  OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
+  ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS
+  SOFTWARE.
+ */
+
+#include <osmscout/util/PolygonCenter.h>
+#include <osmscout/util/Geometry.h>
+#include <osmscout/util/Logger.h>
+
+#include <queue>
+
+namespace osmscout {
+
+namespace {
+// get squared distance from a point to a segment
+double getSegDistSq(const GeoCoord& p,
+                    const Point& a,
+                    const Point& b) {
+  auto x = a.GetLon();
+  auto y = a.GetLat();
+  auto dx = b.GetLon() - x;
+  auto dy = b.GetLat() - y;
+
+  if (dx != 0 || dy != 0) {
+
+    auto t = ((p.GetLon() - x) * dx + (p.GetLat() - y) * dy) / (dx * dx + dy * dy);
+
+    if (t > 1) {
+      x = b.GetLon();
+      y = b.GetLat();
+    } else if (t > 0) {
+      x += dx * t;
+      y += dy * t;
+    }
+  }
+
+  dx = p.GetLon() - x;
+  dy = p.GetLat() - y;
+
+  return dx * dx + dy * dy;
+}
+
+using Ring = std::vector<Point>;
+using Polygon = std::vector<const Ring*>;
+
+// signed distance from point to polygon outline (negative if point is outside)
+auto pointToPolygonDist(const GeoCoord& point, const Polygon& polygon) {
+  bool inside = false;
+  auto minDistSq = std::numeric_limits<double>::infinity();
+
+  for (const auto& ring : polygon) {
+    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.GetLat() > point.GetLat()) != (b.GetLat() > point.GetLat()) &&
+          (point.GetLon() < (b.GetLon() - a.GetLon()) * (point.GetLat() - a.GetLat()) / (b.GetLat() - a.GetLat()) + a.GetLon())) {
+        inside = !inside;
+      }
+
+      minDistSq = std::min(minDistSq, getSegDistSq(point, a, b));
+    }
+  }
+
+  return (inside ? 1 : -1) * std::sqrt(minDistSq);
+}
+
+struct Cell {
+  Cell(const GeoCoord& c_, double h_, const Polygon& polygon)
+      : c(c_),
+        h(h_),
+        d(pointToPolygonDist(c, polygon)),
+        max(d + h * std::sqrt(2))
+  {}
+
+  GeoCoord 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 (of first top-level outer ring)
+Cell getCentroidCell(const Polygon& polygon) {
+  assert(!polygon.empty());
+  double area = 0;
+  GeoCoord c { 0, 0 };
+  const auto &ring = (*polygon[0]);
+
+  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.GetLon() * b.GetLat() - b.GetLon() * a.GetLat();
+    c.Set(
+    c.GetLat() + (a.GetLat() + b.GetLat()) * f,
+    c.GetLon() + (a.GetLon() + b.GetLon()) * f);
+    area += f * 3;
+  }
+
+  assert(!ring.empty());
+  return Cell(area == 0 ? ring[0].GetCoord() : GeoCoord(c.GetLat() / area, c.GetLon() / area),
+              0,
+              polygon);
+}
+
+}
+
+GeoCoord PolygonCenter(const Polygon &polygon, const GeoBox &envelope, double precision)
+{
+  const GeoCoord size{envelope.GetHeight(), envelope.GetWidth()};
+
+  const double cellSize = std::min(size.GetLat(), size.GetLon());
+  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<Cell, std::vector<Cell>, decltype(compareMax)>;
+  Queue cellQueue(compareMax);
+
+  if (cellSize == 0) {
+    return envelope.GetMinCoord();
+  }
+
+  // cover polygon with initial cells
+  for (double lon = envelope.GetMinLon(); lon < envelope.GetMaxLon(); lon += cellSize) {
+    for (double lat = envelope.GetMinLat(); lat < envelope.GetMaxLat(); lat += cellSize) {
+      cellQueue.push(Cell({lat + h, lon + h}, h, polygon));
+    }
+  }
+
+  // take centroid as the first best guess
+  auto bestCell = getCentroidCell(polygon);
+
+  // second guess: bounding box centroid
+  Cell bboxCell(envelope.GetCenter(), 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;
+      log.Debug() << "found best " << ::round(1e4 * cell.d) / 1e4 << " after " << numProbes << " probes";
+    }
+
+    // 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.GetLat() - h, cell.c.GetLon() - h}, h, polygon));
+    cellQueue.push(Cell({cell.c.GetLat() - h, cell.c.GetLon() + h}, h, polygon));
+    cellQueue.push(Cell({cell.c.GetLat() + h, cell.c.GetLon() - h}, h, polygon));
+    cellQueue.push(Cell({cell.c.GetLat() + h, cell.c.GetLon() + h}, h, polygon));
+    numProbes += 4;
+  }
+
+  log.Debug() << "num probes: " << numProbes;
+  log.Debug() << "best distance: " << bestCell.d;
+
+  return bestCell.c;
+}
+
+GeoCoord PolygonCenter(const Area &area, double precision)
+{
+  Polygon polygon;
+  area.VisitRings([&](size_t, const Area::Ring &ring, const TypeInfoRef&) -> bool {
+    polygon.push_back(&ring.nodes);
+    return false; // visit just top-level outer rings
+  });
+  return PolygonCenter(polygon, area.GetBoundingBox(), precision);
+}
+
+GeoCoord PolygonCenter(const std::vector<Point>& ring, double precision)
+{
+  Polygon polygon;
+  polygon.push_back(&ring);
+  GeoBox bbox;
+  GetBoundingBox(ring, bbox);
+  return PolygonCenter(polygon, bbox, precision);
+}
+}