QGIS API Documentation 3.34.0-Prizren (ffbdd678812)
Loading...
Searching...
No Matches
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 QgsChunkNode *node = new QgsChunkNode( QgsChunkNodeId( 0, 0, 0, 0 ), bbox, error );
180 node->setRefinementProcess( mSymbol->renderAsTriangles() ? Qgis::TileRefinementProcess::Replacement : Qgis::TileRefinementProcess::Additive );
181 return node;
182}
183
184QVector<QgsChunkNode *> QgsPointCloudLayerChunkLoaderFactory::createChildren( QgsChunkNode *node ) const
185{
186 QVector<QgsChunkNode *> children;
187 const QgsChunkNodeId nodeId = node->tileId();
188 const QgsAABB bbox = node->bbox();
189 const float childError = node->error() / 2;
190 float xc = bbox.xCenter(), yc = bbox.yCenter(), zc = bbox.zCenter();
191
192 for ( int i = 0; i < 8; ++i )
193 {
194 int dx = i & 1, dy = !!( i & 2 ), dz = !!( i & 4 );
195 const QgsChunkNodeId childId( nodeId.d + 1, nodeId.x * 2 + dx, nodeId.y * 2 + dy, nodeId.z * 2 + dz );
196
197 if ( !mPointCloudIndex->hasNode( IndexedPointCloudNode( childId.d, childId.x, childId.y, childId.z ) ) )
198 continue;
199 if ( !mExtent.isEmpty() &&
200 !mPointCloudIndex->nodeMapExtent( IndexedPointCloudNode( childId.d, childId.x, childId.y, childId.z ) ).intersects( mExtent ) )
201 continue;
202
203 // the Y and Z coordinates below are intentionally flipped, because
204 // in chunk node IDs the X,Y axes define horizontal plane,
205 // while in our 3D scene the X,Z axes define the horizontal plane
206 const float chXMin = dx ? xc : bbox.xMin;
207 const float chXMax = dx ? bbox.xMax : xc;
208 // Z axis: values are increasing to the south
209 const float chZMin = !dy ? zc : bbox.zMin;
210 const float chZMax = !dy ? bbox.zMax : zc;
211 const float chYMin = dz ? yc : bbox.yMin;
212 const float chYMax = dz ? bbox.yMax : yc;
213 QgsChunkNode *child = new QgsChunkNode( childId, QgsAABB( chXMin, chYMin, chZMin, chXMax, chYMax, chZMax ), childError, node );
214 child->setRefinementProcess( mSymbol->renderAsTriangles() ? Qgis::TileRefinementProcess::Replacement : Qgis::TileRefinementProcess::Additive );
215 children << child;
216 }
217 return children;
218}
219
221
222
223QgsAABB nodeBoundsToAABB( QgsPointCloudDataBounds nodeBounds, QgsVector3D offset, QgsVector3D scale, const Qgs3DMapSettings &map, const QgsCoordinateTransform &coordinateTransform, double zValueOffset )
224{
225 QgsVector3D extentMin3D( nodeBounds.xMin() * scale.x() + offset.x(), nodeBounds.yMin() * scale.y() + offset.y(), nodeBounds.zMin() * scale.z() + offset.z() + zValueOffset );
226 QgsVector3D extentMax3D( nodeBounds.xMax() * scale.x() + offset.x(), nodeBounds.yMax() * scale.y() + offset.y(), nodeBounds.zMax() * scale.z() + offset.z() + zValueOffset );
227 QgsCoordinateTransform extentTransform = coordinateTransform;
228 extentTransform.setBallparkTransformsAreAppropriate( true );
229 try
230 {
231 extentMin3D = extentTransform.transform( extentMin3D );
232 extentMax3D = extentTransform.transform( extentMax3D );
233 }
234 catch ( QgsCsException & )
235 {
236 QgsDebugError( QStringLiteral( "Error transforming node bounds coordinate" ) );
237 }
238 const QgsVector3D worldExtentMin3D = Qgs3DUtils::mapToWorldCoordinates( extentMin3D, map.origin() );
239 const QgsVector3D worldExtentMax3D = Qgs3DUtils::mapToWorldCoordinates( extentMax3D, map.origin() );
240 QgsAABB rootBbox( worldExtentMin3D.x(), worldExtentMin3D.y(), worldExtentMin3D.z(),
241 worldExtentMax3D.x(), worldExtentMax3D.y(), worldExtentMax3D.z() );
242 return rootBbox;
243}
244
245
246QgsPointCloudLayerChunkedEntity::QgsPointCloudLayerChunkedEntity( QgsPointCloudIndex *pc, const Qgs3DMapSettings &map,
247 const QgsCoordinateTransform &coordinateTransform, QgsPointCloud3DSymbol *symbol,
248 float maximumScreenSpaceError, bool showBoundingBoxes,
249 double zValueScale, double zValueOffset,
250 int pointBudget )
251 : QgsChunkedEntity( maximumScreenSpaceError,
252 new QgsPointCloudLayerChunkLoaderFactory( map, coordinateTransform, pc, symbol, zValueScale, zValueOffset, pointBudget ), true, pointBudget )
253{
254 setShowBoundingBoxes( showBoundingBoxes );
255}
256
257QgsPointCloudLayerChunkedEntity::~QgsPointCloudLayerChunkedEntity()
258{
259 // cancel / wait for jobs
260 cancelActiveJobs();
261}
262
263QVector<QgsRayCastingUtils::RayHit> QgsPointCloudLayerChunkedEntity::rayIntersection( const QgsRayCastingUtils::Ray3D &ray, const QgsRayCastingUtils::RayCastContext &context ) const
264{
265 QVector<QgsRayCastingUtils::RayHit> result;
266 QgsPointCloudLayerChunkLoaderFactory *factory = static_cast<QgsPointCloudLayerChunkLoaderFactory *>( mChunkLoaderFactory );
267
268 // transform ray
269 const QgsVector3D rayOriginMapCoords = factory->mMap.worldToMapCoordinates( ray.origin() );
270 const QgsVector3D pointMapCoords = factory->mMap.worldToMapCoordinates( ray.origin() + ray.origin().length() * ray.direction().normalized() );
271 QgsVector3D rayDirectionMapCoords = pointMapCoords - rayOriginMapCoords;
272 rayDirectionMapCoords.normalize();
273
274 const int screenSizePx = std::max( context.screenSize.height(), context.screenSize.width() );
275
276 const QgsPointCloud3DSymbol *symbol = factory->mSymbol.get();
277 // Symbol can be null in case of no rendering enabled
278 if ( !symbol )
279 return result;
280 const double pointSize = symbol->pointSize();
281
282 // We're using the angle as a tolerance, effectively meaning we're fetching points intersecting a cone.
283 // This may be revisited to use a cylinder instead, if the balance between near/far points does not scale
284 // well with different point sizes, screen sizes and fov values.
285 const double limitAngle = 2. * pointSize / screenSizePx * factory->mMap.fieldOfView();
286
287 // adjust ray to elevation properties
288 const QgsVector3D adjustedRayOrigin = QgsVector3D( rayOriginMapCoords.x(), rayOriginMapCoords.y(), ( rayOriginMapCoords.z() - factory->mZValueOffset ) / factory->mZValueScale );
289 QgsVector3D adjustedRayDirection = QgsVector3D( rayDirectionMapCoords.x(), rayDirectionMapCoords.y(), rayDirectionMapCoords.z() / factory->mZValueScale );
290 adjustedRayDirection.normalize();
291
292 QgsPointCloudIndex *index = factory->mPointCloudIndex;
293
294 const QgsPointCloudAttributeCollection attributeCollection = index->attributes();
295 QgsPointCloudRequest request;
296 request.setAttributes( attributeCollection );
297
298 double minDist = -1.;
299 const QList<QgsChunkNode *> activeNodes = this->activeNodes();
300 for ( QgsChunkNode *node : activeNodes )
301 {
302 const QgsChunkNodeId id = node->tileId();
303 const IndexedPointCloudNode n( id.d, id.x, id.y, id.z );
304
305 if ( !index->hasNode( n ) )
306 continue;
307
308 if ( !QgsRayCastingUtils::rayBoxIntersection( ray, node->bbox() ) )
309 continue;
310
311 std::unique_ptr<QgsPointCloudBlock> block( index->nodeData( n, request ) );
312 if ( !block )
313 continue;
314
315 const QgsVector3D blockScale = block->scale();
316 const QgsVector3D blockOffset = block->offset();
317
318 const char *ptr = block->data();
319 const QgsPointCloudAttributeCollection blockAttributes = block->attributes();
320 const std::size_t recordSize = blockAttributes.pointRecordSize();
321 int xOffset = 0, yOffset = 0, zOffset = 0;
322 const QgsPointCloudAttribute::DataType xType = blockAttributes.find( QStringLiteral( "X" ), xOffset )->type();
323 const QgsPointCloudAttribute::DataType yType = blockAttributes.find( QStringLiteral( "Y" ), yOffset )->type();
324 const QgsPointCloudAttribute::DataType zType = blockAttributes.find( QStringLiteral( "Z" ), zOffset )->type();
325 for ( int i = 0; i < block->pointCount(); ++i )
326 {
327 double x, y, z;
328 QgsPointCloudAttribute::getPointXYZ( ptr, i, recordSize, xOffset, xType, yOffset, yType, zOffset, zType, blockScale, blockOffset, x, y, z );
329 const QgsVector3D point( x, y, z );
330
331 // check whether point is in front of the ray
332 // similar to QgsRay3D::isInFront(), but using doubles
333 QgsVector3D vectorToPoint = point - adjustedRayOrigin;
334 vectorToPoint.normalize();
335 if ( QgsVector3D::dotProduct( vectorToPoint, adjustedRayDirection ) < 0.0 )
336 continue;
337
338 // calculate the angle between the point and the projected point
339 // similar to QgsRay3D::angleToPoint(), but using doubles
340 const QgsVector3D projPoint = adjustedRayOrigin + adjustedRayDirection * QgsVector3D::dotProduct( point - adjustedRayOrigin, adjustedRayDirection );
341 const QgsVector3D v1 = projPoint - adjustedRayOrigin ;
342 const QgsVector3D v2 = point - projPoint;
343 double angle = std::atan2( v2.length(), v1.length() ) * 180 / M_PI;
344 if ( angle > limitAngle )
345 continue;
346
347 const double dist = rayOriginMapCoords.distance( point );
348
349 if ( minDist < 0 || dist < minDist )
350 {
351 minDist = dist;
352 }
353 else if ( context.singleResult )
354 {
355 continue;
356 }
357
358 // Note : applying elevation properties is done in fromPointCloudIdentificationToIdentifyResults
359 QVariantMap pointAttr = QgsPointCloudAttribute::getAttributeMap( ptr, i * recordSize, blockAttributes );
360 pointAttr[ QStringLiteral( "X" ) ] = x;
361 pointAttr[ QStringLiteral( "Y" ) ] = y;
362 pointAttr[ QStringLiteral( "Z" ) ] = z;
363
364 const QgsVector3D worldPoint = factory->mMap.mapToWorldCoordinates( point );
365 QgsRayCastingUtils::RayHit hit( dist, worldPoint.toVector3D(), FID_NULL, pointAttr );
366 if ( context.singleResult )
367 result.clear();
368 result.append( hit );
369 }
370 }
371 return result;
372}
Represents a indexed point cloud node in octree.
The Qgis class provides global constants for use throughout the application.
Definition qgis.h:54
@ Replacement
When tile is refined then its children should be used in place of itself.
@ Reverse
Reverse/inverse transform (from destination to source)
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,...
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
Transform the point from the source CRS to the destination CRS.
Custom exception class for Coordinate Reference System related exceptions.
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 std::unique_ptr< 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.
Class for storage of 3D vectors similar to QVector3D, with the difference that it uses double precisi...
Definition qgsvector3d.h:32
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.
static double dotProduct(const QgsVector3D &v1, const QgsVector3D &v2)
Returns the dot product of two vectors.
double x() const
Returns X coordinate.
Definition qgsvector3d.h:49
void normalize()
Normalizes the current vector in place.
double length() const
Returns the length of the vector.
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)
#define FID_NULL
#define QgsDebugMsgLevel(str, level)
Definition qgslogger.h:39
#define QgsDebugError(str)
Definition qgslogger.h:38