1#pragma once
2
3#include <mapbox/geometry/polygon.hpp>
4#include <mapbox/geometry/envelope.hpp>
5#include <mapbox/geometry/point.hpp>
6#include <mapbox/geometry/point_arithmetic.hpp>
7
8#include <algorithm>
9#include <cmath>
10#include <iostream>
11#include <queue>
12
13namespace mapbox {
14
15namespace detail {
16
17// get squared distance from a point to a segment
18template <class T>
19T getSegDistSq(const geometry::point<T>& p,
20 const geometry::point<T>& a,
21 const geometry::point<T>& b) {
22 auto x = a.x;
23 auto y = a.y;
24 auto dx = b.x - x;
25 auto dy = b.y - y;
26
27 if (dx != 0 || dy != 0) {
28
29 auto t = ((p.x - x) * dx + (p.y - y) * dy) / (dx * dx + dy * dy);
30
31 if (t > 1) {
32 x = b.x;
33 y = b.y;
34
35 } else if (t > 0) {
36 x += dx * t;
37 y += dy * t;
38 }
39 }
40
41 dx = p.x - x;
42 dy = p.y - y;
43
44 return dx * dx + dy * dy;
45}
46
47// signed distance from point to polygon outline (negative if point is outside)
48template <class T>
49auto pointToPolygonDist(const geometry::point<T>& point, const geometry::polygon<T>& polygon) {
50 bool inside = false;
51 auto minDistSq = std::numeric_limits<double>::infinity();
52
53 for (const auto& ring : polygon) {
54 for (std::size_t i = 0, len = ring.size(), j = len - 1; i < len; j = i++) {
55 const auto& a = ring[i];
56 const auto& b = ring[j];
57
58 if ((a.y > point.y) != (b.y > point.y) &&
59 (point.x < (b.x - a.x) * (point.y - a.y) / (b.y - a.y) + a.x)) inside = !inside;
60
61 minDistSq = std::min(minDistSq, getSegDistSq(point, a, b));
62 }
63 }
64
65 return (inside ? 1 : -1) * std::sqrt(x: minDistSq);
66}
67
68template <class T>
69struct Cell {
70 Cell(const geometry::point<T>& c_, T h_, const geometry::polygon<T>& polygon)
71 : c(c_),
72 h(h_),
73 d(pointToPolygonDist(c, polygon)),
74 max(d + h * std::sqrt(x: 2))
75 {}
76
77 geometry::point<T> c; // cell center
78 T h; // half the cell size
79 T d; // distance from cell center to polygon
80 T max; // max distance to polygon within a cell
81};
82
83// get polygon centroid
84template <class T>
85Cell<T> getCentroidCell(const geometry::polygon<T>& polygon) {
86 T area = 0;
87 geometry::point<T> c { 0, 0 };
88 const auto& ring = polygon.at(0);
89
90 for (std::size_t i = 0, len = ring.size(), j = len - 1; i < len; j = i++) {
91 const geometry::point<T>& a = ring[i];
92 const geometry::point<T>& b = ring[j];
93 auto f = a.x * b.y - b.x * a.y;
94 c.x += (a.x + b.x) * f;
95 c.y += (a.y + b.y) * f;
96 area += f * 3;
97 }
98
99 return Cell<T>(area == 0 ? ring.at(0) : c / area, 0, polygon);
100}
101
102} // namespace detail
103
104template <class T>
105geometry::point<T> polylabel(const geometry::polygon<T>& polygon, T precision = 1, bool debug = false) {
106 using namespace detail;
107
108 // find the bounding box of the outer ring
109 const geometry::box<T> envelope = geometry::envelope(polygon.at(0));
110
111 const geometry::point<T> size {
112 envelope.max.x - envelope.min.x,
113 envelope.max.y - envelope.min.y
114 };
115
116 const T cellSize = std::min(size.x, size.y);
117 T h = cellSize / 2;
118
119 // a priority queue of cells in order of their "potential" (max distance to polygon)
120 auto compareMax = [] (const Cell<T>& a, const Cell<T>& b) {
121 return a.max < b.max;
122 };
123 using Queue = std::priority_queue<Cell<T>, std::vector<Cell<T>>, decltype(compareMax)>;
124 Queue cellQueue(compareMax);
125
126 if (cellSize == 0) {
127 return envelope.min;
128 }
129
130 // cover polygon with initial cells
131 for (T x = envelope.min.x; x < envelope.max.x; x += cellSize) {
132 for (T y = envelope.min.y; y < envelope.max.y; y += cellSize) {
133 cellQueue.push(Cell<T>({x + h, y + h}, h, polygon));
134 }
135 }
136
137 // take centroid as the first best guess
138 auto bestCell = getCentroidCell(polygon);
139
140 // special case for rectangular polygons
141 Cell<T> bboxCell(envelope.min + size / 2.0, 0, polygon);
142 if (bboxCell.d > bestCell.d) {
143 bestCell = bboxCell;
144 }
145
146 auto numProbes = cellQueue.size();
147 while (!cellQueue.empty()) {
148 // pick the most promising cell from the queue
149 auto cell = cellQueue.top();
150 cellQueue.pop();
151
152 // update the best cell if we found a better one
153 if (cell.d > bestCell.d) {
154 bestCell = cell;
155 if (debug) std::cout << "found best " << ::round(x: 1e4 * cell.d) / 1e4 << " after " << numProbes << " probes" << std::endl;
156 }
157
158 // do not drill down further if there's no chance of a better solution
159 if (cell.max - bestCell.d <= precision) continue;
160
161 // split the cell into four cells
162 h = cell.h / 2;
163 cellQueue.push(Cell<T>({cell.c.x - h, cell.c.y - h}, h, polygon));
164 cellQueue.push(Cell<T>({cell.c.x + h, cell.c.y - h}, h, polygon));
165 cellQueue.push(Cell<T>({cell.c.x - h, cell.c.y + h}, h, polygon));
166 cellQueue.push(Cell<T>({cell.c.x + h, cell.c.y + h}, h, polygon));
167 numProbes += 4;
168 }
169
170 if (debug) {
171 std::cout << "num probes: " << numProbes << std::endl;
172 std::cout << "best distance: " << bestCell.d << std::endl;
173 }
174
175 return bestCell.c;
176}
177
178} // namespace mapbox
179

source code of qtlocation/src/3rdparty/mapbox-gl-native/deps/polylabel/1.0.3/include/mapbox/polylabel.hpp