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