From c1a18e1eac21a28e109b36925dfbf600c27939eb Mon Sep 17 00:00:00 2001 From: zmiao Date: Wed, 2 Jun 2021 12:29:25 +0300 Subject: [PATCH] Add slim version of cheap-ruler impl --- .../expression/definitions/distance.js | 5 +- src/style-spec/util/geometry_util.js | 177 +++++++++++++++++- 2 files changed, 178 insertions(+), 4 deletions(-) diff --git a/src/style-spec/expression/definitions/distance.js b/src/style-spec/expression/definitions/distance.js index 0c912976e04..4cc59de41b1 100644 --- a/src/style-spec/expression/definitions/distance.js +++ b/src/style-spec/expression/definitions/distance.js @@ -15,9 +15,8 @@ import type { GeoJSONPolygon, GeoJSONMultiPolygon } from "@mapbox/geojson-types"; -import {classifyRings, updateBBox, boxWithinBox, pointWithinPolygon, segmentIntersectSegment} from '../../util/geometry_util.js'; -import CheapRuler from "cheap-ruler"; -import type {GeometryPoint} from "cheap-ruler"; +import {classifyRings, updateBBox, boxWithinBox, pointWithinPolygon, segmentIntersectSegment, CheapRuler} from '../../util/geometry_util.js'; +import type {GeometryPoint} from '../../util/geometry_util.js'; import Point from "@mapbox/point-geometry"; import Queue from "tinyqueue"; diff --git a/src/style-spec/util/geometry_util.js b/src/style-spec/util/geometry_util.js index ca568d8cd47..5698b378404 100644 --- a/src/style-spec/util/geometry_util.js +++ b/src/style-spec/util/geometry_util.js @@ -3,7 +3,8 @@ import quickselect from 'quickselect'; import type Point from '@mapbox/point-geometry'; -import type {GeometryPoint} from "cheap-ruler"; + +export type GeometryPoint = [number, number] | [number, number, number]; /** * Returns the signed area for the polygon ring. Postive areas are exterior rings and * have a clockwise winding. Negative areas are interior rings and have a counter clockwise @@ -145,3 +146,177 @@ export function segmentIntersectSegment(a: GeometryPoint, b: GeometryPoint, c: G if (twoSided(a, b, c, d) && twoSided(c, d, a, b)) return true; return false; } + +// normalize a degree value into [-180..180] range +function wrap(deg) { + while (deg < -180) deg += 360; + while (deg > 180) deg -= 360; + return deg; +} + +const factors = { + kilometers: 1, + miles: 1000 / 1609.344, + nauticalmiles: 1000 / 1852, + meters: 1000, + metres: 1000, + yards: 1000 / 0.9144, + feet: 1000 / 0.3048, + inches: 1000 / 0.0254 +}; + +// Values that define WGS84 ellipsoid model of the Earth +const RE = 6378.137; // equatorial radius +const FE = 1 / 298.257223563; // flattening + +const E2 = FE * (2 - FE); +const RAD = Math.PI / 180; + +/** + * A collection of very fast approximations to common geodesic measurements. Useful for performance-sensitive code that measures things on a city scale. + * This is the slim version of https://github.com/mapbox/cheap-ruler + * @param {number} lat latitude + * @param {string} [units='kilometers'] + * @returns {CheapRuler} + * @example + * const ruler = cheapRuler(35.05, 'miles'); + * //=ruler + */ +export class CheapRuler { + kx: number; + ky: number; + /** + * Creates a ruler instance for very fast approximations to common geodesic measurements around a certain latitude. + * + * @param {number} lat latitude + * @param {string} [units='kilometers'] + * @returns {CheapRuler} + * @example + * const ruler = cheapRuler(35.05, 'miles'); + * //=ruler + */ + constructor(lat: number, units: string) { + if (lat === undefined) throw new Error('No latitude given.'); + if (units && !factors[units]) throw new Error(`Unknown unit ${units}. Use one of: ${Object.keys(factors).join(', ')}`); + + // Curvature formulas from https://en.wikipedia.org/wiki/Earth_radius#Meridional + const m = RAD * RE * (units ? factors[units] : 1); + const coslat = Math.cos(lat * RAD); + const w2 = 1 / (1 - E2 * (1 - coslat * coslat)); + const w = Math.sqrt(w2); + + // multipliers for converting longitude and latitude degrees into distance + this.kx = m * w * coslat; // based on normal radius of curvature + this.ky = m * w * w2 * (1 - E2); // based on meridonal radius of curvature + } + + /** + * Given two points of the form [longitude, latitude], returns the distance. + * + * @param {GeometryPoint} a point [longitude, latitude] + * @param {GeometryPoint} b point [longitude, latitude] + * @returns {number} distance + * @example + * const distance = ruler.distance([30.5, 50.5], [30.51, 50.49]); + * //=distance + */ + distance(a: GeometryPoint, b: GeometryPoint): number { + const dx = wrap(a[0] - b[0]) * this.kx; + const dy = (a[1] - b[1]) * this.ky; + return Math.sqrt(dx * dx + dy * dy); + } + + /** + * Returns the distance from a point `p` to a line segment `a` to `b`. + * + * @param {GeometryPoint} p point [longitude, latitude] + * @param {GeometryPoint} a segment point 1 [longitude, latitude] + * @param {GeometryPoint} b segment point 2 [longitude, latitude] + * @returns {number} distance + * @example + * const distance = ruler.pointToSegmentDistance([-67.04, 50.5], [-67.05, 50.57], [-67.03, 50.54]); + * //=distance + */ + pointToSegmentDistance(p: GeometryPoint, a: GeometryPoint, b: GeometryPoint): number { + let [x, y] = a; + let dx = wrap(b[0] - x) * this.kx; + let dy = (b[1] - y) * this.ky; + let t = 0; + + if (dx !== 0 || dy !== 0) { + t = (wrap(p[0] - x) * this.kx * dx + (p[1] - y) * this.ky * dy) / (dx * dx + dy * dy); + + if (t > 1) { + x = b[0]; + y = b[1]; + + } else if (t > 0) { + x += (dx / this.kx) * t; + y += (dy / this.ky) * t; + } + } + + dx = wrap(p[0] - x) * this.kx; + dy = (p[1] - y) * this.ky; + + return Math.sqrt(dx * dx + dy * dy); + } + + /** + * Returns an object of the form {point, index, t}, where point is closest point on the line + * from the given point, index is the start index of the segment with the closest point, + * and t is a parameter from 0 to 1 that indicates where the closest point is on that segment. + * + * @param {Array} line + * @param {GeometryPoint} p point [longitude, latitude] + * @returns {Object} {point, index, t} + * @example + * const point = ruler.pointOnLine(line, [-67.04, 50.5]).point; + * //=point + */ + pointOnLine(line: Array, p: GeometryPoint): Object { + let minDist = Infinity; + let minX = Infinity, minY = Infinity, minI = Infinity, minT = Infinity; + + for (let i = 0; i < line.length - 1; i++) { + + let x = line[i][0]; + let y = line[i][1]; + let dx = wrap(line[i + 1][0] - x) * this.kx; + let dy = (line[i + 1][1] - y) * this.ky; + let t = 0; + + if (dx !== 0 || dy !== 0) { + t = (wrap(p[0] - x) * this.kx * dx + (p[1] - y) * this.ky * dy) / (dx * dx + dy * dy); + + if (t > 1) { + x = line[i + 1][0]; + y = line[i + 1][1]; + + } else if (t > 0) { + x += (dx / this.kx) * t; + y += (dy / this.ky) * t; + } + } + + dx = wrap(p[0] - x) * this.kx; + dy = (p[1] - y) * this.ky; + + const sqDist = dx * dx + dy * dy; + if (sqDist < minDist) { + minDist = sqDist; + minX = x; + minY = y; + minI = i; + minT = t; + } + } + + return { + point: [minX, minY], + index: minI, + t: Math.max(0, Math.min(1, minT)) + }; + } + +}