LOD support for tile coverage (#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
+};