1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
|
#pragma once
// Original source from https://github.com/mapbox/polylabel, licensed
// under ISC.
//
// Adapted to use Boost Geometry instead of MapBox's geometry library.
//
// Changes:
// - Default precision changed from 1 to 0.00001.
// @mourner has some comments about what a reasonable precision value is
// for latitude/longitude coordinates, see
// https://github.com/mapbox/polylabel/issues/68#issuecomment-694906027
// https://github.com/mapbox/polylabel/issues/103#issuecomment-1516623862
//
// Possible future changes:
// - Port the change described in https://github.com/mapbox/polylabel/issues/33,
// implemented in Mapnik's Java renderer, to C++.
// - But see counterexample: https://github.com/mapbox/polylabel/pull/63
// @Fil also proposes an alternative approach there.
// - Pick precision as a function of the input geometry.
#include "geom.h"
#include <algorithm>
#include <cmath>
#include <iostream>
#include <queue>
namespace mapbox {
namespace detail {
// get squared distance from a point to a segment
double getSegDistSq(const Point& p,
const Point& a,
const Point& b) {
auto x = a.get<0>();
auto y = a.get<1>();
auto dx = b.get<0>() - x;
auto dy = b.get<1>() - y;
if (dx != 0 || dy != 0) {
auto t = ((p.get<0>() - x) * dx + (p.get<1>() - y) * dy) / (dx * dx + dy * dy);
if (t > 1) {
x = b.get<0>();
y = b.get<1>();
} else if (t > 0) {
x += dx * t;
y += dy * t;
}
}
dx = p.get<0>() - x;
dy = p.get<1>() - y;
return dx * dx + dy * dy;
}
// signed distance from point to polygon outline (negative if point is outside)
auto pointToPolygonDist(const Point& point, const Polygon& polygon) {
bool inside = false;
auto minDistSq = std::numeric_limits<double>::infinity();
{
const auto& ring = polygon.outer();
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.get<1>() > point.get<1>()) != (b.get<1>() > point.get<1>()) &&
(point.get<0>() < (b.get<0>() - a.get<0>()) * (point.get<1>() - a.get<1>()) / (b.get<1>() - a.get<1>()) + a.get<0>())) inside = !inside;
minDistSq = std::min(minDistSq, getSegDistSq(point, a, b));
}
}
for (const auto& ring : polygon.inners()) {
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.get<1>() > point.get<1>()) != (b.get<1>() > point.get<1>()) &&
(point.get<0>() < (b.get<0>() - a.get<0>()) * (point.get<1>() - a.get<1>()) / (b.get<1>() - a.get<1>()) + a.get<0>())) inside = !inside;
minDistSq = std::min(minDistSq, getSegDistSq(point, a, b));
}
}
return (inside ? 1 : -1) * std::sqrt(minDistSq);
}
struct Cell {
Cell(const Point& c_, double h_, const Polygon& polygon)
: c(c_),
h(h_),
d(pointToPolygonDist(c, polygon)),
max(d + h * std::sqrt(2))
{}
Point 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
Cell getCentroidCell(const Polygon& polygon) {
double area = 0;
double cx = 0, cy = 0;
const auto& ring = polygon.outer();
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.get<0>() * b.get<1>() - b.get<0>() * a.get<1>();
cx += (a.get<0>() + b.get<0>()) * f;
cy += (a.get<1>() + b.get<1>()) * f;
area += f * 3;
}
Point c { cx, cy };
return Cell(area == 0 ? ring.at(0) : Point { cx / area, cy / area }, 0, polygon);
}
} // namespace detail
Point polylabel(const Polygon& polygon, double precision = 0.00001, bool debug = false) {
using namespace detail;
// find the bounding box of the outer ring
Box envelope;
geom::envelope(polygon.outer(), envelope);
const Point size {
envelope.max_corner().get<0>() - envelope.min_corner().get<0>(),
envelope.max_corner().get<1>() - envelope.min_corner().get<1>()
};
const double cellSize = std::min(size.get<0>(), size.get<1>());
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.min_corner();
}
// cover polygon with initial cells
for (double x = envelope.min_corner().get<0>(); x < envelope.max_corner().get<0>(); x += cellSize) {
for (double y = envelope.min_corner().get<1>(); y < envelope.max_corner().get<1>(); y += cellSize) {
cellQueue.push(Cell({x + h, y + h}, h, polygon));
}
}
// take centroid as the first best guess
auto bestCell = getCentroidCell(polygon);
// second guess: bounding box centroid
Cell bboxCell(
Point {
envelope.min_corner().get<0>() + size.get<0>() / 2.0,
envelope.min_corner().get<1>() + size.get<1>() / 2.0
}, 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;
if (debug) std::cout << "found best " << ::round(1e4 * cell.d) / 1e4 << " after " << numProbes << " probes" << std::endl;
}
// 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.get<0>() - h, cell.c.get<1>() - h}, h, polygon));
cellQueue.push(Cell({cell.c.get<0>() + h, cell.c.get<1>() - h}, h, polygon));
cellQueue.push(Cell({cell.c.get<0>() - h, cell.c.get<1>() + h}, h, polygon));
cellQueue.push(Cell({cell.c.get<0>() + h, cell.c.get<1>() + h}, h, polygon));
numProbes += 4;
}
if (debug) {
std::cout << "num probes: " << numProbes << std::endl;
std::cout << "best distance: " << bestCell.d << std::endl;
}
return bestCell.c;
}
} // namespace mapbox
|
