QGIS API Documentation 3.31.0-Master (7eca95a70e)
qgspointcloudlayerchunkloader_p.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgspointcloudlayerchunkloader_p.cpp
3 --------------------------------------
4 Date : October 2020
5 Copyright : (C) 2020 by Peter Petrik
6 Email : zilolv dot sk at gmail dot com
7 ***************************************************************************
8 * *
9 * This program is free software; you can redistribute it and/or modify *
10 * it under the terms of the GNU General Public License as published by *
11 * the Free Software Foundation; either version 2 of the License, or *
12 * (at your option) any later version. *
13 * *
14 ***************************************************************************/
15
17
18#include "qgs3dutils.h"
20#include "qgschunknode_p.h"
21#include "qgslogger.h"
22#include "qgspointcloudindex.h"
24#include "qgseventtracing.h"
25
26
31
32#include <QtConcurrent>
33#if QT_VERSION < QT_VERSION_CHECK(6, 0, 0)
34#include <Qt3DRender/QAttribute>
35#else
36#include <Qt3DCore/QAttribute>
37#endif
38#include <Qt3DRender/QTechnique>
39#include <Qt3DRender/QShaderProgram>
40#include <Qt3DRender/QGraphicsApiFilter>
41#include <QPointSize>
42
44
45
47
48QgsPointCloudLayerChunkLoader::QgsPointCloudLayerChunkLoader( const QgsPointCloudLayerChunkLoaderFactory *factory, QgsChunkNode *node, std::unique_ptr< QgsPointCloud3DSymbol > symbol,
49 const QgsCoordinateTransform &coordinateTransform, double zValueScale, double zValueOffset )
50 : QgsChunkLoader( node )
51 , mFactory( factory )
52 , mContext( factory->mMap, coordinateTransform, std::move( symbol ), zValueScale, zValueOffset )
53{
54
55 QgsPointCloudIndex *pc = mFactory->mPointCloudIndex;
56 mContext.setAttributes( pc->attributes() );
57
58 const QgsChunkNodeId nodeId = node->tileId();
59 const IndexedPointCloudNode pcNode( nodeId.d, nodeId.x, nodeId.y, nodeId.z );
60
61 Q_ASSERT( pc->hasNode( pcNode ) );
62
63 QgsDebugMsgLevel( QStringLiteral( "loading entity %1" ).arg( node->tileId().text() ), 2 );
64
65 if ( mContext.symbol()->symbolType() == QLatin1String( "single-color" ) )
66 mHandler.reset( new QgsSingleColorPointCloud3DSymbolHandler() );
67 else if ( mContext.symbol()->symbolType() == QLatin1String( "color-ramp" ) )
68 mHandler.reset( new QgsColorRampPointCloud3DSymbolHandler() );
69 else if ( mContext.symbol()->symbolType() == QLatin1String( "rgb" ) )
70 mHandler.reset( new QgsRGBPointCloud3DSymbolHandler() );
71 else if ( mContext.symbol()->symbolType() == QLatin1String( "classification" ) )
72 {
73 mHandler.reset( new QgsClassificationPointCloud3DSymbolHandler() );
74 const QgsClassificationPointCloud3DSymbol *classificationSymbol = dynamic_cast<const QgsClassificationPointCloud3DSymbol *>( mContext.symbol() );
75 mContext.setFilteredOutCategories( classificationSymbol->getFilteredOutCategories() );
76 }
77
78 //
79 // this will be run in a background thread
80 //
81 mFutureWatcher = new QFutureWatcher<void>( this );
82 connect( mFutureWatcher, &QFutureWatcher<void>::finished, this, &QgsChunkQueueJob::finished );
83
84 const QgsAABB bbox = node->bbox();
85 const QFuture<void> future = QtConcurrent::run( [pc, pcNode, bbox, this]
86 {
87 const QgsEventTracing::ScopedEvent e( QStringLiteral( "3D" ), QStringLiteral( "PC chunk load" ) );
88
89 if ( mContext.isCanceled() )
90 {
91 QgsDebugMsgLevel( QStringLiteral( "canceled" ), 2 );
92 return;
93 }
94
95 mHandler->processNode( pc, pcNode, mContext );
96 if ( mContext.symbol()->renderAsTriangles() )
97 mHandler->triangulate( pc, pcNode, mContext, bbox );
98 } );
99
100 // emit finished() as soon as the handler is populated with features
101 mFutureWatcher->setFuture( future );
102}
103
104QgsPointCloudLayerChunkLoader::~QgsPointCloudLayerChunkLoader()
105{
106 if ( mFutureWatcher && !mFutureWatcher->isFinished() )
107 {
108 disconnect( mFutureWatcher, &QFutureWatcher<void>::finished, this, &QgsChunkQueueJob::finished );
109 mContext.cancelRendering();
110 mFutureWatcher->waitForFinished();
111 }
112}
113
114void QgsPointCloudLayerChunkLoader::cancel()
115{
116 mContext.cancelRendering();
117}
118
119Qt3DCore::QEntity *QgsPointCloudLayerChunkLoader::createEntity( Qt3DCore::QEntity *parent )
120{
121 QgsPointCloudIndex *pc = mFactory->mPointCloudIndex;
122 const QgsChunkNodeId nodeId = mNode->tileId();
123 const IndexedPointCloudNode pcNode( nodeId.d, nodeId.x, nodeId.y, nodeId.z );
124 Q_ASSERT( pc->hasNode( pcNode ) );
125
126 Qt3DCore::QEntity *entity = new Qt3DCore::QEntity( parent );
127 mHandler->finalize( entity, mContext );
128 return entity;
129}
130
132
133
134QgsPointCloudLayerChunkLoaderFactory::QgsPointCloudLayerChunkLoaderFactory( const Qgs3DMapSettings &map, const QgsCoordinateTransform &coordinateTransform, QgsPointCloudIndex *pc, QgsPointCloud3DSymbol *symbol,
135 double zValueScale, double zValueOffset, int pointBudget )
136 : mMap( map )
137 , mCoordinateTransform( coordinateTransform )
138 , mPointCloudIndex( pc )
139 , mZValueScale( zValueScale )
140 , mZValueOffset( zValueOffset )
141 , mPointBudget( pointBudget )
142{
143 mSymbol.reset( symbol );
144
145 try
146 {
147 mExtent = mCoordinateTransform.transformBoundingBox( mMap.extent(), Qgis::TransformDirection::Reverse );
148 }
149 catch ( const QgsCsException & )
150 {
151 // bad luck, can't reproject for some reason
152 QgsDebugError( QStringLiteral( "Transformation of extent failed." ) );
153 }
154}
155
156QgsChunkLoader *QgsPointCloudLayerChunkLoaderFactory::createChunkLoader( QgsChunkNode *node ) const
157{
158 const QgsChunkNodeId id = node->tileId();
159
160 Q_ASSERT( mPointCloudIndex->hasNode( IndexedPointCloudNode( id.d, id.x, id.y, id.z ) ) );
161 QgsPointCloud3DSymbol *symbol = static_cast< QgsPointCloud3DSymbol * >( mSymbol->clone() );
162 return new QgsPointCloudLayerChunkLoader( this, node, std::unique_ptr< QgsPointCloud3DSymbol >( symbol ), mCoordinateTransform, mZValueScale, mZValueOffset );
163}
164
165int QgsPointCloudLayerChunkLoaderFactory::primitivesCount( QgsChunkNode *node ) const
166{
167 const QgsChunkNodeId id = node->tileId();
168 const IndexedPointCloudNode n( id.d, id.x, id.y, id.z );
169 Q_ASSERT( mPointCloudIndex->hasNode( n ) );
170 return mPointCloudIndex->nodePointCount( n );
171}
172
173QgsAABB nodeBoundsToAABB( QgsPointCloudDataBounds nodeBounds, QgsVector3D offset, QgsVector3D scale, const Qgs3DMapSettings &map, const QgsCoordinateTransform &coordinateTransform, double zValueOffset );
174
175QgsChunkNode *QgsPointCloudLayerChunkLoaderFactory::createRootNode() const
176{
177 const QgsAABB bbox = nodeBoundsToAABB( mPointCloudIndex->nodeBounds( IndexedPointCloudNode( 0, 0, 0, 0 ) ), mPointCloudIndex->offset(), mPointCloudIndex->scale(), mMap, mCoordinateTransform, mZValueOffset );
178 const float error = mPointCloudIndex->nodeError( IndexedPointCloudNode( 0, 0, 0, 0 ) );
179 return new QgsChunkNode( QgsChunkNodeId( 0, 0, 0, 0 ), bbox, error );
180}
181
182QVector<QgsChunkNode *> QgsPointCloudLayerChunkLoaderFactory::createChildren( QgsChunkNode *node ) const
183{
184 QVector<QgsChunkNode *> children;
185 const QgsChunkNodeId nodeId = node->tileId();
186 const QgsAABB bbox = node->bbox();
187 const float childError = node->error() / 2;
188 float xc = bbox.xCenter(), yc = bbox.yCenter(), zc = bbox.zCenter();
189
190 for ( int i = 0; i < 8; ++i )
191 {
192 int dx = i & 1, dy = !!( i & 2 ), dz = !!( i & 4 );
193 const QgsChunkNodeId childId( nodeId.d + 1, nodeId.x * 2 + dx, nodeId.y * 2 + dy, nodeId.z * 2 + dz );
194
195 if ( !mPointCloudIndex->hasNode( IndexedPointCloudNode( childId.d, childId.x, childId.y, childId.z ) ) )
196 continue;
197 if ( !mExtent.isEmpty() &&
198 !mPointCloudIndex->nodeMapExtent( IndexedPointCloudNode( childId.d, childId.x, childId.y, childId.z ) ).intersects( mExtent ) )
199 continue;
200
201 // the Y and Z coordinates below are intentionally flipped, because
202 // in chunk node IDs the X,Y axes define horizontal plane,
203 // while in our 3D scene the X,Z axes define the horizontal plane
204 const float chXMin = dx ? xc : bbox.xMin;
205 const float chXMax = dx ? bbox.xMax : xc;
206 // Z axis: values are increasing to the south
207 const float chZMin = !dy ? zc : bbox.zMin;
208 const float chZMax = !dy ? bbox.zMax : zc;
209 const float chYMin = dz ? yc : bbox.yMin;
210 const float chYMax = dz ? bbox.yMax : yc;
211 children << new QgsChunkNode( childId, QgsAABB( chXMin, chYMin, chZMin, chXMax, chYMax, chZMax ), childError, node );
212 }
213 return children;
214}
215
217
218
219QgsAABB nodeBoundsToAABB( QgsPointCloudDataBounds nodeBounds, QgsVector3D offset, QgsVector3D scale, const Qgs3DMapSettings &map, const QgsCoordinateTransform &coordinateTransform, double zValueOffset )
220{
221 QgsVector3D extentMin3D( nodeBounds.xMin() * scale.x() + offset.x(), nodeBounds.yMin() * scale.y() + offset.y(), nodeBounds.zMin() * scale.z() + offset.z() + zValueOffset );
222 QgsVector3D extentMax3D( nodeBounds.xMax() * scale.x() + offset.x(), nodeBounds.yMax() * scale.y() + offset.y(), nodeBounds.zMax() * scale.z() + offset.z() + zValueOffset );
223 QgsCoordinateTransform extentTransform = coordinateTransform;
224 extentTransform.setBallparkTransformsAreAppropriate( true );
225 try
226 {
227 extentMin3D = extentTransform.transform( extentMin3D );
228 extentMax3D = extentTransform.transform( extentMax3D );
229 }
230 catch ( QgsCsException & )
231 {
232 QgsDebugError( QStringLiteral( "Error transforming node bounds coordinate" ) );
233 }
234 const QgsVector3D worldExtentMin3D = Qgs3DUtils::mapToWorldCoordinates( extentMin3D, map.origin() );
235 const QgsVector3D worldExtentMax3D = Qgs3DUtils::mapToWorldCoordinates( extentMax3D, map.origin() );
236 QgsAABB rootBbox( worldExtentMin3D.x(), worldExtentMin3D.y(), worldExtentMin3D.z(),
237 worldExtentMax3D.x(), worldExtentMax3D.y(), worldExtentMax3D.z() );
238 return rootBbox;
239}
240
241
242QgsPointCloudLayerChunkedEntity::QgsPointCloudLayerChunkedEntity( QgsPointCloudIndex *pc, const Qgs3DMapSettings &map,
243 const QgsCoordinateTransform &coordinateTransform, QgsPointCloud3DSymbol *symbol,
244 float maximumScreenSpaceError, bool showBoundingBoxes,
245 double zValueScale, double zValueOffset,
246 int pointBudget )
247 : QgsChunkedEntity( maximumScreenSpaceError,
248 new QgsPointCloudLayerChunkLoaderFactory( map, coordinateTransform, pc, symbol, zValueScale, zValueOffset, pointBudget ), true, pointBudget )
249{
250 setUsingAdditiveStrategy( !symbol->renderAsTriangles() );
251 setShowBoundingBoxes( showBoundingBoxes );
252}
253
254QgsPointCloudLayerChunkedEntity::~QgsPointCloudLayerChunkedEntity()
255{
256 // cancel / wait for jobs
257 cancelActiveJobs();
258}
259
260QVector<QgsRayCastingUtils::RayHit> QgsPointCloudLayerChunkedEntity::rayIntersection( const QgsRayCastingUtils::Ray3D &ray, const QgsRayCastingUtils::RayCastContext &context ) const
261{
262 QVector<QgsRayCastingUtils::RayHit> result;
263 QgsPointCloudLayerChunkLoaderFactory *factory = static_cast<QgsPointCloudLayerChunkLoaderFactory *>( mChunkLoaderFactory );
264
265 // transform ray
266 const QgsVector3D rayOriginMapCoords = factory->mMap.worldToMapCoordinates( ray.origin() );
267 const QgsVector3D pointMapCoords = factory->mMap.worldToMapCoordinates( ray.origin() + ray.origin().length() * ray.direction().normalized() );
268 QgsVector3D rayDirectionMapCoords = pointMapCoords - rayOriginMapCoords;
269 rayDirectionMapCoords.normalize();
270
271 const int screenSizePx = std::max( context.screenSize.height(), context.screenSize.width() );
272
273 const QgsPointCloud3DSymbol *symbol = factory->mSymbol.get();
274 // Symbol can be null in case of no rendering enabled
275 if ( !symbol )
276 return result;
277 const double pointSize = symbol->pointSize();
278
279 // We're using the angle as a tolerance, effectively meaning we're fetching points intersecting a cone.
280 // This may be revisited to use a cylinder instead, if the balance between near/far points does not scale
281 // well with different point sizes, screen sizes and fov values.
282 const double limitAngle = 2. * pointSize / screenSizePx * factory->mMap.fieldOfView();
283
284 // adjust ray to elevation properties
285 const QgsVector3D adjustedRayOrigin = QgsVector3D( rayOriginMapCoords.x(), rayOriginMapCoords.y(), ( rayOriginMapCoords.z() - factory->mZValueOffset ) / factory->mZValueScale );
286 QgsVector3D adjustedRayDirection = QgsVector3D( rayDirectionMapCoords.x(), rayDirectionMapCoords.y(), rayDirectionMapCoords.z() / factory->mZValueScale );
287 adjustedRayDirection.normalize();
288
289 QgsPointCloudIndex *index = factory->mPointCloudIndex;
290
291 const QgsPointCloudAttributeCollection attributeCollection = index->attributes();
292 QgsPointCloudRequest request;
293 request.setAttributes( attributeCollection );
294
295 double minDist = -1.;
296 const QList<QgsChunkNode *> activeNodes = this->activeNodes();
297 for ( QgsChunkNode *node : activeNodes )
298 {
299 const QgsChunkNodeId id = node->tileId();
300 const IndexedPointCloudNode n( id.d, id.x, id.y, id.z );
301
302 if ( !index->hasNode( n ) )
303 continue;
304
305 if ( !QgsRayCastingUtils::rayBoxIntersection( ray, node->bbox() ) )
306 continue;
307
308 std::unique_ptr<QgsPointCloudBlock> block( index->nodeData( n, request ) );
309 if ( !block )
310 continue;
311
312 const QgsVector3D blockScale = block->scale();
313 const QgsVector3D blockOffset = block->offset();
314
315 const char *ptr = block->data();
316 const QgsPointCloudAttributeCollection blockAttributes = block->attributes();
317 const std::size_t recordSize = blockAttributes.pointRecordSize();
318 int xOffset = 0, yOffset = 0, zOffset = 0;
319 const QgsPointCloudAttribute::DataType xType = blockAttributes.find( QStringLiteral( "X" ), xOffset )->type();
320 const QgsPointCloudAttribute::DataType yType = blockAttributes.find( QStringLiteral( "Y" ), yOffset )->type();
321 const QgsPointCloudAttribute::DataType zType = blockAttributes.find( QStringLiteral( "Z" ), zOffset )->type();
322 for ( int i = 0; i < block->pointCount(); ++i )
323 {
324 double x, y, z;
325 QgsPointCloudAttribute::getPointXYZ( ptr, i, recordSize, xOffset, xType, yOffset, yType, zOffset, zType, blockScale, blockOffset, x, y, z );
326 const QgsVector3D point( x, y, z );
327
328 // check whether point is in front of the ray
329 // similar to QgsRay3D::isInFront(), but using doubles
330 QgsVector3D vectorToPoint = point - adjustedRayOrigin;
331 vectorToPoint.normalize();
332 if ( QgsVector3D::dotProduct( vectorToPoint, adjustedRayDirection ) < 0.0 )
333 continue;
334
335 // calculate the angle between the point and the projected point
336 // similar to QgsRay3D::angleToPoint(), but using doubles
337 const QgsVector3D projPoint = adjustedRayOrigin + adjustedRayDirection * QgsVector3D::dotProduct( point - adjustedRayOrigin, adjustedRayDirection );
338 const QgsVector3D v1 = projPoint - adjustedRayOrigin ;
339 const QgsVector3D v2 = point - projPoint;
340 double angle = std::atan2( v2.length(), v1.length() ) * 180 / M_PI;
341 if ( angle > limitAngle )
342 continue;
343
344 const double dist = rayOriginMapCoords.distance( point );
345
346 if ( minDist < 0 || dist < minDist )
347 {
348 minDist = dist;
349 }
350 else if ( context.singleResult )
351 {
352 continue;
353 }
354
355 // Note : applying elevation properties is done in fromPointCloudIdentificationToIdentifyResults
356 QVariantMap pointAttr = QgsPointCloudAttribute::getAttributeMap( ptr, i * recordSize, blockAttributes );
357 pointAttr[ QStringLiteral( "X" ) ] = x;
358 pointAttr[ QStringLiteral( "Y" ) ] = y;
359 pointAttr[ QStringLiteral( "Z" ) ] = z;
360
361 const QgsVector3D worldPoint = factory->mMap.mapToWorldCoordinates( point );
362 QgsRayCastingUtils::RayHit hit( dist, worldPoint.toVector3D(), FID_NULL, pointAttr );
363 if ( context.singleResult )
364 result.clear();
365 result.append( hit );
366 }
367 }
368 return result;
369}
Represents a indexed point cloud node in octree.
QgsVector3D origin() const
Returns coordinates in map CRS at which 3D scene has origin (0,0,0)
static QgsVector3D mapToWorldCoordinates(const QgsVector3D &mapCoords, const QgsVector3D &origin)
Converts map coordinates to 3D world coordinates (applies offset and turns (x,y,z) into (x,...
Definition: qgs3dutils.cpp:551
3
Definition: qgsaabb.h:34
float yMax
Definition: qgsaabb.h:88
float xMax
Definition: qgsaabb.h:87
float xCenter() const
Returns center in X axis.
Definition: qgsaabb.h:50
float xMin
Definition: qgsaabb.h:84
float zMax
Definition: qgsaabb.h:89
float yMin
Definition: qgsaabb.h:85
float yCenter() const
Returns center in Y axis.
Definition: qgsaabb.h:52
float zMin
Definition: qgsaabb.h:86
float zCenter() const
Returns center in Z axis.
Definition: qgsaabb.h:54
QgsPointCloudCategoryList getFilteredOutCategories() const
Gets the list of categories of the classification that should not be rendered.
Class for doing transforms between two map coordinate systems.
void setBallparkTransformsAreAppropriate(bool appropriate)
Sets whether approximate "ballpark" results are appropriate for this coordinate transform.
QgsPointXY transform(const QgsPointXY &point, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward) const SIP_THROW(QgsCsException)
Transform the point from the source CRS to the destination CRS.
Custom exception class for Coordinate Reference System related exceptions.
Definition: qgsexception.h:67
bool renderAsTriangles() const
Returns whether points are triangulated to render solid surface.
float pointSize() const
Returns the point size of the point cloud.
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.
DataType
Systems of unit measurement.
static void getPointXYZ(const char *ptr, int i, std::size_t pointRecordSize, int xOffset, QgsPointCloudAttribute::DataType xType, int yOffset, QgsPointCloudAttribute::DataType yType, int zOffset, QgsPointCloudAttribute::DataType zType, const QgsVector3D &indexScale, const QgsVector3D &indexOffset, double &x, double &y, double &z)
Retrieves the x, y, z values for the point at index i.
static QVariantMap getAttributeMap(const char *data, std::size_t recordOffset, const QgsPointCloudAttributeCollection &attributeCollection)
Retrieves all the attributes of a point.
DataType type() const
Returns the data type.
Represents packaged data bounds.
qint32 xMax() const
Returns x max.
qint32 xMin() const
Returns x min.
qint32 yMax() const
Returns y max.
qint32 zMax() const
Returns z max.
qint32 yMin() const
Returns y min.
qint32 zMin() const
Returns z min.
Represents a indexed point clouds data in octree.
virtual bool hasNode(const IndexedPointCloudNode &n) const
Returns whether the octree contain given node.
virtual QgsPointCloudBlock * nodeData(const IndexedPointCloudNode &n, const QgsPointCloudRequest &request)=0
Returns node data block.
void setAttributes(const QgsPointCloudAttributeCollection &attributes)
Sets native attributes of the data.
QgsPointCloudAttributeCollection attributes() const
Returns all attributes that are stored in the file.
Point cloud data request.
void setAttributes(const QgsPointCloudAttributeCollection &attributes)
Set attributes filter in the request.
double y() const
Returns Y coordinate.
Definition: qgsvector3d.h:51
double z() const
Returns Z coordinate.
Definition: qgsvector3d.h:53
QVector3D toVector3D() const
Converts the current object to QVector3D.
Definition: qgsvector3d.h:169
static double dotProduct(const QgsVector3D &v1, const QgsVector3D &v2)
Returns the dot product of two vectors.
Definition: qgsvector3d.h:99
double x() const
Returns X coordinate.
Definition: qgsvector3d.h:49
void normalize()
Normalizes the current vector in place.
Definition: qgsvector3d.h:119
double length() const
Returns the length of the vector.
Definition: qgsvector3d.h:113
double ANALYSIS_EXPORT angle(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4)
Calculates the angle between two segments (in 2 dimension, z-values are ignored)
Definition: MathUtils.cpp:786
#define FID_NULL
Definition: qgsfeatureid.h:29
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
#define QgsDebugError(str)
Definition: qgslogger.h:38