chore(pointcloud): Add adaptive size mode

Author: ketourneau

packages/Debug/src/PointCloudDebug.js

@@ -276,7 +276,7 @@ export default {
         styleUI.add(layer, 'opacity', 0, 1).name('Layer opacity').onChange(update);
         styleUI.add(layer, 'pointSize', 0, 15).name('Point size').onChange(update);
         if (layer.material.sizeMode != undefined && view.camera.camera3D.isPerspectiveCamera) {
-            styleUI.add(layer.material, 'sizeAttenuation').name('Size attenuation')
+            styleUI.add(layer.material, 'sizeMode', PNTS_SIZE_MODE).name('Size mode')
                 .onChange(update);
             styleUI.add(layer.material, 'minAttenuatedSize', 0, 15).name('Min size')
                 .onChange((value) => {

packages/Main/src/Layer/PointCloudLayer.ts

@@ -1,6 +1,6 @@
 import * as THREE from 'three';
 import GeometryLayer from 'Layer/GeometryLayer';
-import PointsMaterial, { PNTS_MODE } from 'Renderer/PointsMaterial';
+import PointsMaterial, { PNTS_MODE, PNTS_SIZE_MODE } from 'Renderer/PointsMaterial';
 import Picking from 'Core/Picking';
 
 import type PointCloudNode from 'Core/PointCloudNode';
@@ -506,6 +506,58 @@ abstract class PointCloudLayer<S extends PointCloudSource = PointCloudSource>
                 this.dispatchEvent({ type: 'dispose-model', scene: obj, tile: obj.userData.node });
             }
         }
+
+        // @ts-expect-error PointsMaterial is not typed yet
+        if (this.material.sizeMode === PNTS_SIZE_MODE.ADAPTIVE) {
+            const nodes = this.getNodes(this.group.children);
+            const visibilityTextureData = this.computeVisibilityTextureData(nodes);
+
+            // @ts-expect-error PointsMaterial is not typed yet
+            const vnt = this.material.visibleNodes;
+            const data = vnt.image.data;
+            data.set(visibilityTextureData.data);
+            vnt.needsUpdate = true;
+
+            // Use natBox (native octree space) for octreeSize, as the
+            // octree hierarchy is defined in the source CRS coordinate system
+            const rootSize = this.root!.voxelOBB.natBox.getSize(new THREE.Vector3());
+            const octreeSize = Math.max(rootSize.x, rootSize.y, rootSize.z);
+
+            for (const pts of this.group.children) {
+                const node = pts.userData.node;
+                const depth = node.depth;
+                const nodeStartOffset = visibilityTextureData.offsets.get(node.voxelKey);
+                const octreeSpacing = node.source.spacing;
+
+                // Compute the bounding box min of the node's octree cell
+                // in the local space of the THREE.Points geometry.
+                // We must use voxelOBB.natBox.min (the octree cell boundary
+                // in source CRS) rather than geomBBox.min (which only covers
+                // the actual points, not the full octree cell).
+                const natMin = node.voxelOBB.natBox.min;
+                const origin = node.origin;
+                const bboxMin = new THREE.Vector3(
+                    natMin.x - origin.x,
+                    natMin.y - origin.y,
+                    natMin.z - origin.z,
+                );
+
+                pts.onBeforeRender = (_renderer, _scene, _camera, _geometry, material) => {
+                    // @ts-expect-error Material is not typed yet
+                    material.uniforms.nodeStartOffset.value = nodeStartOffset;
+                    // @ts-expect-error Material is not typed yet
+                    material.uniforms.octreeSize.value = octreeSize;
+                    // @ts-expect-error Material is not typed yet
+                    material.uniforms.octreeSpacing.value = octreeSpacing;
+                    // @ts-expect-error Material is not typed yet
+                    material.uniforms.nodeDepth.value = depth;
+                    // @ts-expect-error Material is not typed yet
+                    material.uniforms.nodeBBoxMin.value.copy(bboxMin);
+                    // @ts-expect-error Material is not typed yet
+                    material.uniformsNeedUpdate = true;
+                };
+            }
+        }
     }
 
     // @ts-expect-error Layer and Picking are not typed yet
@@ -529,6 +581,99 @@ abstract class PointCloudLayer<S extends PointCloudSource = PointCloudSource>
             }
         }
     }
+
+    getNodes(children: THREE.Object3D[]) {
+        const nodes = new Set<PointCloudNode>();
+        // Sometimes child node is loaded before parent node, so we need to
+        // make sure to add all the parent nodes of the visible nodes
+        const collectNodeParentRecursively = (node: PointCloudNode): void => {
+            if (nodes.has(node)) {
+                return;
+            }
+            nodes.add(node);
+            if (node.parent) {
+                collectNodeParentRecursively(node.parent);
+            }
+        };
+        children.forEach(child => collectNodeParentRecursively(child.userData.node));
+        return Array.from(nodes);
+    }
+
+    // Encoding the octree hierarchy in breadth-first order
+    // into a texture for adaptive point size rendering
+    // Explanation p36: https://www.cg.tuwien.ac.at/research/publications/2016/SCHUETZ-2016-POT/SCHUETZ-2016-POT-thesis.pdf
+    computeVisibilityTextureData(nodes: PointCloudNode[]) {
+        // sort by level and hierarchy order
+        const sort = function sortNodes(a: PointCloudNode, b: PointCloudNode) {
+            if (a.depth !== b.depth) { return a.depth - b.depth; }
+            // @ts-expect-error PointCloudNode has x properties
+            if (a.x !== b.x) { return a.x - b.x; }
+            // @ts-expect-error PointCloudNode has y properties
+            if (a.y !== b.y) { return a.y - b.y; }
+            // @ts-expect-error PointCloudNode has z properties
+            return a.z - b.z;
+        };
+        // breadth-first order
+        const orderedNodes = nodes.toSorted(sort);
+
+        const data = new Uint8Array(orderedNodes.length * 4);
+        const visibleNodeTextureOffsets = new Map<string, number>();
+        const offsetsToChild: number[] = new Array(orderedNodes.length).fill(Infinity);
+
+        // Helper function to get octree child index from node
+        const getChildIndex = (node: PointCloudNode): number => {
+            if (!node.parent) {
+                return 0;
+            }
+            const parent = node.parent;
+            // @ts-expect-error PointCloudNode has x properties
+            const dx = node.x - parent.x * 2;
+            // @ts-expect-error PointCloudNode has y properties
+            const dy = node.y - parent.y * 2;
+            // @ts-expect-error PointCloudNode has z properties
+            const dz = node.z - parent.z * 2;
+            // Octree child index (Potree convention): 4*x + 2*y + z
+            return 4 * dx + 2 * dy + dz;
+        };
+
+        for (let nodeIndex = 0; nodeIndex < orderedNodes.length; nodeIndex++) {
+            const node = orderedNodes[nodeIndex];
+            // @ts-expect-error PointCloudNode has voxelKey properties
+            visibleNodeTextureOffsets.set(node.voxelKey, nodeIndex);
+
+            if (node.parent) {
+                const childIndex = getChildIndex(node);
+                // @ts-expect-error PointCloudNode has voxelKey properties
+                const parentIndex = visibleNodeTextureOffsets.get(node.parent.voxelKey);
+
+                if (parentIndex === undefined) {
+                    continue;
+                }
+
+                const parentOffsetToChild = nodeIndex - parentIndex;
+                const offsetToFirstChild =
+                    Math.min(offsetsToChild[parentIndex], parentOffsetToChild);
+                offsetsToChild[parentIndex] = offsetToFirstChild;
+
+                // The 8 bits of the red value indicate
+                // which of the children are visible
+                data[parentIndex * 4] = data[parentIndex * 4] | (1 << childIndex);
+                // Offset to child is stored on 2 bytes,
+                // so it can support up to 65536 nodes per subtree.
+                // The green channel contains the relative offset
+                // to the node’s first child (8 most significant bits)
+                data[parentIndex * 4 + 1] = offsetToFirstChild >> 8;
+                // The blue channel contains the relative offset
+                // to the node’s first child (8 least significant bits)
+                data[parentIndex * 4 + 2] = offsetToFirstChild % 256;
+            }
+        }
+
+        return {
+            data,
+            offsets: visibleNodeTextureOffsets,
+        };
+    }
 }
 
 export default PointCloudLayer;

packages/Main/src/Renderer/PointsMaterial.js

@@ -25,6 +25,7 @@ export const PNTS_SHAPE = {
 export const PNTS_SIZE_MODE = {
     VALUE: 0,
     ATTENUATED: 1,
+    ADAPTIVE: 2,
 };
 
 const white = new THREE.Color(1.0,  1.0,  1.0);
@@ -245,6 +246,13 @@ class PointsMaterial extends THREE.ShaderMaterial {
         CommonMaterial.setUniformProperty(this, 'gamma', gamma);
         CommonMaterial.setUniformProperty(this, 'ambientBoost', ambientBoost);
 
+        // Adaptive point size uniforms
+        CommonMaterial.setUniformProperty(this, 'octreeSpacing', 1.0);
+        CommonMaterial.setUniformProperty(this, 'octreeSize', 1.0);
+        CommonMaterial.setUniformProperty(this, 'nodeDepth', 0.0);
+        CommonMaterial.setUniformProperty(this, 'nodeStartOffset', 0.0);
+        CommonMaterial.setUniformProperty(this, 'nodeBBoxMin', new THREE.Vector3());
+
         // add classification texture to apply classification lut.
         const data = new Uint8Array(256 * 4);
         const texture = new THREE.DataTexture(data, 256, 1, THREE.RGBAFormat);
@@ -267,6 +275,12 @@ class PointsMaterial extends THREE.ShaderMaterial {
         textureVisi.magFilter = THREE.NearestFilter;
         CommonMaterial.setUniformProperty(this, 'visibilityTexture', textureVisi);
 
+        const dataNodes = new Uint8Array(2048 * 4);
+        const visibleNodesTexture = new THREE.DataTexture(dataNodes, 2048, 1, THREE.RGBAFormat);
+        visibleNodesTexture.needsUpdate = true;
+        visibleNodesTexture.magFilter = THREE.NearestFilter;
+        CommonMaterial.setUniformProperty(this, 'visibleNodes', visibleNodesTexture);
+
         // Classification and other discrete values scheme
         this.classificationScheme = classificationScheme;
         this.discreteScheme = discreteScheme;

packages/Main/src/Renderer/Shader/PointsVS.glsl

@@ -31,6 +31,14 @@ uniform int sizeMode;
 uniform float minAttenuatedSize;
 uniform float maxAttenuatedSize;
 
+// Adaptive point size uniforms
+uniform sampler2D visibleNodes;
+uniform float octreeSpacing;
+uniform float octreeSize;
+uniform float nodeDepth;
+uniform float nodeStartOffset;
+uniform vec3 nodeBBoxMin;
+
 attribute vec4 unique_id;
 attribute float intensity;
 attribute float classification;
@@ -40,6 +48,61 @@ attribute float returnNumber;
 attribute float numberOfReturns;
 attribute float scanAngle;
 
+// Adaptive point size calculation functions (from Potree)
+/**
+ * number of 1-bits up to inclusive index position
+ * number is treated as if it were an integer in the range 0-255
+ *
+ */
+int numberOfOnes(int number, int index) {
+    int numOnes = 0;
+    int tmp = 128;
+    for (int i = 7; i >= 0; i--) {
+        if (number >= tmp) {
+            number = number - tmp;
+            if (i <= index) {
+                numOnes++;
+            }
+        }
+        tmp = tmp / 2;
+    }
+    return numOnes;
+}
+
+float getLOD() {
+    // Transform position from local space (relative to node origin)
+    // to octree space (relative to node bbox min)
+    vec3 pos = position - nodeBBoxMin;
+    vec3 offset = vec3(0.0, 0.0, 0.0);
+    int iOffset = int(nodeStartOffset);
+    float depth = nodeDepth;
+    for (float i = 0.0; i <= 30.0; i++) {
+        float nodeSizeAtLevel = octreeSize / pow(2.0, i + nodeDepth);
+
+        vec3 index3d = (pos - offset) / nodeSizeAtLevel;
+        index3d = floor(index3d + 0.5);
+        int index = int(round(4.0 * index3d.x + 2.0 * index3d.y + index3d.z));
+
+        vec4 value = texture2D(visibleNodes, vec2(float(iOffset) / 2048.0, 0.0));
+        int mask = int(round(value.r * 255.0));
+        bool childNodeExist = bool(((mask >> index) & 1) != 0);
+
+        if (childNodeExist) {
+            int greenChannelOffset = int(round(value.g * 255.0)) * 256;
+            int blueChannelOffset = int(round(value.b * 255.0));
+            int childIndexOffset = numberOfOnes(mask, index - 1);
+            int totalOffset = greenChannelOffset + blueChannelOffset + childIndexOffset;
+            iOffset = iOffset + totalOffset;
+            depth++;
+        } else {
+            // no more visible child nodes at this position
+            return depth;
+        }
+        offset = offset + (vec3(1.0, 1.0, 1.0) * nodeSizeAtLevel * 0.5) * index3d;
+    }
+    return depth;
+}
+
 void main() {
     vec2 uv = vec2(classification/255., 0.5);
 
@@ -117,13 +180,21 @@ void main() {
 
     gl_PointSize = size;
 
-    if (sizeMode == PNTS_SIZE_MODE_ATTENUATED) {
-        bool isPerspective = isPerspectiveMatrix(projectionMatrix);
+    bool isPerspective = isPerspectiveMatrix(projectionMatrix);
+    float depthAttenuationFactor = scale / -mvPosition.z;
 
+    if (sizeMode == PNTS_SIZE_MODE_ATTENUATED) {
         if (isPerspective) {
-            gl_PointSize *= scale / -mvPosition.z;
+            gl_PointSize *= depthAttenuationFactor;
             gl_PointSize = clamp(gl_PointSize, minAttenuatedSize, maxAttenuatedSize);
         }
+    } else if (sizeMode == PNTS_SIZE_MODE_ADAPTIVE) {
+        if (isPerspective) {
+            float r = octreeSpacing * 1.7;
+            float pointSizeAttenuation = pow(2.0, getLOD());
+            float worldSpaceSize = size * r / pointSizeAttenuation;
+            gl_PointSize = worldSpaceSize * depthAttenuationFactor;
+        }
     }
 
 #include <logdepthbuf_vertex>