24#include <QMutexLocker>
25#include <QJsonDocument>
36#include "qgspointcloudexpression.h"
38#include "lazperf/lazperf.hpp"
39#include "lazperf/readers.hpp"
40#include "lazperf/vlr.hpp"
44#define PROVIDER_KEY QStringLiteral( "copc" )
45#define PROVIDER_DESCRIPTION QStringLiteral( "COPC point cloud provider" )
47QgsCopcPointCloudIndex::QgsCopcPointCloudIndex() =
default;
49QgsCopcPointCloudIndex::~QgsCopcPointCloudIndex() =
default;
51std::unique_ptr<QgsPointCloudIndex> QgsCopcPointCloudIndex::clone()
const
53 QgsCopcPointCloudIndex *clone =
new QgsCopcPointCloudIndex;
54 QMutexLocker locker( &mHierarchyMutex );
55 copyCommonProperties( clone );
56 return std::unique_ptr<QgsPointCloudIndex>( clone );
59void QgsCopcPointCloudIndex::load(
const QString &fileName )
62 mCopcFile.open( QgsLazDecoder::toNativePath( fileName ), std::ios::binary );
64 if ( !mCopcFile.is_open() || !mCopcFile.good() )
66 mError = tr(
"Unable to open %1 for reading" ).arg( fileName );
72 mIsValid = mLazInfo->isValid();
75 mIsValid = loadSchema( *mLazInfo.get() );
79 mError = tr(
"Unable to recognize %1 as a LAZ file: \"%2\"" ).arg( fileName, mLazInfo->error() );
86bool QgsCopcPointCloudIndex::loadSchema(
QgsLazInfo &lazInfo )
88 QByteArray copcInfoVlrData = lazInfo.
vlrData( QStringLiteral(
"copc" ), 1 );
89 if ( copcInfoVlrData.isEmpty() )
91 mError = tr(
"Invalid COPC file" );
94 mCopcInfoVlr.fill( copcInfoVlrData.data(), copcInfoVlrData.size() );
96 mScale = lazInfo.
scale();
97 mOffset = lazInfo.
offset();
103 mExtent.
set( minCoords.
x(), minCoords.
y(), maxCoords.
x(), maxCoords.
y() );
104 mZMin = minCoords.
z();
105 mZMax = maxCoords.
z();
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;
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()
125 double calculatedSpan = nodeMapExtent( root() ).width() / mCopcInfoVlr.spacing;
126 mSpan = calculatedSpan;
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 );
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 );
143 return std::unique_ptr<QgsPointCloudBlock>( cached );
146 const bool found = fetchNodeHierarchy( n );
149 mHierarchyMutex.lock();
150 int pointCount = mHierarchy.value( n );
151 auto [blockOffset, blockSize] = mHierarchyNodePos.value( n );
152 mHierarchyMutex.unlock();
157 QgsPointCloudExpression filterExpression = mFilterExpression;
159 requestAttributes.
extend( attributes(), filterExpression.referencedAttributes() );
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 );
167 QgsDebugError( QStringLiteral(
"Could not read file %1" ).arg( mUri ) );
172 std::unique_ptr<QgsPointCloudBlock> decoded = QgsLazDecoder::decompressCopc( rawBlockData, *mLazInfo.get(), pointCount, requestAttributes, filterExpression, filterRect );
173 storeNodeDataToCache( decoded.get(), n, request );
187 return mLazInfo->crs();
190qint64 QgsCopcPointCloudIndex::pointCount()
const
192 return mLazInfo->pointCount();
195bool QgsCopcPointCloudIndex::loadHierarchy()
197 fetchHierarchyPage( mCopcInfoVlr.root_hier_offset, mCopcInfoVlr.root_hier_size );
203 if ( mLazInfo->version() != qMakePair<uint8_t, uint8_t>( 1, 4 ) )
210 QByteArray statisticsEvlrData = fetchCopcStatisticsEvlrData();
211 if ( !statisticsEvlrData.isEmpty() )
213 QgsMessageLog::logMessage( tr(
"Can't write statistics to \"%1\": file already contains COPC statistics!" ).arg( mUri ) );
217 lazperf::evlr_header statsEvlrHeader;
218 statsEvlrHeader.user_id =
"qgis";
219 statsEvlrHeader.record_id = 0;
220 statsEvlrHeader.description =
"Contains calculated statistics";
222 statsEvlrHeader.data_length = statsJson.size();
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() )
231 lazperf::header14 header = mLazInfo->header();
232 header.evlr_count = header.evlr_count + 1;
234 header.write( copcFile );
237 copcFile.seekg( 0, std::ios::end );
239 statsEvlrHeader.write( copcFile );
240 copcFile.write( statsJson.data(), statsEvlrHeader.data_length );
248 mCopcFile.open( QgsLazDecoder::toNativePath( mUri ), std::ios::binary );
254 QByteArray statisticsEvlrData = fetchCopcStatisticsEvlrData();
256 if ( statisticsEvlrData.isEmpty() )
264bool QgsCopcPointCloudIndex::isValid()
const
271 QMutexLocker locker( &mHierarchyMutex );
273 QVector<IndexedPointCloudNode> ancestors;
275 while ( !mHierarchy.contains( foundRoot ) )
277 ancestors.push_front( foundRoot );
280 ancestors.push_front( foundRoot );
283 auto hierarchyIt = mHierarchy.constFind( n );
284 if ( hierarchyIt == mHierarchy.constEnd() )
286 int nodesCount = *hierarchyIt;
287 if ( nodesCount < 0 )
289 auto hierarchyNodePos = mHierarchyNodePos.constFind( n );
290 mHierarchyMutex.unlock();
291 fetchHierarchyPage( hierarchyNodePos->first, hierarchyNodePos->second );
292 mHierarchyMutex.lock();
295 return mHierarchy.contains( n );
298void QgsCopcPointCloudIndex::fetchHierarchyPage( uint64_t offset, uint64_t byteSize )
const
300 mCopcFile.seekg( offset );
301 std::unique_ptr<char []> data(
new char[ byteSize ] );
302 mCopcFile.read( data.get(), byteSize );
304 populateHierarchy( data.get(), byteSize );
307void QgsCopcPointCloudIndex::populateHierarchy(
const char *hierarchyPageData, uint64_t byteSize )
const
325 QMutexLocker locker( &mHierarchyMutex );
327 for ( uint64_t i = 0; i < byteSize; i +=
sizeof( CopcEntry ) )
329 const CopcEntry *entry =
reinterpret_cast<const CopcEntry *
>( hierarchyPageData + i );
331 mHierarchy[nodeId] = entry->pointCount;
332 mHierarchyNodePos.insert( nodeId, QPair<uint64_t, int32_t>( entry->offset, entry->byteSize ) );
338 return fetchNodeHierarchy( n );
341QList<IndexedPointCloudNode> QgsCopcPointCloudIndex::nodeChildren(
const IndexedPointCloudNode &n )
const
343 fetchNodeHierarchy( n );
345 mHierarchyMutex.lock();
347 auto hierarchyIt = mHierarchy.constFind( n );
348 Q_ASSERT( hierarchyIt != mHierarchy.constEnd() );
349 QList<IndexedPointCloudNode> lst;
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();
357 for (
int i = 0; i < 8; ++i )
359 int dx = i & 1, dy = !!( i & 2 ), dz = !!( i & 4 );
361 if ( fetchNodeHierarchy( n2 ) && mHierarchy[n] >= 0 )
367void QgsCopcPointCloudIndex::copyCommonProperties( QgsCopcPointCloudIndex *destination )
const
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 ) );
381QByteArray QgsCopcPointCloudIndex::fetchCopcStatisticsEvlrData()
383 uint64_t offset = mLazInfo->firstEvlrOffset();
384 uint32_t evlrCount = mLazInfo->evlrCount();
386 QByteArray statisticsEvlrData;
388 for ( uint32_t i = 0; i < evlrCount; ++i )
390 lazperf::evlr_header header;
391 mCopcFile.seekg( offset );
393 mCopcFile.read( buffer, 60 );
394 header.fill( buffer, 60 );
397 if ( header.user_id ==
"qgis" && header.record_id == 0 )
399 statisticsEvlrData = QByteArray( header.data_length, Qt::Initialization::Uninitialized );
400 mCopcFile.read( statisticsEvlrData.data(), header.data_length );
404 offset += 60 + header.data_length;
407 return statisticsEvlrData;
Represents a indexed point cloud node in octree.
IndexedPointCloudNode parentNode() const
Returns the parent of the node.
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...
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...
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)
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.
Class for storage of 3D vectors similar to QVector3D, with the difference that it uses double precisi...
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)