1 #include <mbgl/util/tile_cover_impl.hpp>
2 #include <mbgl/util/tile_coordinate.hpp>
3 
4 #include <functional>
5 #include <cmath>
6 #include <assert.h>
7 #include <limits.h>
8 #include <algorithm>
9 
10 namespace mbgl {
11 namespace util {
12 
13 using PointList = std::vector<Point<double>>;
14 
15 struct TileSpan {
16     int32_t xmin, xmax;
17     bool winding;
18 };
19 
20 // Reorder a ring of points such that it starts at a point with a local minimum y-coordinate
start_list_on_local_minimum(PointList & points)21 void start_list_on_local_minimum(PointList& points) {
22     auto prev_pt = std::prev(points.end(), 2);
23     auto pt = points.begin();
24     auto next_pt = std::next(pt);
25     while (pt != points.end()) {
26         if ((pt->y <= prev_pt->y) &&
27             (pt->y < next_pt->y)) {
28             break;
29         }
30         prev_pt = pt;
31         pt++;
32         next_pt++;
33         if (next_pt == points.end()) { next_pt = std::next(points.begin()); }
34     }
35     if (pt == points.end())
36         return;
37     //Re-close linear rings with first_pt = last_pt
38     if (points.back() == points.front()) {
39         points.pop_back();
40     }
41     std::rotate(points.begin(), pt, points.end());
42     points.push_back(*points.begin());
43 }
44 
45 //Create a bound towards a local maximum point, starting from pt.
46 // Traverse from current pt until the next pt changes y-direction, and copy
47 // all points from start to end (inclusive) into a Bound.
create_bound_towards_maximum(PointList & points,PointList::iterator & pt)48 Bound create_bound_towards_maximum(PointList& points, PointList::iterator& pt) {
49     if (std::distance(pt, points.end()) < 2) { return {}; }
50 
51     const auto begin = pt;
52     auto next_pt = std::next(begin);
53     while (pt->y <= next_pt->y) {
54         pt++;
55         next_pt++;
56         if (next_pt == points.end()) { pt++; break; }
57     }
58 
59     const auto pt_distance = std::distance(begin, next_pt);
60     if (pt_distance < 2) {
61         return {};
62     }
63 
64     Bound bnd;
65     bnd.points.reserve(static_cast<std::size_t>(std::distance(begin, next_pt)));
66     std::copy(begin, next_pt, std::back_inserter(bnd.points));
67     bnd.winding = true;
68     return bnd;
69 }
70 
71 //Create a bound towards a local minimum point, starting from pt.
72 // Traverse from current pt until the next pt changes y-direction, and copy
73 // all points from start to end (inclusive) into a Bound.
create_bound_towards_minimum(PointList & points,PointList::iterator & pt)74 Bound create_bound_towards_minimum(PointList& points, PointList::iterator& pt) {
75     if (std::distance(pt, points.end()) < 2) { return {}; }
76 
77     auto begin = pt;
78     auto next_pt = std::next(begin);
79     while (pt->y > next_pt->y) {
80         pt++;
81         next_pt++;
82         if (next_pt == points.end()) { pt++; break; }
83     }
84 
85     const auto pt_distance = std::distance(begin, next_pt);
86     if (pt_distance < 2) {
87         return {};
88     }
89     Bound bnd;
90     bnd.points.reserve(static_cast<std::size_t>(std::distance(begin, next_pt)));
91     //For bounds that start at a max, reverse copy so that all bounds start at a min
92     std::reverse_copy(begin, next_pt, std::back_inserter(bnd.points));
93     bnd.winding = false;
94     return bnd;
95 }
96 
97 // Given a set of points (ring or list) representing a shape, compute a set of
98 // Bounds, where each Bound represents edges going from a local minima to a local
99 // maxima point. The BoundsMap is an edge table indexed on the starting Y-tile
100 // of each Bound.
build_bounds_map(PointList & points,uint32_t maxTile,BoundsMap & et,bool closed=false)101 void build_bounds_map(PointList& points, uint32_t maxTile, BoundsMap& et, bool closed = false) {
102     if (points.size() < 2) return;
103     //While traversing closed rings, start the bounds at a local minimum.
104     // (For linestrings the starting point is always a local maxima/minima)
105     if (closed) {
106         start_list_on_local_minimum(points);
107     }
108 
109     auto pointsIter = points.begin();
110     while (pointsIter != points.end()) {
111         Bound to_max = create_bound_towards_maximum(points, pointsIter);
112         Bound to_min = create_bound_towards_minimum(points, pointsIter);
113 
114         if (to_max.points.size() >= 2) {
115             // Projections may result in values beyond the bounds, clamp to max tile coordinates
116             const auto y = static_cast<uint32_t>(std::floor(clamp(to_max.points.front().y, 0.0, (double)maxTile)));
117             et[y].push_back(to_max);
118         }
119         if (to_min.points.size() >= 2) {
120             const auto y = static_cast<uint32_t>(std::floor(clamp(to_min.points.front().y, 0.0, (double)maxTile)));
121             et[y].push_back(to_min);
122         }
123     }
124     assert(pointsIter == points.end());
125 }
126 
update_span(TileSpan & xp,double x)127 void update_span(TileSpan& xp, double x) {
128     xp.xmin = std::min(xp.xmin, static_cast<int32_t>(std::floor(x)));
129     xp.xmax = std::max(xp.xmax, static_cast<int32_t>(std::ceil(x)));
130 }
131 
132 // Use the active bounds, an accumulation of all bounds that enter the y tile row,
133 // or start in that row.
134 // Iterate all points of a bound until it exits the row (or ends) and compute the
135 // set of X tiles it spans across. The winding direction of the bound is also
136 // captured for each span to later fill tiles between bounds for polygons
scan_row(uint32_t y,Bounds & activeBounds)137 std::vector<TileSpan> scan_row(uint32_t y, Bounds& activeBounds) {
138     std::vector<TileSpan> tile_range;
139     tile_range.reserve(activeBounds.size());
140 
141     for(Bound& b: activeBounds) {
142         TileSpan xp = { INT_MAX, 0, b.winding };
143         double x;
144         const auto numEdges = b.points.size() - 1;
145         while (b.currentPoint < numEdges) {
146             x = b.interpolate(y);
147             update_span(xp, x);
148 
149             // If this edge ends beyond the current row, find the x-intercept where
150             // it exits the row
151             auto& p1 = b.points[b.currentPoint + 1];
152             if (p1.y > y+1) {
153                 x = b.interpolate(y+1);
154                 update_span(xp, x);
155                 break;
156             } else if (b.currentPoint == numEdges - 1) {
157                 // For last edge, consider x-intercept at the end of the edge.
158                 x = p1.x;
159                 update_span(xp, x);
160             }
161             b.currentPoint++;
162         }
163         tile_range.push_back(xp);
164     }
165     // Erase bounds in the active table whose current edge ends inside this row,
166     // or there are no more edges
167     auto bound = activeBounds.begin();
168     while (bound != activeBounds.end()) {
169         if ( bound->currentPoint == bound->points.size() - 1 &&
170             bound->points[bound->currentPoint].y <= y+1) {
171             bound = activeBounds.erase(bound);
172         } else {
173             bound++;
174         }
175     }
176     // Sort the X-extents of each crossing bound by x_min, x_max
177     std::sort(tile_range.begin(), tile_range.end(), [] (TileSpan& a, TileSpan& b) {
178         return std::tie(a.xmin, a.xmax) < std::tie(b.xmin, b.xmax);
179     });
180 
181     return tile_range;
182 }
183 
184 struct BuildBoundsMap {
185     int32_t zoom;
186     bool project = false;
BuildBoundsMapmbgl::util::BuildBoundsMap187     BuildBoundsMap(int32_t z, bool p): zoom(z), project(p) {}
188 
buildTablembgl::util::BuildBoundsMap189     void buildTable(const std::vector<Point<double>>& points, BoundsMap& et, bool closed = false) const {
190         PointList projectedPoints;
191         if (project) {
192             projectedPoints.reserve(points.size());
193             for(const auto&p : points) {
194                 projectedPoints.push_back(
195                     Projection::project(LatLng{ p.y, p.x }, zoom));
196             }
197         } else {
198             projectedPoints.insert(projectedPoints.end(), points.begin(), points.end());
199         }
200         build_bounds_map(projectedPoints, 1 << zoom, et, closed);
201     }
202 
buildPolygonTablembgl::util::BuildBoundsMap203     void buildPolygonTable(const Polygon<double>& polygon, BoundsMap& et) const {
204         for(const auto&ring : polygon) {
205             buildTable(ring, et, true);
206         }
207     }
operator ()mbgl::util::BuildBoundsMap208     BoundsMap operator()(const Point<double>&p) const {
209         Bound bnd;
210         auto point = p;
211         if (project) {
212             point = Projection::project(LatLng{p.y, p.x}, zoom);
213         }
214         bnd.points.insert(bnd.points.end(), 2, point);
215         bnd.winding = false;
216         BoundsMap et;
217         const auto y = static_cast<uint32_t>(std::floor(clamp(point.y, 0.0, (double)(1 << zoom))));
218         et[y].push_back(bnd);
219         return et;
220     }
221 
operator ()mbgl::util::BuildBoundsMap222     BoundsMap operator()(const MultiPoint<double>& points) const {
223         BoundsMap et;
224         for (const Point<double>& p: points) {
225             Bound bnd;
226             auto point = p;
227             if (project) {
228                 point = Projection::project(LatLng{p.y, p.x}, zoom);
229             }
230             bnd.points.insert(bnd.points.end(), 2, point);
231             bnd.winding = false;
232             const auto y = static_cast<uint32_t>(std::floor(clamp(point.y, 0.0, (double)(1 << zoom))));
233             et[y].push_back(bnd);
234         }
235         return et;
236     }
237 
operator ()mbgl::util::BuildBoundsMap238     BoundsMap operator()(const LineString<double>& lines) const {
239         BoundsMap et;
240         buildTable(lines, et);
241         return et;
242     }
243 
operator ()mbgl::util::BuildBoundsMap244     BoundsMap operator()(const MultiLineString<double>& lines) const {
245         BoundsMap et;
246         for(const auto&line : lines) {
247             buildTable(line, et);
248         }
249         return et;
250     }
251 
operator ()mbgl::util::BuildBoundsMap252     BoundsMap operator()(const Polygon<double>& polygon) const {
253         BoundsMap et;
254         buildPolygonTable(polygon, et);
255         return et;
256     }
257 
operator ()mbgl::util::BuildBoundsMap258     BoundsMap operator()(const MultiPolygon<double>& polygons) const {
259         BoundsMap et;
260         for(const auto& polygon: polygons) {
261             buildPolygonTable(polygon, et);
262         }
263         return et;
264     }
265 
operator ()mbgl::util::BuildBoundsMap266     BoundsMap operator()(const mapbox::geometry::geometry_collection<double>&) const {
267         return {};
268     }
269 };
270 
Impl(int32_t z,const Geometry<double> & geom,bool project)271 TileCover::Impl::Impl(int32_t z, const Geometry<double>& geom, bool project)
272  : zoom(z) {
273     ToFeatureType toFeatureType;
274     isClosed = apply_visitor(toFeatureType, geom) == FeatureType::Polygon;
275 
276     BuildBoundsMap toBoundsMap(z, project);
277     boundsMap = apply_visitor(toBoundsMap, geom);
278     if (boundsMap.size() == 0) return;
279 
280     //Iniitalize the active edge table, and current row span
281     currentBounds = boundsMap.begin();
282     tileY = 0;
283     nextRow();
284     if (tileXSpans.empty()) return;
285     tileX = tileXSpans.front().first;
286 }
287 
288 // Aggregate all Bounds that start in or enter into the next tileY row. Multi-geoms
289 // may have discontinuity in the BoundMap, so skip forward to the next tileY row
290 // when the current/next row has no more bounds in it.
291 // Use scan_row to generate the tileX spans. Merge spans to avoid duplicate tiles
292 // in TileCoverImpl::next(). For closed geometry, use the non-zero rule to expand
293 // (fill) tiles between pairs of spans.
nextRow()294 void TileCover::Impl::nextRow() {
295     // Update activeBounds for next row
296     if (currentBounds != boundsMap.end()) {
297         if (activeBounds.size() == 0 && currentBounds->first > tileY) {
298             //For multi-geoms: use the next row with an edge table starting point
299             tileY = currentBounds->first;
300         }
301         if (tileY == currentBounds->first) {
302             std::move(currentBounds->second.begin(), currentBounds->second.end(),
303                 std::back_inserter(activeBounds));
304             currentBounds++;
305         }
306     }
307     //Scan the active bounds and update currentRange with x_min, x_max pairs
308     auto xps = util::scan_row(tileY, activeBounds);
309     if (xps.size() == 0) {
310         return;
311     }
312 
313     auto x_min = xps[0].xmin;
314     auto x_max = xps[0].xmax;
315     int32_t nzRule = xps[0].winding ? 1 : -1;
316     for (size_t i = 1; i < xps.size(); i++) {
317         auto xp = xps[i];
318         if (!(isClosed && nzRule != 0)) {
319             if (xp.xmin > x_max && xp.xmax >= x_max) {
320                 tileXSpans.emplace(x_min, x_max);
321                 x_min = xp.xmin;
322             }
323         }
324         nzRule += xp.winding ? 1 : -1;
325         x_max = std::max(x_min, xp.xmax);
326     }
327     tileXSpans.emplace(x_min, x_max);
328 }
329 
hasNext() const330 bool TileCover::Impl::hasNext() const {
331     return (!tileXSpans.empty()
332             && tileX < tileXSpans.front().second
333             && tileY < (1u << zoom));
334 }
335 
next()336 optional<UnwrappedTileID> TileCover::Impl::next() {
337     if (!hasNext()) return {};
338 
339     const auto x = tileX;
340     const auto y = tileY;
341     tileX++;
342     if (tileX >= tileXSpans.front().second) {
343         tileXSpans.pop();
344         if (tileXSpans.empty()) {
345             tileY++;
346             nextRow();
347         }
348         if (!tileXSpans.empty()) {
349             tileX = tileXSpans.front().first;
350         }
351     }
352     return UnwrappedTileID(zoom, x, y);
353 }
354 
355 } // namespace util
356 } // namespace mbgl
357