Fix constraining logic for alternate projections

debug/index.html

@@ -20,10 +20,11 @@
 
 var map = window.map = new mapboxgl.Map({
     container: 'map',
-    zoom: 12.5,
-    center: [-122.4194, 37.7749],
+    zoom: 3,
+    center: [-100, 40],
     style: 'mapbox://styles/mapbox/streets-v11',
-    hash: true
+    hash: true,
+    // projection: 'albers'
 });
 
 </script>

src/data/load_geometry.js

@@ -4,7 +4,9 @@ import {warnOnce, clamp} from '../util/util.js';
 
 import EXTENT from './extent.js';
 import {lngFromMercatorX, latFromMercatorY} from '../geo/mercator_coordinate.js';
+import resample from '../geo/projection/resample.js';
 import Point from '@mapbox/point-geometry';
+
 import type {CanonicalTileID} from '../source/tile_id.js';
 import type {TileTransform} from '../geo/projection/tile_transform.js';
 
@@ -28,12 +30,6 @@ function clampPoint(point: Point) {
     return point;
 }
 
-function pointToLineDist(px, py, ax, ay, bx, by) {
-    const dx = ax - bx;
-    const dy = ay - by;
-    return Math.abs((ay - py) * dx - (ax - px) * dy) / Math.hypot(dx, dy);
-}
-
 // a subset of VectorTileGeometry
 type FeatureWithGeometry = {
     extent: number;
@@ -51,64 +47,29 @@ export default function loadGeometry(feature: FeatureWithGeometry, canonical?: C
     const featureExtent = feature.extent;
     const scale = EXTENT / featureExtent;
     const projection = tileTransform ? tileTransform.projection : undefined;
-    let z2;
     const isMercator = !projection || projection.name === 'mercator';
-    if (canonical && !isMercator) {
-        z2 = Math.pow(2, canonical.z);
-    }
 
     function reproject(p) {
         if (isMercator || !canonical || !tileTransform || !projection) {
-            return new Point(Math.round(p.x * scale), Math.round(p.y * scale));
+            return clampPoint(new Point(Math.round(p.x * scale), Math.round(p.y * scale)));
         } else {
+            const z2 = 1 << canonical.z;
             const lng = lngFromMercatorX((canonical.x + p.x / featureExtent) / z2);
             const lat = latFromMercatorY((canonical.y + p.y / featureExtent) / z2);
             const {x, y} = projection.project(lng, lat);
-            return new Point(
+            return clampPoint(new Point(
                 Math.round((x * tileTransform.scale - tileTransform.x) * EXTENT),
                 Math.round((y * tileTransform.scale - tileTransform.y) * EXTENT)
-            );
-        }
-    }
-
-    function addResampled(resampled, startMerc, endMerc, startProj, endProj) {
-        const midMerc = new Point(
-            (startMerc.x + endMerc.x) / 2,
-            (startMerc.y + endMerc.y) / 2);
-        const midProj = reproject(midMerc);
-        const err = pointToLineDist(midProj.x, midProj.y, startProj.x, startProj.y, endProj.x, endProj.y);
-
-        if (err >= 1) {
-            // we're very unlikely to hit max call stack exceeded here,
-            // but we might want to safeguard against it in the future
-            addResampled(resampled, startMerc, midMerc, startProj, midProj);
-            addResampled(resampled, midMerc, endMerc, midProj, endProj);
-        } else {
-            resampled.push(clampPoint(endProj));
+            ));
         }
     }
 
     const geometry = feature.loadGeometry();
 
-    for (let r = 0; r < geometry.length; r++) {
-        const ring = geometry[r];
-        const resampled = [];
-
-        for (let i = 0, prevMerc, prevProj; i < ring.length; i++) {
-            const pointMerc = ring[i];
-            const pointProj = reproject(ring[i]);
-
-            if (projection && projection.name !== 'mercator' && prevMerc && prevProj && feature.type !== 1) {
-                addResampled(resampled, prevMerc, pointMerc, prevProj, pointProj);
-            } else {
-                resampled.push(clampPoint(pointProj));
-            }
-
-            prevMerc = pointMerc;
-            prevProj = pointProj;
-        }
-
-        geometry[r] = resampled;
+    for (let i = 0; i < geometry.length; i++) {
+        geometry[i] = !isMercator && feature.type !== 1 ?
+            resample(geometry[i], reproject, 1) :
+            geometry[i].map(reproject);
     }
 
     return geometry;

src/geo/projection/adjustments.js

@@ -3,6 +3,7 @@
 import LngLat from '../lng_lat.js';
 import MercatorCoordinate from '../mercator_coordinate.js';
 import {mat4, mat2} from 'gl-matrix';
+import {clamp} from '../../util/util.js';
 import type {Projection} from './index.js';
 import type Transform from '../transform.js';
 
@@ -37,22 +38,30 @@ function getInterpolationT(transform: Transform) {
     const rangeAdjustment = Math.log(size / 1024) / Math.LN2;
     const zoomA = range[0] + rangeAdjustment;
     const zoomB = range[1] + rangeAdjustment;
-    const t = Math.max(0, Math.min(1, (transform.zoom - zoomA) / (zoomB - zoomA)));
-
+    const t = clamp((transform.zoom - zoomA) / (zoomB - zoomA), 0, 1);
     return t;
 }
 
+// approx. kilometers per longitude degree at equator
+const offset = 1 / 40000;
+
 /*
  * Calculates the scale difference between Mercator and the given projection at a certain location.
  */
 function getZoomAdjustment(projection: Projection, loc: LngLat) {
-    const loc2 = new LngLat(loc.lng + 360 / 40000, loc.lat);
-    const p1 = projection.project(loc.lng, loc.lat);
+    const loc1 = new LngLat(loc.lng - 180 * offset, loc.lat);
+    const loc2 = new LngLat(loc.lng + 180 * offset, loc.lat);
+
+    const p1 = projection.project(loc1.lng, loc1.lat);
     const p2 = projection.project(loc2.lng, loc2.lat);
 
-    const m1 = MercatorCoordinate.fromLngLat(loc);
+    const m1 = MercatorCoordinate.fromLngLat(loc1);
     const m2 = MercatorCoordinate.fromLngLat(loc2);
 
+    // make sure we operate within mercator space for adjustments (they can go over for other projections)
+    m1.y = clamp(m1.y, -1, 1);
+    m2.y = clamp(m2.y, -1, 1);
+
     const pdx = p2.x - p1.x;
     const pdy = p2.y - p1.y;
     const mdx = m2.x - m1.x;
@@ -65,41 +74,47 @@ function getZoomAdjustment(projection: Projection, loc: LngLat) {
 
 function getShearAdjustment(projection, zoom, loc, interpT, withoutRotation?: boolean) {
 
-    // create a location a tiny amount (~1km) east of the given location
-    const loc2 = new LngLat(loc.lng + 360 / 40000, loc.lat);
+    // create two locations a tiny amount (~1km) east and west of the given location
+    const locw = new LngLat(loc.lng - 180 * offset, loc.lat);
+    const loce = new LngLat(loc.lng + 180 * offset, loc.lat);
 
-    const p1 = projection.project(loc.lng, loc.lat);
-    const p2 = projection.project(loc2.lng, loc2.lat);
+    const pw = projection.project(locw.lng, locw.lat);
+    const pe = projection.project(loce.lng, loce.lat);
 
-    const pdx = p2.x - p1.x;
-    const pdy = p2.y - p1.y;
+    const pdx = pe.x - pw.x;
+    const pdy = pe.y - pw.y;
 
     // Calculate how much the map would need to be rotated to make east-west in
     // projected coordinates be left-right
     const angleAdjust = -Math.atan(pdy / pdx) / Math.PI * 180;
 
-    // Find the projected coordinates of two locations, one slightly north and one slightly east.
+    // Pick a location identical to the original one except for poles to make sure we're within mercator bounds
+    const mc2 = MercatorCoordinate.fromLngLat(loc);
+    mc2.y = clamp(mc2.y, -1 + offset, 1 - offset);
+    const loc2 = mc2.toLngLat();
+    const p2 = projection.project(loc2.lng, loc2.lat);
+
+    // Find the projected coordinates of two locations, one slightly south and one slightly east.
     // Then calculate the transform that would make the projected coordinates of the two locations be:
     // - equal distances from the original location
     // - perpendicular to one another
     //
     // Only the position of the coordinate to the north is adjusted.
     // The coordinate to the east stays where it is.
-    const mc3 = MercatorCoordinate.fromLngLat(loc);
-    const offset = 1 / 40000;
+    const mc3 = MercatorCoordinate.fromLngLat(loc2);
     mc3.x += offset;
     const loc3 = mc3.toLngLat();
     const p3 = projection.project(loc3.lng, loc3.lat);
-    const pdx3 = p3.x - p1.x;
-    const pdy3 = p3.y - p1.y;
+    const pdx3 = p3.x - p2.x;
+    const pdy3 = p3.y - p2.y;
     const delta3 = rotate(pdx3, pdy3, angleAdjust);
 
-    const mc4 = MercatorCoordinate.fromLngLat(loc);
+    const mc4 = MercatorCoordinate.fromLngLat(loc2);
     mc4.y += offset;
     const loc4 = mc4.toLngLat();
     const p4 = projection.project(loc4.lng, loc4.lat);
-    const pdx4 = p4.x - p1.x;
-    const pdy4 = p4.y - p1.y;
+    const pdx4 = p4.x - p2.x;
+    const pdy4 = p4.y - p2.y;
     const delta4 = rotate(pdx4, pdy4, angleAdjust);
 
     const scale = Math.abs(delta3.x) / Math.abs(delta4.y);

src/geo/projection/resample.js

@@ -0,0 +1,49 @@
+// @flow
+
+import Point from '@mapbox/point-geometry';
+
+function pointToLineDist(px, py, ax, ay, bx, by) {
+    const dx = ax - bx;
+    const dy = ay - by;
+    return Math.abs((ay - py) * dx - (ax - px) * dy) / Math.hypot(dx, dy);
+}
+
+function addResampled(resampled, startMerc, endMerc, startProj, endProj, reproject, tolerance) {
+    const midMerc = new Point(
+        (startMerc.x + endMerc.x) / 2,
+        (startMerc.y + endMerc.y) / 2);
+
+    const midProj = reproject(midMerc);
+    const err = pointToLineDist(midProj.x, midProj.y, startProj.x, startProj.y, endProj.x, endProj.y);
+
+    // if reprojected midPoint is too far from geometric midpoint, recurse into two halves
+    if (err >= tolerance) {
+        // we're very unlikely to hit max call stack exceeded here,
+        // but we might want to safeguard against it in the future
+        addResampled(resampled, startMerc, midMerc, startProj, midProj, reproject, tolerance);
+        addResampled(resampled, midMerc, endMerc, midProj, endProj, reproject, tolerance);
+
+    } else { // otherwise, just add the point
+        resampled.push(endProj);
+    }
+}
+
+export default function resample(line: Array<Point>, reproject: (Point) => Point, tolerance: number): Array<Point> {
+    const resampled = [];
+    let prevMerc, prevProj;
+
+    for (const pointMerc of line) {
+        const pointProj = reproject(pointMerc);
+
+        if (prevMerc && prevProj) {
+            addResampled(resampled, prevMerc, pointMerc, prevProj, pointProj, reproject, tolerance);
+        } else {
+            resampled.push(pointProj);
+        }
+
+        prevMerc = pointMerc;
+        prevProj = pointProj;
+    }
+
+    return resampled;
+}

src/geo/projection/tile_transform.js

@@ -1,6 +1,5 @@
 // @flow
 import {lngFromMercatorX, latFromMercatorY} from '../mercator_coordinate.js';
-import {number as interpolate} from '../../style-spec/util/interpolate.js';
 import type {Projection} from './index.js';
 
 export type TileTransform = {
@@ -14,37 +13,66 @@ export type TileTransform = {
 
 export default function tileTransform(id: Object, projection: Projection) {
     const s = Math.pow(2, -id.z);
+
     const x1 = (id.x) * s;
     const x2 = (id.x + 1) * s;
     const y1 = (id.y) * s;
     const y2 = (id.y + 1) * s;
+
+    if (projection.name === 'mercator') {
+        return {scale: 1 << id.z, x: 0, y: 0, x2: 1, y2: 1, projection};
+    }
+
     const lng1 = lngFromMercatorX(x1);
     const lng2 = lngFromMercatorX(x2);
     const lat1 = latFromMercatorY(y1);
     const lat2 = latFromMercatorY(y2);
-    let minX = Infinity;
-    let minY = Infinity;
-    let maxX = -Infinity;
-    let maxY = -Infinity;
-
-    const n = 2;
-    for (let i = 0; i <= n; i++) {
-        const lng = lngFromMercatorX(interpolate(x1, x2, i / n));
-        const lat = latFromMercatorY(interpolate(y1, y2, i / n));
-
-        const p0 = projection.project(lng, lat1);
-        const p1 = projection.project(lng, lat2);
-        const p2 = projection.project(lng1, lat);
-        const p3 = projection.project(lng2, lat);
-
-        minX = Math.min(minX, p0.x, p1.x, p2.x, p3.x);
-        maxX = Math.max(maxX, p0.x, p1.x, p2.x, p3.x);
-        minY = Math.min(minY, p0.y, p1.y, p2.y, p3.y);
-        maxY = Math.max(maxY, p0.y, p1.y, p2.y, p3.y);
+
+    const p0 = projection.project(lng1, lat1);
+    const p1 = projection.project(lng2, lat1);
+    const p2 = projection.project(lng2, lat2);
+    const p3 = projection.project(lng1, lat2);
+
+    let minX = Math.min(p0.x, p1.x, p2.x, p3.x);
+    let minY = Math.min(p0.y, p1.y, p2.y, p3.y);
+    let maxX = Math.max(p0.x, p1.x, p2.x, p3.x);
+    let maxY = Math.max(p0.y, p1.y, p2.y, p3.y);
+
+    // we pick an error threshold for calculating the bbox that balances between performance and precision
+    const maxErr = s / 16;
+
+    function processSegment(pa, pb, ax, ay, bx, by) {
+        const mx = (ax + bx) / 2;
+        const my = (ay + by) / 2;
+
+        const pm = projection.project(lngFromMercatorX(mx), latFromMercatorY(my));
+        const err = Math.max(0, minX - pm.x, minY - pm.y, pm.x - maxX, pm.y - maxY);
+
+        minX = Math.min(minX, pm.x);
+        maxX = Math.max(maxX, pm.x);
+        minY = Math.min(minY, pm.y);
+        maxY = Math.max(maxY, pm.y);
+
+        if (err > maxErr) {
+            processSegment(pa, pm, ax, ay, mx, my);
+            processSegment(pm, pb, mx, my, bx, by);
+        }
     }
 
+    processSegment(p0, p1, x1, y1, x2, y1);
+    processSegment(p1, p2, x2, y1, x2, y2);
+    processSegment(p2, p3, x2, y2, x1, y2);
+    processSegment(p3, p0, x1, y2, x1, y1);
+
+    // extend the bbox by max error to make sure coords don't go past tile extent
+    minX -= maxErr;
+    minY -= maxErr;
+    maxX += maxErr;
+    maxY += maxErr;
+
     const max = Math.max(maxX - minX, maxY - minY);
     const scale = 1 / max;
+
     return {
         scale,
         x: minX * scale,