[core] util::tileCover optimization: three scans and no duplicates handling.

benchmark/util/tilecover.benchmark.cpp

@@ -22,7 +22,8 @@ static void TileCoverPitchedViewport(benchmark::State& state) {
     Transform transform;
     transform.resize({ 512, 512 });
     // slightly offset center so that tile order is better defined
-    transform.jumpTo(CameraOptions().withCenter(LatLng { 0.1, -0.1 }).withZoom(8.0).withBearing(5.0).withPitch(40.0));
+    transform.jumpTo(CameraOptions().withCenter(LatLng { 0.1, -0.1 }).withPadding(EdgeInsets { 376, 0, 0, 0 })
+                                    .withZoom(8.0).withBearing(5.0).withPitch(60.0));
 
     std::size_t length = 0;
     while (state.KeepRunning()) {

include/mbgl/util/projection.hpp

@@ -98,7 +98,7 @@ class Projection {
         return Point<double> {
             util::LONGITUDE_MAX + latLng.longitude(),
             util::LONGITUDE_MAX - util::RAD2DEG * std::log(std::tan(M_PI / 4 + latitude * M_PI / util::DEGREES_MAX))
-        } * worldSize / util::DEGREES_MAX;
+        } * (worldSize / util::DEGREES_MAX);
     }
 };
 

src/mbgl/map/transform.cpp

@@ -138,9 +138,6 @@ void Transform::easeTo(const CameraOptions& camera, const AnimationOptions& anim
         if (bearing != startBearing) {
             state.bearing = util::wrap(util::interpolate(startBearing, bearing, t), -M_PI, M_PI);
         }
-        if (pitch != startPitch) {
-            state.pitch = util::interpolate(startPitch, pitch, t);
-        }
         if (padding != startEdgeInsets) {
             // Interpolate edge insets
             state.edgeInsets = {
@@ -150,6 +147,10 @@ void Transform::easeTo(const CameraOptions& camera, const AnimationOptions& anim
                 util::interpolate(startEdgeInsets.right(), padding.right(), t)
             };
         }
+        auto maxPitch = getMaxPitchForEdgeInsets(state.edgeInsets);
+        if (pitch != startPitch || maxPitch < startPitch) {
+            state.pitch = std::min(maxPitch, util::interpolate(startPitch, pitch, t));
+        }
     }, duration);
 }
 
@@ -302,9 +303,6 @@ void Transform::flyTo(const CameraOptions &camera, const AnimationOptions &anima
         if (bearing != startBearing) {
             state.bearing = util::wrap(util::interpolate(startBearing, bearing, k), -M_PI, M_PI);
         }
-        if (pitch != startPitch) {
-            state.pitch = util::interpolate(startPitch, pitch, k);
-        }
         if (padding != startEdgeInsets) {
             // Interpolate edge insets
             state.edgeInsets = {
@@ -314,6 +312,10 @@ void Transform::flyTo(const CameraOptions &camera, const AnimationOptions &anima
                 util::interpolate(startEdgeInsets.right(), padding.right(), us)
             };
         }
+        auto maxPitch = getMaxPitchForEdgeInsets(state.edgeInsets);
+        if (pitch != startPitch || maxPitch < startPitch) {
+            state.pitch = std::min(maxPitch, util::interpolate(startPitch, pitch, k));
+        }
     }, duration);
 }
 
@@ -576,4 +578,16 @@ LatLng Transform::screenCoordinateToLatLng(const ScreenCoordinate& point, LatLng
     return state.screenCoordinateToLatLng(flippedPoint, wrapMode);
 }
 
+double Transform::getMaxPitchForEdgeInsets(const EdgeInsets &insets) const
+{
+    double centerOffsetY = 0.5 * (insets.top() - insets.bottom()); // See TransformState::getCenterOffset.
+
+    const auto height = state.size.height;
+    assert(height);
+    // See TransformState::fov description: fov = 2 * arctan((height / 2) / (height * 1.5)).
+    const double fovAboveCenter = std::atan((0.666666 + 0.02) * (0.5 + centerOffsetY / height));
+    return M_PI * 0.5 - fovAboveCenter;   // 0.02 added ^^^^ to prevent parallel ground to viewport clipping plane.
+    // e.g. Maximum pitch of 60 degrees is when perspective center's offset from the top is 84% of screen height.
+}
+
 } // namespace mbgl

src/mbgl/map/transform.hpp

@@ -115,6 +115,9 @@ class Transform : private util::noncopyable {
                          std::function<void(double)>,
                          const Duration&);
 
+    // We don't want to show horizon: limit max pitch based on edge insets.
+    double getMaxPitchForEdgeInsets(const EdgeInsets &insets) const;
+
     TimePoint transitionStart;
     Duration transitionDuration;
     std::function<bool(const TimePoint)> transitionFrameFn;

src/mbgl/map/transform_state.cpp

@@ -36,18 +36,18 @@ void TransformState::getProjMatrix(mat4& projMatrix, uint16_t nearZ, bool aligne
     const double cameraToCenterDistance = getCameraToCenterDistance();
     auto offset = getCenterOffset();
 
-    // Find the distance from the viewport center point
-    // [width/2 + offset.x, height/2 + offset.y] to the top edge, to point
-    // [width/2 + offset.x, 0] in Z units, using the law of sines.
+    // Find the Z distance from the viewport center point
+    // [width/2 + offset.x, height/2 + offset.y] to the top edge; to point
+    // [width/2 + offset.x, 0] in Z units.
     // 1 Z unit is equivalent to 1 horizontal px at the center of the map
     // (the distance between[width/2, height/2] and [width/2 + 1, height/2])
-    const double fovAboveCenter = getFieldOfView() * (0.5 + offset.y / size.height);
-    const double groundAngle = M_PI / 2.0 + getPitch();
-    const double aboveCenterSurfaceDistance = std::sin(fovAboveCenter) * cameraToCenterDistance / std::sin(M_PI - groundAngle - fovAboveCenter);
-
-
+    // See https://github.com/mapbox/mapbox-gl-native/pull/15195 for details.
+    // See TransformState::fov description: fov = 2 * arctan((height / 2) / (height * 1.5)).
+    const double tanFovAboveCenter = (size.height * 0.5 + offset.y) / (size.height * 1.5);
+    const double tanMultiple = tanFovAboveCenter * std::tan(getPitch());
+    assert(tanMultiple < 1);
     // Calculate z distance of the farthest fragment that should be rendered.
-    const double furthestDistance = std::cos(M_PI / 2 - getPitch()) * aboveCenterSurfaceDistance + cameraToCenterDistance;
+    const double furthestDistance = cameraToCenterDistance / (1 - tanMultiple);
     // Add a bit extra to avoid precision problems when a fragment's distance is exactly `furthestDistance`
     const double farZ = furthestDistance * 1.01;
 
@@ -64,7 +64,7 @@ void TransformState::getProjMatrix(mat4& projMatrix, uint16_t nearZ, bool aligne
     const bool flippedY = viewportMode == ViewportMode::FlippedY;
     matrix::scale(projMatrix, projMatrix, 1.0, flippedY ? 1 : -1, 1);
 
-    matrix::translate(projMatrix, projMatrix, 0, 0, -getCameraToCenterDistance());
+    matrix::translate(projMatrix, projMatrix, 0, 0, -cameraToCenterDistance);
 
     using NO = NorthOrientation;
     switch (getNorthOrientation()) {

src/mbgl/util/tile_cover.cpp

@@ -39,15 +39,15 @@ static void scanSpans(edge e0, edge e1, int32_t ymin, int32_t ymax, ScanLine sca
     double y1 = ::fmin(ymax, std::ceil(e1.y1));
 
     // sort edges by x-coordinate
-    if ((e0.x0 == e1.x0 && e0.y0 == e1.y0) ?
-        (e0.x0 + e1.dy / e0.dy * e0.dx < e1.x1) :
-        (e0.x1 - e1.dy / e0.dy * e0.dx < e1.x0)) {
+    double m0 = e0.dx / e0.dy;
+    double m1 = e1.dx / e1.dy;
+    double ySort = e0.y0 == e1.y0 ? std::min(e0.y1, e1.y1) : std::max(e0.y0, e1.y0);
+    if (e0.x0 - (e0.y0 - ySort) * m0 < e1.x0 - (e1.y0 - ySort) * m1) {
         std::swap(e0, e1);
+        std::swap(m0, m1);
     }
 
     // scan lines!
-    double m0 = e0.dx / e0.dy;
-    double m1 = e1.dx / e1.dy;
     double d0 = e0.dx > 0; // use y + 1 to compute x0
     double d1 = e1.dx < 0; // use y + 1 to compute x1
     for (int32_t y = y0; y < y1; y++) {
@@ -57,22 +57,6 @@ static void scanSpans(edge e0, edge e1, int32_t ymin, int32_t ymax, ScanLine sca
     }
 }
 
-// scan-line conversion
-static void scanTriangle(const Point<double>& a, const Point<double>& b, const Point<double>& c, int32_t ymin, int32_t ymax, ScanLine& scanLine) {
-    edge ab = edge(a, b);
-    edge bc = edge(b, c);
-    edge ca = edge(c, a);
-
-    // sort edges by y-length
-    if (ab.dy > bc.dy) { std::swap(ab, bc); }
-    if (ab.dy > ca.dy) { std::swap(ab, ca); }
-    if (bc.dy > ca.dy) { std::swap(bc, ca); }
-
-    // scan span! scan span!
-    if (ab.dy) scanSpans(ca, ab, ymin, ymax, scanLine);
-    if (bc.dy) scanSpans(ca, bc, ymin, ymax, scanLine);
-}
-
 } // namespace
 
 namespace util {
@@ -85,7 +69,7 @@ std::vector<UnwrappedTileID> tileCover(const Point<double>& tl,
                                        const Point<double>& bl,
                                        const Point<double>& c,
                                        int32_t z) {
-    const int32_t tiles = 1 << z;
+    const int32_t tiles = (1 << z) + 1;
 
     struct ID {
         int32_t x, y;
@@ -96,30 +80,80 @@ std::vector<UnwrappedTileID> tileCover(const Point<double>& tl,
 
     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 - c.x, dy = y + 0.5 - c.y;
-                t.emplace_back(ID{ x, y, dx * dx + dy * dy });
-            }
+        for (x = x0; x < x1; ++x) {
+            const auto dx = x + 0.5 - c.x, dy = y + 0.5 - c.y;
+            t.emplace_back(ID { x, y, dx * dx + dy * dy });
         }
     };
 
-    // Divide the screen up in two triangles and scan each of them:
-    // \---+
-    // | \ |
-    // +---\.
-    scanTriangle(tl, tr, br, 0, tiles, scanLine);
-    scanTriangle(br, bl, tl, 0, tiles, scanLine);
+    std::vector<Point<double>> bounds = {tl, tr, br, bl};
+    while (bounds[0].y > min(min(bounds[1].y, bounds[2].y), bounds[3].y)) {
+        std::rotate(bounds.begin(), bounds.begin() + 1, bounds.end());
+    }
+    /*
+    Keeping the clockwise winding order (abcd), we rotated convex quadrilateral
+    angles in such way that angle a (bounds[0]) is on top):
+              a
+            /   \
+           /     b
+          /      |
+         /       c
+        /  ....
+       / ..
+      d
+    This is an example: we handle also cases where d.y < c.y, d.y < b.y etc.
+    Split the scan to tree steps:
+              a
+            /   \     (1)
+           /     b
+     -----------------
+          /      |    (2)
+         /       c
+     -----------------
+        /  ....
+       / ..           (3)
+      d
+    */
+    edge ab = edge(bounds[0], bounds[1]);
+    edge ad = edge(bounds[0], bounds[3]);
+
+    // Scan (1).
+    int32_t ymin = std::floor(bounds[0].y);
+    if (bounds[3].y < bounds[1].y) { std::swap(ab, ad); }
+    int32_t ymax = std::ceil(ab.y1);
+    if (ab.dy) {
+        scanSpans(ad, ab, std::max(0, ymin), std::min(tiles, ymax), scanLine);
+        ymin = ymax;
+    }
+
+    // Scan (2).
+    // yCutLower is c or d, whichever is with lower y value.
+    float yCutLower = min(bounds[2].y, ad.y1);
+    ymax = std::ceil(yCutLower);
+
+    // bc is edge opposite of ad.
+    edge bc = bounds[3].y < bounds[1].y ? edge(bounds[3], bounds[2]) : edge(bounds[1], bounds[2]);
+    if (bc.dy) {
+        scanSpans(ad, bc, std::max(0, ymin), std::min(tiles, ymax), scanLine);
+        ymin = ymax;
+    } else {
+        ymin = std::floor(yCutLower);
+    }
+
+    // Scan (3) - the triangle at the bottom.
+    if (ad.y1 < bc.y1) { std::swap(ad, bc); }
+    ymax = std::ceil(ad.y1);
+    bc = edge({ bc.x1, bc.y1 }, { ad.x1, ad.y1 });
+    if (bc.dy) { scanSpans(ad, bc, std::max(0, ymin), std::min(tiles, ymax), scanLine); }
 
     // Sort first by distance, then by x/y.
     std::sort(t.begin(), t.end(), [](const ID& a, const ID& b) {
         return std::tie(a.sqDist, a.x, a.y) < std::tie(b.sqDist, b.x, 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());
+    assert(t.end() == std::unique(t.begin(), t.end(), [](const ID& a, const ID& b) {
+        return a.x == b.x && a.y == b.y;
+    })); // no duplicates.
 
     std::vector<UnwrappedTileID> result;
     for (const auto& id : t) {

test/map/transform.test.cpp

@@ -9,6 +9,7 @@ using namespace mbgl;
 
 TEST(Transform, InvalidZoom) {
     Transform transform;
+    transform.resize({1, 1});
 
     ASSERT_DOUBLE_EQ(0, transform.getLatLng().latitude());
     ASSERT_DOUBLE_EQ(0, transform.getLatLng().longitude());
@@ -56,6 +57,7 @@ TEST(Transform, InvalidZoom) {
 
 TEST(Transform, InvalidBearing) {
     Transform transform;
+    transform.resize({1, 1});
 
     ASSERT_DOUBLE_EQ(0, transform.getLatLng().latitude());
     ASSERT_DOUBLE_EQ(0, transform.getLatLng().longitude());
@@ -78,6 +80,7 @@ TEST(Transform, InvalidBearing) {
 
 TEST(Transform, IntegerZoom) {
     Transform transform;
+    transform.resize({1, 1});
 
     auto checkIntegerZoom = [&transform](uint8_t zoomInt, double zoom) {
         transform.jumpTo(CameraOptions().withZoom(zoom));

test/util/tile_cover.test.cpp

@@ -43,6 +43,41 @@ TEST(TileCover, Pitch) {
               util::tileCover(transform.getState(), 2));
 }
 
+TEST(TileCover, PitchOverAllowedByContentInsets) {
+    Transform transform;
+    transform.resize({ 512, 512 });
+
+    transform.jumpTo(CameraOptions().withCenter(LatLng { 0.1, -0.1 }).withPadding(EdgeInsets { 376, 0, 0, 0 })
+                                    .withZoom(8.0).withBearing(45.0).withPitch(60.0));
+    // Top padding of 376 leads to capped pitch. See Transform::getMaxPitchForEdgeInsets.
+    EXPECT_LE(transform.getPitch() + 0.001, util::DEG2RAD * 60);
+
+    EXPECT_EQ((std::vector<UnwrappedTileID>{
+        { 3, 4, 3 }, { 3, 3, 3 }, { 3, 4, 4 }, { 3, 3, 4 }, { 3, 4, 2 }, { 3, 5, 3 }, { 3, 5, 2 }
+    }),
+              util::tileCover(transform.getState(), 3));
+}
+
+TEST(TileCover, PitchWithLargerResultSet) {
+    Transform transform;
+    transform.resize({ 1024, 768 });
+
+    // The values used here triggered the regression with left and right edge
+    // selection in tile_cover.cpp scanSpans.
+    transform.jumpTo(CameraOptions().withCenter(LatLng { 0.1, -0.1 }).withPadding(EdgeInsets { 400, 0, 0, 0 })
+                                    .withZoom(5).withBearing(-142.2630000003529176).withPitch(60.0));
+
+    auto cover = util::tileCover(transform.getState(), 5);
+    // Returned vector has above 100 elements, we check first 16 as there is a
+    // plan to return lower LOD for distant tiles.
+    EXPECT_EQ((std::vector<UnwrappedTileID> {
+        { 5, 15, 16 }, { 5, 15, 17 }, { 5, 14, 16 }, { 5, 14, 17 },
+        { 5, 16, 16 }, { 5, 16, 17 }, { 5, 15, 15 }, { 5, 14, 15 },
+        { 5, 15, 18 }, { 5, 14, 18 }, { 5, 16, 15 }, { 5, 13, 16 },
+        { 5, 13, 17 }, { 5, 16, 18 }, { 5, 13, 18 }, { 5, 15, 19 }
+    }), (std::vector<UnwrappedTileID> { cover.begin(), cover.begin() + 16}) );
+}
+
 TEST(TileCover, WorldZ1) {
     EXPECT_EQ((std::vector<UnwrappedTileID>{
         { 1, 0, 0 }, { 1, 0, 1 }, { 1, 1, 0 }, { 1, 1, 1 },