LOD support for tile coverage (mapbox#8975)

* LOD support for tile coverage

src/geo/transform.js

@@ -6,12 +6,11 @@ import MercatorCoordinate, {mercatorXfromLng, mercatorYfromLat, mercatorZfromAlt
 import Point from '@mapbox/point-geometry';
 import {wrap, clamp} from '../util/util';
 import {number as interpolate} from '../style-spec/util/interpolate';
-import tileCover from '../util/tile_cover';
-import {UnwrappedTileID} from '../source/tile_id';
 import EXTENT from '../data/extent';
-import {vec4, mat4, mat2} from 'gl-matrix';
+import {vec4, mat4, mat2, vec2} from 'gl-matrix';
+import {Aabb, Frustum} from '../util/primitives.js';
 
-import type {OverscaledTileID, CanonicalTileID} from '../source/tile_id';
+import {UnwrappedTileID, OverscaledTileID, CanonicalTileID} from '../source/tile_id';
 
 /**
  * A single transform, generally used for a single tile to be
@@ -34,6 +33,7 @@ class Transform {
     cameraToCenterDistance: number;
     mercatorMatrix: Array<number>;
     projMatrix: Float64Array;
+    invProjMatrix: Float64Array;
     alignedProjMatrix: Float64Array;
     pixelMatrix: Float64Array;
     pixelMatrixInverse: Float64Array;
@@ -278,15 +278,90 @@ class Transform {
 
         const centerCoord = MercatorCoordinate.fromLngLat(this.center);
         const numTiles = Math.pow(2, z);
-        const centerPoint = new Point(numTiles * centerCoord.x - 0.5, numTiles * centerCoord.y - 0.5);
-        const cornerCoords = [
-            this.pointCoordinate(new Point(0, 0)),
-            this.pointCoordinate(new Point(this.width, 0)),
-            this.pointCoordinate(new Point(this.width, this.height)),
-            this.pointCoordinate(new Point(0, this.height))
-        ];
-        return tileCover(z, cornerCoords, options.reparseOverscaled ? actualZ : z, this._renderWorldCopies)
-            .sort((a, b) => centerPoint.dist(a.canonical) - centerPoint.dist(b.canonical));
+        const centerPoint = [numTiles * centerCoord.x, numTiles * centerCoord.y, 0];
+        const cameraFrustum = Frustum.fromInvProjectionMatrix(this.invProjMatrix, this.worldSize, z);
+
+        // No change of LOD behavior for pitch lower than 60: return only tile ids from the requested zoom level
+        let minZoom = options.minzoom || 0;
+        if (this.pitch <= 60.0)
+            minZoom = z;
+
+        // There should always be a certain number of maximum zoom level tiles surrounding the center location
+        const radiusOfMaxLvlLodInTiles = 3;
+
+        const newRootTile = (wrap: number): any => {
+            return {
+                // All tiles are on zero elevation plane => z difference is zero
+                aabb: new Aabb([wrap * numTiles, 0, 0], [(wrap + 1) * numTiles, numTiles, 0]),
+                zoom: 0,
+                x: 0,
+                y: 0,
+                wrap,
+                fullyVisible: false
+            };
+        };
+
+        // Do a depth-first traversal to find visible tiles and proper levels of detail
+        const stack = [];
+        const result = [];
+        const maxZoom = z;
+        const overscaledZ = options.reparseOverscaled ? actualZ : z;
+
+        if (this._renderWorldCopies) {
+            // Render copy of the globe thrice on both sides
+            for (let i = 1; i <= 3; i++) {
+                stack.push(newRootTile(-i));
+                stack.push(newRootTile(i));
+            }
+        }
+
+        stack.push(newRootTile(0));
+
+        while (stack.length > 0) {
+            const it = stack.pop();
+            const x = it.x;
+            const y = it.y;
+            let fullyVisible = it.fullyVisible;
+
+            // Visibility of a tile is not required if any of its ancestor if fully inside the frustum
+            if (!fullyVisible) {
+                const intersectResult = it.aabb.intersects(cameraFrustum);
+
+                if (intersectResult === 0)
+                    continue;
+
+                fullyVisible = intersectResult === 2;
+            }
+
+            const distanceX = it.aabb.distanceX(centerPoint);
+            const distanceY = it.aabb.distanceY(centerPoint);
+            const longestDim = Math.max(Math.abs(distanceX), Math.abs(distanceY));
+
+            // We're using distance based heuristics to determine if a tile should be split into quadrants or not.
+            // radiusOfMaxLvlLodInTiles defines that there's always a certain number of maxLevel tiles next to the map center.
+            // Using the fact that a parent node in quadtree is twice the size of its children (per dimension)
+            // we can define distance thresholds for each relative level:
+            // f(k) = offset + 2 + 4 + 8 + 16 + ... + 2^k. This is the same as "offset+2^(k+1)-2"
+            const distToSplit = radiusOfMaxLvlLodInTiles + (1 << (maxZoom - it.zoom)) - 2;
+
+            // Have we reached the target depth or is the tile too far away to be any split further?
+            if (it.zoom === maxZoom || (longestDim > distToSplit && it.zoom >= minZoom)) {
+                result.push({
+                    tileID: new OverscaledTileID(it.zoom === maxZoom ? overscaledZ : it.zoom, it.wrap, it.zoom, x, y),
+                    distanceSq: vec2.sqrLen([centerPoint[0] - 0.5 - x, centerPoint[1] - 0.5 - y])
+                });
+                continue;
+            }
+
+            for (let i = 0; i < 4; i++) {
+                const childX = (x << 1) + (i % 2);
+                const childY = (y << 1) + (i >> 1);
+
+                stack.push({aabb: it.aabb.quadrant(i), zoom: it.zoom + 1, x: childX, y: childY, wrap: it.wrap, fullyVisible});
+            }
+        }
+
+        return result.sort((a, b) => a.distanceSq - b.distanceSq).map(a => a.tileID);
     }
 
     resize(width: number, height: number) {
@@ -549,7 +624,7 @@ class Transform {
         // (the distance between[width/2, height/2] and [width/2 + 1, height/2])
         const halfFov = this._fov / 2;
         const groundAngle = Math.PI / 2 + this._pitch;
-        const topHalfSurfaceDistance = Math.sin(halfFov) * this.cameraToCenterDistance / Math.sin(Math.PI - groundAngle - halfFov);
+        const topHalfSurfaceDistance = Math.sin(halfFov) * this.cameraToCenterDistance / Math.sin(clamp(Math.PI - groundAngle - halfFov, 0.01, Math.PI - 0.01));
         const point = this.point;
         const x = point.x, y = point.y;
 
@@ -585,6 +660,7 @@ class Transform {
         mat4.scale(m, m, [1, 1, mercatorZfromAltitude(1, this.center.lat) * this.worldSize, 1]);
 
         this.projMatrix = m;
+        this.invProjMatrix = mat4.invert([], this.projMatrix);
 
         // Make a second projection matrix that is aligned to a pixel grid for rendering raster tiles.
         // We're rounding the (floating point) x/y values to achieve to avoid rendering raster images to fractional

src/util/primitives.js

@@ -0,0 +1,145 @@
+// @flow
+
+import {vec3, vec4} from 'gl-matrix';
+import assert from 'assert';
+
+class Frustum {
+    points: Array<Array<number>>;
+    planes: Array<Array<number>>;
+
+    constructor(points_: Array<Array<number>>, planes_: Array<Array<number>>) {
+        this.points = points_;
+        this.planes = planes_;
+    }
+
+    static fromInvProjectionMatrix(invProj: Float64Array, worldSize: number, zoom: number): Frustum {
+        const clipSpaceCorners = [
+            [-1, 1, -1, 1],
+            [ 1, 1, -1, 1],
+            [ 1, -1, -1, 1],
+            [-1, -1, -1, 1],
+            [-1, 1, 1, 1],
+            [ 1, 1, 1, 1],
+            [ 1, -1, 1, 1],
+            [-1, -1, 1, 1]
+        ];
+
+        const scale = Math.pow(2, zoom);
+
+        // Transform frustum corner points from clip space to tile space
+        const frustumCoords = clipSpaceCorners
+            .map(v => vec4.transformMat4([], v, invProj))
+            .map(v => vec4.scale([], v, 1.0 / v[3] / worldSize * scale));
+
+        const frustumPlanePointIndices = [
+            [0, 1, 2],  // near
+            [6, 5, 4],  // far
+            [0, 3, 7],  // left
+            [2, 1, 5],  // right
+            [3, 2, 6],  // bottom
+            [0, 4, 5]   // top
+        ];
+
+        const frustumPlanes = frustumPlanePointIndices.map((p: Array<number>) => {
+            const a = vec3.sub([], frustumCoords[p[0]], frustumCoords[p[1]]);
+            const b = vec3.sub([], frustumCoords[p[2]], frustumCoords[p[1]]);
+            const n = vec3.normalize([], vec3.cross([], a, b));
+            const d = -vec3.dot(n, frustumCoords[p[1]]);
+            return n.concat(d);
+        });
+
+        return new Frustum(frustumCoords, frustumPlanes);
+    }
+}
+
+class Aabb {
+    min: vec3;
+    max: vec3;
+    center: vec3;
+
+    constructor(min_: vec3, max_: vec3) {
+        this.min = min_;
+        this.max = max_;
+        this.center = vec3.scale([], vec3.add([], this.min, this.max), 0.5);
+    }
+
+    quadrant(index: number): Aabb {
+        const split = [(index % 2) === 0, index < 2];
+        const qMin = vec3.clone(this.min);
+        const qMax = vec3.clone(this.max);
+        for (let axis = 0; axis < split.length; axis++) {
+            qMin[axis] = split[axis] ? this.min[axis] : this.center[axis];
+            qMax[axis] = split[axis] ? this.center[axis] : this.max[axis];
+        }
+        // Elevation is always constant, hence quadrant.max.z = this.max.z
+        qMax[2] = this.max[2];
+        return new Aabb(qMin, qMax);
+    }
+
+    distanceX(point: Array<number>): number {
+        const pointOnAabb = Math.max(Math.min(this.max[0], point[0]), this.min[0]);
+        return pointOnAabb - point[0];
+    }
+
+    distanceY(point: Array<number>): number {
+        const pointOnAabb = Math.max(Math.min(this.max[1], point[1]), this.min[1]);
+        return pointOnAabb - point[1];
+    }
+
+    // Performs a frustum-aabb intersection test. Returns 0 if there's no intersection,
+    // 1 if shapes are intersecting and 2 if the aabb if fully inside the frustum.
+    intersects(frustum: Frustum): number {
+        // Execute separating axis test between two convex objects to find intersections
+        // Each frustum plane together with 3 major axes define the separating axes
+        // Note: test only 4 points as both min and max points have equal elevation
+        assert(this.min[2] === 0 && this.max[2] === 0);
+
+        const aabbPoints = [
+            [this.min[0], this.min[1], 0.0, 1],
+            [this.max[0], this.min[1], 0.0, 1],
+            [this.max[0], this.max[1], 0.0, 1],
+            [this.min[0], this.max[1], 0.0, 1]
+        ];
+
+        let fullyInside = true;
+
+        for (let p = 0; p < frustum.planes.length; p++) {
+            const plane = frustum.planes[p];
+            let pointsInside = 0;
+
+            for (let i = 0; i < aabbPoints.length; i++) {
+                pointsInside += vec4.dot(plane, aabbPoints[i]) >= 0;
+            }
+
+            if (pointsInside === 0)
+                return 0;
+
+            if (pointsInside !== aabbPoints.length)
+                fullyInside = false;
+        }
+
+        if (fullyInside)
+            return 2;
+
+        for (let axis = 0; axis < 3; axis++) {
+            let projMin = Number.MAX_VALUE;
+            let projMax = -Number.MAX_VALUE;
+
+            for (let p = 0; p < frustum.points.length; p++) {
+                const projectedPoint = frustum.points[p][axis] - this.min[axis];
+
+                projMin = Math.min(projMin, projectedPoint);
+                projMax = Math.max(projMax, projectedPoint);
+            }
+
+            if (projMax < 0 || projMin > this.max[axis] - this.min[axis])
+                return 0;
+        }
+
+        return 1;
+    }
+}
+export {
+    Aabb,
+    Frustum
+};