[core] add tileCover overload for polygons

It is now possible to calculate the tile coverage of arbitrary polygons (including ones with holes). The implementation treats (0,0) as the polygon's center.

3 files changed, 344 insertions, 0 deletions

diff --git a/src/mbgl/util/tile_cover.cpp b/src/mbgl/util/tile_cover.cpp
index 5487fb269c..01de65f8c2 100644
--- a/src/mbgl/util/tile_cover.cpp
+++ b/src/mbgl/util/tile_cover.cpp
@@ -1,10 +1,25 @@
 #include <mbgl/util/tile_cover.hpp>
 #include <mbgl/util/constants.hpp>
 #include <mbgl/util/interpolate.hpp>
+#include <mbgl/util/geometry.hpp>
 #include <mbgl/map/transform_state.hpp>
 
+#include <mapbox/earcut.hpp>
+
 #include <functional>
 
+namespace mapbox {
+namespace util {
+    template <> struct nth<0, mbgl::Point<double>> {
+        static double get(const mbgl::Point<double>& t) { return t.x; };
+    };
+    
+    template <> struct nth<1, mbgl::Point<double>> {
+        static double get(const mbgl::Point<double>& t) { return t.y; };
+    };
+} // namespace util
+} // namespace mapbox
+
 namespace mbgl {
 
 namespace {
@@ -24,6 +39,11 @@ struct edge {
         y1 = b.y;
         dx = b.x - a.x;
         dy = b.y - a.y;
+        
+        assert(std::isfinite(x0));
+        assert(std::isfinite(x1));
+        assert(std::isfinite(y0));
+        assert(std::isfinite(y1));
     }
 };
 
@@ -74,6 +94,22 @@ static void scanTriangle(const Point<double>& a, const Point<double>& b, const P
 namespace util {
 
 namespace {
+    
+std::pair<size_t, size_t> ringAndPointIndicesFromFlatIndex(const Polygon<double>& polygon, size_t flatIndex) {
+    size_t ringIndex = 0;
+    while (flatIndex >= polygon[ringIndex].size()) {
+        flatIndex -= polygon[ringIndex].size();
+        ++ringIndex;
+    }
+    assert(ringIndex < polygon.size());
+    assert(flatIndex < polygon[ringIndex].size());
+    return std::make_pair(ringIndex, flatIndex);
+}
+    
+const Point<double>& pointAtFlatIndex(const Polygon<double>& polygon, size_t flatIndex) {
+    std::pair<size_t, size_t> ringAndPointIndices = ringAndPointIndicesFromFlatIndex(polygon, flatIndex);
+    return polygon[ringAndPointIndices.first][ringAndPointIndices.second];
+}
 
 std::vector<UnwrappedTileID> tileCover(const Point<double>& tl,
                                        const Point<double>& tr,
@@ -124,6 +160,70 @@ std::vector<UnwrappedTileID> tileCover(const Point<double>& tl,
     }
     return result;
 }
+    
+std::vector<UnwrappedTileID> tileCover(const Polygon<double> &poly, const Point<double> &center, int32_t z) {
+    const int32_t tiles = 1 << z;
+    
+    struct ID {
+        int32_t x, y;
+        double sqDist;
+    };
+    
+    std::vector<ID> t;
+    
+    auto scanLine = [&](int32_t x0, int32_t x1, int32_t y) {
+        int32_t x;
+        if (y >= 0 && y <= tiles) {
+            for (x = x0; x < x1; ++x) {
+                const auto dx = x + 0.5 - center.x, dy = y + 0.5 - center.y;
+                t.emplace_back(ID{ x, y, dx * dx + dy * dy });
+            }
+        }
+    };
+    
+    std::vector<uint32_t> indices = mapbox::earcut(poly);
+    assert(indices.size() % 3 == 0);
+    assert(indices.size() > 0);
+    
+    for (size_t i = 0; i < indices.size(); i += 3) {
+        // This function is sloooow but I want to just get things working before figuring out how to optimize it
+        Point<double> a = pointAtFlatIndex(poly, indices[i]);
+        Point<double> b = pointAtFlatIndex(poly, indices[i + 1]);
+        Point<double> c = pointAtFlatIndex(poly, indices[i + 2]);
+        scanTriangle(a, b, c, 0, tiles, scanLine);
+    }
+    
+    // Sort first by distance, then by x/y.
+    std::sort(t.begin(), t.end(), [](const ID& a, const ID& b) {
+        return (a.sqDist != b.sqDist) ? (a.sqDist < b.sqDist)
+        : ((a.x != b.x) ? (a.x < b.x) : (a.y < b.y));
+    });
+    
+    // Erase duplicate tile IDs (they typically occur at the common side of both triangles).
+    t.erase(std::unique(t.begin(), t.end(), [](const ID& a, const ID& b) {
+        return a.x == b.x && a.y == b.y;
+    }), t.end());
+    
+    std::vector<UnwrappedTileID> result;
+    for (const auto& id : t) {
+        result.emplace_back(z, id.x, id.y);
+    }
+    return result;
+}
+    
+Polygon<double> fromLatLng(const Polygon<double>& polygon, int32_t z) {
+    TransformState state;
+    Polygon<double> result;
+    for (const LinearRing<double>& ring : polygon) {
+        result.emplace_back();
+        for (const Point<double>& point : ring) {
+            result.back().push_back(TileCoordinate::fromLatLng(state, z, { point.y, point.x }).p);
+            assert(std::isfinite(result.back().back().x));
+            assert(std::isfinite(result.back().back().y));
+        }
+    }
+    return result;
+}
 
 } // namespace
 
@@ -156,6 +256,12 @@ std::vector<UnwrappedTileID> tileCover(const LatLngBounds& bounds_, int32_t z) {
         TileCoordinate::fromLatLng(state, z, bounds.center()).p,
         z);
 }
+    
+std::vector<UnwrappedTileID> tileCover(const Polygon<double>& polygon, int32_t z) {
+    if (polygon.size() == 0) { return {}; }
+    Point<double> center;
+    return tileCover(fromLatLng(polygon, z), center, z);
+}
 
 std::vector<UnwrappedTileID> tileCover(const TransformState& state, int32_t z) {
     const double w = state.getWidth();
diff --git a/src/mbgl/util/tile_cover.hpp b/src/mbgl/util/tile_cover.hpp
index 2d32d8bf41..e0a8c2a380 100644
--- a/src/mbgl/util/tile_cover.hpp
+++ b/src/mbgl/util/tile_cover.hpp
@@ -3,6 +3,7 @@
 #include <mbgl/tile/tile_id.hpp>
 #include <mbgl/style/types.hpp>
 #include <mbgl/util/tile_coordinate.hpp>
+#include <mbgl/util/geometry.hpp>
 
 #include <vector>
 
@@ -16,6 +17,7 @@ namespace util {
 int32_t coveringZoomLevel(double z, SourceType type, uint16_t tileSize);
 
 std::vector<UnwrappedTileID> tileCover(const TransformState&, int32_t z);
+std::vector<UnwrappedTileID> tileCover(const Polygon<double>&, int32_t z);
 std::vector<UnwrappedTileID> tileCover(const LatLngBounds&, int32_t z);
 
 } // namespace util
diff --git a/test/util/tile_cover.cpp b/test/util/tile_cover.cpp
index cc183509d9..9a0a748cb4 100644
--- a/test/util/tile_cover.cpp
+++ b/test/util/tile_cover.cpp
@@ -6,6 +6,29 @@
 
 using namespace mbgl;
 
+namespace {
+    
+bool containsTile(const std::vector<UnwrappedTileID>& tiles, uint8_t z, uint32_t x, uint32_t y) {
+    UnwrappedTileID target(z, x, y);
+    return std::find(tiles.cbegin(), tiles.cend(), target) != tiles.cend();
+}
+
+Polygon<double> toLatLng(const Polygon<double>& polygon, double zoom) {
+    const TransformState state;
+    Polygon<double> result;
+    for (const LinearRing<double>& ring : polygon) {
+        result.emplace_back();
+        for (const Point<double>& point : ring) {
+            TileCoordinate coord{point, zoom};
+            LatLng latLng = coord.toLatLng(state);
+            result.back().emplace_back(latLng.longitude, latLng.latitude);
+        }
+    }
+    return result;
+}
+    
+}
+
 TEST(TileCover, Empty) {
     EXPECT_EQ((std::vector<UnwrappedTileID>{}), util::tileCover(LatLngBounds::empty(), 0));
 }
@@ -81,3 +104,216 @@ TEST(TileCover, SanFranciscoZ0Wrapped) {
     EXPECT_EQ((std::vector<UnwrappedTileID>{ { 0, 1, 0 } }),
               util::tileCover(sanFranciscoWrapped, 0));
 }
+
+TEST(TileCover, PolygonAboveArcticCircleZ1) {
+    Polygon<double> polygon { {
+        { 0.0, 87.0 },
+        { 30.0, 87.0 },
+        { 40.0, 88.0 },
+        { -10.0, 88.0 }
+    } };
+    
+    std::vector<UnwrappedTileID> tiles = util::tileCover(polygon, 1);
+    
+    EXPECT_EQ(0, tiles.size());
+}
+
+
+TEST(TileCover, SquarePolygonZ1) {
+    // A simple square that crosses tile boundaries
+    Polygon<double> polygon { {
+        { 0.5, 0.5 },
+        { 1.5, 0.5 },
+        { 1.5, 1.5 },
+        { 0.5, 1.5 }
+    } };
+    Polygon<double> latLng = toLatLng(polygon, 1);
+    
+    std::vector<UnwrappedTileID> tiles = util::tileCover(latLng, 1);
+    
+    EXPECT_EQ(4, tiles.size());
+    EXPECT_TRUE(containsTile(tiles, 1, 0, 0));
+    EXPECT_TRUE(containsTile(tiles, 1, 0, 1));
+    EXPECT_TRUE(containsTile(tiles, 1, 1, 0));
+    EXPECT_TRUE(containsTile(tiles, 1, 1, 1));
+}
+
+TEST(TileCover, UShapeZ2) {
+    // A block "U"
+    Polygon<double> polygon { {
+        { 0.1, 0.1 },
+        { 2.9, 0.1 },
+        { 2.9, 2.9 },
+        { 2.1, 2.9 },
+        { 2.1, 0.9 },
+        { 0.9, 0.9 },
+        { 0.9, 2.9 },
+        { 0.1, 2.9 }
+    } };
+    Polygon<double> latLng = toLatLng(polygon, 2);
+    
+    std::vector<UnwrappedTileID> tiles = util::tileCover(latLng, 2);
+    
+    EXPECT_EQ(7, tiles.size());
+    EXPECT_TRUE(containsTile(tiles, 2, 0, 0));
+    EXPECT_TRUE(containsTile(tiles, 2, 0, 1));
+    EXPECT_TRUE(containsTile(tiles, 2, 0, 2));
+    EXPECT_TRUE(containsTile(tiles, 2, 1, 0));
+    EXPECT_TRUE(containsTile(tiles, 2, 2, 0));
+    EXPECT_TRUE(containsTile(tiles, 2, 2, 1));
+    EXPECT_TRUE(containsTile(tiles, 2, 2, 2));
+}
+
+TEST(TileCover, DiagonalEdgeZ2) {
+    // A concave shape that crosses tile boundaries on diagonals
+    Polygon<double> polygon { {
+        { 2.9, 0.9 },
+        { 3.1, 2.0 },
+        { 2.9, 3.1 },
+        { 3.6, 3.1 },
+        { 3.6, 0.9 }
+    } };
+    Polygon<double> latLng = toLatLng(polygon, 2);
+    
+    std::vector<UnwrappedTileID> tiles = util::tileCover(latLng, 2);
+    
+    EXPECT_EQ(8, tiles.size());
+    EXPECT_TRUE(containsTile(tiles, 2, 2, 0));
+    EXPECT_TRUE(containsTile(tiles, 2, 2, 1));
+    EXPECT_TRUE(containsTile(tiles, 2, 2, 2));
+    EXPECT_TRUE(containsTile(tiles, 2, 2, 3));
+    EXPECT_TRUE(containsTile(tiles, 2, 3, 0));
+    EXPECT_TRUE(containsTile(tiles, 2, 3, 1));
+    EXPECT_TRUE(containsTile(tiles, 2, 3, 2));
+    EXPECT_TRUE(containsTile(tiles, 2, 3, 3));
+}
+
+TEST(TileCover, HoleContainedInTileZ1) {
+    // A shape with a hole that is entirely contained within a single tile
+    Polygon<double> polygon {
+        {
+            { 0.2, 0.2 },
+            { 0.8, 0.2 },
+            { 0.8, 0.8 },
+            { 0.2, 0.8 }
+        },
+        {
+            { 0.4, 0.4 },
+            { 0.4, 0.6 },
+            { 0.6, 0.6 },
+            { 0.6, 0.4 }
+        }
+    };
+    Polygon<double> latLng = toLatLng(polygon, 1);
+    
+    std::vector<UnwrappedTileID> tiles = util::tileCover(latLng, 1);
+    
+    EXPECT_EQ(1, tiles.size());
+    EXPECT_TRUE(containsTile(tiles, 1, 0, 0));
+}
+
+TEST(TileCover, HoleSpanningTileBoundaryZ2) {
+    // A shape with a hole that crosses tile boundaries, but is not so large
+    // that it causes any tiles to actually be excluded
+    Polygon<double> polygon {
+        {
+            { 0.2, 0.2 },
+            { 1.8, 0.2 },
+            { 1.8, 1.8 },
+            { 0.2, 1.8 }
+        },
+        {
+            { 0.4, 0.4 },
+            { 0.4, 1.6 },
+            { 1.6, 1.6 },
+            { 1.6, 0.4 }
+        }
+    };
+    Polygon<double> latLng = toLatLng(polygon, 2);
+    
+    std::vector<UnwrappedTileID> tiles = util::tileCover(latLng, 2);
+    
+    EXPECT_EQ(4, tiles.size());
+    EXPECT_TRUE(containsTile(tiles, 2, 0, 0));
+    EXPECT_TRUE(containsTile(tiles, 2, 0, 1));
+    EXPECT_TRUE(containsTile(tiles, 2, 1, 0));
+    EXPECT_TRUE(containsTile(tiles, 2, 1, 1));
+}
+
+TEST(TileCover, HoleLargeEnoughToExcludeTileZ2) {
+    // A shape with a hole that is large enough to cause tiles to actually be
+    // excluded.
+    Polygon<double> polygon {
+        {
+            { 0.2, 0.2 },
+            { 3.8, 0.2 },
+            { 3.8, 3.8 },
+            { 0.2, 3.8 }
+        },
+        {
+            { 0.4, 0.4 },
+            { 0.4, 3.6 },
+            { 3.6, 3.6 },
+            { 3.6, 0.4 }
+        }
+    };
+    Polygon<double> latLng = toLatLng(polygon, 2);
+    
+    std::vector<UnwrappedTileID> tiles = util::tileCover(latLng, 2);
+    
+    EXPECT_EQ(12, tiles.size());
+    EXPECT_TRUE(containsTile(tiles, 2, 0, 0));
+    EXPECT_TRUE(containsTile(tiles, 2, 0, 1));
+    EXPECT_TRUE(containsTile(tiles, 2, 0, 2));
+    EXPECT_TRUE(containsTile(tiles, 2, 0, 3));
+    EXPECT_TRUE(containsTile(tiles, 2, 3, 0));
+    EXPECT_TRUE(containsTile(tiles, 2, 3, 1));
+    EXPECT_TRUE(containsTile(tiles, 2, 3, 2));
+    EXPECT_TRUE(containsTile(tiles, 2, 3, 3));
+    EXPECT_TRUE(containsTile(tiles, 2, 1, 0));
+    EXPECT_TRUE(containsTile(tiles, 2, 2, 0));
+    EXPECT_TRUE(containsTile(tiles, 2, 1, 3));
+    EXPECT_TRUE(containsTile(tiles, 2, 2, 3));
+}
+
+TEST(TileCover, PolygonBeyondAntimeridianZ0) {
+    Polygon<double> latLng { {
+        { 190, 10 },
+        { 200, 10 },
+        { 200, 30 },
+        { 190, 30 }
+    } };
+    
+    std::vector<UnwrappedTileID> tiles = util::tileCover(latLng, 0);
+    
+    EXPECT_EQ(1, tiles.size());
+    EXPECT_EQ(0, tiles[0].canonical.z);
+    EXPECT_EQ(0, tiles[0].canonical.x);
+    EXPECT_EQ(0, tiles[0].canonical.y);
+    EXPECT_EQ(1, tiles[0].wrap);
+}
+
+TEST(TileCover, PolygonSpanningAntimeridianZ1) {
+    Polygon<double> latLng { {
+        { 170, -20 },
+        { 190, -20 },
+        { 190, 20 },
+        { 170, 20 }
+    } };
+    
+    std::vector<UnwrappedTileID> tiles = util::tileCover(latLng, 1);
+    
+    EXPECT_EQ(4, tiles.size());
+    EXPECT_TRUE(std::find_if(tiles.cbegin(), tiles.cend(), [](const UnwrappedTileID& tile) {
+        return tile.canonical.x == 0 && tile.canonical.y == 0 && tile.wrap == 1;
+    }) != tiles.end());
+    EXPECT_TRUE(std::find_if(tiles.cbegin(), tiles.cend(), [](const UnwrappedTileID& tile) {
+        return tile.canonical.x == 1 && tile.canonical.y == 0 && tile.wrap == 0;
+    }) != tiles.end());
+    EXPECT_TRUE(std::find_if(tiles.cbegin(), tiles.cend(), [](const UnwrappedTileID& tile) {
+        return tile.canonical.x == 0 && tile.canonical.y == 1 && tile.wrap == 1;
+    }) != tiles.end());
+    EXPECT_TRUE(std::find_if(tiles.cbegin(), tiles.cend(), [](const UnwrappedTileID& tile) {
+        return tile.canonical.x == 1 && tile.canonical.y == 1 && tile.wrap == 0;
+    }) != tiles.end());
+}