QGIS API Documentation 3.40.0-Bratislava (b56115d8743)
Loading...
Searching...
No Matches
qgseptpointcloudindex.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgspointcloudindex.cpp
3 --------------------
4 begin : October 2020
5 copyright : (C) 2020 by Peter Petrik
6 email : zilolv 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#include <QFile>
20#include <QFileInfo>
21#include <QDir>
22#include <QJsonArray>
23#include <QJsonDocument>
24#include <QJsonObject>
25#include <QTime>
26#include <QtDebug>
27#include <QQueue>
28
29#include "qgseptdecoder.h"
30#include "qgslazdecoder.h"
34#include "qgslogger.h"
35#include "qgspointcloudexpression.h"
36
38
39#define PROVIDER_KEY QStringLiteral( "ept" )
40#define PROVIDER_DESCRIPTION QStringLiteral( "EPT point cloud provider" )
41
42QgsEptPointCloudIndex::QgsEptPointCloudIndex() = default;
43
44QgsEptPointCloudIndex::~QgsEptPointCloudIndex() = default;
45
46std::unique_ptr<QgsPointCloudIndex> QgsEptPointCloudIndex::clone() const
47{
48 QgsEptPointCloudIndex *clone = new QgsEptPointCloudIndex;
49 QMutexLocker locker( &mHierarchyMutex );
50 copyCommonProperties( clone );
51 return std::unique_ptr<QgsPointCloudIndex>( clone );
52}
53
54void QgsEptPointCloudIndex::load( const QString &fileName )
55{
56 mUri = fileName;
57 QFile f( fileName );
58 if ( !f.open( QIODevice::ReadOnly ) )
59 {
60 mError = tr( "Unable to open %1 for reading" ).arg( fileName );
61 mIsValid = false;
62 return;
63 }
64
65 const QDir directory = QFileInfo( fileName ).absoluteDir();
66 mDirectory = directory.absolutePath();
67
68 const QByteArray dataJson = f.readAll();
69 bool success = loadSchema( dataJson );
70
71 if ( success )
72 {
73 // try to import the metadata too!
74 QFile manifestFile( mDirectory + QStringLiteral( "/ept-sources/manifest.json" ) );
75 if ( manifestFile.open( QIODevice::ReadOnly ) )
76 {
77 const QByteArray manifestJson = manifestFile.readAll();
78 loadManifest( manifestJson );
79 }
80 }
81
82 if ( success )
83 {
84 success = loadHierarchy();
85 }
86
87 mIsValid = success;
88}
89
90void QgsEptPointCloudIndex::loadManifest( const QByteArray &manifestJson )
91{
92 QJsonParseError err;
93 // try to import the metadata too!
94 const QJsonDocument manifestDoc = QJsonDocument::fromJson( manifestJson, &err );
95 if ( err.error == QJsonParseError::NoError )
96 {
97 const QJsonArray manifestArray = manifestDoc.array();
98 // TODO how to handle multiple?
99 if ( ! manifestArray.empty() )
100 {
101 const QJsonObject sourceObject = manifestArray.at( 0 ).toObject();
102 const QString metadataPath = sourceObject.value( QStringLiteral( "metadataPath" ) ).toString();
103 QFile metadataFile( mDirectory + QStringLiteral( "/ept-sources/" ) + metadataPath );
104 if ( metadataFile.open( QIODevice::ReadOnly ) )
105 {
106 const QByteArray metadataJson = metadataFile.readAll();
107 const QJsonDocument metadataDoc = QJsonDocument::fromJson( metadataJson, &err );
108 if ( err.error == QJsonParseError::NoError )
109 {
110 const QJsonObject metadataObject = metadataDoc.object().value( QStringLiteral( "metadata" ) ).toObject();
111 if ( !metadataObject.empty() )
112 {
113 const QJsonObject sourceMetadata = metadataObject.constBegin().value().toObject();
114 mOriginalMetadata = sourceMetadata.toVariantMap();
115 }
116 }
117 }
118 }
119 }
120}
121
122bool QgsEptPointCloudIndex::loadSchema( const QByteArray &dataJson )
123{
124 QJsonParseError err;
125 const QJsonDocument doc = QJsonDocument::fromJson( dataJson, &err );
126 if ( err.error != QJsonParseError::NoError )
127 return false;
128 const QJsonObject result = doc.object();
129 mDataType = result.value( QLatin1String( "dataType" ) ).toString(); // "binary" or "laszip"
130 if ( mDataType != QLatin1String( "laszip" ) && mDataType != QLatin1String( "binary" ) && mDataType != QLatin1String( "zstandard" ) )
131 return false;
132
133 const QString hierarchyType = result.value( QLatin1String( "hierarchyType" ) ).toString(); // "json" or "gzip"
134 if ( hierarchyType != QLatin1String( "json" ) )
135 return false;
136
137 mSpan = result.value( QLatin1String( "span" ) ).toInt();
138 mPointCount = result.value( QLatin1String( "points" ) ).toDouble();
139
140 // WKT
141 const QJsonObject srs = result.value( QLatin1String( "srs" ) ).toObject();
142 mWkt = srs.value( QLatin1String( "wkt" ) ).toString();
143
144 // rectangular
145 const QJsonArray bounds = result.value( QLatin1String( "bounds" ) ).toArray();
146 if ( bounds.size() != 6 )
147 return false;
148
149 const QJsonArray boundsConforming = result.value( QLatin1String( "boundsConforming" ) ).toArray();
150 if ( boundsConforming.size() != 6 )
151 return false;
152 mExtent.set( boundsConforming[0].toDouble(), boundsConforming[1].toDouble(),
153 boundsConforming[3].toDouble(), boundsConforming[4].toDouble() );
154 mZMin = boundsConforming[2].toDouble();
155 mZMax = boundsConforming[5].toDouble();
156
157 const QJsonArray schemaArray = result.value( QLatin1String( "schema" ) ).toArray();
159
160 for ( const QJsonValue &schemaItem : schemaArray )
161 {
162 const QJsonObject schemaObj = schemaItem.toObject();
163 const QString name = schemaObj.value( QLatin1String( "name" ) ).toString();
164 const QString type = schemaObj.value( QLatin1String( "type" ) ).toString();
165
166 const int size = schemaObj.value( QLatin1String( "size" ) ).toInt();
167
168 if ( name == QLatin1String( "ClassFlags" ) && size == 1 )
169 {
170 attributes.push_back( QgsPointCloudAttribute( QStringLiteral( "Synthetic" ), QgsPointCloudAttribute::UChar ) );
171 attributes.push_back( QgsPointCloudAttribute( QStringLiteral( "KeyPoint" ), QgsPointCloudAttribute::UChar ) );
172 attributes.push_back( QgsPointCloudAttribute( QStringLiteral( "Withheld" ), QgsPointCloudAttribute::UChar ) );
173 attributes.push_back( QgsPointCloudAttribute( QStringLiteral( "Overlap" ), QgsPointCloudAttribute::UChar ) );
174 }
175 else if ( type == QLatin1String( "float" ) && ( size == 4 ) )
176 {
178 }
179 else if ( type == QLatin1String( "float" ) && ( size == 8 ) )
180 {
182 }
183 else if ( size == 1 )
184 {
186 }
187 else if ( type == QLatin1String( "unsigned" ) && size == 2 )
188 {
190 }
191 else if ( size == 2 )
192 {
194 }
195 else if ( size == 4 )
196 {
198 }
199 else
200 {
201 // unknown attribute type
202 return false;
203 }
204
205 double scale = 1.f;
206 if ( schemaObj.contains( QLatin1String( "scale" ) ) )
207 scale = schemaObj.value( QLatin1String( "scale" ) ).toDouble();
208
209 double offset = 0.f;
210 if ( schemaObj.contains( QLatin1String( "offset" ) ) )
211 offset = schemaObj.value( QLatin1String( "offset" ) ).toDouble();
212
213 if ( name == QLatin1String( "X" ) )
214 {
215 mOffset.set( offset, mOffset.y(), mOffset.z() );
216 mScale.set( scale, mScale.y(), mScale.z() );
217 }
218 else if ( name == QLatin1String( "Y" ) )
219 {
220 mOffset.set( mOffset.x(), offset, mOffset.z() );
221 mScale.set( mScale.x(), scale, mScale.z() );
222 }
223 else if ( name == QLatin1String( "Z" ) )
224 {
225 mOffset.set( mOffset.x(), mOffset.y(), offset );
226 mScale.set( mScale.x(), mScale.y(), scale );
227 }
228
229 // store any metadata stats which are present for the attribute
230 AttributeStatistics stats;
231 bool foundStats = false;
232 if ( schemaObj.contains( QLatin1String( "count" ) ) )
233 {
234 stats.count = schemaObj.value( QLatin1String( "count" ) ).toInt();
235 foundStats = true;
236 }
237 if ( schemaObj.contains( QLatin1String( "minimum" ) ) )
238 {
239 stats.minimum = schemaObj.value( QLatin1String( "minimum" ) ).toDouble();
240 foundStats = true;
241 }
242 if ( schemaObj.contains( QLatin1String( "maximum" ) ) )
243 {
244 stats.maximum = schemaObj.value( QLatin1String( "maximum" ) ).toDouble();
245 foundStats = true;
246 }
247 if ( schemaObj.contains( QLatin1String( "count" ) ) )
248 {
249 stats.mean = schemaObj.value( QLatin1String( "mean" ) ).toDouble();
250 foundStats = true;
251 }
252 if ( schemaObj.contains( QLatin1String( "stddev" ) ) )
253 {
254 stats.stDev = schemaObj.value( QLatin1String( "stddev" ) ).toDouble();
255 foundStats = true;
256 }
257 if ( schemaObj.contains( QLatin1String( "variance" ) ) )
258 {
259 stats.variance = schemaObj.value( QLatin1String( "variance" ) ).toDouble();
260 foundStats = true;
261 }
262 if ( foundStats )
263 mMetadataStats.insert( name, stats );
264
265 if ( schemaObj.contains( QLatin1String( "counts" ) ) )
266 {
267 QMap< int, int > classCounts;
268 const QJsonArray counts = schemaObj.value( QLatin1String( "counts" ) ).toArray();
269 for ( const QJsonValue &count : counts )
270 {
271 const QJsonObject countObj = count.toObject();
272 classCounts.insert( countObj.value( QLatin1String( "value" ) ).toInt(), countObj.value( QLatin1String( "count" ) ).toInt() );
273 }
274 mAttributeClasses.insert( name, classCounts );
275 }
276 }
277 setAttributes( attributes );
278
279 // save mRootBounds
280
281 // bounds (cube - octree volume)
282 const double xmin = bounds[0].toDouble();
283 const double ymin = bounds[1].toDouble();
284 const double zmin = bounds[2].toDouble();
285 const double xmax = bounds[3].toDouble();
286 const double ymax = bounds[4].toDouble();
287 const double zmax = bounds[5].toDouble();
288
289 mRootBounds = QgsPointCloudDataBounds(
290 ( xmin - mOffset.x() ) / mScale.x(),
291 ( ymin - mOffset.y() ) / mScale.y(),
292 ( zmin - mOffset.z() ) / mScale.z(),
293 ( xmax - mOffset.x() ) / mScale.x(),
294 ( ymax - mOffset.y() ) / mScale.y(),
295 ( zmax - mOffset.z() ) / mScale.z()
296 );
297
298
299#ifdef QGIS_DEBUG
300 double dx = xmax - xmin, dy = ymax - ymin, dz = zmax - zmin;
301 QgsDebugMsgLevel( QStringLiteral( "lvl0 node size in CRS units: %1 %2 %3" ).arg( dx ).arg( dy ).arg( dz ), 2 ); // all dims should be the same
302 QgsDebugMsgLevel( QStringLiteral( "res at lvl0 %1" ).arg( dx / mSpan ), 2 );
303 QgsDebugMsgLevel( QStringLiteral( "res at lvl1 %1" ).arg( dx / mSpan / 2 ), 2 );
304 QgsDebugMsgLevel( QStringLiteral( "res at lvl2 %1 with node size %2" ).arg( dx / mSpan / 4 ).arg( dx / 4 ), 2 );
305#endif
306
307 return true;
308}
309
310std::unique_ptr<QgsPointCloudBlock> QgsEptPointCloudIndex::nodeData( const IndexedPointCloudNode &n, const QgsPointCloudRequest &request )
311{
312 if ( QgsPointCloudBlock *cached = getNodeDataFromCache( n, request ) )
313 {
314 return std::unique_ptr<QgsPointCloudBlock>( cached );
315 }
316
317 mHierarchyMutex.lock();
318 const bool found = mHierarchy.contains( n );
319 mHierarchyMutex.unlock();
320 if ( !found )
321 return nullptr;
322
323 // we need to create a copy of the expression to pass to the decoder
324 // as the same QgsPointCloudExpression object mighgt be concurrently
325 // used on another thread, for example in a 3d view
326 QgsPointCloudExpression filterExpression = mFilterExpression;
327 QgsPointCloudAttributeCollection requestAttributes = request.attributes();
328 requestAttributes.extend( attributes(), filterExpression.referencedAttributes() );
329 QgsRectangle filterRect = request.filterRect();
330
331 std::unique_ptr<QgsPointCloudBlock> decoded;
332 if ( mDataType == QLatin1String( "binary" ) )
333 {
334 const QString filename = QStringLiteral( "%1/ept-data/%2.bin" ).arg( mDirectory, n.toString() );
335 decoded = QgsEptDecoder::decompressBinary( filename, attributes(), requestAttributes, scale(), offset(), filterExpression, filterRect );
336 }
337 else if ( mDataType == QLatin1String( "zstandard" ) )
338 {
339 const QString filename = QStringLiteral( "%1/ept-data/%2.zst" ).arg( mDirectory, n.toString() );
340 decoded = QgsEptDecoder::decompressZStandard( filename, attributes(), request.attributes(), scale(), offset(), filterExpression, filterRect );
341 }
342 else if ( mDataType == QLatin1String( "laszip" ) )
343 {
344 const QString filename = QStringLiteral( "%1/ept-data/%2.laz" ).arg( mDirectory, n.toString() );
345 decoded = QgsLazDecoder::decompressLaz( filename, requestAttributes, filterExpression, filterRect );
346 }
347
348 storeNodeDataToCache( decoded.get(), n, request );
349 return decoded;
350}
351
352QgsPointCloudBlockRequest *QgsEptPointCloudIndex::asyncNodeData( const IndexedPointCloudNode &n, const QgsPointCloudRequest &request )
353{
354 Q_UNUSED( n );
355 Q_UNUSED( request );
356 Q_ASSERT( false );
357 return nullptr; // unsupported
358}
359
360QgsCoordinateReferenceSystem QgsEptPointCloudIndex::crs() const
361{
363}
364
365qint64 QgsEptPointCloudIndex::pointCount() const
366{
367 return mPointCount;
368}
369
370bool QgsEptPointCloudIndex::hasStatisticsMetadata() const
371{
372 return !mMetadataStats.isEmpty();
373}
374
375QVariant QgsEptPointCloudIndex::metadataStatistic( const QString &attribute, Qgis::Statistic statistic ) const
376{
377 if ( !mMetadataStats.contains( attribute ) )
378 return QVariant();
379
380 const AttributeStatistics &stats = mMetadataStats[ attribute ];
381 switch ( statistic )
382 {
384 return stats.count >= 0 ? QVariant( stats.count ) : QVariant();
385
387 return std::isnan( stats.mean ) ? QVariant() : QVariant( stats.mean );
388
390 return std::isnan( stats.stDev ) ? QVariant() : QVariant( stats.stDev );
391
393 return stats.minimum;
394
396 return stats.maximum;
397
399 return stats.minimum.isValid() && stats.maximum.isValid() ? QVariant( stats.maximum.toDouble() - stats.minimum.toDouble() ) : QVariant();
400
414 return QVariant();
415 }
416 return QVariant();
417}
418
419QVariantList QgsEptPointCloudIndex::metadataClasses( const QString &attribute ) const
420{
421 QVariantList classes;
422 const QMap< int, int > values = mAttributeClasses.value( attribute );
423 for ( auto it = values.constBegin(); it != values.constEnd(); ++it )
424 {
425 classes << it.key();
426 }
427 return classes;
428}
429
430QVariant QgsEptPointCloudIndex::metadataClassStatistic( const QString &attribute, const QVariant &value, Qgis::Statistic statistic ) const
431{
432 if ( statistic != Qgis::Statistic::Count )
433 return QVariant();
434
435 const QMap< int, int > values = mAttributeClasses.value( attribute );
436 if ( !values.contains( value.toInt() ) )
437 return QVariant();
438 return values.value( value.toInt() );
439}
440
441bool QgsEptPointCloudIndex::loadHierarchy()
442{
443 QQueue<QString> queue;
444 queue.enqueue( QStringLiteral( "0-0-0-0" ) );
445 while ( !queue.isEmpty() )
446 {
447 const QString filename = QStringLiteral( "%1/ept-hierarchy/%2.json" ).arg( mDirectory, queue.dequeue() );
448 QFile fH( filename );
449 if ( !fH.open( QIODevice::ReadOnly ) )
450 {
451 QgsDebugMsgLevel( QStringLiteral( "unable to read hierarchy from file %1" ).arg( filename ), 2 );
452 mError = QStringLiteral( "unable to read hierarchy from file %1" ).arg( filename );
453 return false;
454 }
455
456 const QByteArray dataJsonH = fH.readAll();
457 QJsonParseError errH;
458 const QJsonDocument docH = QJsonDocument::fromJson( dataJsonH, &errH );
459 if ( errH.error != QJsonParseError::NoError )
460 {
461 QgsDebugMsgLevel( QStringLiteral( "QJsonParseError when reading hierarchy from file %1" ).arg( filename ), 2 );
462 mError = QStringLiteral( "QJsonParseError when reading hierarchy from file %1" ).arg( filename );
463 return false;
464 }
465
466 const QJsonObject rootHObj = docH.object();
467 for ( auto it = rootHObj.constBegin(); it != rootHObj.constEnd(); ++it )
468 {
469 const QString nodeIdStr = it.key();
470 const int nodePointCount = it.value().toInt();
471 if ( nodePointCount < 0 )
472 {
473 queue.enqueue( nodeIdStr );
474 }
475 else
476 {
478 mHierarchyMutex.lock();
479 mHierarchy[nodeId] = nodePointCount;
480 mHierarchyMutex.unlock();
481 }
482 }
483 }
484 return true;
485}
486
487bool QgsEptPointCloudIndex::isValid() const
488{
489 return mIsValid;
490}
491
492void QgsEptPointCloudIndex::copyCommonProperties( QgsEptPointCloudIndex *destination ) const
493{
495
496 // QgsEptPointCloudIndex specific fields
497 destination->mIsValid = mIsValid;
498 destination->mDataType = mDataType;
499 destination->mDirectory = mDirectory;
500 destination->mWkt = mWkt;
501 destination->mPointCount = mPointCount;
502 destination->mMetadataStats = mMetadataStats;
503 destination->mAttributeClasses = mAttributeClasses;
504 destination->mOriginalMetadata = mOriginalMetadata;
505}
506
Represents a indexed point cloud node in octree.
static IndexedPointCloudNode fromString(const QString &str)
Creates node from string.
QString toString() const
Encode node to string.
Statistic
Available generic statistics.
Definition qgis.h:5432
@ StDev
Standard deviation of values.
@ FirstQuartile
First quartile.
@ Mean
Mean of values.
@ Median
Median of values.
@ Max
Max of values.
@ Min
Min of values.
@ First
First value.
@ Range
Range of values (max - min)
@ Sum
Sum of values.
@ Minority
Minority of values.
@ CountMissing
Number of missing (null) values.
@ All
All statistics.
@ Majority
Majority of values.
@ Variety
Variety (count of distinct) values.
@ StDevSample
Sample standard deviation of values.
@ Last
Last value.
@ ThirdQuartile
Third quartile.
@ InterQuartileRange
Inter quartile range (IQR)
This class represents a coordinate reference system (CRS).
static QgsCoordinateReferenceSystem fromWkt(const QString &wkt)
Creates a CRS from a WKT spatial ref sys definition string.
Collection of point cloud attributes.
void push_back(const QgsPointCloudAttribute &attribute)
Adds extra attribute.
void extend(const QgsPointCloudAttributeCollection &otherCollection, const QSet< QString > &matchingNames)
Adds specific missing attributes from another QgsPointCloudAttributeCollection.
Attribute for point cloud data pair of name and size in bytes.
@ UShort
Unsigned short int 2 bytes.
@ UChar
Unsigned char 1 byte.
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.
A rectangle specified with double values.
#define QgsDebugMsgLevel(str, level)
Definition qgslogger.h:39