1 /***************************************************************************
2 qgspointclouddataprovider.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
18 #include "qgis.h"
19 #include "qgspointclouddataprovider.h"
20 #include "qgspointcloudindex.h"
21 #include "qgsgeometry.h"
22 #include "qgspointcloudrequest.h"
23 #include "qgsgeometryengine.h"
24 #include <mutex>
25 #include <QDebug>
26 #include <QtMath>
27
28 #include <QtConcurrent/QtConcurrentMap>
29
QgsPointCloudDataProvider(const QString & uri,const QgsDataProvider::ProviderOptions & options,QgsDataProvider::ReadFlags flags)30 QgsPointCloudDataProvider::QgsPointCloudDataProvider(
31 const QString &uri,
32 const QgsDataProvider::ProviderOptions &options,
33 QgsDataProvider::ReadFlags flags )
34 : QgsDataProvider( uri, options, flags )
35 {
36 }
37
38 QgsPointCloudDataProvider::~QgsPointCloudDataProvider() = default;
39
capabilities() const40 QgsPointCloudDataProvider::Capabilities QgsPointCloudDataProvider::capabilities() const
41 {
42 return QgsPointCloudDataProvider::NoCapabilities;
43 }
44
hasValidIndex() const45 bool QgsPointCloudDataProvider::hasValidIndex() const
46 {
47 QgsPointCloudIndex *lIndex = index();
48 return lIndex && lIndex->isValid();
49 }
50
polygonBounds() const51 QgsGeometry QgsPointCloudDataProvider::polygonBounds() const
52 {
53 return QgsGeometry::fromRect( extent() );
54 }
55
originalMetadata() const56 QVariantMap QgsPointCloudDataProvider::originalMetadata() const
57 {
58 return QVariantMap();
59 }
60
createRenderer(const QVariantMap &) const61 QgsPointCloudRenderer *QgsPointCloudDataProvider::createRenderer( const QVariantMap & ) const
62 {
63 return nullptr;
64 }
65
lasClassificationCodes()66 QMap<int, QString> QgsPointCloudDataProvider::lasClassificationCodes()
67 {
68 static QMap< int, QString > sCodes
69 {
70 {0, QStringLiteral( "Created, Never Classified" )},
71 {1, QStringLiteral( "Unclassified" )},
72 {2, QStringLiteral( "Ground" )},
73 {3, QStringLiteral( "Low Vegetation" )},
74 {4, QStringLiteral( "Medium Vegetation" )},
75 {5, QStringLiteral( "High Vegetation" )},
76 {6, QStringLiteral( "Building" )},
77 {7, QStringLiteral( "Low Point (Low Noise)" )},
78 {8, QStringLiteral( "Reserved" )},
79 {9, QStringLiteral( "Water" )},
80 {10, QStringLiteral( "Rail" )},
81 {11, QStringLiteral( "Road Surface" )},
82 {12, QStringLiteral( "Reserved" )},
83 {13, QStringLiteral( "Wire - Guard (Shield)" )},
84 {14, QStringLiteral( "Wire - Conductor (Phase)" )},
85 {15, QStringLiteral( "Transmission Tower" )},
86 {16, QStringLiteral( "Wire-Structure Connector (Insulator)" )},
87 {17, QStringLiteral( "Bridge Deck" )},
88 {18, QStringLiteral( "High Noise" )},
89 };
90
91 static std::once_flag initialized;
92 std::call_once( initialized, [ = ]( )
93 {
94 for ( int i = 19; i <= 63; ++i )
95 sCodes.insert( i, QStringLiteral( "Reserved" ) );
96 for ( int i = 64; i <= 255; ++i )
97 sCodes.insert( i, QStringLiteral( "User Definable" ) );
98 } );
99
100 return sCodes;
101 }
102
translatedLasClassificationCodes()103 QMap<int, QString> QgsPointCloudDataProvider::translatedLasClassificationCodes()
104 {
105 static QMap< int, QString > sCodes
106 {
107 {0, QObject::tr( "Created, Never Classified" )},
108 {1, QObject::tr( "Unclassified" )},
109 {2, QObject::tr( "Ground" )},
110 {3, QObject::tr( "Low Vegetation" )},
111 {4, QObject::tr( "Medium Vegetation" )},
112 {5, QObject::tr( "High Vegetation" )},
113 {6, QObject::tr( "Building" )},
114 {7, QObject::tr( "Low Point (Noise)" )},
115 {8, QObject::tr( "Reserved" )},
116 {9, QObject::tr( "Water" )},
117 {10, QObject::tr( "Rail" )},
118 {11, QObject::tr( "Road Surface" )},
119 {12, QObject::tr( "Reserved" )},
120 {13, QObject::tr( "Wire - Guard (Shield)" )},
121 {14, QObject::tr( "Wire - Conductor (Phase)" )},
122 {15, QObject::tr( "Transmission Tower" )},
123 {16, QObject::tr( "Wire-Structure Connector (Insulator)" )},
124 {17, QObject::tr( "Bridge Deck" )},
125 {18, QObject::tr( "High Noise" )},
126 };
127
128 static std::once_flag initialized;
129 std::call_once( initialized, [ = ]( )
130 {
131 for ( int i = 19; i <= 63; ++i )
132 sCodes.insert( i, QObject::tr( "Reserved" ) );
133 for ( int i = 64; i <= 255; ++i )
134 sCodes.insert( i, QObject::tr( "User Definable" ) );
135 } );
136
137 return sCodes;
138 }
139
dataFormatIds()140 QMap<int, QString> QgsPointCloudDataProvider::dataFormatIds()
141 {
142 static const QMap< int, QString > sCodes
143 {
144 {0, QStringLiteral( "No color or time stored" )},
145 {1, QStringLiteral( "Time is stored" )},
146 {2, QStringLiteral( "Color is stored" )},
147 {3, QStringLiteral( "Color and time are stored" )},
148 {6, QStringLiteral( "Time is stored" )},
149 {7, QStringLiteral( "Time and color are stored)" )},
150 {8, QStringLiteral( "Time, color and near infrared are stored" )},
151 };
152
153 return sCodes;
154 }
155
translatedDataFormatIds()156 QMap<int, QString> QgsPointCloudDataProvider::translatedDataFormatIds()
157 {
158 static const QMap< int, QString > sCodes
159 {
160 {0, QObject::tr( "No color or time stored" )},
161 {1, QObject::tr( "Time is stored" )},
162 {2, QObject::tr( "Color is stored" )},
163 {3, QObject::tr( "Color and time are stored" )},
164 {6, QObject::tr( "Time is stored" )},
165 {7, QObject::tr( "Time and color are stored)" )},
166 {8, QObject::tr( "Time, color and near infrared are stored" )},
167 };
168
169 return sCodes;
170 }
171
metadataStatistic(const QString &,QgsStatisticalSummary::Statistic) const172 QVariant QgsPointCloudDataProvider::metadataStatistic( const QString &, QgsStatisticalSummary::Statistic ) const
173 {
174 return QVariant();
175 }
176
metadataClasses(const QString &) const177 QVariantList QgsPointCloudDataProvider::metadataClasses( const QString & ) const
178 {
179 return QVariantList();
180 }
181
metadataClassStatistic(const QString &,const QVariant &,QgsStatisticalSummary::Statistic) const182 QVariant QgsPointCloudDataProvider::metadataClassStatistic( const QString &, const QVariant &, QgsStatisticalSummary::Statistic ) const
183 {
184 return QVariant();
185 }
186
187 struct MapIndexedPointCloudNode
188 {
189 typedef QVector<QMap<QString, QVariant>> result_type;
190
MapIndexedPointCloudNodeMapIndexedPointCloudNode191 MapIndexedPointCloudNode( QgsPointCloudRequest &request, const QgsVector3D &indexScale, const QgsVector3D &indexOffset,
192 const QgsGeometry &extentGeometry, const QgsDoubleRange &zRange, QgsPointCloudIndex *index, int pointsLimit )
193 : mRequest( request ), mIndexScale( indexScale ), mIndexOffset( indexOffset ), mExtentGeometry( extentGeometry ), mZRange( zRange ), mIndex( index ), mPointsLimit( pointsLimit )
194 { }
195
operator ()MapIndexedPointCloudNode196 QVector<QVariantMap> operator()( IndexedPointCloudNode n )
197 {
198 QVector<QVariantMap> acceptedPoints;
199 std::unique_ptr<QgsPointCloudBlock> block( mIndex->nodeData( n, mRequest ) );
200
201 if ( !block || pointsCount == mPointsLimit )
202 return acceptedPoints;
203
204 const char *ptr = block->data();
205 const QgsPointCloudAttributeCollection blockAttributes = block->attributes();
206 const std::size_t recordSize = blockAttributes.pointRecordSize();
207 int xOffset = 0, yOffset = 0, zOffset = 0;
208 const QgsPointCloudAttribute::DataType xType = blockAttributes.find( QStringLiteral( "X" ), xOffset )->type();
209 const QgsPointCloudAttribute::DataType yType = blockAttributes.find( QStringLiteral( "Y" ), yOffset )->type();
210 const QgsPointCloudAttribute::DataType zType = blockAttributes.find( QStringLiteral( "Z" ), zOffset )->type();
211 std::unique_ptr< QgsGeometryEngine > extentEngine( QgsGeometry::createGeometryEngine( mExtentGeometry.constGet() ) );
212 extentEngine->prepareGeometry();
213 for ( int i = 0; i < block->pointCount() && pointsCount < mPointsLimit; ++i )
214 {
215 double x, y, z;
216 QgsPointCloudAttribute::getPointXYZ( ptr, i, recordSize, xOffset, xType, yOffset, yType, zOffset, zType, block->scale(), block->offset(), x, y, z );
217 QgsPoint point( x, y );
218
219 if ( mZRange.contains( z ) && extentEngine->contains( &point ) )
220 {
221 QVariantMap pointAttr = QgsPointCloudAttribute::getAttributeMap( ptr, i * recordSize, blockAttributes );
222 pointAttr[ QStringLiteral( "X" ) ] = x;
223 pointAttr[ QStringLiteral( "Y" ) ] = y;
224 pointAttr[ QStringLiteral( "Z" ) ] = z;
225 pointsCount++;
226 acceptedPoints.push_back( pointAttr );
227 }
228 }
229 return acceptedPoints;
230 }
231
232 QgsPointCloudRequest &mRequest;
233 QgsVector3D mIndexScale;
234 QgsVector3D mIndexOffset;
235 const QgsGeometry &mExtentGeometry;
236 const QgsDoubleRange &mZRange;
237 QgsPointCloudIndex *mIndex = nullptr;
238 int mPointsLimit;
239 int pointsCount = 0;
240 };
241
identify(double maxError,const QgsGeometry & extentGeometry,const QgsDoubleRange & extentZRange,int pointsLimit)242 QVector<QVariantMap> QgsPointCloudDataProvider::identify(
243 double maxError,
244 const QgsGeometry &extentGeometry,
245 const QgsDoubleRange &extentZRange, int pointsLimit )
246 {
247 QVector<QVariantMap> acceptedPoints;
248
249 QgsPointCloudIndex *index = this->index();
250 const IndexedPointCloudNode root = index->root();
251
252 const QgsRectangle rootNodeExtent = index->nodeMapExtent( root );
253 const double rootError = rootNodeExtent.width() / index->span();
254
255 const QVector<IndexedPointCloudNode> nodes = traverseTree( index, root, maxError, rootError, extentGeometry, extentZRange );
256
257 const QgsPointCloudAttributeCollection attributeCollection = index->attributes();
258 QgsPointCloudRequest request;
259 request.setAttributes( attributeCollection );
260
261 acceptedPoints = QtConcurrent::blockingMappedReduced( nodes,
262 MapIndexedPointCloudNode( request, index->scale(), index->offset(), extentGeometry, extentZRange, index, pointsLimit ),
263 qOverload<const QVector<QMap<QString, QVariant>>&>( &QVector<QMap<QString, QVariant>>::append ),
264 QtConcurrent::UnorderedReduce );
265
266 return acceptedPoints;
267 }
268
traverseTree(const QgsPointCloudIndex * pc,IndexedPointCloudNode n,double maxError,double nodeError,const QgsGeometry & extentGeometry,const QgsDoubleRange & extentZRange)269 QVector<IndexedPointCloudNode> QgsPointCloudDataProvider::traverseTree(
270 const QgsPointCloudIndex *pc,
271 IndexedPointCloudNode n,
272 double maxError,
273 double nodeError,
274 const QgsGeometry &extentGeometry,
275 const QgsDoubleRange &extentZRange )
276 {
277 QVector<IndexedPointCloudNode> nodes;
278
279 const QgsDoubleRange nodeZRange = pc->nodeZRange( n );
280 if ( !extentZRange.overlaps( nodeZRange ) )
281 return nodes;
282
283 if ( !extentGeometry.intersects( pc->nodeMapExtent( n ) ) )
284 return nodes;
285
286 nodes.append( n );
287
288 const double childrenError = nodeError / 2.0;
289 if ( childrenError < maxError )
290 return nodes;
291
292 const QList<IndexedPointCloudNode> children = pc->nodeChildren( n );
293 for ( const IndexedPointCloudNode &nn : children )
294 {
295 if ( extentGeometry.intersects( pc->nodeMapExtent( nn ) ) )
296 nodes += traverseTree( pc, nn, maxError, childrenError, extentGeometry, extentZRange );
297 }
298
299 return nodes;
300 }
301