refactor(LASLoader): reproj data during parsing and add elevation att…

…ributs

src/Loader/LASLoader.js

@@ -1,5 +1,6 @@
 import { LazPerf } from 'laz-perf';
 import { Las } from 'copc';
+import proj4 from 'proj4';
 
 /**
  * @typedef {Object} Header - Partial LAS header.
@@ -49,6 +50,11 @@ class LASLoader {
     _parseView(view, options) {
         const colorDepth = options.colorDepth ?? 16;
 
+        const forward = (options.crsIn !== options.crsOut) ?
+            proj4(options.projDefs[options.crsIn], options.projDefs[options.crsOut]).forward :
+            (x => x);
+        const isGeocentric = options.projDefs[options.crsOut].projName === 'geocent';
+
         const getPosition = ['X', 'Y', 'Z'].map(view.getter);
         const getIntensity = view.getter('Intensity');
         const getReturnNumber = view.getter('ReturnNumber');
@@ -60,6 +66,7 @@ class LASLoader {
         const getScanAngle = view.getter('ScanAngle');
 
         const positions = new Float32Array(view.pointCount * 3);
+        const elevations = new Float32Array(view.pointCount);
         const intensities = new Uint16Array(view.pointCount);
         const returnNumbers = new Uint8Array(view.pointCount);
         const numberOfReturns = new Uint8Array(view.pointCount);
@@ -75,17 +82,23 @@ class LASLoader {
         */
         const scanAngles = new Float32Array(view.pointCount);
 
-        // For precision we take the first point that will be use as origin for a local referentiel.
-        const origin = getPosition.map(f => f(0)).map(val => Math.floor(val));
+        // For precision we use the first point to define the origin for a local referentiel.
+        // After projection transformation and only the integer part for simplification.
+        const origin = forward(getPosition.map(f => f(0))).map(val => Math.floor(val));
 
         for (let i = 0; i < view.pointCount; i++) {
             // `getPosition` apply scale and offset transform to the X, Y, Z
             // values. See https://github.com/connormanning/copc.js/blob/master/src/las/extractor.ts.
-            const [x, y, z] = getPosition.map(f => f(i));
+            // we thus apply the projection to get values in the Crs of the view.
+            const point = getPosition.map(f => f(i));
+            const [x, y, z] = forward(point);
             positions[i * 3] = x - origin[0];
             positions[i * 3 + 1] = y - origin[1];
             positions[i * 3 + 2] = z - origin[2];
 
+            elevations[i] = z;
+            // geocentric height to elevation
+            if (isGeocentric) { elevations[i] = point[2]; }
             intensities[i] = getIntensity(i);
             returnNumbers[i] = getReturnNumber(i);
             numberOfReturns[i] = getNumberOfReturns(i);
@@ -115,6 +128,7 @@ class LASLoader {
 
         return {
             position: positions,
+            elevation: elevations,
             intensity: intensities,
             returnNumber: returnNumbers,
             numberOfReturns,
@@ -163,7 +177,12 @@ class LASLoader {
         }, this._initDecoder());
 
         const view = Las.View.create(pointData, header, eb);
-        const attributes = this._parseView(view, { colorDepth });
+        const attributes = this._parseView(view, {
+            colorDepth,
+            crsIn: options.crsIn,
+            crsOut: options.crsOut,
+            projDefs: options.projDefs,
+        });
         return { attributes };
     }
 
@@ -190,7 +209,12 @@ class LASLoader {
         const eb = ebVlr && Las.ExtraBytes.parse(await Las.Vlr.fetch(getter, ebVlr));
 
         const view = Las.View.create(pointData, header, eb);
-        const attributes = this._parseView(view, { colorDepth });
+        const attributes = this._parseView(view, {
+            colorDepth,
+            crsIn: options.crsIn,
+            crsOut: options.crsOut,
+            projDefs: options.projDefs,
+        });
         return {
             header,
             attributes,

src/Parser/LASParser.js

@@ -1,5 +1,6 @@
 import * as THREE from 'three';
 import { spawn, Thread, Transfer } from 'threads';
+import proj4 from 'proj4';
 
 let _lazPerf;
 let _thread;
@@ -24,6 +25,8 @@ function buildBufferGeometry(attributes) {
 
     const positionBuffer = new THREE.BufferAttribute(attributes.position, 3);
     geometry.setAttribute('position', positionBuffer);
+    const elevationBuffer = new THREE.BufferAttribute(attributes.elevation, 1);
+    geometry.setAttribute('elevation', elevationBuffer);
 
     const intensityBuffer = new THREE.BufferAttribute(attributes.intensity, 1);
     geometry.setAttribute('intensity', intensityBuffer);
@@ -105,16 +108,24 @@ export default {
      * `THREE.BufferGeometry`.
      */
     async parseChunk(data, options = {}) {
+        const crsIn = options.in?.crs || 'EPSG:3857';
+        const crsOut = options.out?.crs || crsIn;
+
         const lasLoader = await loader();
         const parsedData = await lasLoader.parseChunk(Transfer(data), {
             pointCount: options.in.pointCount,
             header: options.in.header,
             eb: options.eb,
             colorDepth: options.in.colorDepth,
+            crsIn,
+            crsOut,
+            projDefs: {
+                [crsIn]: proj4.defs(crsIn),
+                [crsOut]: proj4.defs(crsOut),
+            },
         });
 
         const geometry = buildBufferGeometry(parsedData.attributes);
-        geometry.computeBoundingBox();
         return geometry;
     },
 
@@ -128,6 +139,8 @@ export default {
      * @param { 8 | 16 } [options.in.colorDepth] - Color depth (in bits).
      * Defaults to 8 bits for LAS 1.2 and 16 bits for later versions
      * (as mandatory by the specification)
+     * @param {String} [options.in.crs = 'EPSG:3857'] - Crs of the source if any.
+     * @param {String} [options.out.crs = options.in.crs] - Crs of the view if any.
      *
      * @return {Promise} A promise resolving with a `THREE.BufferGeometry`. The
      * header of the file is contained in `userData`.
@@ -137,16 +150,22 @@ export default {
             console.warn("Warning: options 'skip' not supported anymore");
         }
 
-        const input = options.in;
+        const crsIn = options.in?.crs || 'EPSG:3857';
+        const crsOut = options.out?.crs || crsIn;
 
         const lasLoader = await loader();
         const parsedData = await lasLoader.parseFile(Transfer(data), {
-            colorDepth: input?.colorDepth,
+            colorDepth: options.in?.colorDepth,
+            crsIn,
+            crsOut,
+            projDefs: {
+                [crsIn]: proj4.defs(crsIn),
+                [crsOut]: proj4.defs(crsOut),
+            },
         });
 
         const geometry = buildBufferGeometry(parsedData.attributes);
         geometry.userData.header = parsedData.header;
-        geometry.computeBoundingBox();
         return geometry;
     },
 };

src/Renderer/Shader/PointsVS.glsl

@@ -33,6 +33,7 @@ attribute vec4 unique_id;
 attribute float intensity;
 attribute float classification;
 attribute float pointSourceID;
+attribute float elevation;
 
 attribute float returnNumber;
 attribute float numberOfReturns;
@@ -95,8 +96,7 @@ void main() {
             vec2 uv = vec2(i, (1. - i));
             vColor = texture2D(gradientTexture, uv);
         } else if (mode == PNTS_MODE_ELEVATION) {
-            float z = (modelMatrix * vec4(position, 1.0)).z;
-            float i = (z - elevationRange.x) / (elevationRange.y - elevationRange.x);
+            float i = (elevation - elevationRange.x) / (elevationRange.y - elevationRange.x);
             vec2 uv = vec2(i, (1. - i));
             vColor = texture2D(gradientTexture, uv);
         }