24#include <QMutexLocker>
25#include <QJsonDocument>
27#include <qnamespace.h>
42#include "qgspointcloudexpression.h"
44#include "lazperf/vlr.hpp"
49#define PROVIDER_KEY QStringLiteral( "copc" )
50#define PROVIDER_DESCRIPTION QStringLiteral( "COPC point cloud provider" )
52QgsCopcPointCloudIndex::QgsCopcPointCloudIndex() =
default;
54QgsCopcPointCloudIndex::~QgsCopcPointCloudIndex() =
default;
56void QgsCopcPointCloudIndex::load(
const QString &urlString )
60 if ( url.isValid() && ( url.scheme() ==
"http" || url.scheme() ==
"https" ) )
65 mUri = url.toString();
71 mCopcFile.open( QgsLazDecoder::toNativePath( urlString ), std::ios::binary );
72 if ( mCopcFile.fail() )
74 mError = QObject::tr(
"Unable to open %1 for reading" ).arg( urlString );
81 mIsValid = mLazInfo->isValid() && loadSchema( *mLazInfo.get() ) && loadHierarchy();
84 mError = QObject::tr(
"Unable to recognize %1 as a LAZ file: \"%2\"" ).arg( urlString, mLazInfo->error() );
88bool QgsCopcPointCloudIndex::loadSchema(
QgsLazInfo &lazInfo )
90 QByteArray copcInfoVlrData = lazInfo.
vlrData( QStringLiteral(
"copc" ), 1 );
91 if ( copcInfoVlrData.isEmpty() )
93 mError = QObject::tr(
"Invalid COPC file" );
96 mCopcInfoVlr.fill( copcInfoVlrData.data(), copcInfoVlrData.size() );
98 mScale = lazInfo.
scale();
99 mOffset = lazInfo.
offset();
105 mExtent.
set( minCoords.
x(), minCoords.
y(), maxCoords.
x(), maxCoords.
y() );
106 mZMin = minCoords.
z();
107 mZMax = maxCoords.
z();
111 const double xmin = mCopcInfoVlr.center_x - mCopcInfoVlr.halfsize;
112 const double ymin = mCopcInfoVlr.center_y - mCopcInfoVlr.halfsize;
113 const double zmin = mCopcInfoVlr.center_z - mCopcInfoVlr.halfsize;
114 const double xmax = mCopcInfoVlr.center_x + mCopcInfoVlr.halfsize;
115 const double ymax = mCopcInfoVlr.center_y + mCopcInfoVlr.halfsize;
116 const double zmax = mCopcInfoVlr.center_z + mCopcInfoVlr.halfsize;
118 mRootBounds =
QgsBox3D( xmin, ymin, zmin, xmax, ymax, zmax );
121 mSpan = mRootBounds.width() / mCopcInfoVlr.spacing;
124 double dx = xmax - xmin, dy = ymax - ymin, dz = zmax - zmin;
125 QgsDebugMsgLevel( QStringLiteral(
"lvl0 node size in CRS units: %1 %2 %3" ).arg( dx ).arg( dy ).arg( dz ), 2 );
126 QgsDebugMsgLevel( QStringLiteral(
"res at lvl0 %1" ).arg( dx / mSpan ), 2 );
127 QgsDebugMsgLevel( QStringLiteral(
"res at lvl1 %1" ).arg( dx / mSpan / 2 ), 2 );
128 QgsDebugMsgLevel( QStringLiteral(
"res at lvl2 %1 with node size %2" ).arg( dx / mSpan / 4 ).arg( dx / 4 ), 2 );
138 return std::unique_ptr<QgsPointCloudBlock>( cached );
141 std::unique_ptr<QgsPointCloudBlock> block;
144 QByteArray rawBlockData = rawNodeData( n );
145 if ( rawBlockData.isEmpty() )
148 mHierarchyMutex.lock();
149 auto pointCount = mHierarchy.value( n );
150 mHierarchyMutex.unlock();
155 QgsPointCloudExpression filterExpression = request.
ignoreIndexFilterEnabled() ? QgsPointCloudExpression() : mFilterExpression;
157 requestAttributes.
extend( attributes(), filterExpression.referencedAttributes() );
161 block = QgsLazDecoder::decompressCopc( rawBlockData, *mLazInfo.get(), pointCount, requestAttributes, filterExpression, filterRect );
166 std::unique_ptr<QgsPointCloudBlockRequest> blockRequest( asyncNodeData( n, request ) );
174 block = blockRequest->takeBlock();
177 QgsDebugError( QStringLiteral(
"Error downloading node %1 data, error : %2 " ).arg( n.
toString(), blockRequest->errorStr() ) );
180 storeNodeDataToCache( block.get(), n, request );
191 scale(), offset(), mFilterExpression, request.
filterRect() );
194 if ( !fetchNodeHierarchy( n ) )
196 QMutexLocker locker( &mHierarchyMutex );
201 QgsPointCloudExpression filterExpression = request.
ignoreIndexFilterEnabled() ? QgsPointCloudExpression() : mFilterExpression;
203 requestAttributes.
extend( attributes(), filterExpression.referencedAttributes() );
204 auto [ blockOffset, blockSize ] = mHierarchyNodePos.value( n );
205 int pointCount = mHierarchy.value( n );
208 scale(), offset(), filterExpression, request.
filterRect(),
209 blockOffset, blockSize, pointCount, *mLazInfo.get() );
215 const bool found = fetchNodeHierarchy( n );
218 mHierarchyMutex.lock();
219 auto [blockOffset, blockSize] = mHierarchyNodePos.value( n );
220 mHierarchyMutex.unlock();
225 QByteArray rawBlockData( blockSize, Qt::Initialization::Uninitialized );
226 std::ifstream file( QgsLazDecoder::toNativePath( mUri ), std::ios::binary );
227 file.seekg( blockOffset );
228 file.read( rawBlockData.data(), blockSize );
231 QgsDebugError( QStringLiteral(
"Could not read file %1" ).arg( mUri ) );
237 return readRange( blockOffset, blockSize );
242 return mLazInfo->crs();
245qint64 QgsCopcPointCloudIndex::pointCount()
const
247 return mLazInfo->pointCount();
250bool QgsCopcPointCloudIndex::loadHierarchy()
const
252 fetchHierarchyPage( mCopcInfoVlr.root_hier_offset, mCopcInfoVlr.root_hier_size );
264 if ( mLazInfo->version() != qMakePair<uint8_t, uint8_t>( 1, 4 ) )
271 QByteArray statisticsEvlrData = fetchCopcStatisticsEvlrData();
272 if ( !statisticsEvlrData.isEmpty() )
274 QgsMessageLog::logMessage( QObject::tr(
"Can't write statistics to \"%1\": file already contains COPC statistics!" ).arg( mUri ) );
278 lazperf::evlr_header statsEvlrHeader;
279 statsEvlrHeader.user_id =
"qgis";
280 statsEvlrHeader.record_id = 0;
281 statsEvlrHeader.description =
"Contains calculated statistics";
283 statsEvlrHeader.data_length = statsJson.size();
286 QMutexLocker locker( &mFileMutex );
288 std::fstream copcFile;
289 copcFile.open( QgsLazDecoder::toNativePath( mUri ), std::ios_base::binary | std::iostream::in | std::iostream::out );
290 if ( copcFile.is_open() && copcFile.good() )
293 lazperf::header14 header = mLazInfo->header();
294 header.evlr_count = header.evlr_count + 1;
296 header.write( copcFile );
299 copcFile.seekg( 0, std::ios::end );
301 statsEvlrHeader.write( copcFile );
302 copcFile.write( statsJson.data(), statsEvlrHeader.data_length );
310 mCopcFile.open( QgsLazDecoder::toNativePath( mUri ), std::ios::binary );
318 const QByteArray statisticsEvlrData = fetchCopcStatisticsEvlrData();
319 if ( statisticsEvlrData.isEmpty() )
328bool QgsCopcPointCloudIndex::isValid()
const
335 QMutexLocker locker( &mHierarchyMutex );
337 QVector<QgsPointCloudNodeId> ancestors;
339 while ( !mHierarchy.contains( foundRoot ) )
341 ancestors.push_front( foundRoot );
344 ancestors.push_front( foundRoot );
347 auto hierarchyIt = mHierarchy.constFind( n );
348 if ( hierarchyIt == mHierarchy.constEnd() )
350 int nodesCount = *hierarchyIt;
351 if ( nodesCount < 0 )
353 auto hierarchyNodePos = mHierarchyNodePos.constFind( n );
354 mHierarchyMutex.unlock();
355 fetchHierarchyPage( hierarchyNodePos->first, hierarchyNodePos->second );
356 mHierarchyMutex.lock();
359 return mHierarchy.contains( n );
362void QgsCopcPointCloudIndex::fetchHierarchyPage( uint64_t offset, uint64_t byteSize )
const
364 Q_ASSERT( byteSize > 0 );
366 QByteArray data = readRange( offset, byteSize );
367 if ( data.isEmpty() )
370 populateHierarchy( data.constData(), byteSize );
373void QgsCopcPointCloudIndex::populateHierarchy(
const char *hierarchyPageData, uint64_t byteSize )
const
391 QMutexLocker locker( &mHierarchyMutex );
393 for ( uint64_t i = 0; i < byteSize; i +=
sizeof( CopcEntry ) )
395 const CopcEntry *entry =
reinterpret_cast<const CopcEntry *
>( hierarchyPageData + i );
396 const QgsPointCloudNodeId nodeId( entry->key.level, entry->key.x, entry->key.y, entry->key.z );
397 mHierarchy[nodeId] = entry->pointCount;
398 mHierarchyNodePos.insert( nodeId, QPair<uint64_t, int32_t>( entry->offset, entry->byteSize ) );
404 return fetchNodeHierarchy( n );
409 bool nodeFound = fetchNodeHierarchy(
id );
410 Q_ASSERT( nodeFound );
414 QMutexLocker locker( &mHierarchyMutex );
415 pointCount = mHierarchy.value(
id, -1 );
418 QList<QgsPointCloudNodeId> children;
419 children.reserve( 8 );
420 const int d =
id.d() + 1;
421 const int x =
id.x() * 2;
422 const int y =
id.y() * 2;
423 const int z =
id.z() * 2;
425 for (
int i = 0; i < 8; ++i )
427 int dx = i & 1, dy = !!( i & 2 ), dz = !!( i & 4 );
429 bool found = fetchNodeHierarchy( n2 );
431 QMutexLocker locker( &mHierarchyMutex );
432 if ( found && mHierarchy[
id] >= 0 )
433 children.append( n2 );
441QByteArray QgsCopcPointCloudIndex::readRange( uint64_t offset, uint64_t length )
const
445 QMutexLocker locker( &mFileMutex );
447 QByteArray buffer( length, Qt::Initialization::Uninitialized );
448 mCopcFile.seekg( offset );
449 mCopcFile.read( buffer.data(), length );
450 if ( mCopcFile.eof() )
451 QgsDebugError( QStringLiteral(
"Read past end of file (path %1 offset %2 length %3)" ).arg( mUri ).arg( offset ).arg( length ) );
453 QgsDebugError( QStringLiteral(
"Error reading %1" ).arg( mUri ) );
458 QNetworkRequest nr = QNetworkRequest( QUrl( mUri ) );
460 nr.setAttribute( QNetworkRequest::CacheLoadControlAttribute, QNetworkRequest::PreferCache );
461 nr.setAttribute( QNetworkRequest::CacheSaveControlAttribute,
true );
462 QByteArray queryRange = QStringLiteral(
"bytes=%1-%2" ).arg( offset ).arg( offset + length - 1 ).toLocal8Bit();
463 nr.setRawHeader(
"Range", queryRange );
471 if ( reply->error() != QNetworkReply::NoError )
473 QgsDebugError( QStringLiteral(
"Request failed: %1 (offset %1 length %2)" ).arg( mUri ).arg( offset ).arg( length ) );
477 return reply->data();
481QByteArray QgsCopcPointCloudIndex::fetchCopcStatisticsEvlrData()
const
483 uint64_t offset = mLazInfo->firstEvlrOffset();
484 uint32_t evlrCount = mLazInfo->evlrCount();
486 QByteArray statisticsEvlrData;
488 for ( uint32_t i = 0; i < evlrCount; ++i )
490 lazperf::evlr_header header;
492 QByteArray buffer = readRange( offset, 60 );
493 header.fill( buffer.data(), buffer.size() );
495 if ( header.user_id ==
"qgis" && header.record_id == 0 )
497 statisticsEvlrData = readRange( offset + 60, header.data_length );
501 offset += 60 + header.data_length;
504 return statisticsEvlrData;
507void QgsCopcPointCloudIndex::reset()
525 mOriginalMetadata.clear();
528 mHierarchyNodePos.clear();
531QVariantMap QgsCopcPointCloudIndex::extraMetadata()
const
535 { QStringLiteral(
"CopcGpsTimeFlag" ), mLazInfo.get()->header().global_encoding & 1 },
@ Local
Local means the source is a local file on the machine.
@ Remote
Remote means it's loaded through a protocol like HTTP.
virtual QgsPointCloudStatistics metadataStatistics() const
Returns the object containing the statistics metadata extracted from the dataset.
static QgsTileDownloadManager * tileDownloadManager()
Returns the application's tile download manager, used for download of map tiles when rendering.
A 3-dimensional box composed of x, y, z coordinates.
double width() const
Returns the width of the box.
Handles a QgsPointCloudBlockRequest using existing cached QgsPointCloudBlock.
Represents a coordinate reference system (CRS).
Base class for handling loading QgsPointCloudBlock asynchronously from a remote COPC dataset.
Extracts information contained in a LAZ file, such as the public header block and variable length rec...
QgsVector3D maxCoords() const
Returns the maximum coordinate across X, Y and Z axis.
QgsPointCloudAttributeCollection attributes() const
Returns the list of attributes contained in the LAZ file.
QByteArray vlrData(QString userId, int recordId)
Returns the binary data of the variable length record with the user identifier userId and record iden...
static QgsLazInfo fromUrl(QUrl &url)
Static function to create a QgsLazInfo class from a file over network.
QVariantMap toMetadata() const
Returns a map containing various metadata extracted from the LAZ file.
QgsVector3D scale() const
Returns the scale of the points coordinates.
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.
QgsVector3D offset() const
Returns the offset of the points coordinates.
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::MessageLevel::Warning, bool notifyUser=true, const char *file=__builtin_FILE(), const char *function=__builtin_FUNCTION(), int line=__builtin_LINE())
Adds a message to the log instance (and creates it if necessary).
A 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.
void finished()
Emitted when the request processing has finished.
Base class for storing raw data from point cloud nodes.
Represents an indexed point cloud node's position in octree.
QString toString() const
Encode node to string.
QgsPointCloudNodeId parentNode() const
Returns the parent of the node.
Keeps metadata for an indexed point cloud node.
QgsBox3D bounds() const
Returns node's bounding cube in CRS coords.
Point cloud data request.
bool ignoreIndexFilterEnabled() const
Returns whether the request will ignore the point cloud index's filter expression,...
QgsPointCloudAttributeCollection attributes() const
Returns attributes.
QgsRectangle filterRect() const
Returns the rectangle from which points will be taken, in point cloud's crs.
Used to store statistics of a point cloud dataset.
static QgsPointCloudStatistics fromStatisticsJson(const 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.
void finished()
Emitted when the reply has finished (either with a success or with a failure)
A 3D vector (similar to QVector3D) with the difference that it uses double precision instead of singl...
double y() const
Returns Y coordinate.
double z() const
Returns Z coordinate.
double x() const
Returns X coordinate.
void set(double x, double y, double z)
Sets vector coordinates.
#define QgsDebugMsgLevel(str, level)
#define QgsDebugError(str)
#define QgsSetRequestInitiatorClass(request, _class)