QGIS API Documentation 3.37.0-Master (fdefdf9c27f)
qgscopcpointcloudindex.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgscopcpointcloudindex.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
19
20#include <fstream>
21#include <QFile>
22#include <QtDebug>
23#include <QQueue>
24#include <QMutexLocker>
25#include <QJsonDocument>
26#include <QJsonObject>
27
28#include "qgseptdecoder.h"
29#include "qgslazdecoder.h"
33#include "qgslogger.h"
34#include "qgsfeedback.h"
35#include "qgsmessagelog.h"
36#include "qgspointcloudexpression.h"
37
38#include "lazperf/lazperf.hpp"
39#include "lazperf/readers.hpp"
40#include "lazperf/vlr.hpp"
41
43
44#define PROVIDER_KEY QStringLiteral( "copc" )
45#define PROVIDER_DESCRIPTION QStringLiteral( "COPC point cloud provider" )
46
47QgsCopcPointCloudIndex::QgsCopcPointCloudIndex() = default;
48
49QgsCopcPointCloudIndex::~QgsCopcPointCloudIndex() = default;
50
51std::unique_ptr<QgsPointCloudIndex> QgsCopcPointCloudIndex::clone() const
52{
53 QgsCopcPointCloudIndex *clone = new QgsCopcPointCloudIndex;
54 QMutexLocker locker( &mHierarchyMutex );
55 copyCommonProperties( clone );
56 return std::unique_ptr<QgsPointCloudIndex>( clone );
57}
58
59void QgsCopcPointCloudIndex::load( const QString &fileName )
60{
61 mUri = fileName;
62 mCopcFile.open( QgsLazDecoder::toNativePath( fileName ), std::ios::binary );
63
64 if ( !mCopcFile.is_open() || !mCopcFile.good() )
65 {
66 mError = tr( "Unable to open %1 for reading" ).arg( fileName );
67 mIsValid = false;
68 return;
69 }
70
71 mLazInfo.reset( new QgsLazInfo( QgsLazInfo::fromFile( mCopcFile ) ) );
72 mIsValid = mLazInfo->isValid();
73 if ( mIsValid )
74 {
75 mIsValid = loadSchema( *mLazInfo.get() );
76 }
77 if ( !mIsValid )
78 {
79 mError = tr( "Unable to recognize %1 as a LAZ file: \"%2\"" ).arg( fileName, mLazInfo->error() );
80 return;
81 }
82
83 loadHierarchy();
84}
85
86bool QgsCopcPointCloudIndex::loadSchema( QgsLazInfo &lazInfo )
87{
88 QByteArray copcInfoVlrData = lazInfo.vlrData( QStringLiteral( "copc" ), 1 );
89 if ( copcInfoVlrData.isEmpty() )
90 {
91 mError = tr( "Invalid COPC file" );
92 return false;
93 }
94 mCopcInfoVlr.fill( copcInfoVlrData.data(), copcInfoVlrData.size() );
95
96 mScale = lazInfo.scale();
97 mOffset = lazInfo.offset();
98
99 mOriginalMetadata = lazInfo.toMetadata();
100
101 QgsVector3D minCoords = lazInfo.minCoords();
102 QgsVector3D maxCoords = lazInfo.maxCoords();
103 mExtent.set( minCoords.x(), minCoords.y(), maxCoords.x(), maxCoords.y() );
104 mZMin = minCoords.z();
105 mZMax = maxCoords.z();
106
107 setAttributes( lazInfo.attributes() );
108
109 const double xmin = mCopcInfoVlr.center_x - mCopcInfoVlr.halfsize;
110 const double ymin = mCopcInfoVlr.center_y - mCopcInfoVlr.halfsize;
111 const double zmin = mCopcInfoVlr.center_z - mCopcInfoVlr.halfsize;
112 const double xmax = mCopcInfoVlr.center_x + mCopcInfoVlr.halfsize;
113 const double ymax = mCopcInfoVlr.center_y + mCopcInfoVlr.halfsize;
114 const double zmax = mCopcInfoVlr.center_z + mCopcInfoVlr.halfsize;
115
116 mRootBounds = QgsPointCloudDataBounds(
117 ( xmin - mOffset.x() ) / mScale.x(),
118 ( ymin - mOffset.y() ) / mScale.y(),
119 ( zmin - mOffset.z() ) / mScale.z(),
120 ( xmax - mOffset.x() ) / mScale.x(),
121 ( ymax - mOffset.y() ) / mScale.y(),
122 ( zmax - mOffset.z() ) / mScale.z()
123 );
124
125 double calculatedSpan = nodeMapExtent( root() ).width() / mCopcInfoVlr.spacing;
126 mSpan = calculatedSpan;
127
128#ifdef QGIS_DEBUG
129 double dx = xmax - xmin, dy = ymax - ymin, dz = zmax - zmin;
130 QgsDebugMsgLevel( QStringLiteral( "lvl0 node size in CRS units: %1 %2 %3" ).arg( dx ).arg( dy ).arg( dz ), 2 ); // all dims should be the same
131 QgsDebugMsgLevel( QStringLiteral( "res at lvl0 %1" ).arg( dx / mSpan ), 2 );
132 QgsDebugMsgLevel( QStringLiteral( "res at lvl1 %1" ).arg( dx / mSpan / 2 ), 2 );
133 QgsDebugMsgLevel( QStringLiteral( "res at lvl2 %1 with node size %2" ).arg( dx / mSpan / 4 ).arg( dx / 4 ), 2 );
134#endif
135
136 return true;
137}
138
139std::unique_ptr<QgsPointCloudBlock> QgsCopcPointCloudIndex::nodeData( const IndexedPointCloudNode &n, const QgsPointCloudRequest &request )
140{
141 if ( QgsPointCloudBlock *cached = getNodeDataFromCache( n, request ) )
142 {
143 return std::unique_ptr<QgsPointCloudBlock>( cached );
144 }
145
146 const bool found = fetchNodeHierarchy( n );
147 if ( !found )
148 return nullptr;
149 mHierarchyMutex.lock();
150 int pointCount = mHierarchy.value( n );
151 auto [blockOffset, blockSize] = mHierarchyNodePos.value( n );
152 mHierarchyMutex.unlock();
153
154 // we need to create a copy of the expression to pass to the decoder
155 // as the same QgsPointCloudExpression object mighgt be concurrently
156 // used on another thread, for example in a 3d view
157 QgsPointCloudExpression filterExpression = mFilterExpression;
158 QgsPointCloudAttributeCollection requestAttributes = request.attributes();
159 requestAttributes.extend( attributes(), filterExpression.referencedAttributes() );
160
161 QByteArray rawBlockData( blockSize, Qt::Initialization::Uninitialized );
162 std::ifstream file( QgsLazDecoder::toNativePath( mUri ), std::ios::binary );
163 file.seekg( blockOffset );
164 file.read( rawBlockData.data(), blockSize );
165 if ( !file )
166 {
167 QgsDebugError( QStringLiteral( "Could not read file %1" ).arg( mUri ) );
168 return nullptr;
169 }
170 QgsRectangle filterRect = request.filterRect();
171
172 std::unique_ptr<QgsPointCloudBlock> decoded = QgsLazDecoder::decompressCopc( rawBlockData, *mLazInfo.get(), pointCount, requestAttributes, filterExpression, filterRect );
173 storeNodeDataToCache( decoded.get(), n, request );
174 return decoded;
175}
176
177QgsPointCloudBlockRequest *QgsCopcPointCloudIndex::asyncNodeData( const IndexedPointCloudNode &n, const QgsPointCloudRequest &request )
178{
179 Q_UNUSED( n )
180 Q_UNUSED( request )
181 Q_ASSERT( false );
182 return nullptr; // unsupported
183}
184
186{
187 return mLazInfo->crs();
188}
189
190qint64 QgsCopcPointCloudIndex::pointCount() const
191{
192 return mLazInfo->pointCount();
193}
194
195bool QgsCopcPointCloudIndex::loadHierarchy()
196{
197 fetchHierarchyPage( mCopcInfoVlr.root_hier_offset, mCopcInfoVlr.root_hier_size );
198 return true;
199}
200
201bool QgsCopcPointCloudIndex::writeStatistics( QgsPointCloudStatistics &stats )
202{
203 if ( mLazInfo->version() != qMakePair<uint8_t, uint8_t>( 1, 4 ) )
204 {
205 // EVLR isn't supported in the first place
206 QgsMessageLog::logMessage( tr( "Can't write statistics to \"%1\": laz version != 1.4" ).arg( mUri ) );
207 return false;
208 }
209
210 QByteArray statisticsEvlrData = fetchCopcStatisticsEvlrData();
211 if ( !statisticsEvlrData.isEmpty() )
212 {
213 QgsMessageLog::logMessage( tr( "Can't write statistics to \"%1\": file already contains COPC statistics!" ).arg( mUri ) );
214 return false;
215 }
216
217 lazperf::evlr_header statsEvlrHeader;
218 statsEvlrHeader.user_id = "qgis";
219 statsEvlrHeader.record_id = 0;
220 statsEvlrHeader.description = "Contains calculated statistics";
221 QByteArray statsJson = stats.toStatisticsJson();
222 statsEvlrHeader.data_length = statsJson.size();
223
224 // Save the EVLRs to the end of the original file (while erasing the existing EVLRs in the file)
225 mCopcFile.close();
226 std::fstream copcFile;
227 copcFile.open( QgsLazDecoder::toNativePath( mUri ), std::ios_base::binary | std::iostream::in | std::iostream::out );
228 if ( copcFile.is_open() && copcFile.good() )
229 {
230 // Write the new number of EVLRs
231 lazperf::header14 header = mLazInfo->header();
232 header.evlr_count = header.evlr_count + 1;
233 copcFile.seekp( 0 );
234 header.write( copcFile );
235
236 // Append EVLR data to the end
237 copcFile.seekg( 0, std::ios::end );
238
239 statsEvlrHeader.write( copcFile );
240 copcFile.write( statsJson.data(), statsEvlrHeader.data_length );
241 }
242 else
243 {
244 QgsMessageLog::logMessage( tr( "Couldn't open COPC file \"%1\" to write statistics" ).arg( mUri ) );
245 return false;
246 }
247 copcFile.close();
248 mCopcFile.open( QgsLazDecoder::toNativePath( mUri ), std::ios::binary );
249 return true;
250}
251
252QgsPointCloudStatistics QgsCopcPointCloudIndex::readStatistics()
253{
254 QByteArray statisticsEvlrData = fetchCopcStatisticsEvlrData();
255
256 if ( statisticsEvlrData.isEmpty() )
257 {
259 }
260
261 return QgsPointCloudStatistics::fromStatisticsJson( statisticsEvlrData );
262}
263
264bool QgsCopcPointCloudIndex::isValid() const
265{
266 return mIsValid;
267}
268
269bool QgsCopcPointCloudIndex::fetchNodeHierarchy( const IndexedPointCloudNode &n ) const
270{
271 QMutexLocker locker( &mHierarchyMutex );
272
273 QVector<IndexedPointCloudNode> ancestors;
274 IndexedPointCloudNode foundRoot = n;
275 while ( !mHierarchy.contains( foundRoot ) )
276 {
277 ancestors.push_front( foundRoot );
278 foundRoot = foundRoot.parentNode();
279 }
280 ancestors.push_front( foundRoot );
281 for ( IndexedPointCloudNode n : ancestors )
282 {
283 auto hierarchyIt = mHierarchy.constFind( n );
284 if ( hierarchyIt == mHierarchy.constEnd() )
285 return false;
286 int nodesCount = *hierarchyIt;
287 if ( nodesCount < 0 )
288 {
289 auto hierarchyNodePos = mHierarchyNodePos.constFind( n );
290 mHierarchyMutex.unlock();
291 fetchHierarchyPage( hierarchyNodePos->first, hierarchyNodePos->second );
292 mHierarchyMutex.lock();
293 }
294 }
295 return mHierarchy.contains( n );
296}
297
298void QgsCopcPointCloudIndex::fetchHierarchyPage( uint64_t offset, uint64_t byteSize ) const
299{
300 mCopcFile.seekg( offset );
301 std::unique_ptr<char []> data( new char[ byteSize ] );
302 mCopcFile.read( data.get(), byteSize );
303
304 populateHierarchy( data.get(), byteSize );
305}
306
307void QgsCopcPointCloudIndex::populateHierarchy( const char *hierarchyPageData, uint64_t byteSize ) const
308{
309 struct CopcVoxelKey
310 {
311 int32_t level;
312 int32_t x;
313 int32_t y;
314 int32_t z;
315 };
316
317 struct CopcEntry
318 {
319 CopcVoxelKey key;
320 uint64_t offset;
321 int32_t byteSize;
322 int32_t pointCount;
323 };
324
325 QMutexLocker locker( &mHierarchyMutex );
326
327 for ( uint64_t i = 0; i < byteSize; i += sizeof( CopcEntry ) )
328 {
329 const CopcEntry *entry = reinterpret_cast<const CopcEntry *>( hierarchyPageData + i );
330 const IndexedPointCloudNode nodeId( entry->key.level, entry->key.x, entry->key.y, entry->key.z );
331 mHierarchy[nodeId] = entry->pointCount;
332 mHierarchyNodePos.insert( nodeId, QPair<uint64_t, int32_t>( entry->offset, entry->byteSize ) );
333 }
334}
335
336bool QgsCopcPointCloudIndex::hasNode( const IndexedPointCloudNode &n ) const
337{
338 return fetchNodeHierarchy( n );
339}
340
341QList<IndexedPointCloudNode> QgsCopcPointCloudIndex::nodeChildren( const IndexedPointCloudNode &n ) const
342{
343 fetchNodeHierarchy( n );
344
345 mHierarchyMutex.lock();
346
347 auto hierarchyIt = mHierarchy.constFind( n );
348 Q_ASSERT( hierarchyIt != mHierarchy.constEnd() );
349 QList<IndexedPointCloudNode> lst;
350 lst.reserve( 8 );
351 const int d = n.d() + 1;
352 const int x = n.x() * 2;
353 const int y = n.y() * 2;
354 const int z = n.z() * 2;
355 mHierarchyMutex.unlock();
356
357 for ( int i = 0; i < 8; ++i )
358 {
359 int dx = i & 1, dy = !!( i & 2 ), dz = !!( i & 4 );
360 const IndexedPointCloudNode n2( d, x + dx, y + dy, z + dz );
361 if ( fetchNodeHierarchy( n2 ) && mHierarchy[n] >= 0 )
362 lst.append( n2 );
363 }
364 return lst;
365}
366
367void QgsCopcPointCloudIndex::copyCommonProperties( QgsCopcPointCloudIndex *destination ) const
368{
370
371 // QgsCopcPointCloudIndex specific fields
372 destination->mIsValid = mIsValid;
373 destination->mUri = mUri;
374 destination->mCopcFile.open( QgsLazDecoder::toNativePath( mUri ), std::ios::binary );
375 destination->mCopcInfoVlr = mCopcInfoVlr;
376 destination->mHierarchyNodePos = mHierarchyNodePos;
377 destination->mOriginalMetadata = mOriginalMetadata;
378 destination->mLazInfo.reset( new QgsLazInfo( *mLazInfo ) );
379}
380
381QByteArray QgsCopcPointCloudIndex::fetchCopcStatisticsEvlrData()
382{
383 uint64_t offset = mLazInfo->firstEvlrOffset();
384 uint32_t evlrCount = mLazInfo->evlrCount();
385
386 QByteArray statisticsEvlrData;
387
388 for ( uint32_t i = 0; i < evlrCount; ++i )
389 {
390 lazperf::evlr_header header;
391 mCopcFile.seekg( offset );
392 char buffer[60];
393 mCopcFile.read( buffer, 60 );
394 header.fill( buffer, 60 );
395
396 // UserID: "qgis", record id: 0
397 if ( header.user_id == "qgis" && header.record_id == 0 )
398 {
399 statisticsEvlrData = QByteArray( header.data_length, Qt::Initialization::Uninitialized );
400 mCopcFile.read( statisticsEvlrData.data(), header.data_length );
401 break;
402 }
403
404 offset += 60 + header.data_length;
405 }
406
407 return statisticsEvlrData;
408}
409
Represents a indexed point cloud node in octree.
int y() const
Returns y.
int x() const
Returns x.
int d() const
Returns d.
IndexedPointCloudNode parentNode() const
Returns the parent of the node.
int z() const
Returns z.
This class represents a coordinate reference system (CRS).
Class for extracting information contained in LAZ file such as the public header block and variable l...
Definition: qgslazinfo.h:39
QgsVector3D maxCoords() const
Returns the maximum coordinate across X, Y and Z axis.
Definition: qgslazinfo.h:95
QgsPointCloudAttributeCollection attributes() const
Returns the list of attributes contained in the LAZ file.
Definition: qgslazinfo.h:120
QByteArray vlrData(QString userId, int recordId)
Returns the binary data of the variable length record with the user identifier userId and record iden...
Definition: qgslazinfo.cpp:143
QVariantMap toMetadata() const
Returns a map containing various metadata extracted from the LAZ file.
Definition: qgslazinfo.cpp:123
QgsVector3D scale() const
Returns the scale of the points coordinates.
Definition: qgslazinfo.h:77
static QgsLazInfo fromFile(std::ifstream &file)
Static function to create a QgsLazInfo class from a file.
Definition: qgslazinfo.cpp:279
QgsVector3D minCoords() const
Returns the minimum coordinate across X, Y and Z axis.
Definition: qgslazinfo.h:93
QgsVector3D offset() const
Returns the offset of the points coordinates.
Definition: qgslazinfo.h:79
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::MessageLevel::Warning, bool notifyUser=true)
Adds a message to the log instance (and creates it if necessary).
Collection of point cloud attributes.
void extend(const QgsPointCloudAttributeCollection &otherCollection, const QSet< QString > &matchingNames)
Adds specific missing attributes from another QgsPointCloudAttributeCollection.
Base class for handling loading QgsPointCloudBlock asynchronously.
Base class for storing raw data from point cloud nodes.
Represents packaged data bounds.
void copyCommonProperties(QgsPointCloudIndex *destination) const
Copies common properties to the destination index.
Point cloud data request.
QgsPointCloudAttributeCollection attributes() const
Returns attributes.
QgsRectangle filterRect() const
Returns the rectangle from which points will be taken, in point cloud's crs.
Class used to store statistics of a point cloud dataset.
static QgsPointCloudStatistics fromStatisticsJson(QByteArray stats)
Creates a statistics object from the JSON object stats.
QByteArray toStatisticsJson() const
Converts the current statistics object into JSON object.
A rectangle specified with double values.
Definition: qgsrectangle.h:42
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 z() const
Returns Z coordinate.
Definition: qgsvector3d.h:52
double x() const
Returns X coordinate.
Definition: qgsvector3d.h:48
void set(double x, double y, double z)
Sets vector coordinates.
Definition: qgsvector3d.h:73
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
#define QgsDebugError(str)
Definition: qgslogger.h:38
const QgsCoordinateReferenceSystem & crs