From 4a92ed8377324c62b391b489565420ee04c8fa6a Mon Sep 17 00:00:00 2001 From: German Khokhlov Date: Thu, 15 Oct 2026 00:15:20 +0400 Subject: [PATCH] feat(geo): add src/geo, Cesium's vector maths without Cesium An engine-neutral port of the CesiumJS 1.138 core maths gev uses: angle maths, the WGS84 ellipsoid and geodetic coordinates, 2D/3D vectors, 3x3/4x4 matrices, quaternions and local east-north-up / heading-pitch-roll frames. Formulas and operation order match Cesium, so outputs are bit-identical (golden tests pin the centre, pole, -0 and quaternion branches). - Values are plain-data instances of module-private constructors (own enumerable fields only), so geo's hot functions keep their own V8 hidden classes; matrices are plain column-major arrays Cesium accepts. - Cesium's debug argument checks are ported with the same conditions, order and DeveloperError messages, gated on DEBUG (src/geo/debug.js): on in Node and the dev server, folded away by vite build, as Cesium's Build strips them. A release-parity test compares DEBUG=false geo with Build/Cesium. - The import-direction check forbids cesium and @cesium/* in src/geo; CODE-BOUNDARIES.md and THIRD_PARTY_NOTICES.md document the port. --- THIRD_PARTY_NOTICES.md | 24 +- docs/CODE-BOUNDARIES.md | 8 +- scripts/check-import-directions.mjs | 7 +- src/geo/bitExact.test.mjs | 485 +++++++++++++++++++++++++ src/geo/cartesian2.js | 49 +++ src/geo/cartesian2.test.mjs | 19 + src/geo/cartesian3.js | 389 +++++++++++++++++++++ src/geo/cartesian3.test.mjs | 157 +++++++++ src/geo/cartographic.js | 163 +++++++++ src/geo/cartographic.test.mjs | 142 ++++++++ src/geo/check.js | 61 ++++ src/geo/checks.test.mjs | 486 ++++++++++++++++++++++++++ src/geo/debug.js | 12 + src/geo/developerError.js | 13 + src/geo/ellipsoid.js | 180 ++++++++++ src/geo/ellipsoid.test.mjs | 148 ++++++++ src/geo/index.js | 15 + src/geo/index.test.mjs | 22 ++ src/geo/math.js | 133 +++++++ src/geo/math.test.mjs | 121 +++++++ src/geo/matrix.test.mjs | 152 ++++++++ src/geo/matrix3.js | 121 +++++++ src/geo/matrix4.js | 434 +++++++++++++++++++++++ src/geo/quaternion.js | 129 +++++++ src/geo/releaseParity.test.mjs | 446 +++++++++++++++++++++++ src/geo/scaleToGeodeticSurface.js | 142 ++++++++ src/geo/transforms.js | 240 +++++++++++++ src/geo/transforms.test.mjs | 152 ++++++++ src/tooling/importDirections.test.mjs | 26 ++ 29 files changed, 4470 insertions(+), 6 deletions(-) create mode 100644 src/geo/bitExact.test.mjs create mode 100644 src/geo/cartesian2.js create mode 100644 src/geo/cartesian2.test.mjs create mode 100644 src/geo/cartesian3.js create mode 100644 src/geo/cartesian3.test.mjs create mode 100644 src/geo/cartographic.js create mode 100644 src/geo/cartographic.test.mjs create mode 100644 src/geo/check.js create mode 100644 src/geo/checks.test.mjs create mode 100644 src/geo/debug.js create mode 100644 src/geo/developerError.js create mode 100644 src/geo/ellipsoid.js create mode 100644 src/geo/ellipsoid.test.mjs create mode 100644 src/geo/index.js create mode 100644 src/geo/index.test.mjs create mode 100644 src/geo/math.js create mode 100644 src/geo/math.test.mjs create mode 100644 src/geo/matrix.test.mjs create mode 100644 src/geo/matrix3.js create mode 100644 src/geo/matrix4.js create mode 100644 src/geo/quaternion.js create mode 100644 src/geo/releaseParity.test.mjs create mode 100644 src/geo/scaleToGeodeticSurface.js create mode 100644 src/geo/transforms.js create mode 100644 src/geo/transforms.test.mjs diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md index ac186c25b..13246b67b 100644 --- a/THIRD_PARTY_NOTICES.md +++ b/THIRD_PARTY_NOTICES.md @@ -2,8 +2,8 @@ The project depends on third-party npm packages whose own licenses apply. This file records the packages added for the browser-local SDR feature and the -Street Level viewer; the complete resolved dependency inventory remains in -`package-lock.json`. +Street Level viewer, and the CesiumJS code ported into `src/geo/`; the complete +resolved dependency inventory remains in `package-lock.json`. ## Web RTL-SDR @@ -45,3 +45,23 @@ Street Level viewer; the complete resolved dependency inventory remains in - Source: - License: MIT - License text: + +## CesiumJS (portions ported to src/geo) + +- Package: `cesium` 1.138.0 (`@cesium/engine` 22.3.0); the files in `src/geo/` + are adapted from `@cesium/engine/Source/Core` (`Math.js`, `Ellipsoid.js`, + `scaleToGeodeticSurface.js`, `Cartographic.js`, `Cartesian2.js`, + `Cartesian3.js`, `Matrix3.js`, `Matrix4.js`, `Quaternion.js`, + `Transforms.js`, `Check.js`, `DeveloperError.js`) and say so in their first + line. The changes: the classes' static functions became module functions, + vectors, geodetic positions and quaternions became plain-data instances of + module-private constructors (own enumerable fields only, no methods), + matrices plain arrays, and the ellipsoid a frozen object; the debug + argument checks were ported with Cesium's conditions and messages and, like + Cesium's, are active in development and tests and stripped from production + builds; only the subset gev uses was kept. +- Author: Copyright 2011-2024 CesiumJS Contributors +- Source: +- License: Apache License 2.0 +- License text: `node_modules/@cesium/engine/LICENSE.md` (shipped with the + dependency) and diff --git a/docs/CODE-BOUNDARIES.md b/docs/CODE-BOUNDARIES.md index a0a2325ad..47ed184b9 100644 --- a/docs/CODE-BOUNDARIES.md +++ b/docs/CODE-BOUNDARIES.md @@ -32,6 +32,7 @@ export inventory. Use declared exports rather than reaching into internal files. | `src/data/` | Lifecycle, context and feed state; legacy default-layer facades | | `src/layers//` | Source acquisition, records and Cesium resources with explicit controller/renderer owners | | `src/sources/` | Portable protocols and source contracts; request state belongs to each factory instance | +| `src/geo/` | Engine-neutral WGS84 maths, ported from Cesium, never imports Cesium | | `src/services/` | Supplied application operations; scene construction owns caches and cancellation | | `src/voice/` | Portable action schemas/session, common controls, action execution and separate protocol adapters | | `server/providers/` | Node route factories, process-scoped provider caches and shutdown cleanup | @@ -59,9 +60,10 @@ See [application construction](APPLICATION.md) and the Portable source graphs cannot reach application/rendering, Node, Cesium or browser globals. This includes `sources/*`, dedicated `layers/*/source` exports, flight/military/vessel record and ingestion exports, action schemas, the session -interface, lifecycle and feed state. The browser-global rule reserves platform -names such as `document` and `window` in these modules; it is an architectural -check, not a JavaScript sandbox. Common voice controls cannot depend on a +interface, lifecycle and feed state, and every runtime file in `src/geo/` (the +engine-neutral geometry ported from Cesium). The browser-global rule reserves +platform names such as `document` and `window` in these modules; it is an +architectural check, not a JavaScript sandbox. Common voice controls cannot depend on a Realtime protocol implementation. Negative fixtures cover indirect helpers, self-package imports, symlinks and unreachable files. diff --git a/scripts/check-import-directions.mjs b/scripts/check-import-directions.mjs index ac1a3f732..bdbb7ff28 100644 --- a/scripts/check-import-directions.mjs +++ b/scripts/check-import-directions.mjs @@ -19,6 +19,8 @@ const tests = (file) => const source = (file) => file.startsWith('src/sources/') || /^src\/layers\/[^/]+\/(?:source|bundledSource|flowSource)\.js$/.test(file); +// Engine-neutral geometry ported from Cesium; every runtime file is a portable root. +const geo = (file) => file.startsWith('src/geo/'); const renderer = (file) => /^src\/(?:ui|app|standalone|overlays)\//.test(file) || /^src\/layers\/[^/]+\/(?:index|rendering|snapshotRenderer|overlay|presentation|cards|controls)\.js$/.test( @@ -186,6 +188,7 @@ export function checkImportDirections(root) { for (const file of files.filter((file) => file.startsWith('src/'))) browser(file); const roots = new Map(files.filter(source).map((file) => [file, 'source'])); + for (const file of files.filter(geo)) roots.set(file, 'geo'); for (const [key, value] of Object.entries(pkg.exports || {})) { if (portableExport(key) && typeof value === 'string') roots.set(value.replace(/^\.\//, ''), key); @@ -211,7 +214,9 @@ export function checkImportDirections(root) { if ( builtins.has(specifier) || specifier.startsWith('node:') || - /^(?:cesium|vite)(?:\/|$)/.test(specifier) + /^(?:cesium|vite)(?:\/|$)/.test(specifier) || + // Geo must not reach Cesium through its workspace packages either. + (label === 'geo' && /^@cesium\//.test(specifier)) ) report( file, diff --git a/src/geo/bitExact.test.mjs b/src/geo/bitExact.test.mjs new file mode 100644 index 000000000..3628a3c5a --- /dev/null +++ b/src/geo/bitExact.test.mjs @@ -0,0 +1,485 @@ +// Seeded, deterministic bit-for-bit comparison of src/geo against Cesium. +// Every number is compared with Object.is (so -0 and NaN count), a throw on one +// side must be a throw on the other, and `result` aliasing is exercised the +// same way on both sides. +import { test } from 'node:test'; +import assert from 'node:assert/strict'; +import * as Cesium from 'cesium'; +import * as geo from './index.js'; + +const { + cartesian3: c3, + cartographic, + ellipsoid, + matrix3, + matrix4, + quaternion, + transforms, +} = geo; +const C3 = Cesium.Cartesian3; +const CASES = 2000; +const POLAR_RADIUS = 6356752.3142451793; + +// mulberry32: a small seeded generator, so a failure always reproduces. +function generator(seed) { + return () => { + seed = (seed + 0x6d2b79f5) | 0; + let t = Math.imul(seed ^ (seed >>> 15), 1 | seed); + t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t; + return ((t ^ (t >>> 14)) >>> 0) / 4294967296; + }; +} + +const SPECIAL = [ + 0, + -0, + 1e-15, + -1e-15, + 1, + -1, + 6378137, + -6378137, + POLAR_RADIUS, + -POLAR_RADIUS, +]; + +function inputs(seed) { + const random = generator(seed); + const scalar = () => { + if (random() < 0.15) return SPECIAL[Math.floor(random() * SPECIAL.length)]; + const scale = [1e-12, 1, 1e3, 6.4e6, 4.2e7][Math.floor(random() * 5)]; + return (random() - 0.5) * 2 * scale; + }; + const vector = () => ({ x: scalar(), y: scalar(), z: scalar() }); + return { random, scalar, vector }; +} + +const label = (value) => (Object.is(value, -0) ? '-0' : String(value)); + +function sameKeys(ours, theirs, keys, what) { + if (theirs === undefined) { + assert.equal(ours, undefined, `${what}: expected undefined`); + return; + } + assert.notEqual(ours, undefined, `${what}: unexpected undefined`); + for (const k of keys) { + assert.ok( + Object.is(ours[k], theirs[k]), + `${what}.${k}: ${label(ours[k])} vs ${label(theirs[k])}`, + ); + } +} + +const XYZ = ['x', 'y', 'z']; +const XYZW = ['x', 'y', 'z', 'w']; +const GEODETIC = ['longitude', 'latitude', 'height']; +const M3 = Array.from({ length: 9 }, (_, i) => i); +const M4 = Array.from({ length: 16 }, (_, i) => i); + +/** Runs both sides; a throw on one side must be a throw on the other. */ +function compare(ours, theirs, keys, what) { + let a; + let b; + let aThrew = false; + let bThrew = false; + try { + a = ours(); + } catch { + aThrew = true; + } + try { + b = theirs(); + } catch { + bThrew = true; + } + assert.equal( + aThrew, + bThrew, + `${what}: throw mismatch (ours ${aThrew}, Cesium ${bThrew})`, + ); + if (aThrew) return; + if (keys === null) + assert.ok(Object.is(a, b), `${what}: ${label(a)} vs ${label(b)}`); + else sameKeys(a, b, keys, what); +} + +const plain = (v) => ({ x: v.x, y: v.y, z: v.z }); +const cesium = (v) => new C3(v.x, v.y, v.z); + +test('vector operations, including -0 inputs and aliased results, equal Cesium bit for bit', () => { + const { scalar, vector } = inputs(1); + // Directionless vectors first: normalize must throw on both sides. + const DEGENERATE = [ + { x: 0, y: 0, z: 0 }, + { x: -0, y: 0, z: -0 }, + { x: Number.NaN, y: 1, z: 2 }, + ]; + for (let i = 0; i < CASES; i++) { + const a = i < DEGENERATE.length ? DEGENERATE[i] : vector(); + const b = vector(); + const t = scalar(); + for (const op of [ + 'add', + 'subtract', + 'cross', + 'multiplyComponents', + 'midpoint', + ]) { + compare( + () => c3[op](a, b, c3.create()), + () => C3[op](cesium(a), cesium(b), new C3()), + XYZ, + op, + ); + } + for (const op of ['dot', 'distance', 'distanceSquared']) { + compare( + () => c3[op](a, b), + () => C3[op](cesium(a), cesium(b)), + null, + op, + ); + } + compare( + () => c3.magnitude(a), + () => C3.magnitude(cesium(a)), + null, + 'magnitude', + ); + compare( + () => c3.negate(a, c3.create()), + () => C3.negate(cesium(a), new C3()), + XYZ, + 'negate', + ); + compare( + () => c3.multiplyByScalar(a, t, c3.create()), + () => C3.multiplyByScalar(cesium(a), t, new C3()), + XYZ, + 'multiplyByScalar', + ); + compare( + () => c3.divideByScalar(a, t, c3.create()), + () => C3.divideByScalar(cesium(a), t, new C3()), + XYZ, + 'divideByScalar', + ); + compare( + () => c3.equalsEpsilon(a, b, Math.abs(t), Math.abs(t) / 3), + () => + C3.equalsEpsilon(cesium(a), cesium(b), Math.abs(t), Math.abs(t) / 3), + null, + 'equalsEpsilon', + ); + compare( + () => c3.lerp(a, b, t, c3.create()), + () => C3.lerp(cesium(a), cesium(b), t, new C3()), + XYZ, + 'lerp', + ); + + // Aliasing: the result is one of the inputs on both sides. + compare( + () => { + const left = plain(a); + return c3.lerp(left, b, t, left); + }, + () => { + const left = cesium(a); + return C3.lerp(left, cesium(b), t, left); + }, + XYZ, + 'lerp result = start', + ); + compare( + () => { + const right = plain(b); + return c3.lerp(a, right, t, right); + }, + () => { + const right = cesium(b); + return C3.lerp(cesium(a), right, t, right); + }, + XYZ, + 'lerp result = end', + ); + compare( + () => { + const left = plain(a); + return c3.cross(left, b, left); + }, + () => { + const left = cesium(a); + return C3.cross(left, cesium(b), left); + }, + XYZ, + 'cross result = left', + ); + compare( + () => { + const right = plain(b); + return c3.cross(a, right, right); + }, + () => { + const right = cesium(b); + return C3.cross(cesium(a), right, right); + }, + XYZ, + 'cross result = right', + ); + compare( + () => { + const v = plain(a); + return c3.normalize(v, v); + }, + () => { + const v = cesium(a); + return C3.normalize(v, v); + }, + XYZ, + 'normalize in place', + ); + } +}); + +test('geodetic conversions, including the centre, poles and -0 inputs, equal Cesium bit for bit', () => { + const { random, scalar, vector } = inputs(2); + const W = Cesium.Ellipsoid.WGS84; + for (let i = 0; i < CASES; i++) { + const p = vector(); + const cp = cesium(p); + compare( + () => ellipsoid.WGS84.geodeticSurfaceNormal(p), + () => W.geodeticSurfaceNormal(cp), + XYZ, + 'normal', + ); + compare( + () => ellipsoid.WGS84.scaleToGeodeticSurface(p), + () => W.scaleToGeodeticSurface(cp), + XYZ, + 'scale', + ); + compare( + () => ellipsoid.WGS84.cartesianToCartographic(p), + () => W.cartesianToCartographic(cp), + GEODETIC, + 'cartesianToCartographic', + ); + compare( + () => cartographic.fromCartesian(p), + () => Cesium.Cartographic.fromCartesian(cp), + GEODETIC, + 'fromCartesian', + ); + compare( + () => { + const v = plain(p); + return ellipsoid.WGS84.geodeticSurfaceNormal(v, v); + }, + () => { + const v = cesium(p); + return W.geodeticSurfaceNormal(v, v); + }, + XYZ, + 'normal in place', + ); + + const g = { + longitude: random() < 0.1 ? -0 : (random() - 0.5) * 8, + latitude: random() < 0.1 ? Math.PI / 2 : (random() - 0.5) * Math.PI, + height: scalar() / 10, + }; + const cg = new Cesium.Cartographic(g.longitude, g.latitude, g.height); + compare( + () => ellipsoid.WGS84.cartographicToCartesian(g), + () => W.cartographicToCartesian(cg), + XYZ, + 'cartographicToCartesian', + ); + compare( + () => cartographic.toCartesian(g), + () => Cesium.Cartographic.toCartesian(cg), + XYZ, + 'toCartesian', + ); + const lon = (random() - 0.5) * 720; + const lat = (random() - 0.5) * 180; + compare( + () => c3.fromDegrees(lon, lat, g.height), + () => C3.fromDegrees(lon, lat, g.height), + XYZ, + 'fromDegrees', + ); + } +}); + +/** A rotation by angle around axis (Cesium's own construction, so both sides get the same matrix). */ +function rotation(axis, angle) { + const q = Cesium.Quaternion.fromAxisAngle(cesium(axis), angle); + return Array.from(Cesium.Matrix3.fromQuaternion(q, new Cesium.Matrix3())); +} + +test('matrices and quaternions, including aliased results and every fromRotationMatrix branch, equal Cesium', () => { + const { random, scalar, vector } = inputs(3); + const branches = new Set(); + for (let i = 0; i < CASES; i++) { + const a = Array.from({ length: 16 }, () => scalar()); + const b = Array.from({ length: 16 }, () => scalar()); + const ca = Cesium.Matrix4.fromArray(a); + const cb = Cesium.Matrix4.fromArray(b); + for (const op of ['multiply', 'multiplyTransformation']) { + compare( + () => matrix4[op](a, b, matrix4.create()), + () => Cesium.Matrix4[op](ca, cb, new Cesium.Matrix4()), + M4, + op, + ); + compare( + () => { + const left = [...a]; + return matrix4[op](left, b, left); + }, + () => { + const left = ca.clone(); + return Cesium.Matrix4[op](left, cb, left); + }, + M4, + `${op} result = left`, + ); + compare( + () => { + const right = [...b]; + return matrix4[op](a, right, right); + }, + () => { + const right = cb.clone(); + return Cesium.Matrix4[op](ca, right, right); + }, + M4, + `${op} result = right`, + ); + } + compare( + () => { + const m = [...a]; + return matrix4.inverseTransformation(m, m); + }, + () => { + const m = ca.clone(); + return Cesium.Matrix4.inverseTransformation(m, m); + }, + M4, + 'inverseTransformation in place', + ); + const p = vector(); + for (const op of ['multiplyByPoint', 'multiplyByPointAsVector']) { + compare( + () => { + const v = plain(p); + return matrix4[op](a, v, v); + }, + () => { + const v = cesium(p); + return Cesium.Matrix4[op](ca, v, v); + }, + XYZ, + `${op} in place`, + ); + } + + // Rotations near a half turn have trace <= 0, which selects the i = 0, 1, 2 branches. + const axis = [ + { x: 1, y: 0.1 * random(), z: 0.1 * random() }, + { x: 0.1 * random(), y: 1, z: 0.1 * random() }, + { x: 0.1 * random(), y: 0.1 * random(), z: 1 }, + vector(), + ][i % 4]; + const angle = + i % 8 < 4 ? Math.PI - 0.5 * random() : (random() - 0.5) * 2 * Math.PI; + let m; + try { + m = rotation(axis, angle); + } catch { + continue; // a zero axis: Cesium throws, nothing to compare + } + const trace = m[0] + m[4] + m[8]; + branches.add( + trace > 0 + ? 'trace>0' + : m[8] > m[0] && m[8] > m[4] + ? 'i=2' + : m[4] > m[0] + ? 'i=1' + : 'i=0', + ); + compare( + () => quaternion.fromRotationMatrix(m), + () => Cesium.Quaternion.fromRotationMatrix(Cesium.Matrix3.fromArray(m)), + XYZW, + `fromRotationMatrix (trace ${trace})`, + ); + const q = { x: scalar(), y: scalar(), z: scalar(), w: scalar() }; + compare( + () => matrix3.fromQuaternion(q), + () => + Cesium.Matrix3.fromQuaternion( + new Cesium.Quaternion(q.x, q.y, q.z, q.w), + ), + M3, + 'fromQuaternion', + ); + compare( + () => quaternion.fromAxisAngle(axis, angle), + () => Cesium.Quaternion.fromAxisAngle(cesium(axis), angle), + XYZW, + 'fromAxisAngle', + ); + } + assert.deepEqual([...branches].sort(), ['i=0', 'i=1', 'i=2', 'trace>0']); +}); + +test('local frames, including the centre and the exact poles, equal Cesium bit for bit', () => { + const { random, vector } = inputs(4); + const DEGENERATE = [ + { x: 0, y: 0, z: 0 }, + { x: -0, y: -0, z: -0 }, + { x: 0, y: 0, z: POLAR_RADIUS }, + { x: 0, y: 0, z: -POLAR_RADIUS }, + { x: -0, y: 1e-15, z: 7 }, + ]; + for (let i = 0; i < CASES; i++) { + const origin = i < DEGENERATE.length ? DEGENERATE[i] : vector(); + compare( + () => transforms.eastNorthUpToFixedFrame(origin), + () => Cesium.Transforms.eastNorthUpToFixedFrame(cesium(origin)), + M4, + 'eastNorthUpToFixedFrame', + ); + const hpr = { + heading: (random() - 0.5) * 10, + pitch: (random() - 0.5) * 4, + roll: (random() - 0.5) * 4, + }; + compare( + () => + transforms.headingPitchRollToFixedFrame( + origin, + hpr, + undefined, + undefined, + matrix4.create(), + ), + () => + Cesium.Transforms.headingPitchRollToFixedFrame( + cesium(origin), + new Cesium.HeadingPitchRoll(hpr.heading, hpr.pitch, hpr.roll), + undefined, + undefined, + new Cesium.Matrix4(), + ), + M4, + 'headingPitchRollToFixedFrame', + ); + } +}); diff --git a/src/geo/cartesian2.js b/src/geo/cartesian2.js new file mode 100644 index 000000000..44cefb563 --- /dev/null +++ b/src/geo/cartesian2.js @@ -0,0 +1,49 @@ +// Ported from CesiumJS 1.138 (Apache-2.0) — Cartesian2.js +// +// 2D vectors (screen coordinates) as plain-data {x, y} objects. With DEBUG on +// (development and tests; see debug.js) invalid arguments throw the same +// DeveloperError as Cesium's debug checks. +// +// Geo values are not Object.prototype objects, so src/app/stateChannel.js +// copyState and src/director/documentFields.js jsonTree reject them: copy +// the fields into a literal before putting one into product state or a +// scene document. + +import { DEBUG } from './debug.js'; +import { Check } from './check.js'; + +// Values are instances of a module-private constructor rather than object +// literals: every `{x, y}` literal in the process shares one V8 hidden +// class, and code elsewhere that builds such literals would make geo's hot +// functions polymorphic. The instances still look like plain data: own +// enumerable fields only, no prototype methods or getters. +function Vector2(x, y) { + this.x = x; + this.y = y; +} + +/** A vector (Cesium's `new Cartesian2(x, y)`). */ +export function create(x, y) { + return new Vector2(x ?? 0.0, y ?? 0.0); +} + +/** Copies a vector; `undefined` when there is none. */ +export function clone(cartesian, result) { + if (cartesian == null) return undefined; + if (result == null) return create(cartesian.x, cartesian.y); + result.x = cartesian.x; + result.y = cartesian.y; + return result; +} + +/** The distance between two points. */ +export function distance(left, right) { + if (DEBUG) { + Check.typeOf.object('left', left); + Check.typeOf.object('right', right); + } + + const x = left.x - right.x; + const y = left.y - right.y; + return Math.sqrt(x * x + y * y); +} diff --git a/src/geo/cartesian2.test.mjs b/src/geo/cartesian2.test.mjs new file mode 100644 index 000000000..1bf0cb188 --- /dev/null +++ b/src/geo/cartesian2.test.mjs @@ -0,0 +1,19 @@ +// src/geo/cartesian2.test.mjs +import { test } from 'node:test'; +import assert from 'node:assert/strict'; +import * as Cesium from 'cesium'; +import * as c2 from './cartesian2.js'; + +test('cartesian2 create, clone and distance equal Cesium', () => { + const a = c2.create(3, 4); + assert.deepEqual({ ...a }, { x: 3, y: 4 }); + const out = c2.create(); + assert.equal(c2.clone(a, out), out); + assert.equal( + c2.distance(a, c2.create()), + Cesium.Cartesian2.distance( + new Cesium.Cartesian2(3, 4), + new Cesium.Cartesian2(), + ), + ); +}); diff --git a/src/geo/cartesian3.js b/src/geo/cartesian3.js new file mode 100644 index 000000000..18f8ef3ad --- /dev/null +++ b/src/geo/cartesian3.js @@ -0,0 +1,389 @@ +// Ported from CesiumJS 1.138 (Apache-2.0) — Cartesian3.js +// +// 3D vectors as plain-data {x, y, z} objects. Arithmetic takes a required `result` +// (as Cesium does) that may be one of the inputs; creators take an optional +// one. The operation order matches Cesium so the outputs are bit-identical. +// With DEBUG on (development and tests; see debug.js) invalid arguments throw +// the same DeveloperError as Cesium's debug checks, and `normalize` throws on a +// zero or non-finite vector; production builds skip the checks and return NaN +// components there, as release Cesium does. +// +// Geo values are not Object.prototype objects, so src/app/stateChannel.js +// copyState and src/director/documentFields.js jsonTree reject them: copy +// the fields into a literal before putting one into product state or a +// scene document. + +import { DEBUG } from './debug.js'; +import { Check } from './check.js'; +import { DeveloperError } from './developerError.js'; +import { equalsEpsilon as scalarEqualsEpsilon, toRadians } from './math.js'; + +// Values are instances of a module-private constructor rather than object +// literals: every `{x, y, z}` literal in the process shares one V8 hidden +// class, and code elsewhere that builds such literals would make geo's hot +// functions polymorphic. The instances still look like plain data: own +// enumerable fields only, no prototype methods or getters. +function Vector3(x, y, z) { + this.x = x; + this.y = y; + this.z = z; +} + +/** A vector (Cesium's `new Cartesian3(x, y, z)`). */ +export function create(x, y, z) { + return new Vector3(x ?? 0.0, y ?? 0.0, z ?? 0.0); +} + +/** Sets a vector's components. */ +export function fromElements(x, y, z, result) { + if (result == null) return create(x, y, z); + result.x = x; + result.y = y; + result.z = z; + return result; +} + +/** Copies a vector; `undefined` when there is none. */ +export function clone(cartesian, result) { + if (cartesian == null) return undefined; + if (result == null) return create(cartesian.x, cartesian.y, cartesian.z); + result.x = cartesian.x; + result.y = cartesian.y; + result.z = cartesian.z; + return result; +} + +export const ZERO = Object.freeze(create(0.0, 0.0, 0.0)); +export const UNIT_X = Object.freeze(create(1.0, 0.0, 0.0)); +export const UNIT_Y = Object.freeze(create(0.0, 1.0, 0.0)); +export const UNIT_Z = Object.freeze(create(0.0, 0.0, 1.0)); + +export function magnitudeSquared(cartesian) { + if (DEBUG) { + Check.typeOf.object('cartesian', cartesian); + } + + return ( + cartesian.x * cartesian.x + + cartesian.y * cartesian.y + + cartesian.z * cartesian.z + ); +} + +export function magnitude(cartesian) { + return Math.sqrt(magnitudeSquared(cartesian)); +} + +const distanceScratch = create(); + +export function distance(left, right) { + if (DEBUG) { + Check.typeOf.object('left', left); + Check.typeOf.object('right', right); + } + + subtract(left, right, distanceScratch); + return magnitude(distanceScratch); +} + +export function distanceSquared(left, right) { + if (DEBUG) { + Check.typeOf.object('left', left); + Check.typeOf.object('right', right); + } + + subtract(left, right, distanceScratch); + return magnitudeSquared(distanceScratch); +} + +/** The unit vector; with DEBUG on, throws when the vector has no direction (zero or non-finite). */ +export function normalize(cartesian, result) { + if (DEBUG) { + Check.typeOf.object('cartesian', cartesian); + Check.typeOf.object('result', result); + } + + const length = magnitude(cartesian); + result.x = cartesian.x / length; + result.y = cartesian.y / length; + result.z = cartesian.z / length; + if (DEBUG) { + if (isNaN(result.x) || isNaN(result.y) || isNaN(result.z)) { + throw new DeveloperError('normalized result is not a number'); + } + } + return result; +} + +export function dot(left, right) { + if (DEBUG) { + Check.typeOf.object('left', left); + Check.typeOf.object('right', right); + } + + return left.x * right.x + left.y * right.y + left.z * right.z; +} + +export function multiplyComponents(left, right, result) { + if (DEBUG) { + Check.typeOf.object('left', left); + Check.typeOf.object('right', right); + Check.typeOf.object('result', result); + } + + result.x = left.x * right.x; + result.y = left.y * right.y; + result.z = left.z * right.z; + return result; +} + +export function add(left, right, result) { + if (DEBUG) { + Check.typeOf.object('left', left); + Check.typeOf.object('right', right); + Check.typeOf.object('result', result); + } + + result.x = left.x + right.x; + result.y = left.y + right.y; + result.z = left.z + right.z; + return result; +} + +export function subtract(left, right, result) { + if (DEBUG) { + Check.typeOf.object('left', left); + Check.typeOf.object('right', right); + Check.typeOf.object('result', result); + } + + result.x = left.x - right.x; + result.y = left.y - right.y; + result.z = left.z - right.z; + return result; +} + +export function multiplyByScalar(cartesian, scalar, result) { + if (DEBUG) { + Check.typeOf.object('cartesian', cartesian); + Check.typeOf.number('scalar', scalar); + Check.typeOf.object('result', result); + } + + result.x = cartesian.x * scalar; + result.y = cartesian.y * scalar; + result.z = cartesian.z * scalar; + return result; +} + +export function divideByScalar(cartesian, scalar, result) { + if (DEBUG) { + Check.typeOf.object('cartesian', cartesian); + Check.typeOf.number('scalar', scalar); + Check.typeOf.object('result', result); + } + + result.x = cartesian.x / scalar; + result.y = cartesian.y / scalar; + result.z = cartesian.z / scalar; + return result; +} + +export function negate(cartesian, result) { + if (DEBUG) { + Check.typeOf.object('cartesian', cartesian); + Check.typeOf.object('result', result); + } + + result.x = -cartesian.x; + result.y = -cartesian.y; + result.z = -cartesian.z; + return result; +} + +const lerpScratch = create(); + +export function lerp(start, end, t, result) { + if (DEBUG) { + Check.typeOf.object('start', start); + Check.typeOf.object('end', end); + Check.typeOf.number('t', t); + Check.typeOf.object('result', result); + } + + multiplyByScalar(end, t, lerpScratch); + result = multiplyByScalar(start, 1.0 - t, result); + return add(lerpScratch, result, result); +} + +export function cross(left, right, result) { + if (DEBUG) { + Check.typeOf.object('left', left); + Check.typeOf.object('right', right); + Check.typeOf.object('result', result); + } + + const leftX = left.x; + const leftY = left.y; + const leftZ = left.z; + const rightX = right.x; + const rightY = right.y; + const rightZ = right.z; + + const x = leftY * rightZ - leftZ * rightY; + const y = leftZ * rightX - leftX * rightZ; + const z = leftX * rightY - leftY * rightX; + + result.x = x; + result.y = y; + result.z = z; + return result; +} + +export function midpoint(left, right, result) { + if (DEBUG) { + Check.typeOf.object('left', left); + Check.typeOf.object('right', right); + Check.typeOf.object('result', result); + } + + result.x = (left.x + right.x) * 0.5; + result.y = (left.y + right.y) * 0.5; + result.z = (left.z + right.z) * 0.5; + return result; +} + +export function equals(left, right) { + return ( + left === right || + (left != null && + right != null && + left.x === right.x && + left.y === right.y && + left.z === right.z) + ); +} + +export function equalsEpsilon(left, right, relativeEpsilon, absoluteEpsilon) { + return ( + left === right || + (left != null && + right != null && + scalarEqualsEpsilon(left.x, right.x, relativeEpsilon, absoluteEpsilon) && + scalarEqualsEpsilon(left.y, right.y, relativeEpsilon, absoluteEpsilon) && + scalarEqualsEpsilon(left.z, right.z, relativeEpsilon, absoluteEpsilon)) + ); +} + +// WGS84 radii squared, used when no ellipsoid is given (Cesium's +// Cartesian3._ellipsoidRadiiSquared; equal to ellipsoid.WGS84.radiiSquared). +// Cesium follows Ellipsoid.default here; the default is WGS84 because gev never +// changes Ellipsoid.default. +const WGS84_RADII_SQUARED = Object.freeze( + create( + 6378137.0 * 6378137.0, + 6378137.0 * 6378137.0, + 6356752.3142451793 * 6356752.3142451793, + ), +); + +let scratchN = create(); +let scratchK = create(); + +/** The Earth-fixed position of a geodetic one in radians on the ellipsoid (WGS84 by default). */ +export function fromRadians(longitude, latitude, height, ellipsoid, result) { + if (DEBUG) { + Check.typeOf.number('longitude', longitude); + Check.typeOf.number('latitude', latitude); + } + + height = height ?? 0.0; + const radiiSquared = + ellipsoid == null ? WGS84_RADII_SQUARED : ellipsoid.radiiSquared; + + const cosLatitude = Math.cos(latitude); + scratchN.x = cosLatitude * Math.cos(longitude); + scratchN.y = cosLatitude * Math.sin(longitude); + scratchN.z = Math.sin(latitude); + scratchN = normalize(scratchN, scratchN); + + multiplyComponents(radiiSquared, scratchN, scratchK); + const gamma = Math.sqrt(dot(scratchN, scratchK)); + scratchK = divideByScalar(scratchK, gamma, scratchK); + scratchN = multiplyByScalar(scratchN, height, scratchN); + + if (result == null) result = create(); + return add(scratchK, scratchN, result); +} + +/** The Earth-fixed position of a geodetic one in degrees on the ellipsoid (WGS84 by default). */ +export function fromDegrees(longitude, latitude, height, ellipsoid, result) { + if (DEBUG) { + Check.typeOf.number('longitude', longitude); + Check.typeOf.number('latitude', latitude); + } + + longitude = toRadians(longitude); + latitude = toRadians(latitude); + return fromRadians(longitude, latitude, height, ellipsoid, result); +} + +/** Positions from a flat [lon, lat, lon, lat, …] list in degrees, at height 0. */ +export function fromDegreesArray(coordinates, ellipsoid, result) { + if (DEBUG) { + Check.defined('coordinates', coordinates); + if (coordinates.length < 2 || coordinates.length % 2 !== 0) { + throw new DeveloperError( + 'the number of coordinates must be a multiple of 2 and at least 2', + ); + } + } + + const length = coordinates.length; + if (result == null) { + result = new Array(length / 2); + } else { + result.length = length / 2; + } + for (let i = 0; i < length; i += 2) { + const index = i / 2; + result[index] = fromDegrees( + coordinates[i], + coordinates[i + 1], + 0, + ellipsoid, + result[index], + ); + } + return result; +} + +/** Positions from a flat [lon, lat, height, …] list in degrees and metres. */ +export function fromDegreesArrayHeights(coordinates, ellipsoid, result) { + if (DEBUG) { + Check.defined('coordinates', coordinates); + if (coordinates.length < 3 || coordinates.length % 3 !== 0) { + throw new DeveloperError( + 'the number of coordinates must be a multiple of 3 and at least 3', + ); + } + } + + const length = coordinates.length; + if (result == null) { + result = new Array(length / 3); + } else { + result.length = length / 3; + } + for (let i = 0; i < length; i += 3) { + const index = i / 3; + result[index] = fromDegrees( + coordinates[i], + coordinates[i + 1], + coordinates[i + 2], + ellipsoid, + result[index], + ); + } + return result; +} diff --git a/src/geo/cartesian3.test.mjs b/src/geo/cartesian3.test.mjs new file mode 100644 index 000000000..ba661f550 --- /dev/null +++ b/src/geo/cartesian3.test.mjs @@ -0,0 +1,157 @@ +// src/geo/cartesian3.test.mjs +import { test } from 'node:test'; +import assert from 'node:assert/strict'; +import * as Cesium from 'cesium'; +import * as c3 from './cartesian3.js'; + +const C = Cesium.Cartesian3; +const SAMPLES = [ + [1, 2, 3], + [-4.5, 0.25, 1e6], + [0, 0, 0], + [6378137, -1e-3, 42], + [1e-12, -1e-12, 3], +]; +const asCesium = ([x, y, z]) => new C(x, y, z); +const asPlain = ([x, y, z]) => ({ x, y, z }); +const same = (ours, theirs, label) => + assert.deepEqual( + { x: ours.x, y: ours.y, z: ours.z }, + { x: theirs.x, y: theirs.y, z: theirs.z }, + label, + ); + +test('binary vector ops equal Cesium bit for bit', () => { + for (const a of SAMPLES) + for (const b of SAMPLES) { + for (const op of ['add', 'subtract', 'cross', 'multiplyComponents']) { + same( + c3[op](asPlain(a), asPlain(b), c3.create()), + C[op](asCesium(a), asCesium(b), new C()), + `${op}`, + ); + } + assert.equal( + c3.dot(asPlain(a), asPlain(b)), + C.dot(asCesium(a), asCesium(b)), + ); + assert.equal( + c3.distance(asPlain(a), asPlain(b)), + C.distance(asCesium(a), asCesium(b)), + ); + assert.equal( + c3.distanceSquared(asPlain(a), asPlain(b)), + C.distanceSquared(asCesium(a), asCesium(b)), + ); + same( + c3.lerp(asPlain(a), asPlain(b), 0.3, c3.create()), + C.lerp(asCesium(a), asCesium(b), 0.3, new C()), + 'lerp', + ); + same( + c3.midpoint(asPlain(a), asPlain(b), c3.create()), + C.midpoint(asCesium(a), asCesium(b), new C()), + 'midpoint', + ); + assert.equal( + c3.equalsEpsilon(asPlain(a), asPlain(b), 1e-7), + C.equalsEpsilon(asCesium(a), asCesium(b), 1e-7), + ); + } +}); + +test('unary ops and scalars equal Cesium', () => { + for (const a of SAMPLES) { + assert.equal(c3.magnitude(asPlain(a)), C.magnitude(asCesium(a))); + assert.equal( + c3.magnitudeSquared(asPlain(a)), + C.magnitudeSquared(asCesium(a)), + ); + same( + c3.negate(asPlain(a), c3.create()), + C.negate(asCesium(a), new C()), + 'negate', + ); + same( + c3.multiplyByScalar(asPlain(a), 2.5, c3.create()), + C.multiplyByScalar(asCesium(a), 2.5, new C()), + 'mul', + ); + same( + c3.divideByScalar(asPlain(a), 4, c3.create()), + C.divideByScalar(asCesium(a), 4, new C()), + 'div', + ); + if (C.magnitude(asCesium(a)) !== 0) { + same( + c3.normalize(asPlain(a), c3.create()), + C.normalize(asCesium(a), new C()), + 'normalize', + ); + } + } +}); + +test('geodetic constructors equal Cesium', () => { + for (const [lon, lat, h] of [ + [0, 0, 0], + [-73.9857, 40.7484, 381], + [139.76, 35.68, -20], + [180, -90, 1e4], + [12.5, 89.999999, 1], + ]) { + same( + c3.fromDegrees(lon, lat, h), + C.fromDegrees(lon, lat, h), + `fromDegrees ${lon},${lat}`, + ); + const lr = [Cesium.Math.toRadians(lon), Cesium.Math.toRadians(lat), h]; + same(c3.fromRadians(...lr), C.fromRadians(...lr), 'fromRadians'); + } + const flat = [-74, 40.7, -73.9, 40.8, -73.8, 40.9]; + c3.fromDegreesArray(flat).forEach((p, i) => + same(p, C.fromDegreesArray(flat)[i], `array ${i}`), + ); + const heights = [-74, 40.7, 10, -73.9, 40.8, 20]; + c3.fromDegreesArrayHeights(heights).forEach((p, i) => + same(p, C.fromDegreesArrayHeights(heights)[i], `heights ${i}`), + ); +}); + +test('result aliasing works like Cesium and creators reuse result', () => { + const a = { x: 1, y: 2, z: 3 }; + c3.add(a, a, a); + assert.deepEqual({ ...a }, { x: 2, y: 4, z: 6 }); + const out = c3.create(); + assert.equal(c3.fromDegrees(1, 2, 3, undefined, out), out); + assert.ok(Object.isFrozen(c3.ZERO)); +}); + +test('plain results are accepted by Cesium functions (duck typing)', () => { + const p = c3.fromDegrees(-73.9857, 40.7484, 381); + const carto = Cesium.Cartographic.fromCartesian(p); + assert.ok(Math.abs(Cesium.Math.toDegrees(carto.latitude) - 40.7484) < 1e-9); + assert.ok( + Cesium.Cartesian3.equals( + Cesium.Cartesian3.clone(p), + C.fromDegrees(-73.9857, 40.7484, 381), + ), + ); +}); + +test('normalize throws on a zero or NaN vector, as Cesium does when run from source', () => { + for (const [x, y, z] of [ + [0, 0, 0], + [-0, 0, -0], + [Number.NaN, 1, 2], + ]) { + assert.throws( + () => C.normalize(new C(x, y, z), new C()), + /normalized result is not a number/, + ); + assert.throws( + () => c3.normalize({ x, y, z }, c3.create()), + /normalized result is not a number/, + ); + } +}); diff --git a/src/geo/cartographic.js b/src/geo/cartographic.js new file mode 100644 index 000000000..5500a1f94 --- /dev/null +++ b/src/geo/cartographic.js @@ -0,0 +1,163 @@ +// Ported from CesiumJS 1.138 (Apache-2.0) — Cartographic.js +// +// Geodetic positions as plain-data {longitude, latitude, height} objects (radians and +// metres). Every function takes an optional `result` and returns it. With +// DEBUG on (development and tests; see debug.js) invalid arguments throw the +// same DeveloperError as Cesium's debug checks. +// +// Geo values are not Object.prototype objects, so src/app/stateChannel.js +// copyState and src/director/documentFields.js jsonTree reject them: copy +// the fields into a literal before putting one into product state or a +// scene document. + +import { DEBUG } from './debug.js'; +import { Check } from './check.js'; +import { toRadians } from './math.js'; +import * as cartesian3 from './cartesian3.js'; +import { scaleToGeodeticSurface } from './scaleToGeodeticSurface.js'; + +// Values are instances of a module-private constructor rather than object +// literals: every `{longitude, latitude, height}` literal in the process shares one V8 hidden +// class, and code elsewhere that builds such literals would make geo's hot +// functions polymorphic. The instances still look like plain data: own +// enumerable fields only, no prototype methods or getters. +function Geodetic(longitude, latitude, height) { + this.longitude = longitude; + this.latitude = latitude; + this.height = height; +} + +/** A geodetic position (Cesium's `new Cartographic(longitude, latitude, height)`). */ +export function create(longitude, latitude, height) { + return new Geodetic(longitude ?? 0.0, latitude ?? 0.0, height ?? 0.0); +} + +// WGS84 constants used when no ellipsoid is given (Cesium's +// Cartographic._ellipsoidOneOverRadii and friends; bit-equal to +// ellipsoid.WGS84's fields). Local, as in Cesium, so that ellipsoid.js can +// import this module without a cycle. +const WGS84_ONE_OVER_RADII = Object.freeze( + cartesian3.create(1.0 / 6378137.0, 1.0 / 6378137.0, 1.0 / 6356752.3142451793), +); +const WGS84_ONE_OVER_RADII_SQUARED = Object.freeze( + cartesian3.create( + 1.0 / (6378137.0 * 6378137.0), + 1.0 / (6378137.0 * 6378137.0), + 1.0 / (6356752.3142451793 * 6356752.3142451793), + ), +); +// Cesium's CesiumMath.EPSILON1. +const WGS84_CENTER_TOLERANCE_SQUARED = 0.1; + +/** A geodetic position from radians. */ +export function fromRadians(longitude, latitude, height, result) { + if (DEBUG) { + Check.typeOf.number('longitude', longitude); + Check.typeOf.number('latitude', latitude); + } + + height = height ?? 0.0; + if (result == null) return create(longitude, latitude, height); + result.longitude = longitude; + result.latitude = latitude; + result.height = height; + return result; +} + +/** A geodetic position from degrees; the result is in radians. */ +export function fromDegrees(longitude, latitude, height, result) { + if (DEBUG) { + Check.typeOf.number('longitude', longitude); + Check.typeOf.number('latitude', latitude); + } + + longitude = toRadians(longitude); + latitude = toRadians(latitude); + return fromRadians(longitude, latitude, height, result); +} + +const cartesianToCartographicN = cartesian3.create(); +const cartesianToCartographicP = cartesian3.create(); +const cartesianToCartographicH = cartesian3.create(); + +/** + * The geodetic position of an Earth-fixed point, or `undefined` when the point + * is at the centre of the ellipsoid (WGS84 by default). + */ +export function fromCartesian(cartesian, ellipsoid, result) { + // Cesium falls back to Ellipsoid.default; gev never changes it, so WGS84. + // A Cesium Ellipsoid keeps the tolerance in `_centerToleranceSquared`. + const oneOverRadii = + ellipsoid == null ? WGS84_ONE_OVER_RADII : ellipsoid.oneOverRadii; + const oneOverRadiiSquared = + ellipsoid == null + ? WGS84_ONE_OVER_RADII_SQUARED + : ellipsoid.oneOverRadiiSquared; + const centerToleranceSquared = + ellipsoid == null + ? WGS84_CENTER_TOLERANCE_SQUARED + : (ellipsoid.centerToleranceSquared ?? ellipsoid._centerToleranceSquared); + + const p = scaleToGeodeticSurface( + cartesian, + oneOverRadii, + oneOverRadiiSquared, + centerToleranceSquared, + cartesianToCartographicP, + ); + if (p === undefined) return undefined; + + let n = cartesian3.multiplyComponents( + p, + oneOverRadiiSquared, + cartesianToCartographicN, + ); + n = cartesian3.normalize(n, n); + + const h = cartesian3.subtract(cartesian, p, cartesianToCartographicH); + + const longitude = Math.atan2(n.y, n.x); + const latitude = Math.asin(n.z); + const height = + Math.sign(cartesian3.dot(h, cartesian)) * cartesian3.magnitude(h); + + if (result == null) return create(longitude, latitude, height); + result.longitude = longitude; + result.latitude = latitude; + result.height = height; + return result; +} + +/** + * The Earth-fixed position of a geodetic one on the ellipsoid (WGS84 by default, + * where Cesium uses Ellipsoid.default; gev never changes it). + */ +export function toCartesian(cartographic, ellipsoid, result) { + if (DEBUG) { + Check.defined('cartographic', cartographic); + } + + return cartesian3.fromRadians( + cartographic.longitude, + cartographic.latitude, + cartographic.height, + ellipsoid, + result, + ); +} + +/** Copies a geodetic position; `undefined` when there is none. */ +export function clone(cartographic, result) { + if (cartographic == null) return undefined; + if (result == null) { + return create( + cartographic.longitude, + cartographic.latitude, + cartographic.height, + ); + } + result.longitude = cartographic.longitude; + result.latitude = cartographic.latitude; + result.height = cartographic.height; + return result; +} diff --git a/src/geo/cartographic.test.mjs b/src/geo/cartographic.test.mjs new file mode 100644 index 000000000..724312821 --- /dev/null +++ b/src/geo/cartographic.test.mjs @@ -0,0 +1,142 @@ +// src/geo/cartographic.test.mjs +import { test } from 'node:test'; +import assert from 'node:assert/strict'; +import * as Cesium from 'cesium'; +import * as cartographic from './cartographic.js'; + +const close = (a, b, tol, label) => + assert.ok( + Math.abs(a - b) <= tol * Math.max(1, Math.abs(b)), + `${label}: ${a} vs ${b}`, + ); + +test('fromDegrees and fromRadians equal Cesium', () => { + for (const [lon, lat, h] of [ + [0, 0, 0], + [-73.9857, 40.7484, 381], + [139.76, 35.68, -20], + [180, -90, 1e4], + ]) { + const ours = cartographic.fromDegrees(lon, lat, h); + const theirs = Cesium.Cartographic.fromDegrees(lon, lat, h); + assert.deepEqual( + { ...ours }, + { + longitude: theirs.longitude, + latitude: theirs.latitude, + height: theirs.height, + }, + ); + } +}); + +test('fromCartesian and toCartesian equal Cesium, including the centre case', () => { + for (const [x, y, z] of [ + [1334762.2, -4654869.1, 4137287.9], + [6378137, 0, 0], + [0.01, 0, 0], + ]) { + const ours = cartographic.fromCartesian({ x, y, z }); + const theirs = Cesium.Cartographic.fromCartesian( + new Cesium.Cartesian3(x, y, z), + ); + if (theirs === undefined) { + assert.equal(ours, undefined); + continue; + } + close(ours.longitude, theirs.longitude, 1e-12, 'lon'); + close(ours.latitude, theirs.latitude, 1e-12, 'lat'); + close(ours.height, theirs.height, 1e-9, 'h'); + const back = cartographic.toCartesian(ours); + const cBack = Cesium.Cartographic.toCartesian(theirs); + for (const k of ['x', 'y', 'z']) close(back[k], cBack[k], 1e-12, k); + } +}); + +test('result objects are reused', () => { + const out = cartographic.create(); + assert.equal(cartographic.fromDegrees(1, 2, 3, out), out); + assert.equal(cartographic.clone(out).longitude, out.longitude); +}); + +const POLAR_RADIUS = 6356752.3142451793; +const sameFields = (ours, theirs, keys, label) => { + for (const k of keys) { + assert.ok( + Object.is(ours[k], theirs[k]), + `${label}.${k}: ${ours[k]} vs ${theirs[k]}`, + ); + } +}; + +test('fromCartesian is undefined at the exact centre, as in Cesium', () => { + for (const [x, y, z] of [ + [0, 0, 0], + [-0, -0, -0], + ]) { + assert.equal( + Cesium.Cartographic.fromCartesian(new Cesium.Cartesian3(x, y, z)), + undefined, + ); + assert.equal(cartographic.fromCartesian({ x, y, z }), undefined); + } +}); + +test('fromCartesian and toCartesian equal Cesium bit for bit at the poles and the antimeridian', () => { + for (const [x, y, z, label] of [ + [0, 0, POLAR_RADIUS, 'north pole'], + [0, 0, -POLAR_RADIUS, 'south pole'], + [-6378137, 0, 0, 'antimeridian, y = +0'], + [-6378137, -0, 0, 'antimeridian, y = -0'], + [-6378137, -0, -0, 'antimeridian, y = z = -0'], + ]) { + const ours = cartographic.fromCartesian({ x, y, z }); + const theirs = Cesium.Cartographic.fromCartesian( + new Cesium.Cartesian3(x, y, z), + ); + sameFields(ours, theirs, ['longitude', 'latitude', 'height'], label); + sameFields( + cartographic.toCartesian(ours), + Cesium.Cartographic.toCartesian(theirs), + ['x', 'y', 'z'], + `${label} back`, + ); + } + // atan2 keeps the sign of a zero y: +π on one side of the antimeridian, -π on the other. + assert.equal( + cartographic.fromCartesian({ x: -6378137, y: 0, z: 0 }).longitude, + Math.PI, + ); + assert.equal( + cartographic.fromCartesian({ x: -6378137, y: -0, z: 0 }).longitude, + -Math.PI, + ); +}); + +// layers/streetLevel/view.js passes a Cesium Ellipsoid (scene.globe.ellipsoid), +// which has no public `centerToleranceSquared`; geo then reads Cesium's private +// `_centerToleranceSquared`. This pins that coupling to a Cesium private field: +// if a Cesium upgrade renames it, the centre case below stops matching. +test('fromCartesian with a Cesium Ellipsoid equals Cesium, including the centre', () => { + const ellipsoid = Cesium.Ellipsoid.WGS84; + assert.equal(ellipsoid.centerToleranceSquared, undefined); + assert.equal(typeof ellipsoid._centerToleranceSquared, 'number'); + for (const [x, y, z, label] of [ + [0, 0, 0, 'centre'], + [0.1, 0.1, 0.1, 'inside the centre tolerance'], + [6378137, 0, 0, 'equator'], + [1334000, -4654000, 4138000, 'mid-latitude'], + [0, 0, POLAR_RADIUS + 500, 'above the north pole'], + ]) { + const ours = cartographic.fromCartesian({ x, y, z }, ellipsoid); + const theirs = Cesium.Cartographic.fromCartesian( + new Cesium.Cartesian3(x, y, z), + ellipsoid, + ); + if (theirs === undefined) { + assert.equal(ours, undefined, label); + continue; + } + sameFields(ours, theirs, ['longitude', 'latitude', 'height'], label); + } +}); diff --git a/src/geo/check.js b/src/geo/check.js new file mode 100644 index 000000000..16d83f368 --- /dev/null +++ b/src/geo/check.js @@ -0,0 +1,61 @@ +// Ported from CesiumJS 1.138 (Apache-2.0) — Check.js +// +// Cesium's argument checks with the same conditions and message text. Only the +// checks src/geo needs are ported. Callers run them inside `if (DEBUG)` blocks +// (see debug.js), as Cesium runs them inside its debug pragmas, so production +// builds drop the calls and this module with them. + +import { DeveloperError } from './developerError.js'; + +const undefinedMessage = (name) => + `${name} is required, actual value was undefined`; +const typeMessage = (actual, expected, name) => + `Expected ${name} to be typeof ${expected}, actual typeof was ${actual}`; + +/** Throws unless `test` is neither undefined nor null. */ +function defined(name, test) { + if (test === undefined || test === null) { + throw new DeveloperError(undefinedMessage(name)); + } +} + +/** Throws unless `typeof test` is 'number' (NaN passes, as in Cesium). */ +function number(name, test) { + if (typeof test !== 'number') { + throw new DeveloperError(typeMessage(typeof test, 'number', name)); + } +} + +// Attached in a function marked pure, not by top-level assignments, so that a +// production build, where no caller runs a check, drops this module entirely. +function withBounds(check) { + check.lessThanOrEquals = function (name, test, limit) { + check(name, test); + if (test > limit) { + throw new DeveloperError( + `Expected ${name} to be less than or equal to ${limit}, actual value was ${test}`, + ); + } + }; + check.greaterThanOrEquals = function (name, test, limit) { + check(name, test); + if (test < limit) { + throw new DeveloperError( + `Expected ${name} to be greater than or equal to ${limit}, actual value was ${test}`, + ); + } + }; + return check; +} + +/** Throws unless `typeof test` is 'object' (null passes, as in Cesium). */ +function object(name, test) { + if (typeof test !== 'object') { + throw new DeveloperError(typeMessage(typeof test, 'object', name)); + } +} + +export const Check = /* @__PURE__ */ Object.freeze({ + defined, + typeOf: Object.freeze({ number: /* @__PURE__ */ withBounds(number), object }), +}); diff --git a/src/geo/checks.test.mjs b/src/geo/checks.test.mjs new file mode 100644 index 000000000..41f59b0d0 --- /dev/null +++ b/src/geo/checks.test.mjs @@ -0,0 +1,486 @@ +// Argument-check parity with Cesium. gev runs Cesium from source, where its +// debug checks are live, so for any input both sides must throw an error with +// the same name and message, or neither may throw. The inputs are seeded and +// deterministic: each case builds its arguments twice from the same seed, once +// for geo and once for Cesium, so no object is shared between the two sides. +import { test } from 'node:test'; +import assert from 'node:assert/strict'; +import * as Cesium from 'cesium'; +import * as geo from './index.js'; +import { DeveloperError } from './developerError.js'; + +const TRIALS = 150; + +function generator(seed) { + return () => { + seed = (seed + 0x6d2b79f5) | 0; + let t = Math.imul(seed ^ (seed >>> 15), 1 | seed); + t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t; + return ((t ^ (t >>> 14)) >>> 0) / 4294967296; + }; +} + +// Invalid values every argument may receive. NaN is a number to Check.typeOf.number +// and null an object to Check.typeOf.object, so both reach the arithmetic. +const INVALID = [undefined, null, Number.NaN, 'abc', 5, true, {}]; + +const vectorValues = { x: 1.5, y: -2, z: 3 }; +const VALID = { + num: (side, random) => [0, -0, 0.7, -2.5, 1e6][Math.floor(random() * 5)], + index3: (side, random) => Math.floor(random() * 5) - 1, + index4: (side, random) => Math.floor(random() * 6) - 1, + nonNegative: (side, random) => [0, 3, 6378137, -1][Math.floor(random() * 4)], + vec: (side) => + side === 'geo' ? { ...vectorValues } : new Cesium.Cartesian3(1.5, -2, 3), + zeroVec: (side) => + side === 'geo' ? { x: 0, y: 0, z: 0 } : new Cesium.Cartesian3(), + vec4: (side) => + side === 'geo' ? { x: 0, y: 0, z: 0, w: 0 } : new Cesium.Cartesian4(), + carto: (side) => + side === 'geo' + ? { longitude: 0.3, latitude: -0.2, height: 10 } + : new Cesium.Cartographic(0.3, -0.2, 10), + partialCarto: (side) => + side === 'geo' + ? { longitude: 0.3, latitude: -0.2 } + : { longitude: 0.3, latitude: -0.2 }, + quat: (side) => + side === 'geo' + ? { x: 0.1, y: 0.2, z: 0.3, w: 0.9 } + : new Cesium.Quaternion(0.1, 0.2, 0.3, 0.9), + mat3: (side) => + side === 'geo' + ? geo.matrix3.fromRotationZ(0.4) + : Cesium.Matrix3.fromRotationZ(0.4, new Cesium.Matrix3()), + mat4: (side) => + side === 'geo' + ? geo.transforms.eastNorthUpToFixedFrame({ x: 6378137, y: 1, z: 2 }) + : Cesium.Transforms.eastNorthUpToFixedFrame( + new Cesium.Cartesian3(6378137, 1, 2), + ), + hpr: (side) => + side === 'geo' + ? { heading: 0.3, pitch: -0.2, roll: 0.1 } + : new Cesium.HeadingPitchRoll(0.3, -0.2, 0.1), + partialHpr: () => ({ heading: 0.3 }), + ellipsoid: (side) => + side === 'geo' ? geo.ellipsoid.WGS84 : Cesium.Ellipsoid.WGS84, + coords2: (side, random) => + [[1, 2], [1, 2, 3], [], [1, 2, 3, 4]][Math.floor(random() * 4)], + coords3: (side, random) => + [[1, 2, 3], [1, 2], [], [1, 2, 3, 4, 5, 6], [1, 2, 'h']][ + Math.floor(random() * 5) + ], + nanVec: (side) => + side === 'geo' + ? { x: 1, y: Number.NaN, z: 2 } + : new Cesium.Cartesian3(1, Number.NaN, 2), + partialVec: () => ({ x: 1, y: 2 }), + stringHeight: () => '5', +}; + +/** One argument: usually a valid value of its kind, sometimes an invalid one. */ +function argument(kinds, side, random) { + if (random() < 0.4) return INVALID[Math.floor(random() * INVALID.length)]; + const kind = kinds[Math.floor(random() * kinds.length)]; + return VALID[kind](side, random); +} + +const W = Cesium.Ellipsoid.WGS84; +const C3 = Cesium.Cartesian3; +const V = ['vec', 'zeroVec', 'nanVec', 'partialVec']; +const N = ['num']; +const R = ['vec']; + +// [name, geo function, Cesium function, argument kinds per position] +const CASES = [ + ['math.toRadians', geo.math.toRadians, Cesium.Math.toRadians, [N]], + ['math.toDegrees', geo.math.toDegrees, Cesium.Math.toDegrees, [N]], + ['math.clamp', geo.math.clamp, Cesium.Math.clamp, [N, N, N]], + ['math.lerp', geo.math.lerp, Cesium.Math.lerp, [N, N, N]], + ['math.mod', geo.math.mod, Cesium.Math.mod, [N, ['num', 'nonNegative']]], + ['math.zeroToTwoPi', geo.math.zeroToTwoPi, Cesium.Math.zeroToTwoPi, [N]], + [ + 'math.negativePiToPi', + geo.math.negativePiToPi, + Cesium.Math.negativePiToPi, + [N], + ], + [ + 'math.equalsEpsilon', + geo.math.equalsEpsilon, + Cesium.Math.equalsEpsilon, + [N, N, N, N], + ], + + ['cartesian2.clone', geo.cartesian2.clone, Cesium.Cartesian2.clone, [R, R]], + [ + 'cartesian2.distance', + geo.cartesian2.distance, + Cesium.Cartesian2.distance, + [R, R], + ], + + ['cartesian3.clone', geo.cartesian3.clone, C3.clone, [V, R]], + [ + 'cartesian3.fromElements', + geo.cartesian3.fromElements, + C3.fromElements, + [N, N, N, R], + ], + ['cartesian3.magnitude', geo.cartesian3.magnitude, C3.magnitude, [V]], + [ + 'cartesian3.magnitudeSquared', + geo.cartesian3.magnitudeSquared, + C3.magnitudeSquared, + [V], + ], + ['cartesian3.distance', geo.cartesian3.distance, C3.distance, [V, V]], + [ + 'cartesian3.distanceSquared', + geo.cartesian3.distanceSquared, + C3.distanceSquared, + [V, V], + ], + ['cartesian3.normalize', geo.cartesian3.normalize, C3.normalize, [V, R]], + ['cartesian3.dot', geo.cartesian3.dot, C3.dot, [V, V]], + [ + 'cartesian3.multiplyComponents', + geo.cartesian3.multiplyComponents, + C3.multiplyComponents, + [V, V, R], + ], + ['cartesian3.add', geo.cartesian3.add, C3.add, [V, V, R]], + ['cartesian3.subtract', geo.cartesian3.subtract, C3.subtract, [V, V, R]], + [ + 'cartesian3.multiplyByScalar', + geo.cartesian3.multiplyByScalar, + C3.multiplyByScalar, + [V, N, R], + ], + [ + 'cartesian3.divideByScalar', + geo.cartesian3.divideByScalar, + C3.divideByScalar, + [V, N, R], + ], + ['cartesian3.negate', geo.cartesian3.negate, C3.negate, [V, R]], + ['cartesian3.lerp', geo.cartesian3.lerp, C3.lerp, [V, V, N, R]], + ['cartesian3.cross', geo.cartesian3.cross, C3.cross, [V, V, R]], + ['cartesian3.midpoint', geo.cartesian3.midpoint, C3.midpoint, [V, V, R]], + ['cartesian3.equals', geo.cartesian3.equals, C3.equals, [V, V]], + [ + 'cartesian3.equalsEpsilon', + geo.cartesian3.equalsEpsilon, + C3.equalsEpsilon, + [V, V, N, N], + ], + [ + 'cartesian3.fromRadians', + geo.cartesian3.fromRadians, + C3.fromRadians, + [N, N, ['num', 'stringHeight'], ['ellipsoid'], R], + ], + [ + 'cartesian3.fromDegrees', + geo.cartesian3.fromDegrees, + C3.fromDegrees, + [N, N, ['num', 'stringHeight'], ['ellipsoid'], R], + ], + [ + 'cartesian3.fromDegreesArray', + geo.cartesian3.fromDegreesArray, + C3.fromDegreesArray, + [['coords2'], ['ellipsoid']], + ], + [ + 'cartesian3.fromDegreesArrayHeights', + geo.cartesian3.fromDegreesArrayHeights, + C3.fromDegreesArrayHeights, + [['coords3'], ['ellipsoid']], + ], + + [ + 'cartographic.fromRadians', + geo.cartographic.fromRadians, + Cesium.Cartographic.fromRadians, + [N, N, N, ['carto']], + ], + [ + 'cartographic.fromDegrees', + geo.cartographic.fromDegrees, + Cesium.Cartographic.fromDegrees, + [N, N, N, ['carto']], + ], + [ + 'cartographic.fromCartesian', + geo.cartographic.fromCartesian, + Cesium.Cartographic.fromCartesian, + [V, ['ellipsoid'], ['carto']], + ], + [ + 'cartographic.toCartesian', + geo.cartographic.toCartesian, + Cesium.Cartographic.toCartesian, + [['carto', 'partialCarto'], ['ellipsoid'], R], + ], + [ + 'cartographic.clone', + geo.cartographic.clone, + Cesium.Cartographic.clone, + [['carto'], ['carto']], + ], + + [ + 'ellipsoid.create', + geo.ellipsoid.create, + (x, y, z) => new Cesium.Ellipsoid(x, y, z), + [['nonNegative'], ['nonNegative'], ['nonNegative']], + ], + [ + 'WGS84.geodeticSurfaceNormal', + geo.ellipsoid.WGS84.geodeticSurfaceNormal, + (c, r) => W.geodeticSurfaceNormal(c, r), + [V, R], + ], + [ + 'WGS84.geodeticSurfaceNormalCartographic', + geo.ellipsoid.WGS84.geodeticSurfaceNormalCartographic, + (c, r) => W.geodeticSurfaceNormalCartographic(c, r), + [['carto', 'partialCarto'], R], + ], + [ + 'WGS84.scaleToGeodeticSurface', + geo.ellipsoid.WGS84.scaleToGeodeticSurface, + (c, r) => W.scaleToGeodeticSurface(c, r), + [V, R], + ], + [ + 'WGS84.cartographicToCartesian', + geo.ellipsoid.WGS84.cartographicToCartesian, + (c, r) => W.cartographicToCartesian(c, r), + [['carto', 'partialCarto'], R], + ], + [ + 'WGS84.cartesianToCartographic', + geo.ellipsoid.WGS84.cartesianToCartographic, + (c, r) => W.cartesianToCartographic(c, r), + [V, ['carto']], + ], + + [ + 'matrix3.fromQuaternion', + geo.matrix3.fromQuaternion, + Cesium.Matrix3.fromQuaternion, + [['quat'], ['mat3']], + ], + [ + 'matrix3.fromRotationZ', + geo.matrix3.fromRotationZ, + Cesium.Matrix3.fromRotationZ, + [N, ['mat3']], + ], + [ + 'matrix3.multiplyByVector', + geo.matrix3.multiplyByVector, + Cesium.Matrix3.multiplyByVector, + [['mat3'], V, R], + ], + [ + 'matrix3.getColumn', + geo.matrix3.getColumn, + Cesium.Matrix3.getColumn, + [['mat3'], ['index3'], R], + ], + + [ + 'matrix4.clone', + geo.matrix4.clone, + Cesium.Matrix4.clone, + [['mat4'], ['mat4']], + ], + [ + 'matrix4.equals', + geo.matrix4.equals, + Cesium.Matrix4.equals, + [['mat4'], ['mat4']], + ], + [ + 'matrix4.multiply', + geo.matrix4.multiply, + Cesium.Matrix4.multiply, + [['mat4'], ['mat4'], ['mat4']], + ], + [ + 'matrix4.multiplyTransformation', + geo.matrix4.multiplyTransformation, + Cesium.Matrix4.multiplyTransformation, + [['mat4'], ['mat4'], ['mat4']], + ], + [ + 'matrix4.multiplyByPoint', + geo.matrix4.multiplyByPoint, + Cesium.Matrix4.multiplyByPoint, + [['mat4'], V, R], + ], + [ + 'matrix4.multiplyByPointAsVector', + geo.matrix4.multiplyByPointAsVector, + Cesium.Matrix4.multiplyByPointAsVector, + [['mat4'], V, R], + ], + [ + 'matrix4.getTranslation', + geo.matrix4.getTranslation, + Cesium.Matrix4.getTranslation, + [['mat4'], R], + ], + [ + 'matrix4.getMatrix3', + geo.matrix4.getMatrix3, + Cesium.Matrix4.getMatrix3, + [['mat4'], ['mat3']], + ], + [ + 'matrix4.getColumn', + geo.matrix4.getColumn, + Cesium.Matrix4.getColumn, + [['mat4'], ['index4'], ['vec4']], + ], + [ + 'matrix4.inverseTransformation', + geo.matrix4.inverseTransformation, + Cesium.Matrix4.inverseTransformation, + [['mat4'], ['mat4']], + ], + [ + 'matrix4.fromRotationTranslation', + geo.matrix4.fromRotationTranslation, + Cesium.Matrix4.fromRotationTranslation, + [['mat3'], R, ['mat4']], + ], + [ + 'matrix4.fromUniformScale', + geo.matrix4.fromUniformScale, + Cesium.Matrix4.fromUniformScale, + [N, ['mat4']], + ], + + [ + 'quaternion.fromAxisAngle', + geo.quaternion.fromAxisAngle, + Cesium.Quaternion.fromAxisAngle, + [V, N, ['quat']], + ], + [ + 'quaternion.fromRotationMatrix', + geo.quaternion.fromRotationMatrix, + Cesium.Quaternion.fromRotationMatrix, + [['mat3'], ['quat']], + ], + + [ + 'transforms.eastNorthUpToFixedFrame', + geo.transforms.eastNorthUpToFixedFrame, + Cesium.Transforms.eastNorthUpToFixedFrame, + [V, ['ellipsoid'], ['mat4']], + ], + [ + 'transforms.headingPitchRollToFixedFrame', + geo.transforms.headingPitchRollToFixedFrame, + Cesium.Transforms.headingPitchRollToFixedFrame, + [V, ['hpr', 'partialHpr'], ['ellipsoid'], [], ['mat4']], + ], +]; + +function outcome(fn, args) { + try { + fn(...args); + return 'returned'; + } catch (error) { + return `${error?.name}: ${error?.message}`; + } +} + +test('every ported function throws exactly where Cesium does, with the same name and message', () => { + let thrown = 0; + let checked = 0; + const mismatches = []; + CASES.forEach(([name, ours, theirs, kinds], caseIndex) => { + for (let trial = 0; trial < TRIALS; trial++) { + const seed = caseIndex * 100003 + trial; + // Trial 0 passes no arguments at all; the rest draw from the pools. + const build = (side) => { + if (trial === 0) return []; + const random = generator(seed); + return kinds.map((k) => + k.length ? argument(k, side, random) : undefined, + ); + }; + const a = outcome(ours, build('geo')); + const b = outcome(theirs, build('cesium')); + if (a !== b) + mismatches.push(`${name} (seed ${seed}): geo ${a} | Cesium ${b}`); + checked++; + if (b !== 'returned') thrown++; + } + }); + assert.deepEqual(mismatches, []); + // The pools must actually reach the checks. + assert.ok(thrown > checked / 4, `only ${thrown} of ${checked} cases threw`); +}); + +test('the checks throw a DeveloperError with Cesium message text', () => { + const cases = [ + [() => geo.math.toRadians(undefined), 'degrees is required.'], + [() => geo.math.mod(5, 0), 'divisor cannot be 0.'], + [ + () => geo.math.clamp('1', 0, 1), + 'Expected value to be typeof number, actual typeof was string', + ], + [ + () => geo.cartesian3.fromDegrees(undefined, 1), + 'Expected longitude to be typeof number, actual typeof was undefined', + ], + [ + () => geo.cartesian3.add({ x: 0, y: 0, z: 0 }, { x: 0, y: 0, z: 0 }), + 'Expected result to be typeof object, actual typeof was undefined', + ], + [ + () => geo.cartesian3.fromDegreesArray([1]), + 'the number of coordinates must be a multiple of 2 and at least 2', + ], + [ + () => geo.cartesian3.fromDegreesArray(), + 'coordinates is required, actual value was undefined', + ], + [ + () => geo.ellipsoid.create(-1, 1, 1), + 'Expected x to be greater than or equal to 0, actual value was -1', + ], + [ + () => geo.matrix3.getColumn(geo.matrix3.create(), 3, {}), + 'Expected index to be less than or equal to 2, actual value was 3', + ], + [ + () => + geo.transforms.eastNorthUpToFixedFrame({ x: 1, y: Number.NaN, z: 0 }), + 'origin has a NaN component', + ], + [ + () => geo.cartographic.fromCartesian({ x: 1, y: 2, z: 3 }, {}), + 'oneOverRadii is required.', + ], + ]; + for (const [fn, message] of cases) { + assert.throws(fn, (error) => { + assert.ok(error instanceof DeveloperError); + assert.equal(error.name, 'DeveloperError'); + assert.equal(error.message, message); + return true; + }); + } + // NaN is typeof number, so Check.typeOf.number lets it through (as in Cesium). + assert.ok(Number.isNaN(geo.math.clamp(Number.NaN, 0, 1))); + assert.ok(Number.isNaN(Cesium.Math.clamp(Number.NaN, 0, 1))); +}); diff --git a/src/geo/debug.js b/src/geo/debug.js new file mode 100644 index 000000000..b3eed7c23 --- /dev/null +++ b/src/geo/debug.js @@ -0,0 +1,12 @@ +// Whether src/geo runs Cesium's debug argument checks. +// +// Cesium wraps its argument checks in `//>>includeStart('debug')` pragmas: they +// are live in its Source modules (Node tests, the Vite dev server) and stripped +// from Build/Cesium, which vite-plugin-cesium copies into production builds. +// src/geo does the same: every Check call and every debug throw sits in an +// `if (DEBUG)` block. DEBUG is true in Node (no import.meta.env) and in Vite +// dev; `vite build` folds it to false and the minifier drops the blocks, so a +// production build returns what release Cesium returns (NaN, not a throw, for +// a degenerate input). + +export const DEBUG = import.meta.env?.PROD !== true; diff --git a/src/geo/developerError.js b/src/geo/developerError.js new file mode 100644 index 000000000..459a8bacf --- /dev/null +++ b/src/geo/developerError.js @@ -0,0 +1,13 @@ +// Ported from CesiumJS 1.138 (Apache-2.0) — DeveloperError.js +// +// The error Cesium throws for invalid arguments. Its debug argument checks are +// live in Cesium's Source modules (Node tests, the Vite dev server) and +// stripped from production builds; src/geo throws the same name and message +// wherever Cesium does while DEBUG is on (see debug.js). + +export class DeveloperError extends Error { + constructor(message) { + super(message); + this.name = 'DeveloperError'; + } +} diff --git a/src/geo/ellipsoid.js b/src/geo/ellipsoid.js new file mode 100644 index 000000000..2da0bcaaf --- /dev/null +++ b/src/geo/ellipsoid.js @@ -0,0 +1,180 @@ +// Ported from CesiumJS 1.138 (Apache-2.0) — Ellipsoid.js +// +// A quadratic surface x²/a² + y²/b² + z²/c² = 1 as a plain, frozen object. +// Every method takes an optional `result` and returns it (a new object when it +// is absent); `result` may be the input. The algorithms and the operation order +// match Cesium so the outputs are bit-identical. With DEBUG on (development +// and tests; see debug.js) invalid arguments throw the same DeveloperError as +// Cesium's debug checks. + +import { DEBUG } from './debug.js'; +import { Check } from './check.js'; +import { create as createCartographic } from './cartographic.js'; +import { DeveloperError } from './developerError.js'; +import { EPSILON14 } from './math.js'; +import { scaleToGeodeticSurface } from './scaleToGeodeticSurface.js'; +import { + add, + create as vec, + equalsEpsilon, + ZERO, + divideByScalar, + dot, + magnitude, + multiplyByScalar, + multiplyComponents, + normalize, + subtract, +} from './cartesian3.js'; + +/** Cesium's CesiumMath.EPSILON1, the squared norm below which a point counts as the centre. */ +const CENTER_TOLERANCE_SQUARED = 0.1; + +/** + * Builds an ellipsoid with radii `x`, `y`, `z` (Cesium's `new Ellipsoid(x, y, z)`). + * The returned object and its vector fields are frozen (deeply; Cesium freezes + * only the WGS84 object itself, not its radii vectors). + */ +export function create(x, y, z) { + x = x ?? 0.0; + y = y ?? 0.0; + z = z ?? 0.0; + if (DEBUG) { + Check.typeOf.number.greaterThanOrEquals('x', x, 0.0); + Check.typeOf.number.greaterThanOrEquals('y', y, 0.0); + Check.typeOf.number.greaterThanOrEquals('z', z, 0.0); + } + + const radii = Object.freeze(vec(x, y, z)); + const radiiSquared = Object.freeze(vec(x * x, y * y, z * z)); + const radiiToTheFourth = Object.freeze( + vec(x * x * x * x, y * y * y * y, z * z * z * z), + ); + const oneOverRadii = Object.freeze( + vec( + x === 0.0 ? 0.0 : 1.0 / x, + y === 0.0 ? 0.0 : 1.0 / y, + z === 0.0 ? 0.0 : 1.0 / z, + ), + ); + const oneOverRadiiSquared = Object.freeze( + vec( + x === 0.0 ? 0.0 : 1.0 / (x * x), + y === 0.0 ? 0.0 : 1.0 / (y * y), + z === 0.0 ? 0.0 : 1.0 / (z * z), + ), + ); + const centerToleranceSquared = CENTER_TOLERANCE_SQUARED; + + /** The normal to the ellipsoid surface at a geodetic position. */ + function geodeticSurfaceNormalCartographic(cartographic, result) { + if (DEBUG) { + Check.typeOf.object('cartographic', cartographic); + } + + const longitude = cartographic.longitude; + const latitude = cartographic.latitude; + const cosLatitude = Math.cos(latitude); + + const nx = cosLatitude * Math.cos(longitude); + const ny = cosLatitude * Math.sin(longitude); + const nz = Math.sin(latitude); + + if (result == null) result = vec(); + result.x = nx; + result.y = ny; + result.z = nz; + return normalize(result, result); + } + + /** The normal to the ellipsoid surface at a position, or `undefined` at the centre. */ + function geodeticSurfaceNormal(cartesian, result) { + if (DEBUG) { + Check.typeOf.object('cartesian', cartesian); + if (isNaN(cartesian.x) || isNaN(cartesian.y) || isNaN(cartesian.z)) { + throw new DeveloperError('cartesian has a NaN component'); + } + } + + if (equalsEpsilon(cartesian, ZERO, EPSILON14)) { + return undefined; + } + if (result == null) result = vec(); + result = multiplyComponents(cartesian, oneOverRadiiSquared, result); + return normalize(result, result); + } + + /** Scales a position onto the geodetic surface, or `undefined` near the centre. */ + function scaleToSurface(cartesian, result) { + return scaleToGeodeticSurface( + cartesian, + oneOverRadii, + oneOverRadiiSquared, + centerToleranceSquared, + result, + ); + } + + const cartographicToCartesianNormal = vec(); + const cartographicToCartesianK = vec(); + + /** Converts a geodetic position to Earth-fixed Cartesian coordinates. */ + function cartographicToCartesian(cartographic, result) { + const n = cartographicToCartesianNormal; + const k = cartographicToCartesianK; + geodeticSurfaceNormalCartographic(cartographic, n); + multiplyComponents(radiiSquared, n, k); + const gamma = Math.sqrt(dot(n, k)); + divideByScalar(k, gamma, k); + multiplyByScalar(n, cartographic.height, n); + if (result == null) result = vec(); + return add(k, n, result); + } + + const cartesianToCartographicN = vec(); + const cartesianToCartographicP = vec(); + const cartesianToCartographicH = vec(); + + /** Converts an Earth-fixed position to geodetic coordinates, or `undefined` at the centre. */ + function cartesianToCartographic(cartesian, result) { + const p = scaleToSurface(cartesian, cartesianToCartographicP); + if (p === undefined) return undefined; + + const n = geodeticSurfaceNormal(p, cartesianToCartographicN); + const h = subtract(cartesian, p, cartesianToCartographicH); + + const longitude = Math.atan2(n.y, n.x); + const latitude = Math.asin(n.z); + const height = Math.sign(dot(h, cartesian)) * magnitude(h); + + if (result == null) return createCartographic(longitude, latitude, height); + result.longitude = longitude; + result.latitude = latitude; + result.height = height; + return result; + } + + return Object.freeze({ + radii, + radiiSquared, + radiiToTheFourth, + oneOverRadii, + oneOverRadiiSquared, + minimumRadius: Math.min(x, y, z), + maximumRadius: Math.max(x, y, z), + centerToleranceSquared, + geodeticSurfaceNormal, + geodeticSurfaceNormalCartographic, + scaleToGeodeticSurface: scaleToSurface, + cartographicToCartesian, + cartesianToCartographic, + }); +} + +/** + * The WGS84 ellipsoid. For src/geo only: never pass it into a Cesium function + * or constructor (pickEllipsoid, EllipsoidalOccluder, …). It lacks Cesium's + * private fields and methods such as `_centerToleranceSquared` and + * `transformPositionToScaledSpace`; use `Cesium.Ellipsoid.WGS84` there. + */ +export const WGS84 = create(6378137.0, 6378137.0, 6356752.3142451793); diff --git a/src/geo/ellipsoid.test.mjs b/src/geo/ellipsoid.test.mjs new file mode 100644 index 000000000..5389c8a2b --- /dev/null +++ b/src/geo/ellipsoid.test.mjs @@ -0,0 +1,148 @@ +// src/geo/ellipsoid.test.mjs +import { test } from 'node:test'; +import assert from 'node:assert/strict'; +import * as Cesium from 'cesium'; +import { WGS84, create } from './ellipsoid.js'; + +const close = (a, b, tol, label) => + assert.ok( + Math.abs(a - b) <= tol * Math.max(1, Math.abs(b)), + `${label}: ${a} vs ${b}`, + ); +const closeVec = (a, b, tol, label) => { + for (const k of ['x', 'y', 'z']) close(a[k], b[k], tol, `${label}.${k}`); +}; +const POINTS = [ + [6378137, 0, 0], + [0, 6378137, 0], + [0, 0, 6356752.314245], + [1334762.2, -4654869.1, 4137287.9], + [-2694045, -4293642, 3857878], + [1e3, 1e3, 1e3], + [0.01, 0.01, 0.01], + [7e6, -2e6, 1e6], + [4510000, 4510000, 1], + [1, 1, 6400000], +]; + +test('WGS84 constants equal Cesium', () => { + const c = Cesium.Ellipsoid.WGS84; + for (const k of [ + 'radii', + 'radiiSquared', + 'radiiToTheFourth', + 'oneOverRadii', + 'oneOverRadiiSquared', + ]) { + closeVec(WGS84[k], c[k], 0, k); + } + assert.equal(WGS84.minimumRadius, c.minimumRadius); + assert.equal(WGS84.maximumRadius, c.maximumRadius); +}); + +test('surface normals, scaling and conversions equal Cesium', () => { + const c = Cesium.Ellipsoid.WGS84; + for (const [x, y, z] of POINTS) { + const p = { x, y, z }; + const cp = new Cesium.Cartesian3(x, y, z); + closeVec( + WGS84.geodeticSurfaceNormal(p), + c.geodeticSurfaceNormal(cp), + 1e-12, + 'normal', + ); + const ours = WGS84.scaleToGeodeticSurface(p); + const theirs = c.scaleToGeodeticSurface(cp); + if (theirs === undefined) assert.equal(ours, undefined); + else closeVec(ours, theirs, 1e-12, 'scale'); + const carto = WGS84.cartesianToCartographic(p); + const cCarto = c.cartesianToCartographic(cp); + if (cCarto === undefined) { + assert.equal(carto, undefined); + } else { + close(carto.longitude, cCarto.longitude, 1e-12, 'lon'); + close(carto.latitude, cCarto.latitude, 1e-12, 'lat'); + close(carto.height, cCarto.height, 1e-9, 'height'); + closeVec( + WGS84.cartographicToCartesian(carto), + c.cartographicToCartesian(cCarto), + 1e-12, + 'roundtrip', + ); + } + } +}); + +test('result parameters are filled in place and may alias the input', () => { + const p = { x: 1334762.2, y: -4654869.1, z: 4137287.9 }; + const out = { x: 0, y: 0, z: 0 }; + assert.equal(WGS84.geodeticSurfaceNormal(p, out), out); + const alias = { ...p }; + WGS84.geodeticSurfaceNormal(alias, alias); + closeVec(alias, out, 0, 'alias'); +}); + +test('a custom ellipsoid behaves like Cesium', () => { + const ours = create(10, 20, 30); + const theirs = new Cesium.Ellipsoid(10, 20, 30); + const p = { x: 5, y: 7, z: 9 }; + closeVec( + ours.scaleToGeodeticSurface(p), + theirs.scaleToGeodeticSurface(new Cesium.Cartesian3(5, 7, 9)), + 1e-12, + 'custom', + ); +}); + +const sameXyz = (ours, theirs, label) => { + if (theirs === undefined) return assert.equal(ours, undefined, label); + for (const k of ['x', 'y', 'z']) { + assert.ok( + Object.is(ours[k], theirs[k]), + `${label}.${k}: ${ours[k]} vs ${theirs[k]}`, + ); + } +}; + +test('the centre, the poles and the antimeridian equal Cesium bit for bit', () => { + const c = Cesium.Ellipsoid.WGS84; + for (const [x, y, z] of [ + [0, 0, 0], + [-0, 0, -0], + [1e-15, 0, 0], + [0, 0, 6356752.3142451793], + [0, 0, -6356752.3142451793], + [-6378137, 0, 0], + [-6378137, -0, 0], + [-7e6, -0, 1e5], + ]) { + const p = { x, y, z }; + const cp = new Cesium.Cartesian3(x, y, z); + const label = `(${Object.is(x, -0) ? '-0' : x}, ${Object.is(y, -0) ? '-0' : y}, ${z})`; + sameXyz( + WGS84.geodeticSurfaceNormal(p), + c.geodeticSurfaceNormal(cp), + `normal ${label}`, + ); + sameXyz( + WGS84.scaleToGeodeticSurface(p), + c.scaleToGeodeticSurface(cp), + `scale ${label}`, + ); + const ours = WGS84.cartesianToCartographic(p); + const theirs = c.cartesianToCartographic(cp); + if (theirs === undefined) { + assert.equal(ours, undefined, label); + continue; + } + for (const k of ['longitude', 'latitude', 'height']) { + assert.ok( + Object.is(ours[k], theirs[k]), + `${label}.${k}: ${ours[k]} vs ${theirs[k]}`, + ); + } + } + assert.equal(WGS84.geodeticSurfaceNormal({ x: 0, y: 0, z: 0 }), undefined); + assert.equal(WGS84.cartesianToCartographic({ x: 0, y: 0, z: 0 }), undefined); + assert.equal(WGS84.scaleToGeodeticSurface({ x: 0, y: 0, z: 0 }), undefined); +}); diff --git a/src/geo/index.js b/src/geo/index.js new file mode 100644 index 000000000..c86b10401 --- /dev/null +++ b/src/geo/index.js @@ -0,0 +1,15 @@ +// Engine-neutral geometry (WGS84, vectors, matrices, local frames) ported from +// CesiumJS 1.138 (Apache-2.0). Never imports Cesium; see docs/CODE-BOUNDARIES.md. +// +// import * as geo from '../geo/index.js'; +// geo.cartesian3.fromDegrees(lon, lat, height); + +export * as math from './math.js'; +export * as ellipsoid from './ellipsoid.js'; +export * as cartographic from './cartographic.js'; +export * as cartesian2 from './cartesian2.js'; +export * as cartesian3 from './cartesian3.js'; +export * as matrix3 from './matrix3.js'; +export * as matrix4 from './matrix4.js'; +export * as quaternion from './quaternion.js'; +export * as transforms from './transforms.js'; diff --git a/src/geo/index.test.mjs b/src/geo/index.test.mjs new file mode 100644 index 000000000..e15050449 --- /dev/null +++ b/src/geo/index.test.mjs @@ -0,0 +1,22 @@ +// src/geo/index.test.mjs +import { test } from 'node:test'; +import assert from 'node:assert/strict'; +import * as geo from './index.js'; + +test('the public entry exposes every namespace', () => { + for (const name of [ + 'math', + 'ellipsoid', + 'cartographic', + 'cartesian2', + 'cartesian3', + 'matrix3', + 'matrix4', + 'quaternion', + 'transforms', + ]) { + assert.equal(typeof geo[name], 'object', name); + } + assert.equal(typeof geo.cartesian3.fromDegrees, 'function'); + assert.equal(geo.ellipsoid.WGS84.maximumRadius, 6378137); +}); diff --git a/src/geo/math.js b/src/geo/math.js new file mode 100644 index 000000000..3e77f5eef --- /dev/null +++ b/src/geo/math.js @@ -0,0 +1,133 @@ +// Ported from CesiumJS 1.138 (Apache-2.0) — Math.js +// +// Engine-neutral angle maths. Formulas and constants match Cesium.Math bit for +// bit. With DEBUG on (development and tests; see debug.js) invalid arguments +// throw the same DeveloperError as Cesium's debug checks; production builds +// skip the checks, as release Cesium does. + +import { DEBUG } from './debug.js'; +import { Check } from './check.js'; +import { DeveloperError } from './developerError.js'; + +export const EPSILON7 = 0.0000001; +export const EPSILON10 = 0.0000000001; +export const EPSILON12 = 0.000000000001; +export const EPSILON14 = 0.00000000000001; +export const PI_OVER_TWO = Math.PI / 2.0; +export const TWO_PI = 2.0 * Math.PI; +export const RADIANS_PER_DEGREE = Math.PI / 180.0; +export const DEGREES_PER_RADIAN = 180.0 / Math.PI; + +/** Converts degrees to radians. */ +export function toRadians(degrees) { + if (DEBUG) { + if (degrees == null) { + throw new DeveloperError('degrees is required.'); + } + } + return degrees * RADIANS_PER_DEGREE; +} + +/** Converts radians to degrees. */ +export function toDegrees(radians) { + if (DEBUG) { + if (radians == null) { + throw new DeveloperError('radians is required.'); + } + } + return radians * DEGREES_PER_RADIAN; +} + +/** Constrains a value to lie between two values. */ +export function clamp(value, min, max) { + if (DEBUG) { + Check.typeOf.number('value', value); + Check.typeOf.number('min', min); + Check.typeOf.number('max', max); + } + + return value < min ? min : value > max ? max : value; +} + +/** Computes the linear interpolation of two values. */ +export function lerp(p, q, time) { + return (1.0 - time) * p + time * q; +} + +/** The modulo operation that also works for negative dividends. */ +export function mod(m, n) { + if (DEBUG) { + if (m == null) { + throw new DeveloperError('m is required.'); + } + if (n == null) { + throw new DeveloperError('n is required.'); + } + if (n === 0.0) { + throw new DeveloperError('divisor cannot be 0.'); + } + } + if (Math.sign(m) === Math.sign(n) && Math.abs(m) < Math.abs(n)) { + // Early exit if the input does not need to be modded. This avoids + // unnecessary math which could introduce floating point error. + return m; + } + return ((m % n) + n) % n; +} + +/** Produces an angle in the range 0 <= angle <= 2Pi equivalent to the provided angle. */ +export function zeroToTwoPi(angle) { + if (DEBUG) { + if (angle == null) { + throw new DeveloperError('angle is required.'); + } + } + if (angle >= 0 && angle <= TWO_PI) { + // Early exit if the input is already inside the range. This avoids + // unnecessary math which could introduce floating point error. + return angle; + } + const result = mod(angle, TWO_PI); + if (Math.abs(result) < EPSILON14 && Math.abs(angle) > EPSILON14) { + return TWO_PI; + } + return result; +} + +/** Produces an angle in the range -Pi <= angle <= Pi equivalent to the provided angle. */ +export function negativePiToPi(angle) { + if (DEBUG) { + if (angle == null) { + throw new DeveloperError('angle is required.'); + } + } + if (angle >= -Math.PI && angle <= Math.PI) { + // Early exit if the input is already inside the range. This avoids + // unnecessary math which could introduce floating point error. + return angle; + } + return zeroToTwoPi(angle + Math.PI) - Math.PI; +} + +/** + * Whether two values are equal within an absolute or a relative tolerance + * (the absolute one defaults to the relative one). + */ +export function equalsEpsilon(left, right, relativeEpsilon, absoluteEpsilon) { + if (DEBUG) { + if (left == null) { + throw new DeveloperError('left is required.'); + } + if (right == null) { + throw new DeveloperError('right is required.'); + } + } + + relativeEpsilon = relativeEpsilon ?? 0.0; + absoluteEpsilon = absoluteEpsilon ?? relativeEpsilon; + const absDiff = Math.abs(left - right); + return ( + absDiff <= absoluteEpsilon || + absDiff <= relativeEpsilon * Math.max(Math.abs(left), Math.abs(right)) + ); +} diff --git a/src/geo/math.test.mjs b/src/geo/math.test.mjs new file mode 100644 index 000000000..d2cecc266 --- /dev/null +++ b/src/geo/math.test.mjs @@ -0,0 +1,121 @@ +// src/geo/math.test.mjs +import { test } from 'node:test'; +import assert from 'node:assert/strict'; +import * as Cesium from 'cesium'; +import * as geo from './math.js'; + +const ANGLES = [ + 0, + 1e-15, + -1e-15, + 0.5, + -0.5, + Math.PI, + -Math.PI, + Math.PI + 1e-12, + 2 * Math.PI, + -2 * Math.PI, + 2 * Math.PI - 1e-15, + 2 * Math.PI + 1e-15, + 3 * Math.PI, + -3 * Math.PI, + 7.5, + -7.5, + 100, + -100, + 1e6, + -1e6, +]; + +test('constants equal Cesium.Math', () => { + for (const name of [ + 'PI_OVER_TWO', + 'TWO_PI', + 'EPSILON7', + 'EPSILON10', + 'EPSILON12', + 'EPSILON14', + 'RADIANS_PER_DEGREE', + 'DEGREES_PER_RADIAN', + ]) { + assert.equal(geo[name], Cesium.Math[name], name); + } +}); + +test('angle conversions equal Cesium bit for bit', () => { + for (const value of [ + 0, 1, -1, 45, 90, 180, 360, 123.456789, -73.9857, 1e-9, + ]) { + assert.equal(geo.toRadians(value), Cesium.Math.toRadians(value)); + assert.equal(geo.toDegrees(value), Cesium.Math.toDegrees(value)); + } +}); + +test('clamp, lerp and mod equal Cesium', () => { + for (const [v, lo, hi] of [ + [5, 0, 1], + [-5, 0, 1], + [0.5, 0, 1], + [0, 0, 0], + ]) { + assert.equal(geo.clamp(v, lo, hi), Cesium.Math.clamp(v, lo, hi)); + } + for (const t of [0, 0.25, 1, 1.5, -0.5]) { + assert.equal(geo.lerp(2, 10, t), Cesium.Math.lerp(2, 10, t)); + } + for (const [m, n] of [ + [5, 3], + [-5, 3], + [5, -3], + [-5, -3], + [0, 3], + [2, 3], + [-2, 3], + ]) { + assert.equal(geo.mod(m, n), Cesium.Math.mod(m, n), `${m} mod ${n}`); + } +}); + +test('angle normalisation equals Cesium at the edges', () => { + for (const a of ANGLES) { + assert.equal( + geo.zeroToTwoPi(a), + Cesium.Math.zeroToTwoPi(a), + `zeroToTwoPi(${a})`, + ); + assert.equal( + geo.negativePiToPi(a), + Cesium.Math.negativePiToPi(a), + `negativePiToPi(${a})`, + ); + } +}); + +test('the module does not import cesium', async () => { + const { readFile } = await import('node:fs/promises'); + const source = await readFile(new URL('./math.js', import.meta.url), 'utf8'); + assert.doesNotMatch(source, /from ['"]cesium['"]/); +}); + +test('equalsEpsilon equals Cesium, including the default epsilons', () => { + const CASES = [ + [1, 1], + [1, 1 + 1e-8, 1e-7], + [1, 1 + 1e-6, 1e-7], + [1e6, 1e6 + 0.05, 1e-7], + [1e6, 1e6 + 0.5, 1e-7], + [1e6, 1e6 + 0.5, 1e-7, 1], + [0, 1e-15, 1e-14], + [0, 1e-13, 1e-14], + [-3, 3, 0.1, 6], + [Number.NaN, Number.NaN, 1], + [1, 2], + ]; + for (const [l, r, rel, abs] of CASES) { + assert.equal( + geo.equalsEpsilon(l, r, rel, abs), + Cesium.Math.equalsEpsilon(l, r, rel, abs), + `${l} ~ ${r} (${rel}, ${abs})`, + ); + } +}); diff --git a/src/geo/matrix.test.mjs b/src/geo/matrix.test.mjs new file mode 100644 index 000000000..b4db952b1 --- /dev/null +++ b/src/geo/matrix.test.mjs @@ -0,0 +1,152 @@ +// src/geo/matrix.test.mjs +import { test } from 'node:test'; +import assert from 'node:assert/strict'; +import * as Cesium from 'cesium'; +import * as m3 from './matrix3.js'; +import * as m4 from './matrix4.js'; +import * as quat from './quaternion.js'; + +const arr = (m, n) => Array.from({ length: n }, (_, i) => m[i]); +const sameArr = (ours, theirs, n, label) => + assert.deepEqual(arr(ours, n), arr(theirs, n), label); +const close = (a, b, tol, label) => + assert.ok(Math.abs(a - b) <= tol, `${label}: ${a} vs ${b}`); + +const T = Cesium.Matrix4.fromTranslationQuaternionRotationScale( + new Cesium.Cartesian3(10, -20, 30), + Cesium.Quaternion.fromAxisAngle( + Cesium.Cartesian3.normalize( + new Cesium.Cartesian3(1, 2, 3), + new Cesium.Cartesian3(), + ), + 0.7, + ), + new Cesium.Cartesian3(1, 1, 1), +); +const U = Cesium.Matrix4.fromRotationTranslation( + Cesium.Matrix3.fromRotationZ(1.1), + new Cesium.Cartesian3(5, 6, 7), +); +const plain4 = (m) => arr(m, 16); + +test('matrix4 products, points and inverses equal Cesium', () => { + sameArr( + m4.multiply(plain4(T), plain4(U), m4.create()), + Cesium.Matrix4.multiply(T, U, new Cesium.Matrix4()), + 16, + 'multiply', + ); + sameArr( + m4.multiplyTransformation(plain4(T), plain4(U), m4.create()), + Cesium.Matrix4.multiplyTransformation(T, U, new Cesium.Matrix4()), + 16, + 'multiplyTransformation', + ); + const p = new Cesium.Cartesian3(1.5, -2, 3); + const ours = m4.multiplyByPoint( + plain4(T), + { x: 1.5, y: -2, z: 3 }, + { x: 0, y: 0, z: 0 }, + ); + const theirs = Cesium.Matrix4.multiplyByPoint(T, p, new Cesium.Cartesian3()); + assert.deepEqual({ ...ours }, { x: theirs.x, y: theirs.y, z: theirs.z }); + const v = m4.multiplyByPointAsVector( + plain4(T), + { x: 1.5, y: -2, z: 3 }, + { x: 0, y: 0, z: 0 }, + ); + const cv = Cesium.Matrix4.multiplyByPointAsVector( + T, + p, + new Cesium.Cartesian3(), + ); + assert.deepEqual({ ...v }, { x: cv.x, y: cv.y, z: cv.z }); + const inv = m4.inverseTransformation(plain4(T), m4.create()); + const cinv = Cesium.Matrix4.inverseTransformation(T, new Cesium.Matrix4()); + for (let i = 0; i < 16; i++) close(inv[i], cinv[i], 1e-12, `inv ${i}`); + const t = m4.getTranslation(plain4(T), { x: 0, y: 0, z: 0 }); + assert.deepEqual({ ...t }, { x: 10, y: -20, z: 30 }); + sameArr( + m4.getMatrix3(plain4(T), m3.create()), + Cesium.Matrix4.getMatrix3(T, new Cesium.Matrix3()), + 9, + 'getMatrix3', + ); + assert.ok(m4.equals(m4.clone(plain4(T)), plain4(T))); + assert.ok(Object.isFrozen(m4.IDENTITY)); + sameArr( + m4.fromUniformScale(2.5, m4.create()), + Cesium.Matrix4.fromUniformScale(2.5, new Cesium.Matrix4()), + 16, + 'scale', + ); +}); + +test('matrix3, quaternion and rotation-translation equal Cesium', () => { + sameArr( + m3.fromRotationZ(1.1, m3.create()), + Cesium.Matrix3.fromRotationZ(1.1, new Cesium.Matrix3()), + 9, + 'rotZ', + ); + const axis = { + x: 0.2672612419124244, + y: 0.5345224838248488, + z: 0.8017837257372732, + }; + const q = quat.fromAxisAngle(axis, 0.7, quat.create()); + const cq = Cesium.Quaternion.fromAxisAngle( + new Cesium.Cartesian3(axis.x, axis.y, axis.z), + 0.7, + new Cesium.Quaternion(), + ); + assert.deepEqual({ ...q }, { x: cq.x, y: cq.y, z: cq.z, w: cq.w }); + sameArr( + m3.fromQuaternion(q, m3.create()), + Cesium.Matrix3.fromQuaternion(cq, new Cesium.Matrix3()), + 9, + 'fromQuaternion', + ); + const back = quat.fromRotationMatrix( + m3.fromQuaternion(q, m3.create()), + quat.create(), + ); + const cback = Cesium.Quaternion.fromRotationMatrix( + Cesium.Matrix3.fromQuaternion(cq, new Cesium.Matrix3()), + new Cesium.Quaternion(), + ); + for (const k of ['x', 'y', 'z', 'w']) + close(back[k], cback[k], 1e-15, `q.${k}`); + const mv = m3.multiplyByVector( + m3.fromRotationZ(1.1, m3.create()), + { x: 1, y: 2, z: 3 }, + { x: 0, y: 0, z: 0 }, + ); + const cmv = Cesium.Matrix3.multiplyByVector( + Cesium.Matrix3.fromRotationZ(1.1, new Cesium.Matrix3()), + new Cesium.Cartesian3(1, 2, 3), + new Cesium.Cartesian3(), + ); + assert.deepEqual({ ...mv }, { x: cmv.x, y: cmv.y, z: cmv.z }); + const col = m3.getColumn(m3.fromRotationZ(1.1, m3.create()), 1, { + x: 0, + y: 0, + z: 0, + }); + const ccol = Cesium.Matrix3.getColumn( + Cesium.Matrix3.fromRotationZ(1.1, new Cesium.Matrix3()), + 1, + new Cesium.Cartesian3(), + ); + assert.deepEqual({ ...col }, { x: ccol.x, y: ccol.y, z: ccol.z }); + sameArr( + m4.fromRotationTranslation( + m3.fromRotationZ(1.1, m3.create()), + { x: 5, y: 6, z: 7 }, + m4.create(), + ), + U, + 16, + 'fromRotationTranslation', + ); +}); diff --git a/src/geo/matrix3.js b/src/geo/matrix3.js new file mode 100644 index 000000000..09550e6cc --- /dev/null +++ b/src/geo/matrix3.js @@ -0,0 +1,121 @@ +// Ported from CesiumJS 1.138 (Apache-2.0) — Matrix3.js +// +// 3x3 matrices as plain arrays of 9 numbers in Cesium's column-major index +// order (m[column * 3 + row]), so Cesium accepts them wherever it reads +// m[0]..m[8]. Only the subset gev uses is ported. With DEBUG on (development +// and tests; see debug.js) invalid arguments throw the same DeveloperError as +// Cesium's debug checks. + +import { DEBUG } from './debug.js'; +import { Check } from './check.js'; + +/** A matrix of zeros (Cesium's `new Matrix3()`). */ +export function create() { + return [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]; +} + +/** The rotation matrix of a unit quaternion. */ +export function fromQuaternion(quaternion, result) { + if (DEBUG) { + Check.typeOf.object('quaternion', quaternion); + } + + const x2 = quaternion.x * quaternion.x; + const xy = quaternion.x * quaternion.y; + const xz = quaternion.x * quaternion.z; + const xw = quaternion.x * quaternion.w; + const y2 = quaternion.y * quaternion.y; + const yz = quaternion.y * quaternion.z; + const yw = quaternion.y * quaternion.w; + const z2 = quaternion.z * quaternion.z; + const zw = quaternion.z * quaternion.w; + const w2 = quaternion.w * quaternion.w; + + const m00 = x2 - y2 - z2 + w2; + const m01 = 2.0 * (xy - zw); + const m02 = 2.0 * (xz + yw); + + const m10 = 2.0 * (xy + zw); + const m11 = -x2 + y2 - z2 + w2; + const m12 = 2.0 * (yz - xw); + + const m20 = 2.0 * (xz - yw); + const m21 = 2.0 * (yz + xw); + const m22 = -x2 - y2 + z2 + w2; + + if (result == null) result = create(); + result[0] = m00; + result[1] = m10; + result[2] = m20; + result[3] = m01; + result[4] = m11; + result[5] = m21; + result[6] = m02; + result[7] = m12; + result[8] = m22; + return result; +} + +/** A rotation by `angle` radians around the z-axis. */ +export function fromRotationZ(angle, result) { + if (DEBUG) { + Check.typeOf.number('angle', angle); + } + + const cosAngle = Math.cos(angle); + const sinAngle = Math.sin(angle); + + if (result == null) result = create(); + result[0] = cosAngle; + result[1] = sinAngle; + result[2] = 0.0; + result[3] = -sinAngle; + result[4] = cosAngle; + result[5] = 0.0; + result[6] = 0.0; + result[7] = 0.0; + result[8] = 1.0; + return result; +} + +/** The product of a matrix and a column vector; `result` is required. */ +export function multiplyByVector(matrix, cartesian, result) { + if (DEBUG) { + Check.typeOf.object('matrix', matrix); + Check.typeOf.object('cartesian', cartesian); + Check.typeOf.object('result', result); + } + + const vX = cartesian.x; + const vY = cartesian.y; + const vZ = cartesian.z; + + const x = matrix[0] * vX + matrix[3] * vY + matrix[6] * vZ; + const y = matrix[1] * vX + matrix[4] * vY + matrix[7] * vZ; + const z = matrix[2] * vX + matrix[5] * vY + matrix[8] * vZ; + + result.x = x; + result.y = y; + result.z = z; + return result; +} + +/** The column at `index` (0..2) as {x, y, z}; `result` is required. */ +export function getColumn(matrix, index, result) { + if (DEBUG) { + Check.typeOf.object('matrix', matrix); + Check.typeOf.number.greaterThanOrEquals('index', index, 0); + Check.typeOf.number.lessThanOrEquals('index', index, 2); + Check.typeOf.object('result', result); + } + + const startIndex = index * 3; + const x = matrix[startIndex]; + const y = matrix[startIndex + 1]; + const z = matrix[startIndex + 2]; + + result.x = x; + result.y = y; + result.z = z; + return result; +} diff --git a/src/geo/matrix4.js b/src/geo/matrix4.js new file mode 100644 index 000000000..7a98fd8d9 --- /dev/null +++ b/src/geo/matrix4.js @@ -0,0 +1,434 @@ +// Ported from CesiumJS 1.138 (Apache-2.0) — Matrix4.js +// +// 4x4 matrices as plain arrays of 16 numbers in Cesium's column-major index +// order (m[column * 4 + row]), so Cesium accepts them wherever it reads +// m[0]..m[15] (e.g. the Model.modelMatrix setter clones by index). Products +// and getters take a required `result` that may be an input; constructors +// take an optional one. Only the subset gev uses is ported. With DEBUG on +// (development and tests; see debug.js) invalid arguments throw the same +// DeveloperError as Cesium's debug checks. + +import { DEBUG } from './debug.js'; +import { Check } from './check.js'; +import { ZERO } from './cartesian3.js'; + +/** A matrix of zeros (Cesium's `new Matrix4()`). */ +export function create() { + return [ + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, + ]; +} + +export const IDENTITY = Object.freeze([ + 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, + 1.0, +]); + +/** Copies a matrix; `undefined` when there is none. */ +export function clone(matrix, result) { + if (matrix == null) return undefined; + if (result == null) result = create(); + for (let i = 0; i < 16; i++) result[i] = matrix[i]; + return result; +} + +/** Whether two matrices are component-wise equal. */ +export function equals(left, right) { + // Given that most matrices will be transformation matrices, the elements + // are tested in order such that the test is likely to fail as early as possible. + return ( + left === right || + (left != null && + right != null && + // Translation + left[12] === right[12] && + left[13] === right[13] && + left[14] === right[14] && + // Rotation/scale + left[0] === right[0] && + left[1] === right[1] && + left[2] === right[2] && + left[4] === right[4] && + left[5] === right[5] && + left[6] === right[6] && + left[8] === right[8] && + left[9] === right[9] && + left[10] === right[10] && + // Bottom row + left[3] === right[3] && + left[7] === right[7] && + left[11] === right[11] && + left[15] === right[15]) + ); +} + +/** The product of two matrices. */ +export function multiply(left, right, result) { + if (DEBUG) { + Check.typeOf.object('left', left); + Check.typeOf.object('right', right); + Check.typeOf.object('result', result); + } + + const left0 = left[0]; + const left1 = left[1]; + const left2 = left[2]; + const left3 = left[3]; + const left4 = left[4]; + const left5 = left[5]; + const left6 = left[6]; + const left7 = left[7]; + const left8 = left[8]; + const left9 = left[9]; + const left10 = left[10]; + const left11 = left[11]; + const left12 = left[12]; + const left13 = left[13]; + const left14 = left[14]; + const left15 = left[15]; + + const right0 = right[0]; + const right1 = right[1]; + const right2 = right[2]; + const right3 = right[3]; + const right4 = right[4]; + const right5 = right[5]; + const right6 = right[6]; + const right7 = right[7]; + const right8 = right[8]; + const right9 = right[9]; + const right10 = right[10]; + const right11 = right[11]; + const right12 = right[12]; + const right13 = right[13]; + const right14 = right[14]; + const right15 = right[15]; + + const column0Row0 = + left0 * right0 + left4 * right1 + left8 * right2 + left12 * right3; + const column0Row1 = + left1 * right0 + left5 * right1 + left9 * right2 + left13 * right3; + const column0Row2 = + left2 * right0 + left6 * right1 + left10 * right2 + left14 * right3; + const column0Row3 = + left3 * right0 + left7 * right1 + left11 * right2 + left15 * right3; + + const column1Row0 = + left0 * right4 + left4 * right5 + left8 * right6 + left12 * right7; + const column1Row1 = + left1 * right4 + left5 * right5 + left9 * right6 + left13 * right7; + const column1Row2 = + left2 * right4 + left6 * right5 + left10 * right6 + left14 * right7; + const column1Row3 = + left3 * right4 + left7 * right5 + left11 * right6 + left15 * right7; + + const column2Row0 = + left0 * right8 + left4 * right9 + left8 * right10 + left12 * right11; + const column2Row1 = + left1 * right8 + left5 * right9 + left9 * right10 + left13 * right11; + const column2Row2 = + left2 * right8 + left6 * right9 + left10 * right10 + left14 * right11; + const column2Row3 = + left3 * right8 + left7 * right9 + left11 * right10 + left15 * right11; + + const column3Row0 = + left0 * right12 + left4 * right13 + left8 * right14 + left12 * right15; + const column3Row1 = + left1 * right12 + left5 * right13 + left9 * right14 + left13 * right15; + const column3Row2 = + left2 * right12 + left6 * right13 + left10 * right14 + left14 * right15; + const column3Row3 = + left3 * right12 + left7 * right13 + left11 * right14 + left15 * right15; + + result[0] = column0Row0; + result[1] = column0Row1; + result[2] = column0Row2; + result[3] = column0Row3; + result[4] = column1Row0; + result[5] = column1Row1; + result[6] = column1Row2; + result[7] = column1Row3; + result[8] = column2Row0; + result[9] = column2Row1; + result[10] = column2Row2; + result[11] = column2Row3; + result[12] = column3Row0; + result[13] = column3Row1; + result[14] = column3Row2; + result[15] = column3Row3; + return result; +} + +/** The product of two affine transformations (bottom rows assumed [0, 0, 0, 1]). */ +export function multiplyTransformation(left, right, result) { + if (DEBUG) { + Check.typeOf.object('left', left); + Check.typeOf.object('right', right); + Check.typeOf.object('result', result); + } + + const left0 = left[0]; + const left1 = left[1]; + const left2 = left[2]; + const left4 = left[4]; + const left5 = left[5]; + const left6 = left[6]; + const left8 = left[8]; + const left9 = left[9]; + const left10 = left[10]; + const left12 = left[12]; + const left13 = left[13]; + const left14 = left[14]; + + const right0 = right[0]; + const right1 = right[1]; + const right2 = right[2]; + const right4 = right[4]; + const right5 = right[5]; + const right6 = right[6]; + const right8 = right[8]; + const right9 = right[9]; + const right10 = right[10]; + const right12 = right[12]; + const right13 = right[13]; + const right14 = right[14]; + + const column0Row0 = left0 * right0 + left4 * right1 + left8 * right2; + const column0Row1 = left1 * right0 + left5 * right1 + left9 * right2; + const column0Row2 = left2 * right0 + left6 * right1 + left10 * right2; + + const column1Row0 = left0 * right4 + left4 * right5 + left8 * right6; + const column1Row1 = left1 * right4 + left5 * right5 + left9 * right6; + const column1Row2 = left2 * right4 + left6 * right5 + left10 * right6; + + const column2Row0 = left0 * right8 + left4 * right9 + left8 * right10; + const column2Row1 = left1 * right8 + left5 * right9 + left9 * right10; + const column2Row2 = left2 * right8 + left6 * right9 + left10 * right10; + + const column3Row0 = + left0 * right12 + left4 * right13 + left8 * right14 + left12; + const column3Row1 = + left1 * right12 + left5 * right13 + left9 * right14 + left13; + const column3Row2 = + left2 * right12 + left6 * right13 + left10 * right14 + left14; + + result[0] = column0Row0; + result[1] = column0Row1; + result[2] = column0Row2; + result[3] = 0.0; + result[4] = column1Row0; + result[5] = column1Row1; + result[6] = column1Row2; + result[7] = 0.0; + result[8] = column2Row0; + result[9] = column2Row1; + result[10] = column2Row2; + result[11] = 0.0; + result[12] = column3Row0; + result[13] = column3Row1; + result[14] = column3Row2; + result[15] = 1.0; + return result; +} + +/** Transforms a point (w = 1). */ +export function multiplyByPoint(matrix, cartesian, result) { + if (DEBUG) { + Check.typeOf.object('matrix', matrix); + Check.typeOf.object('cartesian', cartesian); + Check.typeOf.object('result', result); + } + + const vX = cartesian.x; + const vY = cartesian.y; + const vZ = cartesian.z; + + const x = matrix[0] * vX + matrix[4] * vY + matrix[8] * vZ + matrix[12]; + const y = matrix[1] * vX + matrix[5] * vY + matrix[9] * vZ + matrix[13]; + const z = matrix[2] * vX + matrix[6] * vY + matrix[10] * vZ + matrix[14]; + + result.x = x; + result.y = y; + result.z = z; + return result; +} + +/** Transforms a direction (w = 0), ignoring the translation. */ +export function multiplyByPointAsVector(matrix, cartesian, result) { + if (DEBUG) { + Check.typeOf.object('matrix', matrix); + Check.typeOf.object('cartesian', cartesian); + Check.typeOf.object('result', result); + } + + const vX = cartesian.x; + const vY = cartesian.y; + const vZ = cartesian.z; + + const x = matrix[0] * vX + matrix[4] * vY + matrix[8] * vZ; + const y = matrix[1] * vX + matrix[5] * vY + matrix[9] * vZ; + const z = matrix[2] * vX + matrix[6] * vY + matrix[10] * vZ; + + result.x = x; + result.y = y; + result.z = z; + return result; +} + +/** The translation part as {x, y, z}. */ +export function getTranslation(matrix, result) { + if (DEBUG) { + Check.typeOf.object('matrix', matrix); + Check.typeOf.object('result', result); + } + + result.x = matrix[12]; + result.y = matrix[13]; + result.z = matrix[14]; + return result; +} + +/** The upper-left 3x3 part. */ +export function getMatrix3(matrix, result) { + if (DEBUG) { + Check.typeOf.object('matrix', matrix); + Check.typeOf.object('result', result); + } + + result[0] = matrix[0]; + result[1] = matrix[1]; + result[2] = matrix[2]; + result[3] = matrix[4]; + result[4] = matrix[5]; + result[5] = matrix[6]; + result[6] = matrix[8]; + result[7] = matrix[9]; + result[8] = matrix[10]; + return result; +} + +/** The column at `index` (0..3) as {x, y, z, w}. */ +export function getColumn(matrix, index, result) { + if (DEBUG) { + Check.typeOf.object('matrix', matrix); + Check.typeOf.number.greaterThanOrEquals('index', index, 0); + Check.typeOf.number.lessThanOrEquals('index', index, 3); + Check.typeOf.object('result', result); + } + + const startIndex = index * 4; + const x = matrix[startIndex]; + const y = matrix[startIndex + 1]; + const z = matrix[startIndex + 2]; + const w = matrix[startIndex + 3]; + + result.x = x; + result.y = y; + result.z = z; + result.w = w; + return result; +} + +/** The inverse of a rotation-plus-translation matrix (no scale or shear). */ +export function inverseTransformation(matrix, result) { + if (DEBUG) { + Check.typeOf.object('matrix', matrix); + Check.typeOf.object('result', result); + } + + // This function is an optimized version of the below 4 lines. + // const rT = Matrix3.transpose(Matrix4.getMatrix3(matrix)); + // const rTN = Matrix3.negate(rT); + // const rTT = Matrix3.multiplyByVector(rTN, Matrix4.getTranslation(matrix)); + // return Matrix4.fromRotationTranslation(rT, rTT, result); + + const matrix0 = matrix[0]; + const matrix1 = matrix[1]; + const matrix2 = matrix[2]; + const matrix4 = matrix[4]; + const matrix5 = matrix[5]; + const matrix6 = matrix[6]; + const matrix8 = matrix[8]; + const matrix9 = matrix[9]; + const matrix10 = matrix[10]; + + const vX = matrix[12]; + const vY = matrix[13]; + const vZ = matrix[14]; + + const x = -matrix0 * vX - matrix1 * vY - matrix2 * vZ; + const y = -matrix4 * vX - matrix5 * vY - matrix6 * vZ; + const z = -matrix8 * vX - matrix9 * vY - matrix10 * vZ; + + result[0] = matrix0; + result[1] = matrix4; + result[2] = matrix8; + result[3] = 0.0; + result[4] = matrix1; + result[5] = matrix5; + result[6] = matrix9; + result[7] = 0.0; + result[8] = matrix2; + result[9] = matrix6; + result[10] = matrix10; + result[11] = 0.0; + result[12] = x; + result[13] = y; + result[14] = z; + result[15] = 1.0; + return result; +} + +/** A transformation from a 3x3 rotation and a translation (zero when absent). */ +export function fromRotationTranslation(rotation, translation, result) { + if (DEBUG) { + Check.typeOf.object('rotation', rotation); + } + + translation = translation ?? ZERO; + if (result == null) result = create(); + result[0] = rotation[0]; + result[1] = rotation[1]; + result[2] = rotation[2]; + result[3] = 0.0; + result[4] = rotation[3]; + result[5] = rotation[4]; + result[6] = rotation[5]; + result[7] = 0.0; + result[8] = rotation[6]; + result[9] = rotation[7]; + result[10] = rotation[8]; + result[11] = 0.0; + result[12] = translation.x; + result[13] = translation.y; + result[14] = translation.z; + result[15] = 1.0; + return result; +} + +/** A uniform scale. */ +export function fromUniformScale(scale, result) { + if (DEBUG) { + Check.typeOf.number('scale', scale); + } + + if (result == null) result = create(); + result[0] = scale; + result[1] = 0.0; + result[2] = 0.0; + result[3] = 0.0; + result[4] = 0.0; + result[5] = scale; + result[6] = 0.0; + result[7] = 0.0; + result[8] = 0.0; + result[9] = 0.0; + result[10] = scale; + result[11] = 0.0; + result[12] = 0.0; + result[13] = 0.0; + result[14] = 0.0; + result[15] = 1.0; + return result; +} diff --git a/src/geo/quaternion.js b/src/geo/quaternion.js new file mode 100644 index 000000000..00234c1c6 --- /dev/null +++ b/src/geo/quaternion.js @@ -0,0 +1,129 @@ +// Ported from CesiumJS 1.138 (Apache-2.0) — Quaternion.js +// +// Rotations as plain-data {x, y, z, w} objects. Only the subset gev uses is ported. +// With DEBUG on (development and tests; see debug.js) invalid arguments throw +// the same DeveloperError as Cesium's debug checks. +// +// Geo values are not Object.prototype objects, so src/app/stateChannel.js +// copyState and src/director/documentFields.js jsonTree reject them: copy +// the fields into a literal before putting one into product state or a +// scene document. + +import { DEBUG } from './debug.js'; +import { Check } from './check.js'; +import { normalize, create as createCartesian3 } from './cartesian3.js'; + +// Values are instances of a module-private constructor rather than object +// literals: every `{x, y, z, w}` literal in the process shares one V8 hidden +// class, and code elsewhere that builds such literals would make geo's hot +// functions polymorphic. The instances still look like plain data: own +// enumerable fields only, no prototype methods or getters. +function Rotation(x, y, z, w) { + this.x = x; + this.y = y; + this.z = z; + this.w = w; +} + +/** A quaternion of zeros by default (Cesium's `new Quaternion(x, y, z, w)`). */ +export function create(x, y, z, w) { + return new Rotation(x ?? 0.0, y ?? 0.0, z ?? 0.0, w ?? 0.0); +} + +let fromAxisAngleScratch = createCartesian3(); + +/** The rotation by `angle` radians around `axis` (normalised here; with DEBUG on, a zero axis throws). */ +export function fromAxisAngle(axis, angle, result) { + if (DEBUG) { + Check.typeOf.object('axis', axis); + Check.typeOf.number('angle', angle); + } + + const halfAngle = angle / 2.0; + const s = Math.sin(halfAngle); + fromAxisAngleScratch = normalize(axis, fromAxisAngleScratch); + + const x = fromAxisAngleScratch.x * s; + const y = fromAxisAngleScratch.y * s; + const z = fromAxisAngleScratch.z * s; + const w = Math.cos(halfAngle); + if (result == null) return create(x, y, z, w); + result.x = x; + result.y = y; + result.z = z; + result.w = w; + return result; +} + +const fromRotationMatrixNext = [1, 2, 0]; +const fromRotationMatrixQuat = new Array(3); + +// Matrix3.getElementIndex(column, row) +const elementIndex = (column, row) => column * 3 + row; + +/** The rotation of a 3x3 rotation matrix. */ +export function fromRotationMatrix(matrix, result) { + if (DEBUG) { + Check.typeOf.object('matrix', matrix); + } + + let root; + let x; + let y; + let z; + let w; + + const m00 = matrix[0]; + const m11 = matrix[4]; + const m22 = matrix[8]; + const trace = m00 + m11 + m22; + + if (trace > 0.0) { + // |w| > 1/2, may as well choose w > 1/2 + root = Math.sqrt(trace + 1.0); // 2w + w = 0.5 * root; + root = 0.5 / root; // 1/(4w) + + x = (matrix[5] - matrix[7]) * root; + y = (matrix[6] - matrix[2]) * root; + z = (matrix[1] - matrix[3]) * root; + } else { + // |w| <= 1/2 + const next = fromRotationMatrixNext; + + let i = 0; + if (m11 > m00) { + i = 1; + } + if (m22 > m00 && m22 > m11) { + i = 2; + } + const j = next[i]; + const k = next[j]; + + root = Math.sqrt( + matrix[elementIndex(i, i)] - + matrix[elementIndex(j, j)] - + matrix[elementIndex(k, k)] + + 1.0, + ); + + const quat = fromRotationMatrixQuat; + quat[i] = 0.5 * root; + root = 0.5 / root; + w = (matrix[elementIndex(k, j)] - matrix[elementIndex(j, k)]) * root; + quat[j] = (matrix[elementIndex(j, i)] + matrix[elementIndex(i, j)]) * root; + quat[k] = (matrix[elementIndex(k, i)] + matrix[elementIndex(i, k)]) * root; + + x = -quat[0]; + y = -quat[1]; + z = -quat[2]; + } + + if (result == null) return create(x, y, z, w); + result.x = x; + result.y = y; + result.z = z; + result.w = w; + return result; +} diff --git a/src/geo/releaseParity.test.mjs b/src/geo/releaseParity.test.mjs new file mode 100644 index 000000000..21faa39e8 --- /dev/null +++ b/src/geo/releaseParity.test.mjs @@ -0,0 +1,446 @@ +// Release parity: with DEBUG off, src/geo returns what release Cesium returns. +// +// Production builds fold geo's DEBUG to false (see debug.js) and ship Cesium's +// Build/Cesium, whose debug checks are stripped. So for a degenerate input that +// throws in development (a zero vector to normalize, a NaN origin, a divisor +// of 0, an index out of range) both must now return the same values, NaN for +// NaN and -0 for -0, or fail with the same error. This file loads a copy of +// src/geo whose debug.js says `DEBUG = false`, and the stripped Cesium build. +import { test, after } from 'node:test'; +import assert from 'node:assert/strict'; +import { createRequire } from 'node:module'; +import { + cpSync, + mkdtempSync, + readdirSync, + rmSync, + writeFileSync, +} from 'node:fs'; +import { tmpdir } from 'node:os'; +import path from 'node:path'; +import { fileURLToPath, pathToFileURL } from 'node:url'; + +const here = path.dirname(fileURLToPath(import.meta.url)); +const copy = mkdtempSync(path.join(tmpdir(), 'geo-release-')); +for (const name of readdirSync(here)) { + if (name.endsWith('.js')) + cpSync(path.join(here, name), path.join(copy, name)); +} +writeFileSync(path.join(copy, 'debug.js'), 'export const DEBUG = false;\n'); +after(() => rmSync(copy, { recursive: true, force: true })); + +const geo = await import(pathToFileURL(path.join(copy, 'index.js')).href); +const { DEBUG } = await import(pathToFileURL(path.join(copy, 'debug.js')).href); +const Cesium = createRequire(import.meta.url)('cesium/Build/Cesium/index.cjs'); + +const W = Cesium.Ellipsoid.WGS84; +const C3 = Cesium.Cartesian3; + +// Builders run once per side so no object is shared between geo and Cesium. +const vec = (x, y, z) => (side) => + side === 'geo' ? geo.cartesian3.create(x, y, z) : new C3(x, y, z); +const carto = (lon, lat, h) => (side) => + side === 'geo' + ? geo.cartographic.create(lon, lat, h) + : new Cesium.Cartographic(lon, lat, h); +const mat3 = (side) => + side === 'geo' + ? geo.matrix3.fromRotationZ(0.4) + : Cesium.Matrix3.fromRotationZ(0.4, new Cesium.Matrix3()); +const mat4 = (side) => + side === 'geo' ? geo.matrix4.create() : new Cesium.Matrix4(); +const nanMat4 = (side) => { + const m = side === 'geo' ? geo.matrix4.create() : new Cesium.Matrix4(); + for (let i = 0; i < 16; i++) + m[i] = i === 5 ? Number.NaN : i % 5 === 0 ? 1 : 0; + return m; +}; +const vec4 = (side) => + side === 'geo' ? { x: 0, y: 0, z: 0, w: 0 } : new Cesium.Cartesian4(); +const quat = (side) => + side === 'geo' ? geo.quaternion.create() : new Cesium.Quaternion(); +const hpr = (h, p, r) => (side) => + side === 'geo' + ? { heading: h, pitch: p, roll: r } + : new Cesium.HeadingPitchRoll(h, p, r); +const ellipsoid = (side) => (side === 'geo' ? geo.ellipsoid.WGS84 : W); +const value = (v) => () => v; +const list = + (...v) => + () => [...v]; + +const ZERO = vec(0, 0, 0); +const NAN = vec(1, Number.NaN, 2); +const INF = vec(Number.POSITIVE_INFINITY, 0, 0); +const R = vec(0, 0, 0); +const NaNum = value(Number.NaN); + +// [name, geo function, Cesium function, argument builders] +const CASES = [ + ['math.toRadians()', geo.math.toRadians, Cesium.Math.toRadians, []], + [ + 'math.toDegrees(null)', + geo.math.toDegrees, + Cesium.Math.toDegrees, + [value(null)], + ], + [ + 'math.clamp(str)', + geo.math.clamp, + Cesium.Math.clamp, + [value('5'), value(0), value(2)], + ], + ['math.mod(5, 0)', geo.math.mod, Cesium.Math.mod, [value(5), value(0)]], + ['math.mod(-5, 0)', geo.math.mod, Cesium.Math.mod, [value(-5), value(0)]], + [ + 'math.mod(undefined, 3)', + geo.math.mod, + Cesium.Math.mod, + [value(undefined), value(3)], + ], + ['math.zeroToTwoPi()', geo.math.zeroToTwoPi, Cesium.Math.zeroToTwoPi, []], + [ + 'math.negativePiToPi(NaN)', + geo.math.negativePiToPi, + Cesium.Math.negativePiToPi, + [NaNum], + ], + [ + 'math.equalsEpsilon(null, 1)', + geo.math.equalsEpsilon, + Cesium.Math.equalsEpsilon, + [value(null), value(1), value(0)], + ], + + [ + 'cartesian3.normalize(zero)', + geo.cartesian3.normalize, + C3.normalize, + [ZERO, R], + ], + [ + 'cartesian3.normalize(NaN)', + geo.cartesian3.normalize, + C3.normalize, + [NAN, R], + ], + [ + 'cartesian3.normalize(Infinity)', + geo.cartesian3.normalize, + C3.normalize, + [INF, R], + ], + ['cartesian3.magnitude(NaN)', geo.cartesian3.magnitude, C3.magnitude, [NAN]], + [ + 'cartesian3.multiplyByScalar(str)', + geo.cartesian3.multiplyByScalar, + C3.multiplyByScalar, + [vec(1, 2, 3), value('2'), R], + ], + [ + 'cartesian3.divideByScalar(0)', + geo.cartesian3.divideByScalar, + C3.divideByScalar, + [vec(1, -2, 0), value(0), R], + ], + [ + 'cartesian3.lerp(NaN t)', + geo.cartesian3.lerp, + C3.lerp, + [vec(1, 2, 3), ZERO, NaNum, R], + ], + [ + 'cartesian3.fromRadians(NaN)', + geo.cartesian3.fromRadians, + C3.fromRadians, + [NaNum, value(0.2), value(10)], + ], + [ + 'cartesian3.fromRadians(str lat)', + geo.cartesian3.fromRadians, + C3.fromRadians, + [value(0.1), value('0.2')], + ], + [ + 'cartesian3.fromDegrees(NaN, NaN)', + geo.cartesian3.fromDegrees, + C3.fromDegrees, + [NaNum, NaNum], + ], + [ + 'cartesian3.fromDegrees(str height)', + geo.cartesian3.fromDegrees, + C3.fromDegrees, + [value(10), value(20), value('5')], + ], + [ + 'cartesian3.fromDegreesArray([])', + geo.cartesian3.fromDegreesArray, + C3.fromDegreesArray, + [list()], + ], + [ + 'cartesian3.fromDegreesArray(odd)', + geo.cartesian3.fromDegreesArray, + C3.fromDegreesArray, + [list(1, 2, 3)], + ], + [ + 'cartesian3.fromDegreesArrayHeights(4)', + geo.cartesian3.fromDegreesArrayHeights, + C3.fromDegreesArrayHeights, + [list(1, 2, 3, 4)], + ], + [ + 'cartesian3.fromDegreesArrayHeights(str)', + geo.cartesian3.fromDegreesArrayHeights, + C3.fromDegreesArrayHeights, + [list(1, 2, 'h')], + ], + + [ + 'cartographic.fromRadians(NaN)', + geo.cartographic.fromRadians, + Cesium.Cartographic.fromRadians, + [NaNum, value(0.1)], + ], + [ + 'cartographic.fromDegrees(str)', + geo.cartographic.fromDegrees, + Cesium.Cartographic.fromDegrees, + [value('10'), value(20)], + ], + [ + 'cartographic.fromCartesian(NaN)', + geo.cartographic.fromCartesian, + Cesium.Cartographic.fromCartesian, + [NAN], + ], + [ + 'cartographic.fromCartesian(zero)', + geo.cartographic.fromCartesian, + Cesium.Cartographic.fromCartesian, + [ZERO], + ], + [ + 'cartographic.fromCartesian(Infinity)', + geo.cartographic.fromCartesian, + Cesium.Cartographic.fromCartesian, + [INF], + ], + [ + 'cartographic.toCartesian(NaN)', + geo.cartographic.toCartesian, + Cesium.Cartographic.toCartesian, + [carto(Number.NaN, 0, 0)], + ], + + [ + 'ellipsoid.create(-1, 2, 3)', + geo.ellipsoid.create, + (x, y, z) => new Cesium.Ellipsoid(x, y, z), + [value(-1), value(2), value(3)], + ], + [ + 'WGS84.geodeticSurfaceNormal(NaN)', + geo.ellipsoid.WGS84.geodeticSurfaceNormal, + (c, r) => W.geodeticSurfaceNormal(c, r), + [NAN, R], + ], + [ + 'WGS84.geodeticSurfaceNormal(Infinity)', + geo.ellipsoid.WGS84.geodeticSurfaceNormal, + (c, r) => W.geodeticSurfaceNormal(c, r), + [INF, R], + ], + [ + 'WGS84.geodeticSurfaceNormalCartographic(NaN)', + geo.ellipsoid.WGS84.geodeticSurfaceNormalCartographic, + (c, r) => W.geodeticSurfaceNormalCartographic(c, r), + [carto(0.1, Number.NaN, 0), R], + ], + [ + 'WGS84.scaleToGeodeticSurface(NaN)', + geo.ellipsoid.WGS84.scaleToGeodeticSurface, + (c, r) => W.scaleToGeodeticSurface(c, r), + [NAN, R], + ], + [ + 'WGS84.cartographicToCartesian(NaN)', + geo.ellipsoid.WGS84.cartographicToCartesian, + (c, r) => W.cartographicToCartesian(c, r), + [carto(Number.NaN, 0.1, 5), R], + ], + [ + 'WGS84.cartesianToCartographic(NaN)', + geo.ellipsoid.WGS84.cartesianToCartographic, + (c, r) => W.cartesianToCartographic(c, r), + [NAN, carto(0, 0, 0)], + ], + + [ + 'matrix3.fromRotationZ(NaN)', + geo.matrix3.fromRotationZ, + Cesium.Matrix3.fromRotationZ, + [NaNum, mat3], + ], + [ + 'matrix3.getColumn(5)', + geo.matrix3.getColumn, + Cesium.Matrix3.getColumn, + [mat3, value(5), R], + ], + [ + 'matrix3.getColumn(-1)', + geo.matrix3.getColumn, + Cesium.Matrix3.getColumn, + [mat3, value(-1), R], + ], + [ + 'matrix4.getColumn(4)', + geo.matrix4.getColumn, + Cesium.Matrix4.getColumn, + [nanMat4, value(4), vec4], + ], + [ + 'matrix4.inverseTransformation(NaN)', + geo.matrix4.inverseTransformation, + Cesium.Matrix4.inverseTransformation, + [nanMat4, mat4], + ], + [ + 'matrix4.multiplyByPoint(NaN)', + geo.matrix4.multiplyByPoint, + Cesium.Matrix4.multiplyByPoint, + [nanMat4, NAN, R], + ], + [ + 'matrix4.fromUniformScale(str)', + geo.matrix4.fromUniformScale, + Cesium.Matrix4.fromUniformScale, + [value('2'), mat4], + ], + + [ + 'quaternion.fromAxisAngle(zero)', + geo.quaternion.fromAxisAngle, + Cesium.Quaternion.fromAxisAngle, + [ZERO, value(0.5), quat], + ], + [ + 'quaternion.fromAxisAngle(NaN angle)', + geo.quaternion.fromAxisAngle, + Cesium.Quaternion.fromAxisAngle, + [vec(0, 0, 1), NaNum, quat], + ], + + [ + 'transforms.eastNorthUpToFixedFrame(NaN)', + geo.transforms.eastNorthUpToFixedFrame, + Cesium.Transforms.eastNorthUpToFixedFrame, + [NAN], + ], + [ + 'transforms.eastNorthUpToFixedFrame(Infinity)', + geo.transforms.eastNorthUpToFixedFrame, + Cesium.Transforms.eastNorthUpToFixedFrame, + [INF], + ], + [ + 'transforms.eastNorthUpToFixedFrame(zero)', + geo.transforms.eastNorthUpToFixedFrame, + Cesium.Transforms.eastNorthUpToFixedFrame, + [ZERO], + ], + [ + 'transforms.headingPitchRollToFixedFrame(NaN origin)', + geo.transforms.headingPitchRollToFixedFrame, + Cesium.Transforms.headingPitchRollToFixedFrame, + [NAN, hpr(0.3, -0.2, 0.1)], + ], + [ + 'transforms.headingPitchRollToFixedFrame(NaN heading)', + geo.transforms.headingPitchRollToFixedFrame, + Cesium.Transforms.headingPitchRollToFixedFrame, + [vec(6378137, 1, 2), hpr(Number.NaN, 0, 0)], + ], + [ + 'transforms.headingPitchRollToFixedFrame(partial)', + geo.transforms.headingPitchRollToFixedFrame, + Cesium.Transforms.headingPitchRollToFixedFrame, + [vec(6378137, 1, 2), (side) => ({ heading: 0.3 })], + ], +]; + +const FIELDS = ['x', 'y', 'z', 'w', 'longitude', 'latitude', 'height']; + +/** The numbers a result holds, in a fixed order, so geo and Cesium values compare. */ +function numbers(result) { + if (result === undefined || result === null || typeof result !== 'object') { + return [result]; + } + if (Array.isArray(result) || 0 in result) { + const out = []; + for (let i = 0; i in result; i++) out.push(numbers(result[i])); + return out.flat(); + } + if ('radii' in result) { + return [ + 'radii', + 'radiiSquared', + 'oneOverRadii', + 'oneOverRadiiSquared', + ].flatMap((key) => numbers(result[key])); + } + return FIELDS.filter((key) => key in result).map((key) => result[key]); +} + +function outcome(fn, args) { + try { + return { values: numbers(fn(...args)) }; + } catch (error) { + return { error: `${error.name}: ${error.message}` }; + } +} + +test('the copy under test has its debug checks off', () => { + assert.equal(DEBUG, false); + assert.ok(!('normalized result is not a number' in Cesium)); +}); + +for (const [name, geoFn, cesiumFn, builders] of CASES) { + test(`release ${name} matches Build/Cesium`, () => { + const g = outcome( + geoFn, + builders.map((build) => build('geo')), + ); + const c = outcome( + cesiumFn, + builders.map((build) => build('cesium')), + ); + if (c.error !== undefined || g.error !== undefined) { + assert.equal(g.error, c.error); + return; + } + assert.equal(g.values.length, c.values.length); + g.values.forEach((v, i) => + assert.ok( + Object.is(v, c.values[i]), + `${name}[${i}]: geo ${v}, Cesium ${c.values[i]}`, + ), + ); + }); +} + +test('release geo returns NaN where development geo throws', () => { + const r = geo.cartesian3.normalize( + geo.cartesian3.create(), + geo.cartesian3.create(), + ); + assert.ok(Number.isNaN(r.x) && Number.isNaN(r.y) && Number.isNaN(r.z)); + const m = geo.transforms.eastNorthUpToFixedFrame( + geo.cartesian3.create(Number.NaN, 0, 0), + ); + assert.ok(Number.isNaN(m[0])); + assert.ok(Number.isNaN(geo.math.mod(5, 0))); +}); diff --git a/src/geo/scaleToGeodeticSurface.js b/src/geo/scaleToGeodeticSurface.js new file mode 100644 index 000000000..feeed0022 --- /dev/null +++ b/src/geo/scaleToGeodeticSurface.js @@ -0,0 +1,142 @@ +// Ported from CesiumJS 1.138 (Apache-2.0) — scaleToGeodeticSurface.js +// +// Scales a position onto an ellipsoid's geodetic surface (Newton iteration). +// Internal to src/geo (not re-exported by index.js): used by ellipsoid.js and, +// as in Cesium, directly by cartographic.fromCartesian. + +import { DEBUG } from './debug.js'; +import { DeveloperError } from './developerError.js'; +import { EPSILON12 } from './math.js'; +import { + clone, + create as vec, + magnitude, + multiplyByScalar, +} from './cartesian3.js'; + +const scaleIntersectionScratch = vec(); +const scaleGradientScratch = vec(); + +/** + * Scales a position to the geodetic surface of the ellipsoid given by + * `oneOverRadii`/`oneOverRadiiSquared`, or `undefined` near the centre. + */ +export function scaleToGeodeticSurface( + cartesian, + oneOverRadii, + oneOverRadiiSquared, + centerToleranceSquared, + result, +) { + if (DEBUG) { + if (cartesian == null) { + throw new DeveloperError('cartesian is required.'); + } + if (oneOverRadii == null) { + throw new DeveloperError('oneOverRadii is required.'); + } + if (oneOverRadiiSquared == null) { + throw new DeveloperError('oneOverRadiiSquared is required.'); + } + if (centerToleranceSquared == null) { + throw new DeveloperError('centerToleranceSquared is required.'); + } + } + + const positionX = cartesian.x; + const positionY = cartesian.y; + const positionZ = cartesian.z; + + const oneOverRadiiX = oneOverRadii.x; + const oneOverRadiiY = oneOverRadii.y; + const oneOverRadiiZ = oneOverRadii.z; + + const x2 = positionX * positionX * oneOverRadiiX * oneOverRadiiX; + const y2 = positionY * positionY * oneOverRadiiY * oneOverRadiiY; + const z2 = positionZ * positionZ * oneOverRadiiZ * oneOverRadiiZ; + + // Compute the squared ellipsoid norm. + const squaredNorm = x2 + y2 + z2; + const ratio = Math.sqrt(1.0 / squaredNorm); + + // As an initial approximation, assume that the radial intersection is the projection point. + const intersection = multiplyByScalar( + cartesian, + ratio, + scaleIntersectionScratch, + ); + + // If the position is near the center, the iteration will not converge. + if (squaredNorm < centerToleranceSquared) { + return !isFinite(ratio) ? undefined : clone(intersection, result); + } + + const oneOverRadiiSquaredX = oneOverRadiiSquared.x; + const oneOverRadiiSquaredY = oneOverRadiiSquared.y; + const oneOverRadiiSquaredZ = oneOverRadiiSquared.z; + + // Use the gradient at the intersection point in place of the true unit normal. + // The difference in magnitude will be absorbed in the multiplier. + const gradient = scaleGradientScratch; + gradient.x = intersection.x * oneOverRadiiSquaredX * 2.0; + gradient.y = intersection.y * oneOverRadiiSquaredY * 2.0; + gradient.z = intersection.z * oneOverRadiiSquaredZ * 2.0; + + // Compute the initial guess at the normal vector multiplier, lambda. + let lambda = + ((1.0 - ratio) * magnitude(cartesian)) / (0.5 * magnitude(gradient)); + let correction = 0.0; + + let func; + let denominator; + let xMultiplier; + let yMultiplier; + let zMultiplier; + let xMultiplier2; + let yMultiplier2; + let zMultiplier2; + let xMultiplier3; + let yMultiplier3; + let zMultiplier3; + + do { + lambda -= correction; + + xMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquaredX); + yMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquaredY); + zMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquaredZ); + + xMultiplier2 = xMultiplier * xMultiplier; + yMultiplier2 = yMultiplier * yMultiplier; + zMultiplier2 = zMultiplier * zMultiplier; + + xMultiplier3 = xMultiplier2 * xMultiplier; + yMultiplier3 = yMultiplier2 * yMultiplier; + zMultiplier3 = zMultiplier2 * zMultiplier; + + func = x2 * xMultiplier2 + y2 * yMultiplier2 + z2 * zMultiplier2 - 1.0; + + // "denominator" here refers to the use of this expression in the velocity and acceleration + // computations in the sections to follow. + denominator = + x2 * xMultiplier3 * oneOverRadiiSquaredX + + y2 * yMultiplier3 * oneOverRadiiSquaredY + + z2 * zMultiplier3 * oneOverRadiiSquaredZ; + + const derivative = -2.0 * denominator; + + correction = func / derivative; + } while (Math.abs(func) > EPSILON12); + + if (result == null) { + return vec( + positionX * xMultiplier, + positionY * yMultiplier, + positionZ * zMultiplier, + ); + } + result.x = positionX * xMultiplier; + result.y = positionY * yMultiplier; + result.z = positionZ * zMultiplier; + return result; +} diff --git a/src/geo/transforms.js b/src/geo/transforms.js new file mode 100644 index 000000000..76101a186 --- /dev/null +++ b/src/geo/transforms.js @@ -0,0 +1,240 @@ +// Ported from CesiumJS 1.138 (Apache-2.0) — Transforms.js, Quaternion.js, Matrix4.js +// +// Local east-north-up frames and heading-pitch-roll orientations as 4x4 +// column-major matrices (see matrix4.js). Outputs match Cesium bit for bit, +// including the degenerate frames at the Earth's centre and at the poles. +// With DEBUG on (development and tests; see debug.js) invalid arguments throw +// the same DeveloperError as Cesium's debug checks. + +import { DEBUG } from './debug.js'; +import { Check } from './check.js'; +import { DeveloperError } from './developerError.js'; +import { EPSILON14, equalsEpsilon as scalarEqualsEpsilon } from './math.js'; +import { + create as createCartesian3, + cross, + equalsEpsilon, + multiplyByScalar, + normalize, + UNIT_X, + UNIT_Y, + UNIT_Z, + ZERO, +} from './cartesian3.js'; +import { WGS84 } from './ellipsoid.js'; +import { create as createMatrix4, multiply } from './matrix4.js'; +import { create as createQuaternion, fromAxisAngle } from './quaternion.js'; + +// Cesium's degeneratePositionLocalFrame for the three ENU axes. +const DEGENERATE_EAST = [0, 1, 0]; +const DEGENERATE_NORTH = [-1, 0, 0]; +const DEGENERATE_UP = [0, 0, 1]; + +const scratchEast = createCartesian3(); +const scratchNorth = createCartesian3(); +const scratchUp = createCartesian3(); + +function unpack(array, result) { + result.x = array[0]; + result.y = array[1]; + result.z = array[2]; + return result; +} + +/** + * The transformation from the local east-north-up frame centred at `origin` + * to the Earth-fixed frame (Cesium's + * `Transforms.localFrameToFixedFrameGenerator('east', 'north')`). + */ +export function eastNorthUpToFixedFrame(origin, ellipsoid, result) { + if (DEBUG) { + if (origin == null) { + throw new DeveloperError('origin is required.'); + } + if (isNaN(origin.x) || isNaN(origin.y) || isNaN(origin.z)) { + throw new DeveloperError('origin has a NaN component'); + } + } + + if (result == null) result = createMatrix4(); + + const east = scratchEast; + const north = scratchNorth; + const up = scratchUp; + + if (equalsEpsilon(origin, ZERO, EPSILON14)) { + // If x, y, and z are zero, use the degenerate local frame, which is a special case + unpack(DEGENERATE_EAST, east); + unpack(DEGENERATE_NORTH, north); + unpack(DEGENERATE_UP, up); + } else if ( + scalarEqualsEpsilon(origin.x, 0.0, EPSILON14) && + scalarEqualsEpsilon(origin.y, 0.0, EPSILON14) + ) { + // If x and y are zero, assume origin is at a pole, which is a special case. + const sign = Math.sign(origin.z); + unpack(DEGENERATE_EAST, east); + unpack(DEGENERATE_NORTH, north); + multiplyByScalar(north, sign, north); + unpack(DEGENERATE_UP, up); + multiplyByScalar(up, sign, up); + } else { + // Cesium falls back to Ellipsoid.default; gev never changes it, so WGS84. + ellipsoid = ellipsoid ?? WGS84; + ellipsoid.geodeticSurfaceNormal(origin, up); + east.x = -origin.y; + east.y = origin.x; + east.z = 0.0; + normalize(east, east); + cross(up, east, north); + } + + result[0] = east.x; + result[1] = east.y; + result[2] = east.z; + result[3] = 0.0; + result[4] = north.x; + result[5] = north.y; + result[6] = north.z; + result[7] = 0.0; + result[8] = up.x; + result[9] = up.y; + result[10] = up.z; + result[11] = 0.0; + result[12] = origin.x; + result[13] = origin.y; + result[14] = origin.z; + result[15] = 1.0; + return result; +} + +// Quaternion.multiply +function multiplyQuaternions(left, right, result) { + const leftX = left.x; + const leftY = left.y; + const leftZ = left.z; + const leftW = left.w; + + const rightX = right.x; + const rightY = right.y; + const rightZ = right.z; + const rightW = right.w; + + const x = leftW * rightX + leftX * rightW + leftY * rightZ - leftZ * rightY; + const y = leftW * rightY - leftX * rightZ + leftY * rightW + leftZ * rightX; + const z = leftW * rightZ + leftX * rightY - leftY * rightX + leftZ * rightW; + const w = leftW * rightW - leftX * rightX - leftY * rightY - leftZ * rightZ; + + result.x = x; + result.y = y; + result.z = z; + result.w = w; + return result; +} + +const scratchHPRAxisQuaternion = createQuaternion(); + +// Quaternion.fromHeadingPitchRoll +function quaternionFromHeadingPitchRoll(headingPitchRoll, result) { + const rollQuaternion = fromAxisAngle( + UNIT_X, + headingPitchRoll.roll, + scratchHPRAxisQuaternion, + ); + const pitchQuaternion = fromAxisAngle( + UNIT_Y, + -headingPitchRoll.pitch, + result, + ); + result = multiplyQuaternions( + pitchQuaternion, + rollQuaternion, + pitchQuaternion, + ); + const headingQuaternion = fromAxisAngle( + UNIT_Z, + -headingPitchRoll.heading, + scratchHPRAxisQuaternion, + ); + return multiplyQuaternions(headingQuaternion, result, result); +} + +// Matrix4.fromTranslationQuaternionRotationScale with a zero translation and a unit scale. +function rotationFromQuaternion(rotation, result) { + const scaleX = 1.0; + const scaleY = 1.0; + const scaleZ = 1.0; + + const x2 = rotation.x * rotation.x; + const xy = rotation.x * rotation.y; + const xz = rotation.x * rotation.z; + const xw = rotation.x * rotation.w; + const y2 = rotation.y * rotation.y; + const yz = rotation.y * rotation.z; + const yw = rotation.y * rotation.w; + const z2 = rotation.z * rotation.z; + const zw = rotation.z * rotation.w; + const w2 = rotation.w * rotation.w; + + const m00 = x2 - y2 - z2 + w2; + const m01 = 2.0 * (xy - zw); + const m02 = 2.0 * (xz + yw); + + const m10 = 2.0 * (xy + zw); + const m11 = -x2 + y2 - z2 + w2; + const m12 = 2.0 * (yz - xw); + + const m20 = 2.0 * (xz - yw); + const m21 = 2.0 * (yz + xw); + const m22 = -x2 - y2 + z2 + w2; + + result[0] = m00 * scaleX; + result[1] = m10 * scaleX; + result[2] = m20 * scaleX; + result[3] = 0.0; + result[4] = m01 * scaleY; + result[5] = m11 * scaleY; + result[6] = m21 * scaleY; + result[7] = 0.0; + result[8] = m02 * scaleZ; + result[9] = m12 * scaleZ; + result[10] = m22 * scaleZ; + result[11] = 0.0; + result[12] = 0.0; + result[13] = 0.0; + result[14] = 0.0; + result[15] = 1.0; + return result; +} + +const scratchHPRQuaternion = createQuaternion(); +const scratchHPRMatrix4 = createMatrix4(); + +/** + * The transformation from a frame rotated by `headingPitchRoll` ({heading, + * pitch, roll} in radians) inside the local frame centred at `origin` to the + * Earth-fixed frame. Same argument order as Cesium: `fixedFrameTransform` + * defaults to the east-north-up frame and `result` is the fifth argument. + */ +export function headingPitchRollToFixedFrame( + origin, + headingPitchRoll, + ellipsoid, + fixedFrameTransform, + result, +) { + // The private helpers below skip Cesium's checks on their internal + // arguments (scratch objects and constants), which can never fail. + if (DEBUG) { + Check.typeOf.object('HeadingPitchRoll', headingPitchRoll); + } + + fixedFrameTransform = fixedFrameTransform ?? eastNorthUpToFixedFrame; + const hprQuaternion = quaternionFromHeadingPitchRoll( + headingPitchRoll, + scratchHPRQuaternion, + ); + const hprMatrix = rotationFromQuaternion(hprQuaternion, scratchHPRMatrix4); + result = fixedFrameTransform(origin, ellipsoid, result); + return multiply(result, hprMatrix, result); +} diff --git a/src/geo/transforms.test.mjs b/src/geo/transforms.test.mjs new file mode 100644 index 000000000..2811d48c5 --- /dev/null +++ b/src/geo/transforms.test.mjs @@ -0,0 +1,152 @@ +// src/geo/transforms.test.mjs +import { test } from 'node:test'; +import assert from 'node:assert/strict'; +import * as Cesium from 'cesium'; +import * as transforms from './transforms.js'; + +const close = (a, b, tol, label) => + assert.ok(Math.abs(a - b) <= tol, `${label}: ${a} vs ${b}`); +const ORIGINS = [ + [-73.9857, 40.7484, 381], + [0, 0, 0], + [179.9, -45, 10], + [10, 89.9999999, 0], + [10, -90, 0], + [0, 90, 0], +]; + +test('eastNorthUpToFixedFrame equals Cesium, including the poles', () => { + for (const [lon, lat, h] of ORIGINS) { + const origin = Cesium.Cartesian3.fromDegrees(lon, lat, h); + const ours = transforms.eastNorthUpToFixedFrame({ + x: origin.x, + y: origin.y, + z: origin.z, + }); + const theirs = Cesium.Transforms.eastNorthUpToFixedFrame(origin); + for (let i = 0; i < 16; i++) + close( + ours[i], + theirs[i], + 1e-9 * Math.max(1, Math.abs(theirs[i])), + `${lon},${lat} [${i}]`, + ); + } +}); + +test('headingPitchRollToFixedFrame equals Cesium', () => { + for (const [lon, lat, h] of ORIGINS.slice(0, 3)) { + const origin = Cesium.Cartesian3.fromDegrees(lon, lat, h); + for (const [hd, p, r] of [ + [0, 0, 0], + [0.7, -0.2, 0.1], + [-2.5, 1.2, -0.4], + ]) { + const ours = transforms.headingPitchRollToFixedFrame( + { x: origin.x, y: origin.y, z: origin.z }, + { heading: hd, pitch: p, roll: r }, + ); + const theirs = Cesium.Transforms.headingPitchRollToFixedFrame( + origin, + new Cesium.HeadingPitchRoll(hd, p, r), + ); + for (let i = 0; i < 16; i++) + close( + ours[i], + theirs[i], + 1e-9 * Math.max(1, Math.abs(theirs[i])), + `hpr [${i}]`, + ); + } + } +}); + +test('headingPitchRollToFixedFrame keeps Cesium argument order (result is the fifth argument)', () => { + const origin = Cesium.Cartesian3.fromDegrees(-73.9857, 40.7484, 381); + const hpr = new Cesium.HeadingPitchRoll(0.7, -0.2, 0.1); + const out = Array.from({ length: 16 }, () => 0); + const ours = transforms.headingPitchRollToFixedFrame( + { x: origin.x, y: origin.y, z: origin.z }, + hpr, + Cesium.Ellipsoid.WGS84, + undefined, + out, + ); + assert.equal(ours, out); + const theirs = Cesium.Transforms.headingPitchRollToFixedFrame( + origin, + hpr, + Cesium.Ellipsoid.WGS84, + undefined, + new Cesium.Matrix4(), + ); + assert.deepEqual( + Array.from(ours), + Array.from({ length: 16 }, (_, i) => theirs[i]), + ); + let calledWith; + const frame = (o, e, r) => { + calledWith = [o, e, r]; + return transforms.eastNorthUpToFixedFrame(o, e, r); + }; + transforms.headingPitchRollToFixedFrame( + { x: origin.x, y: origin.y, z: origin.z }, + hpr, + undefined, + frame, + out, + ); + assert.equal(calledWith[2], out); +}); + +// WGS84 polar radius: {0, 0, ±b} has x = y = 0 exactly, so it takes Cesium's pole branch. +const POLAR_RADIUS = 6356752.3142451793; +const sameMatrix = (ours, theirs, label) => { + for (let i = 0; i < 16; i++) { + assert.ok( + Object.is(ours[i], theirs[i]), + `${label} [${i}]: ${ours[i]} vs ${theirs[i]}`, + ); + } +}; + +test('eastNorthUpToFixedFrame equals Cesium bit for bit at the centre and the exact poles', () => { + for (const [x, y, z, branch] of [ + [0, 0, 0, 'centre'], + [-0, 0, -0, 'centre (-0)'], + [1e-15, -1e-15, 1e-15, 'within EPSILON14 of the centre'], + [0, 0, POLAR_RADIUS, 'north pole'], + [0, 0, -POLAR_RADIUS, 'south pole'], + [1e-15, -0, -5, 'south pole within EPSILON14'], + ]) { + const ours = transforms.eastNorthUpToFixedFrame({ x, y, z }); + const theirs = Cesium.Transforms.eastNorthUpToFixedFrame( + new Cesium.Cartesian3(x, y, z), + ); + sameMatrix(ours, theirs, branch); + const hpr = new Cesium.HeadingPitchRoll(0.3, -0.2, 0.1); + sameMatrix( + transforms.headingPitchRollToFixedFrame({ x, y, z }, hpr), + Cesium.Transforms.headingPitchRollToFixedFrame( + new Cesium.Cartesian3(x, y, z), + hpr, + ), + `hpr ${branch}`, + ); + } + // The centre uses the fixed degenerate frame; the pole frame flips north and up with the sign of z. + const centre = transforms.eastNorthUpToFixedFrame({ x: 0, y: 0, z: 0 }); + assert.deepEqual( + Array.from(centre.slice(0, 11)), + [0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 1], + ); + const south = transforms.eastNorthUpToFixedFrame({ + x: 0, + y: 0, + z: -POLAR_RADIUS, + }); + assert.deepEqual( + Array.from(south.slice(0, 11)), + [0, 1, 0, 0, 1, -0, -0, 0, -0, -0, -1], + ); +}); diff --git a/src/tooling/importDirections.test.mjs b/src/tooling/importDirections.test.mjs index b34fee3b9..0685c181b 100644 --- a/src/tooling/importDirections.test.mjs +++ b/src/tooling/importDirections.test.mjs @@ -82,6 +82,24 @@ for (const [name, files, pattern] of [ }, /reaches browser globals/, ], + [ + 'geo module imports cesium', + { 'src/geo/demo.js': "import 'cesium';" }, + /src\/geo\/demo\.js: geo reaches platform\/rendering dependency: cesium/, + ], + [ + 'geo module imports a Cesium workspace package', + { 'src/geo/demo.js': "import { Cartesian3 } from '@cesium/engine';" }, + /src\/geo\/demo\.js: geo reaches platform\/rendering dependency: @cesium\/engine/, + ], + [ + 'geo module reaches cesium through a helper', + { + 'src/geo/demo.js': "import '../data/helper.js';", + 'src/data/helper.js': "export * from 'cesium';", + }, + /src\/data\/helper\.js: geo reaches platform\/rendering dependency: cesium/, + ], [ 'browser imports Node without prefix', { 'src/unused.mjs': "import 'fs/promises';" }, @@ -149,6 +167,14 @@ for (const [name, files, pattern] of [ }); } +test('geo modules may import each other, and its tests may import cesium', (t) => { + const { root, write } = fixture(t); + write('src/geo/math.js', 'export const half = (value) => value / 2;'); + write('src/geo/demo.js', "export { half } from './math.js';"); + write('src/geo/demo.test.mjs', "import 'cesium'; import './demo.js';"); + assert.doesNotThrow(() => checkImportDirections(root)); +}); + test('self package imports cannot evade portable ownership', (t) => { const { root, write } = fixture(t, { './view': './src/ui/view.js' }); write('src/sources/demo.js', "import '@gods-eye-view/core/view';");