new interior algorithm

include/mapnik/geom_util.hpp

@@ -537,80 +537,6 @@ bool hit_test(PathType & path, double x, double y, double tol)
     return inside;
 }
 
-template <typename PathType>
-bool interior_position(PathType & path, double & x, double & y)
-{
-    // start with the centroid
-    if (!label::centroid(path, x,y))
-        return false;
-
-    // otherwise we find a horizontal line across the polygon and then return the
-    // center of the widest intersection between the polygon and the line.
-
-    std::vector<double> intersections; // only need to store the X as we know the y
-    geometry::point<double> p0, p1, move_to;
-    unsigned command = SEG_END;
-
-    path.rewind(0);
-
-    while (SEG_END != (command = path.vertex(&p0.x, &p0.y)))
-    {
-        switch (command)
-        {
-            case SEG_MOVETO:
-                move_to = p0;
-                break;
-            case SEG_CLOSE:
-                p0 = move_to;
-            case SEG_LINETO:
-                // if the segments overlap
-                if (p0.y == p1.y)
-                {
-                    if (p0.y == y)
-                    {
-                        double xi = (p0.x + p1.x) / 2.0;
-                        intersections.push_back(xi);
-                    }
-                }
-                // if the path segment crosses the bisector
-                else if ((p0.y <= y && p1.y >= y) ||
-                         (p0.y >= y && p1.y <= y))
-                {
-                    // then calculate the intersection
-                    double xi = p0.x;
-                    if (p0.x != p1.x)
-                    {
-                        double m = (p1.y - p0.y) / (p1.x - p0.x);
-                        double c = p0.y - m * p0.x;
-                        xi = (y - c) / m;
-                    }
-
-                    intersections.push_back(xi);
-                }
-                break;
-        }
-        p1 = p0;
-    }
-
-    // no intersections we just return the default
-    if (intersections.empty())
-        return true;
-    std::sort(intersections.begin(), intersections.end());
-    double max_width = 0;
-    for (unsigned ii = 1; ii < intersections.size(); ii += 2)
-    {
-        double xlow = intersections[ii-1];
-        double xhigh = intersections[ii];
-        double width = xhigh - xlow;
-        if (width > max_width)
-        {
-            x = (xlow + xhigh) / 2.0;
-            max_width = width;
-        }
-    }
-    return true;
-}
-
 }}
 
 #endif // MAPNIK_GEOM_UTIL_HPP

include/mapnik/geometry/interior.hpp

@@ -0,0 +1,36 @@
+/*****************************************************************************
+ *
+ * This file is part of Mapnik (c++ mapping toolkit)
+ *
+ * Copyright (C) 2017 Artem Pavlenko
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ *
+ *****************************************************************************/
+
+#ifndef MAPNIK_GEOMETRY_INTERIOR_HPP
+#define MAPNIK_GEOMETRY_INTERIOR_HPP
+
+#include <mapnik/geometry/polygon.hpp>
+#include <mapnik/geometry/point.hpp>
+
+namespace mapnik { namespace geometry {
+
+template <class T>
+point<T> interior(polygon<T> const& polygon, double scale_factor);
+
+} }
+
+#endif // MAPNIK_GEOMETRY_INTERIOR_HPP

include/mapnik/markers_placements/interior.hpp

@@ -26,6 +26,8 @@
 #include <mapnik/markers_placements/point.hpp>
 #include <mapnik/geom_util.hpp>
 #include <mapnik/geometry/geometry_types.hpp>
+#include <mapnik/geometry/interior.hpp>
+#include <mapnik/geometry/polygon_vertex_processor.hpp>
 
 namespace mapnik {
 
@@ -58,11 +60,13 @@ class markers_interior_placement : public markers_point_placement<Locator, Detec
         }
         else
         {
-            if (!label::interior_position(this->locator_, x, y))
-            {
-                this->done_ = true;
-                return false;
-            }
+            geometry::polygon_vertex_processor<double> vertex_processor;
+            vertex_processor.add_path(this->locator_);
+            geometry::point<double> placement = geometry::interior(vertex_processor.polygon_,
+                                                                   this->params_.scale_factor);
+
+            x = placement.x;
+            y = placement.y;
         }
 
         angle = 0;

include/mapnik/renderer_common/process_point_symbolizer.hpp

@@ -30,6 +30,7 @@
 #include <mapnik/marker_cache.hpp>
 #include <mapnik/label_collision_detector.hpp>
 #include <mapnik/geometry/centroid.hpp>
+#include <mapnik/geometry/interior.hpp>
 #include <mapnik/geometry/geometry_type.hpp>
 #include <mapnik/geometry/geometry_types.hpp>
 #include <mapnik/vertex_adapters.hpp>
@@ -79,9 +80,7 @@ void render_point_symbolizer(point_symbolizer const &sym,
         else if (type == mapnik::geometry::geometry_types::Polygon)
         {
             auto const& poly = mapnik::util::get<geometry::polygon<double> >(geometry);
-            geometry::polygon_vertex_adapter<double> va(poly);
-            if (!label::interior_position(va ,pt.x, pt.y))
-                return;
+            pt = geometry::interior(poly, common.scale_factor_);
         }
         else
         {

src/build.py

@@ -173,6 +173,7 @@ def ldconfig(*args,**kwargs):
     geometry/closest_point.cpp
     geometry/reprojection.cpp
     geometry/envelope.cpp
+    geometry/interior.cpp
     geometry/polylabel.cpp
     expression_node.cpp
     expression_string.cpp

src/geometry/interior.cpp

@@ -0,0 +1,226 @@
+/*****************************************************************************
+ *
+ * This file is part of Mapnik (c++ mapping toolkit)
+ *
+ * Copyright (C) 2017 Artem Pavlenko
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ *
+ *****************************************************************************/
+
+#include <mapnik/geometry/interior.hpp>
+#include <mapnik/geometry/envelope.hpp>
+#include <mapnik/geometry/box2d.hpp>
+#include <mapbox/geometry/point_arithmetic.hpp>
+#include <mapnik/geometry/centroid.hpp>
+
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+#include <queue>
+
+namespace mapnik { namespace geometry {
+
+// Interior algorithm is realized as a modification of Polylabel algorithm
+// from https://github.com/mapbox/polylabel.
+// The modification aims to improve visual output by prefering
+// placements closer to centroid.
+
+namespace detail {
+
+// get squared distance from a point to a segment
+template <class T>
+T segment_dist_sq(const point<T>& p,
+                  const point<T>& a,
+                  const point<T>& b)
+{
+    auto x = a.x;
+    auto y = a.y;
+    auto dx = b.x - x;
+    auto dy = b.y - y;
+
+    if (dx != 0 || dy != 0) {
+
+        auto t = ((p.x - x) * dx + (p.y - y) * dy) / (dx * dx + dy * dy);
+
+        if (t > 1) {
+            x = b.x;
+            y = b.y;
+
+        } else if (t > 0) {
+            x += dx * t;
+            y += dy * t;
+        }
+    }
+
+    dx = p.x - x;
+    dy = p.y - y;
+
+    return dx * dx + dy * dy;
+}
+
+// signed distance from point to polygon outline (negative if point is outside)
+template <class T>
+auto point_to_polygon_dist(const point<T>& point, const polygon<T>& polygon)
+{
+    bool inside = false;
+    auto min_dist_sq = 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.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;
+
+            min_dist_sq = std::min(min_dist_sq, segment_dist_sq(point, a, b));
+        }
+    }
+
+    return (inside ? 1 : -1) * std::sqrt(min_dist_sq);
+}
+
+template <class T>
+struct fitness_functor
+{
+    fitness_functor(point<T> const& centroid, point<T> const& polygon_size)
+        : centroid(centroid),
+          max_size(std::max(polygon_size.x, polygon_size.y))
+        {}
+
+    T operator()(const point<T>& cell_center, T distance_polygon) const
+    {
+        if (distance_polygon <= 0)
+        {
+            return distance_polygon;
+        }
+        point<T> d = cell_center - centroid;
+        double distance_centroid = std::sqrt(d.x * d.x + d.y * d.y);
+        return distance_polygon * (1 - distance_centroid / max_size);
+    }
+
+    point<T> centroid;
+    T max_size;
+};
+
+template <class T>
+struct cell
+{
+    template <class FitnessFunc>
+    cell(const point<T>& c_, T h_,
+         const polygon<T>& polygon,
+         const FitnessFunc& ff)
+        : c(c_),
+          h(h_),
+          d(point_to_polygon_dist(c, polygon)),
+          fitness(ff(c, d)),
+          max_fitness(ff(c, d + h * std::sqrt(2)))
+        {}
+
+    point<T> c; // cell center
+    T h; // half the cell size
+    T d; // distance from cell center to polygon
+    T fitness; // fitness of the cell center
+    T max_fitness; // a "potential" of the cell calculated from max distance to polygon within the cell
+};
+
+template <class T>
+point<T> polylabel(const polygon<T>& polygon, T precision = 1)
+{
+    // find the bounding box of the outer ring
+    const box2d<T> bbox = envelope(polygon.at(0));
+    const point<T> size { bbox.width(), bbox.height() };
+
+    const T cell_size = std::min(size.x, size.y);
+    T h = cell_size / 2;
+
+    // a priority queue of cells in order of their "potential" (max distance to polygon)
+    auto compare_func = [] (const cell<T>& a, const cell<T>& b)
+    {
+        return a.max_fitness < b.max_fitness;
+    };
+    using Queue = std::priority_queue<cell<T>, std::vector<cell<T>>, decltype(compare_func)>;
+    Queue queue(compare_func);
+
+    if (cell_size == 0)
+    {
+        return { bbox.minx(), bbox.miny() };
+    }
+
+    point<T> centroid;
+    if (!mapnik::geometry::centroid(polygon, centroid))
+    {
+        auto center = bbox.center();
+        return { center.x, center.y };
+    }
+
+    fitness_functor<T> fitness_func(centroid, size);
+
+    // cover polygon with initial cells
+    for (T x = bbox.minx(); x < bbox.maxx(); x += cell_size)
+    {
+        for (T y = bbox.miny(); y < bbox.maxy(); y += cell_size)
+        {
+            queue.push(cell<T>({x + h, y + h}, h, polygon, fitness_func));
+        }
+    }
+
+    // take centroid as the first best guess
+    auto best_cell = cell<T>(centroid, 0, polygon, fitness_func);
+
+    while (!queue.empty())
+    {
+        // pick the most promising cell from the queue
+        auto current_cell = queue.top();
+        queue.pop();
+
+        // update the best cell if we found a better one
+        if (current_cell.fitness > best_cell.fitness)
+        {
+            best_cell = current_cell;
+        }
+
+        // do not drill down further if there's no chance of a better solution
+        if (current_cell.max_fitness - best_cell.fitness <= precision) continue;
+
+        // split the cell into four cells
+        h = current_cell.h / 2;
+        queue.push(cell<T>({current_cell.c.x - h, current_cell.c.y - h}, h, polygon, fitness_func));
+        queue.push(cell<T>({current_cell.c.x + h, current_cell.c.y - h}, h, polygon, fitness_func));
+        queue.push(cell<T>({current_cell.c.x - h, current_cell.c.y + h}, h, polygon, fitness_func));
+        queue.push(cell<T>({current_cell.c.x + h, current_cell.c.y + h}, h, polygon, fitness_func));
+    }
+
+    return best_cell.c;
+}
+
+} // namespace detail
+
+template <class T>
+point<T> interior(polygon<T> const& polygon, double scale_factor)
+{
+    // This precision has been chosen to work well in the map (viewport) coordinates.
+    double precision = 10.0 * scale_factor;
+    return detail::polylabel(polygon, precision);
+}
+
+template
+point<double> interior(polygon<double> const& polygon, double scale_factor);
+
+} }
+