QGIS API Documentation  3.21.0-Master (909859188c)
qgseptdecoder.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgspointcloudrenderer.cpp
3  --------------------
4  begin : October 2020
5  copyright : (C) 2020 by Peter Petrik
6  email : zilolv 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 "qgseptdecoder.h"
19 #include "qgseptpointcloudindex.h"
20 #include "qgspointcloudattribute.h"
21 #include "qgsvector3d.h"
22 #include "qgsconfig.h"
23 #include "qgslogger.h"
24 
25 #include <QFile>
26 #include <QDir>
27 #include <iostream>
28 #include <memory>
29 #include <cstring>
30 #include <QTemporaryFile>
31 
32 #include <zstd.h>
33 
34 #include "laz-perf/io.hpp"
35 #include "laz-perf/common/common.hpp"
36 
38 
39 template <typename T>
40 bool _storeToStream( char *s, size_t position, QgsPointCloudAttribute::DataType type, T value )
41 {
42  switch ( type )
43  {
45  {
46  const char val = char( value );
47  s[position] = val;
48  break;
49  }
51  {
52  short val = short( value );
53  memcpy( s + position, reinterpret_cast<char * >( &val ), sizeof( short ) );
54  break;
55  }
56 
58  {
59  unsigned short val = static_cast< unsigned short>( value );
60  memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( unsigned short ) );
61  break;
62  }
63 
65  {
66  float val = float( value );
67  memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( float ) );
68  break;
69  }
71  {
72  qint32 val = qint32( value );
73  memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( qint32 ) );
74  break;
75  }
77  {
78  double val = double( value );
79  memcpy( s + position, reinterpret_cast< char * >( &val ), sizeof( double ) );
80  break;
81  }
82  }
83 
84  return true;
85 }
86 
87 bool __serialize( char *data, size_t outputPosition, QgsPointCloudAttribute::DataType outputType,
88  const char *input, QgsPointCloudAttribute::DataType inputType, int inputSize, size_t inputPosition )
89 {
90  if ( outputType == inputType )
91  {
92  memcpy( data + outputPosition, input + inputPosition, inputSize );
93  return true;
94  }
95 
96  switch ( inputType )
97  {
99  {
100  const char val = *( input + inputPosition );
101  return _storeToStream<char>( data, outputPosition, outputType, val );
102  }
104  {
105  const short val = *reinterpret_cast< const short * >( input + inputPosition );
106  return _storeToStream<short>( data, outputPosition, outputType, val );
107  }
109  {
110  const unsigned short val = *reinterpret_cast< const unsigned short * >( input + inputPosition );
111  return _storeToStream<unsigned short>( data, outputPosition, outputType, val );
112  }
114  {
115  const float val = *reinterpret_cast< const float * >( input + inputPosition );
116  return _storeToStream<float>( data, outputPosition, outputType, val );
117  }
119  {
120  const qint32 val = *reinterpret_cast<const qint32 * >( input + inputPosition );
121  return _storeToStream<qint32>( data, outputPosition, outputType, val );
122  }
124  {
125  const double val = *reinterpret_cast< const double * >( input + inputPosition );
126  return _storeToStream<double>( data, outputPosition, outputType, val );
127  }
128  }
129  return true;
130 }
131 
132 // //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
133 
134 QgsPointCloudBlock *_decompressBinary( const QByteArray &dataUncompressed, const QgsPointCloudAttributeCollection &attributes, const QgsPointCloudAttributeCollection &requestedAttributes, const QgsVector3D &scale, const QgsVector3D &offset )
135 {
136  const std::size_t pointRecordSize = attributes.pointRecordSize( );
137  const std::size_t requestedPointRecordSize = requestedAttributes.pointRecordSize();
138  const int count = dataUncompressed.size() / pointRecordSize;
139  QByteArray data;
140  data.resize( requestedPointRecordSize * count );
141  char *destinationBuffer = data.data();
142  const char *s = dataUncompressed.data();
143 
144  const QVector<QgsPointCloudAttribute> requestedAttributesVector = requestedAttributes.attributes();
145 
146  // calculate input attributes and offsets once in advance
147 
148  struct AttributeData
149  {
150  AttributeData( int inputOffset, int inputSize, QgsPointCloudAttribute::DataType inputType, int requestedSize, QgsPointCloudAttribute::DataType requestedType )
151  : inputOffset( inputOffset )
152  , inputSize( inputSize )
153  , inputType( inputType )
154  , requestedSize( requestedSize )
155  , requestedType( requestedType )
156  {}
157 
158  int inputOffset;
159  int inputSize;
161  int requestedSize;
162  QgsPointCloudAttribute::DataType requestedType;
163  };
164 
165  std::vector< AttributeData > attributeData;
166  attributeData.reserve( requestedAttributesVector.size() );
167  for ( const QgsPointCloudAttribute &requestedAttribute : requestedAttributesVector )
168  {
169  int inputAttributeOffset;
170  const QgsPointCloudAttribute *inputAttribute = attributes.find( requestedAttribute.name(), inputAttributeOffset );
171  if ( !inputAttribute )
172  {
173  return nullptr;
174  }
175  attributeData.emplace_back( AttributeData( inputAttributeOffset, inputAttribute->size(), inputAttribute->type(),
176  requestedAttribute.size(), requestedAttribute.type() ) );
177  }
178 
179  // now loop through points
180  size_t outputOffset = 0;
181  for ( int i = 0; i < count; ++i )
182  {
183  for ( const AttributeData &attribute : attributeData )
184  {
185  __serialize( destinationBuffer, outputOffset,
186  attribute.requestedType, s,
187  attribute.inputType, attribute.inputSize, i * pointRecordSize + attribute.inputOffset );
188 
189  outputOffset += attribute.requestedSize;
190  }
191  }
192  return new QgsPointCloudBlock(
193  count,
194  requestedAttributes,
195  data, scale, offset
196  );
197 }
198 
199 QgsPointCloudBlock *QgsEptDecoder::decompressBinary( const QString &filename, const QgsPointCloudAttributeCollection &attributes, const QgsPointCloudAttributeCollection &requestedAttributes, const QgsVector3D &scale, const QgsVector3D &offset )
200 {
201  if ( ! QFile::exists( filename ) )
202  return nullptr;
203 
204  QFile f( filename );
205  const bool r = f.open( QIODevice::ReadOnly );
206  if ( !r )
207  return nullptr;
208 
209  const QByteArray dataUncompressed = f.read( f.size() );
210  return _decompressBinary( dataUncompressed, attributes, requestedAttributes, scale, offset );
211 }
212 
213 QgsPointCloudBlock *QgsEptDecoder::decompressBinary( const QByteArray &data, const QgsPointCloudAttributeCollection &attributes, const QgsPointCloudAttributeCollection &requestedAttributes, const QgsVector3D &scale, const QgsVector3D &offset )
214 {
215  return _decompressBinary( data, attributes, requestedAttributes, scale, offset );
216 }
217 
218 /* *************************************************************************************** */
219 
220 QByteArray decompressZtdStream( const QByteArray &dataCompressed )
221 {
222  // NOTE: this is very primitive implementation because we expect the uncompressed
223  // data will be always less than 10 MB
224 
225  const int MAXSIZE = 10000000;
226  QByteArray dataUncompressed;
227  dataUncompressed.resize( MAXSIZE );
228 
229  ZSTD_DStream *strm = ZSTD_createDStream();
230  ZSTD_initDStream( strm );
231 
232  ZSTD_inBuffer m_inBuf;
233  m_inBuf.src = reinterpret_cast<const void *>( dataCompressed.constData() );
234  m_inBuf.size = dataCompressed.size();
235  m_inBuf.pos = 0;
236 
237  ZSTD_outBuffer outBuf { reinterpret_cast<void *>( dataUncompressed.data() ), MAXSIZE, 0 };
238  const size_t ret = ZSTD_decompressStream( strm, &outBuf, &m_inBuf );
239  Q_ASSERT( !ZSTD_isError( ret ) );
240  Q_ASSERT( outBuf.pos );
241  Q_ASSERT( outBuf.pos < outBuf.size );
242 
243  ZSTD_freeDStream( strm );
244  dataUncompressed.resize( outBuf.pos );
245  return dataUncompressed;
246 }
247 
248 QgsPointCloudBlock *QgsEptDecoder::decompressZStandard( const QString &filename, const QgsPointCloudAttributeCollection &attributes, const QgsPointCloudAttributeCollection &requestedAttributes, const QgsVector3D &scale, const QgsVector3D &offset )
249 {
250  if ( ! QFile::exists( filename ) )
251  return nullptr;
252 
253  QFile f( filename );
254  const bool r = f.open( QIODevice::ReadOnly );
255  if ( !r )
256  return nullptr;
257 
258  const QByteArray dataCompressed = f.readAll();
259  const QByteArray dataUncompressed = decompressZtdStream( dataCompressed );
260  return _decompressBinary( dataUncompressed, attributes, requestedAttributes, scale, offset );
261 }
262 
263 QgsPointCloudBlock *QgsEptDecoder::decompressZStandard( const QByteArray &data, const QgsPointCloudAttributeCollection &attributes, const QgsPointCloudAttributeCollection &requestedAttributes, const QgsVector3D &scale, const QgsVector3D &offset )
264 {
265  const QByteArray dataUncompressed = decompressZtdStream( data );
266  return _decompressBinary( dataUncompressed, attributes, requestedAttributes, scale, offset );
267 }
268 
269 /* *************************************************************************************** */
270 
271 template<typename FileType>
272 QgsPointCloudBlock *__decompressLaz( FileType &file, const QgsPointCloudAttributeCollection &attributes, const QgsPointCloudAttributeCollection &requestedAttributes, const QgsVector3D &_scale, const QgsVector3D &_offset )
273 {
274  Q_UNUSED( attributes );
275  Q_UNUSED( _scale );
276  Q_UNUSED( _offset );
277 
278  if ( ! file.good() )
279  return nullptr;
280 
281 #ifdef QGISDEBUG
282  const auto start = common::tick();
283 #endif
284 
285  laszip::io::reader::basic_file<FileType> f( file );
286 
287  const size_t count = f.get_header().point_count;
288  const QgsVector3D scale( f.get_header().scale.x, f.get_header().scale.y, f.get_header().scale.z );
289  const QgsVector3D offset( f.get_header().offset.x, f.get_header().offset.y, f.get_header().offset.z );
290 
291  QByteArray bufArray( f.get_header().point_record_length, 0 );
292  char *buf = bufArray.data();
293 
294  const size_t requestedPointRecordSize = requestedAttributes.pointRecordSize();
295  QByteArray data;
296  data.resize( requestedPointRecordSize * count );
297  char *dataBuffer = data.data();
298 
299  const QVector<QgsPointCloudAttribute> requestedAttributesVector = requestedAttributes.attributes();
300 
301  std::size_t outputOffset = 0;
302 
303  enum class LazAttribute
304  {
305  X,
306  Y,
307  Z,
308  Classification,
309  Intensity,
310  ReturnNumber,
311  NumberOfReturns,
312  ScanDirectionFlag,
313  EdgeOfFlightLine,
314  ScanAngleRank,
315  UserData,
316  PointSourceId,
317  Red,
318  Green,
319  Blue,
320  MissingOrUnknown
321  };
322 
323  struct RequestedAttributeDetails
324  {
325  RequestedAttributeDetails( LazAttribute attribute, QgsPointCloudAttribute::DataType type, int size )
326  : attribute( attribute )
327  , type( type )
328  , size( size )
329  {}
330 
331  LazAttribute attribute;
333  int size;
334  };
335 
336  std::vector< RequestedAttributeDetails > requestedAttributeDetails;
337  requestedAttributeDetails.reserve( requestedAttributesVector.size() );
338  for ( const QgsPointCloudAttribute &requestedAttribute : requestedAttributesVector )
339  {
340  if ( requestedAttribute.name().compare( QLatin1String( "X" ), Qt::CaseInsensitive ) == 0 )
341  {
342  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::X, requestedAttribute.type(), requestedAttribute.size() ) );
343  }
344  else if ( requestedAttribute.name().compare( QLatin1String( "Y" ), Qt::CaseInsensitive ) == 0 )
345  {
346  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::Y, requestedAttribute.type(), requestedAttribute.size() ) );
347  }
348  else if ( requestedAttribute.name().compare( QLatin1String( "Z" ), Qt::CaseInsensitive ) == 0 )
349  {
350  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::Z, requestedAttribute.type(), requestedAttribute.size() ) );
351  }
352  else if ( requestedAttribute.name().compare( QLatin1String( "Classification" ), Qt::CaseInsensitive ) == 0 )
353  {
354  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::Classification, requestedAttribute.type(), requestedAttribute.size() ) );
355  }
356  else if ( requestedAttribute.name().compare( QLatin1String( "Intensity" ), Qt::CaseInsensitive ) == 0 )
357  {
358  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::Intensity, requestedAttribute.type(), requestedAttribute.size() ) );
359  }
360  else if ( requestedAttribute.name().compare( QLatin1String( "ReturnNumber" ), Qt::CaseInsensitive ) == 0 )
361  {
362  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::ReturnNumber, requestedAttribute.type(), requestedAttribute.size() ) );
363  }
364  else if ( requestedAttribute.name().compare( QLatin1String( "NumberOfReturns" ), Qt::CaseInsensitive ) == 0 )
365  {
366  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::NumberOfReturns, requestedAttribute.type(), requestedAttribute.size() ) );
367  }
368  else if ( requestedAttribute.name().compare( QLatin1String( "ScanDirectionFlag" ), Qt::CaseInsensitive ) == 0 )
369  {
370  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::ScanDirectionFlag, requestedAttribute.type(), requestedAttribute.size() ) );
371  }
372  else if ( requestedAttribute.name().compare( QLatin1String( "EdgeOfFlightLine" ), Qt::CaseInsensitive ) == 0 )
373  {
374  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::EdgeOfFlightLine, requestedAttribute.type(), requestedAttribute.size() ) );
375  }
376  else if ( requestedAttribute.name().compare( QLatin1String( "ScanAngleRank" ), Qt::CaseInsensitive ) == 0 )
377  {
378  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::ScanAngleRank, requestedAttribute.type(), requestedAttribute.size() ) );
379  }
380  else if ( requestedAttribute.name().compare( QLatin1String( "UserData" ), Qt::CaseInsensitive ) == 0 )
381  {
382  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::UserData, requestedAttribute.type(), requestedAttribute.size() ) );
383  }
384  else if ( requestedAttribute.name().compare( QLatin1String( "PointSourceId" ), Qt::CaseInsensitive ) == 0 )
385  {
386  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::PointSourceId, requestedAttribute.type(), requestedAttribute.size() ) );
387  }
388  else if ( requestedAttribute.name().compare( QLatin1String( "Red" ), Qt::CaseInsensitive ) == 0 )
389  {
390  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::Red, requestedAttribute.type(), requestedAttribute.size() ) );
391  }
392  else if ( requestedAttribute.name().compare( QLatin1String( "Green" ), Qt::CaseInsensitive ) == 0 )
393  {
394  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::Green, requestedAttribute.type(), requestedAttribute.size() ) );
395  }
396  else if ( requestedAttribute.name().compare( QLatin1String( "Blue" ), Qt::CaseInsensitive ) == 0 )
397  {
398  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::Blue, requestedAttribute.type(), requestedAttribute.size() ) );
399  }
400  else
401  {
402  // 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
403  requestedAttributeDetails.emplace_back( RequestedAttributeDetails( LazAttribute::MissingOrUnknown, requestedAttribute.type(), requestedAttribute.size() ) );
404  }
405  }
406 
407  for ( size_t i = 0 ; i < count ; i ++ )
408  {
409  f.readPoint( buf ); // read the point out
410  const laszip::formats::las::point10 p = laszip::formats::packers<laszip::formats::las::point10>::unpack( buf );
411  const laszip::formats::las::rgb rgb = laszip::formats::packers<laszip::formats::las::rgb>::unpack( buf + sizeof( laszip::formats::las::point10 ) + sizeof( laszip::formats::las::gpstime ) );
412 
413  for ( const RequestedAttributeDetails &requestedAttribute : requestedAttributeDetails )
414  {
415  switch ( requestedAttribute.attribute )
416  {
417  case LazAttribute::X:
418  _storeToStream<qint32>( dataBuffer, outputOffset, requestedAttribute.type, p.x );
419  break;
420  case LazAttribute::Y:
421  _storeToStream<qint32>( dataBuffer, outputOffset, requestedAttribute.type, p.y );
422  break;
423  case LazAttribute::Z:
424  _storeToStream<qint32>( dataBuffer, outputOffset, requestedAttribute.type, p.z );
425  break;
426  case LazAttribute::Classification:
427  _storeToStream<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, p.classification );
428  break;
429  case LazAttribute::Intensity:
430  _storeToStream<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, p.intensity );
431  break;
432  case LazAttribute::ReturnNumber:
433  _storeToStream<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, p.return_number );
434  break;
435  case LazAttribute::NumberOfReturns:
436  _storeToStream<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, p.number_of_returns_of_given_pulse );
437  break;
438  case LazAttribute::ScanDirectionFlag:
439  _storeToStream<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, p.scan_direction_flag );
440  break;
441  case LazAttribute::EdgeOfFlightLine:
442  _storeToStream<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, p.edge_of_flight_line );
443  break;
444  case LazAttribute::ScanAngleRank:
445  _storeToStream<char>( dataBuffer, outputOffset, requestedAttribute.type, p.scan_angle_rank );
446  break;
447  case LazAttribute::UserData:
448  _storeToStream<unsigned char>( dataBuffer, outputOffset, requestedAttribute.type, p.user_data );
449  break;
450  case LazAttribute::PointSourceId:
451  _storeToStream<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, p.point_source_ID );
452  break;
453  case LazAttribute::Red:
454  _storeToStream<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, rgb.r );
455  break;
456  case LazAttribute::Green:
457  _storeToStream<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, rgb.g );
458  break;
459  case LazAttribute::Blue:
460  _storeToStream<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, rgb.b );
461  break;
462  case LazAttribute::MissingOrUnknown:
463  // just store 0 for unknown/missing attributes
464  _storeToStream<unsigned short>( dataBuffer, outputOffset, requestedAttribute.type, 0 );
465  break;
466  }
467 
468  outputOffset += requestedAttribute.size;
469  }
470  }
471 
472 #ifdef QGISDEBUG
473  const float t = common::since( start );
474  QgsDebugMsgLevel( QStringLiteral( "LAZ-PERF Read through the points in %1 seconds." ).arg( t ), 2 );
475 #endif
477  count,
478  requestedAttributes,
479  data, scale, offset
480  );
481  return block;
482 }
483 
484 QgsPointCloudBlock *QgsEptDecoder::decompressLaz( const QString &filename,
485  const QgsPointCloudAttributeCollection &attributes,
486  const QgsPointCloudAttributeCollection &requestedAttributes,
487  const QgsVector3D &scale, const QgsVector3D &offset )
488 {
489  const QByteArray arr = filename.toUtf8();
490  std::ifstream file( arr.constData(), std::ios::binary );
491 
492  return __decompressLaz<std::ifstream>( file, attributes, requestedAttributes, scale, offset );
493 }
494 
495 QgsPointCloudBlock *QgsEptDecoder::decompressLaz( const QByteArray &byteArrayData,
496  const QgsPointCloudAttributeCollection &attributes,
497  const QgsPointCloudAttributeCollection &requestedAttributes,
498  const QgsVector3D &scale, const QgsVector3D &offset )
499 {
500  std::istringstream file( byteArrayData.toStdString() );
501  return __decompressLaz<std::istringstream>( file, attributes, requestedAttributes, scale, offset );
502 }
503 
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.
int size() const
Returns size of the attribute in bytes.
DataType type() const
Returns the data type.
Base class for storing raw data from point cloud nodes.
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39