Add pole_of_inaccessibility() function in Lua for polygons

Implementation of the "polylabel" algorithm for finding the "pole of
inaccessibility" or center of the maximum inscribed circle of a polygon.
This is accessible from Lua as the pole_of_inaccessibility() function.

The current implementation is only defined for polygons, because it is
unclear how it should behave for other types of geometries. We can
extend this later if needed.

Author: joto

flex-config/labelpoint.lua

@@ -0,0 +1,63 @@
+-- This config example file is released into the Public Domain.
+--
+-- This example shows the use of the centroid() and pole_of_inaccessibility()
+-- functions. For all named polygons several types of "center" points are
+-- calculated that can be used for labelling.
+
+local tables = {}
+
+tables.polygons = osm2pgsql.define_area_table('polygons', {
+    { column = 'name', type = 'text' },
+    { column = 'tags', type = 'jsonb' },
+    { column = 'geom', type = 'polygon', not_null = true },
+    { column = 'centroid', type = 'point', not_null = true },
+    { column = 'poi1', type = 'point', not_null = true },
+    { column = 'poi2', type = 'point', not_null = true },
+})
+
+function add(tags, geom)
+    tables.polygons:insert({
+        name = tags.name,
+        tags = tags,
+        geom = geom,
+        -- The centroid is easy and fast to calculate
+        centroid = geom:centroid(),
+        -- The pole of inaccessibility can take a bit longer. It calculates
+        -- the place in the polygon that's farthest away from the boundary.
+        -- In other words it finds the center of the largest circle that fits
+        -- completely into the polygon. (Note that due to the time required
+        -- to calculate this, only an approximation will be calculated.)
+        poi1 = geom:pole_of_inaccessibility(),
+        -- The function allows a Lua table with options as parameter.
+        -- Currently only a single option "stretch" is defined. If you set
+        -- this to something other than 1 (the default), the functions doesn't
+        -- use a circle but an ellipse stretched horizontally by the specified
+        -- factor when calculating the result. This can make sense, because
+        -- labels are usually written horizontally so you need more space in
+        -- that direction. (Use a value > 1 to stretch horizontally or values
+        -- between 0 and 1 to stretch vertically. The value 0 ist not allowed.)
+        poi2 = geom:pole_of_inaccessibility({ stretch = 3 }),
+    })
+end
+
+function osm2pgsql.process_way(object)
+    if object.is_closed and object.tags.name then
+        local geom = object:as_polygon()
+        add(object.tags, geom)
+    end
+end
+
+function osm2pgsql.process_relation(object)
+    local relation_type = object:grab_tag('type')
+
+    if relation_type == 'multipolygon' and object.tags.name then
+        local geom = object:as_multipolygon()
+        -- The pole_of_inaccessibility() function only works for polygons,
+        -- not multipolygons. So we split up the multipolygons here and
+        -- calculate the pole for each part separately.
+        for g in geom:geometries() do
+            add(object.tags, g)
+        end
+    end
+end
+

src/CMakeLists.txt

@@ -10,6 +10,7 @@ target_sources(osm2pgsql_lib PRIVATE
   geom-box.cpp
   geom-from-osm.cpp
   geom-functions.cpp
+  geom-pole-of-inaccessibility.cpp
   input.cpp
   logging.cpp
   middle.cpp

src/flex-lua-geom.cpp

@@ -9,6 +9,7 @@
 
 #include "flex-lua-geom.hpp"
 #include "geom-functions.hpp"
+#include "geom-pole-of-inaccessibility.hpp"
 #include "lua-utils.hpp"
 
 extern "C"
@@ -167,6 +168,34 @@ static int geom_num_geometries(lua_State *lua_state)
     return 1;
 }
 
+static int geom_pole_of_inaccessibility(lua_State *lua_state)
+{
+    auto const *const input_geometry = unpack_geometry(lua_state);
+
+    double stretch = 1.0;
+    if (lua_gettop(lua_state) > 1) {
+        if (lua_type(lua_state, 2) != LUA_TTABLE) {
+            throw std::runtime_error{
+                "Argument #2 to 'pole_of_inaccessibility' must be a table."};
+        }
+
+        lua_getfield(lua_state, 2, "stretch");
+        if (lua_isnumber(lua_state, -1)) {
+            stretch = lua_tonumber(lua_state, -1);
+            if (stretch <= 0.0) {
+                throw std::runtime_error{"The 'stretch' factor must be > 0."};
+            }
+        } else {
+            throw std::runtime_error{"The 'stretch' factor must be a number."};
+        }
+    }
+
+    auto *geom = create_lua_geometry_object(lua_state);
+    geom::pole_of_inaccessibility(geom, *input_geometry, 0, stretch);
+
+    return 1;
+}
+
 static int geom_segmentize(lua_State *lua_state)
 {
     auto const *const input_geometry = unpack_geometry(lua_state);
@@ -257,6 +286,8 @@ void init_geometry_class(lua_State *lua_state)
     luaX_add_table_func(lua_state, "line_merge", geom_line_merge);
     luaX_add_table_func(lua_state, "reverse", geom_reverse);
     luaX_add_table_func(lua_state, "num_geometries", geom_num_geometries);
+    luaX_add_table_func(lua_state, "pole_of_inaccessibility",
+                        geom_pole_of_inaccessibility);
     luaX_add_table_func(lua_state, "segmentize", geom_segmentize);
     luaX_add_table_func(lua_state, "simplify", geom_simplify);
     luaX_add_table_func(lua_state, "srid", geom_srid);

src/geom-box.hpp

@@ -48,6 +48,11 @@ class box_t
     constexpr double width() const noexcept { return m_max_x - m_min_x; }
     constexpr double height() const noexcept { return m_max_y - m_min_y; }
 
+    constexpr point_t center() const noexcept
+    {
+        return {m_min_x + width() / 2.0, m_min_y + height() / 2.0};
+    }
+
     constexpr point_t min() const noexcept { return {m_min_x, m_min_y}; }
     constexpr point_t max() const noexcept { return {m_max_x, m_max_y}; }
 

src/geom-pole-of-inaccessibility.cpp

@@ -0,0 +1,259 @@
+/**
+ * SPDX-License-Identifier: GPL-2.0-or-later
+ *
+ * This file is part of osm2pgsql (https://osm2pgsql.org/).
+ *
+ * Copyright (C) 2006-2023 by the osm2pgsql developer community.
+ * For a full list of authors see the git log.
+ */
+
+#include "geom-pole-of-inaccessibility.hpp"
+#include "geom-boost-adaptor.hpp"
+#include "geom-box.hpp"
+#include "logging.hpp"
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <iostream>
+#include <queue>
+
+/**
+ * \file
+ *
+ * Implementation of the "Polylabel" algorithm for finding the "pole of
+ * inaccessibility", the internal point most distant from the polygon outline,
+ * or center of the maximum inscribed circle, described in
+ * https://blog.mapbox.com/a-new-algorithm-for-finding-a-visual-center-of-a-polygon-7c77e6492fbc
+ *
+ * Adapted from https://github.com/mapbox/polylabel and
+ * https://github.com/libgeos/geos/blob/main/src/algorithm/construct/MaximumInscribedCircle.cpp
+ * including the change from https://github.com/mapbox/polylabel/issues/82 .
+ *
+ * Forcing the precision to be no smaller than max(width, height) / 1000 of the
+ * envelope makes sure the algorithm terminates in sensible run-time and without
+ * taking too much memory. The value of 1000.0 was taken from the PostGIS
+ * implementation, but unlike the PostGIS implementation you can set a higher
+ * value.
+ */
+
+namespace geom {
+
+/// Get squared distance from a point p to a segment (a, b).
+static double point_to_segment_distance_squared(point_t p, point_t a, point_t b,
+                                                double stretch) noexcept
+{
+    double x = a.x();
+    double y = a.y() * stretch;
+    double dx = b.x() - x;
+    double dy = b.y() * stretch - y;
+
+    if (dx != 0 || dy != 0) {
+        double const t =
+            ((p.x() - x) * dx + (p.y() - y) * dy) / (dx * dx + dy * dy);
+
+        if (t > 1) {
+            x = b.x();
+            y = b.y() * stretch;
+        } else if (t > 0) {
+            x += dx * t;
+            y += dy * t;
+        }
+    }
+
+    dx = p.x() - x;
+    dy = p.y() - y;
+
+    return dx * dx + dy * dy;
+}
+
+/// Get squared distance from a point p to ring.
+static bool point_to_ring_distance_squared(point_t point, ring_t const &ring,
+                                           bool inside, double stretch,
+                                           double *min_dist_squared) noexcept
+{
+    std::size_t const len = ring.size();
+
+    for (std::size_t i = 0, j = len - 1; i < len; j = i++) {
+        auto const &a = ring[i];
+        auto const &b = ring[j];
+
+        if (((a.y() * stretch) > point.y()) !=
+                ((b.y() * stretch) > point.y()) &&
+            (point.x() < (b.x() - a.x()) * (point.y() - (a.y() * stretch)) /
+                                 ((b.y() - a.y()) * stretch) +
+                             a.x())) {
+            inside = !inside;
+        }
+
+        double const d =
+            point_to_segment_distance_squared(point, a, b, stretch);
+        if (d < *min_dist_squared) {
+            *min_dist_squared = d;
+        }
+    }
+
+    return inside;
+}
+
+/**
+ * Signed distance from point to polygon boundary. The result is negative if
+ * the point is outside.
+ */
+static auto point_to_polygon_distance(point_t point, polygon_t const &polygon,
+                                      double stretch)
+{
+    double min_dist_squared = std::numeric_limits<double>::infinity();
+
+    bool inside = point_to_ring_distance_squared(point, polygon.outer(), false,
+                                                 stretch, &min_dist_squared);
+
+    for (auto const &ring : polygon.inners()) {
+        inside = point_to_ring_distance_squared(point, ring, inside, stretch,
+                                                &min_dist_squared);
+    }
+
+    return (inside ? 1 : -1) * std::sqrt(min_dist_squared);
+}
+
+namespace {
+
+struct Cell
+{
+    static constexpr double const SQRT2 = 1.4142135623730951;
+
+    Cell(point_t c, double h, polygon_t const &polygon, double stretch)
+    : center(c), half_size(h),
+      dist(point_to_polygon_distance(center, polygon, stretch)),
+      max(dist + half_size * SQRT2)
+    {}
+
+    point_t center;   // cell center
+    double half_size; // half the cell size
+    double dist;      // distance from cell center to polygon
+    double max;       // max distance to polygon within a cell
+
+    friend bool operator<(Cell const &a, Cell const &b) noexcept
+    {
+        return a.max < b.max;
+    }
+};
+
+} // anonymous namespace
+
+static Cell make_centroid_cell(polygon_t const &polygon, double stretch)
+{
+    point_t centroid{0, 0};
+    boost::geometry::centroid(polygon, centroid);
+    centroid.set_y(stretch * centroid.y());
+    return {centroid, 0, polygon, stretch};
+}
+
+point_t pole_of_inaccessibility(const polygon_t &polygon, double precision,
+                                double stretch)
+{
+    assert(stretch > 0);
+
+    box_t const envelope = geom::envelope(polygon);
+
+    double const min_precision =
+        std::max(envelope.width(), envelope.height()) / 1000.0;
+    if (min_precision > precision) {
+        precision = min_precision;
+    }
+
+    box_t const stretched_envelope{envelope.min_x(), envelope.min_y() * stretch,
+                                   envelope.max_x(),
+                                   envelope.max_y() * stretch};
+
+    double const cell_size =
+        std::min(stretched_envelope.width(), stretched_envelope.height());
+    if (cell_size == 0) {
+        return envelope.min();
+    }
+
+    std::priority_queue<Cell, std::vector<Cell>> cell_queue;
+
+    // cover polygon with initial cells
+    {
+        double const h = cell_size / 2.0;
+        for (double x = stretched_envelope.min_x();
+             x < stretched_envelope.max_x(); x += cell_size) {
+            for (double y = stretched_envelope.min_y();
+                 y < stretched_envelope.max_y(); y += cell_size) {
+                cell_queue.emplace(point_t{x + h, y + h}, h, polygon, stretch);
+            }
+        }
+    }
+
+    // take centroid as the first best guess
+    auto best_cell = make_centroid_cell(polygon, stretch);
+
+    // second guess: bounding box centroid
+    Cell const bbox_cell{stretched_envelope.center(), 0, polygon, stretch};
+    if (bbox_cell.dist > best_cell.dist) {
+        best_cell = bbox_cell;
+    }
+
+    auto num_probes = cell_queue.size();
+    while (!cell_queue.empty()) {
+        // pick the most promising cell from the queue
+        auto cell = cell_queue.top();
+        cell_queue.pop();
+
+        // update the best cell if we found a better one
+        if (cell.dist > best_cell.dist) {
+            best_cell = cell;
+            log_debug("polyline: found best {} after {} probes",
+                      ::round(1e4 * cell.dist) / 1e4, num_probes);
+        }
+
+        // do not drill down further if there's no chance of a better solution
+        if (cell.max - best_cell.dist <= precision) {
+            continue;
+        }
+
+        // split the cell into four cells
+        auto const h = cell.half_size / 2.0;
+        auto const center = cell.center;
+
+        for (auto dy : {-h, h}) {
+            for (auto dx : {-h, h}) {
+                Cell c{point_t{center.x() + dx, center.y() + dy}, h, polygon,
+                       stretch};
+                if (c.max > best_cell.dist) {
+                    cell_queue.push(c);
+                }
+            }
+        }
+
+        num_probes += 4;
+    }
+
+    log_debug("polyline: num probes: {}", num_probes);
+    log_debug("polyline: best distance: {}", best_cell.dist);
+
+    return {best_cell.center.x(), best_cell.center.y() / stretch};
+}
+
+void pole_of_inaccessibility(geometry_t *output, geometry_t const &input,
+                             double precision, double stretch)
+{
+    if (input.is_polygon()) {
+        output->set<geom::point_t>() = pole_of_inaccessibility(
+            input.get<geom::polygon_t>(), precision, stretch);
+        output->set_srid(input.srid());
+    } else {
+        output->reset();
+    }
+}
+
+geometry_t pole_of_inaccessibility(geometry_t const &input, double precision,
+                                   double stretch)
+{
+    geometry_t geom;
+    pole_of_inaccessibility(&geom, input, precision, stretch);
+    return geom;
+}
+
+} // namespace geom