QGIS API Documentation 3.41.0-Master (fda2aa46e9a)
Loading...
Searching...
No Matches
qgslazdecoder.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgslazdecoder.cpp
3 --------------------
4 begin : March 2022
5 copyright : (C) 2022 by Belgacem Nedjima
6 email : belgacem dot nedjima at gmail dot com
7 ***************************************************************************/
8
9/***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17
18#include "qgslazdecoder.h"
20#include "qgsvector3d.h"
21#include "qgsconfig.h"
22#include "qgslogger.h"
23#include "qgslazinfo.h"
24#include "qgspointcloudexpression.h"
25
26#include <QFile>
27#include <QDir>
28#include <QElapsedTimer>
29#include <QTemporaryFile>
30#include <iostream>
31#include <memory>
32#include <cstring>
33#include <string>
34
35#include <zstd.h>
36
37#include "lazperf/las.hpp"
38#include "lazperf/readers.hpp"
39
40#if defined(_MSC_VER)
41#ifndef UNICODE
42#define UNICODE
43#endif
44#include <locale>
45#include <codecvt>
46#endif
47
49
50template <typename T>
51bool lazStoreToStream_( char *s, size_t position, QgsPointCloudAttribute::DataType type, T value )
52{
53 switch ( type )
54 {
56 {
57 const char val = char( value );
58 s[position] = val;
59 break;
60 }
62 {
63 const unsigned char val = ( unsigned char )( value );
64 s[position] = val;
65 break;
66 }
67
69 {
70 short val = short( value );
71 memcpy( s + position, reinterpret_cast<char * >( &val ), sizeof( short ) );
72 break;
73 }
75 {
76 unsigned short val = static_cast< unsigned short>( value );
77 memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( unsigned short ) );
78 break;
79 }
80
82 {
83 qint32 val = qint32( value );
84 memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( qint32 ) );
85 break;
86 }
88 {
89 quint32 val = quint32( value );
90 memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( quint32 ) );
91 break;
92 }
93
95 {
96 qint64 val = qint64( value );
97 memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( qint64 ) );
98 break;
99 }
101 {
102 quint64 val = quint64( value );
103 memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( quint64 ) );
104 break;
105 }
106
108 {
109 float val = float( value );
110 memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( float ) );
111 break;
112 }
114 {
115 double val = double( value );
116 memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( double ) );
117 break;
118 }
119 }
120
121 return true;
122}
123
124bool lazSerialize_( char *data, size_t outputPosition, QgsPointCloudAttribute::DataType outputType,
125 const char *input, QgsPointCloudAttribute::DataType inputType, int inputSize, size_t inputPosition )
126{
127 if ( outputType == inputType )
128 {
129 memcpy( data + outputPosition, input + inputPosition, inputSize );
130 return true;
131 }
132
133 switch ( inputType )
134 {
136 {
137 const char val = *( input + inputPosition );
138 return lazStoreToStream_<char>( data, outputPosition, outputType, val );
139 }
141 {
142 const unsigned char val = *( input + inputPosition );
143 return lazStoreToStream_<unsigned char>( data, outputPosition, outputType, val );
144 }
146 {
147 const short val = *reinterpret_cast< const short * >( input + inputPosition );
148 return lazStoreToStream_<short>( data, outputPosition, outputType, val );
149 }
151 {
152 const unsigned short val = *reinterpret_cast< const unsigned short * >( input + inputPosition );
153 return lazStoreToStream_<unsigned short>( data, outputPosition, outputType, val );
154 }
156 {
157 const qint32 val = *reinterpret_cast<const qint32 * >( input + inputPosition );
158 return lazStoreToStream_<qint32>( data, outputPosition, outputType, val );
159 }
161 {
162 const quint32 val = *reinterpret_cast<const quint32 * >( input + inputPosition );
163 return lazStoreToStream_<quint32>( data, outputPosition, outputType, val );
164 }
166 {
167 const qint64 val = *reinterpret_cast<const qint64 * >( input + inputPosition );
168 return lazStoreToStream_<qint64>( data, outputPosition, outputType, val );
169 }
171 {
172 const quint64 val = *reinterpret_cast<const quint64 * >( input + inputPosition );
173 return lazStoreToStream_<quint64>( data, outputPosition, outputType, val );
174 }
176 {
177 const float val = *reinterpret_cast< const float * >( input + inputPosition );
178 return lazStoreToStream_<float>( data, outputPosition, outputType, val );
179 }
181 {
182 const double val = *reinterpret_cast< const double * >( input + inputPosition );
183 return lazStoreToStream_<double>( data, outputPosition, outputType, val );
184 }
185 }
186 return true;
187}
188
189// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
190
191std::vector< QgsLazDecoder::RequestedAttributeDetails > prepareRequestedAttributeDetails_( const QgsPointCloudAttributeCollection &requestedAttributes, QVector<QgsLazInfo::ExtraBytesAttributeDetails> &extrabytesAttr )
192{
193 const QVector<QgsPointCloudAttribute> requestedAttributesVector = requestedAttributes.attributes();
194
195 std::vector< QgsLazDecoder::RequestedAttributeDetails > requestedAttributeDetails;
196 requestedAttributeDetails.reserve( requestedAttributesVector.size() );
197
198 for ( const QgsPointCloudAttribute &requestedAttribute : requestedAttributesVector )
199 {
200 if ( requestedAttribute.name().compare( QLatin1String( "X" ), Qt::CaseInsensitive ) == 0 )
201 {
202 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::X, requestedAttribute.type(), requestedAttribute.size() ) );
203 }
204 else if ( requestedAttribute.name().compare( QLatin1String( "Y" ), Qt::CaseInsensitive ) == 0 )
205 {
206 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Y, requestedAttribute.type(), requestedAttribute.size() ) );
207 }
208 else if ( requestedAttribute.name().compare( QLatin1String( "Z" ), Qt::CaseInsensitive ) == 0 )
209 {
210 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Z, requestedAttribute.type(), requestedAttribute.size() ) );
211 }
212 else if ( requestedAttribute.name().compare( QLatin1String( "Classification" ), Qt::CaseInsensitive ) == 0 )
213 {
214 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Classification, requestedAttribute.type(), requestedAttribute.size() ) );
215 }
216 else if ( requestedAttribute.name().compare( QLatin1String( "Intensity" ), Qt::CaseInsensitive ) == 0 )
217 {
218 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Intensity, requestedAttribute.type(), requestedAttribute.size() ) );
219 }
220 else if ( requestedAttribute.name().compare( QLatin1String( "ReturnNumber" ), Qt::CaseInsensitive ) == 0 )
221 {
222 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::ReturnNumber, requestedAttribute.type(), requestedAttribute.size() ) );
223 }
224 else if ( requestedAttribute.name().compare( QLatin1String( "NumberOfReturns" ), Qt::CaseInsensitive ) == 0 )
225 {
226 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::NumberOfReturns, requestedAttribute.type(), requestedAttribute.size() ) );
227 }
228 else if ( requestedAttribute.name().compare( QLatin1String( "ScanDirectionFlag" ), Qt::CaseInsensitive ) == 0 )
229 {
230 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::ScanDirectionFlag, requestedAttribute.type(), requestedAttribute.size() ) );
231 }
232 else if ( requestedAttribute.name().compare( QLatin1String( "EdgeOfFlightLine" ), Qt::CaseInsensitive ) == 0 )
233 {
234 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::EdgeOfFlightLine, requestedAttribute.type(), requestedAttribute.size() ) );
235 }
236 else if ( requestedAttribute.name().compare( QLatin1String( "ScanAngleRank" ), Qt::CaseInsensitive ) == 0 )
237 {
238 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::ScanAngleRank, requestedAttribute.type(), requestedAttribute.size() ) );
239 }
240 else if ( requestedAttribute.name().compare( QLatin1String( "UserData" ), Qt::CaseInsensitive ) == 0 )
241 {
242 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::UserData, requestedAttribute.type(), requestedAttribute.size() ) );
243 }
244 else if ( requestedAttribute.name().compare( QLatin1String( "PointSourceId" ), Qt::CaseInsensitive ) == 0 )
245 {
246 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::PointSourceId, requestedAttribute.type(), requestedAttribute.size() ) );
247 }
248 else if ( requestedAttribute.name().compare( QLatin1String( "GpsTime" ), Qt::CaseInsensitive ) == 0 )
249 {
250 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::GpsTime, requestedAttribute.type(), requestedAttribute.size() ) );
251 }
252 else if ( requestedAttribute.name().compare( QLatin1String( "Red" ), Qt::CaseInsensitive ) == 0 )
253 {
254 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Red, requestedAttribute.type(), requestedAttribute.size() ) );
255 }
256 else if ( requestedAttribute.name().compare( QLatin1String( "Green" ), Qt::CaseInsensitive ) == 0 )
257 {
258 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Green, requestedAttribute.type(), requestedAttribute.size() ) );
259 }
260 else if ( requestedAttribute.name().compare( QLatin1String( "Blue" ), Qt::CaseInsensitive ) == 0 )
261 {
262 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Blue, requestedAttribute.type(), requestedAttribute.size() ) );
263 }
264 else if ( requestedAttribute.name().compare( QLatin1String( "ScannerChannel" ), Qt::CaseInsensitive ) == 0 )
265 {
266 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::ScannerChannel, requestedAttribute.type(), requestedAttribute.size() ) );
267 }
268 else if ( requestedAttribute.name().compare( QLatin1String( "Synthetic" ), Qt::CaseInsensitive ) == 0 )
269 {
270 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Synthetic, requestedAttribute.type(), requestedAttribute.size() ) );
271 }
272 else if ( requestedAttribute.name().compare( QLatin1String( "KeyPoint" ), Qt::CaseInsensitive ) == 0 )
273 {
274 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::KeyPoint, requestedAttribute.type(), requestedAttribute.size() ) );
275 }
276 else if ( requestedAttribute.name().compare( QLatin1String( "Withheld" ), Qt::CaseInsensitive ) == 0 )
277 {
278 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Withheld, requestedAttribute.type(), requestedAttribute.size() ) );
279 }
280 else if ( requestedAttribute.name().compare( QLatin1String( "Overlap" ), Qt::CaseInsensitive ) == 0 )
281 {
282 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Overlap, requestedAttribute.type(), requestedAttribute.size() ) );
283 }
284 else if ( requestedAttribute.name().compare( QLatin1String( "Infrared" ), Qt::CaseInsensitive ) == 0 )
285 {
286 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::NIR, requestedAttribute.type(), requestedAttribute.size() ) );
287 }
288 else
289 {
290 bool foundAttr = false;
291 for ( QgsLazInfo::ExtraBytesAttributeDetails &eba : extrabytesAttr )
292 {
293 if ( requestedAttribute.name().compare( eba.attribute.trimmed() ) == 0 )
294 {
295 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::ExtraBytes, eba.type, eba.size, eba.offset ) );
296 foundAttr = true;
297 break;
298 }
299 }
300 if ( !foundAttr )
301 {
302 // this can possibly happen -- e.g. if a style built using a different point cloud format references an attribute which isn't available from the laz file
303 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::MissingOrUnknown, requestedAttribute.type(), requestedAttribute.size() ) );
304 }
305 }
306 }
307 return requestedAttributeDetails;
308}
309
310bool decodePoint( char *buf, int lasPointFormat, char *dataBuffer, std::size_t &outputOffset, std::vector< QgsLazDecoder::RequestedAttributeDetails > &requestedAttributeDetails )
311{
312 lazperf::las::point10 p10;
313 lazperf::las::gpstime gps;
314 lazperf::las::rgb rgb;
315 lazperf::las::nir14 nir;
316 lazperf::las::point14 p14;
317
318 // Does the point record start with the common fields for formats introduced
319 // in the LAS 1.4 spec?
320 const bool isLas14 = ( lasPointFormat == 6 || lasPointFormat == 7 || lasPointFormat == 8 || lasPointFormat == 9 || lasPointFormat == 10 );
321
322 switch ( lasPointFormat )
323 {
324 // LAS 1.2 file support
325 case 0: // base
326 p10.unpack( buf );
327 break;
328 case 1: // base + gps time
329 p10.unpack( buf );
330 gps.unpack( buf + sizeof( lazperf::las::point10 ) );
331 break;
332 case 2: // base + rgb
333 p10.unpack( buf );
334 rgb.unpack( buf + sizeof( lazperf::las::point10 ) );
335 break;
336 case 3: // base + gps time + rgb
337 p10.unpack( buf );
338 gps.unpack( buf + sizeof( lazperf::las::point10 ) );
339 rgb.unpack( buf + sizeof( lazperf::las::point10 ) + sizeof( lazperf::las::gpstime ) );
340 break;
341
342 // LAS 1.4 file support
343 case 6: // base (includes gps time)
344 p14.unpack( buf );
345 break;
346 case 7: // base + rgb
347 p14.unpack( buf );
348 rgb.unpack( buf + sizeof( lazperf::las::point14 ) );
349 break;
350 case 8: // base + rgb + nir
351 p14.unpack( buf );
352 rgb.unpack( buf + sizeof( lazperf::las::point14 ) );
353 nir.unpack( buf + sizeof( lazperf::las::point14 ) + sizeof( lazperf::las::rgb ) );
354 break;
355
356 default:
357 Q_ASSERT( false ); // must not happen - we checked earlier that the format is supported
358 return false;
359 }
360
361 for ( const QgsLazDecoder::RequestedAttributeDetails &requestedAttribute : requestedAttributeDetails )
362 {
363 switch ( requestedAttribute.attribute )
364 {
365 case QgsLazDecoder::LazAttribute::X:
366 lazStoreToStream_<qint32>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.x() : p10.x );
367 break;
368 case QgsLazDecoder::LazAttribute::Y:
369 lazStoreToStream_<qint32>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.y() : p10.y );
370 break;
371 case QgsLazDecoder::LazAttribute::Z:
372 lazStoreToStream_<qint32>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.z() : p10.z );
373 break;
374 case QgsLazDecoder::LazAttribute::Classification:
375 {
376 if ( isLas14 )
377 {
378 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, p14.classification() );
379 }
380 else
381 {
382 // p10 format encoded "Overlap" as Classification=12, so in that case we set Classification=0 (Never classified) and will set Overlap=1 a few lines below
383 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, ( p10.classification & 0x1F ) == 12 ? 0 : p10.classification & 0x1F );
384 }
385 break;
386 }
387 case QgsLazDecoder::LazAttribute::Intensity:
388 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.intensity() : p10.intensity );
389 break;
390 case QgsLazDecoder::LazAttribute::ReturnNumber:
391 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.returnNum() : p10.return_number );
392 break;
393 case QgsLazDecoder::LazAttribute::NumberOfReturns:
394 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.numReturns() : p10.number_of_returns_of_given_pulse );
395 break;
396 case QgsLazDecoder::LazAttribute::ScanDirectionFlag:
397 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.scanDirFlag() : p10.scan_direction_flag );
398 break;
399 case QgsLazDecoder::LazAttribute::EdgeOfFlightLine:
400 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.eofFlag() : p10.edge_of_flight_line );
401 break;
402 case QgsLazDecoder::LazAttribute::ScanAngleRank:
403 lazStoreToStream_<float>( dataBuffer, outputOffset, requestedAttribute.type,
404 isLas14
405 // Formats from LAS 1.4 spec store the angle in 0.006 degree increments
406 ? p14.scanAngle() * 0.006f
407 // Older formats store integer values
408 : p10.scan_angle_rank );
409 break;
410 case QgsLazDecoder::LazAttribute::UserData:
411 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.userData() : p10.user_data );
412 break;
413 case QgsLazDecoder::LazAttribute::PointSourceId:
414 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.pointSourceID() : p10.point_source_ID );
415 break;
416 case QgsLazDecoder::LazAttribute::GpsTime:
417 // lazperf internally stores gps value as int64 field, but in fact it is a double value
418 lazStoreToStream_<double>( dataBuffer, outputOffset, requestedAttribute.type,
419 isLas14 ? p14.gpsTime() : *reinterpret_cast<const double *>( reinterpret_cast<const void *>( &gps.value ) ) );
420 break;
421 case QgsLazDecoder::LazAttribute::Red:
422 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, rgb.r );
423 break;
424 case QgsLazDecoder::LazAttribute::Green:
425 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, rgb.g );
426 break;
427 case QgsLazDecoder::LazAttribute::Blue:
428 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, rgb.b );
429 break;
430 case QgsLazDecoder::LazAttribute::ScannerChannel:
431 lazStoreToStream_<char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? char( p14.scannerChannel() ) : 0 );
432 break;
433 case QgsLazDecoder::LazAttribute::Synthetic:
434 lazStoreToStream_<char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? char( ( p14.classFlags() >> 0 ) & 0x01 ) : char( ( p10.classification >> 5 ) & 0x01 ) );
435 break;
436 case QgsLazDecoder::LazAttribute::KeyPoint:
437 lazStoreToStream_<char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? char( ( p14.classFlags() >> 1 ) & 0x01 ) : char( ( p10.classification >> 6 ) & 0x01 ) );
438 break;
439 case QgsLazDecoder::LazAttribute::Withheld:
440 lazStoreToStream_<char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? char( ( p14.classFlags() >> 2 ) & 0x01 ) : char( ( p10.classification >> 7 ) & 0x01 ) );
441 break;
442 case QgsLazDecoder::LazAttribute::Overlap:
443 {
444 if ( isLas14 )
445 {
446 lazStoreToStream_<char>( dataBuffer, outputOffset, requestedAttribute.type, char( ( p14.classFlags() >> 3 ) & 0x01 ) );
447 }
448 else
449 {
450 // p10 format encoded "Overlap" as Classification=12, so in that case we set Overlap=1 (we have already set Classification=0)
451 lazStoreToStream_<char>( dataBuffer, outputOffset, requestedAttribute.type, ( p10.classification & 0x1F ) == 12 ? 1 : 0 );
452 }
453 break;
454 }
455 case QgsLazDecoder::LazAttribute::NIR:
456 {
457 if ( lasPointFormat == 8 || lasPointFormat == 10 )
458 {
459 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, nir.val );
460 }
461 else
462 {
463 // just store 0 for missing attributes
464 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, 0 );
465 }
466 break;
467
468 }
469 case QgsLazDecoder::LazAttribute::ExtraBytes:
470 {
471 switch ( requestedAttribute.type )
472 {
474 lazStoreToStream_<char>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<char * >( &buf[requestedAttribute.offset] ) );
475 break;
477 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<unsigned char * >( &buf[requestedAttribute.offset] ) );
478 break;
480 lazStoreToStream_<qint16>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<qint16 * >( &buf[requestedAttribute.offset] ) );
481 break;
483 lazStoreToStream_<quint16>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<quint16 * >( &buf[requestedAttribute.offset] ) );
484 break;
486 lazStoreToStream_<qint32>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<qint32 * >( &buf[requestedAttribute.offset] ) );
487 break;
489 lazStoreToStream_<quint32>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<quint32 * >( &buf[requestedAttribute.offset] ) );
490 break;
492 lazStoreToStream_<qint64>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<qint64 * >( &buf[requestedAttribute.offset] ) );
493 break;
495 lazStoreToStream_<quint64>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<quint64 * >( &buf[requestedAttribute.offset] ) );
496 break;
498 lazStoreToStream_<float>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<float * >( &buf[requestedAttribute.offset] ) );
499 break;
501 lazStoreToStream_<double>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<double * >( &buf[requestedAttribute.offset] ) );
502 break;
503 }
504 }
505 break;
506 case QgsLazDecoder::LazAttribute::MissingOrUnknown:
507 // just store 0 for unknown/missing attributes
508 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, 0 );
509 break;
510 }
511
512 outputOffset += requestedAttribute.size;
513 }
514 return true;
515}
516
517template<typename FileType>
518std::unique_ptr<QgsPointCloudBlock> decompressLaz_( FileType &file, const QgsPointCloudAttributeCollection &requestedAttributes, QgsPointCloudExpression &filterExpression, QgsRectangle &filterRect )
519{
520 if ( ! file.good() )
521 return nullptr;
522
523#ifdef QGISDEBUG
524 QElapsedTimer t;
525 t.start();
526#endif
527
528 // lazperf may throw exceptions
529 try
530 {
531 lazperf::reader::generic_file f( file );
532
533
534 // output file formats from entwine/untwine:
535 // - older versions write LAZ 1.2 files with point formats 0, 1, 2 or 3
536 // - newer versions write LAZ 1.4 files with point formats 6, 7 or 8
537
538 int lasPointFormat = f.header().pointFormat();
539 if ( lasPointFormat != 0 && lasPointFormat != 1 && lasPointFormat != 2 && lasPointFormat != 3 &&
540 lasPointFormat != 6 && lasPointFormat != 7 && lasPointFormat != 8 )
541 {
542 QgsDebugError( QStringLiteral( "Unexpected point format record (%1) - only 0, 1, 2, 3, 6, 7, 8 are supported" ).arg( lasPointFormat ) );
543 return nullptr;
544 }
545
546 const size_t count = f.header().point_count;
547 const QgsVector3D scale( f.header().scale.x, f.header().scale.y, f.header().scale.z );
548 const QgsVector3D offset( f.header().offset.x, f.header().offset.y, f.header().offset.z );
549
550 QByteArray bufArray( f.header().point_record_length, 0 );
551 char *buf = bufArray.data();
552
553 const size_t requestedPointRecordSize = requestedAttributes.pointRecordSize();
554 QByteArray data;
555 data.resize( requestedPointRecordSize * count );
556 char *dataBuffer = data.data();
557
558 std::size_t outputOffset = 0;
559
560 std::unique_ptr< QgsPointCloudBlock > block = std::make_unique< QgsPointCloudBlock >(
561 count,
562 requestedAttributes,
563 data, scale, offset
564 );
565
566 int skippedPoints = 0;
567 const bool filterIsValid = filterExpression.isValid();
568 if ( !filterExpression.prepare( block.get() ) && filterIsValid )
569 {
570 // skip processing if the expression cannot be prepared
571 block->setPointCount( 0 );
572 return block;
573 }
574
575 int xAttributeOffset, yAttributeOffset;
576 const QgsPointCloudAttribute *attributeX = nullptr;
577 const QgsPointCloudAttribute *attributeY = nullptr;
578 const bool hasFilterRect = !filterRect.isEmpty();
579 if ( hasFilterRect )
580 {
581 attributeX = requestedAttributes.find( QLatin1String( "X" ), xAttributeOffset );
582 attributeY = requestedAttributes.find( QLatin1String( "Y" ), yAttributeOffset );
583 filterRect.setXMinimum( ( filterRect.xMinimum() - offset.x() ) / scale.x() );
584 filterRect.setXMaximum( ( filterRect.xMaximum() - offset.x() ) / scale.x() );
585 filterRect.setYMinimum( ( filterRect.yMinimum() - offset.y() ) / scale.y() );
586 filterRect.setYMaximum( ( filterRect.yMaximum() - offset.y() ) / scale.y() );
587 }
588
589 std::vector<char> rawExtrabytes = f.vlrData( "LASF_Spec", 4 );
590 QVector<QgsLazInfo::ExtraBytesAttributeDetails> extrabyteAttributesDetails = QgsLazInfo::parseExtrabytes( rawExtrabytes.data(), rawExtrabytes.size(), f.header().point_record_length );
591 std::vector< QgsLazDecoder::RequestedAttributeDetails > requestedAttributeDetails = prepareRequestedAttributeDetails_( requestedAttributes, extrabyteAttributesDetails );
592
593 for ( size_t i = 0 ; i < count ; i ++ )
594 {
595 f.readPoint( buf ); // read the point out
596
597 bool skipThisPoint = !decodePoint( buf, lasPointFormat, dataBuffer, outputOffset, requestedAttributeDetails );
598
599 // check if point needs to be filtered out
600 if ( !skipThisPoint && hasFilterRect && attributeX && attributeY )
601 {
602 const double x = attributeX->convertValueToDouble( dataBuffer + outputOffset - requestedPointRecordSize + xAttributeOffset );
603 const double y = attributeY->convertValueToDouble( dataBuffer + outputOffset - requestedPointRecordSize + yAttributeOffset );
604 if ( !filterRect.contains( x, y ) )
605 skipThisPoint = true;
606 }
607 if ( !skipThisPoint && filterIsValid )
608 {
609 // we're always evaluating the last written point in the buffer
610 double eval = filterExpression.evaluate( i - skippedPoints );
611 if ( !eval || std::isnan( eval ) )
612 skipThisPoint = true;
613 }
614 if ( skipThisPoint )
615 {
616 // if the point is filtered out, rewind the offset so the next point is written over it
617 outputOffset -= requestedPointRecordSize;
618 ++skippedPoints;
619 }
620 }
621
622#ifdef QGISDEBUG
623 QgsDebugMsgLevel( QStringLiteral( "LAZ-PERF Read through the points in %1 seconds." ).arg( t.elapsed() / 1000. ), 2 );
624#endif
625 block->setPointCount( count - skippedPoints );
626 return block;
627 }
628 catch ( std::exception &e )
629 {
630 QgsDebugError( "Error decompressing laz file: " + QString::fromLatin1( e.what() ) );
631 return nullptr;
632 }
633}
634
635std::unique_ptr<QgsPointCloudBlock> QgsLazDecoder::decompressLaz( const QString &filename,
636 const QgsPointCloudAttributeCollection &requestedAttributes,
637 QgsPointCloudExpression &filterExpression, QgsRectangle &filterRect )
638{
639 std::ifstream file( toNativePath( filename ), std::ios::binary );
640
641 return decompressLaz_<std::ifstream>( file, requestedAttributes, filterExpression, filterRect );
642}
643
644std::unique_ptr<QgsPointCloudBlock> QgsLazDecoder::decompressLaz( const QByteArray &byteArrayData,
645 const QgsPointCloudAttributeCollection &requestedAttributes,
646 QgsPointCloudExpression &filterExpression, QgsRectangle &filterRect )
647{
648 std::istringstream file( byteArrayData.toStdString() );
649 return decompressLaz_<std::istringstream>( file, requestedAttributes, filterExpression, filterRect );
650}
651
652std::unique_ptr<QgsPointCloudBlock> QgsLazDecoder::decompressCopc( const QByteArray &data, QgsLazInfo &lazInfo, int32_t pointCount, const QgsPointCloudAttributeCollection &requestedAttributes, QgsPointCloudExpression &filterExpression, QgsRectangle &filterRect )
653{
654 // COPC only supports point formats 6, 7 and 8
655 int lasPointFormat = lazInfo.pointFormat();
656 if ( lasPointFormat != 6 && lasPointFormat != 7 && lasPointFormat != 8 )
657 {
658 QgsDebugError( QStringLiteral( "Unexpected point format record (%1) - only 6, 7, 8 are supported for COPC format" ).arg( lasPointFormat ) );
659 return nullptr;
660 }
661
662 std::unique_ptr<char []> decodedData( new char[ lazInfo.pointRecordLength() ] );
663
664 lazperf::reader::chunk_decompressor decompressor( lasPointFormat, lazInfo.extrabytesCount(), data.data() );
665
666 const size_t requestedPointRecordSize = requestedAttributes.pointRecordSize();
667 QByteArray blockData;
668 blockData.resize( requestedPointRecordSize * pointCount );
669 char *dataBuffer = blockData.data();
670
671 std::size_t outputOffset = 0;
672
673 QVector<QgsLazInfo::ExtraBytesAttributeDetails> extrabyteAttributesDetails = lazInfo.extrabytes();
674 std::vector< RequestedAttributeDetails > requestedAttributeDetails = prepareRequestedAttributeDetails_( requestedAttributes, extrabyteAttributesDetails );
675 std::unique_ptr< QgsPointCloudBlock > block = std::make_unique< QgsPointCloudBlock >(
676 pointCount, requestedAttributes,
677 blockData, lazInfo.scale(), lazInfo.offset()
678 );
679
680 int skippedPoints = 0;
681 const bool filterIsValid = filterExpression.isValid();
682 if ( !filterExpression.prepare( block.get() ) && filterIsValid )
683 {
684 // skip processing if the expression cannot be prepared
685 block->setPointCount( 0 );
686 return block;
687 }
688
689 int xAttributeOffset, yAttributeOffset;
690 const QgsPointCloudAttribute *attributeX = nullptr;
691 const QgsPointCloudAttribute *attributeY = nullptr;
692 const bool hasFilterRect = !filterRect.isEmpty();
693 if ( hasFilterRect )
694 {
695 attributeX = requestedAttributes.find( QLatin1String( "X" ), xAttributeOffset );
696 attributeY = requestedAttributes.find( QLatin1String( "Y" ), yAttributeOffset );
697 filterRect.setXMinimum( ( filterRect.xMinimum() - lazInfo.offset().x() ) / lazInfo.scale().x() );
698 filterRect.setXMaximum( ( filterRect.xMaximum() - lazInfo.offset().x() ) / lazInfo.scale().x() );
699 filterRect.setYMinimum( ( filterRect.yMinimum() - lazInfo.offset().y() ) / lazInfo.scale().y() );
700 filterRect.setYMaximum( ( filterRect.yMaximum() - lazInfo.offset().y() ) / lazInfo.scale().y() );
701 }
702 for ( int i = 0 ; i < pointCount; ++i )
703 {
704 decompressor.decompress( decodedData.get() );
705 char *buf = decodedData.get();
706
707 bool skipThisPoint = !decodePoint( buf, lasPointFormat, dataBuffer, outputOffset, requestedAttributeDetails );
708
709 // check if point needs to be filtered out
710 if ( !skipThisPoint && hasFilterRect && attributeX && attributeY )
711 {
712 const double x = attributeX->convertValueToDouble( dataBuffer + outputOffset - requestedPointRecordSize + xAttributeOffset );
713 const double y = attributeY->convertValueToDouble( dataBuffer + outputOffset - requestedPointRecordSize + yAttributeOffset );
714 if ( !filterRect.contains( x, y ) )
715 skipThisPoint = true;
716 }
717 if ( !skipThisPoint && filterIsValid )
718 {
719 // we're always evaluating the last written point in the buffer
720 double eval = filterExpression.evaluate( i - skippedPoints );
721 if ( !eval || std::isnan( eval ) )
722 skipThisPoint = true;
723 }
724 if ( skipThisPoint )
725 {
726 // if the point is filtered out, rewind the offset so the next point is written over it
727 outputOffset -= requestedPointRecordSize;
728 ++skippedPoints;
729 }
730 }
731
732 block->setPointCount( pointCount - skippedPoints );
733 return block;
734}
735
736#if defined(_MSC_VER)
737std::wstring QgsLazDecoder::toNativePath( const QString &filename )
738{
739 std::wstring_convert< std::codecvt_utf8_utf16< wchar_t > > converter;
740 return converter.from_bytes( filename.toStdString() );
741}
742#else
743std::string QgsLazDecoder::toNativePath( const QString &filename )
744{
745 return filename.toStdString();
746}
747#endif
748
Class for extracting information contained in LAZ file such as the public header block and variable l...
Definition qgslazinfo.h:39
int extrabytesCount() const
Returns the number of extrabytes contained in the LAZ dataset.
Definition qgslazinfo.h:103
QgsVector3D scale() const
Returns the scale of the points coordinates.
Definition qgslazinfo.h:77
int pointFormat() const
Returns the point format of the point records contained in the LAZ file.
Definition qgslazinfo.h:85
QVector< ExtraBytesAttributeDetails > extrabytes() const
Returns the list of extrabytes contained in the LAZ file.
Definition qgslazinfo.h:123
QgsVector3D offset() const
Returns the offset of the points coordinates.
Definition qgslazinfo.h:79
int pointRecordLength() const
Returns the length of each point record in bytes.
Definition qgslazinfo.h:101
static QVector< ExtraBytesAttributeDetails > parseExtrabytes(char *rawData, int length, int pointRecordLength)
Static function to parse the raw extrabytes VLR into a list of recognizable extrabyte attributes.
Collection of point cloud attributes.
int pointRecordSize() const
Returns total size of record.
const QgsPointCloudAttribute * find(const QString &attributeName, int &offset) const
Finds the attribute with the name.
QVector< QgsPointCloudAttribute > attributes() const
Returns all attributes.
Attribute for point cloud data pair of name and size in bytes.
DataType
Systems of unit measurement.
@ UShort
Unsigned short int 2 bytes.
@ UChar
Unsigned char 1 byte.
@ UInt32
Unsigned int32 4 bytes.
@ UInt64
Unsigned int64 8 bytes.
double convertValueToDouble(const char *ptr) const
Returns the attribute's value as a double for data pointed to by ptr.
A rectangle specified with double values.
bool contains(const QgsRectangle &rect) const
Returns true when rectangle contains other rectangle.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
void setYMinimum(double y)
Set the minimum y value.
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
void setXMinimum(double x)
Set the minimum x value.
double xMaximum() const
Returns the x maximum value (right side of rectangle).
double yMaximum() const
Returns the y maximum value (top side of rectangle).
void setYMaximum(double y)
Set the maximum y value.
void setXMaximum(double x)
Set the maximum x value.
bool isEmpty() const
Returns true if the rectangle has no area.
Class for storage of 3D vectors similar to QVector3D, with the difference that it uses double precisi...
Definition qgsvector3d.h:31
double y() const
Returns Y coordinate.
Definition qgsvector3d.h:50
double x() const
Returns X coordinate.
Definition qgsvector3d.h:48
#define QgsDebugMsgLevel(str, level)
Definition qgslogger.h:39
#define QgsDebugError(str)
Definition qgslogger.h:38