[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.

12 files changed, 886 insertions, 55 deletions

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 <benchmark/benchmark.h>
+#include <mbgl/util/geo.hpp>
+#include <mbgl/util/geometry.hpp>
+#include <mbgl/util/tile_coordinate.hpp>
+#include <mbgl/util/tile_cover.hpp>
+#include <mbgl/map/transform.hpp>
+
+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<double>{
+    {
+        {-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<double> project(const LatLng& latLng, uint8_t zoom) {
-        return project_(latLng, std::pow(2.0, zoom));
+    //Returns point on tile
+    static Point<double> project(const LatLng& latLng, int32_t zoom) {
+        return project_(latLng, 1 << zoom);
     }
 
     static LatLng unproject(const Point<double>& 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<double> project(const LatLng& point, uint8_t zoom, uint16_t tileSize) {
-        const double t2z = tileSize * std::pow(2, zoom);
-        Point<double> 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<double> project_(const LatLng& latLng, double worldSize) {
         return Point<double> {
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<uint8_t> 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 <mbgl/programs/programs.hpp>
 #include <mbgl/programs/background_program.hpp>
 #include <mbgl/util/tile_cover.hpp>
+#include <mbgl/map/transform_state.hpp>
 
 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 <mbgl/util/constants.hpp>
 #include <mbgl/util/interpolate.hpp>
 #include <mbgl/map/transform_state.hpp>
+#include <mbgl/util/tile_cover_impl.hpp>
+#include <mbgl/util/tile_coordinate.hpp>
 
 #include <functional>
+#include <list>
 
 namespace mbgl {
 
 namespace {
 
+using ScanLine = const std::function<void(int32_t x0, int32_t x1, int32_t y)>;
+
 // 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<void(int32_t x0, int32_t x1, int32_t y)>;
-
 // 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<UnwrappedTileID> 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<UnwrappedTileID> tileCover(const TransformState& state, int32_t z) {
         z);
 }
 
+std::vector<UnwrappedTileID> tileCover(const Geometry<double>& geometry, int32_t z) {
+    std::vector<UnwrappedTileID> 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<double>& 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<double> p({ {sw, nw, ne, se, sw} });
+    impl = std::make_unique<TileCover::Impl>(z, p, false);
+}
+
+TileCover::TileCover(const Geometry<double>& geom, int32_t z, bool project/* = true*/)
+ : impl( std::make_unique<TileCover::Impl>(z, geom, project)) {
+}
+
+TileCover::~TileCover() {
+
+}
+
+optional<UnwrappedTileID> 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 <mbgl/tile/tile_id.hpp>
 #include <mbgl/style/types.hpp>
-#include <mbgl/util/tile_coordinate.hpp>
+#include <mbgl/util/geometry.hpp>
+#include <mbgl/util/optional.hpp>
 
 #include <vector>
+#include <memory>
 
 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<double>&, int32_t z, bool project = true);
+    ~TileCover();
+
+    optional<UnwrappedTileID> next();
+    bool hasNext();
+
+private:
+    class Impl;
+    std::unique_ptr<Impl> impl;
+};
+
 int32_t coveringZoomLevel(double z, style::SourceType type, uint16_t tileSize);
 
 std::vector<UnwrappedTileID> tileCover(const TransformState&, int32_t z);
 std::vector<UnwrappedTileID> tileCover(const LatLngBounds&, int32_t z);
+std::vector<UnwrappedTileID> tileCover(const Geometry<double>&, 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<double>&, 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 <mbgl/util/tile_cover_impl.hpp>
+#include <mbgl/util/tile_coordinate.hpp>
+
+#include <functional>
+#include <cmath>
+#include <assert.h>
+#include <limits.h>
+#include <algorithm>
+
+namespace mbgl {
+namespace util {
+
+using PointList = std::vector<Point<double>>;
+
+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::size_t>(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::size_t>(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<uint32_t>(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<uint32_t>(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<int32_t>(std::floor(x)));
+    xp.xmax = std::max(xp.xmax, static_cast<int32_t>(std::ceil(x)));
+}
+
+//Build a vector of X tile-coordinates spanned by each bound.
+std::vector<TileSpan> scan_row(uint32_t y, Bounds& aet) {
+    std::vector<TileSpan> 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<Point<double>>& 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<double>& polygon, BoundsMap& et) const {
+        for(const auto&ring : polygon) {
+            buildTable(ring, et, true);
+        }
+    }
+    BoundsMap operator()(const Point<double>&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<uint32_t>(std::floor(clamp(point.y, 0.0, (double)(1 << zoom))));
+        et[y].push_back(bnd);
+        return et;
+    }
+
+    BoundsMap operator()(const MultiPoint<double>& points) const {
+        BoundsMap et;
+        for (const Point<double>& 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<uint32_t>(std::floor(clamp(point.y, 0.0, (double)(1 << zoom))));
+            et[y].push_back(bnd);
+        }
+        return et;
+    }
+
+    BoundsMap operator()(const LineString<double>& lines) const {
+        BoundsMap et;
+        buildTable(lines, et);
+        return et;
+    }
+
+    BoundsMap operator()(const MultiLineString<double>& lines) const {
+        BoundsMap et;
+        for(const auto&line : lines) {
+            buildTable(line, et);
+        }
+        return et;
+    }
+
+    BoundsMap operator()(const Polygon<double>& polygon) const {
+        BoundsMap et;
+        buildPolygonTable(polygon, et);
+        return et;
+    }
+
+    BoundsMap operator()(const MultiPolygon<double>& polygons) const {
+        BoundsMap et;
+        for(const auto& polygon: polygons) {
+            buildPolygonTable(polygon, et);
+        }
+        return et;
+    }
+
+    BoundsMap operator()(const mapbox::geometry::geometry_collection<double>&) const {
+        return {};
+    }
+};
+
+TileCover::Impl::Impl(int32_t z, const Geometry<double>& 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<UnwrappedTileID> 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 <mbgl/util/tile_cover.hpp>
+#include <mbgl/util/geometry.hpp>
+#include <mbgl/util/optional.hpp>
+
+#include <vector>
+#include <map>
+#include <queue>
+
+namespace mbgl {
+
+class TransformState;
+class LatLngBounds;
+
+namespace util {
+
+struct Bound;
+
+using Bounds = std::vector<Bound>;
+using BoundsMap = std::map<uint32_t, Bounds>;
+
+// A chain of points from a local minimum to a local maximum. `winding` indicates
+//  the direction of the original geometry.
+struct Bound {
+    std::vector<Point<double>> 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<double>& geom, bool project = true);
+    ~Impl() = default;
+
+    optional<UnwrappedTileID> next();
+    bool hasNext() const;
+
+private:
+    using TileSpans = std::queue<std::pair<int32_t, int32_t>>;
+
+    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 <mbgl/util/geo.hpp>
 #include <mbgl/map/transform.hpp>
 
+#include <algorithm>
+
 #include <gtest/gtest.h>
 
 using namespace mbgl;
@@ -22,8 +24,8 @@ TEST(TileCover, Antarctic) {
 
 TEST(TileCover, WorldZ0) {
     EXPECT_EQ((std::vector<UnwrappedTileID>{
-                  { 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<UnwrappedTileID>{
-                  { 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<UnwrappedTileID>{
-                  { 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<UnwrappedTileID> 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<UnwrappedTileID> t;
+    while(tc.hasNext()) {
+        t.push_back(*tc.next());
+    };
+    EXPECT_EQ((std::vector<UnwrappedTileID>{
+        { 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<UnwrappedTileID>{
-                  { 0, 0, 0 },
-              }),
+        { 0, 0, 0 },
+    }),
               util::tileCover(sanFrancisco, 0));
 }
 
 TEST(TileCover, SanFranciscoZ10) {
     EXPECT_EQ((std::vector<UnwrappedTileID>{
-                  { 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<double>{ -122.5744, 37.6609 };
+
+    EXPECT_EQ((std::vector<UnwrappedTileID>{ {2 ,0 ,1 } }),
+              util::tileCover(point, 2));
+}
+
+TEST(TileCover, GeomMultiPoint) {
+    auto points = MultiPoint<double>{ { -122.5, 37.76 }, { -122.4, 37.76} };
+
+    EXPECT_EQ((std::vector<UnwrappedTileID>{
+        {20, 167480, 405351}, {20, 167772, 405351} }),
+              util::tileCover(points, 20));
+}
+
+TEST(TileCover, GeomLineZ10) {
+    auto lineCover = util::tileCover(LineString<double>{
+        {-121.49368286132812,38.57903714667459},
+        {-122.4422836303711,37.773157169570695}
+    }, 10);
+    EXPECT_EQ((std::vector<UnwrappedTileID>{
+        { 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<double>{
+        {-179.93342914581299,38.892101707724315},
+        {-180.02394485473633,38.89203490311832}
+    }, 10);
+    EXPECT_EQ((std::vector<UnwrappedTileID>{ { 10, -1, 391 }, { 10, 0, 391 } }),
+              lineCover);
+    
+    lineCover = util::tileCover(LineString<double>{
+        {179.93342914581299,38.892101707724315},
+        {180.02394485473633,38.89203490311832}
+    }, 10);
+    EXPECT_EQ((std::vector<UnwrappedTileID>{ { 10, 1023, 391 }, { 10, 1024, 391 } }),
+              lineCover);
+}
+
+TEST(TileCover, GeomMultiLineString) {
+    auto geom = MultiLineString<double>{
+        { { -122.5, 37.76 }, { -122.4, 37.76} },
+        { { -122.5, 37.72 }, { -122.4, 37.72} } };
+
+    EXPECT_EQ((std::vector<UnwrappedTileID>{
+        {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<double>{
+        {
+            {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<double>{
+        {{
+            {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<double>{
+        {
+            {-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<UnwrappedTileID>{
+        { 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<double>{ -122.5744, 97.6609 };
+    EXPECT_THROW(util::tileCover(point, 2), std::domain_error);
+
+    auto badPoly = Polygon<double> { { {1.0,  35.0} } };
+    EXPECT_EQ((std::vector<UnwrappedTileID>{ }), util::tileCover(badPoly, 16));
+
+    //Should handle open polygons.
+    badPoly = Polygon<double> { { {1.0,  34.2}, {1.0, 34.4}, {0.5, 34.3} } };
+    EXPECT_EQ((std::vector<UnwrappedTileID>{
+        { 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)));