diff --git a/examples/jsm/OGC3DTilesHelper.js b/examples/jsm/OGC3DTilesHelper.js index babf74596a..ce0306bad6 100644 --- a/examples/jsm/OGC3DTilesHelper.js +++ b/examples/jsm/OGC3DTilesHelper.js @@ -47,7 +47,7 @@ function zoomToSphere(view, tile, sphere) { const distance = radius * Math.tan(fov * 2); return { - coord: new Coordinates('EPSG:4978', center), + coord: new Coordinates('EPSG:4978').setFromVector3(center), range: distance + radius, }; } diff --git a/src/Controls/GlobeControls.js b/src/Controls/GlobeControls.js index 76b5901a2a..68ae25518d 100644 --- a/src/Controls/GlobeControls.js +++ b/src/Controls/GlobeControls.js @@ -755,7 +755,7 @@ class GlobeControls extends THREE.EventDispatcher { const range = this.getRange(point); if (point && range > this.minDistance) { return this.lookAtCoordinate({ - coord: new Coordinates('EPSG:4978', point), + coord: new Coordinates('EPSG:4978').setFromVector3(point), range: range * (event.direction === 'out' ? 1 / 0.6 : 0.6), time: 1500, }); @@ -1028,7 +1028,9 @@ class GlobeControls extends THREE.EventDispatcher { */ getCameraCoordinate() { - return new Coordinates('EPSG:4978', this.camera.position).as('EPSG:4326'); + return new Coordinates('EPSG:4978') + .setFromVector3(this.camera.position) + .as('EPSG:4326'); } /** @@ -1216,7 +1218,9 @@ class GlobeControls extends THREE.EventDispatcher { return; } - return new Coordinates('EPSG:4978', pickedPosition).as('EPSG:4326'); + return new Coordinates('EPSG:4978') + .setFromVector3(pickedPosition) + .as('EPSG:4326'); } } diff --git a/src/Core/Geographic/Coordinates.js b/src/Core/Geographic/Coordinates.js deleted file mode 100644 index d0934066ea..0000000000 --- a/src/Core/Geographic/Coordinates.js +++ /dev/null @@ -1,332 +0,0 @@ -import * as THREE from 'three'; -import proj4 from 'proj4'; -import * as CRS from 'Core/Geographic/Crs'; -import Ellipsoid from 'Core/Math/Ellipsoid'; - -proj4.defs('EPSG:4978', '+proj=geocent +datum=WGS84 +units=m +no_defs'); - -const ellipsoid = new Ellipsoid(); -const projectionCache = {}; - -const v0 = new THREE.Vector3(); -const v1 = new THREE.Vector3(); - -let coord0; -let coord1; - -function proj4cache(crsIn, crsOut) { - if (!projectionCache[crsIn]) { - projectionCache[crsIn] = {}; - } - - if (!projectionCache[crsIn][crsOut]) { - projectionCache[crsIn][crsOut] = proj4(crsIn, crsOut); - } - - return projectionCache[crsIn][crsOut]; -} - -/** - * A Coordinates object, defined by a [crs](http://inspire.ec.europa.eu/theme/rs) - * and three values. These values are accessible through `x`, `y` and `z`, - * although it can also be accessible through `latitude`, `longitude` and - * `altitude`. To change a value, prefer the `set()` method below. - * - * `EPSG:4978` and `EPSG:4326` are supported by default. To use another CRS, - * you have to declare it with `proj4`. You can find most projections and their - * proj4 code at [epsg.io](https://epsg.io/). - * - * @example - * new Coordinates('EPSG:4978', 20885167, 849862, 23385912); //Geocentric coordinates - * new Coordinates('EPSG:4326', 2.33, 48.24, 24999549); //Geographic coordinates - * - * @example - * // Declare EPSG:3946 with proj4 - * itowns.proj4.defs('EPSG:3946', '+proj=lcc +lat_1=45.25 +lat_2=46.75 +lat_0=46 +lon_0=3 +x_0=1700000 +y_0=5200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'); - * - * @property {boolean} isCoordinates - Used to checkout whether this coordinates - * is a Coordinates. Default is true. You should not change this, as it is used - * internally for optimisation. - * @property {string} crs - A supported crs by default in - * [`proj4js`](https://github.com/proj4js/proj4js#named-projections), or an - * added crs to `proj4js` (using `proj4.defs`). Note that `EPSG:4978` is also - * supported by default in itowns. - * @property {number} x - The first value of the coordinate. - * @property {number} y - The second value of the coordinate. - * @property {number} z - The third value of the coordinate. - * @property {number} latitude - The first value of the coordinate. - * @property {number} longitude - The second value of the coordinate. - * @property {number} altitude - The third value of the coordinate. - * @property {THREE.Vector3} geodesicNormal - The geodesic normal of the - * coordinate. - */ -class Coordinates { - /** - * @param {string} crs - A supported Coordinate Reference System. - * @param {number|Array|Coordinates|THREE.Vector3} [v0=0] - - * x or longitude value, or a more complex one: it can be an array of three - * numbers, being x/lon, y/lat, z/alt, or it can be `THREE.Vector3`. It can - * also simply be a Coordinates. - * @param {number} [v1=0] - y or latitude value. - * @param {number} [v2=0] - z or altitude value. - */ - constructor(crs, v0 = 0, v1 = 0, v2 = 0) { - this.isCoordinates = true; - - CRS.isValid(crs); - this.crs = crs; - - // Storing the coordinates as is, not in arrays, as it is - // slower (see https://jsbench.me/40jumfag6g/1) - this.x = 0; - this.y = 0; - this.z = 0; - - // Normal - this._normal = new THREE.Vector3(); - - if (v0.length > 0) { - this.setFromArray(v0); - } else if (v0.isVector3 || v0.isCoordinates) { - this.setFromVector3(v0); - } else { - this.setFromValues(v0, v1, v2); - } - - this._normalNeedsUpdate = true; - } - - /** - * Sets the Coordinate Reference System. - * @param {String} crs Coordinate Reference System (e.g. 'EPSG:4978') - */ - setCrs(crs) { - CRS.isValid(crs); - this.crs = crs; - } - - /** - * Set the values of this Coordinates. - * - * @param {number} [v0=0] - x or longitude value. - * @param {number} [v1=0] - y or latitude value. - * @param {number} [v2=0] - z or altitude value. - * - * @return {Coordinates} This Coordinates. - */ - setFromValues(v0 = 0, v1 = 0, v2 = 0) { - this.x = v0 == undefined ? 0 : v0; - this.y = v1 == undefined ? 0 : v1; - this.z = v2 == undefined ? 0 : v2; - - this._normalNeedsUpdate = true; - return this; - } - - /** - * Set the values of this Coordinates from an array. - * - * @param {Array} array - An array of number to assign to the - * Coordinates. - * @param {number} [offset] - Optional offset into the array. - * - * @return {Coordinates} This Coordinates. - */ - setFromArray(array, offset = 0) { - return this.setFromValues(array[offset], array[offset + 1], array[offset + 2]); - } - - /** - * Set the values of this Coordinates from a `THREE.Vector3` or an `Object` - * having `x/y/z` properties, like a `Coordinates`. - * - * @param {THREE.Vector3|Coordinates} v0 - The object to read the values - * from. - * - * @return {Coordinates} This Coordinates. - */ - setFromVector3(v0) { - return this.setFromValues(v0.x, v0.y, v0.z); - } - - /** - * Returns a new Coordinates with the same values as this one. It will - * instantiate a new Coordinates with the same CRS as this one. - * - * @return {Coordinates} The target with its new coordinates. - */ - clone() { - return new Coordinates(this.crs, this); - } - - /** - * Copies the values of the passed Coordinates to this one. The CRS is - * however not copied. - * - * @param {Coordinates} src - The source to copy from. - * - * @return {Coordinates} This Coordinates. - */ - copy(src) { - this.crs = src.crs; - return this.setFromVector3(src); - } - - get longitude() { - return this.x; - } - - get latitude() { - return this.y; - } - - get altitude() { - return this.z; - } - - set altitude(value) { - this.z = value; - } - - get geodesicNormal() { - if (this._normalNeedsUpdate) { - this._normalNeedsUpdate = false; - - if (CRS.is4326(this.crs)) { - ellipsoid.geodeticSurfaceNormalCartographic(this, this._normal); - } else if (this.crs == 'EPSG:4978') { - ellipsoid.geodeticSurfaceNormal(this, this._normal); - } else { - this._normal.set(0, 0, 1); - } - } - - return this._normal; - } - - /** - * Return this Coordinates values into a `THREE.Vector3`. - * - * @param {THREE.Vector3} [target] - The target to put the values in. If not - * specified, a new vector will be created. - * - * @return {THREE.Vector3} - */ - toVector3(target = new THREE.Vector3()) { - return target.copy(this); - } - - /** - * Copy values coordinates to array - * - * @param {number[]} array - array to store this vector to. If this is not - * provided a new array will be created. - * @param {number} [offset=0] - optional offset into the array. - * - * @return {number[]} Returns an array [x, y, z], or copies x, y and z into - * the provided array. - */ - toArray(array = [], offset = 0) { - return THREE.Vector3.prototype.toArray.call(this, array, offset); - } - - /** - * Calculate planar distance between this coordinates and `coord`. - * Planar distance is the straight-line euclidean distance calculated in a 2D cartesian coordinate system. - * - * @param {Coordinates} coord The coordinate - * @return {number} planar distance - * - */ - planarDistanceTo(coord) { - this.toVector3(v0).setZ(0); - coord.toVector3(v1).setZ(0); - return v0.distanceTo(v1); - } - - /** - * Calculate geodetic distance between this coordinates and `coord`. - * **Geodetic distance** is calculated in an ellispoid space as the shortest distance - * across the curved surface of the world. - * - * => As the crow flies/ Orthodromy - * - * @param {Coordinates} coord The coordinate - * @return {number} geodetic distance - * - */ - geodeticDistanceTo(coord) { - this.as('EPSG:4326', coord0); - coord.as('EPSG:4326', coord1); - return ellipsoid.geodesicDistance(coord0, coord1); - } - - /** - * Calculate earth euclidean distance between this coordinates and `coord`. - * - * @param {Coordinates} coord The coordinate - * @return {number} earth euclidean distance - * - */ - spatialEuclideanDistanceTo(coord) { - this.as('EPSG:4978', coord0).toVector3(v0); - coord.as('EPSG:4978', coord1).toVector3(v1); - return v0.distanceTo(v1); - } - - /** - * Multiplies this `coordinates` (with an implicit 1 in the 4th dimension) and `mat`. - * - * @param {THREE.Matrix4} mat The matrix. - * @return {Coordinates} return this object. - */ - applyMatrix4(mat) { - return THREE.Vector3.prototype.applyMatrix4.call(this, mat); - } - - /** - * Returns coordinates in the wanted [CRS](http://inspire.ec.europa.eu/theme/rs). - * - * @param {string} crs - The CRS to convert the Coordinates into. - * @param {Coordinates} [target] - The target to put the converted - * Coordinates into. If not specified a new one will be created. - * - * @return {Coordinates} - The resulting Coordinates after the conversion. - * - * @example - * const position = { longitude: 2.33, latitude: 48.24, altitude: 24999549 }; - * const coords = new Coordinates('EPSG:4326', position.longitude, position.latitude, position.altitude); // Geographic system - * const coordinates = coords.as('EPSG:4978'); // Geocentric system - * - * @example - * const position = { x: 20885167, y: 849862, z: 23385912 }; - * const coords = new Coordinates('EPSG:4978', position.x, position.y, position.z); // Geocentric system - * const coordinates = coords.as('EPSG:4326'); // Geographic system - * - * @example - * new Coordinates('EPSG:4326', longitude: 2.33, latitude: 48.24, altitude: 24999549).as('EPSG:4978'); // Geocentric system - * - * @example - * new Coordinates('EPSG:4978', x: 20885167, y: 849862, z: 23385912).as('EPSG:4326'); // Geographic system - */ - as(crs, target = new Coordinates(crs)) { - if (this.crs == crs) { - target.copy(this); - } else { - if (CRS.is4326(this.crs) && crs == 'EPSG:3857') { - this.y = THREE.MathUtils.clamp(this.y, -89.999999, 89.999999); - } - - target.setFromArray(proj4cache(this.crs, crs).forward([this.x, this.y, this.z])); - } - - target.crs = crs; - - return target; - } -} - -coord0 = new Coordinates('EPSG:4326', 0, 0, 0); -coord1 = new Coordinates('EPSG:4326', 0, 0, 0); - -export default Coordinates; diff --git a/src/Core/Geographic/Coordinates.ts b/src/Core/Geographic/Coordinates.ts new file mode 100644 index 0000000000..ca134ee415 --- /dev/null +++ b/src/Core/Geographic/Coordinates.ts @@ -0,0 +1,364 @@ +import * as THREE from 'three'; +import proj4 from 'proj4'; +import * as CRS from 'Core/Geographic/Crs'; +import Ellipsoid from 'Core/Math/Ellipsoid'; + +import type { ProjectionLike } from './Crs'; + +const ellipsoid = new Ellipsoid(); +const projectionCache: Record> = {}; + +const v0 = new THREE.Vector3(); +const v1 = new THREE.Vector3(); + +let coord0: Coordinates; +let coord1: Coordinates; + +export interface CoordinatesLike { + readonly crs: string; + readonly x: number; + readonly y: number; + readonly z: number; +} + +function proj4cache(crsIn: string, crsOut: string): proj4.Converter { + if (!projectionCache[crsIn]) { + projectionCache[crsIn] = {}; + } + + if (!projectionCache[crsIn][crsOut]) { + projectionCache[crsIn][crsOut] = proj4(crsIn, crsOut); + } + + return projectionCache[crsIn][crsOut]; +} + +/** + * A class representing a geographic or geocentric coordinate. + * + * A coordinate is defined by a [CRS](http://inspire.ec.europa.eu/theme/rs) + * (Coordinate Reference System) and a 3-dimensional vector `(x, y, z)`. + * For geocentric projections, it is recommended to use the `latitude`, + * `longitude` and `altitude` aliases to refer to vector components. + * + * To change a value, prefer the use of the `set*` methods. + * + * By default, the `EPSG:4978` and `EPSG:4326` projections are supported. To use + * a different projection, it must have been declared previously with `proj4`. + * A comprehensive list of projections and their corresponding proj4 string can + * be found at [epsg.io](https://epsg.io/). + * + * @example Geocentric coordinates + * ```js + * new Coordinates('EPSG:4978', 20885167, 849862, 23385912); + * ``` + * + * @example Geographic coordinates + * ```js + * new Coordinates('EPSG:4326', 2.33, 48.24, 24999549); + * ``` + * + * @example Defining the EPSG:2154 projection with proj4 + * ```js + * proj4.defs('EPSG:2154', `+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 + * +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m + * +no_defs +type=crs`); + * ``` + */ +class Coordinates { + /** + * Read-only flag to check if a given object is of type `Coordinates`. + */ + readonly isCoordinates: boolean; + /** + * A default or user-defined CRS (see {@link ProjectionLike}). + */ + crs: ProjectionLike; + + /** The x value (or longitude) of this coordinate. */ + x: number; + /** The y value (or latitude) of this coordinate. */ + y: number; + /** The z value (or altitude) of this coordinate. */ + z: number; + + private _normal: THREE.Vector3; + private _normalNeedsUpdate: boolean; + + /** + * @param crs - A default or user-defined CRS (see {@link ProjectionLike}). + * @param x - x or longitude value. + * @param y - y or latitude value. + * @param z - z or altitude value. + */ + constructor(crs: ProjectionLike, x: number = 0, y: number = 0, z: number = 0) { + this.isCoordinates = true; + + CRS.isValid(crs); + this.crs = crs; + + // Storing the coordinates as is, not in arrays, as it is + // slower (see https://jsbench.me/40jumfag6g/1) + this.x = 0; + this.y = 0; + this.z = 0; + + // Normal + this._normal = new THREE.Vector3(); + + // eslint-disable-next-line @typescript-eslint/no-explicit-any + if ((x as any).length > 0) { // deepscan-disable-line + console.warn( + 'Deprecated Coordinates#constructor(string, number[]),', + 'use `new Coordinates(string).setFromArray(number[])` instead.', + ); + // eslint-disable-next-line @typescript-eslint/no-explicit-any + this.setFromArray(x as any); + // eslint-disable-next-line @typescript-eslint/no-explicit-any + } else if ((x as any).isVector3 || (x as any).isCoordinates) { + console.warn( + 'Deprecated Coordinates#constructor(string, Vector3),', + 'use `new Coordinates(string).setFromVector3(Vector3)` instead.', + ); + // eslint-disable-next-line @typescript-eslint/no-explicit-any + this.setFromVector3(x as any); + } else { + this.setFromValues(x, y, z); + } + + this._normalNeedsUpdate = true; + } + + /** + * Sets the Coordinate Reference System. + * @param crs - Coordinate Reference System (e.g. 'EPSG:4978') + */ + setCrs(crs: ProjectionLike): this { + CRS.isValid(crs); + this.crs = crs; + return this; + } + + /** + * Sets the x, y and z components of this coordinate. + * + * @param x - x or longitude value. + * @param y - y or latitude value. + * @param z - z or altitude value. + */ + setFromValues(x: number = 0, y: number = 0, z: number = 0): this { + this.x = x; + this.y = y; + this.z = z; + + this._normalNeedsUpdate = true; + return this; + } + + /** + * Sets the coordinates's {@link Coordinates#x | x} component to + * `array[offset + 0]`, {@link Coordinates#y | y} component to + * `array[offset + 1]` and {@link Coordinates#z | z} component to + * `array[offset + 2]`. + * + * @param array - The source array. + * @param offset - Optional offset into the array. Default is 0. + */ + setFromArray(array: number[], offset: number = 0): this { + return this.setFromValues( + array[offset], + array[offset + 1], + array[offset + 2], + ); + } + + /** + * Sets the `(x, y, z)` vector of this coordinate from a 3-dimensional + * vector-like object. This object shall have both `x`, `y` and `z` + * properties. + * + * @param v - The source object. + */ + setFromVector3(v: THREE.Vector3Like): this { + return this.setFromValues(v.x, v.y, v.z); + } + + /** + * Returns a new coordinate with the same `(x, y, z)` vector and crs as this + * one. + */ + clone(): Coordinates { + return new Coordinates(this.crs, this.x, this.y, this.z); + } + + /** + * Copies the `(x, y, z)` vector components and crs of the passed coordinate + * to this coordinate. + * + * @param src - The source coordinate to copy from. + */ + copy(src: CoordinatesLike): this { + this.crs = src.crs; + return this.setFromVector3(src); + } + + get longitude() { + return this.x; + } + + get latitude() { + return this.y; + } + + get altitude() { + return this.z; + } + + set altitude(value) { + this.z = value; + } + + /** + * The geodesic normal of the coordinate. + */ + get geodesicNormal() { + if (this._normalNeedsUpdate) { + this._normalNeedsUpdate = false; + + if (CRS.is4326(this.crs)) { + ellipsoid.geodeticSurfaceNormalCartographic(this, this._normal); + } else if (this.crs == 'EPSG:4978') { + ellipsoid.geodeticSurfaceNormal(this, this._normal); + } else { + this._normal.set(0, 0, 1); + } + } + + return this._normal; + } + + /** + * Copies the `x`, `y` and `z` components into the provided `THREE.Vector3`. + * + * @param target - An object to store this vector to. If this is not + * specified, a new vector will be created. + * + * @returns A vector `(x, y, z)`, or copies x, y and z into the provided + * vector. + */ + toVector3(target: THREE.Vector3 = new THREE.Vector3()): THREE.Vector3 { + return target.copy(this); + } + + /** + * Copies the `x`, `y` and `z` components into the provided array. + * + * @param array - An array to store this vector to. If this is not + * provided a new array will be created. + * @param offset - An optional offset into the array. + * + * @returns An array [x, y, z], or copies x, y and z into the provided + * array. + */ + toArray(array: number[] = [], offset: number = 0): ArrayLike { + return THREE.Vector3.prototype.toArray.call(this, array, offset); + } + + /** + * Computes the planar distance from this coordinates to `coord`. + * **Planar distance** is the straight-line euclidean distance calculated in + * a 2D cartesian coordinate system. + */ + planarDistanceTo(coord: Coordinates): number { + this.toVector3(v0).setZ(0); + coord.toVector3(v1).setZ(0); + return v0.distanceTo(v1); + } + + /** + * Computes the geodetic distance from this coordinates to `coord`. + * **Geodetic distance** is calculated in an ellipsoid space as the shortest + * distance across the curved surface of the ellipsoid. + */ + geodeticDistanceTo(coord: Coordinates): number { + this.as('EPSG:4326', coord0); + coord.as('EPSG:4326', coord1); + return ellipsoid.geodesicDistance(coord0, coord1); + } + + /** + * Computes the euclidean distance from this coordinates to `coord` in a + * WGS84 projection. + * + * @param coord - The coordinate + * @returns earth euclidean distance + */ + spatialEuclideanDistanceTo(coord: Coordinates): number { + this.as('EPSG:4978', coord0).toVector3(v0); + coord.as('EPSG:4978', coord1).toVector3(v1); + return v0.distanceTo(v1); + } + + /** + * Multiplies this coordinate (with an implicit 1 in the 4th dimension) + * by `mat`, and divides by perspective. + * + * @param mat - The matrix. + */ + applyMatrix4(mat: THREE.Matrix4): this { + THREE.Vector3.prototype.applyMatrix4.call(this, mat); + return this; + } + + /** + * Projects this coordinate to the specified + * [CRS](http://inspire.ec.europa.eu/theme/rs). + * + * @param crs - The target CRS to which the coordinate will be converted. + * @param target - The target to store the projected coordinate. If this not + * provided a new coordinate will be created. + * + * @returns The coordinate projected into the specified CRS. + * + * @example Conversion from a geographic to a geocentric reference system + * ```js + * const geographicCoords = new Coordinates('EPSG:4326', + * 2.33, // longitude + * 48.24, // latitude + * 24999549, // altitude + * ); + * const geocentricCoords = geographicCoords.as('EPSG:4978'); + * ``` + * + * @example Conversion from a geocentric to a geographic reference system + * ```js + * const geocentricCoords = new Coordinates('EPSG:4978', + * 20885167, // x + * 849862, // y + * 23385912, // z + * ); + * const geographicCoords = geocentricCoords.as('EPSG:4326'); + * ``` + */ + as(crs: ProjectionLike, target = new Coordinates(crs)): Coordinates { + if (this.crs == crs) { + target.copy(this); + } else { + if (CRS.is4326(this.crs) && crs == 'EPSG:3857') { + this.y = THREE.MathUtils.clamp(this.y, -89.999999, 89.999999); + } + + target.setFromArray(proj4cache(this.crs, crs) + .forward([this.x, this.y, this.z])); + } + + target.crs = crs; + + return target; + } +} + +coord0 = new Coordinates('EPSG:4326', 0, 0, 0); +coord1 = new Coordinates('EPSG:4326', 0, 0, 0); + +export default Coordinates; diff --git a/src/Core/Math/Ellipsoid.js b/src/Core/Math/Ellipsoid.js deleted file mode 100644 index 8c0ea17221..0000000000 --- a/src/Core/Math/Ellipsoid.js +++ /dev/null @@ -1,175 +0,0 @@ -import * as THREE from 'three'; -import proj4 from 'proj4'; -import Coordinates from 'Core/Geographic/Coordinates'; - -export const ellipsoidSizes = new THREE.Vector3( - proj4.WGS84.a, - proj4.WGS84.a, - proj4.WGS84.b); - -const normal = new THREE.Vector3(); - -class Ellipsoid { - constructor(size = ellipsoidSizes) { - this.size = new THREE.Vector3(); - this._radiiSquared = new THREE.Vector3(); - this._invRadiiSquared = new THREE.Vector3(); - this.eccentricity = 0; - - this.setSize(size); - } - - geodeticSurfaceNormal(cartesian, target = new THREE.Vector3()) { - return cartesian.toVector3(target).multiply(this._invRadiiSquared).normalize(); - } - - geodeticSurfaceNormalCartographic(coordCarto, target = new THREE.Vector3()) { - const longitude = THREE.MathUtils.degToRad(coordCarto.longitude); - const latitude = THREE.MathUtils.degToRad(coordCarto.latitude); - const cosLatitude = Math.cos(latitude); - - return target.set(cosLatitude * Math.cos(longitude), - cosLatitude * Math.sin(longitude), - Math.sin(latitude)); - } - - setSize(size) { - this.size.set(size.x, size.y, size.z); - - this._radiiSquared.multiplyVectors(size, size); - - this._invRadiiSquared.x = (size.x == 0) ? 0 : 1 / this._radiiSquared.x; - this._invRadiiSquared.y = (size.y == 0) ? 0 : 1 / this._radiiSquared.y; - this._invRadiiSquared.z = (size.z == 0) ? 0 : 1 / this._radiiSquared.z; - - this.eccentricity = Math.sqrt(this._radiiSquared.x - this._radiiSquared.z) / this.size.x; - } - - cartographicToCartesian(coordCarto, target = new THREE.Vector3()) { - normal.copy(coordCarto.geodesicNormal); - - target.multiplyVectors(this._radiiSquared, normal); - - const gamma = Math.sqrt(normal.dot(target)); - - target.divideScalar(gamma); - - normal.multiplyScalar(coordCarto.altitude); - - return target.add(normal); - } - - /** - * Convert cartesian coordinates to geographic according to the current ellipsoid of revolution. - * @param {Object} position - The coordinate to convert - * @param {number} position.x - * @param {number} position.y - * @param {number} position.z - * @param {Coordinate} [target] coordinate to copy result - * @returns {Coordinate} an object describing the coordinates on the reference ellipsoid, angles are in degree - */ - cartesianToCartographic(position, target = new Coordinates('EPSG:4326', 0, 0, 0)) { - // for details, see for example http://www.linz.govt.nz/data/geodetic-system/coordinate-conversion/geodetic-datum-conversions/equations-used-datum - // TODO the following is only valable for oblate ellipsoid of revolution. do we want to support triaxial ellipsoid? - const R = Math.sqrt(position.x * position.x + position.y * position.y + position.z * position.z); - const a = this.size.x; // x - const b = this.size.z; // z - const e = Math.abs((a * a - b * b) / (a * a)); - const f = 1 - Math.sqrt(1 - e); - const rsqXY = Math.sqrt(position.x * position.x + position.y * position.y); - - const theta = Math.atan2(position.y, position.x); - const nu = Math.atan(position.z / rsqXY * ((1 - f) + e * a / R)); - - const sinu = Math.sin(nu); - const cosu = Math.cos(nu); - - const phi = Math.atan((position.z * (1 - f) + e * a * sinu * sinu * sinu) / ((1 - f) * (rsqXY - e * a * cosu * cosu * cosu))); - - const h = (rsqXY * Math.cos(phi)) + position.z * Math.sin(phi) - a * Math.sqrt(1 - e * Math.sin(phi) * Math.sin(phi)); - - return target.setFromValues(THREE.MathUtils.radToDeg(theta), THREE.MathUtils.radToDeg(phi), h); - } - - cartographicToCartesianArray(coordCartoArray) { - const cartesianArray = []; - for (let i = 0; i < coordCartoArray.length; i++) { - cartesianArray.push(this.cartographicToCartesian(coordCartoArray[i])); - } - - return cartesianArray; - } - - intersection(ray) { - const EPSILON = 0.0001; - const O_C = ray.origin; - const dir = ray.direction; - // normalizeVector( dir ); - - const a = - ((dir.x * dir.x) * this._invRadiiSquared.x) + ((dir.y * dir.y) * this._invRadiiSquared.y) + ((dir.z * dir.z) * this._invRadiiSquared.z); - - const b = - ((2 * O_C.x * dir.x) * this._invRadiiSquared.x) + ((2 * O_C.y * dir.y) * this._invRadiiSquared.y) + ((2 * O_C.z * dir.z) * this._invRadiiSquared.z); - const c = - ((O_C.x * O_C.x) * this._invRadiiSquared.x) + ((O_C.y * O_C.y) * this._invRadiiSquared.y) + ((O_C.z * O_C.z) * this._invRadiiSquared.z) - 1; - - let d = ((b * b) - (4 * a * c)); - if (d < 0 || a === 0 || b === 0 || c === 0) { return false; } - - d = Math.sqrt(d); - - const t1 = (-b + d) / (2 * a); - const t2 = (-b - d) / (2 * a); - - if (t1 <= EPSILON && t2 <= EPSILON) { return false; } // both intersections are behind the ray origin - // var back = (t1 <= EPSILON || t2 <= EPSILON); // If only one intersection (t>0) then we are inside the ellipsoid and the intersection is at the back of the ellipsoid - let t = 0; - if (t1 <= EPSILON) { t = t2; } else - if (t2 <= EPSILON) { t = t1; } else { t = (t1 < t2) ? t1 : t2; } - - if (t < EPSILON) { return false; } // Too close to intersection - - const inter = new THREE.Vector3(); - - inter.addVectors(ray.origin, dir.clone().setLength(t)); - - return inter; - } - - computeDistance(coordCarto1, coordCarto2) { - console.warn('computeDistance is renamed to geodesicDistance'); - this.geodesicDistance(coordCarto1, coordCarto2); - } - - /** - * Calculate the geodesic distance, between coordCarto1 and coordCarto2. - * It's most short distance on ellipsoid surface between coordCarto1 and coordCarto2. - * It's called orthodromy. - * - * @param {Coordinates} coordCarto1 The coordinate carto 1 - * @param {Coordinates} coordCarto2 The coordinate carto 2 - * @return {number} The orthodromic distance between the two given coordinates. - */ - geodesicDistance(coordCarto1, coordCarto2) { - // The formula uses the distance on approximated sphere, - // with the nearest local radius of curvature of the ellipsoid - // https://geodesie.ign.fr/contenu/fichiers/Distance_longitude_latitude.pdf - const longitude1 = THREE.MathUtils.degToRad(coordCarto1.longitude); - const latitude1 = THREE.MathUtils.degToRad(coordCarto1.latitude); - const longitude2 = THREE.MathUtils.degToRad(coordCarto2.longitude); - const latitude2 = THREE.MathUtils.degToRad(coordCarto2.latitude); - - const distRad = Math.acos(Math.sin(latitude1) * Math.sin(latitude2) + Math.cos(latitude1) * Math.cos(latitude2) * Math.cos(longitude2 - longitude1)); - - const e = this.eccentricity; - const latMoy = (latitude1 + latitude2) * 0.5; - const es = (e * Math.sin(latMoy)) ** 2; - const rho = this.size.x * (1 - e ** 2) / ((1 - es) ** (3 / 2)); - const N = this.size.x / Math.sqrt(1 - es); - - return distRad * Math.sqrt(rho * N); - } -} - -export default Ellipsoid; diff --git a/src/Core/Math/Ellipsoid.ts b/src/Core/Math/Ellipsoid.ts new file mode 100644 index 0000000000..67ef10d405 --- /dev/null +++ b/src/Core/Math/Ellipsoid.ts @@ -0,0 +1,250 @@ +import * as THREE from 'three'; +import proj4 from 'proj4'; +import Coordinates from '../Geographic/Coordinates'; + +/** + * Length of the semi-axes of the WGS84 ellipsoid. + * @internal + */ +export const ellipsoidSizes = new THREE.Vector3( + proj4.WGS84.a, + proj4.WGS84.a, + proj4.WGS84.b); + +const normal = new THREE.Vector3(); + +class Ellipsoid { + /** + * Length of the semi-axes of the ellipsoid. + */ + size: THREE.Vector3; + /** + * Eccentricity of the ellipsoid. + */ + eccentricity: number; + + private _radiiSquared: THREE.Vector3; + private _invRadiiSquared: THREE.Vector3; + + /** + * @param size - Length of the semi-axes of the ellipsoid. Defaults to those + * defined by the WGS84 ellipsoid. + */ + constructor(size: THREE.Vector3 = ellipsoidSizes) { + this.size = new THREE.Vector3(); + this._radiiSquared = new THREE.Vector3(); + this._invRadiiSquared = new THREE.Vector3(); + this.eccentricity = 0; + + this.setSize(size); + } + + /** + * Computes the normal vector to an ellipsoid at the given cartesian + * coordinate `(x, y, z)`. + * + * @param cartesian - The given cartesian coordinate. + * @param target - An object to store this vector to. If this is not + * specified, a new vector will be created. + */ + geodeticSurfaceNormal( + cartesian: Coordinates, + target = new THREE.Vector3(), + ): THREE.Vector3 { + return cartesian.toVector3(target).multiply(this._invRadiiSquared).normalize(); + } + + /** + * Computes the normal vector to an ellipsoid at the given geographic + * coordinate `(longitude, latitude, altitude)`. + * + * @param coordCarto - The given geographic coordinate. + * @param target - An object to store this vector to. If this is not + * specified, a new vector will be created. + */ + geodeticSurfaceNormalCartographic( + coordCarto: Coordinates, + target = new THREE.Vector3(), + ): THREE.Vector3 { + const longitude = THREE.MathUtils.degToRad(coordCarto.longitude); + const latitude = THREE.MathUtils.degToRad(coordCarto.latitude); + const cosLatitude = Math.cos(latitude); + + return target.set(cosLatitude * Math.cos(longitude), + cosLatitude * Math.sin(longitude), + Math.sin(latitude)); + } + + /** + * Sets the length of the semi-axes of this ellipsoid from a 3-dimensional + * vector-like object. The object shall have both `x`, `y` and `z` + * properties. + * + * @param size - The source vector. + */ + setSize(size: THREE.Vector3Like): this { + this.size.set(size.x, size.y, size.z); + + this._radiiSquared.multiplyVectors(size, size); + + this._invRadiiSquared.x = (size.x === 0) ? 0 : 1 / this._radiiSquared.x; + this._invRadiiSquared.y = (size.y === 0) ? 0 : 1 / this._radiiSquared.y; + this._invRadiiSquared.z = (size.z === 0) ? 0 : 1 / this._radiiSquared.z; + + this.eccentricity = Math.sqrt(this._radiiSquared.x - this._radiiSquared.z) / this.size.x; + + return this; + } + + cartographicToCartesian( + coordCarto: Coordinates, + target = new THREE.Vector3(), + ): THREE.Vector3 { + normal.copy(coordCarto.geodesicNormal); + + target.multiplyVectors(this._radiiSquared, normal); + + const gamma = Math.sqrt(normal.dot(target)); + + target.divideScalar(gamma); + + normal.multiplyScalar(coordCarto.altitude); + + return target.add(normal); + } + + /** + * Convert cartesian coordinates to geographic according to the current + * ellipsoid of revolution. + * @param position - The coordinate to convert + * @param target - coordinate to copy result + * @returns an object describing the coordinates on the reference ellipsoid, + * angles are in degree + */ + cartesianToCartographic( + position: THREE.Vector3Like, + target = new Coordinates('EPSG:4326', 0, 0, 0), + ): Coordinates { + // for details, see for example http://www.linz.govt.nz/data/geodetic-system/coordinate-conversion/geodetic-datum-conversions/equations-used-datum + // TODO the following is only valable for oblate ellipsoid of + // revolution. do we want to support triaxial ellipsoid? + const R = Math.sqrt(position.x * position.x + + position.y * position.y + + position.z * position.z); + const a = this.size.x; // x + const b = this.size.z; // z + const e = Math.abs((a * a - b * b) / (a * a)); + const f = 1 - Math.sqrt(1 - e); + const rsqXY = Math.sqrt(position.x * position.x + position.y * position.y); + + const theta = Math.atan2(position.y, position.x); + const nu = Math.atan(position.z / rsqXY * ((1 - f) + e * a / R)); + + const sinu = Math.sin(nu); + const cosu = Math.cos(nu); + + const phi = Math.atan( + (position.z * (1 - f) + e * a * sinu * sinu * sinu) / + ((1 - f) * (rsqXY - e * a * cosu * cosu * cosu))); + + const h = (rsqXY * Math.cos(phi)) + + position.z * Math.sin(phi) - + a * Math.sqrt(1 - e * Math.sin(phi) * Math.sin(phi)); + + return target.setFromValues( + THREE.MathUtils.radToDeg(theta), + THREE.MathUtils.radToDeg(phi), + h, + ); + } + + cartographicToCartesianArray(coordCartoArray: Coordinates[]): THREE.Vector3[] { + const cartesianArray = []; + for (let i = 0; i < coordCartoArray.length; i++) { + cartesianArray.push(this.cartographicToCartesian(coordCartoArray[i])); + } + + return cartesianArray; + } + + intersection(ray: THREE.Ray): THREE.Vector3 | false { + const EPSILON = 0.0001; + const O_C = ray.origin; + const dir = ray.direction; + // normalizeVector( dir ); + + const a = + ((dir.x * dir.x) * this._invRadiiSquared.x) + + ((dir.y * dir.y) * this._invRadiiSquared.y) + + ((dir.z * dir.z) * this._invRadiiSquared.z); + + const b = + ((2 * O_C.x * dir.x) * this._invRadiiSquared.x) + + ((2 * O_C.y * dir.y) * this._invRadiiSquared.y) + + ((2 * O_C.z * dir.z) * this._invRadiiSquared.z); + + const c = + ((O_C.x * O_C.x) * this._invRadiiSquared.x) + + ((O_C.y * O_C.y) * this._invRadiiSquared.y) + + ((O_C.z * O_C.z) * this._invRadiiSquared.z) - 1; + + let d = ((b * b) - (4 * a * c)); + if (d < 0 || a === 0 || b === 0 || c === 0) { return false; } + + d = Math.sqrt(d); + + const t1 = (-b + d) / (2 * a); + const t2 = (-b - d) / (2 * a); + + if (t1 <= EPSILON && t2 <= EPSILON) { + // both intersections are behind the ray origin + return false; + } + + let t = 0; + if (t1 <= EPSILON) { t = t2; } else + if (t2 <= EPSILON) { t = t1; } else { t = (t1 < t2) ? t1 : t2; } + + if (t < EPSILON) { return false; } // Too close to intersection + + const inter = new THREE.Vector3(); + + inter.addVectors(ray.origin, dir.clone().setLength(t)); + + return inter; + } + + /** + * Calculate the geodesic distance, between coordCarto1 and coordCarto2. + * It's most short distance on ellipsoid surface between coordCarto1 and + * coordCarto2. + * It's called orthodromy. + * + * @param coordCarto1 - The coordinate carto 1 + * @param coordCarto2 - The coordinate carto 2 + * @returns The orthodromic distance between the two given coordinates. + */ + geodesicDistance(coordCarto1: Coordinates, coordCarto2: Coordinates): number { + // The formula uses the distance on approximated sphere, + // with the nearest local radius of curvature of the ellipsoid + // https://geodesie.ign.fr/contenu/fichiers/Distance_longitude_latitude.pdf + const longitude1 = THREE.MathUtils.degToRad(coordCarto1.longitude); + const latitude1 = THREE.MathUtils.degToRad(coordCarto1.latitude); + const longitude2 = THREE.MathUtils.degToRad(coordCarto2.longitude); + const latitude2 = THREE.MathUtils.degToRad(coordCarto2.latitude); + + const distRad = Math.acos( + Math.sin(latitude1) * Math.sin(latitude2) + + Math.cos(latitude1) * Math.cos(latitude2) * Math.cos(longitude2 - longitude1)); + + const e = this.eccentricity; + const latMoy = (latitude1 + latitude2) * 0.5; + const es = (e * Math.sin(latMoy)) ** 2; + const rho = this.size.x * (1 - e ** 2) / ((1 - es) ** (3 / 2)); + const N = this.size.x / Math.sqrt(1 - es); + + return distRad * Math.sqrt(rho * N); + } +} + +export default Ellipsoid; diff --git a/src/Renderer/Camera.js b/src/Renderer/Camera.js index 848936cf2b..99a29b9ee4 100644 --- a/src/Renderer/Camera.js +++ b/src/Renderer/Camera.js @@ -188,7 +188,9 @@ class Camera { * @return {Coordinates} Coordinates object holding camera's position. */ position(crs) { - return new Coordinates(this.crs, this.camera3D.position).as(crs || this.crs); + return new Coordinates(this.crs) + .setFromVector3(this.camera3D.position) + .as(crs || this.crs); } /** diff --git a/src/Utils/CameraUtils.js b/src/Utils/CameraUtils.js index 5bccc8ece6..a586a242e0 100644 --- a/src/Utils/CameraUtils.js +++ b/src/Utils/CameraUtils.js @@ -154,7 +154,10 @@ class CameraRig extends THREE.Object3D { // set rig's objects transformation from camera's position and target's position setFromPositions(view, cameraPosition) { - this.setTargetFromCoordinate(view, new Coordinates(view.referenceCrs, targetPosition)); + this.setTargetFromCoordinate( + view, + new Coordinates(view.referenceCrs).setFromVector3(targetPosition), + ); this.target.rotation.set(0, 0, 0); this.updateMatrixWorld(true); this.camera.position.copy(cameraPosition); diff --git a/test/unit/ellipsoid.js b/test/unit/ellipsoid.js index 877dbe3be3..88f5c23fe9 100644 --- a/test/unit/ellipsoid.js +++ b/test/unit/ellipsoid.js @@ -5,7 +5,20 @@ import Ellipsoid from 'Core/Math/Ellipsoid'; describe('Ellipsoid', function () { const c1 = new Coordinates('EPSG:4326', 0, 0, 0); const ellipsoid = new Ellipsoid(); - it('geodeticSurfaceNormalCartographic', () => { + + it('geodeticSurfaceNormal', function () { + c1.setFromValues(6378137, 0, 0); + const v = ellipsoid.geodeticSurfaceNormal(c1); + assert.equal(v.x, 1); + assert.equal(v.y, 0); + assert.equal(v.z, 0); + c1.x = -6378137; + ellipsoid.geodeticSurfaceNormal(c1, v); + assert.equal(v.x, -1); + }); + + it('geodeticSurfaceNormalCartographic', function () { + c1.setFromValues(0, 0, 0); const v = ellipsoid.geodeticSurfaceNormalCartographic(c1); assert.equal(v.x, 1); assert.equal(v.y, 0); @@ -15,13 +28,14 @@ describe('Ellipsoid', function () { assert.equal(v.x, -1); }); - it('cartographicToCartesian', () => { - c1.x = 0; + it('cartographicToCartesian', function () { + c1.setFromValues(0, 0, 0); const v = ellipsoid.cartographicToCartesian(c1); assert.equal(v.x, ellipsoid.size.x); }); - it('cartesianToCartographic', () => { + it('cartesianToCartographic', function () { + c1.setFromValues(0, 0, 0); const altitude = 2000; const v = ellipsoid.cartographicToCartesian(c1); v.x += altitude; @@ -29,14 +43,15 @@ describe('Ellipsoid', function () { assert.equal(c1.z, altitude); }); - it('cartographicToCartesianArray', () => { - c1.z = 0; + it('cartographicToCartesianArray', function () { + c1.setFromValues(0, 0, 0); const a = ellipsoid.cartographicToCartesianArray([c1]); assert.equal(a.length, 1); assert.equal(a[0].x, ellipsoid.size.x); }); - it('geodesic distance', () => { + it('geodesic distance', function () { + c1.setFromValues(0, 0, 0); const a = ellipsoid.geodesicDistance(c1, c1); assert.equal(a, 0); const c2 = new Coordinates('EPSG:4326', 180, 0, 0); diff --git a/utils/debug/Debug.js b/utils/debug/Debug.js index 97420a456b..55361cf5b6 100644 --- a/utils/debug/Debug.js +++ b/utils/debug/Debug.js @@ -177,7 +177,9 @@ function Debug(view, datDebugTool, chartDivContainer) { const size = { x: g.width * ratio, y: g.height * ratio }; debugCamera.aspect = size.x / size.y; const camera = view.camera3D; - const coord = new Coordinates(view.referenceCrs, camera.position).as(tileLayer.extent.crs); + const coord = new Coordinates(view.referenceCrs) + .setFromVector3(camera.position) + .as(tileLayer.extent.crs); const extent = view.tileLayer.info.displayed.extent; displayedTilesObb.setFromExtent(extent); displayedTilesObbHelper.visible = true;