From 4344812ef5ef3b5b67e9f313e0c3cc87b3e95c46 Mon Sep 17 00:00:00 2001 From: Asheem Mamoowala Date: Thu, 26 Apr 2018 07:43:38 -0700 Subject: [core] Streaming TileCover for polygonal regions (#11267) A per-tile streaming algorithm for tile cover on points, lines, and polygons. Works for individual zoom levels, and not zoom ranges. --- benchmark/util/tilecover.benchmark.cpp | 96 ++++++ cmake/benchmark-files.cmake | 1 + cmake/core-files.cmake | 2 + include/mbgl/util/projection.hpp | 16 +- platform/default/mbgl/storage/offline.cpp | 2 +- .../renderer/layers/render_background_layer.cpp | 1 + src/mbgl/util/tile_cover.cpp | 104 ++++-- src/mbgl/util/tile_cover.hpp | 24 +- src/mbgl/util/tile_cover_impl.cpp | 365 +++++++++++++++++++++ src/mbgl/util/tile_cover_impl.hpp | 90 +++++ test/util/tile_cover.test.cpp | 232 ++++++++++++- test/util/tile_range.test.cpp | 8 +- 12 files changed, 886 insertions(+), 55 deletions(-) create mode 100644 benchmark/util/tilecover.benchmark.cpp create mode 100644 src/mbgl/util/tile_cover_impl.cpp create mode 100644 src/mbgl/util/tile_cover_impl.hpp diff --git a/benchmark/util/tilecover.benchmark.cpp b/benchmark/util/tilecover.benchmark.cpp new file mode 100644 index 0000000000..186de6f216 --- /dev/null +++ b/benchmark/util/tilecover.benchmark.cpp @@ -0,0 +1,96 @@ +#include +#include +#include +#include +#include +#include + +using namespace mbgl; + +static const LatLngBounds sanFrancisco = + LatLngBounds::hull({ 37.6609, -122.5744 }, { 37.8271, -122.3204 }); + +static void TileCountBounds(benchmark::State& state) { + std::size_t length = 0; + while (state.KeepRunning()) { + auto count = util::tileCount(sanFrancisco, 10); + length += count; + } +} + +static void TileCoverPitchedViewport(benchmark::State& state) { + Transform transform; + transform.resize({ 512, 512 }); + // slightly offset center so that tile order is better defined + transform.setLatLng({ 0.1, -0.1 }); + transform.setZoom(8); + transform.setAngle(5.0); + transform.setPitch(40.0 * M_PI / 180.0); + + std::size_t length = 0; + while (state.KeepRunning()) { + auto tiles = util::tileCover(transform.getState(), 8); + length += tiles.size(); + } +} + +static void TileCoverBounds(benchmark::State& state) { + std::size_t length = 0; + while (state.KeepRunning()) { + auto tiles = util::tileCover(sanFrancisco, 8); + length += tiles.size(); + } +} + +static const auto geomPolygon = Polygon{ + { + {-122.5143814086914,37.779127216982424}, + {-122.50811576843262,37.72721239056709}, + {-122.50313758850099,37.70820178063929}, + {-122.3938751220703,37.707454835665274}, + {-122.37567901611328,37.70663997801684}, + {-122.36297607421874,37.71343018466285}, + {-122.354736328125,37.727280276860036}, + {-122.36469268798828,37.73868429065797}, + {-122.38014221191408,37.75442980295571}, + {-122.38391876220702,37.78753873820529}, + {-122.35919952392578,37.8065289741725}, + {-122.35679626464844,37.820632846207864}, + {-122.3712158203125,37.835276322922695}, + {-122.3818588256836,37.82958198283902}, + {-122.37190246582031,37.80788523279169}, + {-122.38735198974608,37.791337175930686}, + {-122.40966796874999,37.812767557570204}, + {-122.46425628662108,37.807071480609274}, + {-122.46803283691405,37.810326435534755}, + {-122.47901916503906,37.81168262440736}, + {-122.48966217041016,37.78916666399649}, + {-122.50579833984375,37.78781006166096}, + {-122.5143814086914,37.779127216982424} + } +}; + +static void TileCoverPolygon(benchmark::State& state) { + std::size_t length = 0; + + while (state.KeepRunning()) { + auto tiles = util::tileCover(geomPolygon, 8); + length += tiles.size(); + } +} + +static void TileCountPolygon(benchmark::State& state) { + std::size_t length = 0; + + while (state.KeepRunning()) { + auto tiles = util::tileCount(geomPolygon, 16); + length += tiles; + } +} + +BENCHMARK(TileCountBounds); +BENCHMARK(TileCountPolygon); +BENCHMARK(TileCoverPitchedViewport); +BENCHMARK(TileCoverBounds); +BENCHMARK(TileCoverPolygon); + diff --git a/cmake/benchmark-files.cmake b/cmake/benchmark-files.cmake index fdafcf30f9..21547eb9c4 100644 --- a/cmake/benchmark-files.cmake +++ b/cmake/benchmark-files.cmake @@ -23,5 +23,6 @@ set(MBGL_BENCHMARK_FILES # util benchmark/util/dtoa.benchmark.cpp + benchmark/util/tilecover.benchmark.cpp ) diff --git a/cmake/core-files.cmake b/cmake/core-files.cmake index 0853097630..d4fec34bf5 100644 --- a/cmake/core-files.cmake +++ b/cmake/core-files.cmake @@ -753,6 +753,8 @@ set(MBGL_CORE_FILES src/mbgl/util/tile_coordinate.hpp src/mbgl/util/tile_cover.cpp src/mbgl/util/tile_cover.hpp + src/mbgl/util/tile_cover_impl.cpp + src/mbgl/util/tile_cover_impl.hpp src/mbgl/util/tile_range.hpp src/mbgl/util/tiny_sdf.cpp src/mbgl/util/tiny_sdf.hpp diff --git a/include/mbgl/util/projection.hpp b/include/mbgl/util/projection.hpp index b4a34521a4..65a84d8dc2 100644 --- a/include/mbgl/util/projection.hpp +++ b/include/mbgl/util/projection.hpp @@ -78,8 +78,9 @@ public: return project_(latLng, worldSize(scale)); } - static Point project(const LatLng& latLng, uint8_t zoom) { - return project_(latLng, std::pow(2.0, zoom)); + //Returns point on tile + static Point project(const LatLng& latLng, int32_t zoom) { + return project_(latLng, 1 << zoom); } static LatLng unproject(const Point& p, double scale, LatLng::WrapMode wrapMode = LatLng::Unwrapped) { @@ -90,16 +91,7 @@ public: wrapMode }; } - - // Project lat, lon to point in a zoom-dependent world size - static Point project(const LatLng& point, uint8_t zoom, uint16_t tileSize) { - const double t2z = tileSize * std::pow(2, zoom); - Point pt = project_(point, t2z); - // Flip y coordinate - auto x = ::round(std::min(pt.x, t2z)); - auto y = ::round(std::min(t2z - pt.y, t2z)); - return { x, y }; - } + private: static Point project_(const LatLng& latLng, double worldSize) { return Point { diff --git a/platform/default/mbgl/storage/offline.cpp b/platform/default/mbgl/storage/offline.cpp index 7670790be9..598a0b182b 100644 --- a/platform/default/mbgl/storage/offline.cpp +++ b/platform/default/mbgl/storage/offline.cpp @@ -43,7 +43,7 @@ uint64_t OfflineTilePyramidRegionDefinition::tileCount(style::SourceType type, u const Range clampedZoomRange = coveringZoomRange(type, tileSize, zoomRange); unsigned long result = 0;; for (uint8_t z = clampedZoomRange.min; z <= clampedZoomRange.max; z++) { - result += util::tileCount(bounds, z, tileSize); + result += util::tileCount(bounds, z); } return result; diff --git a/src/mbgl/renderer/layers/render_background_layer.cpp b/src/mbgl/renderer/layers/render_background_layer.cpp index aebc4cc9aa..44c3fffb6c 100644 --- a/src/mbgl/renderer/layers/render_background_layer.cpp +++ b/src/mbgl/renderer/layers/render_background_layer.cpp @@ -7,6 +7,7 @@ #include #include #include +#include namespace mbgl { diff --git a/src/mbgl/util/tile_cover.cpp b/src/mbgl/util/tile_cover.cpp index 39b562d811..488e6b88ce 100644 --- a/src/mbgl/util/tile_cover.cpp +++ b/src/mbgl/util/tile_cover.cpp @@ -2,13 +2,18 @@ #include #include #include +#include +#include #include +#include namespace mbgl { namespace { +using ScanLine = const std::function; + // Taken from polymaps src/Layer.js // https://github.com/simplegeo/polymaps/blob/master/src/Layer.js#L333-L383 struct edge { @@ -27,8 +32,6 @@ struct edge { } }; -using ScanLine = const std::function; - // scan-line conversion static void scanSpans(edge e0, edge e1, int32_t ymin, int32_t ymax, ScanLine scanLine) { double y0 = ::fmax(ymin, std::floor(e1.y0)); @@ -147,11 +150,11 @@ std::vector tileCover(const LatLngBounds& bounds_, int32_t z) { { std::min(bounds_.north(), util::LATITUDE_MAX), bounds_.east() }); return tileCover( - TileCoordinate::fromLatLng(z, bounds.northwest()).p, - TileCoordinate::fromLatLng(z, bounds.northeast()).p, - TileCoordinate::fromLatLng(z, bounds.southeast()).p, - TileCoordinate::fromLatLng(z, bounds.southwest()).p, - TileCoordinate::fromLatLng(z, bounds.center()).p, + Projection::project(bounds.northwest(), z), + Projection::project(bounds.northeast(), z), + Projection::project(bounds.southeast(), z), + Projection::project(bounds.southwest(), z), + Projection::project(bounds.center(), z), z); } @@ -169,25 +172,80 @@ std::vector tileCover(const TransformState& state, int32_t z) { z); } +std::vector tileCover(const Geometry& geometry, int32_t z) { + std::vector result; + TileCover tc(geometry, z, true); + while (tc.hasNext()) { + result.push_back(*tc.next()); + }; + + return result; +} + // Taken from https://github.com/mapbox/sphericalmercator#xyzbbox-zoom-tms_style-srs // Computes the projected tiles for the lower left and upper right points of the bounds // and uses that to compute the tile cover count -uint64_t tileCount(const LatLngBounds& bounds, uint8_t zoom, uint16_t tileSize_){ - - auto sw = Projection::project(bounds.southwest().wrapped(), zoom, tileSize_); - auto ne = Projection::project(bounds.northeast().wrapped(), zoom, tileSize_); - - auto x1 = floor(sw.x/ tileSize_); - auto x2 = floor((ne.x - 1) / tileSize_); - auto y1 = floor(sw.y/ tileSize_); - auto y2 = floor((ne.y - 1) / tileSize_); - - auto minX = ::fmax(std::min(x1, x2), 0); - auto maxX = std::max(x1, x2); - auto minY = (std::pow(2, zoom) - 1) - std::max(y1, y2); - auto maxY = (std::pow(2, zoom) - 1) - ::fmax(std::min(y1, y2), 0); - - return (maxX - minX + 1) * (maxY - minY + 1); +uint64_t tileCount(const LatLngBounds& bounds, uint8_t zoom){ + if (zoom == 0) { + return 1; + } + auto sw = Projection::project(bounds.southwest(), zoom); + auto ne = Projection::project(bounds.northeast(), zoom); + auto maxTile = std::pow(2.0, zoom); + auto x1 = floor(sw.x); + auto x2 = ceil(ne.x) - 1; + auto y1 = util::clamp(floor(sw.y), 0.0, maxTile - 1); + auto y2 = util::clamp(floor(ne.y), 0.0, maxTile - 1); + + auto dx = x1 > x2 ? (maxTile - x1) + x2 : x2 - x1; + auto dy = y1 - y2; + return (dx + 1) * (dy + 1); +} + +uint64_t tileCount(const Geometry& geometry, uint8_t z) { + uint64_t tileCount = 0; + + TileCover tc(geometry, z, true); + while (tc.next()) { + tileCount++; + }; + return tileCount; +} + +TileCover::TileCover(const LatLngBounds&bounds_, int32_t z) { + LatLngBounds bounds = LatLngBounds::hull( + { std::max(bounds_.south(), -util::LATITUDE_MAX), bounds_.west() }, + { std::min(bounds_.north(), util::LATITUDE_MAX), bounds_.east() }); + + if (bounds.isEmpty() || + bounds.south() > util::LATITUDE_MAX || + bounds.north() < -util::LATITUDE_MAX) { + bounds = LatLngBounds::world(); + } + + auto sw = Projection::project(bounds.southwest(), z); + auto ne = Projection::project(bounds.northeast(), z); + auto se = Projection::project(bounds.southeast(), z); + auto nw = Projection::project(bounds.northwest(), z); + + Polygon p({ {sw, nw, ne, se, sw} }); + impl = std::make_unique(z, p, false); +} + +TileCover::TileCover(const Geometry& geom, int32_t z, bool project/* = true*/) + : impl( std::make_unique(z, geom, project)) { +} + +TileCover::~TileCover() { + +} + +optional TileCover::next() { + return impl->next(); +} + +bool TileCover::hasNext() { + return impl->hasNext(); } } // namespace util diff --git a/src/mbgl/util/tile_cover.hpp b/src/mbgl/util/tile_cover.hpp index b2098b59b8..c953d764d2 100644 --- a/src/mbgl/util/tile_cover.hpp +++ b/src/mbgl/util/tile_cover.hpp @@ -2,9 +2,11 @@ #include #include -#include +#include +#include #include +#include namespace mbgl { @@ -13,13 +15,31 @@ class LatLngBounds; namespace util { +// Helper class to stream tile-cover results per row +class TileCover { +public: + TileCover(const LatLngBounds&, int32_t z); + // When project == true, projects the geometry points to tile coordinates + TileCover(const Geometry&, int32_t z, bool project = true); + ~TileCover(); + + optional next(); + bool hasNext(); + +private: + class Impl; + std::unique_ptr impl; +}; + int32_t coveringZoomLevel(double z, style::SourceType type, uint16_t tileSize); std::vector tileCover(const TransformState&, int32_t z); std::vector tileCover(const LatLngBounds&, int32_t z); +std::vector tileCover(const Geometry&, int32_t z); // Compute only the count of tiles needed for tileCover -uint64_t tileCount(const LatLngBounds&, uint8_t z, uint16_t tileSize); +uint64_t tileCount(const LatLngBounds&, uint8_t z); +uint64_t tileCount(const Geometry&, uint8_t z); } // namespace util } // namespace mbgl diff --git a/src/mbgl/util/tile_cover_impl.cpp b/src/mbgl/util/tile_cover_impl.cpp new file mode 100644 index 0000000000..b3fc07f7dd --- /dev/null +++ b/src/mbgl/util/tile_cover_impl.cpp @@ -0,0 +1,365 @@ +#include +#include + +#include +#include +#include +#include +#include + +namespace mbgl { +namespace util { + +using PointList = std::vector>; + +struct TileSpan { + int32_t xmin, xmax; + bool winding; +}; + + +// Find the first local minimum going forward in the list. +void start_list_on_local_minimum(PointList& points) { + auto prev_pt = std::prev(points.end(), 2); + auto pt = points.begin(); + auto next_pt = std::next(pt); + while (pt != points.end()) { + if ((pt->y <= prev_pt->y) && + (pt->y < next_pt->y)) { + break; + } + prev_pt = pt; + pt++; + next_pt++; + if (next_pt == points.end()) { next_pt = std::next(points.begin()); } + } + //Re-close linear rings with first_pt = last_pt + if (points.back() == points.front()) { + points.pop_back(); + } + std::rotate(points.begin(), pt, points.end()); + points.push_back(*points.begin()); +} + +//Create a bound towards a local maximum point, starting from pt. +Bound create_bound_towards_maximum(PointList& points, PointList::iterator& pt) { + if (std::distance(pt, points.end()) < 2) { return {}; } + if (std::distance(pt, points.end()) == 2) { + Bound bnd; + if (pt->y < std::next(pt)->y) { + std::copy(pt, points.end(), std::back_inserter(bnd.points)); + bnd.winding = true; + } + else { + std::reverse_copy(pt, points.end(), std::back_inserter(bnd.points)); + bnd.winding = false; + } + pt = points.end(); + return bnd; + } + const auto begin = pt; + auto prev_pt = pt == points.begin() ? std::prev(points.end(), 2) : std::prev(pt); + auto next_pt = std::next(pt) == points.end() ? std::next(points.begin()) : std::next(pt); + while (pt != points.end()) { + if ((pt->y >= prev_pt->y) && + (pt->y > next_pt->y )) { + break; + } + prev_pt = pt; + pt++; + next_pt++; + if (next_pt == points.end()) { next_pt = std::next(points.begin()); } + } + + Bound bnd; + if (std::next(pt) == points.end()) { next_pt = points.end(); pt++; }; + bnd.points.reserve(static_cast(std::distance(begin, next_pt))); + std::copy(begin, next_pt, std::back_inserter(bnd.points)); + bnd.winding = true; + return bnd; +} + +//Create a bound towards a local minimum point, starting from pt. +Bound create_bound_towards_minimum(PointList& points, PointList::iterator& pt) { + if (std::distance(pt, points.end()) < 2) { return {}; } + if (std::distance(pt, points.end()) == 2) { + Bound bnd; + if (pt->y < std::next(pt)->y) { + std::copy(pt, points.end(), std::back_inserter(bnd.points)); + bnd.winding = true; + } + else { + std::reverse_copy(pt, points.end(), std::back_inserter(bnd.points)); + bnd.winding = false; + } + pt = points.end(); + return bnd; + } + auto begin = pt; + auto prev_pt = pt == points.begin() ? std::prev(points.end(), 2) : std::prev(pt); + auto next_pt = std::next(pt) == points.end() ? std::next(points.begin()) : std::next(pt); + while (pt != points.end()) { + if ((pt->y <= prev_pt->y) && + (pt->y < next_pt->y)) { + break; + } + prev_pt = pt; + pt++; + next_pt++; + if (next_pt == points.end()) { next_pt = std::next(points.begin()); } + } + + Bound bnd; + if (std::next(pt) == points.end()) { next_pt = points.end(); pt++; }; + bnd.points.reserve(static_cast(std::distance(begin, next_pt))); + //For bounds that start at a max, reverse copy so that all bounds start at a min + std::reverse_copy(begin, next_pt, std::back_inserter(bnd.points)); + bnd.winding = false; + return bnd; +} + +//Build a map of bounds and their starting Y tile coordinate. +void build_bounds_map(PointList& points, uint32_t maxTile, BoundsMap& et, bool closed = false) { + if (points.size() < 2) return; + //While traversing closed rings, start the bounds at a local minimum + if (closed) { + start_list_on_local_minimum(points); + } + + auto pointsIter = points.begin(); + while (pointsIter != points.end()) { + Bound to_max = create_bound_towards_maximum(points, pointsIter); + Bound to_min = create_bound_towards_minimum(points, pointsIter); + + if (to_max.points.size() > 0) { + // Projections may result in values beyond the bounds, clamp to max tile coordinates + const auto y = static_cast(std::floor(clamp(to_max.points.front().y, 0.0, (double)maxTile))); + et[y].push_back(to_max); + } + if (to_min.points.size() > 0) { + const auto y = static_cast(std::floor(clamp(to_min.points.front().y, 0.0, (double)maxTile))); + et[y].push_back(to_min); + } + } + assert(pointsIter == points.end()); +} + +void update_span(TileSpan& xp, double x) { + xp.xmin = std::min(xp.xmin, static_cast(std::floor(x))); + xp.xmax = std::max(xp.xmax, static_cast(std::ceil(x))); +} + +//Build a vector of X tile-coordinates spanned by each bound. +std::vector scan_row(uint32_t y, Bounds& aet) { + std::vector tile_range; + tile_range.reserve(aet.size()); + + for(Bound& b: aet) { + TileSpan xp = { INT_MAX, 0, b.winding }; + double x; + const auto numEdges = b.points.size() - 1; + assert(numEdges >= 1); + while (b.currentPoint < numEdges) { + x = b.interpolate(y); + update_span(xp, x); + + // If this edge ends beyond the current row, find the x-intercept where + // it exits the row + auto& p1 = b.points[b.currentPoint + 1]; + if (p1.y > y+1) { + x = b.interpolate(y+1); + update_span(xp, x); + break; + } else if(b.currentPoint == numEdges - 1) { + // For last edge, consider x-intercept at the end of the edge. + x = p1.x; + update_span(xp, x); + } + b.currentPoint++; + } + tile_range.push_back(xp); + } + // Erase bounds in the active table whose current edge ends inside this row, + // or there are no more edges + auto bound = aet.begin(); + while (bound != aet.end()) { + if ( bound->currentPoint == bound->points.size() - 1 && + bound->points[bound->currentPoint].y <= y+1) { + bound = aet.erase(bound); + } else { + bound++; + } + } + // Sort the X-extents of each crossing bound by x_min, x_max + std::sort(tile_range.begin(), tile_range.end(), [] (TileSpan& a, TileSpan& b) { + return std::tie(a.xmin, a.xmax) < std::tie(b.xmin, b.xmax); + }); + + return tile_range; +} + +struct BuildBoundsMap { + int32_t zoom; + bool project = false; + BuildBoundsMap(int32_t z, bool p): zoom(z), project(p) {} + + void buildTable(const std::vector>& points, BoundsMap& et, bool closed = false) const { + PointList projectedPoints; + if (project) { + projectedPoints.reserve(points.size()); + for(const auto&p : points) { + projectedPoints.push_back( + Projection::project(LatLng{ p.y, p.x }, zoom)); + } + } else { + projectedPoints.insert(projectedPoints.end(), points.begin(), points.end()); + } + build_bounds_map(projectedPoints, 1 << zoom, et, closed); + } + + void buildPolygonTable(const Polygon& polygon, BoundsMap& et) const { + for(const auto&ring : polygon) { + buildTable(ring, et, true); + } + } + BoundsMap operator()(const Point&p) const { + Bound bnd; + auto point = p; + if(project) { + point = Projection::project(LatLng{p.y, p.x}, zoom); + } + bnd.points.insert(bnd.points.end(), 2, point); + bnd.winding = false; + BoundsMap et; + const auto y = static_cast(std::floor(clamp(point.y, 0.0, (double)(1 << zoom)))); + et[y].push_back(bnd); + return et; + } + + BoundsMap operator()(const MultiPoint& points) const { + BoundsMap et; + for (const Point& p: points) { + Bound bnd; + auto point = p; + if(project) { + point = Projection::project(LatLng{p.y, p.x}, zoom); + } + bnd.points.insert(bnd.points.end(), 2, point); + bnd.winding = false; + const auto y = static_cast(std::floor(clamp(point.y, 0.0, (double)(1 << zoom)))); + et[y].push_back(bnd); + } + return et; + } + + BoundsMap operator()(const LineString& lines) const { + BoundsMap et; + buildTable(lines, et); + return et; + } + + BoundsMap operator()(const MultiLineString& lines) const { + BoundsMap et; + for(const auto&line : lines) { + buildTable(line, et); + } + return et; + } + + BoundsMap operator()(const Polygon& polygon) const { + BoundsMap et; + buildPolygonTable(polygon, et); + return et; + } + + BoundsMap operator()(const MultiPolygon& polygons) const { + BoundsMap et; + for(const auto& polygon: polygons) { + buildPolygonTable(polygon, et); + } + return et; + } + + BoundsMap operator()(const mapbox::geometry::geometry_collection&) const { + return {}; + } +}; + +TileCover::Impl::Impl(int32_t z, const Geometry& geom, bool project) + : zoom(z) { + ToFeatureType toFeatureType; + isClosed = apply_visitor(toFeatureType, geom) == FeatureType::Polygon; + + BuildBoundsMap toBoundsMap(z, project); + boundsMap = apply_visitor(toBoundsMap, geom); + if (boundsMap.size() == 0) return; + + //Iniitalize the active edge table, and current row span + currentBounds = boundsMap.begin(); + tileY = 0; + nextRow(); + if (tileXSpans.empty()) return; + tileX = tileXSpans.front().first; +} + +void TileCover::Impl::nextRow() { + // Update AET for next row + if (currentBounds != boundsMap.end()) { + if (activeBounds.size() == 0 && currentBounds->first > tileY) { + //For multi-geoms: use the next row with an edge table starting point + tileY = currentBounds->first; + } + if (tileY == currentBounds->first) { + + std::move(currentBounds->second.begin(), currentBounds->second.end(), std::back_inserter(activeBounds)); + currentBounds++; + } + } + //Scan aet and update currenRange with x_min, x_max pairs + auto xps = util::scan_row(tileY, activeBounds); + if (xps.size() == 0) { + return; + } + + auto x_min = xps[0].xmin; + auto x_max = xps[0].xmax; + int32_t nzRule = xps[0].winding ? 1 : -1; + for (size_t i = 1; i < xps.size(); i++) { + auto xp = xps[i]; + if (!(isClosed && nzRule != 0)) { + if (xp.xmin > x_max && xp.xmax >= x_max) { + tileXSpans.emplace(x_min, x_max); + x_min = xp.xmin; + } + } + nzRule += xp.winding ? 1 : -1; + x_max = std::max(x_min, xp.xmax); + } + tileXSpans.emplace(x_min, x_max); +} + +bool TileCover::Impl::hasNext() const { + return (!tileXSpans.empty() && tileX < tileXSpans.front().second && tileY < (1u << zoom)); +} + +optional TileCover::Impl::next() { + if (!hasNext()) return {}; + + const auto x = tileX; + const auto y = tileY; + tileX++; + if (tileX >= tileXSpans.front().second) { + tileXSpans.pop(); + if(tileXSpans.empty()) { + tileY++; + nextRow(); + } + if (!tileXSpans.empty()) { + tileX = tileXSpans.front().first; + } + } + return UnwrappedTileID(zoom, x, y); +} + +} // namespace util +} // namespace mbgl diff --git a/src/mbgl/util/tile_cover_impl.hpp b/src/mbgl/util/tile_cover_impl.hpp new file mode 100644 index 0000000000..7c16718984 --- /dev/null +++ b/src/mbgl/util/tile_cover_impl.hpp @@ -0,0 +1,90 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include + +namespace mbgl { + +class TransformState; +class LatLngBounds; + +namespace util { + +struct Bound; + +using Bounds = std::vector; +using BoundsMap = std::map; + +// A chain of points from a local minimum to a local maximum. `winding` indicates +// the direction of the original geometry. +struct Bound { + std::vector> points; + size_t currentPoint = 0; + bool winding = false; + + Bound() = default; + Bound(const Bound& rhs) { + points = rhs.points; + currentPoint = rhs.currentPoint; + winding = rhs.winding; + } + Bound& operator=(Bound&& rhs) { + points = std::move(rhs.points); + currentPoint = rhs.currentPoint; + winding = rhs.winding; + return *this; + } + + // Compute the interpolated x coordinate at y for the current edge + double interpolate(uint32_t y) { + const auto& p0 = points[currentPoint]; + const auto& p1 = points[currentPoint + 1]; + + const auto dx = p1.x - p0.x; + const auto dy = p1.y - p0.y; + auto x = p0.x; + if (dx == 0) { + return x; + } else if (dy == 0){ + return y <= p0.y ? p0.x : p1.x; + } + if (y < p0.y) return x; + if (y > p1.y) return p1.x; + x = (dx / dy) * (y - p0.y) + p0.x; + return x; + } +}; + +class TileCover::Impl { +public: + Impl(int32_t z, const Geometry& geom, bool project = true); + ~Impl() = default; + + optional next(); + bool hasNext() const; + +private: + using TileSpans = std::queue>; + + void nextRow(); + + const int32_t zoom; + bool isClosed; + + BoundsMap boundsMap; + BoundsMap::iterator currentBounds; + // List of bounds that begin at or before `tileY` + Bounds activeBounds; + + TileSpans tileXSpans; + uint32_t tileY; + int32_t tileX; +}; + +} // namespace util +} // namespace mbgl diff --git a/test/util/tile_cover.test.cpp b/test/util/tile_cover.test.cpp index 933c18b5ea..7defa761af 100644 --- a/test/util/tile_cover.test.cpp +++ b/test/util/tile_cover.test.cpp @@ -2,6 +2,8 @@ #include #include +#include + #include using namespace mbgl; @@ -22,8 +24,8 @@ TEST(TileCover, Antarctic) { TEST(TileCover, WorldZ0) { EXPECT_EQ((std::vector{ - { 0, 0, 0 }, - }), + { 0, 0, 0 }, + }), util::tileCover(LatLngBounds::world(), 0)); } @@ -37,15 +39,15 @@ TEST(TileCover, Pitch) { transform.setPitch(40.0 * M_PI / 180.0); EXPECT_EQ((std::vector{ - { 2, 1, 2 }, { 2, 1, 1 }, { 2, 2, 2 }, { 2, 2, 1 }, { 2, 3, 2 } - }), + { 2, 1, 2 }, { 2, 1, 1 }, { 2, 2, 2 }, { 2, 2, 1 }, { 2, 3, 2 } + }), util::tileCover(transform.getState(), 2)); } TEST(TileCover, WorldZ1) { EXPECT_EQ((std::vector{ - { 1, 0, 0 }, { 1, 0, 1 }, { 1, 1, 0 }, { 1, 1, 1 }, - }), + { 1, 0, 0 }, { 1, 0, 1 }, { 1, 1, 0 }, { 1, 1, 1 }, + }), util::tileCover(LatLngBounds::world(), 1)); } @@ -59,21 +61,44 @@ TEST(TileCover, SingletonZ1) { util::tileCover(LatLngBounds::singleton({ 0, 0 }), 1)); } +TEST(TileCoverStream, Arctic) { + auto bounds = LatLngBounds::hull({ 84, -180 }, { 70, 180 }); + auto zoom = 3; + util::TileCover tc(bounds, zoom); + auto results = util::tileCover(bounds, zoom); + std::vector t; + while(tc.hasNext()) { + t.push_back(*tc.next()); + }; + EXPECT_EQ(t.size(), results.size()); +} + +TEST(TileCoverStream, WorldZ1) { + auto zoom = 1; + util::TileCover tc(LatLngBounds::world(), zoom); + std::vector t; + while(tc.hasNext()) { + t.push_back(*tc.next()); + }; + EXPECT_EQ((std::vector{ + { 1, 0, 0 }, { 1, 1, 0 }, { 1, 0, 1 }, { 1, 1, 1 }, + }), t); +} + static const LatLngBounds sanFrancisco = LatLngBounds::hull({ 37.6609, -122.5744 }, { 37.8271, -122.3204 }); TEST(TileCover, SanFranciscoZ0) { EXPECT_EQ((std::vector{ - { 0, 0, 0 }, - }), + { 0, 0, 0 }, + }), util::tileCover(sanFrancisco, 0)); } TEST(TileCover, SanFranciscoZ10) { EXPECT_EQ((std::vector{ - { 10, 163, 395 }, { 10, 163, 396 }, { 10, 164, 395 }, { 10, 164, 396 }, - - }), + { 10, 163, 395 }, { 10, 163, 396 }, { 10, 164, 395 }, { 10, 164, 396 }, + }), util::tileCover(sanFrancisco, 10)); } @@ -85,11 +110,192 @@ TEST(TileCover, SanFranciscoZ0Wrapped) { util::tileCover(sanFranciscoWrapped, 0)); } +TEST(TileCover, GeomPoint) { + auto point = Point{ -122.5744, 37.6609 }; + + EXPECT_EQ((std::vector{ {2 ,0 ,1 } }), + util::tileCover(point, 2)); +} + +TEST(TileCover, GeomMultiPoint) { + auto points = MultiPoint{ { -122.5, 37.76 }, { -122.4, 37.76} }; + + EXPECT_EQ((std::vector{ + {20, 167480, 405351}, {20, 167772, 405351} }), + util::tileCover(points, 20)); +} + +TEST(TileCover, GeomLineZ10) { + auto lineCover = util::tileCover(LineString{ + {-121.49368286132812,38.57903714667459}, + {-122.4422836303711,37.773157169570695} + }, 10); + EXPECT_EQ((std::vector{ + { 10, 166, 392}, {10, 165, 393}, {10, 166, 393}, + {10, 164, 394}, {10, 165, 394}, {10,163,395}, {10, 164, 395} + }),lineCover); + +} + +TEST(TileCover, WrappedGeomLineZ10) { + auto lineCover = util::tileCover(LineString{ + {-179.93342914581299,38.892101707724315}, + {-180.02394485473633,38.89203490311832} + }, 10); + EXPECT_EQ((std::vector{ { 10, -1, 391 }, { 10, 0, 391 } }), + lineCover); + + lineCover = util::tileCover(LineString{ + {179.93342914581299,38.892101707724315}, + {180.02394485473633,38.89203490311832} + }, 10); + EXPECT_EQ((std::vector{ { 10, 1023, 391 }, { 10, 1024, 391 } }), + lineCover); +} + +TEST(TileCover, GeomMultiLineString) { + auto geom = MultiLineString{ + { { -122.5, 37.76 }, { -122.4, 37.76} }, + { { -122.5, 37.72 }, { -122.4, 37.72} } }; + + EXPECT_EQ((std::vector{ + {14, 2616, 6333}, {14, 2617, 6333}, {14, 2618, 6333}, + {14, 2619, 6333}, {14, 2620, 6333}, {14, 2621, 6333}, + {14, 2616, 6335}, {14, 2617, 6335}, {14, 2618, 6335}, + {14, 2619, 6335}, {14, 2620, 6335}, {14, 2621, 6335}}), + util::tileCover(geom, 14)); +} + +TEST(TileCover, GeomPolygon) { + auto polygon = Polygon{ + { + {5.09765625,53.067626642387374}, + {2.373046875,43.389081939117496}, + {-4.74609375,48.45835188280866}, + {-1.494140625,37.09023980307208}, + {22.587890625,36.24427318493909}, + {31.640625,46.13417004624326}, + {17.841796875,54.7246201949245}, + {5.09765625,53.067626642387374}, + },{ + {19.6875,49.66762782262194}, + {22.8515625,43.51668853502906}, + {13.623046875,45.089035564831036}, + {16.34765625,39.095962936305476}, + {5.185546875,41.244772343082076}, + {8.701171874999998,50.233151832472245}, + {19.6875,49.66762782262194} + } + }; + + auto results = util::tileCover(polygon, 8); + + EXPECT_NE(std::find(results.begin(), results.end(), UnwrappedTileID{8, 134, 87}), results.end()); + EXPECT_NE(std::find(results.begin(), results.end(), UnwrappedTileID{8, 139, 87}), results.end()); + // Should have a hole + EXPECT_EQ(std::find(results.begin(), results.end(), UnwrappedTileID{8, 136, 87}), results.end()); +} + +TEST(TileCover, GeomMultiPolygon) { + auto multiPolygon = MultiPolygon{ + {{ + {5.09765625,53.067626642387374}, + {2.373046875,43.389081939117496}, + {-4.74609375,48.45835188280866}, + {-1.494140625,37.09023980307208}, + {22.587890625,36.24427318493909}, + {31.640625,46.13417004624326}, + {17.841796875,54.7246201949245}, + {5.09765625,53.067626642387374}, + }},{{ + {59.150390625,45.460130637921004}, + {65.126953125,41.11246878918088}, + {69.169921875,47.45780853075031}, + {63.896484375,50.064191736659104}, + {59.150390625,45.460130637921004} + }} + }; + auto results = util::tileCover(multiPolygon, 8); + + EXPECT_EQ(424u, results.size()); + EXPECT_NE(std::find(results.begin(), results.end(), UnwrappedTileID{8, 139, 87}), results.end()); + EXPECT_NE(std::find(results.begin(), results.end(), UnwrappedTileID{8, 136, 87}), results.end()); + EXPECT_NE(std::find(results.begin(), results.end(), UnwrappedTileID{8, 174, 94}), results.end()); +} + +TEST(TileCover, GeomSanFranciscoPoly) { + auto sanFranciscoGeom = Polygon{ + { + {-122.5143814086914,37.779127216982424}, + {-122.50811576843262,37.72721239056709}, + {-122.50313758850099,37.70820178063929}, + {-122.3938751220703,37.707454835665274}, + {-122.37567901611328,37.70663997801684}, + {-122.36297607421874,37.71343018466285}, + {-122.354736328125,37.727280276860036}, + {-122.36469268798828,37.73868429065797}, + {-122.38014221191408,37.75442980295571}, + {-122.38391876220702,37.78753873820529}, + {-122.35919952392578,37.8065289741725}, + {-122.35679626464844,37.820632846207864}, + {-122.3712158203125,37.835276322922695}, + {-122.3818588256836,37.82958198283902}, + {-122.37190246582031,37.80788523279169}, + {-122.38735198974608,37.791337175930686}, + {-122.40966796874999,37.812767557570204}, + {-122.46425628662108,37.807071480609274}, + {-122.46803283691405,37.810326435534755}, + {-122.47901916503906,37.81168262440736}, + {-122.48966217041016,37.78916666399649}, + {-122.50579833984375,37.78781006166096}, + {-122.5143814086914,37.779127216982424} + } + }; + + EXPECT_EQ((std::vector{ + { 12, 654, 1582 }, { 12, 655, 1582 }, + { 12, 654, 1583 }, { 12, 655, 1583 }, + { 12, 654, 1584 }, { 12, 655, 1584 } + }), util::tileCover(sanFranciscoGeom, 12)); +} + +TEST(TileCover, GeomInvalid) { + auto point = Point{ -122.5744, 97.6609 }; + EXPECT_THROW(util::tileCover(point, 2), std::domain_error); + + auto badPoly = Polygon { { {1.0, 35.0} } }; + EXPECT_EQ((std::vector{ }), util::tileCover(badPoly, 16)); + + //Should handle open polygons. + badPoly = Polygon { { {1.0, 34.2}, {1.0, 34.4}, {0.5, 34.3} } }; + EXPECT_EQ((std::vector{ + { 10, 513, 407 }, { 10, 514, 407}, + { 10, 513, 408 }, { 10, 514, 408} + }), util::tileCover(badPoly, 10)); +} + + +TEST(TileCount, World) { + EXPECT_EQ(1u, util::tileCount(LatLngBounds::world(), 0)); + EXPECT_EQ(4u, util::tileCount(LatLngBounds::world(), 1)); +} + TEST(TileCount, SanFranciscoZ10) { - EXPECT_EQ(4u, util::tileCount(sanFrancisco, 10, util::tileSize)); + EXPECT_EQ(4u, util::tileCount(sanFrancisco, 10)); +} + +TEST(TileCount, SanFranciscoWrappedZ10) { + EXPECT_EQ(4u, util::tileCount(sanFranciscoWrapped, 10)); } TEST(TileCount, SanFranciscoZ22) { - EXPECT_EQ(7254450u, util::tileCount(sanFrancisco, 22, util::tileSize)); + EXPECT_EQ(7254450u, util::tileCount(sanFrancisco, 22)); } +TEST(TileCount, BoundsCrossingAntimeridian) { + auto crossingBounds = LatLngBounds::hull({-20.9615, -214.309}, {19.477, -155.830}); + + EXPECT_EQ(1u, util::tileCount(crossingBounds, 0)); + EXPECT_EQ(4u, util::tileCount(crossingBounds, 3)); + EXPECT_EQ(8u, util::tileCount(crossingBounds, 4)); +} diff --git a/test/util/tile_range.test.cpp b/test/util/tile_range.test.cpp index c4c37c74d7..83ac5096a5 100644 --- a/test/util/tile_range.test.cpp +++ b/test/util/tile_range.test.cpp @@ -52,14 +52,14 @@ TEST(TileRange, ContainsWrappedBounds) { TEST(TileRange, ContainsBoundsCrossingAntimeridian) { { - auto cossingBounds = LatLngBounds::hull({-20.9615, -214.309}, {19.477, -155.830}); - auto range = util::TileRange::fromLatLngBounds(cossingBounds, 1); + auto crossingBounds = LatLngBounds::hull({-20.9615, -214.309}, {19.477, -155.830}); + auto range = util::TileRange::fromLatLngBounds(crossingBounds, 1); EXPECT_TRUE(range.contains(CanonicalTileID(1, 1, 1))); EXPECT_TRUE(range.contains(CanonicalTileID(1, 0, 0))); } { - auto cossingBounds = LatLngBounds::hull({-20.9615, -214.309}, {19.477, -155.830}); - auto range = util::TileRange::fromLatLngBounds(cossingBounds, 6); + auto crossingBounds = LatLngBounds::hull({-20.9615, -214.309}, {19.477, -155.830}); + auto range = util::TileRange::fromLatLngBounds(crossingBounds, 6); EXPECT_FALSE(range.contains(CanonicalTileID(6, 55, 34))); EXPECT_FALSE(range.contains(CanonicalTileID(6, 5, 28))); EXPECT_TRUE(range.contains(CanonicalTileID(6, 63, 28))); -- cgit v1.2.1