Repository: NASA-AMMOS/3DTilesRendererJS
Pull request #483 (2024-02-16 11:36:37)

GlobeControls limit Frustum with horizon for performance
Currently for simplicity GlobeControls camera far is far enough to include the entire earth resulting in unnecessary tiles being loaded.

With these changes we limit the frustum to the distance of the horizon based on our current latitude and elevation.

This seems to work well but some tiles are being loaded while they don't intersect the frustum.

I edited the googleMapsExample.js to visualize it.

I'm not sure where is the issue but I'm pretty sure it's either the frustum or the OOB.
This particular example might help identifying the issue, when the camera is zoomed in on the north pole we get false positive on the intersectsFrustum check. 
If we look at it from the side, it's considered like an intersection as it seems like the frustum far limit does intersect with the obbs.
![wrongFrustumCheck](https://github.com/NASA-AMMOS/3DTilesRendererJS/assets/7174039/f309ce0f-8af7-4d28-a49d-6a841c0d1f68)
But if we look at it from above, we can see that the frustum doesn't.
![wrongFrustumCheckFromAbove](https://github.com/NASA-AMMOS/3DTilesRendererJS/assets/7174039/c47d7258-eca8-430f-a307-cb8c7ec635a5)

Those false positive unload as soon as on the side view, the frustum move outside of those obb.

gkjohnson: Awesome! I'll take a deeper look at this next week. It would be nice to have some before / after stats on the number of tiles loaded and number of tiles rendered after the changes you've made so far.

> I'm not sure where is the issue but I'm pretty sure it's either the frustum or the OOB.
> This particular example might help identifying the issue, when the camera is zoomed in on the north pole we get false positive on the intersectsFrustum check.

I took a look and the issue is in the OBB frustum intersect code. In fact three.js' AABB frustum intersect code that the OBB logic is based on isn't correct, either. It will mark bounds as in the frustum even when they're clearly outside of it since it's only checking if the box is completely on one side of any plane which misses the common case of a box that intersects multiple planes if it's large. Here's a demo:

<img width="226" alt="image" src="https://github.com/NASA-AMMOS/3DTilesRendererJS/assets/734200/bb5e3474-2d71-4c77-ba81-7fc6e23f0648">

https://jsfiddle.net/89rtuaqp/2/

[This post](https://iquilezles.org/articles/frustumcorrect/) discusses the problem and a correct application of separating axis theorem. I'll make a three.js issue. We can see what they say - we may have to make our own extended frustum class that stores enough data for performing a more correct OBB intersection test.


gkjohnson: https://github.com/mrdoob/three.js/issues/27756 may take a bit to resolve so I think it's best for us to modify our current implementation to be more robust, but I think that can happen in another PR. This can be done by storing the frustum corners as points on the Frustum so they can be tested against the OBB planes as referenced in the inigo post.
AlaricBaraou: Good timing, I was working on a custom Frustum but it's hard to test it with the tiles logic, I might try in a simpler environment.

I added a class AdvFrustum that calculate the points of the Frustum to use them for different operations as suggested in the article you linked. But I think there are still case that are not handled or it might be the logic of the nested tiles etc that I don't really grasp at the moment.

Also the latest commit take care of "Google Maps Example: Add ability to specify height in hash" #482 
AlaricBaraou: a simpler example to play with it 

https://jsfiddle.net/0sb3gden/30/

The extended Frustum works well with axis aligned box but not with OBB at the moment.
AlaricBaraou: I was about to write that it's good now but apparently now we have some false negative.

Looking good here
![frustumOk](https://github.com/NASA-AMMOS/3DTilesRendererJS/assets/7174039/3cfccb54-a869-4cbe-88ca-15cb5db3ddb4)

But here it discard some tiles that are in obvious intersections. 
![frustumNOk](https://github.com/NASA-AMMOS/3DTilesRendererJS/assets/7174039/636b02c9-9e88-40cd-a1b8-380f2b629b1f)

I`ll try to reproduce that one in jsfiddle

AlaricBaraou: I addressed all points, just a note about this issue https://github.com/NASA-AMMOS/3DTilesRendererJS/issues/482 
It can't be init under a certain height, at least in the area I was testing it.
I added a comment to check later
gkjohnson: > It can't be init under a certain height, at least in the area I was testing it.

Likely what's happening is that the coarse levels of detail on the model load and have a much higher height than the lower LoDs which pushes the camera up due to the "GlobeControls" update function. I think it's okay to say that this fixes #482.

Patch:
diff --git a/example/googleMapsExample.js b/example/googleMapsExample.js
index 5faebad..e0ba6db 100644
--- a/example/googleMapsExample.js
+++ b/example/googleMapsExample.js
@@ -146,7 +146,10 @@ function updateHash() {
 
 	res.lat *= MathUtils.RAD2DEG;
 	res.lon *= MathUtils.RAD2DEG;
-	window.history.replaceState( undefined, undefined, `#${ res.lat.toFixed( 4 ) },${ res.lon.toFixed( 4 ) }` );
+
+	const elevation = WGS84_ELLIPSOID.getPositionElevation( vec );
+
+	window.history.replaceState( undefined, undefined, `#${ res.lat.toFixed( 4 ) },${ res.lon.toFixed( 4 ) },${ elevation.toFixed( 4 ) }` );
 
 }
 
@@ -154,17 +157,18 @@ function initFromHash() {
 
 	const hash = window.location.hash.replace( /^#/, '' );
 	const tokens = hash.split( /,/g ).map( t => parseFloat( t ) );
-	if ( tokens.length !== 2 || tokens.findIndex( t => Number.isNaN( t ) ) !== - 1 ) {
+	if ( tokens.length < 3 || tokens.findIndex( t => Number.isNaN( t ) ) !== - 1 ) {
 
 		return;
 
 	}
 
-	const [ lat, lon ] = tokens;
-	WGS84_ELLIPSOID.getCartographicToPosition( lat * MathUtils.DEG2RAD, lon * MathUtils.DEG2RAD, 0, camera.position );
+	//todo it looks like we can't init under a certain height, I assume it has something to do with the first tile loaded and collision
+	const [ lat, lon, height ] = tokens;
+	WGS84_ELLIPSOID.getCartographicToPosition( lat * MathUtils.DEG2RAD, lon * MathUtils.DEG2RAD, height, camera.position );
 
 	tiles.group.updateMatrixWorld();
-	camera.position.applyMatrix4( tiles.group.matrixWorld ).multiplyScalar( 2 );
+	camera.position.applyMatrix4( tiles.group.matrixWorld );
 	camera.lookAt( 0, 0, 0 );
 
 }
diff --git a/src/three/TilesRenderer.js b/src/three/TilesRenderer.js
index ddaf80d..f5a7bf9 100644
--- a/src/three/TilesRenderer.js
+++ b/src/three/TilesRenderer.js
@@ -16,6 +16,7 @@ import {
 import { raycastTraverse, raycastTraverseFirstHit } from './raycastTraverse.js';
 import { readMagicBytes } from '../utilities/readMagicBytes.js';
 import { TileBoundingVolume } from './math/TileBoundingVolume.js';
+import { ExtendedFrustum } from './math/ExtendedFrustum.js';
 
 const INITIAL_FRUSTUM_CULLED = Symbol( 'INITIAL_FRUSTUM_CULLED' );
 const tempMat = new Matrix4();
@@ -369,7 +370,7 @@ export class TilesRenderer extends TilesRendererBase {
 
 			cameraInfo.push( {
 
-				frustum: new Frustum(),
+				frustum: new ExtendedFrustum(),
 				isOrthographic: false,
 				sseDenominator: - 1, // used if isOrthographic:false
 				position: new Vector3(),
diff --git a/src/three/controls/GlobeControls.js b/src/three/controls/GlobeControls.js
index 051af81..c07419a 100644
--- a/src/three/controls/GlobeControls.js
+++ b/src/three/controls/GlobeControls.js
@@ -10,6 +10,7 @@ import { makeRotateAroundPoint } from './utils.js';
 
 const _invMatrix = new Matrix4();
 const _rotMatrix = new Matrix4();
+const _pos = new Vector3();
 const _vec = new Vector3();
 const _center = new Vector3();
 const _up = new Vector3();
@@ -18,6 +19,7 @@ const _right = new Vector3();
 const _targetRight = new Vector3();
 const _globalUp = new Vector3();
 const _quaternion = new Quaternion();
+const _latLon = {};
 
 const _pointer = new Vector2();
 const _prevPointer = new Vector2();
@@ -152,7 +154,20 @@ export class GlobeControls extends EnvironmentControls {
 		// update the projection matrix
 		const largestDistance = Math.max( ...ellipsoid.radius );
 		camera.near = Math.max( 1, distanceToCenter - largestDistance * 1.25 );
-		camera.far = distanceToCenter + largestDistance + 0.1;
+
+		// update the far plane to the horizon distance
+		const invMatrix = _invMatrix.copy( tilesGroup.matrixWorld ).invert();
+		_pos.copy( camera.position ).applyMatrix4( invMatrix );
+		ellipsoid.getPositionToCartographic( _pos, _latLon );
+		const elevation = ellipsoid.getPositionElevation( _pos );
+
+		// due to the level of precision of the tiles / sphere
+		// we can get an elevation that will be higher than the actual elevation.
+		// when we end up below the surface of the ellipsoid the elevation is returned as positive instead of negative
+		// resulting in a horizon distance that is too large.
+		const horizonDistance = ellipsoid.calculateHorizonDistance( _latLon.lat, elevation );
+		camera.far = horizonDistance + 0.1;
+
 		camera.updateProjectionMatrix();
 
 	}
diff --git a/src/three/math/Ellipsoid.js b/src/three/math/Ellipsoid.js
index 89bc8a6..6a03f8c 100644
--- a/src/three/math/Ellipsoid.js
+++ b/src/three/math/Ellipsoid.js
@@ -1,4 +1,4 @@
-import { Vector3, Spherical } from 'three';
+import { Vector3, Spherical, MathUtils } from 'three';
 import { swapToGeoFrame, latitudeToSphericalPhi } from './GeoUtils.js';
 
 const _spherical = new Spherical();
@@ -183,9 +183,9 @@ export class Ellipsoid {
 			// "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;
+				x2 * xMultiplier3 * invRadiusSqX +
+				y2 * yMultiplier3 * invRadiusSqY +
+				z2 * zMultiplier3 * invRadiusSqZ;
 
 			const derivative = - 2.0 * denominator;
 			correction = func / derivative;
@@ -200,4 +200,39 @@ export class Ellipsoid {
 
 	}
 
+	calculateHorizonDistance( latitude, elevation ) {
+
+		// from https://aty.sdsu.edu/explain/atmos_refr/horizon.html
+		// OG = sqrt ( 2 R h + h2 ) .
+		const effectiveRadius = this.calculateEffectiveRadius( latitude );
+		return Math.sqrt( 2 * effectiveRadius * elevation + elevation ** 2 );
+
+	}
+
+	calculateEffectiveRadius( latitude ) {
+
+		// This radius represents the distance from the center of the ellipsoid to the surface along the normal at the given latitude.
+		// from https://en.wikipedia.org/wiki/Earth_radius#Prime_vertical
+		// N = a / sqrt(1 - e^2 * sin^2(phi))
+		const semiMajorAxis = this.radius.x;
+		const semiMinorAxis = this.radius.z;
+		const eSquared = 1 - ( semiMinorAxis ** 2 / semiMajorAxis ** 2 );
+		const phi = latitude * MathUtils.DEG2RAD;
+
+		const sinPhiSquared = Math.sin( phi ) ** 2;
+		const N = semiMajorAxis / Math.sqrt( 1 - eSquared * sinPhiSquared );
+		return N;
+
+	}
+
+	getPositionElevation( pos ) {
+
+		this.getPositionToSurfacePoint( pos, _vec3 );
+
+		const elevation = _vec3.distanceTo( pos );
+
+		return elevation;
+
+	}
+
 }
diff --git a/src/three/math/ExtendedFrustum.js b/src/three/math/ExtendedFrustum.js
new file mode 100644
index 0000000..23121bc
--- /dev/null
+++ b/src/three/math/ExtendedFrustum.js
@@ -0,0 +1,65 @@
+import { Frustum, Matrix3, Vector3 } from 'three';
+
+const _mat3 = new Matrix3();
+
+// Solve a system of equations to find the point where the three planes intersect
+function findIntersectionPoint( plane1, plane2, plane3, target ) {
+
+	// Create the matrix A using the normals of the planes as rows
+	const A = _mat3.set(
+		plane1.normal.x, plane1.normal.y, plane1.normal.z,
+		plane2.normal.x, plane2.normal.y, plane2.normal.z,
+		plane3.normal.x, plane3.normal.y, plane3.normal.z
+	);
+
+	// Create the vector B using the constants of the planes
+	target.set( - plane1.constant, - plane2.constant, - plane3.constant );
+
+	// Solve for X by applying the inverse matrix to B
+	target.applyMatrix3( A.invert() );
+
+	return target;
+
+}
+
+class ExtendedFrustum extends Frustum {
+
+	constructor() {
+
+		super();
+		this.points = Array( 8 ).fill().map( () => new Vector3() );
+
+	}
+
+	setFromProjectionMatrix( m, coordinateSystem ) {
+
+		super.setFromProjectionMatrix( m, coordinateSystem );
+		this.calculateFrustumPoints();
+
+	}
+
+	calculateFrustumPoints() {
+
+		const { planes, points } = this;
+		const planeIntersections = [
+			[ planes[ 0 ], planes[ 3 ], planes[ 4 ] ], // Near top left
+			[ planes[ 1 ], planes[ 3 ], planes[ 4 ] ], // Near top right
+			[ planes[ 0 ], planes[ 2 ], planes[ 4 ] ], // Near bottom left
+			[ planes[ 1 ], planes[ 2 ], planes[ 4 ] ], // Near bottom right
+			[ planes[ 0 ], planes[ 3 ], planes[ 5 ] ], // Far top left
+			[ planes[ 1 ], planes[ 3 ], planes[ 5 ] ], // Far top right
+			[ planes[ 0 ], planes[ 2 ], planes[ 5 ] ], // Far bottom left
+			[ planes[ 1 ], planes[ 2 ], planes[ 5 ] ], // Far bottom right
+		];
+
+		planeIntersections.forEach( ( planes, index ) => {
+
+			findIntersectionPoint( planes[ 0 ], planes[ 1 ], planes[ 2 ], points[ index ] );
+
+		} );
+
+	}
+
+}
+
+export { ExtendedFrustum };
diff --git a/src/three/math/OBB.js b/src/three/math/OBB.js
index 45279ae..86a0269 100644
--- a/src/three/math/OBB.js
+++ b/src/three/math/OBB.js
@@ -1,4 +1,8 @@
-import { Matrix4, Box3, Vector3 } from 'three';
+import { Matrix4, Box3, Vector3, Plane } from 'three';
+
+const _worldMin = new Vector3();
+const _worldMax = new Vector3();
+const _norm = new Vector3();
 
 export class OBB {
 
@@ -8,6 +12,7 @@ export class OBB {
 		this.transform = transform.clone();
 		this.inverseTransform = new Matrix4();
 		this.points = new Array( 8 ).fill().map( () => new Vector3() );
+		this.planes = new Array( 6 ).fill().map( () => new Plane() );
 
 	}
 
@@ -37,6 +42,27 @@ export class OBB {
 
 		}
 
+		this.updatePlanes();
+
+	}
+
+	updatePlanes() {
+
+		_worldMin.copy( this.box.min ).applyMatrix4( this.transform );
+		_worldMax.copy( this.box.max ).applyMatrix4( this.transform );
+
+		_norm.set( 0, 0, 1 ).transformDirection( this.transform );
+		this.planes[ 0 ].setFromNormalAndCoplanarPoint( _norm, _worldMin );
+		this.planes[ 1 ].setFromNormalAndCoplanarPoint( _norm, _worldMax ).negate();
+
+		_norm.set( 0, 1, 0 ).transformDirection( this.transform );
+		this.planes[ 2 ].setFromNormalAndCoplanarPoint( _norm, _worldMin );
+		this.planes[ 3 ].setFromNormalAndCoplanarPoint( _norm, _worldMax ).negate();
+
+		_norm.set( 1, 0, 0 ).transformDirection( this.transform );
+		this.planes[ 4 ].setFromNormalAndCoplanarPoint( _norm, _worldMin );
+		this.planes[ 5 ].setFromNormalAndCoplanarPoint( _norm, _worldMax ).negate();
+
 	}
 
 	// based on three.js' Box3 "intersects frustum" function
@@ -64,6 +90,27 @@ export class OBB {
 
 		}
 
+		// do the opposite check using the obb planes to avoid false positives
+		for ( let i = 0; i < 6; i ++ ) {
+
+			const plane = this.planes[ i ];
+			let maxDistance = - Infinity;
+			for ( let j = 0; j < 8; j ++ ) {
+
+				const v = frustum.points[ j ];
+				const dist = plane.distanceToPoint( v );
+				maxDistance = maxDistance < dist ? dist : maxDistance;
+
+			}
+
+			if ( maxDistance < 0 ) {
+
+				return false;
+
+			}
+
+		}
+
 		return true;
 
 	}


Test patch:
diff --git a/test/Ellipsoid.test.js b/test/Ellipsoid.test.js
index 040e659..f6a8015 100644
--- a/test/Ellipsoid.test.js
+++ b/test/Ellipsoid.test.js
@@ -178,6 +178,37 @@ describe( 'Ellipsoid', () => {
 
 	} );
 
+	it( 'should match the expected elevation.', () => {
+
+		//ellipsoid rotation
+		const northPos = new Vector3( 0, 0, WGS84_HEIGHT );
+		expect( wgsEllipse.getPositionElevation( northPos ) ).toBeCloseTo( 0, 1e-6 );
+		const southPos = new Vector3( 0, 0, - WGS84_HEIGHT );
+		expect( wgsEllipse.getPositionElevation( southPos ) ).toBeCloseTo( 0, 1e-6 );
+		const xPos = new Vector3( WGS84_RADIUS, 0, 0 );
+		expect( wgsEllipse.getPositionElevation( xPos ) ).toBeCloseTo( 0, 1e-6 );
+		const mxPos = new Vector3( - WGS84_RADIUS, 0, 0 );
+		expect( wgsEllipse.getPositionElevation( mxPos ) ).toBeCloseTo( 0, 1e-6 );
+		const zPos = new Vector3( 0, WGS84_RADIUS, 0 );
+		expect( wgsEllipse.getPositionElevation( zPos ) ).toBeCloseTo( 0, 1e-6 );
+		const mzPos = new Vector3( 0, - WGS84_RADIUS, 0 );
+		expect( wgsEllipse.getPositionElevation( mzPos ) ).toBeCloseTo( 0, 1e-6 );
+
+		const northPos100 = new Vector3( 0, 0, WGS84_HEIGHT + 100 );
+		expect( wgsEllipse.getPositionElevation( northPos100 ) ).toBeCloseTo( 100, 1e-6 );
+		const southPos100 = new Vector3( 0, 0, - WGS84_HEIGHT + 100 );
+		expect( wgsEllipse.getPositionElevation( southPos100 ) ).toBeCloseTo( 100, 1e-6 );
+		const xPos100 = new Vector3( WGS84_RADIUS + 100, 0, 0 );
+		expect( wgsEllipse.getPositionElevation( xPos100 ) ).toBeCloseTo( 100, 1e-6 );
+		const mxPos100 = new Vector3( - WGS84_RADIUS + 100, 0, 0 );
+		expect( wgsEllipse.getPositionElevation( mxPos100 ) ).toBeCloseTo( 100, 1e-6 );
+		const zPos100 = new Vector3( 0, WGS84_RADIUS + 100, 0 );
+		expect( wgsEllipse.getPositionElevation( zPos100 ) ).toBeCloseTo( 100, 1e-6 );
+		const mzPos100 = new Vector3( 0, - WGS84_RADIUS + 100, 0 );
+		expect( wgsEllipse.getPositionElevation( mzPos100 ) ).toBeCloseTo( 100, 1e-6 );
+
+	} );
+
 } );
 
 describe( 'EllipsoidRegion', () => {