From 59c7573e2fad263847dd8993868839076d0cb883 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schn=C3=BCrer?= Date: Thu, 26 Dec 2024 19:18:49 +0000 Subject: [PATCH 01/13] ported methods: - mod - zeroToTwoPi - negativePiToPi --- modules/core/src/index.ts | 3 ++ modules/core/src/lib/common.ts | 50 ++++++++++++++++++++++++++++++++++ 2 files changed, 53 insertions(+) diff --git a/modules/core/src/index.ts b/modules/core/src/index.ts index 373c167c..16ae8c34 100644 --- a/modules/core/src/index.ts +++ b/modules/core/src/index.ts @@ -58,6 +58,9 @@ export { // math.gl "GLSL"-style functions radians, degrees, + mod, + zeroToTwoPi, + negativePiToPi, sin, cos, tan, diff --git a/modules/core/src/lib/common.ts b/modules/core/src/lib/common.ts index c571f93f..f3fc79d2 100644 --- a/modules/core/src/lib/common.ts +++ b/modules/core/src/lib/common.ts @@ -7,6 +7,7 @@ import type {NumericArray} from '@math.gl/types'; import type {MathArray} from '../classes/base/math-array'; +import {EPSILON14, TWO_PI, PI} from './math-utils'; const RADIANS_TO_DEGREES = (1 / Math.PI) * 180; const DEGREES_TO_RADIANS = (1 / 180) * Math.PI; @@ -122,6 +123,55 @@ export function degrees( return map(radians, (radians) => radians * RADIANS_TO_DEGREES, result); } + +/** + * The modulo operation that also works for negative dividends. + * + * @param {number} m The dividend. + * @param {number} n The divisor. + * @returns {number} The remainder. + */ +export function mod(m: number, n: number): number { + if (Math.sign(m) === Math.sign(n) && Math.abs(m) < Math.abs(n)) { + return m; + } + + return ((m % n) + n) % n; +} + +/** + * Produces an angle in the range 0 <= angle <= 2Pi which is equivalent to the provided angle. + * + * @param {number} angle in radians + * @returns {number} The angle in the range [0, TWO_PI]. + */ +export function zeroToTwoPi(angle: number): number { + if (angle >= 0 && angle <= TWO_PI) { + return angle; + } + const remainder = mod(angle, TWO_PI); + if ( + Math.abs(remainder) < EPSILON14 && + Math.abs(angle) > EPSILON14 + ) { + return TWO_PI; + } + return remainder; +} + +/** + * Produces an angle in the range -Pi <= angle <= Pi which is equivalent to the provided angle. + * + * @param {number} angle in radians + * @returns {number} The angle in the range [-PI, PI]. + */ +export function negativePiToPi(angle: number): number { + if (angle >= -PI && angle <= PI) { + return angle; + } + return zeroToTwoPi(angle + PI) - PI; +} + /** * "GLSL equivalent" of `Math.sin`: Works on single values and vectors * @deprecated From 647c759b352dcb517d8062afbd170998772dc0ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schn=C3=BCrer?= Date: Thu, 26 Dec 2024 19:19:22 +0000 Subject: [PATCH 02/13] partly ported class 'EllipsoidTangentPlane' from Cesium --- .../src/lib/ellipsoid-tangent-plane.ts | 102 ++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100644 modules/culling/src/lib/ellipsoid-tangent-plane.ts diff --git a/modules/culling/src/lib/ellipsoid-tangent-plane.ts b/modules/culling/src/lib/ellipsoid-tangent-plane.ts new file mode 100644 index 00000000..4230abc1 --- /dev/null +++ b/modules/culling/src/lib/ellipsoid-tangent-plane.ts @@ -0,0 +1,102 @@ +// math.gl +// SPDX-License-Identifier: MIT and Apache-2.0 +// Copyright (c) vis.gl contributors + +// This file is derived from the Cesium library under Apache 2 license +// See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md + +import {Vector2, Vector3, Matrix4} from '@math.gl/core'; +// @ts-ignore +// eslint-disable-next-line import/no-unresolved +import {Ellipsoid} from '@math.gl/geospatial'; +import {Plane} from './plane'; +import {Ray} from './ray'; + +const scratchOrigin = new Vector3(); +const scratchCart3 = new Vector3(); +const scratchEastNorthUp = new Matrix4(); +const scratchPlane = new Plane(); + +const scratchProjectPointOntoPlaneRay = new Ray(); +const scratchProjectPointOntoPlaneCartesian3 = new Vector3(); +const scratchDirection = new Vector3(); + +/** A plane tangent to the WGS84 ellipsoid at the provided origin */ +export class EllipsoidTangentPlane { + private _origin: Vector3; + private _xAxis: Vector3; + private _yAxis: Vector3; + private _plane: Plane; + + /** + * Creates a new plane tangent to the WGS84 ellipsoid at the provided origin. + * If origin is not on the surface of the ellipsoid, it's surface projection will be used. + * + * @param {Cartesian3} origin The point on the surface of the ellipsoid where the tangent plane touches. + */ + constructor(origin: number[]) { + origin = Ellipsoid.WGS84.scaleToGeodeticSurface(origin, scratchOrigin); + + const eastNorthUp = Ellipsoid.WGS84.eastNorthUpToFixedFrame(origin, scratchEastNorthUp); + + this._origin = origin as Vector3; + this._xAxis = new Vector3(scratchCart3.from(eastNorthUp.getColumn(0))); + this._yAxis = new Vector3(scratchCart3.from(eastNorthUp.getColumn(1))); + const normal = new Vector3(scratchCart3.from(eastNorthUp.getColumn(2))); + + this._plane = scratchPlane.fromPointNormal(origin, normal); + } + + /** + * Computes the projection of the provided 3D position onto the 2D plane, along the plane normal. + * + * @param {Vector3} cartesian The point to project. + * @param {Vector2} [result] The object onto which to store the result. + * @returns {Vector2} The modified result parameter or a new Cartesian2 instance if none was provided. + */ + projectPointToNearestOnPlane (cartesian: Vector3, result?: Vector2): Vector2 { + if (!result) + result = new Vector2(); + + const plane = this._plane; + + const ray = scratchProjectPointOntoPlaneRay; + scratchProjectPointOntoPlaneRay.origin = cartesian; + scratchProjectPointOntoPlaneRay.direction = scratchDirection.copy(plane.normal); + + let intersectionPoint = plane.intersectWithRay(ray, scratchProjectPointOntoPlaneCartesian3); + + if (!intersectionPoint) { + ray.direction = ray.direction.negate(); + intersectionPoint = plane.intersectWithRay(ray, scratchProjectPointOntoPlaneCartesian3); + } + + const v = intersectionPoint.subtract(this._origin); + const x = this._xAxis.dot(v); + const y = this._yAxis.dot(v); + + result.x = x; + result.y = y; + return result; + } + + get plane() { + return this._plane; + } + + get origin() { + return this._origin; + } + + get xAxis() { + return this._xAxis; + } + + get yAxis() { + return this._yAxis; + } + + get zAxis() { + return this._plane.normal; + } +} \ No newline at end of file From 16ef1f8443583c68b46354a232325e5f66da0f78 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schn=C3=BCrer?= Date: Thu, 26 Dec 2024 19:20:25 +0000 Subject: [PATCH 03/13] ported constant from Cesium: - PI --- modules/core/src/lib/math-utils.ts | 1 + 1 file changed, 1 insertion(+) diff --git a/modules/core/src/lib/math-utils.ts b/modules/core/src/lib/math-utils.ts index f28ea7c5..a523b10d 100644 --- a/modules/core/src/lib/math-utils.ts +++ b/modules/core/src/lib/math-utils.ts @@ -29,4 +29,5 @@ export const PI_OVER_TWO = Math.PI / 2; export const PI_OVER_FOUR = Math.PI / 4; export const PI_OVER_SIX = Math.PI / 6; +export const PI = Math.PI; export const TWO_PI = Math.PI * 2; From 30da0750b7deadeecd679d232862ebe32b78a921 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schn=C3=BCrer?= Date: Thu, 26 Dec 2024 19:21:01 +0000 Subject: [PATCH 04/13] ported methods from Cesium: - multiplyByVector - multiplyByScale --- modules/core/src/classes/matrix3.ts | 51 +++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/modules/core/src/classes/matrix3.ts b/modules/core/src/classes/matrix3.ts index 2c622286..c2b77637 100644 --- a/modules/core/src/classes/matrix3.ts +++ b/modules/core/src/classes/matrix3.ts @@ -5,6 +5,7 @@ import {NumericArray, NumericArray9} from '@math.gl/types'; import {Matrix} from './base/matrix'; +import {Vector3} from './vector3'; import {checkVector} from '../lib/validators'; import {vec4_transformMat3} from '../lib/gl-matrix-extras'; @@ -202,6 +203,56 @@ export class Matrix3 extends Matrix { return this.check(); } + /** + * Computes the product of this matrix and a column vector. + * + * @param {Vector3} cartesian The column. + * @param {Vector3} result The object onto which to store the result. + * @returns {Vector3} The modified result parameter. + */ + multiplyByVector(cartesian: Vector3, result?: Vector3): Vector3 { + if (!result) + result = new Vector3() + + const vX = cartesian.x; + const vY = cartesian.y; + const vZ = cartesian.z; + + const x = this[0] * vX + this[3] * vY + this[6] * vZ; + const y = this[1] * vX + this[4] * vY + this[7] * vZ; + const z = this[2] * vX + this[5] * vY + this[8] * vZ; + + result.x = x; + result.y = y; + result.z = z; + + return result; + } + + /** + * Computes the product of this matrix times a (non-uniform) scale, as if the scale were a scale matrix. + * + * @param {Vector3} scale The non-uniform scale on the right-hand side. + * @param {Matrix3} result The object onto which to store the result. + * @returns {Matrix3} The modified result parameter. + */ + multiplyByScale(scale: Vector3, result?: Matrix3): Matrix3 { + if (!result) + result = new Matrix3() + + result[0] = this[0] * scale.x; + result[1] = this[1] * scale.x; + result[2] = this[2] * scale.x; + result[3] = this[3] * scale.y; + result[4] = this[4] * scale.y; + result[5] = this[5] * scale.y; + result[6] = this[6] * scale.z; + result[7] = this[7] * scale.z; + result[8] = this[8] * scale.z; + + return result; + } + rotate(radians: number): this { mat3_rotate(this, this, radians); return this.check(); From 88cd8d6e8f6f0a97745b617d80823d4dcbc497cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schn=C3=BCrer?= Date: Thu, 26 Dec 2024 19:22:02 +0000 Subject: [PATCH 05/13] ported methods from Cesium: - fromRegion - fromPlaneExtents --- .../bounding-volumes/oriented-bounding-box.ts | 255 +++++++++++++++++- 1 file changed, 253 insertions(+), 2 deletions(-) diff --git a/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts b/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts index d5fa7c91..5ba1a82b 100644 --- a/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts +++ b/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts @@ -6,10 +6,16 @@ // See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md import {NumericArray} from '@math.gl/types'; -import {Vector3, Matrix3, Matrix4, Quaternion} from '@math.gl/core'; +import {Vector2, Vector3, Matrix3, Matrix4, Quaternion, + degrees, _MathUtils} from '@math.gl/core'; + +// @ts-ignore +// eslint-disable-next-line import/no-unresolved +import {Ellipsoid, Rectangle} from '@math.gl/geospatial'; import type {BoundingVolume} from './bounding-volume'; import {BoundingSphere} from './bounding-sphere'; -import type {Plane} from '../plane'; +import {Plane} from '../plane'; +import {EllipsoidTangentPlane} from '../ellipsoid-tangent-plane' import {INTERSECTION} from '../../constants'; const scratchVector3 = new Vector3(); @@ -20,6 +26,36 @@ const scratchVectorW = new Vector3(); const scratchCorner = new Vector3(); const scratchToCenter = new Vector3(); +const scratchTangentPoint = new Vector3(); +const scratchPerimeterCartographicNC = new Vector3(); +const scratchPerimeterCartographicNW = new Vector3(); +const scratchPerimeterCartographicCW = new Vector3(); +const scratchPerimeterCartographicSW = new Vector3(); +const scratchPerimeterCartographicSC = new Vector3(); +const scratchPerimeterCartesianNC = new Vector3(); +const scratchPerimeterCartesianNW = new Vector3(); +const scratchPerimeterCartesianCW = new Vector3(); +const scratchPerimeterCartesianSW = new Vector3(); +const scratchPerimeterCartesianSC = new Vector3(); +const scratchPerimeterProjectedNC = new Vector2(); +const scratchPerimeterProjectedNW = new Vector2(); +const scratchPerimeterProjectedCW = new Vector2(); +const scratchPerimeterProjectedSW = new Vector2(); +const scratchPerimeterProjectedSC = new Vector2(); + +const scratchPlane = new Plane(); +const scratchPlaneOrigin = new Vector3(); +const scratchPlaneNormal = new Vector3(); +const scratchPlaneXAxis = new Vector3(); +const scratchHorizonCartesian = new Vector3(); +const scratchHorizonProjected = new Vector2(); +const scratchMaxY = new Vector3(); +const scratchMinY = new Vector3(); +const scratchZ = new Vector3(); + +const VECTOR3_UNIT_X = new Vector3(1, 0, 0); +const VECTOR3_UNIT_Z = new Vector3(0, 0, 1); + const MATRIX3 = { COLUMN0ROW0: 0, COLUMN0ROW1: 1, @@ -98,6 +134,178 @@ export class OrientedBoundingBox implements BoundingVolume { return this; } + /** + * Computes an OrientedBoundingBox that bounds a region on the surface of the WGS84 ellipsoid. + * There are no guarantees about the orientation of the bounding box. + * + * @param {number[]} region The cartographic region ([west, south, east, north, minimum height, maximum height]) + * on the surface of the ellipsoid. + * @returns {OrientedBoundingBox} The modified result parameter or a new OrientedBoundingBox instance if none was provided. + */ + // eslint-disable-next-line max-statements + fromRegion(region: number[]): OrientedBoundingBox { + const [west, south, east, north, minimumHeight, maximumHeight] = region; + + const northDeg = degrees(north); + const southDeg = degrees(south); + + let maxX: number, maxY: number, maxZ: number, minX: number, minY: number, minZ: number, plane: Plane; + + const rectangle = new Rectangle(west, south, east, north); + const tangentPoint = Rectangle.center(rectangle, scratchTangentPoint); + const tangentPointCartographic = new Vector3([degrees(tangentPoint.x), degrees(tangentPoint.y), 0.0]); + + const lonCenter = tangentPoint.x; + const lonCenterDeg = tangentPointCartographic.x; + + if (rectangle.width <= _MathUtils.PI) { + const westDeg = degrees(west); + + const tangentPoint = Ellipsoid.WGS84.cartographicToCartesian(tangentPointCartographic); + const ellipsoidTangentPlane = new EllipsoidTangentPlane(tangentPoint); + + const latCenter = southDeg < 0.0 && northDeg > 0.0 ? 0.0 : tangentPointCartographic.y; + + const perimeterCartographicNC = scratchPerimeterCartographicNC.copy([lonCenterDeg, northDeg, maximumHeight]); + const perimeterCartographicNW = scratchPerimeterCartographicNW.copy([westDeg, northDeg, maximumHeight]); + const perimeterCartographicCW = scratchPerimeterCartographicCW.copy([westDeg, latCenter, maximumHeight]); + const perimeterCartographicSW = scratchPerimeterCartographicSW.copy([westDeg, southDeg, maximumHeight]); + const perimeterCartographicSC = scratchPerimeterCartographicSC.copy([lonCenterDeg, southDeg, maximumHeight]); + + const perimeterCartesianNC = Ellipsoid.WGS84.cartographicToCartesian( + perimeterCartographicNC, + scratchPerimeterCartesianNC, + ); + let perimeterCartesianNW = Ellipsoid.WGS84.cartographicToCartesian( + perimeterCartographicNW, + scratchPerimeterCartesianNW, + ); + const perimeterCartesianCW = Ellipsoid.WGS84.cartographicToCartesian( + perimeterCartographicCW, + scratchPerimeterCartesianCW, + ); + let perimeterCartesianSW = Ellipsoid.WGS84.cartographicToCartesian( + perimeterCartographicSW, + scratchPerimeterCartesianSW, + ); + const perimeterCartesianSC = Ellipsoid.WGS84.cartographicToCartesian( + perimeterCartographicSC, + scratchPerimeterCartesianSC, + ); + + const perimeterProjectedNC = ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianNC, + scratchPerimeterProjectedNC, + ); + const perimeterProjectedNW = ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianNW, + scratchPerimeterProjectedNW, + ); + const perimeterProjectedCW = ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianCW, + scratchPerimeterProjectedCW, + ); + const perimeterProjectedSW = ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianSW, + scratchPerimeterProjectedSW, + ); + const perimeterProjectedSC = ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianSC, + scratchPerimeterProjectedSC, + ); + + minX = Math.min( + perimeterProjectedNW.x, + perimeterProjectedCW.x, + perimeterProjectedSW.x, + ); + maxX = -minX; + + maxY = Math.max(perimeterProjectedNW.y, perimeterProjectedNC.y); + minY = Math.min(perimeterProjectedSW.y, perimeterProjectedSC.y); + + perimeterCartographicNW.z = perimeterCartographicSW.z = + minimumHeight; + perimeterCartesianNW = Ellipsoid.WGS84.cartographicToCartesian( + perimeterCartographicNW, + scratchPerimeterCartesianNW, + ); + perimeterCartesianSW = Ellipsoid.WGS84.cartographicToCartesian( + perimeterCartographicSW, + scratchPerimeterCartesianSW, + ); + + plane = ellipsoidTangentPlane.plane; + minZ = Math.min( + plane.getPointDistance(perimeterCartesianNW), + plane.getPointDistance(perimeterCartesianSW), + ); + maxZ = maximumHeight; + + return fromPlaneExtents( + ellipsoidTangentPlane.origin, + ellipsoidTangentPlane.xAxis, + ellipsoidTangentPlane.yAxis, + ellipsoidTangentPlane.zAxis, + minX, + maxX, + minY, + maxY, + minZ, + maxZ, + this + ); + } + + const eastDeg = degrees(east); + + const fullyAboveEquator = south > 0.0; + const fullyBelowEquator = north < 0.0; + const latitudeNearestToEquator = fullyAboveEquator ? southDeg : fullyBelowEquator ? northDeg : 0.0; + + const planeOrigin = Ellipsoid.WGS84.cartographicToCartesian( + [lonCenterDeg, latitudeNearestToEquator, maximumHeight], scratchPlaneOrigin + ); + planeOrigin.z = 0.0; + const isPole = Math.abs(planeOrigin.x) < _MathUtils.EPSILON10 && Math.abs(planeOrigin.y) < _MathUtils.EPSILON10; + const planeNormal = !isPole ? scratchPlaneNormal.copy(planeOrigin).normalize() : VECTOR3_UNIT_X; + const planeYAxis = VECTOR3_UNIT_Z; + const planeXAxis = scratchPlaneXAxis.copy(planeNormal).cross(planeYAxis); + plane = scratchPlane.fromPointNormal(planeOrigin, planeNormal); + + const horizonCartesian = Ellipsoid.WGS84.cartographicToCartesian([ + degrees(lonCenter + _MathUtils.PI_OVER_TWO), latitudeNearestToEquator, maximumHeight], scratchHorizonCartesian + ); + const projectedPoint = plane.projectPointOntoPlane(horizonCartesian, scratchHorizonProjected) as Vector3; + maxX = projectedPoint.dot(planeXAxis); + minX = -maxX; + + maxY = Ellipsoid.WGS84.cartographicToCartesian( + [0.0, northDeg, fullyBelowEquator ? minimumHeight : maximumHeight], scratchMaxY).z; + minY = Ellipsoid.WGS84.cartographicToCartesian( + [0.0, southDeg, fullyAboveEquator ? minimumHeight : maximumHeight], scratchMinY).z; + + const farZ = Ellipsoid.WGS84.cartographicToCartesian( + [eastDeg, latitudeNearestToEquator, maximumHeight], scratchZ); + + minZ = plane.getPointDistance(farZ); + maxZ = 0.0; + + return fromPlaneExtents( + planeOrigin, + planeXAxis, + planeYAxis, + planeNormal, + minX, + maxX, + minY, + maxY, + minZ, + maxZ, + this + ); + } + /** Duplicates a OrientedBoundingBox instance. */ clone(): OrientedBoundingBox { return new OrientedBoundingBox(this.center, this.halfAxes); @@ -350,3 +558,46 @@ export class OrientedBoundingBox implements BoundingVolume { throw new Error('not implemented'); } } + +const scratchScale = new Vector3(); +const scratchHalfAxes = new Matrix3(); + +/** helper function for OrientedBoundingBox.fromRegion() */ +// eslint-disable-next-line max-params +function fromPlaneExtents( + planeOrigin: Vector3, + planeXAxis: Vector3, + planeYAxis: Vector3, + planeZAxis: Vector3, + minimumX: number, + maximumX: number, + minimumY: number, + maximumY: number, + minimumZ: number, + maximumZ: number, + result: OrientedBoundingBox +) { + if (!result) + result = new OrientedBoundingBox() + + const halfAxes = new Matrix3(); + halfAxes.setColumn(0, planeXAxis); + halfAxes.setColumn(1, planeYAxis); + halfAxes.setColumn(2, planeZAxis); + + let centerOffset = scratchOffset; + centerOffset.x = (minimumX + maximumX) / 2.0; + centerOffset.y = (minimumY + maximumY) / 2.0; + centerOffset.z = (minimumZ + maximumZ) / 2.0; + centerOffset = halfAxes.multiplyByVector(centerOffset, scratchOffset); + + const scale = scratchScale; + scale.x = (maximumX - minimumX) / 2.0; + scale.y = (maximumY - minimumY) / 2.0; + scale.z = (maximumZ - minimumZ) / 2.0; + + result.center = planeOrigin.add(centerOffset); + result.halfAxes = halfAxes.multiplyByScale(scale, scratchHalfAxes); + + return result; +} From d83cd381a961ae54ac0ef110ca9e19162124c30d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schn=C3=BCrer?= Date: Thu, 26 Dec 2024 19:22:22 +0000 Subject: [PATCH 06/13] partly ported class 'Rectangle' from Cesium --- modules/geospatial/src/rectangle.ts | 68 +++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) create mode 100644 modules/geospatial/src/rectangle.ts diff --git a/modules/geospatial/src/rectangle.ts b/modules/geospatial/src/rectangle.ts new file mode 100644 index 00000000..402c14cb --- /dev/null +++ b/modules/geospatial/src/rectangle.ts @@ -0,0 +1,68 @@ +// math.gl +// SPDX-License-Identifier: MIT and Apache-2.0 +// Copyright (c) vis.gl contributors + +// This file is derived from the Cesium library under Apache 2 license +// See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md + +import {Vector3, negativePiToPi, _MathUtils} from '@math.gl/core' + +/** A two dimensional region specified as longitude and latitude coordinates. */ +export class Rectangle { + west: number; + south: number; + east: number; + north: number; + + /** + * Creates a new two dimensional region specified as longitude and latitude coordinates. + * + * @param {number} [west=0.0] The westernmost longitude, in radians, in the range [-Pi, Pi]. + * @param {number} [south=0.0] The southernmost latitude, in radians, in the range [-Pi/2, Pi/2]. + * @param {number} [east=0.0] The easternmost longitude, in radians, in the range [-Pi, Pi]. + * @param {number} [north=0.0] The northernmost latitude, in radians, in the range [-Pi/2, Pi/2]. + */ + constructor(west: number, south: number, east: number, north: number) { + this.west = west; + this.south = south; + this.east = east; + this.north = north; + } + + /** + * Computes the center of a rectangle. + * + * @param {Rectangle} rectangle The rectangle for which to find the center + * @param {Cartographic} [result] The object onto which to store the result. + * @returns {Cartographic} The modified result parameter or a new Cartographic instance if none was provided. + */ + static center(rectangle: Rectangle, result?: Vector3) { + if (!result) + result = new Vector3() + + let east = rectangle.east; + const west = rectangle.west; + + if (east < west) { + east += _MathUtils.TWO_PI; + } + + const longitude = negativePiToPi((west + east) * 0.5); + const latitude = (rectangle.south + rectangle.north) * 0.5; + + result.x = longitude; + result.y = latitude; + result.z = 0.0; + return result; + }; + + get width() { + let east = this.east; + const west = this.west; + + if (east < west) { + east += _MathUtils.TWO_PI; + } + return east - west; + } +} \ No newline at end of file From ef97e409c9f1febecc78d15047352feef112dba2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schn=C3=BCrer?= Date: Thu, 26 Dec 2024 19:22:57 +0000 Subject: [PATCH 07/13] partly ported class 'Ray' from Cesium --- modules/culling/src/lib/ray.ts | 35 ++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 modules/culling/src/lib/ray.ts diff --git a/modules/culling/src/lib/ray.ts b/modules/culling/src/lib/ray.ts new file mode 100644 index 00000000..2ba4518b --- /dev/null +++ b/modules/culling/src/lib/ray.ts @@ -0,0 +1,35 @@ +// math.gl +// SPDX-License-Identifier: MIT and Apache-2.0 +// Copyright (c) vis.gl contributors + +// This file is derived from the Cesium library under Apache 2 license +// See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md + +import {Vector3} from '@math.gl/core'; + +/* Represents a ray that extends infinitely from the provided origin in the provided direction. */ +export class Ray { + origin: Vector3; + direction: Vector3; + + /** + * Creates a new ray that extends infinitely from the provided origin in the provided direction. + * + * @param {Vector3} [origin=Vector3] The origin of the ray. + * @param {Vector3} [direction=Vector3] The direction of the ray. + */ + constructor(origin?: Vector3, direction?: Vector3) { + if (origin) + origin = origin.clone(); + else + origin = new Vector3(); + + if (direction) + direction = direction.clone().normalize(); + else + direction = new Vector3(); + + this.origin = origin + this.direction = direction; + } +} \ No newline at end of file From 815f081b1f683b8cabea4355fca6bf7ff3ea643f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schn=C3=BCrer?= Date: Thu, 26 Dec 2024 19:23:43 +0000 Subject: [PATCH 08/13] ported method from Cesium: - intersectWithRay (originally in IntersectionTests) --- modules/culling/src/lib/plane.ts | 33 +++++++++++++++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/modules/culling/src/lib/plane.ts b/modules/culling/src/lib/plane.ts index 3547d72d..e004a454 100644 --- a/modules/culling/src/lib/plane.ts +++ b/modules/culling/src/lib/plane.ts @@ -6,7 +6,8 @@ // See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md /* eslint-disable */ -import {Vector3, equals, assert, NumericArray} from '@math.gl/core'; +import {Vector3, equals, assert, NumericArray, _MathUtils} from '@math.gl/core'; +import {Ray} from './ray'; const scratchPosition = new Vector3(); const scratchNormal = new Vector3(); @@ -86,4 +87,34 @@ export class Plane { return scratchPoint.subtract(scaledNormal).to(result); } + + /** + * Computes the intersection of a ray and this plane. + * + * @param {Ray} ray The ray. + * @param {Vector3} [result] The object onto which to store the result. + * @returns {Vector3} The intersection point or undefined if there is no intersections. + */ + intersectWithRay(ray: Ray, result?: Vector3): Vector3 { + if (!result) + result = new Vector3() + + const origin = ray.origin; + const direction = ray.direction; + const normal = this.normal; + const denominator = normal.dot(direction); + + if (Math.abs(denominator) < _MathUtils.EPSILON15) { + return undefined; + } + + const t = (-this.distance - normal.dot(origin)) / denominator; + + if (t < 0) { + return undefined; + } + + result = result.copy(direction).multiplyByScalar(t); + return origin.add(result); + } } From 22238fdd499d42d501a796edbbb5316fed270c98 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schn=C3=BCrer?= Date: Thu, 26 Dec 2024 19:24:46 +0000 Subject: [PATCH 09/13] export ported classes --- modules/culling/src/index.ts | 2 ++ modules/geospatial/src/index.ts | 1 + 2 files changed, 3 insertions(+) diff --git a/modules/culling/src/index.ts b/modules/culling/src/index.ts index 5cbc162b..319b2b9a 100644 --- a/modules/culling/src/index.ts +++ b/modules/culling/src/index.ts @@ -9,6 +9,8 @@ export {BoundingSphere} from './lib/bounding-volumes/bounding-sphere'; export {OrientedBoundingBox} from './lib/bounding-volumes/oriented-bounding-box'; export {CullingVolume} from './lib/culling-volume'; export {Plane} from './lib/plane'; +export {Ray} from './lib/ray'; +export {EllipsoidTangentPlane} from './lib/ellipsoid-tangent-plane'; export {PerspectiveOffCenterFrustum as _PerspectiveOffCenterFrustum} from './lib/perspective-off-center-frustum'; export {PerspectiveFrustum as _PerspectiveFrustum} from './lib/perspective-frustum'; diff --git a/modules/geospatial/src/index.ts b/modules/geospatial/src/index.ts index 93a582ac..a02d1447 100644 --- a/modules/geospatial/src/index.ts +++ b/modules/geospatial/src/index.ts @@ -3,4 +3,5 @@ // Copyright (c) vis.gl contributors export {Ellipsoid} from './ellipsoid/ellipsoid'; +export {Rectangle} from './rectangle'; export {isWGS84} from './type-utils'; From 87fd0efe7da04d25b86cc39d2be16a0b2fa9e4bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schn=C3=BCrer?= Date: Thu, 26 Dec 2024 19:25:05 +0000 Subject: [PATCH 10/13] added geospatial module for Rectangle calculations --- modules/culling/package.json | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/modules/culling/package.json b/modules/culling/package.json index b8ca8e4c..c350b139 100644 --- a/modules/culling/package.json +++ b/modules/culling/package.json @@ -43,7 +43,8 @@ ], "dependencies": { "@math.gl/core": "4.1.0-alpha.9", - "@math.gl/types": "4.1.0-alpha.9" + "@math.gl/types": "4.1.0-alpha.9", + "@math.gl/geospatial": "4.1.0-alpha.9" }, "gitHead": "e1a95300cb225a90da6e90333d4adf290f7ba501" } From 1794fd52cea4efcf4cf7b6c1e7dea1e47c6c259a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schn=C3=BCrer?= Date: Mon, 30 Dec 2024 22:55:05 +0000 Subject: [PATCH 11/13] fix: reuse OBBs properties --- .../src/lib/bounding-volumes/oriented-bounding-box.ts | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts b/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts index 5ba1a82b..1c0414cd 100644 --- a/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts +++ b/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts @@ -560,7 +560,6 @@ export class OrientedBoundingBox implements BoundingVolume { } const scratchScale = new Vector3(); -const scratchHalfAxes = new Matrix3(); /** helper function for OrientedBoundingBox.fromRegion() */ // eslint-disable-next-line max-params @@ -577,10 +576,8 @@ function fromPlaneExtents( maximumZ: number, result: OrientedBoundingBox ) { - if (!result) - result = new OrientedBoundingBox() - - const halfAxes = new Matrix3(); + const center = result.center; + const halfAxes = result.halfAxes; halfAxes.setColumn(0, planeXAxis); halfAxes.setColumn(1, planeYAxis); halfAxes.setColumn(2, planeZAxis); @@ -596,8 +593,9 @@ function fromPlaneExtents( scale.y = (maximumY - minimumY) / 2.0; scale.z = (maximumZ - minimumZ) / 2.0; + planeOrigin = center.copy(planeOrigin); result.center = planeOrigin.add(centerOffset); - result.halfAxes = halfAxes.multiplyByScale(scale, scratchHalfAxes); + result.halfAxes = halfAxes.multiplyByScale(scale, halfAxes); return result; } From 80e7cb36bd5c4069a37051a144e24f3e2533801c Mon Sep 17 00:00:00 2001 From: Ib Green Date: Sat, 4 Jan 2025 14:23:21 -0500 Subject: [PATCH 12/13] --amend --- modules/core/src/classes/matrix3.ts | 6 +- modules/core/src/index.ts | 5 +- modules/core/src/lib/common.ts | 49 ++- modules/culling/package.json | 3 +- modules/culling/src/index.ts | 1 - .../bounding-volumes/oriented-bounding-box.ts | 247 +------------- .../src/lib/ellipsoid-tangent-plane.ts | 102 ------ modules/culling/src/lib/plane.ts | 3 +- modules/culling/src/lib/ray.ts | 38 +-- modules/culling/test/lib/ray.spec.ts | 16 + modules/geospatial/package.json | 1 + .../ellipsoid-transform.ts | 0 .../scale-to-geodetic-surface.ts | 0 .../geospatial/src/ellipsoid-tangent-plane.ts | 100 ++++++ .../src/{ellipsoid => }/ellipsoid.ts | 10 +- modules/geospatial/src/index.ts | 5 +- modules/geospatial/src/lng-lat-rectangle.ts | 67 ++++ .../geospatial/src/make-obb-from-region.ts | 321 ++++++++++++++++++ modules/geospatial/src/rectangle.ts | 68 ---- .../geospatial/test/lng-lat-rectangle.spec.ts | 7 + yarn.lock | 3 +- 21 files changed, 579 insertions(+), 473 deletions(-) delete mode 100644 modules/culling/src/lib/ellipsoid-tangent-plane.ts create mode 100644 modules/culling/test/lib/ray.spec.ts rename modules/geospatial/src/{ellipsoid/helpers => ellipsoid-helpers}/ellipsoid-transform.ts (100%) rename modules/geospatial/src/{ellipsoid/helpers => ellipsoid-helpers}/scale-to-geodetic-surface.ts (100%) create mode 100644 modules/geospatial/src/ellipsoid-tangent-plane.ts rename modules/geospatial/src/{ellipsoid => }/ellipsoid.ts (96%) create mode 100644 modules/geospatial/src/lng-lat-rectangle.ts create mode 100644 modules/geospatial/src/make-obb-from-region.ts delete mode 100644 modules/geospatial/src/rectangle.ts create mode 100644 modules/geospatial/test/lng-lat-rectangle.spec.ts diff --git a/modules/core/src/classes/matrix3.ts b/modules/core/src/classes/matrix3.ts index c2b77637..1b8f04d5 100644 --- a/modules/core/src/classes/matrix3.ts +++ b/modules/core/src/classes/matrix3.ts @@ -211,8 +211,7 @@ export class Matrix3 extends Matrix { * @returns {Vector3} The modified result parameter. */ multiplyByVector(cartesian: Vector3, result?: Vector3): Vector3 { - if (!result) - result = new Vector3() + if (!result) result = new Vector3(); const vX = cartesian.x; const vY = cartesian.y; @@ -237,8 +236,7 @@ export class Matrix3 extends Matrix { * @returns {Matrix3} The modified result parameter. */ multiplyByScale(scale: Vector3, result?: Matrix3): Matrix3 { - if (!result) - result = new Matrix3() + if (!result) result = new Matrix3(); result[0] = this[0] * scale.x; result[1] = this[1] * scale.x; diff --git a/modules/core/src/index.ts b/modules/core/src/index.ts index 16ae8c34..ec5932ab 100644 --- a/modules/core/src/index.ts +++ b/modules/core/src/index.ts @@ -48,6 +48,8 @@ export { // math.gl global utility methods config, configure, + safeMod, + normalizeAngle, formatValue, isArray, clone, @@ -58,9 +60,6 @@ export { // math.gl "GLSL"-style functions radians, degrees, - mod, - zeroToTwoPi, - negativePiToPi, sin, cos, tan, diff --git a/modules/core/src/lib/common.ts b/modules/core/src/lib/common.ts index f3fc79d2..c32c91f8 100644 --- a/modules/core/src/lib/common.ts +++ b/modules/core/src/lib/common.ts @@ -123,15 +123,14 @@ export function degrees( return map(radians, (radians) => radians * RADIANS_TO_DEGREES, result); } - /** * The modulo operation that also works for negative dividends. * - * @param {number} m The dividend. - * @param {number} n The divisor. - * @returns {number} The remainder. + * @param m The dividend. + * @param n The divisor. + * @returns The remainder. */ -export function mod(m: number, n: number): number { +export function safeMod(m: number, n: number): number { if (Math.sign(m) === Math.sign(n) && Math.abs(m) < Math.abs(n)) { return m; } @@ -139,21 +138,39 @@ export function mod(m: number, n: number): number { return ((m % n) + n) % n; } +/** + * Produces an angle restricted to its equivalent in a normalized range + * + * @param angle in radians + * @param range 'zero-to-two-pi' - in the range 0 <= angle <= 2PI, | 'negative-pi-to-pi' - -Pi <= angle <= Pi + * @returns The angle in the range [0, TWO_PI] or [-PI, PI].. + */ +export function normalizeAngle( + angle: number, + range: 'zero-to-two-pi' | 'negative-pi-to-pi' +): number { + switch (range) { + case 'negative-pi-to-pi': + return negativePiToPi(angle); + case 'zero-to-two-pi': + return zeroToTwoPi(angle); + default: + return angle; + } +} + /** * Produces an angle in the range 0 <= angle <= 2Pi which is equivalent to the provided angle. * - * @param {number} angle in radians - * @returns {number} The angle in the range [0, TWO_PI]. + * @param angle in radians + * @returns The angle in the range [0, TWO_PI]. */ -export function zeroToTwoPi(angle: number): number { +function zeroToTwoPi(angle: number): number { if (angle >= 0 && angle <= TWO_PI) { return angle; } - const remainder = mod(angle, TWO_PI); - if ( - Math.abs(remainder) < EPSILON14 && - Math.abs(angle) > EPSILON14 - ) { + const remainder = safeMod(angle, TWO_PI); + if (Math.abs(remainder) < EPSILON14 && Math.abs(angle) > EPSILON14) { return TWO_PI; } return remainder; @@ -162,10 +179,10 @@ export function zeroToTwoPi(angle: number): number { /** * Produces an angle in the range -Pi <= angle <= Pi which is equivalent to the provided angle. * - * @param {number} angle in radians - * @returns {number} The angle in the range [-PI, PI]. + * @param angle in radians + * @returns The angle in the range [-PI, PI]. */ -export function negativePiToPi(angle: number): number { +function negativePiToPi(angle: number): number { if (angle >= -PI && angle <= PI) { return angle; } diff --git a/modules/culling/package.json b/modules/culling/package.json index c350b139..b8ca8e4c 100644 --- a/modules/culling/package.json +++ b/modules/culling/package.json @@ -43,8 +43,7 @@ ], "dependencies": { "@math.gl/core": "4.1.0-alpha.9", - "@math.gl/types": "4.1.0-alpha.9", - "@math.gl/geospatial": "4.1.0-alpha.9" + "@math.gl/types": "4.1.0-alpha.9" }, "gitHead": "e1a95300cb225a90da6e90333d4adf290f7ba501" } diff --git a/modules/culling/src/index.ts b/modules/culling/src/index.ts index 319b2b9a..53d6282c 100644 --- a/modules/culling/src/index.ts +++ b/modules/culling/src/index.ts @@ -10,7 +10,6 @@ export {OrientedBoundingBox} from './lib/bounding-volumes/oriented-bounding-box' export {CullingVolume} from './lib/culling-volume'; export {Plane} from './lib/plane'; export {Ray} from './lib/ray'; -export {EllipsoidTangentPlane} from './lib/ellipsoid-tangent-plane'; export {PerspectiveOffCenterFrustum as _PerspectiveOffCenterFrustum} from './lib/perspective-off-center-frustum'; export {PerspectiveFrustum as _PerspectiveFrustum} from './lib/perspective-frustum'; diff --git a/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts b/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts index 1c0414cd..ea015c99 100644 --- a/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts +++ b/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts @@ -6,8 +6,7 @@ // See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md import {NumericArray} from '@math.gl/types'; -import {Vector2, Vector3, Matrix3, Matrix4, Quaternion, - degrees, _MathUtils} from '@math.gl/core'; +import {Vector3, Matrix3, Matrix4, Quaternion, _MathUtils} from '@math.gl/core'; // @ts-ignore // eslint-disable-next-line import/no-unresolved @@ -15,7 +14,6 @@ import {Ellipsoid, Rectangle} from '@math.gl/geospatial'; import type {BoundingVolume} from './bounding-volume'; import {BoundingSphere} from './bounding-sphere'; import {Plane} from '../plane'; -import {EllipsoidTangentPlane} from '../ellipsoid-tangent-plane' import {INTERSECTION} from '../../constants'; const scratchVector3 = new Vector3(); @@ -26,36 +24,6 @@ const scratchVectorW = new Vector3(); const scratchCorner = new Vector3(); const scratchToCenter = new Vector3(); -const scratchTangentPoint = new Vector3(); -const scratchPerimeterCartographicNC = new Vector3(); -const scratchPerimeterCartographicNW = new Vector3(); -const scratchPerimeterCartographicCW = new Vector3(); -const scratchPerimeterCartographicSW = new Vector3(); -const scratchPerimeterCartographicSC = new Vector3(); -const scratchPerimeterCartesianNC = new Vector3(); -const scratchPerimeterCartesianNW = new Vector3(); -const scratchPerimeterCartesianCW = new Vector3(); -const scratchPerimeterCartesianSW = new Vector3(); -const scratchPerimeterCartesianSC = new Vector3(); -const scratchPerimeterProjectedNC = new Vector2(); -const scratchPerimeterProjectedNW = new Vector2(); -const scratchPerimeterProjectedCW = new Vector2(); -const scratchPerimeterProjectedSW = new Vector2(); -const scratchPerimeterProjectedSC = new Vector2(); - -const scratchPlane = new Plane(); -const scratchPlaneOrigin = new Vector3(); -const scratchPlaneNormal = new Vector3(); -const scratchPlaneXAxis = new Vector3(); -const scratchHorizonCartesian = new Vector3(); -const scratchHorizonProjected = new Vector2(); -const scratchMaxY = new Vector3(); -const scratchMinY = new Vector3(); -const scratchZ = new Vector3(); - -const VECTOR3_UNIT_X = new Vector3(1, 0, 0); -const VECTOR3_UNIT_Z = new Vector3(0, 0, 1); - const MATRIX3 = { COLUMN0ROW0: 0, COLUMN0ROW1: 1, @@ -134,178 +102,6 @@ export class OrientedBoundingBox implements BoundingVolume { return this; } - /** - * Computes an OrientedBoundingBox that bounds a region on the surface of the WGS84 ellipsoid. - * There are no guarantees about the orientation of the bounding box. - * - * @param {number[]} region The cartographic region ([west, south, east, north, minimum height, maximum height]) - * on the surface of the ellipsoid. - * @returns {OrientedBoundingBox} The modified result parameter or a new OrientedBoundingBox instance if none was provided. - */ - // eslint-disable-next-line max-statements - fromRegion(region: number[]): OrientedBoundingBox { - const [west, south, east, north, minimumHeight, maximumHeight] = region; - - const northDeg = degrees(north); - const southDeg = degrees(south); - - let maxX: number, maxY: number, maxZ: number, minX: number, minY: number, minZ: number, plane: Plane; - - const rectangle = new Rectangle(west, south, east, north); - const tangentPoint = Rectangle.center(rectangle, scratchTangentPoint); - const tangentPointCartographic = new Vector3([degrees(tangentPoint.x), degrees(tangentPoint.y), 0.0]); - - const lonCenter = tangentPoint.x; - const lonCenterDeg = tangentPointCartographic.x; - - if (rectangle.width <= _MathUtils.PI) { - const westDeg = degrees(west); - - const tangentPoint = Ellipsoid.WGS84.cartographicToCartesian(tangentPointCartographic); - const ellipsoidTangentPlane = new EllipsoidTangentPlane(tangentPoint); - - const latCenter = southDeg < 0.0 && northDeg > 0.0 ? 0.0 : tangentPointCartographic.y; - - const perimeterCartographicNC = scratchPerimeterCartographicNC.copy([lonCenterDeg, northDeg, maximumHeight]); - const perimeterCartographicNW = scratchPerimeterCartographicNW.copy([westDeg, northDeg, maximumHeight]); - const perimeterCartographicCW = scratchPerimeterCartographicCW.copy([westDeg, latCenter, maximumHeight]); - const perimeterCartographicSW = scratchPerimeterCartographicSW.copy([westDeg, southDeg, maximumHeight]); - const perimeterCartographicSC = scratchPerimeterCartographicSC.copy([lonCenterDeg, southDeg, maximumHeight]); - - const perimeterCartesianNC = Ellipsoid.WGS84.cartographicToCartesian( - perimeterCartographicNC, - scratchPerimeterCartesianNC, - ); - let perimeterCartesianNW = Ellipsoid.WGS84.cartographicToCartesian( - perimeterCartographicNW, - scratchPerimeterCartesianNW, - ); - const perimeterCartesianCW = Ellipsoid.WGS84.cartographicToCartesian( - perimeterCartographicCW, - scratchPerimeterCartesianCW, - ); - let perimeterCartesianSW = Ellipsoid.WGS84.cartographicToCartesian( - perimeterCartographicSW, - scratchPerimeterCartesianSW, - ); - const perimeterCartesianSC = Ellipsoid.WGS84.cartographicToCartesian( - perimeterCartographicSC, - scratchPerimeterCartesianSC, - ); - - const perimeterProjectedNC = ellipsoidTangentPlane.projectPointToNearestOnPlane( - perimeterCartesianNC, - scratchPerimeterProjectedNC, - ); - const perimeterProjectedNW = ellipsoidTangentPlane.projectPointToNearestOnPlane( - perimeterCartesianNW, - scratchPerimeterProjectedNW, - ); - const perimeterProjectedCW = ellipsoidTangentPlane.projectPointToNearestOnPlane( - perimeterCartesianCW, - scratchPerimeterProjectedCW, - ); - const perimeterProjectedSW = ellipsoidTangentPlane.projectPointToNearestOnPlane( - perimeterCartesianSW, - scratchPerimeterProjectedSW, - ); - const perimeterProjectedSC = ellipsoidTangentPlane.projectPointToNearestOnPlane( - perimeterCartesianSC, - scratchPerimeterProjectedSC, - ); - - minX = Math.min( - perimeterProjectedNW.x, - perimeterProjectedCW.x, - perimeterProjectedSW.x, - ); - maxX = -minX; - - maxY = Math.max(perimeterProjectedNW.y, perimeterProjectedNC.y); - minY = Math.min(perimeterProjectedSW.y, perimeterProjectedSC.y); - - perimeterCartographicNW.z = perimeterCartographicSW.z = - minimumHeight; - perimeterCartesianNW = Ellipsoid.WGS84.cartographicToCartesian( - perimeterCartographicNW, - scratchPerimeterCartesianNW, - ); - perimeterCartesianSW = Ellipsoid.WGS84.cartographicToCartesian( - perimeterCartographicSW, - scratchPerimeterCartesianSW, - ); - - plane = ellipsoidTangentPlane.plane; - minZ = Math.min( - plane.getPointDistance(perimeterCartesianNW), - plane.getPointDistance(perimeterCartesianSW), - ); - maxZ = maximumHeight; - - return fromPlaneExtents( - ellipsoidTangentPlane.origin, - ellipsoidTangentPlane.xAxis, - ellipsoidTangentPlane.yAxis, - ellipsoidTangentPlane.zAxis, - minX, - maxX, - minY, - maxY, - minZ, - maxZ, - this - ); - } - - const eastDeg = degrees(east); - - const fullyAboveEquator = south > 0.0; - const fullyBelowEquator = north < 0.0; - const latitudeNearestToEquator = fullyAboveEquator ? southDeg : fullyBelowEquator ? northDeg : 0.0; - - const planeOrigin = Ellipsoid.WGS84.cartographicToCartesian( - [lonCenterDeg, latitudeNearestToEquator, maximumHeight], scratchPlaneOrigin - ); - planeOrigin.z = 0.0; - const isPole = Math.abs(planeOrigin.x) < _MathUtils.EPSILON10 && Math.abs(planeOrigin.y) < _MathUtils.EPSILON10; - const planeNormal = !isPole ? scratchPlaneNormal.copy(planeOrigin).normalize() : VECTOR3_UNIT_X; - const planeYAxis = VECTOR3_UNIT_Z; - const planeXAxis = scratchPlaneXAxis.copy(planeNormal).cross(planeYAxis); - plane = scratchPlane.fromPointNormal(planeOrigin, planeNormal); - - const horizonCartesian = Ellipsoid.WGS84.cartographicToCartesian([ - degrees(lonCenter + _MathUtils.PI_OVER_TWO), latitudeNearestToEquator, maximumHeight], scratchHorizonCartesian - ); - const projectedPoint = plane.projectPointOntoPlane(horizonCartesian, scratchHorizonProjected) as Vector3; - maxX = projectedPoint.dot(planeXAxis); - minX = -maxX; - - maxY = Ellipsoid.WGS84.cartographicToCartesian( - [0.0, northDeg, fullyBelowEquator ? minimumHeight : maximumHeight], scratchMaxY).z; - minY = Ellipsoid.WGS84.cartographicToCartesian( - [0.0, southDeg, fullyAboveEquator ? minimumHeight : maximumHeight], scratchMinY).z; - - const farZ = Ellipsoid.WGS84.cartographicToCartesian( - [eastDeg, latitudeNearestToEquator, maximumHeight], scratchZ); - - minZ = plane.getPointDistance(farZ); - maxZ = 0.0; - - return fromPlaneExtents( - planeOrigin, - planeXAxis, - planeYAxis, - planeNormal, - minX, - maxX, - minY, - maxY, - minZ, - maxZ, - this - ); - } - /** Duplicates a OrientedBoundingBox instance. */ clone(): OrientedBoundingBox { return new OrientedBoundingBox(this.center, this.halfAxes); @@ -558,44 +354,3 @@ export class OrientedBoundingBox implements BoundingVolume { throw new Error('not implemented'); } } - -const scratchScale = new Vector3(); - -/** helper function for OrientedBoundingBox.fromRegion() */ -// eslint-disable-next-line max-params -function fromPlaneExtents( - planeOrigin: Vector3, - planeXAxis: Vector3, - planeYAxis: Vector3, - planeZAxis: Vector3, - minimumX: number, - maximumX: number, - minimumY: number, - maximumY: number, - minimumZ: number, - maximumZ: number, - result: OrientedBoundingBox -) { - const center = result.center; - const halfAxes = result.halfAxes; - halfAxes.setColumn(0, planeXAxis); - halfAxes.setColumn(1, planeYAxis); - halfAxes.setColumn(2, planeZAxis); - - let centerOffset = scratchOffset; - centerOffset.x = (minimumX + maximumX) / 2.0; - centerOffset.y = (minimumY + maximumY) / 2.0; - centerOffset.z = (minimumZ + maximumZ) / 2.0; - centerOffset = halfAxes.multiplyByVector(centerOffset, scratchOffset); - - const scale = scratchScale; - scale.x = (maximumX - minimumX) / 2.0; - scale.y = (maximumY - minimumY) / 2.0; - scale.z = (maximumZ - minimumZ) / 2.0; - - planeOrigin = center.copy(planeOrigin); - result.center = planeOrigin.add(centerOffset); - result.halfAxes = halfAxes.multiplyByScale(scale, halfAxes); - - return result; -} diff --git a/modules/culling/src/lib/ellipsoid-tangent-plane.ts b/modules/culling/src/lib/ellipsoid-tangent-plane.ts deleted file mode 100644 index 4230abc1..00000000 --- a/modules/culling/src/lib/ellipsoid-tangent-plane.ts +++ /dev/null @@ -1,102 +0,0 @@ -// math.gl -// SPDX-License-Identifier: MIT and Apache-2.0 -// Copyright (c) vis.gl contributors - -// This file is derived from the Cesium library under Apache 2 license -// See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md - -import {Vector2, Vector3, Matrix4} from '@math.gl/core'; -// @ts-ignore -// eslint-disable-next-line import/no-unresolved -import {Ellipsoid} from '@math.gl/geospatial'; -import {Plane} from './plane'; -import {Ray} from './ray'; - -const scratchOrigin = new Vector3(); -const scratchCart3 = new Vector3(); -const scratchEastNorthUp = new Matrix4(); -const scratchPlane = new Plane(); - -const scratchProjectPointOntoPlaneRay = new Ray(); -const scratchProjectPointOntoPlaneCartesian3 = new Vector3(); -const scratchDirection = new Vector3(); - -/** A plane tangent to the WGS84 ellipsoid at the provided origin */ -export class EllipsoidTangentPlane { - private _origin: Vector3; - private _xAxis: Vector3; - private _yAxis: Vector3; - private _plane: Plane; - - /** - * Creates a new plane tangent to the WGS84 ellipsoid at the provided origin. - * If origin is not on the surface of the ellipsoid, it's surface projection will be used. - * - * @param {Cartesian3} origin The point on the surface of the ellipsoid where the tangent plane touches. - */ - constructor(origin: number[]) { - origin = Ellipsoid.WGS84.scaleToGeodeticSurface(origin, scratchOrigin); - - const eastNorthUp = Ellipsoid.WGS84.eastNorthUpToFixedFrame(origin, scratchEastNorthUp); - - this._origin = origin as Vector3; - this._xAxis = new Vector3(scratchCart3.from(eastNorthUp.getColumn(0))); - this._yAxis = new Vector3(scratchCart3.from(eastNorthUp.getColumn(1))); - const normal = new Vector3(scratchCart3.from(eastNorthUp.getColumn(2))); - - this._plane = scratchPlane.fromPointNormal(origin, normal); - } - - /** - * Computes the projection of the provided 3D position onto the 2D plane, along the plane normal. - * - * @param {Vector3} cartesian The point to project. - * @param {Vector2} [result] The object onto which to store the result. - * @returns {Vector2} The modified result parameter or a new Cartesian2 instance if none was provided. - */ - projectPointToNearestOnPlane (cartesian: Vector3, result?: Vector2): Vector2 { - if (!result) - result = new Vector2(); - - const plane = this._plane; - - const ray = scratchProjectPointOntoPlaneRay; - scratchProjectPointOntoPlaneRay.origin = cartesian; - scratchProjectPointOntoPlaneRay.direction = scratchDirection.copy(plane.normal); - - let intersectionPoint = plane.intersectWithRay(ray, scratchProjectPointOntoPlaneCartesian3); - - if (!intersectionPoint) { - ray.direction = ray.direction.negate(); - intersectionPoint = plane.intersectWithRay(ray, scratchProjectPointOntoPlaneCartesian3); - } - - const v = intersectionPoint.subtract(this._origin); - const x = this._xAxis.dot(v); - const y = this._yAxis.dot(v); - - result.x = x; - result.y = y; - return result; - } - - get plane() { - return this._plane; - } - - get origin() { - return this._origin; - } - - get xAxis() { - return this._xAxis; - } - - get yAxis() { - return this._yAxis; - } - - get zAxis() { - return this._plane.normal; - } -} \ No newline at end of file diff --git a/modules/culling/src/lib/plane.ts b/modules/culling/src/lib/plane.ts index e004a454..2d3edc67 100644 --- a/modules/culling/src/lib/plane.ts +++ b/modules/culling/src/lib/plane.ts @@ -96,8 +96,7 @@ export class Plane { * @returns {Vector3} The intersection point or undefined if there is no intersections. */ intersectWithRay(ray: Ray, result?: Vector3): Vector3 { - if (!result) - result = new Vector3() + if (!result) result = new Vector3(); const origin = ray.origin; const direction = ray.direction; diff --git a/modules/culling/src/lib/ray.ts b/modules/culling/src/lib/ray.ts index 2ba4518b..a2c19969 100644 --- a/modules/culling/src/lib/ray.ts +++ b/modules/culling/src/lib/ray.ts @@ -9,27 +9,23 @@ import {Vector3} from '@math.gl/core'; /* Represents a ray that extends infinitely from the provided origin in the provided direction. */ export class Ray { - origin: Vector3; - direction: Vector3; + origin: Vector3; + direction: Vector3; - /** - * Creates a new ray that extends infinitely from the provided origin in the provided direction. - * - * @param {Vector3} [origin=Vector3] The origin of the ray. - * @param {Vector3} [direction=Vector3] The direction of the ray. - */ - constructor(origin?: Vector3, direction?: Vector3) { - if (origin) - origin = origin.clone(); - else - origin = new Vector3(); + /** + * Creates a new ray that extends infinitely from the provided origin in the provided direction. + * + * @param [origin=Vector3] The origin of the ray. + * @param [direction=Vector3] The direction of the ray. + */ + constructor(origin?: Vector3, direction?: Vector3) { + if (origin) origin = origin.clone(); + else origin = new Vector3(); - if (direction) - direction = direction.clone().normalize(); - else - direction = new Vector3(); + if (direction) direction = direction.clone().normalize(); + else direction = new Vector3(); - this.origin = origin - this.direction = direction; - } -} \ No newline at end of file + this.origin = origin; + this.direction = direction; + } +} diff --git a/modules/culling/test/lib/ray.spec.ts b/modules/culling/test/lib/ray.spec.ts new file mode 100644 index 00000000..609fbb94 --- /dev/null +++ b/modules/culling/test/lib/ray.spec.ts @@ -0,0 +1,16 @@ +// This file is derived from the Cesium math library under Apache 2 license +// See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md + +/* eslint-disable */ +import test from 'tape-promise/tape'; +// import {tapeEquals} from 'test/utils/tape-assertions'; + +import {Vector3, _MathUtils} from '@math.gl/core'; +import {Ray} from '@math.gl/culling'; + +test('Ray#constructs', (t) => { + const ray = new Ray(new Vector3(1, 0, 0)); + t.ok(ray); + t.end(); +}); + diff --git a/modules/geospatial/package.json b/modules/geospatial/package.json index d421f1a0..7dcd409f 100644 --- a/modules/geospatial/package.json +++ b/modules/geospatial/package.json @@ -40,6 +40,7 @@ ], "dependencies": { "@math.gl/core": "4.1.0-alpha.9", + "@math.gl/culling": "4.1.0-alpha.9", "@math.gl/types": "4.1.0-alpha.9" }, "gitHead": "e1a95300cb225a90da6e90333d4adf290f7ba501" diff --git a/modules/geospatial/src/ellipsoid/helpers/ellipsoid-transform.ts b/modules/geospatial/src/ellipsoid-helpers/ellipsoid-transform.ts similarity index 100% rename from modules/geospatial/src/ellipsoid/helpers/ellipsoid-transform.ts rename to modules/geospatial/src/ellipsoid-helpers/ellipsoid-transform.ts diff --git a/modules/geospatial/src/ellipsoid/helpers/scale-to-geodetic-surface.ts b/modules/geospatial/src/ellipsoid-helpers/scale-to-geodetic-surface.ts similarity index 100% rename from modules/geospatial/src/ellipsoid/helpers/scale-to-geodetic-surface.ts rename to modules/geospatial/src/ellipsoid-helpers/scale-to-geodetic-surface.ts diff --git a/modules/geospatial/src/ellipsoid-tangent-plane.ts b/modules/geospatial/src/ellipsoid-tangent-plane.ts new file mode 100644 index 00000000..a13d5320 --- /dev/null +++ b/modules/geospatial/src/ellipsoid-tangent-plane.ts @@ -0,0 +1,100 @@ +// math.gl +// SPDX-License-Identifier: MIT and Apache-2.0 +// Copyright (c) vis.gl contributors + +// This file is derived from the Cesium library under Apache 2 license +// See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md + +import {Vector2, Vector3, Matrix4} from '@math.gl/core'; +// @ts-ignore +// eslint-disable-next-line import/no-unresolved +import {Ellipsoid} from '@math.gl/geospatial'; +import {Plane,Ray} from '@math.gl/culling'; + +const scratchOrigin = new Vector3(); +const scratchCart3 = new Vector3(); +const scratchEastNorthUp = new Matrix4(); +const scratchPlane = new Plane(); + +const scratchProjectPointOntoPlaneRay = new Ray(); +const scratchProjectPointOntoPlaneCartesian3 = new Vector3(); +const scratchDirection = new Vector3(); + +/** A plane tangent to the WGS84 ellipsoid at the provided origin */ +export class EllipsoidTangentPlane { + private _origin: Vector3; + private _xAxis: Vector3; + private _yAxis: Vector3; + private _plane: Plane; + + /** + * Creates a new plane tangent to the WGS84 ellipsoid at the provided origin. + * If origin is not on the surface of the ellipsoid, it's surface projection will be used. + * + * @param {Cartesian3} origin The point on the surface of the ellipsoid where the tangent plane touches. + */ + constructor(origin: number[]) { + origin = Ellipsoid.WGS84.scaleToGeodeticSurface(origin, scratchOrigin); + + const eastNorthUp = Ellipsoid.WGS84.eastNorthUpToFixedFrame(origin, scratchEastNorthUp); + + this._origin = origin as Vector3; + this._xAxis = new Vector3(scratchCart3.from(eastNorthUp.getColumn(0))); + this._yAxis = new Vector3(scratchCart3.from(eastNorthUp.getColumn(1))); + const normal = new Vector3(scratchCart3.from(eastNorthUp.getColumn(2))); + + this._plane = scratchPlane.fromPointNormal(origin, normal); + } + + /** + * Computes the projection of the provided 3D position onto the 2D plane, along the plane normal. + * + * @param {Vector3} cartesian The point to project. + * @param {Vector2} [result] The object onto which to store the result. + * @returns {Vector2} The modified result parameter or a new Cartesian2 instance if none was provided. + */ + projectPointToNearestOnPlane(cartesian: Vector3, result?: Vector2): Vector2 { + if (!result) result = new Vector2(); + + const plane = this._plane; + + const ray = scratchProjectPointOntoPlaneRay; + scratchProjectPointOntoPlaneRay.origin = cartesian; + scratchProjectPointOntoPlaneRay.direction = scratchDirection.copy(plane.normal); + + let intersectionPoint = plane.intersectWithRay(ray, scratchProjectPointOntoPlaneCartesian3); + + if (!intersectionPoint) { + ray.direction = ray.direction.negate(); + intersectionPoint = plane.intersectWithRay(ray, scratchProjectPointOntoPlaneCartesian3); + } + + const v = intersectionPoint.subtract(this._origin); + const x = this._xAxis.dot(v); + const y = this._yAxis.dot(v); + + result.x = x; + result.y = y; + return result; + } + + get plane() { + return this._plane; + } + + get origin() { + return this._origin; + } + + get xAxis() { + return this._xAxis; + } + + get yAxis() { + return this._yAxis; + } + + get zAxis() { + return this._plane.normal; + } +} diff --git a/modules/geospatial/src/ellipsoid/ellipsoid.ts b/modules/geospatial/src/ellipsoid.ts similarity index 96% rename from modules/geospatial/src/ellipsoid/ellipsoid.ts rename to modules/geospatial/src/ellipsoid.ts index fee125e3..73b535ff 100644 --- a/modules/geospatial/src/ellipsoid/ellipsoid.ts +++ b/modules/geospatial/src/ellipsoid.ts @@ -8,12 +8,12 @@ /* eslint-disable */ import {Vector3, Matrix4, assert, equals, _MathUtils, NumericArray, vec3} from '@math.gl/core'; -import {WGS84_RADIUS_X, WGS84_RADIUS_Y, WGS84_RADIUS_Z} from '../constants'; -import {fromCartographicToRadians, toCartographicFromRadians} from '../type-utils'; +import {WGS84_RADIUS_X, WGS84_RADIUS_Y, WGS84_RADIUS_Z} from './constants'; +import {fromCartographicToRadians, toCartographicFromRadians} from './type-utils'; -import type {AxisDirection} from './helpers/ellipsoid-transform'; -import {localFrameToFixedFrame} from './helpers/ellipsoid-transform'; -import {scaleToGeodeticSurface} from './helpers/scale-to-geodetic-surface'; +import type {AxisDirection} from './ellipsoid-helpers/ellipsoid-transform'; +import {localFrameToFixedFrame} from './ellipsoid-helpers/ellipsoid-transform'; +import {scaleToGeodeticSurface} from './ellipsoid-helpers/scale-to-geodetic-surface'; const scratchVector = new Vector3(); const scratchNormal = new Vector3(); diff --git a/modules/geospatial/src/index.ts b/modules/geospatial/src/index.ts index a02d1447..09fabb3f 100644 --- a/modules/geospatial/src/index.ts +++ b/modules/geospatial/src/index.ts @@ -2,6 +2,7 @@ // SPDX-License-Identifier: MIT // Copyright (c) vis.gl contributors -export {Ellipsoid} from './ellipsoid/ellipsoid'; -export {Rectangle} from './rectangle'; +export {Ellipsoid} from './ellipsoid'; +export {EllipsoidTangentPlane} from './ellipsoid-tangent-plane'; +export {LngLatRectangle} from './lng-lat-rectangle'; export {isWGS84} from './type-utils'; diff --git a/modules/geospatial/src/lng-lat-rectangle.ts b/modules/geospatial/src/lng-lat-rectangle.ts new file mode 100644 index 00000000..8e4d4a06 --- /dev/null +++ b/modules/geospatial/src/lng-lat-rectangle.ts @@ -0,0 +1,67 @@ +// math.gl +// SPDX-License-Identifier: MIT and Apache-2.0 +// Copyright (c) vis.gl contributors + +// This file is derived from the Cesium library under Apache 2 license +// See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md + +import {Vector3, normalizeAngle, _MathUtils} from '@math.gl/core'; + +/** A two dimensional region specified as longitude and latitude coordinates. */ +export class LngLatRectangle { + west: number; + south: number; + east: number; + north: number; + + /** + * Creates a new two dimensional region specified as longitude and latitude coordinates. + * + * @param west The westernmost longitude, in radians, in the range [-Pi, Pi]. + * @param south The southernmost latitude, in radians, in the range [-Pi/2, Pi/2]. + * @param east The easternmost longitude, in radians, in the range [-Pi, Pi]. + * @param north The northernmost latitude, in radians, in the range [-Pi/2, Pi/2]. + */ + constructor(west: number, south: number, east: number, north: number) { + this.west = west; + this.south = south; + this.east = east; + this.north = north; + } + + /** + * Computes the center of a rectangle. + * + * @param rectangle The rectangle for which to find the center + * @param [result] The object onto which to store the result. + * @returns The modified result parameter or a new Cartographic instance if none was provided. + */ + static center(rectangle: LngLatRectangle, result?: Vector3): Vector3 { + if (!result) result = new Vector3(); + + let east = rectangle.east; + const west = rectangle.west; + + if (east < west) { + east += _MathUtils.TWO_PI; + } + + const longitude = normalizeAngle((west + east) * 0.5, 'negative-pi-to-pi'); + const latitude = (rectangle.south + rectangle.north) * 0.5; + + result.x = longitude; + result.y = latitude; + result.z = 0.0; + return result; + } + + get width() { + let east = this.east; + const west = this.west; + + if (east < west) { + east += _MathUtils.TWO_PI; + } + return east - west; + } +} diff --git a/modules/geospatial/src/make-obb-from-region.ts b/modules/geospatial/src/make-obb-from-region.ts new file mode 100644 index 00000000..a11f6b39 --- /dev/null +++ b/modules/geospatial/src/make-obb-from-region.ts @@ -0,0 +1,321 @@ +// math.gl +// SPDX-License-Identifier: MIT and Apache-2.0 +// Copyright (c) vis.gl contributors + +// This file is derived from the Cesium math library under Apache 2 license +// See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md + +import { Vector2, Vector3, degrees, _MathUtils } from "@math.gl/core"; +import { OrientedBoundingBox, Plane } from "@math.gl/culling"; +// @ts-ignore +// eslint-disable-next-line import/no-unresolved +import { LngLatRectangle } from "./lng-lat-rectangle"; +import { Ellipsoid } from "./ellipsoid"; +import { EllipsoidTangentPlane } from "./ellipsoid-tangent-plane"; + +const scratchOffset = new Vector3(); + +const scratchTangentPoint = new Vector3(); +const scratchPerimeterCartographicNC = new Vector3(); +const scratchPerimeterCartographicNW = new Vector3(); +const scratchPerimeterCartographicCW = new Vector3(); +const scratchPerimeterCartographicSW = new Vector3(); +const scratchPerimeterCartographicSC = new Vector3(); +const scratchPerimeterCartesianNC = new Vector3(); +const scratchPerimeterCartesianNW = new Vector3(); +const scratchPerimeterCartesianCW = new Vector3(); +const scratchPerimeterCartesianSW = new Vector3(); +const scratchPerimeterCartesianSC = new Vector3(); +const scratchPerimeterProjectedNC = new Vector2(); +const scratchPerimeterProjectedNW = new Vector2(); +const scratchPerimeterProjectedCW = new Vector2(); +const scratchPerimeterProjectedSW = new Vector2(); +const scratchPerimeterProjectedSC = new Vector2(); + +const scratchPlane = new Plane(); +const scratchPlaneOrigin = new Vector3(); +const scratchPlaneNormal = new Vector3(); +const scratchPlaneXAxis = new Vector3(); +const scratchHorizonCartesian = new Vector3(); +const scratchHorizonProjected = new Vector2(); +const scratchMaxY = new Vector3(); +const scratchMinY = new Vector3(); +const scratchZ = new Vector3(); + +const VECTOR3_UNIT_X = new Vector3(1, 0, 0); +const VECTOR3_UNIT_Z = new Vector3(0, 0, 1); + +/** + * Computes an OrientedBoundingBox that bounds a region on the surface of the WGS84 ellipsoid. + * There are no guarantees about the orientation of the bounding box. + * + * @param {number[]} region The cartographic region ([west, south, east, north, minimum height, maximum height]) + * on the surface of the ellipsoid. + * @returns {OrientedBoundingBox} The modified result parameter or a new OrientedBoundingBox instance if none was provided. + */ +// eslint-disable-next-line max-statements +export function makeOBBfromRegion(region: number[]): OrientedBoundingBox { + const obb = new OrientedBoundingBox(); + + const [west, south, east, north, minimumHeight, maximumHeight] = region; + + const northDeg = degrees(north); + const southDeg = degrees(south); + + let maxX: number; + let maxY: number; + let maxZ: number; + let minX: number; + let minY: number; + let minZ: number; + let plane: Plane; + + const rectangle = new LngLatRectangle(west, south, east, north); + const tangentPoint = LngLatRectangle.center(rectangle, scratchTangentPoint); + const tangentPointCartographic = new Vector3([ + degrees(tangentPoint.x), + degrees(tangentPoint.y), + 0.0, + ]); + + const lonCenter = tangentPoint.x; + const lonCenterDeg = tangentPointCartographic.x; + + if (rectangle.width <= _MathUtils.PI) { + const westDeg = degrees(west); + + const tangentPoint = Ellipsoid.WGS84.cartographicToCartesian( + tangentPointCartographic + ); + const ellipsoidTangentPlane = new EllipsoidTangentPlane(tangentPoint); + + const latCenter = + southDeg < 0.0 && northDeg > 0.0 ? 0.0 : tangentPointCartographic.y; + + const perimeterCartographicNC = scratchPerimeterCartographicNC.copy([ + lonCenterDeg, + northDeg, + maximumHeight, + ]); + const perimeterCartographicNW = scratchPerimeterCartographicNW.copy([ + westDeg, + northDeg, + maximumHeight, + ]); + const perimeterCartographicCW = scratchPerimeterCartographicCW.copy([ + westDeg, + latCenter, + maximumHeight, + ]); + const perimeterCartographicSW = scratchPerimeterCartographicSW.copy([ + westDeg, + southDeg, + maximumHeight, + ]); + const perimeterCartographicSC = scratchPerimeterCartographicSC.copy([ + lonCenterDeg, + southDeg, + maximumHeight, + ]); + + const perimeterCartesianNC = Ellipsoid.WGS84.cartographicToCartesian( + perimeterCartographicNC, + scratchPerimeterCartesianNC + ); + let perimeterCartesianNW = Ellipsoid.WGS84.cartographicToCartesian( + perimeterCartographicNW, + scratchPerimeterCartesianNW + ); + const perimeterCartesianCW = Ellipsoid.WGS84.cartographicToCartesian( + perimeterCartographicCW, + scratchPerimeterCartesianCW + ); + let perimeterCartesianSW = Ellipsoid.WGS84.cartographicToCartesian( + perimeterCartographicSW, + scratchPerimeterCartesianSW + ); + const perimeterCartesianSC = Ellipsoid.WGS84.cartographicToCartesian( + perimeterCartographicSC, + scratchPerimeterCartesianSC + ); + + const perimeterProjectedNC = + ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianNC, + scratchPerimeterProjectedNC + ); + const perimeterProjectedNW = + ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianNW, + scratchPerimeterProjectedNW + ); + const perimeterProjectedCW = + ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianCW, + scratchPerimeterProjectedCW + ); + const perimeterProjectedSW = + ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianSW, + scratchPerimeterProjectedSW + ); + const perimeterProjectedSC = + ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianSC, + scratchPerimeterProjectedSC + ); + + minX = Math.min( + perimeterProjectedNW.x, + perimeterProjectedCW.x, + perimeterProjectedSW.x + ); + maxX = -minX; + + maxY = Math.max(perimeterProjectedNW.y, perimeterProjectedNC.y); + minY = Math.min(perimeterProjectedSW.y, perimeterProjectedSC.y); + + perimeterCartographicNW.z = perimeterCartographicSW.z = minimumHeight; + perimeterCartesianNW = Ellipsoid.WGS84.cartographicToCartesian( + perimeterCartographicNW, + scratchPerimeterCartesianNW + ); + perimeterCartesianSW = Ellipsoid.WGS84.cartographicToCartesian( + perimeterCartographicSW, + scratchPerimeterCartesianSW + ); + + plane = ellipsoidTangentPlane.plane; + minZ = Math.min( + plane.getPointDistance(perimeterCartesianNW), + plane.getPointDistance(perimeterCartesianSW) + ); + maxZ = maximumHeight; + + return fromPlaneExtents( + ellipsoidTangentPlane.origin, + ellipsoidTangentPlane.xAxis, + ellipsoidTangentPlane.yAxis, + ellipsoidTangentPlane.zAxis, + minX, + maxX, + minY, + maxY, + minZ, + maxZ, + obb + ); + } + + const eastDeg = degrees(east); + + const fullyAboveEquator = south > 0.0; + const fullyBelowEquator = north < 0.0; + const latitudeNearestToEquator = fullyAboveEquator + ? southDeg + : fullyBelowEquator + ? northDeg + : 0.0; + + const planeOrigin = Ellipsoid.WGS84.cartographicToCartesian( + [lonCenterDeg, latitudeNearestToEquator, maximumHeight], + scratchPlaneOrigin + ); + planeOrigin.z = 0.0; + const isPole = + Math.abs(planeOrigin.x) < _MathUtils.EPSILON10 && + Math.abs(planeOrigin.y) < _MathUtils.EPSILON10; + const planeNormal = !isPole + ? scratchPlaneNormal.copy(planeOrigin).normalize() + : VECTOR3_UNIT_X; + const planeYAxis = VECTOR3_UNIT_Z; + const planeXAxis = scratchPlaneXAxis.copy(planeNormal).cross(planeYAxis); + plane = scratchPlane.fromPointNormal(planeOrigin, planeNormal); + + const horizonCartesian = Ellipsoid.WGS84.cartographicToCartesian( + [ + degrees(lonCenter + _MathUtils.PI_OVER_TWO), + latitudeNearestToEquator, + maximumHeight, + ], + scratchHorizonCartesian + ); + const projectedPoint = plane.projectPointOntoPlane( + horizonCartesian, + scratchHorizonProjected + ) as Vector3; + maxX = projectedPoint.dot(planeXAxis); + minX = -maxX; + + maxY = Ellipsoid.WGS84.cartographicToCartesian( + [0.0, northDeg, fullyBelowEquator ? minimumHeight : maximumHeight], + scratchMaxY + ).z; + minY = Ellipsoid.WGS84.cartographicToCartesian( + [0.0, southDeg, fullyAboveEquator ? minimumHeight : maximumHeight], + scratchMinY + ).z; + + const farZ = Ellipsoid.WGS84.cartographicToCartesian( + [eastDeg, latitudeNearestToEquator, maximumHeight], + scratchZ + ); + + minZ = plane.getPointDistance(farZ); + maxZ = 0.0; + + return fromPlaneExtents( + planeOrigin, + planeXAxis, + planeYAxis, + planeNormal, + minX, + maxX, + minY, + maxY, + minZ, + maxZ, + obb + ); +} + +const scratchScale = new Vector3(); + +/** helper function for OrientedBoundingBox.fromRegion() */ +// eslint-disable-next-line max-params +function fromPlaneExtents( + planeOrigin: Vector3, + planeXAxis: Vector3, + planeYAxis: Vector3, + planeZAxis: Vector3, + minimumX: number, + maximumX: number, + minimumY: number, + maximumY: number, + minimumZ: number, + maximumZ: number, + result: OrientedBoundingBox +) { + const center = result.center; + const halfAxes = result.halfAxes; + halfAxes.setColumn(0, planeXAxis); + halfAxes.setColumn(1, planeYAxis); + halfAxes.setColumn(2, planeZAxis); + + let centerOffset = scratchOffset; + centerOffset.x = (minimumX + maximumX) / 2.0; + centerOffset.y = (minimumY + maximumY) / 2.0; + centerOffset.z = (minimumZ + maximumZ) / 2.0; + centerOffset = halfAxes.multiplyByVector(centerOffset, scratchOffset); + + const scale = scratchScale; + scale.x = (maximumX - minimumX) / 2.0; + scale.y = (maximumY - minimumY) / 2.0; + scale.z = (maximumZ - minimumZ) / 2.0; + + planeOrigin = center.copy(planeOrigin); + result.center = planeOrigin.add(centerOffset); + result.halfAxes = halfAxes.multiplyByScale(scale, halfAxes); + + return result; +} diff --git a/modules/geospatial/src/rectangle.ts b/modules/geospatial/src/rectangle.ts deleted file mode 100644 index 402c14cb..00000000 --- a/modules/geospatial/src/rectangle.ts +++ /dev/null @@ -1,68 +0,0 @@ -// math.gl -// SPDX-License-Identifier: MIT and Apache-2.0 -// Copyright (c) vis.gl contributors - -// This file is derived from the Cesium library under Apache 2 license -// See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md - -import {Vector3, negativePiToPi, _MathUtils} from '@math.gl/core' - -/** A two dimensional region specified as longitude and latitude coordinates. */ -export class Rectangle { - west: number; - south: number; - east: number; - north: number; - - /** - * Creates a new two dimensional region specified as longitude and latitude coordinates. - * - * @param {number} [west=0.0] The westernmost longitude, in radians, in the range [-Pi, Pi]. - * @param {number} [south=0.0] The southernmost latitude, in radians, in the range [-Pi/2, Pi/2]. - * @param {number} [east=0.0] The easternmost longitude, in radians, in the range [-Pi, Pi]. - * @param {number} [north=0.0] The northernmost latitude, in radians, in the range [-Pi/2, Pi/2]. - */ - constructor(west: number, south: number, east: number, north: number) { - this.west = west; - this.south = south; - this.east = east; - this.north = north; - } - - /** - * Computes the center of a rectangle. - * - * @param {Rectangle} rectangle The rectangle for which to find the center - * @param {Cartographic} [result] The object onto which to store the result. - * @returns {Cartographic} The modified result parameter or a new Cartographic instance if none was provided. - */ - static center(rectangle: Rectangle, result?: Vector3) { - if (!result) - result = new Vector3() - - let east = rectangle.east; - const west = rectangle.west; - - if (east < west) { - east += _MathUtils.TWO_PI; - } - - const longitude = negativePiToPi((west + east) * 0.5); - const latitude = (rectangle.south + rectangle.north) * 0.5; - - result.x = longitude; - result.y = latitude; - result.z = 0.0; - return result; - }; - - get width() { - let east = this.east; - const west = this.west; - - if (east < west) { - east += _MathUtils.TWO_PI; - } - return east - west; - } -} \ No newline at end of file diff --git a/modules/geospatial/test/lng-lat-rectangle.spec.ts b/modules/geospatial/test/lng-lat-rectangle.spec.ts new file mode 100644 index 00000000..fda35cbc --- /dev/null +++ b/modules/geospatial/test/lng-lat-rectangle.spec.ts @@ -0,0 +1,7 @@ +import test from 'tape-promise/tape'; +import {LngLatRectangle} from '@math.gl/geospatial'; + +test('LngLatRectangle', (t) => { + t.ok(new LngLatRectangle(0, 0, 0, 0)); + t.end(); +}); diff --git a/yarn.lock b/yarn.lock index 84b13758..a7e3a340 100644 --- a/yarn.lock +++ b/yarn.lock @@ -812,6 +812,7 @@ __metadata: resolution: "@math.gl/culling@workspace:modules/culling" dependencies: "@math.gl/core": "npm:4.1.0-alpha.9" + "@math.gl/geospatial": "npm:4.1.0-alpha.9" "@math.gl/types": "npm:4.1.0-alpha.9" languageName: unknown linkType: soft @@ -844,7 +845,7 @@ __metadata: languageName: unknown linkType: soft -"@math.gl/geospatial@workspace:modules/geospatial": +"@math.gl/geospatial@npm:4.1.0-alpha.9, @math.gl/geospatial@workspace:modules/geospatial": version: 0.0.0-use.local resolution: "@math.gl/geospatial@workspace:modules/geospatial" dependencies: From 80496c602dad4f3f81f69bc3dfcd3cacf0266c20 Mon Sep 17 00:00:00 2001 From: Ib Green Date: Sat, 4 Jan 2025 14:31:36 -0500 Subject: [PATCH 13/13] wip --- .../bounding-volumes/oriented-bounding-box.ts | 5 +- modules/culling/test/lib/ray.spec.ts | 1 - .../geospatial/src/ellipsoid-tangent-plane.ts | 6 +- .../geospatial/src/make-obb-from-region.ts | 90 ++++++++----------- modules/geospatial/tsconfig.json | 3 +- 5 files changed, 41 insertions(+), 64 deletions(-) diff --git a/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts b/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts index ea015c99..2e512952 100644 --- a/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts +++ b/modules/culling/src/lib/bounding-volumes/oriented-bounding-box.ts @@ -6,11 +6,8 @@ // See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md import {NumericArray} from '@math.gl/types'; -import {Vector3, Matrix3, Matrix4, Quaternion, _MathUtils} from '@math.gl/core'; +import {Vector3, Matrix3, Matrix4, Quaternion} from '@math.gl/core'; -// @ts-ignore -// eslint-disable-next-line import/no-unresolved -import {Ellipsoid, Rectangle} from '@math.gl/geospatial'; import type {BoundingVolume} from './bounding-volume'; import {BoundingSphere} from './bounding-sphere'; import {Plane} from '../plane'; diff --git a/modules/culling/test/lib/ray.spec.ts b/modules/culling/test/lib/ray.spec.ts index 609fbb94..daf6a394 100644 --- a/modules/culling/test/lib/ray.spec.ts +++ b/modules/culling/test/lib/ray.spec.ts @@ -13,4 +13,3 @@ test('Ray#constructs', (t) => { t.ok(ray); t.end(); }); - diff --git a/modules/geospatial/src/ellipsoid-tangent-plane.ts b/modules/geospatial/src/ellipsoid-tangent-plane.ts index a13d5320..de7c3d75 100644 --- a/modules/geospatial/src/ellipsoid-tangent-plane.ts +++ b/modules/geospatial/src/ellipsoid-tangent-plane.ts @@ -6,10 +6,8 @@ // See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md import {Vector2, Vector3, Matrix4} from '@math.gl/core'; -// @ts-ignore -// eslint-disable-next-line import/no-unresolved -import {Ellipsoid} from '@math.gl/geospatial'; -import {Plane,Ray} from '@math.gl/culling'; +import {Plane, Ray} from '@math.gl/culling'; +import {Ellipsoid} from './ellipsoid'; const scratchOrigin = new Vector3(); const scratchCart3 = new Vector3(); diff --git a/modules/geospatial/src/make-obb-from-region.ts b/modules/geospatial/src/make-obb-from-region.ts index a11f6b39..ec52a800 100644 --- a/modules/geospatial/src/make-obb-from-region.ts +++ b/modules/geospatial/src/make-obb-from-region.ts @@ -5,13 +5,13 @@ // This file is derived from the Cesium math library under Apache 2 license // See LICENSE.md and https://github.com/AnalyticalGraphicsInc/cesium/blob/master/LICENSE.md -import { Vector2, Vector3, degrees, _MathUtils } from "@math.gl/core"; -import { OrientedBoundingBox, Plane } from "@math.gl/culling"; +import {Vector2, Vector3, degrees, _MathUtils} from '@math.gl/core'; +import {OrientedBoundingBox, Plane} from '@math.gl/culling'; // @ts-ignore // eslint-disable-next-line import/no-unresolved -import { LngLatRectangle } from "./lng-lat-rectangle"; -import { Ellipsoid } from "./ellipsoid"; -import { EllipsoidTangentPlane } from "./ellipsoid-tangent-plane"; +import {LngLatRectangle} from './lng-lat-rectangle'; +import {Ellipsoid} from './ellipsoid'; +import {EllipsoidTangentPlane} from './ellipsoid-tangent-plane'; const scratchOffset = new Vector3(); @@ -75,7 +75,7 @@ export function makeOBBfromRegion(region: number[]): OrientedBoundingBox { const tangentPointCartographic = new Vector3([ degrees(tangentPoint.x), degrees(tangentPoint.y), - 0.0, + 0.0 ]); const lonCenter = tangentPoint.x; @@ -84,38 +84,35 @@ export function makeOBBfromRegion(region: number[]): OrientedBoundingBox { if (rectangle.width <= _MathUtils.PI) { const westDeg = degrees(west); - const tangentPoint = Ellipsoid.WGS84.cartographicToCartesian( - tangentPointCartographic - ); + const tangentPoint = Ellipsoid.WGS84.cartographicToCartesian(tangentPointCartographic); const ellipsoidTangentPlane = new EllipsoidTangentPlane(tangentPoint); - const latCenter = - southDeg < 0.0 && northDeg > 0.0 ? 0.0 : tangentPointCartographic.y; + const latCenter = southDeg < 0.0 && northDeg > 0.0 ? 0.0 : tangentPointCartographic.y; const perimeterCartographicNC = scratchPerimeterCartographicNC.copy([ lonCenterDeg, northDeg, - maximumHeight, + maximumHeight ]); const perimeterCartographicNW = scratchPerimeterCartographicNW.copy([ westDeg, northDeg, - maximumHeight, + maximumHeight ]); const perimeterCartographicCW = scratchPerimeterCartographicCW.copy([ westDeg, latCenter, - maximumHeight, + maximumHeight ]); const perimeterCartographicSW = scratchPerimeterCartographicSW.copy([ westDeg, southDeg, - maximumHeight, + maximumHeight ]); const perimeterCartographicSC = scratchPerimeterCartographicSC.copy([ lonCenterDeg, southDeg, - maximumHeight, + maximumHeight ]); const perimeterCartesianNC = Ellipsoid.WGS84.cartographicToCartesian( @@ -139,37 +136,28 @@ export function makeOBBfromRegion(region: number[]): OrientedBoundingBox { scratchPerimeterCartesianSC ); - const perimeterProjectedNC = - ellipsoidTangentPlane.projectPointToNearestOnPlane( - perimeterCartesianNC, - scratchPerimeterProjectedNC - ); - const perimeterProjectedNW = - ellipsoidTangentPlane.projectPointToNearestOnPlane( - perimeterCartesianNW, - scratchPerimeterProjectedNW - ); - const perimeterProjectedCW = - ellipsoidTangentPlane.projectPointToNearestOnPlane( - perimeterCartesianCW, - scratchPerimeterProjectedCW - ); - const perimeterProjectedSW = - ellipsoidTangentPlane.projectPointToNearestOnPlane( - perimeterCartesianSW, - scratchPerimeterProjectedSW - ); - const perimeterProjectedSC = - ellipsoidTangentPlane.projectPointToNearestOnPlane( - perimeterCartesianSC, - scratchPerimeterProjectedSC - ); - - minX = Math.min( - perimeterProjectedNW.x, - perimeterProjectedCW.x, - perimeterProjectedSW.x + const perimeterProjectedNC = ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianNC, + scratchPerimeterProjectedNC + ); + const perimeterProjectedNW = ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianNW, + scratchPerimeterProjectedNW + ); + const perimeterProjectedCW = ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianCW, + scratchPerimeterProjectedCW ); + const perimeterProjectedSW = ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianSW, + scratchPerimeterProjectedSW + ); + const perimeterProjectedSC = ellipsoidTangentPlane.projectPointToNearestOnPlane( + perimeterCartesianSC, + scratchPerimeterProjectedSC + ); + + minX = Math.min(perimeterProjectedNW.x, perimeterProjectedCW.x, perimeterProjectedSW.x); maxX = -minX; maxY = Math.max(perimeterProjectedNW.y, perimeterProjectedNC.y); @@ -225,19 +213,13 @@ export function makeOBBfromRegion(region: number[]): OrientedBoundingBox { const isPole = Math.abs(planeOrigin.x) < _MathUtils.EPSILON10 && Math.abs(planeOrigin.y) < _MathUtils.EPSILON10; - const planeNormal = !isPole - ? scratchPlaneNormal.copy(planeOrigin).normalize() - : VECTOR3_UNIT_X; + const planeNormal = !isPole ? scratchPlaneNormal.copy(planeOrigin).normalize() : VECTOR3_UNIT_X; const planeYAxis = VECTOR3_UNIT_Z; const planeXAxis = scratchPlaneXAxis.copy(planeNormal).cross(planeYAxis); plane = scratchPlane.fromPointNormal(planeOrigin, planeNormal); const horizonCartesian = Ellipsoid.WGS84.cartographicToCartesian( - [ - degrees(lonCenter + _MathUtils.PI_OVER_TWO), - latitudeNearestToEquator, - maximumHeight, - ], + [degrees(lonCenter + _MathUtils.PI_OVER_TWO), latitudeNearestToEquator, maximumHeight], scratchHorizonCartesian ); const projectedPoint = plane.projectPointOntoPlane( diff --git a/modules/geospatial/tsconfig.json b/modules/geospatial/tsconfig.json index b6d00108..d57a67bf 100644 --- a/modules/geospatial/tsconfig.json +++ b/modules/geospatial/tsconfig.json @@ -9,6 +9,7 @@ }, "references": [ {"path": "../types"}, - {"path": "../core"} + {"path": "../core"}, + {"path": "../culling"} ] } \ No newline at end of file