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