QGIS API Documentation 3.43.0-Master (ac9f54ad1f7)
qgspointcloudstatscalculator.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgspointcloudstatscalculator.cpp
3 --------------------
4 begin : April 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#include "moc_qgspointcloudstatscalculator.cpp"
20
22
23#include "qgspointcloudindex.h"
24#include "qgsmessagelog.h"
27
29
30#include "qgsfeedback.h"
32
33#include <QQueue>
34#include <QtConcurrent/QtConcurrentMap>
35
37{
40
41 StatsProcessor( QgsPointCloudIndex index, QgsPointCloudRequest request, QgsFeedback *feedback, double progressValue )
42 : mIndex( index ), mRequest( request ), mFeedback( feedback ), mProgressValue( progressValue )
43 {
44 }
45
46 StatsProcessor( const StatsProcessor &processor )
47 : mIndex( processor.mIndex ), mRequest( processor.mRequest ), mFeedback( processor.mFeedback ), mProgressValue( processor.mProgressValue )
48 {
49 }
50
52 {
53 mIndex = rhs.mIndex;
54 mRequest = rhs.mRequest;
55 mFeedback = rhs.mFeedback;
56 mProgressValue = rhs.mProgressValue;
57 return *this;
58 }
59
61 {
62 QgsPointCloudNode node = mIndex.getNode( nodeId );
63 if ( node.pointCount() < 1 )
65
66 std::unique_ptr<QgsPointCloudBlock> block = nullptr;
68 {
69 block = mIndex.nodeData( nodeId, mRequest );
70 }
71 else
72 {
73 QgsPointCloudBlockRequest *request = mIndex.asyncNodeData( nodeId, mRequest );
74 if ( request == nullptr )
75 {
76 QgsDebugError( QStringLiteral( "Unable to calculate statistics for node %1: Got nullptr async request" ).arg( nodeId.toString() ) );
78 }
79 QEventLoop loop;
80 QObject::connect( request, &QgsPointCloudBlockRequest::finished, &loop, &QEventLoop::quit );
81 QObject::connect( mFeedback, &QgsFeedback::canceled, &loop, &QEventLoop::quit );
82 loop.exec();
83 if ( !mFeedback->isCanceled() )
84 {
85 block = request->takeBlock();
86 if ( !block )
87 {
88 QgsMessageLog::logMessage( QObject::tr( "Unable to calculate statistics for node %1, error: \"%2\"" ).arg( nodeId.toString(), request->errorStr() ) );
89 }
90 }
91 }
92
93 if ( !block )
94 {
95 updateFeedback();
97 }
98
99 const QgsPointCloudAttributeCollection attributesCollection = block->attributes();
100 const QVector<QgsPointCloudAttribute> attributes = attributesCollection.attributes();
101 const char *ptr = block->data();
102 int count = block->pointCount();
103 int recordSize = attributesCollection.pointRecordSize();
104
105 QMap<QString, QgsPointCloudAttributeStatistics> statsMap;
106 for ( const QgsPointCloudAttribute &attribute : attributes )
107 {
109 summary.minimum = std::numeric_limits<double>::max();
110 summary.maximum = std::numeric_limits<double>::lowest();
111 summary.count = 0;
112 summary.mean = 0;
113 summary.stDev = std::numeric_limits<double>::quiet_NaN();
114 summary.classCount.clear();
115 statsMap[ attribute.name() ] = std::move( summary );
116 }
117
118 QVector<int> attributeOffsetVector;
119 QSet<int> classifiableAttributesOffsetSet;
120 for ( const QgsPointCloudAttribute &attribute : attributes )
121 {
122 int attributeOffset = 0;
123 attributesCollection.find( attribute.name(), attributeOffset );
124 attributeOffsetVector.push_back( attributeOffset );
125 if ( attribute.name() == QLatin1String( "ScannerChannel" ) ||
126 attribute.name() == QLatin1String( "ReturnNumber" ) ||
127 attribute.name() == QLatin1String( "NumberOfReturns" ) ||
128 attribute.name() == QLatin1String( "ScanDirectionFlag" ) ||
129 attribute.name() == QLatin1String( "Classification" ) ||
130 attribute.name() == QLatin1String( "EdgeOfFlightLine" ) ||
131 attribute.name() == QLatin1String( "PointSourceId" ) ||
132 attribute.name() == QLatin1String( "Synthetic" ) ||
133 attribute.name() == QLatin1String( "KeyPoint" ) ||
134 attribute.name() == QLatin1String( "Withheld" ) ||
135 attribute.name() == QLatin1String( "Overlap" ) )
136 {
137 classifiableAttributesOffsetSet.insert( attributeOffset );
138 }
139 }
140
141 QVector<double> sumOfSquares( attributes.size(), 0 );
142 QVector<double> partialMean( attributes.size(), 0 );
143 for ( int i = 0; i < count; ++i )
144 {
145 for ( int j = 0; j < attributes.size(); ++j )
146 {
147 if ( mFeedback->isCanceled() )
148 {
150 }
151 QString attributeName = attributes.at( j ).name();
152 QgsPointCloudAttribute::DataType attributeType = attributes.at( j ).type();
153
154 double attributeValue = 0;
155 int attributeOffset = attributeOffsetVector[ j ];
156
157 QgsPointCloudAttributeStatistics &stats = statsMap[ attributeName ];
158 QgsPointCloudRenderContext::getAttribute( ptr, i * recordSize + attributeOffset, attributeType, attributeValue );
159 stats.minimum = std::min( stats.minimum, attributeValue );
160 stats.maximum = std::max( stats.maximum, attributeValue );
161 stats.mean += attributeValue / count;
162 stats.count++;
163
164 // Single pass stdev
165 const double delta = attributeValue - partialMean[j];
166 partialMean[j] += delta / stats.count;
167 sumOfSquares[j] += delta * ( attributeValue - partialMean[j] );
168
169 if ( classifiableAttributesOffsetSet.contains( attributeOffset ) )
170 {
171 stats.classCount[( int )attributeValue ]++;
172 }
173 }
174 }
175
176 for ( int j = 0; j < attributes.size(); ++j )
177 {
178 const QString attributeName = attributes.at( j ).name();
179 QgsPointCloudAttributeStatistics &stats = statsMap[ attributeName ];
180 stats.stDev = std::sqrt( sumOfSquares[j] / ( stats.count - 1.0 ) );
181 }
182
183 updateFeedback();
184 return QgsPointCloudStatistics( count, statsMap );
185 }
186 private:
187 QgsPointCloudIndex mIndex;
188 QgsPointCloudRequest mRequest;
189 QgsFeedback *mFeedback = nullptr;
190 double mProgressValue = 0.0;
191
192 void updateFeedback()
193 {
194 QMutexLocker locker( &sStatsProcessorFeedbackMutex );
195 mFeedback->setProgress( mFeedback->progress() + mProgressValue );
196 }
197};
198
200
206
207bool QgsPointCloudStatsCalculator::calculateStats( QgsFeedback *feedback, const QVector<QgsPointCloudAttribute> &attributes, qint64 pointsLimit )
208{
209 if ( !mIndex.isValid() )
210 {
211 QgsMessageLog::logMessage( QObject::tr( "Unable to calculate statistics of an invalid index" ) );
212 return false;
213 }
214 mRequest.setAttributes( attributes );
215
216 qint64 pointCount = 0;
217 QVector<QgsPointCloudNodeId> nodes;
218 QQueue<QgsPointCloudNodeId> queue;
219 queue.push_back( mIndex.root() );
220 while ( !queue.empty() )
221 {
222 QgsPointCloudNode node = mIndex.getNode( queue.front() );
223 queue.pop_front();
224 if ( !mProcessedNodes.contains( node.id() ) )
225 pointCount += node.pointCount();
226 if ( pointsLimit != -1 && pointCount > pointsLimit )
227 break;
228 if ( !mProcessedNodes.contains( node.id() ) )
229 {
230 nodes.push_back( node.id() );
231 mProcessedNodes.insert( node.id() );
232 }
233 for ( const QgsPointCloudNodeId &child : node.children() )
234 {
235 queue.push_back( child );
236 }
237 }
238
239 feedback->setProgress( 0 );
240
241 QVector<QgsPointCloudStatistics> list = QtConcurrent::blockingMapped( nodes, StatsProcessor( mIndex, mRequest, feedback, 100.0 / ( double )nodes.size() ) );
242
243 for ( QgsPointCloudStatistics &s : list )
244 {
245 mStats.combineWith( s );
246 }
247 return !feedback->isCanceled() && mStats.sampledPointsCount() != 0;
248}
@ Local
Local means the source is a local file on the machine.
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition qgsfeedback.h:44
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
void canceled()
Internal routines can connect to this signal if they use event loop.
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:61
double progress() const
Returns the current progress reported by the feedback object.
Definition qgsfeedback.h:77
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.
int pointRecordSize() const
Returns total size of record.
const QgsPointCloudAttribute * find(const QString &attributeName, int &offset) const
Finds the attribute with the name.
QVector< QgsPointCloudAttribute > attributes() const
Returns all attributes.
Attribute for point cloud data pair of name and size in bytes.
DataType
Systems of unit measurement.
Base class for handling loading QgsPointCloudBlock asynchronously.
QString errorStr()
Returns the error message string of the request.
void finished()
Emitted when the request processing has finished.
std::unique_ptr< QgsPointCloudBlock > takeBlock()
Returns the requested block.
Smart pointer for QgsAbstractPointCloudIndex.
QgsPointCloudBlockRequest * asyncNodeData(const QgsPointCloudNodeId &n, const QgsPointCloudRequest &request)
Returns a handle responsible for loading a node data block.
bool isValid() const
Returns whether index is loaded and valid.
std::unique_ptr< QgsPointCloudBlock > nodeData(const QgsPointCloudNodeId &n, const QgsPointCloudRequest &request)
Returns node data block.
QgsPointCloudNodeId root() const
Returns root node of the index.
QgsPointCloudNode getNode(const QgsPointCloudNodeId &id) const
Returns object for a given node.
Qgis::PointCloudAccessType accessType() const
Returns the access type of the data If the access type is Remote, data will be fetched from an HTTP s...
Represents an indexed point cloud node's position in octree.
QString toString() const
Encode node to string.
Keeps metadata for an indexed point cloud node.
QList< QgsPointCloudNodeId > children() const
Returns IDs of child nodes.
qint64 pointCount() const
Returns number of points contained in node data.
QgsPointCloudNodeId id() const
Returns node's ID (unique in index)
static void getAttribute(const char *data, std::size_t offset, QgsPointCloudAttribute::DataType type, T &value)
Retrieves the attribute value from data at the specified offset, where type indicates the original da...
Point cloud data request.
void setAttributes(const QgsPointCloudAttributeCollection &attributes)
Set attributes filter in the request.
Used to store statistics of a point cloud dataset.
void combineWith(const QgsPointCloudStatistics &stats)
Merges the current statistics with the statistics from stats.
int sampledPointsCount() const
Returns the number of points used to calculate the statistics.
QgsPointCloudStatsCalculator(QgsPointCloudIndex index)
Constructor.
bool calculateStats(QgsFeedback *feedback, const QVector< QgsPointCloudAttribute > &attributes, qint64 pointsLimit=-1)
Calculates the statistics of given attributes attributes up to new pointsLimit points Note: the alrea...
#define QgsDebugError(str)
Definition qgslogger.h:40
Stores statistics of one attribute of a point cloud dataset.
static QMutex sStatsProcessorFeedbackMutex
StatsProcessor(const StatsProcessor &processor)
QgsPointCloudStatistics result_type
QgsPointCloudStatistics operator()(QgsPointCloudNodeId nodeId)
StatsProcessor & operator=(const StatsProcessor &rhs)
StatsProcessor(QgsPointCloudIndex index, QgsPointCloudRequest request, QgsFeedback *feedback, double progressValue)