QGIS API Documentation 3.29.0-Master (ade4f0cf0f)
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 mFileName = 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
139QgsPointCloudBlock *QgsCopcPointCloudIndex::nodeData( const IndexedPointCloudNode &n, const QgsPointCloudRequest &request )
140{
141 const bool found = fetchNodeHierarchy( n );
142 if ( !found )
143 return nullptr;
144 mHierarchyMutex.lock();
145 int pointCount = mHierarchy.value( n );
146 auto [blockOffset, blockSize] = mHierarchyNodePos.value( n );
147 mHierarchyMutex.unlock();
148
149 // we need to create a copy of the expression to pass to the decoder
150 // as the same QgsPointCloudExpression object mighgt be concurrently
151 // used on another thread, for example in a 3d view
152 QgsPointCloudExpression filterExpression = mFilterExpression;
153 QgsPointCloudAttributeCollection requestAttributes = request.attributes();
154 requestAttributes.extend( attributes(), filterExpression.referencedAttributes() );
155
156 QByteArray rawBlockData( blockSize, Qt::Initialization::Uninitialized );
157 std::ifstream file( QgsLazDecoder::toNativePath( mFileName ), std::ios::binary );
158 file.seekg( blockOffset );
159 file.read( rawBlockData.data(), blockSize );
160 if ( !file )
161 {
162 QgsDebugMsg( QStringLiteral( "Could not read file %1" ).arg( mFileName ) );
163 return nullptr;
164 }
165
166 return QgsLazDecoder::decompressCopc( rawBlockData, *mLazInfo.get(), pointCount, requestAttributes, filterExpression );
167}
168
169QgsPointCloudBlockRequest *QgsCopcPointCloudIndex::asyncNodeData( const IndexedPointCloudNode &n, const QgsPointCloudRequest &request )
170{
171 Q_UNUSED( n )
172 Q_UNUSED( request )
173 Q_ASSERT( false );
174 return nullptr; // unsupported
175}
176
178{
179 return mLazInfo->crs();
180}
181
182qint64 QgsCopcPointCloudIndex::pointCount() const
183{
184 return mLazInfo->pointCount();
185}
186
187bool QgsCopcPointCloudIndex::loadHierarchy()
188{
189 QMutexLocker locker( &mHierarchyMutex );
190 fetchHierarchyPage( mCopcInfoVlr.root_hier_offset, mCopcInfoVlr.root_hier_size );
191 return true;
192}
193
194bool QgsCopcPointCloudIndex::writeStatistics( QgsPointCloudStatistics &stats )
195{
196 if ( mLazInfo->version() != qMakePair<uint8_t, uint8_t>( 1, 4 ) )
197 {
198 // EVLR isn't supported in the first place
199 QgsMessageLog::logMessage( tr( "Can't write statistics to \"%1\": laz version != 1.4" ).arg( mFileName ) );
200 return false;
201 }
202
203 QByteArray statisticsEvlrData = fetchCopcStatisticsEvlrData();
204 if ( !statisticsEvlrData.isEmpty() )
205 {
206 QgsMessageLog::logMessage( tr( "Can't write statistics to \"%1\": file already contains COPC statistics!" ).arg( mFileName ) );
207 return false;
208 }
209
210 lazperf::evlr_header statsEvlrHeader;
211 statsEvlrHeader.user_id = "qgis";
212 statsEvlrHeader.record_id = 0;
213 statsEvlrHeader.description = "Contains calculated statistics";
214 QByteArray statsJson = stats.toStatisticsJson();
215 statsEvlrHeader.data_length = statsJson.size();
216
217 // Save the EVLRs to the end of the original file (while erasing the exisitng EVLRs in the file)
218 mCopcFile.close();
219 std::fstream copcFile;
220 copcFile.open( QgsLazDecoder::toNativePath( mFileName ), std::ios_base::binary | std::iostream::in | std::iostream::out );
221 if ( copcFile.is_open() && copcFile.good() )
222 {
223 // Write the new number of EVLRs
224 lazperf::header14 header = mLazInfo->header();
225 header.evlr_count = header.evlr_count + 1;
226 copcFile.seekp( 0 );
227 header.write( copcFile );
228
229 // Append EVLR data to the end
230 copcFile.seekg( 0, std::ios::end );
231
232 statsEvlrHeader.write( copcFile );
233 copcFile.write( statsJson.data(), statsEvlrHeader.data_length );
234 }
235 else
236 {
237 QgsMessageLog::logMessage( tr( "Couldn't open COPC file \"%1\" to write statistics" ).arg( mFileName ) );
238 return false;
239 }
240 copcFile.close();
241 mCopcFile.open( QgsLazDecoder::toNativePath( mFileName ), std::ios::binary );
242 return true;
243}
244
245QgsPointCloudStatistics QgsCopcPointCloudIndex::readStatistics()
246{
247 QByteArray statisticsEvlrData = fetchCopcStatisticsEvlrData();
248
249 if ( statisticsEvlrData.isEmpty() )
250 {
252 }
253
254 return QgsPointCloudStatistics::fromStatisticsJson( statisticsEvlrData );
255}
256
257bool QgsCopcPointCloudIndex::isValid() const
258{
259 return mIsValid;
260}
261
262bool QgsCopcPointCloudIndex::fetchNodeHierarchy( const IndexedPointCloudNode &n ) const
263{
264 QMutexLocker locker( &mHierarchyMutex );
265
266 QVector<IndexedPointCloudNode> ancestors;
267 IndexedPointCloudNode foundRoot = n;
268 while ( !mHierarchy.contains( foundRoot ) )
269 {
270 ancestors.push_front( foundRoot );
271 foundRoot = foundRoot.parentNode();
272 }
273 ancestors.push_front( foundRoot );
274 for ( IndexedPointCloudNode n : ancestors )
275 {
276 auto hierarchyIt = mHierarchy.constFind( n );
277 if ( hierarchyIt == mHierarchy.constEnd() )
278 return false;
279 int nodesCount = *hierarchyIt;
280 if ( nodesCount < 0 )
281 {
282 auto hierarchyNodePos = mHierarchyNodePos.constFind( n );
283 fetchHierarchyPage( hierarchyNodePos->first, hierarchyNodePos->second );
284 }
285 }
286 return true;
287}
288
289void QgsCopcPointCloudIndex::fetchHierarchyPage( uint64_t offset, uint64_t byteSize ) const
290{
291 mCopcFile.seekg( offset );
292 std::unique_ptr<char []> data( new char[ byteSize ] );
293 mCopcFile.read( data.get(), byteSize );
294
295 struct CopcVoxelKey
296 {
297 int32_t level;
298 int32_t x;
299 int32_t y;
300 int32_t z;
301 };
302
303 struct CopcEntry
304 {
305 CopcVoxelKey key;
306 uint64_t offset;
307 int32_t byteSize;
308 int32_t pointCount;
309 };
310
311 for ( uint64_t i = 0; i < byteSize; i += sizeof( CopcEntry ) )
312 {
313 CopcEntry *entry = reinterpret_cast<CopcEntry *>( data.get() + i );
314 const IndexedPointCloudNode nodeId( entry->key.level, entry->key.x, entry->key.y, entry->key.z );
315 mHierarchy[nodeId] = entry->pointCount;
316 mHierarchyNodePos.insert( nodeId, QPair<uint64_t, int32_t>( entry->offset, entry->byteSize ) );
317 }
318}
319
320bool QgsCopcPointCloudIndex::hasNode( const IndexedPointCloudNode &n ) const
321{
322 fetchNodeHierarchy( n );
323 mHierarchyMutex.lock();
324
325 auto it = mHierarchy.constFind( n );
326 const bool found = it != mHierarchy.constEnd() && ( *it ) >= 0;
327 mHierarchyMutex.unlock();
328 return found;
329}
330
331QList<IndexedPointCloudNode> QgsCopcPointCloudIndex::nodeChildren( const IndexedPointCloudNode &n ) const
332{
333 fetchNodeHierarchy( n );
334
335 mHierarchyMutex.lock();
336
337 auto hierarchyIt = mHierarchy.constFind( n );
338 Q_ASSERT( hierarchyIt != mHierarchy.constEnd() );
339 QList<IndexedPointCloudNode> lst;
340 lst.reserve( 8 );
341 const int d = n.d() + 1;
342 const int x = n.x() * 2;
343 const int y = n.y() * 2;
344 const int z = n.z() * 2;
345 mHierarchyMutex.unlock();
346
347 for ( int i = 0; i < 8; ++i )
348 {
349 int dx = i & 1, dy = !!( i & 2 ), dz = !!( i & 4 );
350 const IndexedPointCloudNode n2( d, x + dx, y + dy, z + dz );
351 if ( fetchNodeHierarchy( n2 ) && mHierarchy[n] >= 0 )
352 lst.append( n2 );
353 }
354 return lst;
355}
356
357void QgsCopcPointCloudIndex::copyCommonProperties( QgsCopcPointCloudIndex *destination ) const
358{
360
361 // QgsCopcPointCloudIndex specific fields
362 destination->mIsValid = mIsValid;
363 destination->mFileName = mFileName;
364 destination->mCopcFile.open( QgsLazDecoder::toNativePath( mFileName ), std::ios::binary );
365 destination->mCopcInfoVlr = mCopcInfoVlr;
366 destination->mHierarchyNodePos = mHierarchyNodePos;
367 destination->mOriginalMetadata = mOriginalMetadata;
368 destination->mLazInfo.reset( new QgsLazInfo( *mLazInfo ) );
369}
370
371QByteArray QgsCopcPointCloudIndex::fetchCopcStatisticsEvlrData()
372{
373 uint64_t offset = mLazInfo->firstEvlrOffset();
374 uint32_t evlrCount = mLazInfo->evlrCount();
375
376 QByteArray statisticsEvlrData;
377
378 for ( uint32_t i = 0; i < evlrCount; ++i )
379 {
380 lazperf::evlr_header header;
381 mCopcFile.seekg( offset );
382 char buffer[60];
383 mCopcFile.read( buffer, 60 );
384 header.fill( buffer, 60 );
385
386 // UserID: "qgis", record id: 0
387 if ( header.user_id == "qgis" && header.record_id == 0 )
388 {
389 statisticsEvlrData = QByteArray( header.data_length, Qt::Initialization::Uninitialized );
390 mCopcFile.read( statisticsEvlrData.data(), header.data_length );
391 break;
392 }
393
394 offset += 60 + header.data_length;
395 }
396
397 return statisticsEvlrData;
398}
399
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:141
QVariantMap toMetadata() const
Returns a map containing various metadata extracted from the LAZ file.
Definition: qgslazinfo.cpp:121
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:274
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.
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.
double y() const
Returns Y coordinate.
Definition: qgsvector3d.h:51
double z() const
Returns Z coordinate.
Definition: qgsvector3d.h:53
double x() const
Returns X coordinate.
Definition: qgsvector3d.h:49
void set(double x, double y, double z)
Sets vector coordinates.
Definition: qgsvector3d.h:56
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
const QgsCoordinateReferenceSystem & crs