QGIS API Documentation 3.30.0-'s-Hertogenbosch (f186b8efe0)
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#define UNICODE
42#include <locale>
43#include <codecvt>
44#endif
45
47
48template <typename T>
49bool lazStoreToStream_( char *s, size_t position, QgsPointCloudAttribute::DataType type, T value )
50{
51 switch ( type )
52 {
54 {
55 const char val = char( value );
56 s[position] = val;
57 break;
58 }
60 {
61 const unsigned char val = ( unsigned char )( value );
62 s[position] = val;
63 break;
64 }
65
67 {
68 short val = short( value );
69 memcpy( s + position, reinterpret_cast<char * >( &val ), sizeof( short ) );
70 break;
71 }
73 {
74 unsigned short val = static_cast< unsigned short>( value );
75 memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( unsigned short ) );
76 break;
77 }
78
80 {
81 qint32 val = qint32( value );
82 memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( qint32 ) );
83 break;
84 }
86 {
87 quint32 val = quint32( value );
88 memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( quint32 ) );
89 break;
90 }
91
93 {
94 qint64 val = qint64( value );
95 memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( qint64 ) );
96 break;
97 }
99 {
100 quint64 val = quint64( value );
101 memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( quint64 ) );
102 break;
103 }
104
106 {
107 float val = float( value );
108 memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( float ) );
109 break;
110 }
112 {
113 double val = double( value );
114 memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( double ) );
115 break;
116 }
117 }
118
119 return true;
120}
121
122bool lazSerialize_( char *data, size_t outputPosition, QgsPointCloudAttribute::DataType outputType,
123 const char *input, QgsPointCloudAttribute::DataType inputType, int inputSize, size_t inputPosition )
124{
125 if ( outputType == inputType )
126 {
127 memcpy( data + outputPosition, input + inputPosition, inputSize );
128 return true;
129 }
130
131 switch ( inputType )
132 {
134 {
135 const char val = *( input + inputPosition );
136 return lazStoreToStream_<char>( data, outputPosition, outputType, val );
137 }
139 {
140 const unsigned char val = *( input + inputPosition );
141 return lazStoreToStream_<unsigned char>( data, outputPosition, outputType, val );
142 }
144 {
145 const short val = *reinterpret_cast< const short * >( input + inputPosition );
146 return lazStoreToStream_<short>( data, outputPosition, outputType, val );
147 }
149 {
150 const unsigned short val = *reinterpret_cast< const unsigned short * >( input + inputPosition );
151 return lazStoreToStream_<unsigned short>( data, outputPosition, outputType, val );
152 }
154 {
155 const qint32 val = *reinterpret_cast<const qint32 * >( input + inputPosition );
156 return lazStoreToStream_<qint32>( data, outputPosition, outputType, val );
157 }
159 {
160 const quint32 val = *reinterpret_cast<const quint32 * >( input + inputPosition );
161 return lazStoreToStream_<quint32>( data, outputPosition, outputType, val );
162 }
164 {
165 const qint64 val = *reinterpret_cast<const qint64 * >( input + inputPosition );
166 return lazStoreToStream_<qint64>( data, outputPosition, outputType, val );
167 }
169 {
170 const quint64 val = *reinterpret_cast<const quint64 * >( input + inputPosition );
171 return lazStoreToStream_<quint64>( data, outputPosition, outputType, val );
172 }
174 {
175 const float val = *reinterpret_cast< const float * >( input + inputPosition );
176 return lazStoreToStream_<float>( data, outputPosition, outputType, val );
177 }
179 {
180 const double val = *reinterpret_cast< const double * >( input + inputPosition );
181 return lazStoreToStream_<double>( data, outputPosition, outputType, val );
182 }
183 }
184 return true;
185}
186
187// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
188
189std::vector< QgsLazDecoder::RequestedAttributeDetails > prepareRequestedAttributeDetails_( const QgsPointCloudAttributeCollection &requestedAttributes, QVector<QgsLazInfo::ExtraBytesAttributeDetails> &extrabytesAttr )
190{
191 const QVector<QgsPointCloudAttribute> requestedAttributesVector = requestedAttributes.attributes();
192
193 std::vector< QgsLazDecoder::RequestedAttributeDetails > requestedAttributeDetails;
194 requestedAttributeDetails.reserve( requestedAttributesVector.size() );
195 for ( const QgsPointCloudAttribute &requestedAttribute : requestedAttributesVector )
196 {
197 if ( requestedAttribute.name().compare( QLatin1String( "X" ), Qt::CaseInsensitive ) == 0 )
198 {
199 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::X, requestedAttribute.type(), requestedAttribute.size() ) );
200 }
201 else if ( requestedAttribute.name().compare( QLatin1String( "Y" ), Qt::CaseInsensitive ) == 0 )
202 {
203 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Y, requestedAttribute.type(), requestedAttribute.size() ) );
204 }
205 else if ( requestedAttribute.name().compare( QLatin1String( "Z" ), Qt::CaseInsensitive ) == 0 )
206 {
207 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Z, requestedAttribute.type(), requestedAttribute.size() ) );
208 }
209 else if ( requestedAttribute.name().compare( QLatin1String( "Classification" ), Qt::CaseInsensitive ) == 0 )
210 {
211 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Classification, requestedAttribute.type(), requestedAttribute.size() ) );
212 }
213 else if ( requestedAttribute.name().compare( QLatin1String( "Intensity" ), Qt::CaseInsensitive ) == 0 )
214 {
215 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Intensity, requestedAttribute.type(), requestedAttribute.size() ) );
216 }
217 else if ( requestedAttribute.name().compare( QLatin1String( "ReturnNumber" ), Qt::CaseInsensitive ) == 0 )
218 {
219 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::ReturnNumber, requestedAttribute.type(), requestedAttribute.size() ) );
220 }
221 else if ( requestedAttribute.name().compare( QLatin1String( "NumberOfReturns" ), Qt::CaseInsensitive ) == 0 )
222 {
223 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::NumberOfReturns, requestedAttribute.type(), requestedAttribute.size() ) );
224 }
225 else if ( requestedAttribute.name().compare( QLatin1String( "ScanDirectionFlag" ), Qt::CaseInsensitive ) == 0 )
226 {
227 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::ScanDirectionFlag, requestedAttribute.type(), requestedAttribute.size() ) );
228 }
229 else if ( requestedAttribute.name().compare( QLatin1String( "EdgeOfFlightLine" ), Qt::CaseInsensitive ) == 0 )
230 {
231 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::EdgeOfFlightLine, requestedAttribute.type(), requestedAttribute.size() ) );
232 }
233 else if ( requestedAttribute.name().compare( QLatin1String( "ScanAngleRank" ), Qt::CaseInsensitive ) == 0 )
234 {
235 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::ScanAngleRank, requestedAttribute.type(), requestedAttribute.size() ) );
236 }
237 else if ( requestedAttribute.name().compare( QLatin1String( "UserData" ), Qt::CaseInsensitive ) == 0 )
238 {
239 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::UserData, requestedAttribute.type(), requestedAttribute.size() ) );
240 }
241 else if ( requestedAttribute.name().compare( QLatin1String( "PointSourceId" ), Qt::CaseInsensitive ) == 0 )
242 {
243 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::PointSourceId, requestedAttribute.type(), requestedAttribute.size() ) );
244 }
245 else if ( requestedAttribute.name().compare( QLatin1String( "GpsTime" ), Qt::CaseInsensitive ) == 0 )
246 {
247 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::GpsTime, requestedAttribute.type(), requestedAttribute.size() ) );
248 }
249 else if ( requestedAttribute.name().compare( QLatin1String( "Red" ), Qt::CaseInsensitive ) == 0 )
250 {
251 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Red, requestedAttribute.type(), requestedAttribute.size() ) );
252 }
253 else if ( requestedAttribute.name().compare( QLatin1String( "Green" ), Qt::CaseInsensitive ) == 0 )
254 {
255 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Green, requestedAttribute.type(), requestedAttribute.size() ) );
256 }
257 else if ( requestedAttribute.name().compare( QLatin1String( "Blue" ), Qt::CaseInsensitive ) == 0 )
258 {
259 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::Blue, requestedAttribute.type(), requestedAttribute.size() ) );
260 }
261 else if ( requestedAttribute.name().compare( QLatin1String( "ScannerChannel" ), Qt::CaseInsensitive ) == 0 )
262 {
263 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::ScannerChannel, requestedAttribute.type(), requestedAttribute.size() ) );
264 }
265 else if ( requestedAttribute.name().compare( QLatin1String( "ClassificationFlags" ), Qt::CaseInsensitive ) == 0 )
266 {
267 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::ClassificationFlags, requestedAttribute.type(), requestedAttribute.size() ) );
268 }
269 else if ( requestedAttribute.name().compare( QLatin1String( "Infrared" ), Qt::CaseInsensitive ) == 0 )
270 {
271 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::NIR, requestedAttribute.type(), requestedAttribute.size() ) );
272 }
273 else
274 {
275 bool foundAttr = false;
276 for ( QgsLazInfo::ExtraBytesAttributeDetails &eba : extrabytesAttr )
277 {
278 if ( requestedAttribute.name().compare( eba.attribute.trimmed() ) == 0 )
279 {
280 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::ExtraBytes, eba.type, eba.size, eba.offset ) );
281 foundAttr = true;
282 break;
283 }
284 }
285 if ( !foundAttr )
286 {
287 // 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
288 requestedAttributeDetails.emplace_back( QgsLazDecoder::RequestedAttributeDetails( QgsLazDecoder::LazAttribute::MissingOrUnknown, requestedAttribute.type(), requestedAttribute.size() ) );
289 }
290 }
291 }
292 return requestedAttributeDetails;
293}
294
295void decodePoint( char *buf, int lasPointFormat, char *dataBuffer, std::size_t &outputOffset, std::vector< QgsLazDecoder::RequestedAttributeDetails > &requestedAttributeDetails )
296{
297 lazperf::las::point10 p10;
298 lazperf::las::gpstime gps;
299 lazperf::las::rgb rgb;
300 lazperf::las::nir14 nir;
301 lazperf::las::point14 p14;
302
303 bool isLas14 = ( lasPointFormat == 6 || lasPointFormat == 7 || lasPointFormat == 8 );
304
305 switch ( lasPointFormat )
306 {
307 // LAS 1.2 file support
308 case 0: // base
309 p10.unpack( buf );
310 break;
311 case 1: // base + gps time
312 p10.unpack( buf );
313 gps.unpack( buf + sizeof( lazperf::las::point10 ) );
314 break;
315 case 2: // base + rgb
316 p10.unpack( buf );
317 rgb.unpack( buf + sizeof( lazperf::las::point10 ) );
318 break;
319 case 3: // base + gps time + rgb
320 p10.unpack( buf );
321 gps.unpack( buf + sizeof( lazperf::las::point10 ) );
322 rgb.unpack( buf + sizeof( lazperf::las::point10 ) + sizeof( lazperf::las::gpstime ) );
323 break;
324
325 // LAS 1.4 file support
326 case 6: // base (includes gps time)
327 p14.unpack( buf );
328 break;
329 case 7: // base + rgb
330 p14.unpack( buf );
331 rgb.unpack( buf + sizeof( lazperf::las::point14 ) );
332 break;
333 case 8: // base + rgb + nir
334 p14.unpack( buf );
335 rgb.unpack( buf + sizeof( lazperf::las::point14 ) );
336 nir.unpack( buf + sizeof( lazperf::las::point14 ) + sizeof( lazperf::las::rgb ) );
337 break;
338
339 default:
340 Q_ASSERT( false ); // must not happen - we checked earlier that the format is supported
341 }
342
343 for ( const QgsLazDecoder::RequestedAttributeDetails &requestedAttribute : requestedAttributeDetails )
344 {
345 switch ( requestedAttribute.attribute )
346 {
347 case QgsLazDecoder::LazAttribute::X:
348 lazStoreToStream_<qint32>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.x() : p10.x );
349 break;
350 case QgsLazDecoder::LazAttribute::Y:
351 lazStoreToStream_<qint32>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.y() : p10.y );
352 break;
353 case QgsLazDecoder::LazAttribute::Z:
354 lazStoreToStream_<qint32>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.z() : p10.z );
355 break;
356 case QgsLazDecoder::LazAttribute::Classification:
357 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.classification() : p10.classification );
358 break;
359 case QgsLazDecoder::LazAttribute::Intensity:
360 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.intensity() : p10.intensity );
361 break;
362 case QgsLazDecoder::LazAttribute::ReturnNumber:
363 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.returnNum() : p10.return_number );
364 break;
365 case QgsLazDecoder::LazAttribute::NumberOfReturns:
366 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.numReturns() : p10.number_of_returns_of_given_pulse );
367 break;
368 case QgsLazDecoder::LazAttribute::ScanDirectionFlag:
369 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.scanDirFlag() : p10.scan_direction_flag );
370 break;
371 case QgsLazDecoder::LazAttribute::EdgeOfFlightLine:
372 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.eofFlag() : p10.edge_of_flight_line );
373 break;
374 case QgsLazDecoder::LazAttribute::ScanAngleRank:
375 lazStoreToStream_<char>( dataBuffer, outputOffset, requestedAttribute.type, char( isLas14 ? p14.scanAngle() : p10.scan_angle_rank ) );
376 break;
377 case QgsLazDecoder::LazAttribute::UserData:
378 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.userData() : p10.user_data );
379 break;
380 case QgsLazDecoder::LazAttribute::PointSourceId:
381 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, isLas14 ? p14.pointSourceID() : p10.point_source_ID );
382 break;
383 case QgsLazDecoder::LazAttribute::GpsTime:
384 // lazperf internally stores gps value as int64 field, but in fact it is a double value
385 lazStoreToStream_<double>( dataBuffer, outputOffset, requestedAttribute.type,
386 isLas14 ? p14.gpsTime() : *reinterpret_cast<const double *>( reinterpret_cast<const void *>( &gps.value ) ) );
387 break;
388 case QgsLazDecoder::LazAttribute::Red:
389 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, rgb.r );
390 break;
391 case QgsLazDecoder::LazAttribute::Green:
392 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, rgb.g );
393 break;
394 case QgsLazDecoder::LazAttribute::Blue:
395 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, rgb.b );
396 break;
397 case QgsLazDecoder::LazAttribute::ScannerChannel:
398 lazStoreToStream_<char>( dataBuffer, outputOffset, requestedAttribute.type, char( p14.scannerChannel() ) );
399 break;
400 case QgsLazDecoder::LazAttribute::ClassificationFlags:
401 lazStoreToStream_<char>( dataBuffer, outputOffset, requestedAttribute.type, char( p14.classFlags() ) );
402 break;
403 case QgsLazDecoder::LazAttribute::NIR:
404 {
405 if ( lasPointFormat == 8 || lasPointFormat == 10 )
406 {
407 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, nir.val );
408 }
409 else
410 {
411 // just store 0 for missing attributes
412 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, 0 );
413 }
414 break;
415
416 }
417 case QgsLazDecoder::LazAttribute::ExtraBytes:
418 {
419 switch ( requestedAttribute.type )
420 {
422 lazStoreToStream_<char>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<char * >( &buf[requestedAttribute.offset] ) );
423 break;
425 lazStoreToStream_<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<unsigned char * >( &buf[requestedAttribute.offset] ) );
426 break;
428 lazStoreToStream_<qint16>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<qint16 * >( &buf[requestedAttribute.offset] ) );
429 break;
431 lazStoreToStream_<quint16>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<quint16 * >( &buf[requestedAttribute.offset] ) );
432 break;
434 lazStoreToStream_<qint32>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<qint32 * >( &buf[requestedAttribute.offset] ) );
435 break;
437 lazStoreToStream_<quint32>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<quint32 * >( &buf[requestedAttribute.offset] ) );
438 break;
440 lazStoreToStream_<qint64>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<qint64 * >( &buf[requestedAttribute.offset] ) );
441 break;
443 lazStoreToStream_<quint64>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<quint64 * >( &buf[requestedAttribute.offset] ) );
444 break;
446 lazStoreToStream_<float>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<float * >( &buf[requestedAttribute.offset] ) );
447 break;
449 lazStoreToStream_<double>( dataBuffer, outputOffset, requestedAttribute.type, *reinterpret_cast<double * >( &buf[requestedAttribute.offset] ) );
450 break;
451 }
452 }
453 break;
454 case QgsLazDecoder::LazAttribute::MissingOrUnknown:
455 // just store 0 for unknown/missing attributes
456 lazStoreToStream_<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, 0 );
457 break;
458 }
459
460 outputOffset += requestedAttribute.size;
461 }
462}
463
464template<typename FileType>
465QgsPointCloudBlock *decompressLaz_( FileType &file, const QgsPointCloudAttributeCollection &requestedAttributes, QgsPointCloudExpression &filterExpression, QgsRectangle &filterRect )
466{
467 if ( ! file.good() )
468 return nullptr;
469
470#ifdef QGISDEBUG
471 QElapsedTimer t;
472 t.start();
473#endif
474
475 // lazperf may throw exceptions
476 try
477 {
478 lazperf::reader::generic_file f( file );
479
480
481 // output file formats from entwine/untwine:
482 // - older versions write LAZ 1.2 files with point formats 0, 1, 2 or 3
483 // - newer versions write LAZ 1.4 files with point formats 6, 7 or 8
484
485 int lasPointFormat = f.header().pointFormat();
486 if ( lasPointFormat != 0 && lasPointFormat != 1 && lasPointFormat != 2 && lasPointFormat != 3 &&
487 lasPointFormat != 6 && lasPointFormat != 7 && lasPointFormat != 8 )
488 {
489 QgsDebugMsg( QStringLiteral( "Unexpected point format record (%1) - only 0, 1, 2, 3, 6, 7, 8 are supported" ).arg( lasPointFormat ) );
490 return nullptr;
491 }
492
493 const size_t count = f.header().point_count;
494 const QgsVector3D scale( f.header().scale.x, f.header().scale.y, f.header().scale.z );
495 const QgsVector3D offset( f.header().offset.x, f.header().offset.y, f.header().offset.z );
496
497 QByteArray bufArray( f.header().point_record_length, 0 );
498 char *buf = bufArray.data();
499
500 const size_t requestedPointRecordSize = requestedAttributes.pointRecordSize();
501 QByteArray data;
502 data.resize( requestedPointRecordSize * count );
503 char *dataBuffer = data.data();
504
505 std::size_t outputOffset = 0;
506
507 std::unique_ptr< QgsPointCloudBlock > block = std::make_unique< QgsPointCloudBlock >(
508 count,
509 requestedAttributes,
510 data, scale, offset
511 );
512
513 int skippedPoints = 0;
514 const bool filterIsValid = filterExpression.isValid();
515 if ( !filterExpression.prepare( block.get() ) && filterIsValid )
516 {
517 // skip processing if the expression cannot be prepared
518 block->setPointCount( 0 );
519 return block.release();
520 }
521
522 int xAttributeOffset, yAttributeOffset;
523 const QgsPointCloudAttribute *attributeX = nullptr;
524 const QgsPointCloudAttribute *attributeY = nullptr;
525 const bool hasFilterRect = !filterRect.isEmpty();
526 if ( hasFilterRect )
527 {
528 attributeX = requestedAttributes.find( QLatin1String( "X" ), xAttributeOffset );
529 attributeY = requestedAttributes.find( QLatin1String( "Y" ), yAttributeOffset );
530 filterRect.setXMinimum( ( filterRect.xMinimum() - offset.x() ) / scale.x() );
531 filterRect.setXMaximum( ( filterRect.xMaximum() - offset.x() ) / scale.x() );
532 filterRect.setYMinimum( ( filterRect.yMinimum() - offset.y() ) / scale.y() );
533 filterRect.setYMaximum( ( filterRect.yMaximum() - offset.y() ) / scale.y() );
534 }
535
536 std::vector<char> rawExtrabytes = f.vlrData( "LASF_Spec", 4 );
537 QVector<QgsLazInfo::ExtraBytesAttributeDetails> extrabyteAttributesDetails = QgsLazInfo::parseExtrabytes( rawExtrabytes.data(), rawExtrabytes.size(), f.header().point_record_length );
538 std::vector< QgsLazDecoder::RequestedAttributeDetails > requestedAttributeDetails = prepareRequestedAttributeDetails_( requestedAttributes, extrabyteAttributesDetails );
539
540 for ( size_t i = 0 ; i < count ; i ++ )
541 {
542 f.readPoint( buf ); // read the point out
543
544 decodePoint( buf, lasPointFormat, dataBuffer, outputOffset, requestedAttributeDetails );
545
546 // check if point needs to be filtered out
547 bool skipThisPoint = false;
548 if ( hasFilterRect && attributeX && attributeY )
549 {
550 const double x = attributeX->convertValueToDouble( dataBuffer + outputOffset - requestedPointRecordSize + xAttributeOffset );
551 const double y = attributeY->convertValueToDouble( dataBuffer + outputOffset - requestedPointRecordSize + yAttributeOffset );
552 if ( !filterRect.contains( x, y ) )
553 skipThisPoint = true;
554 }
555 if ( !skipThisPoint && filterIsValid )
556 {
557 // we're always evaluating the last written point in the buffer
558 double eval = filterExpression.evaluate( i - skippedPoints );
559 if ( !eval || std::isnan( eval ) )
560 skipThisPoint = true;
561 }
562 if ( skipThisPoint )
563 {
564 // if the point is filtered out, rewind the offset so the next point is written over it
565 outputOffset -= requestedPointRecordSize;
566 ++skippedPoints;
567 }
568 }
569
570#ifdef QGISDEBUG
571 QgsDebugMsgLevel( QStringLiteral( "LAZ-PERF Read through the points in %1 seconds." ).arg( t.elapsed() / 1000. ), 2 );
572#endif
573 block->setPointCount( count - skippedPoints );
574 return block.release();
575 }
576 catch ( std::exception &e )
577 {
578 QgsDebugMsg( "Error decompressing laz file: " + QString::fromLatin1( e.what() ) );
579 return nullptr;
580 }
581}
582
583QgsPointCloudBlock *QgsLazDecoder::decompressLaz( const QString &filename,
584 const QgsPointCloudAttributeCollection &requestedAttributes,
585 QgsPointCloudExpression &filterExpression, QgsRectangle &filterRect )
586{
587 std::ifstream file( toNativePath( filename ), std::ios::binary );
588
589 return decompressLaz_<std::ifstream>( file, requestedAttributes, filterExpression, filterRect );
590}
591
592QgsPointCloudBlock *QgsLazDecoder::decompressLaz( const QByteArray &byteArrayData,
593 const QgsPointCloudAttributeCollection &requestedAttributes,
594 QgsPointCloudExpression &filterExpression, QgsRectangle &filterRect )
595{
596 std::istringstream file( byteArrayData.toStdString() );
597 return decompressLaz_<std::istringstream>( file, requestedAttributes, filterExpression, filterRect );
598}
599
600QgsPointCloudBlock *QgsLazDecoder::decompressCopc( const QByteArray &data, QgsLazInfo &lazInfo, int32_t pointCount, const QgsPointCloudAttributeCollection &requestedAttributes, QgsPointCloudExpression &filterExpression, QgsRectangle &filterRect )
601{
602 // COPC only supports point formats 6, 7 and 8
603 int lasPointFormat = lazInfo.pointFormat();
604 if ( lasPointFormat != 6 && lasPointFormat != 7 && lasPointFormat != 8 )
605 {
606 QgsDebugMsg( QStringLiteral( "Unexpected point format record (%1) - only 6, 7, 8 are supported for COPC format" ).arg( lasPointFormat ) );
607 return nullptr;
608 }
609
610 std::unique_ptr<char []> decodedData( new char[ lazInfo.pointRecordLength() ] );
611
612 lazperf::reader::chunk_decompressor decompressor( lasPointFormat, lazInfo.extrabytesCount(), data.data() );
613
614 const size_t requestedPointRecordSize = requestedAttributes.pointRecordSize();
615 QByteArray blockData;
616 blockData.resize( requestedPointRecordSize * pointCount );
617 char *dataBuffer = blockData.data();
618
619 std::size_t outputOffset = 0;
620
621 QVector<QgsLazInfo::ExtraBytesAttributeDetails> extrabyteAttributesDetails = lazInfo.extrabytes();
622 std::vector< RequestedAttributeDetails > requestedAttributeDetails = prepareRequestedAttributeDetails_( requestedAttributes, extrabyteAttributesDetails );
623 std::unique_ptr< QgsPointCloudBlock > block = std::make_unique< QgsPointCloudBlock >(
624 pointCount, requestedAttributes,
625 blockData, lazInfo.scale(), lazInfo.offset()
626 );
627
628 int skippedPoints = 0;
629 const bool filterIsValid = filterExpression.isValid();
630 if ( !filterExpression.prepare( block.get() ) && filterIsValid )
631 {
632 // skip processing if the expression cannot be prepared
633 block->setPointCount( 0 );
634 return block.release();
635 }
636
637 int xAttributeOffset, yAttributeOffset;
638 const QgsPointCloudAttribute *attributeX = nullptr;
639 const QgsPointCloudAttribute *attributeY = nullptr;
640 const bool hasFilterRect = !filterRect.isEmpty();
641 if ( hasFilterRect )
642 {
643 attributeX = requestedAttributes.find( QLatin1String( "X" ), xAttributeOffset );
644 attributeY = requestedAttributes.find( QLatin1String( "Y" ), yAttributeOffset );
645 filterRect.setXMinimum( ( filterRect.xMinimum() - lazInfo.offset().x() ) / lazInfo.scale().x() );
646 filterRect.setXMaximum( ( filterRect.xMaximum() - lazInfo.offset().x() ) / lazInfo.scale().x() );
647 filterRect.setYMinimum( ( filterRect.yMinimum() - lazInfo.offset().y() ) / lazInfo.scale().y() );
648 filterRect.setYMaximum( ( filterRect.yMaximum() - lazInfo.offset().y() ) / lazInfo.scale().y() );
649 }
650 for ( int i = 0 ; i < pointCount; ++i )
651 {
652 decompressor.decompress( decodedData.get() );
653 char *buf = decodedData.get();
654
655 decodePoint( buf, lasPointFormat, dataBuffer, outputOffset, requestedAttributeDetails );
656
657 // check if point needs to be filtered out
658 bool skipThisPoint = false;
659
660 if ( hasFilterRect && attributeX && attributeY )
661 {
662 const double x = attributeX->convertValueToDouble( dataBuffer + outputOffset - requestedPointRecordSize + xAttributeOffset );
663 const double y = attributeY->convertValueToDouble( dataBuffer + outputOffset - requestedPointRecordSize + yAttributeOffset );
664 if ( !filterRect.contains( x, y ) )
665 skipThisPoint = true;
666 }
667 if ( !skipThisPoint && filterIsValid )
668 {
669 // we're always evaluating the last written point in the buffer
670 double eval = filterExpression.evaluate( i - skippedPoints );
671 if ( !eval || std::isnan( eval ) )
672 skipThisPoint = true;
673 }
674 if ( skipThisPoint )
675 {
676 // if the point is filtered out, rewind the offset so the next point is written over it
677 outputOffset -= requestedPointRecordSize;
678 ++skippedPoints;
679 }
680 }
681
682 block->setPointCount( pointCount - skippedPoints );
683 return block.release();
684}
685
686#if defined(_MSC_VER)
687std::wstring QgsLazDecoder::toNativePath( const QString &filename )
688{
689 std::wstring_convert< std::codecvt_utf8_utf16< wchar_t > > converter;
690 return converter.from_bytes( filename.toStdString() );
691}
692#else
693std::string QgsLazDecoder::toNativePath( const QString &filename )
694{
695 return filename.toStdString();
696}
697#endif
698
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.
Definition: qgslazinfo.cpp:206
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.
Base class for storing raw data from point cloud nodes.
A rectangle specified with double values.
Definition: qgsrectangle.h:42
double yMaximum() const SIP_HOLDGIL
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:193
double xMaximum() const SIP_HOLDGIL
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:183
double xMinimum() const SIP_HOLDGIL
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:188
double yMinimum() const SIP_HOLDGIL
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:198
void setYMinimum(double y) SIP_HOLDGIL
Set the minimum y value.
Definition: qgsrectangle.h:161
void setXMaximum(double x) SIP_HOLDGIL
Set the maximum x value.
Definition: qgsrectangle.h:156
void setXMinimum(double x) SIP_HOLDGIL
Set the minimum x value.
Definition: qgsrectangle.h:151
void setYMaximum(double y) SIP_HOLDGIL
Set the maximum y value.
Definition: qgsrectangle.h:166
bool isEmpty() const
Returns true if the rectangle is empty.
Definition: qgsrectangle.h:469
bool contains(const QgsRectangle &rect) const SIP_HOLDGIL
Returns true when rectangle contains other rectangle.
Definition: qgsrectangle.h:363
double y() const
Returns Y coordinate.
Definition: qgsvector3d.h:51
double x() const
Returns X coordinate.
Definition: qgsvector3d.h:49
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
#define QgsDebugMsg(str)
Definition: qgslogger.h:38