Source: mercator-gl.js

///////////////////////////////////////////////////////////////////////////////////
// The MIT License (MIT)
//
// Copyright (c) 2019 Tarek Sherif
//
// Permission is hereby granted, free of charge, to any person obtaining a copy of
// this software and associated documentation files (the "Software"), to deal in
// the Software without restriction, including without limitation the rights to
// use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
// the Software, and to permit persons to whom the Software is furnished to do so,
// subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
// COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
// IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
// CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
///////////////////////////////////////////////////////////////////////////////////

/**
    @module JavaScript
*/

// GLSL projection code from deck.gl https://github.com/uber/deck.gl
// Used under MIT licence

// transformMat4 from gl-matrix https://github.com/toji/gl-matrix/blob/master/src/vec4.js
// Used under MIT licence

import {PROJECTION_GLSL} from "./projection-glsl.js";

const PI = Math.PI;
const PI_4 = PI / 4;
const DEGREES_TO_RADIANS = PI / 180;
const TILE_SIZE = 512;
const EARTH_CIRCUMFERENCE = 40.03e6;


// High-precision for intermediate calculations
let tempCenter64 = new Float64Array(4);

// Low-precision for uniforms and to avoid instability
let lngLat32 = new Float32Array(2);

/**
    Create a Float32Array containing the low-part of 64-bit floats. The
    low parts can be passed to GLSL projection functions to improve numerical
    stability.

    @function
    @param {Float32Array} lnglat Array containing longitude latitude data (assumed to be in adjacent elements).
    @param {number} offset Offset of the first longitude value from the start of the array.
    @param {number} stride Number of elements between each pair of longitude/latitude values.
    @return {Float32Array} Array containing the low parts of input longitude/latitude data.
*/
export function highPrecisionLngLat(lngLat, offset = 0, stride = 2) {
    let numElements = Math.ceil((lngLat.length - offset) / stride);
    let precisionData = new Float32Array(numElements * 2);
    for (let i = 0; i < numElements; ++i) {
        let lli = offset + i * stride;
        let pi = i * 2;

        precisionData[pi]     = lngLat[lli]     - Math.fround(lngLat[lli]);
        precisionData[pi + 1] = lngLat[lli + 1] - Math.fround(lngLat[lli + 1]);
    }

    return precisionData;
}

/**
    Inject MercatorGL projection functions in to GLSL vertex shader source string.

    @function
    @param {string} vsSource Vertex shader source code. Can be GLSL 1 or 3.
    @return {string} Source code with the projection functions injected.
*/
export function injectMercatorGLSL(vsSource) {
    let versionMatch = vsSource.match(/#version \d+(\s+es)?\s*\n/);
    let versionLine = versionMatch ? versionMatch[0] : "";

    return vsSource.replace(versionLine, versionLine + PROJECTION_GLSL);
}

/**
    Allocate arrays to store projection uniform values. Will add uniforms
    to existing object, if provided.

    @function
    @param {object} [uniforms] Optional object containing existing uniforms.
    @return {object} Object containing projection uniforms.
*/
export function allocateMercatorUniforms(uniforms = {}) {
    uniforms.mercator_gl_lngLatCenter = new Float32Array(2);
    uniforms.mercator_gl_angleDerivatives = new Float32Array(3);
    uniforms.mercator_gl_clipCenter = new Float32Array(4);
    uniforms.mercator_gl_viewProjectionMatrix = new Float32Array(16);
    uniforms.mercator_gl_scale = 0;

    return uniforms;
}

/**
    Update projection uniforms based on current mercator settings and provided
    view-projection matrix. The view projection matrix is expected to project from
    a 512x512x<elevation meters> mercator space to clip space.

    @function
    @param {object} uniforms Object containing existing uniforms.
    @param {array} lngLat Current camera center as a 2-element longitude/latitude array.
    @param {number} zoom Current mercator zoom level.
    @param {array} viewProjectionMatrix 4x4 matrix to transform from mercator space to clip space.
    @return {object} Object containing updated uniforms.
*/
export function updateMercatorUniforms(uniforms, lngLat, zoom, viewProjectionMatrix) {
    let longitude = lngLat[0];
    let latitude = lngLat[1];

    uniforms.mercator_gl_scale = Math.pow(2, zoom);

    uniforms.mercator_gl_lngLatCenter[0] = longitude;
    uniforms.mercator_gl_lngLatCenter[1] = latitude;

    angleDerivatives(uniforms.mercator_gl_angleDerivatives, latitude, zoom);

    lngLat32[0] = longitude;
    lngLat32[1] = latitude;
    lngLatToClip(uniforms.mercator_gl_clipCenter, lngLat32, zoom, viewProjectionMatrix);

    uniforms.mercator_gl_viewProjectionMatrix.set(viewProjectionMatrix);
 
    return uniforms;
}

/**
    Mercator pixels/meter at the current latitude and zoom level (using 512x512 tile size).

    @function
    @param {number} latitude Current latitude of the camera center.
    @param {number} zoom Current mercator zoom level.
    @return {number} Number of pixels per meter.
*/
export function pixelsPerMeter(latitude, zoom) {
    let scale = Math.pow(2, zoom);
    let latCosine = Math.cos(latitude * DEGREES_TO_RADIANS);
    
    // Number of pixels occupied by one meter around current lat/lon
    return scale * TILE_SIZE / EARTH_CIRCUMFERENCE / latCosine;
}

/**
    Mercator pixels/degree of longitude and latitude at the current latitude and zoom level (using 512x512 tile size).

    @function
    @param {array} out Two-element array into which the results will be written.
    @param {number} latitude Current latitude of the camera center.
    @param {number} zoom Current mercator zoom level.
    @return {array} The "out" array provided, containing the results.
*/
export function pixelsPerDegree(out, latitude, zoom) {
    let scale = Math.pow(2, zoom);
    let latCosine = Math.cos(latitude * DEGREES_TO_RADIANS);
    
    // Number of pixels occupied by one degree around current lat/lon
    out[0] = scale * TILE_SIZE / 360;
    out[1] = out[0] / latCosine;

    return out;
}

/**
    Transform a longitude/latitude point to a 4d homegeneous coordinate in Mercator space.

    @function
    @param {array} out 4-element array into which the results will be written.
    @param {array} lngLat Longitude/latitude coordinate.
    @param {number} zoom Current mercator zoom level.
    @return {array} The "out" array provided, containing the results.
*/
export function lngLatToMercator(out, lngLat, zoom) {
    let longitude = lngLat[0];
    let latitude = lngLat[1];
    let scale = Math.pow(2, zoom);

    let lambda2 = longitude * DEGREES_TO_RADIANS;
    let phi2 = latitude * DEGREES_TO_RADIANS;
    let x = scale * TILE_SIZE * (lambda2 + PI) / (2 * PI);
    let y = scale * TILE_SIZE * (PI - Math.log(Math.tan(PI_4 + phi2 * 0.5))) / (2 * PI);

    out[0] = x;
    out[1] = y;
    out[2] = 0;
    out[3] = 1;

    return out;
}

/**
    Transform a homegeneous 4d representation of a mercator point on the plane
    (x, y, 0, 1) to clip space.

    @function
    @param {array} out 4-element array into which the results will be written.
    @param {array} mercatorPosition Homegeneous 4d coordinate on the mercator plane.
    @param {array} viewProjectionMatrix 4x4 matrix to transform from mercator space to clip space.
    @return {array} The "out" array provided, containing the results.
*/
export function mercatorToClip(out, mercatorPosition, viewProjectionMatrix) {
    transformMat4(out, mercatorPosition, viewProjectionMatrix);

    return out;
}

/**
    Transform a 2d longitude/latitude point to clip space to clip space.

    @function
    @param {array} out 4-element array into which the results will be written.
    @param {array} lngLat Longitude/latitude coordinate.
    @param {array} viewProjectionMatrix 4x4 matrix to transform from mercator space to clip space.
    @return {array} The "out" array provided, containing the results.
*/
export function lngLatToClip(out, lngLat, zoom, viewProjectionMatrix) {
    lngLatToMercator(tempCenter64, lngLat, zoom);
    mercatorToClip(out, tempCenter64, viewProjectionMatrix);

    return out;
}

function angleDerivatives(out, latitude, zoom, latCosine2) {
    let latCosine = Math.cos(latitude * DEGREES_TO_RADIANS);
    let latCosine2 = DEGREES_TO_RADIANS * Math.tan(latitude * DEGREES_TO_RADIANS) / latCosine;

    pixelsPerDegree(out, latitude, zoom);
    out[2] = out[0] * latCosine2 / 2;
}

function transformMat4(out, a, m) {
    let x = a[0];
    let y = a[1];
    let z = a[2];
    let w = a[3];
    out[0] = m[0] * x + m[4] * y + m[8] * z + m[12] * w;
    out[1] = m[1] * x + m[5] * y + m[9] * z + m[13] * w;
    out[2] = m[2] * x + m[6] * y + m[10] * z + m[14] * w;
    out[3] = m[3] * x + m[7] * y + m[11] * z + m[15] * w;
    
    return out;
}