New interior algorithm

include/mapnik/geom_util.hpp

@@ -518,80 +518,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,35 @@
+/*****************************************************************************
+ *
+ * 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.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/geometry/polygon_vertex_processor.hpp

@@ -0,0 +1,71 @@
+/*****************************************************************************
+ *
+ * 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_POLYGON_VERTEX_PROCESSOR_HPP
+#define MAPNIK_GEOMETRY_POLYGON_VERTEX_PROCESSOR_HPP
+
+// geometry
+#include <mapnik/geometry.hpp>
+
+namespace mapnik { namespace geometry {
+
+template <typename T>
+struct polygon_vertex_processor
+{
+    template <typename Path>
+    void add_path(Path & path)
+    {
+        point<T> p;
+        unsigned cmd;
+        linear_ring<T> ring;
+        bool exterior = true;
+        while ((cmd = path.vertex(&p.x, &p.y)) != SEG_END)
+        {
+            switch (cmd)
+            {
+                case SEG_MOVETO:
+                case SEG_LINETO:
+                    ring.emplace_back(p);
+                    break;
+                case SEG_CLOSE:
+                    ring.emplace_back(ring.front());
+                    if (exterior)
+                    {
+                        polygon_.exterior_ring = std::move(ring);
+                        exterior = false;
+                    }
+                    else
+                    {
+                        polygon_.interior_rings.emplace_back(std::move(ring));
+                    }
+                    ring = linear_ring<T>();
+                    break;
+            }
+        }
+    }
+
+    polygon<T> polygon_;
+};
+
+} }
+
+#endif // MAPNIK_GEOMETRY_POLYGON_VERTEX_PROCESSOR_HPP

include/mapnik/markers_placements/basic.hpp

@@ -46,6 +46,7 @@ struct markers_placement_params
     bool allow_overlap;
     bool avoid_edges;
     direction_enum direction;
+    double scale_factor;
 };
 
 class markers_basic_placement : util::noncopyable

include/mapnik/markers_placements/interior.hpp

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

include/mapnik/renderer_common/process_point_symbolizer.hpp

@@ -32,6 +32,7 @@
 #include <mapnik/geometry_centroid.hpp>
 #include <mapnik/geometry_type.hpp>
 #include <mapnik/geometry_types.hpp>
+#include <mapnik/geometry/interior.hpp>
 #include <mapnik/vertex_adapters.hpp>
 #include <mapnik/geom_util.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);
+            pt = geometry::interior(poly, common.scale_factor_);
-            if (!label::interior_position(va ,pt.x, pt.y))
-                return;
         }
         else
         {

src/build.py

@@ -169,6 +169,7 @@ source = Split(
     datasource_cache_static.cpp
     debug.cpp
     geometry_reprojection.cpp
+    geometry/interior.cpp
     expression_node.cpp
     expression_string.cpp
     expression.cpp

src/geometry/interior.cpp

@@ -0,0 +1,234 @@
+/*****************************************************************************
+ *
+ * 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/box2d.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;
+}
+
+template <class T>
+void point_to_ring_dist(point<T> const& point, linear_ring<T> const& ring,
+                        bool & inside, double & min_dist_sq)
+{
+    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));
+    }
+}
+
+// signed distance from point to polygon outline (negative if point is outside)
+template <class T>
+double point_to_polygon_dist(const point<T>& point, const polygon<T>& polygon)
+{
+    bool inside = false;
+    double min_dist_sq = std::numeric_limits<double>::infinity();
+
+    point_to_ring_dist(point, polygon.exterior_ring, inside, min_dist_sq);
+
+    for (const auto& ring : polygon.interior_rings)
+    {
+        point_to_ring_dist(point, ring, inside, min_dist_sq);
+    }
+
+    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()(point<T> const& cell_center, T distance_polygon) const
+    {
+        if (distance_polygon <= 0)
+        {
+            return distance_polygon;
+        }
+        point<T> d(cell_center.x - centroid.x, cell_center.y - centroid.y);
+        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.exterior_ring);
+    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);
+
+} }
+

src/renderer_common/render_markers_symbolizer.cpp

@@ -228,7 +228,8 @@ markers_dispatch_params::markers_dispatch_params(box2d<double> const& size,
         get<value_double, keys::max_error>(sym, feature, vars),
         get<value_bool, keys::allow_overlap>(sym, feature, vars),
         get<value_bool, keys::avoid_edges>(sym, feature, vars),
-        get<direction_enum, keys::direction>(sym, feature, vars)}
+        get<direction_enum, keys::direction>(sym, feature, vars),
+        scale}
     , placement_method(get<marker_placement_enum, keys::markers_placement_type>(sym, feature, vars))
     , ignore_placement(get<value_bool, keys::ignore_placement>(sym, feature, vars))
     , snap_to_pixels(snap)

src/text/symbolizer_helpers.cpp

@@ -31,6 +31,7 @@
 #include <mapnik/geometry.hpp>
 #include <mapnik/geometry_type.hpp>
 #include <mapnik/geometry_centroid.hpp>
+#include <mapnik/geometry/interior.hpp>
 #include <mapnik/vertex_processor.hpp>
 #include <mapnik/geom_util.hpp>
 #include <mapnik/parse_path.hpp>
@@ -40,6 +41,10 @@
 #include <mapnik/text/placement_finder_impl.hpp>
 #include <mapnik/text/placements/base.hpp>
 #include <mapnik/text/placements/dummy.hpp>
+#include <mapnik/geometry_transform.hpp>
+#include <mapnik/geometry_strategy.hpp>
+#include <mapnik/proj_strategy.hpp>
+#include <mapnik/view_strategy.hpp>
 
 namespace mapnik {
 namespace geometry {
@@ -288,9 +293,16 @@ void base_symbolizer_helper::initialize_points() const
             }
             else if (how_placed == INTERIOR_PLACEMENT && type == geometry::geometry_types::Polygon)
             {
-                auto const& poly = mapnik::util::get<geometry::polygon<double> >(geom);
+                auto const& poly = util::get<geometry::polygon<double>>(geom);
-                geometry::polygon_vertex_adapter<double> va(poly);
+                proj_transform backwart_transform(prj_trans_.dest(), prj_trans_.source());
-                success = label::interior_position(va, label_x, label_y);
+                view_strategy vs(t_);
+                proj_strategy ps(backwart_transform);
+                using transform_group_type = geometry::strategy_group<proj_strategy, view_strategy>;
+                transform_group_type transform_group(ps, vs);
+                geometry::polygon<double> tranformed_poly(geometry::transform<double>(poly, transform_group));
+                geometry::point<double> pt = geometry::interior(tranformed_poly, scale_factor_);
+                points_.emplace_back(pt.x, pt.y);
+                continue;
             }
             else
             {