Repository: NASA-AMMOS/3DTilesRendererJS
Pull request #349 (2023-06-06 04:20:36)

More google tiles updates
Related to #331 

- Cleanup
- Add Ellipsoid function for retrieving lat lon from position
- Jump to hash lat lon

**TODO**
- ~Ensure Ellipsoid does not have to "swap frame"~
- Add tests


Patch:
diff --git a/example/googleMapsExample.js b/example/googleMapsExample.js
index 822d8d0..9e9e652 100644
--- a/example/googleMapsExample.js
+++ b/example/googleMapsExample.js
@@ -5,10 +5,7 @@ import {
 	AmbientLight,
 	WebGLRenderer,
 	PerspectiveCamera,
-	Vector3,
-	sRGBEncoding,
 	Raycaster,
-	Vector2,
 	Box3,
 } from 'three';
 import { GLTFLoader } from 'three/examples/jsm/loaders/GLTFLoader.js';
@@ -19,7 +16,6 @@ import Stats from 'three/examples/jsm/libs/stats.module.js';
 
 import { MapControls } from './src/lib/MapControls.js';
 import { MapsTilesCredits } from './src/MapsTilesCredits.js';
-import { GeoCoord } from './src/GeoCoord.js';
 import * as GeoUtils from './src/GeoUtils.js';
 
 const apiOrigin = 'https://tile.googleapis.com';
@@ -27,13 +23,8 @@ const apiOrigin = 'https://tile.googleapis.com';
 let camera, controls, scene, renderer, tiles, credits;
 let statsContainer, stats;
 
-const ellipsoid = GeoUtils.WGS84_ELLIPSOID;
-
 const raycaster = new Raycaster();
 raycaster.firstHitOnly = true;
-const pointer = new Vector2();
-const deltaTarget = new Vector3();
-const geocoord = new GeoCoord();
 
 const params = {
 
@@ -171,6 +162,8 @@ function reinstantiateTiles() {
 			tiles.group.rotation.x = - Math.PI / 2;
 
 			setupTiles();
+			initFromHash();
+			setInterval( updateHash, 100 );
 
 		} )
 		.catch( err => {
@@ -190,7 +183,6 @@ function init() {
 	renderer.setPixelRatio( window.devicePixelRatio );
 	renderer.setSize( window.innerWidth, window.innerHeight );
 	renderer.setClearColor( 0x151c1f );
-	renderer.outputEncoding = sRGBEncoding;
 
 	document.body.appendChild( renderer.domElement );
 
@@ -276,8 +268,6 @@ function onWindowResize() {
 
 function updateControls() {
 
-	// TODO: add tilting of the camera - keep a frame at the surface to retain the camera orientation
-
 	const raycaster = new Raycaster();
 	raycaster.ray.origin.copy( controls.target ).normalize().multiplyScalar( GeoUtils.WGS84_RADIUS * 1.5 );
 	raycaster.ray.direction.copy( raycaster.ray.origin ).normalize().multiplyScalar( - 1 );
@@ -309,16 +299,35 @@ function updateControls() {
 
 }
 
-window.setLatLon = function setLatLon( lat, lon ) {
+function updateHash() {
 
-	geocoord.lat = lat;
-	geocoord.lon = lon;
+	const res = {};
+	const pos = controls.target.clone();
+	GeoUtils.swapToGeoFrame( pos );
 
-	// geocoord.getVector3( controls.target );
-	ellipsoid.getCartographicToPosition( lat, lon, 0, controls.target );
-	controls.target.applyMatrix4( tiles.group.matrixWorld );
+	GeoUtils.WGS84_ELLIPSOID.getPositionToCartographic( pos, res );
+	window.history.replaceState( undefined, undefined, `#${ res.lat.toFixed( 4 ) },${ res.lon.toFixed( 4 ) }` );
 
-};
+}
+
+function initFromHash() {
+
+	const hash = window.location.hash.replace( /^#/, '' );
+	const tokens = hash.split( /,/g ).map( t => parseFloat( t ) );
+	console.log( tokens );
+	if ( tokens.length !== 2 || tokens.findIndex( t => Number.isNaN( t ) ) !== - 1 ) {
+
+		return;
+
+	}
+
+	const [ lat, lon ] = tokens;
+	GeoUtils.WGS84_ELLIPSOID.getCartographicToPosition( lat, lon, 0, controls.target );
+	GeoUtils.swapToThreeFrame( controls.target );
+
+	updateControls();
+
+}
 
 function animate() {
 
@@ -412,7 +421,11 @@ function render() {
 
 		const mat = tiles.group.matrixWorld.clone().invert();
 		const vec = camera.position.clone().applyMatrix4( mat );
-		document.getElementById( 'credits' ).innerText = geocoord.fromVector3( vec ).toString() + '\n' + credits.getCredits();
+
+
+		const res = {};
+		GeoUtils.WGS84_ELLIPSOID.getPositionToCartographic( vec, res );
+		document.getElementById( 'credits' ).innerText = GeoUtils.toLatLonString( res.lat, res.lon ) + '\n' + credits.getCredits();
 
 	}
 
diff --git a/example/src/GeoCoord.js b/example/src/GeoCoord.js
deleted file mode 100644
index e7d0daf..0000000
--- a/example/src/GeoCoord.js
+++ /dev/null
@@ -1,97 +0,0 @@
-import { MathUtils, Spherical, Vector3 } from 'three';
-import { swapToGeoFrame, swapToThreeFrame, sphericalPhiToLatitude, latitudeToSphericalPhi } from './GeoUtils.js';
-
-function getHoursMinutesSeconds( value, pos = 'E', neg = 'W' ) {
-
-	const direction = value < 0 ? neg : pos;
-	value = Math.abs( value );
-
-	const hours = ~ ~ value;
-
-	const minDec = ( value - hours ) * 60;
-	const minutes = ~ ~ minDec;
-
-	const secDec = ( minDec - minutes ) * 60;
-	const seconds = ~ ~ secDec;
-
-	return `${ hours }° ${ minutes }' ${ seconds }" ${ direction }`
-
-
-}
-
-// lat = [ -90, 90 ], lon = [ -180, 180 ]
-const _vec = new Vector3();
-const _spherical = new Spherical();
-export class GeoCoord {
-
-	constructor( lat = 0, lon = 0 ) {
-
-		this.lat = lat;
-		this.lon = lon;
-
-	}
-
-	fromSpherical( spherical ) {
-
-		this.lon = spherical.theta;
-		this.lat = sphericalPhiToLatitude( spherical.phi );
-		return this;
-
-	}
-
-	getSpherical( target ) {
-
-		target.theta = this.lon;
-		target.phi = latitudeToSphericalPhi( this.lat );
-		return this;
-
-	}
-
-	fromVector3( vector ) {
-
-		swapToThreeFrame( vector );
-
-		_spherical.setFromVector3( vector );
-		this.fromSpherical( _spherical );
-		return this;
-
-	}
-
-	getVector3( target ) {
-
-		this.getSpherical( _spherical );
-		target.setFromSpherical( _spherical );
-
-		swapToGeoFrame( target );
-
-		return this;
-
-	}
-
-	normalize() {
-
-		this.getVector3( _vec );
-		this.fromVector3( _vec );
-		return this;
-
-	}
-
-	toLatitudeString() {
-
-		return getHoursMinutesSeconds( MathUtils.RAD2DEG * this.lat, 'N', 'S' );
-
-	}
-
-	toLongitudeString() {
-
-		return getHoursMinutesSeconds( MathUtils.RAD2DEG * this.lon, 'E', 'W' );
-
-	}
-
-	toString() {
-
-		return `${ this.toLatitudeString() }, ${ this.toLongitudeString() }`;
-
-	}
-
-}
diff --git a/example/src/GeoUtils.js b/example/src/GeoUtils.js
index 584d91e..3be7d37 100644
--- a/example/src/GeoUtils.js
+++ b/example/src/GeoUtils.js
@@ -1,5 +1,10 @@
+import { Spherical, Vector3, MathUtils } from 'three';
 import { Ellipsoid } from '../../';
 
+const _spherical = new Spherical();
+const _vec = new Vector3();
+const _geoResults = {};
+
 export const WGS84_RADIUS = 6378137;
 
 export const WGS84_FLATTENING = 1 / 298.257223563;
@@ -37,3 +42,53 @@ export function latitudeToSphericalPhi( latitude ) {
 	return - latitude + Math.PI / 2;
 
 }
+
+export function correctGeoCoordWrap( lat, lon, target = {} ) {
+
+	_spherical.theta = lon;
+	_spherical.phi = latitudeToSphericalPhi( lat );
+	_vec.setFromSpherical( _spherical );
+
+	_spherical.setFromVector3( _vec );
+	target.lat = sphericalPhiToLatitude( _spherical.phi );
+	target.lon = _spherical.theta;
+	return target;
+
+}
+
+function toHoursMinutesSecondsString( value, pos = 'E', neg = 'W' ) {
+
+	const direction = value < 0 ? neg : pos;
+	value = Math.abs( value );
+
+	const hours = ~ ~ value;
+
+	const minDec = ( value - hours ) * 60;
+	const minutes = ~ ~ minDec;
+
+	const secDec = ( minDec - minutes ) * 60;
+	const seconds = ~ ~ secDec;
+
+	return `${ hours }° ${ minutes }' ${ seconds }" ${ direction }`
+
+}
+
+export function toLatLonString( lat, lon, decimalFormat = false ) {
+
+	const result = correctGeoCoordWrap( lat, lon, _geoResults );
+	let latString, lonString;
+	if ( decimalFormat ) {
+
+		latString =  `${ ( MathUtils.RAD2DEG * result.lat ).toFixed( 4 ) }°`;
+		lonString =  `${ ( MathUtils.RAD2DEG * result.lon ).toFixed( 4 ) }°`;
+
+	} else {
+
+		latString = toHoursMinutesSecondsString( MathUtils.RAD2DEG * result.lat, 'N', 'S' );
+		lonString = toHoursMinutesSecondsString( MathUtils.RAD2DEG * result.lon, 'E', 'W' );
+
+	}
+
+	return `${ latString } ${ lonString }`;
+
+}
diff --git a/src/three/math/Ellipsoid.d.ts b/src/three/math/Ellipsoid.d.ts
index 0cb9fd5..8493a04 100644
--- a/src/three/math/Ellipsoid.d.ts
+++ b/src/three/math/Ellipsoid.d.ts
@@ -10,6 +10,7 @@ export class Ellipsoid {
 
 	constructor( x: Number, y: Number, z: Number );
 	getCartographicToPosition( lat: Number, lon: Number, height: Number, target: Vector3 ): Vector3;
+	getPositionToCartographic( pos: Vector3, target: Object ): Vector3;
 	getCartographicToNormal( lat: Number, lon: Number, target: Vector3 ): Vector3;
 	getPositionToNormal( pos: Vector3, target: Vector3 ): Vector3;
 
diff --git a/src/three/math/Ellipsoid.js b/src/three/math/Ellipsoid.js
index a08ae0b..901348d 100644
--- a/src/three/math/Ellipsoid.js
+++ b/src/three/math/Ellipsoid.js
@@ -45,6 +45,10 @@ export function latitudeToSphericalPhi( latitude ) {
 const _spherical = new Spherical();
 const _norm = new Vector3();
 const _vec = new Vector3();
+const _vec2 = new Vector3();
+
+const EPSILON12 = 1e-12;
+const CENTER_EPS = 0.1;
 
 export class Ellipsoid {
 
@@ -56,13 +60,11 @@ export class Ellipsoid {
 
 	getCartographicToPosition( lat, lon, height, target ) {
 
-		// https://github.com/CesiumGS/cesium/blob/main/Source/Core/Ellipsoid.js#L396
-		const radius = this.radius;
-		_spherical.set( 1, latitudeToSphericalPhi( lat ), lon );
-		_norm.setFromSpherical( _spherical ).normalize();
-
-		swapFrame( _norm );
+		// From Cesium function Ellipsoid.cartographicToCartesian
+		// https://github.com/CesiumGS/cesium/blob/665ec32e813d5d6fe906ec3e87187f6c38ed5e49/packages/engine/Source/Core/Ellipsoid.js#L396
+		this.getCartographicToNormal( lat, lon, _norm );
 
+		const radius = this.radius;
 		_vec.copy( _norm );
 		_vec.x *= radius.x ** 2;
 		_vec.y *= radius.y ** 2;
@@ -75,11 +77,28 @@ export class Ellipsoid {
 
 	}
 
+	getPositionToCartographic( pos, target ) {
+
+		// From Cesium function Ellipsoid.cartesianToCartographic
+		// https://github.com/CesiumGS/cesium/blob/665ec32e813d5d6fe906ec3e87187f6c38ed5e49/packages/engine/Source/Core/Ellipsoid.js#L463
+		this.getPositionToSurfacePoint( pos, _vec );
+		this.getPositionToNormal( pos, _norm );
+
+		const heightDelta = _vec2.subVectors( pos, _vec );
+
+		target.lon = Math.atan2( _norm.y, _norm.x );
+		target.lat = Math.asin( _norm.z );
+		target.height = Math.sign( heightDelta.dot( pos ) ) * heightDelta.length();
+		return target;
+
+	}
+
 	getCartographicToNormal( lat, lon, target ) {
 
-		_spherical.set( 1, ( - lat + Math.PI / 2 ), lon );
+		_spherical.set( 1, latitudeToSphericalPhi( lat ), lon );
 		target.setFromSpherical( _spherical ).normalize();
 
+		// swap frame from the three.js frame to the geo coord frame
 		swapFrame( target );
 		return target;
 
@@ -100,12 +119,81 @@ export class Ellipsoid {
 
 	getPositionToSurfacePoint( pos, target ) {
 
-		const normal = this.getPositionToNormal( pos, target );
-		normal.x *= this.radius.x;
-		normal.y *= this.radius.y;
-		normal.z *= this.radius.z;
+		// From Cesium function Ellipsoid.scaleToGeodeticSurface
+		// https://github.com/CesiumGS/cesium/blob/d11b746e5809ac115fcff65b7b0c6bdfe81dcf1c/packages/engine/Source/Core/scaleToGeodeticSurface.js#L25
+		const radius = this.radius;
+		const invRadiusSqX = 1 / ( radius.x ** 2 );
+		const invRadiusSqY = 1 / ( radius.y ** 2 );
+		const invRadiusSqZ = 1 / ( radius.z ** 2 );
+
+		const x2 = pos.x * pos.x * invRadiusSqX;
+		const y2 = pos.y * pos.y * invRadiusSqY;
+		const z2 = pos.z * pos.z * invRadiusSqZ;
 
-		return target;
+		// 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 = _vec.copy( pos ).multiplyScalar( ratio );
+		if ( squaredNorm < CENTER_EPS ) {
+
+			return ! isFinite( ratio ) ? null : target.copy( intersection );
+
+		}
+
+		// 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 = _vec2.set(
+			intersection.x * invRadiusSqX * 2.0,
+			intersection.y * invRadiusSqY * 2.0,
+			intersection.z * invRadiusSqZ * 2.0
+		);
+
+		// Compute the initial guess at the normal vector multiplier, lambda.
+		let lambda = ( 1.0 - ratio ) * pos.length() / ( 0.5 * gradient.length() );
+		let correction = 0.0;
+
+		let func, denominator;
+		let xMultiplier, yMultiplier, zMultiplier;
+		let xMultiplier2, yMultiplier2, zMultiplier2;
+		let xMultiplier3, yMultiplier3, zMultiplier3;
+
+		do {
+
+			lambda -= correction;
+
+			xMultiplier = 1.0 / ( 1.0 + lambda * invRadiusSqX );
+			yMultiplier = 1.0 / ( 1.0 + lambda * invRadiusSqY );
+			zMultiplier = 1.0 / ( 1.0 + lambda * invRadiusSqZ );
+
+			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 * invRadiusSqX +
+			  	y2 * yMultiplier3 * invRadiusSqY +
+			  	z2 * zMultiplier3 * invRadiusSqZ;
+
+			const derivative = - 2.0 * denominator;
+			correction = func / derivative;
+
+		} while ( Math.abs( func ) > EPSILON12 );
+
+		return target.set(
+			pos.x * xMultiplier,
+			pos.y * yMultiplier,
+			pos.z * zMultiplier
+		);
 
 	}
 


Test patch:
diff --git a/test/Ellipsoid.test.js b/test/Ellipsoid.test.js
index 47d79ab..6c2a820 100644
--- a/test/Ellipsoid.test.js
+++ b/test/Ellipsoid.test.js
@@ -4,23 +4,55 @@ import { EllipsoidRegion } from '../src/three/math/EllipsoidRegion.js';
 import { Ellipsoid } from '../src/three/math/Ellipsoid.js';
 import { WGS84_HEIGHT, WGS84_RADIUS } from '../src/base/constants.js';
 
-function compareCartesianVector( cart, vec, EPSILON = 1e-13 ) {
+function epsCompare( a, b, EPSILON = 1e-10 ) {
 
-	expect( Math.abs( cart.x - vec.x ) ).toBeLessThanOrEqual( EPSILON );
-	expect( Math.abs( cart.y - vec.y ) ).toBeLessThanOrEqual( EPSILON );
-	expect( Math.abs( cart.z - vec.z ) ).toBeLessThanOrEqual( EPSILON );
+	expect( Math.abs( a - b ) ).toBeLessThanOrEqual( EPSILON );
+
+}
+
+function compareCartesianVector( cart, vec, EPSILON = 1e-10 ) {
+
+	epsCompare( cart.x, vec.x, EPSILON );
+	epsCompare( cart.y, vec.y, EPSILON );
+	epsCompare( cart.z, vec.z, EPSILON );
+
+}
+
+function compareCoord( a, b, EPSILON = 1e-7 ) {
+
+	epsCompare( a, b, EPSILON );
+
+}
+
+function randomRange( min, max ) {
+
+	const delta = max - min;
+	return min + Math.random() * delta;
+
+}
+
+function randomLat() {
+
+	return randomRange( - Math.PI / 2, Math.PI / 2 );
+
+}
+
+function randomLon() {
+
+	return randomRange( - Math.PI, Math.PI );
 
 }
 
 describe( 'Ellipsoid', () => {
 
-	let v, c;
+	let v, c, c2;
 	let c_unitEllipse, unitEllipse;
 	let c_wgsEllipse, wgsEllipse;
 	let norm, norm2, cnorm, cnorm2;
 	beforeEach( () => {
 
 		c = new Cesium.Cartesian3();
+		c2 = new Cesium.Cartesian3();
 		v = new Vector3();
 
 		norm = new Vector3();
@@ -36,28 +68,80 @@ describe( 'Ellipsoid', () => {
 
 	} );
 
+	it( 'should convert between lat lon and position consistently.', () => {
+
+		const POS_EPSILON = 1e-6;
+		const SURF_EPSILON = 1e-5;
+		const LAT_LON_EPSILON = 1e-4;
+		for ( let i = 0; i < 100; i ++ ) {
+
+			const lat = randomLat();
+			const lon = randomLon();
+			const height = randomRange( - 1e5, 1e5 );
+			const cart = new Cesium.Cartographic( lon, lat, height );
+
+			wgsEllipse.getCartographicToPosition( lat, lon, height, v );
+			c_wgsEllipse.cartographicToCartesian( cart, c );
+
+			compareCartesianVector( c, v, POS_EPSILON );
+
+			const result = {};
+			wgsEllipse.getPositionToCartographic( v, result );
+			c_wgsEllipse.scaleToGeodeticSurface( c, c2 );
+
+			compareCoord( result.lat, lat, LAT_LON_EPSILON );
+			compareCoord( result.lon, lon, LAT_LON_EPSILON );
+			compareCoord( result.height, height, SURF_EPSILON );
+
+		}
+
+	} );
+
+	it( 'should match the surface points.', () => {
+
+		const LOCAL_EPSILON = 1e-8;
+		for ( let i = 0; i < 100; i ++ ) {
+
+			v.random().normalize().multiplyScalar( WGS84_RADIUS * 2 );
+			c.x = v.x;
+			c.y = v.y;
+			c.z = v.z;
+
+			wgsEllipse.getPositionToSurfacePoint( v, norm );
+			c_wgsEllipse.scaleToGeodeticSurface( c, cnorm );
+
+			compareCartesianVector( cnorm, norm, LOCAL_EPSILON );
+
+		}
+
+	} );
+
 	it( 'should match Cesium Unit Sphere results.', () => {
 
 		for ( let i = 0; i < 100; i ++ ) {
 
-			const lat = Math.random();
-			const lon = Math.random();
-			const height = Math.random();
-			const cart = new Cesium.Cartographic( lon, lat, 0 );
+			const lat = randomLat();
+			const lon = randomLon();
+			const height = randomRange( - 1e5, 1e5 );
+			const cart = new Cesium.Cartographic( lon, lat, height );
 
 			// test positions
-			c_unitEllipse.cartographicToCartesian( new Cesium.Cartographic( lon, lat, height ), c );
+			c_unitEllipse.cartographicToCartesian( cart, c );
 			unitEllipse.getCartographicToPosition( lat, lon, height, v );
 
 			compareCartesianVector( c, v );
 
-			// test normals
-			c_unitEllipse.geodeticSurfaceNormal( c, cnorm2 );
-			c_unitEllipse.geodeticSurfaceNormalCartographic( cart, cnorm );
+			// test pos to normal
+			c_unitEllipse.geodeticSurfaceNormal( c, cnorm );
 			unitEllipse.getPositionToNormal( v, norm );
 
 			compareCartesianVector( cnorm, norm );
-			compareCartesianVector( cnorm, cnorm2 );
+
+			// test cart to normal
+			c_wgsEllipse.geodeticSurfaceNormalCartographic( cart, cnorm );
+			wgsEllipse.getCartographicToNormal( lat, lon, norm );
+
+			compareCartesianVector( cnorm, norm );
 
 		}
 
@@ -68,10 +152,9 @@ describe( 'Ellipsoid', () => {
 		const LOCAL_EPSILON = 1e-6;
 		for ( let i = 0; i < 100; i ++ ) {
 
-			const lat = Math.random();
-			const lon = Math.random();
-			const height = ( Math.random() - 0.5 ) * 1000;
-
+			const lat = randomLat();
+			const lon = randomLon();
+			const height = randomRange( - 1e5, 1e5 );
 			const cart = new Cesium.Cartographic( lon, lat, 0 );
 
 			// test positions
@@ -80,15 +163,17 @@ describe( 'Ellipsoid', () => {
 
 			compareCartesianVector( c, v, LOCAL_EPSILON );
 
-			// test normals
-			c_wgsEllipse.geodeticSurfaceNormal( c, cnorm2 );
+			// test pos to normal
+			c_unitEllipse.geodeticSurfaceNormal( c, cnorm );
+			unitEllipse.getPositionToNormal( v, norm );
+
+			compareCartesianVector( cnorm, norm );
+
+			// test cart to normal
 			c_wgsEllipse.geodeticSurfaceNormalCartographic( cart, cnorm );
-			wgsEllipse.getPositionToNormal( v, norm );
-			wgsEllipse.getCartographicToNormal( lat, lon, norm2 );
+			wgsEllipse.getCartographicToNormal( lat, lon, norm );
 
-			compareCartesianVector( cnorm, norm, LOCAL_EPSILON );
-			compareCartesianVector( cnorm, cnorm2, LOCAL_EPSILON );
-			compareCartesianVector( cnorm, norm2, LOCAL_EPSILON );
+			compareCartesianVector( cnorm, norm );
 
 		}