refactor(core): unify SSE computation methods

Until this commit each geometry type had its own SSE computation methods.
The goal of this commit is to provide a single method based on:
  - bounding-box (BB)
  - geometricError (in meters)
  - distance

The computation is simple:
  - transform camera in BB's basis
  - compute distance camera - BB
  - project geometricError on screen at distance

The last step produce a screen-space-error in pixels.

Geometric error values depend on the data type:
  - 3dtiles: the tileset defines a geometricError per node so we can use this value
  - pointcloud: PotreeConverter defines a spacing, which is the average distance
between points in each node. We use: geometricError = spacing / 2^depth
  - tiles: this commit adapted the existing formula to use consistent units

src/Core/Prefab/GlobeView.js

@@ -159,7 +159,6 @@ export function createGlobeLayer(id, options) {
         if (SubdivisionControl.hasEnoughTexturesToSubdivide(context, layer, node)) {
             return globeSubdivisionControl(2,
                 options.maxSubdivisionLevel || 18,
-                options.sseSubdivisionThreshold || 1.0,
                 options.maxDeltaElevationLevel || 4)(context, layer, node);
         }
         return false;

src/Core/Scheduler/Providers/PointCloudProvider.js

@@ -83,7 +83,15 @@ function parseOctree(layer, hierarchyStepSize, root) {
                         const myname = childname.substr(root.name.length);
                         url = `${root.baseurl}/${myname}`;
                     }
-                    const item = { numPoints: n, childrenBitField: c, children: [], name: childname, baseurl: url, bbox: bounds };
+                    const item = {
+                        numPoints: n,
+                        childrenBitField: c,
+                        children: [],
+                        name: childname,
+                        baseurl: url,
+                        bbox: bounds,
+                        geometricError: layer.metadata.spacing / Math.pow(2, childname.length),
+                    };
                     snode.children.push(item);
                     stack.push(item);
                 }
@@ -153,7 +161,7 @@ export default {
         layer.octreeDepthLimit = layer.octreeDepthLimit || -1;
         layer.pointBudget = layer.pointBudget || 15000000;
         layer.pointSize = layer.pointSize === 0 || !isNaN(layer.pointSize) ? layer.pointSize : 4;
-        layer.overdraw = layer.overdraw || 2;
+        layer.sseThreshold = layer.sseThreshold || 5;
         layer.type = 'geometry';
 
         // default update methods
@@ -202,7 +210,12 @@ export default {
             return parseOctree(
                     layer,
                     layer.metadata.hierarchyStepSize,
-                    { baseurl: `${layer.url}/${cloud.octreeDir}/r`, name: '', bbox });
+                {
+                    baseurl: `${layer.url}/${cloud.octreeDir}/r`,
+                    name: '',
+                    bbox,
+                    geometricError: layer.metadata.spacing,
+                });
         }).then((root) => {
             // eslint-disable-next-line no-console
             console.log('LAYER metadata:', root);

src/Core/Scheduler/Providers/TileProvider.js

@@ -18,6 +18,7 @@ function preprocessDataLayer(layer, view, scheduler) {
 
     layer.level0Nodes = [];
     layer.onTileCreated = layer.onTileCreated || (() => {});
+    layer.sseThreshold = layer.sseThreshold || 5.0;
 
     const promises = [];
 

src/Core/ScreenSpaceError.js

@@ -0,0 +1,88 @@
+import * as THREE from 'three';
+
+const v = new THREE.Vector3();
+const m = new THREE.Matrix4();
+const m2 = new THREE.Matrix4();
+const m3 = new THREE.Matrix4();
+
+function compute(vector, matrix, camera, distance, _3d) {
+    const basis = [
+        new THREE.Vector3(vector.x, 0, 0),
+        new THREE.Vector3(0, vector.y, 0),
+    ];
+
+    if (_3d) {
+        basis.push(new THREE.Vector3(0, 0, vector.z));
+    }
+
+    m2.identity();
+    m2.extractRotation(matrix);
+    m3.identity();
+    m3.extractRotation(camera.camera3D.matrixWorldInverse);
+
+    for (const b of basis) {
+        // Apply rotation
+        b.applyMatrix4(m2);
+        // Apply inverse camera rotation
+        b.applyMatrix4(m3);
+        // Move at 'distance' from camera
+        b.z += -distance;
+        // Project on screen
+        b.applyMatrix4(camera.camera3D.projectionMatrix);
+        // cancel z component
+        b.z = 0;
+        b.x = b.x * camera.width * 0.5;
+        b.y = b.y * camera.height * 0.5;
+    }
+
+    const lengthsq = basis.map(b => b.lengthSq());
+    const min = Math.min.apply(Math, lengthsq);
+    return Math.sqrt(min);
+}
+
+function findBox3Distance(camera, box3, matrix) {
+    // TODO: can be cached
+    m.getInverse(matrix);
+    // Move camera position in box3 basis
+    // (we don't transform box3 to camera basis because box3 are AABB)
+    const pt = camera.camera3D.position
+        .clone(v).applyMatrix4(m);
+    // Compute distance between the camera / box3
+    return box3.distanceToPoint(pt);
+}
+
+function computeSizeFromGeometricError(box3, geometricError) {
+    const size = box3.getSize();
+    // Build a vector with the same ratio than box3,
+    // and with the biggest component being geometricError
+    size.multiplyScalar(geometricError /
+        Math.max(size.x, Math.max(size.y, size.z)));
+    return size;
+}
+
+export default {
+    MODE_2D: 1,
+
+    MODE_3D: 2,
+
+    computeFromBox3(camera, box3, matrix, geometricError, mode) {
+        const distance = findBox3Distance(camera, box3, matrix);
+
+        if (distance <= geometricError) {
+            return {
+                sse: Infinity,
+                distance,
+            };
+        }
+
+        const size = computeSizeFromGeometricError(box3, geometricError);
+
+        const sse = compute(size, matrix, camera, distance, mode == this.MODE_3D);
+
+        return {
+            sse,
+            distance,
+            size,
+        };
+    },
+};

src/Core/TileMesh.js

@@ -8,7 +8,7 @@ import * as THREE from 'three';
 import LayeredMaterial from '../Renderer/LayeredMaterial';
 import { l_ELEVATION } from '../Renderer/LayeredMaterialConstants';
 import RendererConstant from '../Renderer/RendererConstant';
-import OGCWebServiceHelper, { SIZE_TEXTURE_TILE } from './Scheduler/Providers/OGCWebServiceHelper';
+import OGCWebServiceHelper from './Scheduler/Providers/OGCWebServiceHelper';
 
 function TileMesh(geometry, params) {
     // Constructor
@@ -142,9 +142,7 @@ TileMesh.prototype.setBBoxZ = function setBBoxZ(min, max) {
 };
 
 TileMesh.prototype.updateGeometricError = function updateGeometricError() {
-    // The geometric error is calculated to have a correct texture display.
-    // For the projection of a texture's texel to be less than or equal to one pixel
-    this.geometricError = this.boundingSphere.radius / SIZE_TEXTURE_TILE;
+    this.geometricError = this.boundingSphere.radius;
 };
 
 TileMesh.prototype.setTexturesLayer = function setTexturesLayer(textures, layerType, layerId) {

src/Main.js

@@ -24,3 +24,4 @@ export { default as GeoJSON2Features } from './Renderer/ThreeExtended/GeoJSON2Fe
 export { default as FeaturesUtils } from './Renderer/ThreeExtended/FeaturesUtils';
 export { CONTROL_EVENTS } from './Renderer/ThreeExtended/GlobeControls';
 export { default as DEMUtils } from './utils/DEMUtils';
+export { default as ScreenSpaceError } from './Core/ScreenSpaceError';

src/Process/3dTilesProcessing.js

@@ -1,4 +1,5 @@
 import * as THREE from 'three';
+import ScreenSpaceError from '../Core/ScreenSpaceError';
 
 function requestNewTile(view, scheduler, geometryLayer, metadata, parent, redraw) {
     const command = {
@@ -225,25 +226,30 @@ const worldPosition = new THREE.Vector3();
 function computeNodeSSE(camera, node) {
     node.distance = 0;
     if (node.boundingVolume.region) {
-        worldPosition.setFromMatrixPosition(node.boundingVolume.region.matrixWorld);
-        cameraLocalPosition.copy(camera.camera3D.position).sub(worldPosition);
-        node.distance = node.boundingVolume.region.box3D.distanceToPoint(cameraLocalPosition);
-    } else if (node.boundingVolume.box) {
-        worldPosition.setFromMatrixPosition(node.matrixWorld);
-        cameraLocalPosition.copy(camera.camera3D.position).sub(worldPosition);
-        node.distance = node.boundingVolume.box.distanceToPoint(cameraLocalPosition);
-    } else if (node.boundingVolume.sphere) {
+        return ScreenSpaceError.computeFromBox3(
+            camera,
+            node.boundingVolume.region.box3D,
+            node.boundingVolume.region.matrixWorld,
+            node.geometricError,
+            ScreenSpaceError.MODE_2D);
+    }
+    if (node.boundingVolume.box) {
+        return ScreenSpaceError.computeFromBox3(
+            camera,
+            node.boundingVolume.box,
+            node.matrixWorld,
+            node.geometricError,
+            ScreenSpaceError.MODE_3D);
+    }
+    if (node.boundingVolume.sphere) {
+        // TODO USE http://iquilezles.org/www/articles/sphereproj/sphereproj.htm
         worldPosition.setFromMatrixPosition(node.matrixWorld);
         cameraLocalPosition.copy(camera.camera3D.position).sub(worldPosition);
-        node.distance = Math.max(0.0, node.boundingVolume.sphere.distanceToPoint(cameraLocalPosition));
+        const distance = Math.max(0.0, node.boundingVolume.sphere.distanceToPoint(cameraLocalPosition));
+        return camera.preSSE * (node.geometricError / distance);
     } else {
         return Infinity;
     }
-    if (node.distance === 0) {
-        // This test is needed in case geometricError = distance = 0
-        return Infinity;
-    }
-    return camera.preSSE * (node.geometricError / node.distance);
 }
 
 export function init3dTilesLayer(view, scheduler, layer) {
@@ -312,6 +318,8 @@ export function process3dTilesNode(cullingTest, subdivisionTest) {
                 });
             }
             return returnValue;
+        } else {
+            node.sse = Infinity;
         }
 
         markForDeletion(layer, node);
@@ -324,6 +332,6 @@ export function $3dTilesSubdivisionControl(context, layer, node) {
     if (layer.tileIndex.index[node.tileId].children === undefined) {
         return false;
     }
-    const sse = computeNodeSSE(context.camera, node);
-    return sse > layer.sseThreshold;
+    node.sse = computeNodeSSE(context.camera, node);
+    return node.sse.sse > layer.sseThreshold;
 }

src/Process/GlobeTileProcessing.js

@@ -3,6 +3,7 @@ import MathExt from '../Core/Math/MathExtended';
 import { UNIT, ellipsoidSizes } from '../Core/Geographic/Coordinates';
 import { SIZE_TEXTURE_TILE } from '../Core/Scheduler/Providers/OGCWebServiceHelper';
 import Extent from '../Core/Geographic/Extent';
+import ScreenSpaceError from '../Core/ScreenSpaceError';
 
 const cV = new THREE.Vector3();
 let vhMagnitudeSquared;
@@ -11,18 +12,6 @@ let SSE_SUBDIVISION_THRESHOLD;
 
 const worldToScaledEllipsoid = new THREE.Matrix4();
 
-function _preSSE(view) {
-    const canvasSize = view.mainLoop.gfxEngine.getWindowSize();
-    const hypotenuse = canvasSize.length();
-    const radAngle = view.camera.camera3D.fov * Math.PI / 180;
-
-     // TODO: not correct -> see new preSSE
-    // const HFOV = 2.0 * Math.atan(Math.tan(radAngle * 0.5) / context.camera.ratio);
-    const HYFOV = 2.0 * Math.atan(Math.tan(radAngle * 0.5) * hypotenuse / canvasSize.x);
-
-    return hypotenuse * (2.0 * Math.tan(HYFOV * 0.5));
-}
-
 export function preGlobeUpdate(context, layer) {
     // We're going to use the method described here:
     //    https://cesiumjs.org/2013/04/25/Horizon-culling/
@@ -40,9 +29,6 @@ export function preGlobeUpdate(context, layer) {
     cV.copy(context.camera.camera3D.position).applyMatrix4(worldToScaledEllipsoid);
     vhMagnitudeSquared = cV.lengthSq() - 1.0;
 
-    // pre-sse
-    context.camera.preSSE = _preSSE(context.view);
-
     const elevationLayers = context.view.getLayers((l, a) => a && a.id == layer.id && l.type == 'elevation');
     context.maxElevationLevel = -1;
     for (const e of elevationLayers) {
@@ -92,30 +78,7 @@ export function globeCulling(minLevelForHorizonCulling) {
     };
 }
 
-const v = new THREE.Vector3();
-function computeNodeSSE(camera, node) {
-    v.setFromMatrixScale(node.matrixWorld);
-    const boundingSphereCenter = node.boundingSphere.center.clone().applyMatrix4(node.matrixWorld);
-    const distance = Math.max(
-        0.0,
-        camera.camera3D.position.distanceTo(boundingSphereCenter) - node.boundingSphere.radius * v.x);
-
-    // Removed because is false computation, it doesn't consider the altitude of node
-    // Added small oblique weight (distance is not enough, tile orientation is needed)
-    /*
-    var altiW = node.bbox.top() === 10000 ? 0. : node.bbox.bottom() / 10000.;
-    var dotProductW = Math.min(altiW + Math.abs(this.camera3D.getWorldDirection().dot(node.centerSphere.clone().normalize())), 1.);
-    if (this.camera3D.position.length() > 6463300) dotProductW = 1;
-    var SSE = Math.sqrt(dotProductW) * this.preSSE * (node.geometricError / distance);
-    */
-
-    // TODO: node.geometricError is computed using a hardcoded 18 level
-    // The computation of node.geometricError is surely false
-    return camera.preSSE * (node.geometricError * v.x) / distance;
-}
-
-export function globeSubdivisionControl(minLevel, maxLevel, sseThreshold, maxDeltaElevationLevel) {
-    SSE_SUBDIVISION_THRESHOLD = sseThreshold;
+export function globeSubdivisionControl(minLevel, maxLevel, maxDeltaElevationLevel) {
     return function _globeSubdivisionControl(context, layer, node) {
         if (node.level < minLevel) {
             return true;
@@ -133,9 +96,14 @@ export function globeSubdivisionControl(minLevel, maxLevel, sseThreshold, maxDel
             return false;
         }
 
-        const sse = computeNodeSSE(context.camera, node);
-
-        return SSE_SUBDIVISION_THRESHOLD < sse;
+        node.sse = ScreenSpaceError.computeFromBox3(
+            context.camera,
+            node.OBB().box3D,
+            node.OBB().matrixWorld,
+            node.geometricError,
+            ScreenSpaceError.MODE_2D);
+        node.sse.offset = SIZE_TEXTURE_TILE;
+        return node.sse.sse > (SIZE_TEXTURE_TILE + layer.sseThreshold);
     };
 }
 
@@ -164,6 +132,7 @@ export function globeSchemeTileWMTS(type) {
 }
 
 export function computeTileZoomFromDistanceCamera(distance, view) {
+    // TODO fixme
     const sizeEllipsoid = ellipsoidSizes().x;
     const preSinus = SIZE_TEXTURE_TILE * (SSE_SUBDIVISION_THRESHOLD * 0.5) / view.camera.preSSE / sizeEllipsoid;
 
@@ -182,6 +151,7 @@ export function computeTileZoomFromDistanceCamera(distance, view) {
 }
 
 export function computeDistanceCameraFromTileZoom(zoom, view) {
+    // TODO fixme
     const delta = Math.PI / Math.pow(2, zoom);
     const circleChord = 2.0 * ellipsoidSizes().x * Math.sin(delta * 0.5);
     const radius = circleChord * 0.5;

src/Process/PlanarTileProcessing.js

@@ -1,3 +1,6 @@
+import ScreenSpaceError from '../Core/ScreenSpaceError';
+import { SIZE_TEXTURE_TILE } from '../Core/Scheduler/Providers/OGCWebServiceHelper';
+
 function frustumCullingOBB(node, camera) {
     return camera.isBox3Visible(node.OBB().box3D, node.OBB().matrixWorld);
 }
@@ -6,22 +9,6 @@ export function planarCulling(node, camera) {
     return !frustumCullingOBB(node, camera);
 }
 
-function _isTileBigOnScreen(camera, node) {
-    const onScreen = camera.box3SizeOnScreen(node.OBB().box3D, node.matrixWorld);
-
-    // onScreen.x/y/z are [-1, 1] so divide by 2
-    // (so x = 0.4 means the object width on screen is 40% of the total screen width)
-    const dim = {
-        x: 0.5 * (onScreen.max.x - onScreen.min.x),
-        y: 0.5 * (onScreen.max.y - onScreen.min.y),
-    };
-
-    // subdivide if on-screen width (and resp. height) is bigger than 30% of the screen width (resp. height)
-    // TODO: the 30% value is arbitrary and needs to be configurable by the user
-    // TODO: we might want to use texture resolution here as well
-    return (dim.x >= 0.3 && dim.y >= 0.3);
-}
-
 export function prePlanarUpdate(context, layer) {
     const elevationLayers = context.view.getLayers((l, a) => a && a.id == layer.id && l.type == 'elevation');
     context.maxElevationLevel = -1;
@@ -48,7 +35,14 @@ export function planarSubdivisionControl(maxLevel, maxDeltaElevationLevel) {
             (node.level - currentElevationLevel) >= maxDeltaElevationLevel) {
             return false;
         }
-
-        return _isTileBigOnScreen(context.camera, node);
+        node.sse = ScreenSpaceError.computeFromBox3(
+            context.camera,
+            node.OBB().box3D,
+            node.OBB().matrixWorld,
+            node.geometricError,
+            ScreenSpaceError.MODE_2D);
+        node.sse.sse = Math.max(0, node.sse.sse - SIZE_TEXTURE_TILE);
+        node.sse.offset = SIZE_TEXTURE_TILE;
+        return node.sse.sse > layer.sseThreshold;
     };
 }