GlobeControls limit Frustum with horizon for performance (#483)

* GlobeControls limit Frustum with horizon

* add custom Frustum

* fix obb intersectsFrustum false positive

* revert mistake change

* fix obb planes

* document new ellippsoid methods

* make obb more readable

* Ellipsoid fix getPositionElevation

* add comment to explain horizon edge cases

* clean PR

* Small style changes

---------

Co-authored-by: Garrett Johnson <[email protected]>

example/googleMapsExample.js

@@ -146,25 +146,29 @@ 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 ) }` );
 
 }
 
 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 );
 
 }

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(),

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();
 
 	}

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;
+
+	}
+
 }

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 };

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/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', () => {