From e32988ae5e79911c2a5790ef65f44989b27ffc3f Mon Sep 17 00:00:00 2001 From: Sameh Khamis Date: Thu, 5 Nov 2026 14:28:10 -0800 Subject: [PATCH] Match Cesium's fromDegrees exactly and keep toEcef in the stub fromRadians divides by gamma as Cesium does, so geodetic to ECEF is now bit for bit Cesium's, and the test checks it exactly. Its tests also cover three height bands, points inside the Earth and NaN input. toEcef reads the renderer contract's Position, so it moves to the stub's maths. The ellipsoid docs pair each conversion with Cesium's. --- src/geo/ellipsoid.js | 21 ++++++++++--------- src/geo/geo.test.mjs | 45 ++++++++++++++++++++++++++++++++++++++--- src/geo/vector.js | 1 + src/render/stub/math.js | 9 +++++++++ 4 files changed, 64 insertions(+), 12 deletions(-) diff --git a/src/geo/ellipsoid.js b/src/geo/ellipsoid.js index 45d67da0d..dbff51e3f 100644 --- a/src/geo/ellipsoid.js +++ b/src/geo/ellipsoid.js @@ -4,10 +4,20 @@ * (ECEF) conversions, surface normals and projection onto the surface. The * formulas follow the default renderer's step for step, so they agree with * it to floating-point noise. + * + * Into ECEF, positions are arguments (`fromDegrees(lon, lat, height)`, + * `fromRadians` likewise); out of it, objects (`toLonLatH` gives + * `{lon, lat, height}` in degrees, `toCartographicRadians` the same in + * radians). Cesium's `Cartographic.fromDegrees(lon, lat, h)` is plain + * `{lon: toRadians(lon), lat: toRadians(lat), height: h}`; for its ECEF + * point, use `fromDegrees`. + * + * NaN in gives NaN fields out, never a throw (Cesium's debug build throws). + * The Earth's centre, and points so near it that their squares round to + * zero, answer null, as Cesium answers undefined. */ import { add, - clone, dot, magnitude, normalize, @@ -58,7 +68,7 @@ export function fromRadians(lon, lat, height = 0) { const r = WGS84.radiiSquared; let k = vec(r.x * n.x, r.y * n.y, r.z * n.z); const gamma = Math.sqrt(dot(n, k)); - k = scale(k, 1 / gamma); + k = vec(k.x / gamma, k.y / gamma, k.z / gamma); return add(k, scale(n, height)); } @@ -66,13 +76,6 @@ export function fromRadians(lon, lat, height = 0) { export const fromDegrees = (lon, lat, height = 0) => fromRadians(toRadians(lon), toRadians(lat), height); -/** A position in either contract form, as ECEF. */ -export function toEcef(position) { - if ('lon' in position && 'lat' in position) - return fromDegrees(position.lon, position.lat, position.height ?? 0); - return clone(position); -} - /** The point on the ellipsoid along the geodetic normal, or null at the centre. */ export function scaleToGeodeticSurface(p) { const { x: positionX, y: positionY, z: positionZ } = p; diff --git a/src/geo/geo.test.mjs b/src/geo/geo.test.mjs index 9e0ad1a58..e0fcde931 100644 --- a/src/geo/geo.test.mjs +++ b/src/geo/geo.test.mjs @@ -29,12 +29,20 @@ function points(count = 500) { { lon: 180, lat: 0, height: -400 }, { lon: -97.7431, lat: 30.2672, height: 150 }, ]; - for (let i = 0; i < count; i++) + // Each height band gets its share: near the ground, low orbit, far out. + const bands = [ + [-500, 1e4], + [1e4, 1e6], + [1e6, 4e7], + ]; + for (let i = 0; i < count; i++) { + const [lo, hi] = bands[i % bands.length]; out.push({ lon: next() * 360 - 180, lat: next() * 180 - 90, - height: next() * 4e7 - 1e4, + height: lo + next() * (hi - lo), }); + } return out; } @@ -60,7 +68,7 @@ test('geodetic to ECEF matches the default renderer', () => { for (const p of points()) { const ours = fromDegrees(p.lon, p.lat, p.height); const theirs = Cesium.Cartesian3.fromDegrees(p.lon, p.lat, p.height); - closeVec(ours, theirs, 1e-6, `${p.lon},${p.lat},${p.height}`); + closeVec(ours, theirs, 0, `${p.lon},${p.lat},${p.height}`); } }); @@ -76,6 +84,37 @@ test('ECEF to geodetic matches the default renderer', () => { close(ours.lat, theirs.latitude, 1e-12, `${where} lat`); close(ours.height, theirs.height, 1e-6, `${where} height`); } + // Inside the Earth: within the centre cutoff (the first four) and between + // it and the surface (the last). + for (const p of [ + { x: 1e6, y: 0, z: 0 }, + { x: 1, y: 2, z: 3 }, + { x: -2e5, y: 3e5, z: 1e5 }, + { x: 1e6, y: 1e6, z: 5e5 }, + { x: 2e6, y: 1.5e6, z: 1.5e6 }, + ]) { + const ours = toCartographicRadians(p); + const theirs = Cesium.Cartographic.fromCartesian( + new Cesium.Cartesian3(p.x, p.y, p.z), + ); + const where = `${p.x},${p.y},${p.z}`; + close(ours.lon, theirs.longitude, 1e-12, `${where} lon`); + close(ours.lat, theirs.latitude, 1e-12, `${where} lat`); + close(ours.height, theirs.height, 1e-6, `${where} height`); + } + for (const p of [ + { x: NaN, y: 0, z: 0 }, + { x: 7e6, y: NaN, z: 0 }, + { x: 7e6, y: 0, z: NaN }, + ]) { + const nan = toLonLatH(p); + assert.ok( + Number.isNaN(nan.lon) && + Number.isNaN(nan.lat) && + Number.isNaN(nan.height), + JSON.stringify(nan), + ); + } assert.equal(toCartographicRadians({ x: 0, y: 0, z: 0 }), null); assert.equal( Cesium.Cartographic.fromCartesian(Cesium.Cartesian3.ZERO), diff --git a/src/geo/vector.js b/src/geo/vector.js index e020dd5cc..a5b7400fc 100644 --- a/src/geo/vector.js +++ b/src/geo/vector.js @@ -18,6 +18,7 @@ export const dot = (a, b) => a.x * b.x + a.y * b.y + a.z * b.z; export const cross = (a, b) => vec(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x); export const magnitude = (a) => Math.sqrt(dot(a, a)); +/** Straight-line (chord) distance, not along the surface. */ export const distance = (a, b) => magnitude(sub(a, b)); export function normalize(a) { const m = magnitude(a); diff --git a/src/render/stub/math.js b/src/render/stub/math.js index 4fa6665b1..896309eca 100644 --- a/src/render/stub/math.js +++ b/src/render/stub/math.js @@ -14,6 +14,8 @@ import { clamp, equalsEpsilon, negate, + clone, + fromDegrees, rotate, scale, vec, @@ -22,6 +24,13 @@ import { export * from '../../geo/index.js'; +/** A contract `Position` (`Ecef` or `LonLatH`, degrees), as ECEF. */ +export function toEcef(position) { + if ('lon' in position && 'lat' in position) + return fromDegrees(position.lon, position.lat, position.height ?? 0); + return clone(position); +} + const EPSILON3 = 1e-3; const EPSILON10 = 1e-10;