QGIS API Documentation  3.23.0-Master (22c16f2067)
qgsdemterraintileloader_p.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsdemterraintileloader_p.cpp
3  --------------------------------------
4  Date : July 2017
5  Copyright : (C) 2017 by Martin Dobias
6  Email : wonder 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 "qgs3dmapsettings.h"
19 #include "qgschunknode_p.h"
20 #include "qgsdemterraingenerator.h"
22 #include "qgseventtracing.h"
24 #include "qgsterrainentity_p.h"
26 #include "qgsterraintileentity_p.h"
27 #include "qgsterraingenerator.h"
28 
29 #include <Qt3DRender/QGeometryRenderer>
30 #include <Qt3DCore/QTransform>
31 #include <QMutexLocker>
32 
34 
35 static void _heightMapMinMax( const QByteArray &heightMap, float &zMin, float &zMax )
36 {
37  const float *zBits = ( const float * ) heightMap.constData();
38  int zCount = heightMap.count() / sizeof( float );
39  bool first = true;
40 
41  zMin = zMax = std::numeric_limits<float>::quiet_NaN();
42  for ( int i = 0; i < zCount; ++i )
43  {
44  float z = zBits[i];
45  if ( std::isnan( z ) )
46  continue;
47  if ( first )
48  {
49  zMin = zMax = z;
50  first = false;
51  }
52  zMin = std::min( zMin, z );
53  zMax = std::max( zMax, z );
54  }
55 }
56 
57 
58 QgsDemTerrainTileLoader::QgsDemTerrainTileLoader( QgsTerrainEntity *terrain, QgsChunkNode *node, QgsTerrainGenerator *terrainGenerator )
59  : QgsTerrainTileLoader( terrain, node )
60  , mResolution( 0 )
61 {
62 
63  QgsDemHeightMapGenerator *heightMapGenerator = nullptr;
64  if ( terrainGenerator->type() == QgsTerrainGenerator::Dem )
65  {
66  QgsDemTerrainGenerator *generator = static_cast<QgsDemTerrainGenerator *>( terrainGenerator );
67  heightMapGenerator = generator->heightMapGenerator();
68  mSkirtHeight = generator->skirtHeight();
69  }
70  else if ( terrainGenerator->type() == QgsTerrainGenerator::Online )
71  {
72  QgsOnlineTerrainGenerator *generator = static_cast<QgsOnlineTerrainGenerator *>( terrainGenerator );
73  heightMapGenerator = generator->heightMapGenerator();
74  mSkirtHeight = generator->skirtHeight();
75  }
76  else
77  Q_ASSERT( false );
78 
79  // get heightmap asynchronously
80  connect( heightMapGenerator, &QgsDemHeightMapGenerator::heightMapReady, this, &QgsDemTerrainTileLoader::onHeightMapReady );
81  mHeightMapJobId = heightMapGenerator->render( node->tileId() );
82  mResolution = heightMapGenerator->resolution();
83 }
84 
85 Qt3DCore::QEntity *QgsDemTerrainTileLoader::createEntity( Qt3DCore::QEntity *parent )
86 {
87  float zMin, zMax;
88  _heightMapMinMax( mHeightMap, zMin, zMax );
89 
90  if ( std::isnan( zMin ) || std::isnan( zMax ) )
91  {
92  // no data available for this tile
93  return nullptr;
94  }
95 
96  const Qgs3DMapSettings &map = terrain()->map3D();
97  QgsChunkNodeId nodeId = mNode->tileId();
98  QgsRectangle extent = map.terrainGenerator()->tilingScheme().tileToExtent( nodeId );
99  double x0 = extent.xMinimum() - map.origin().x();
100  double y0 = extent.yMinimum() - map.origin().y();
101  double side = extent.width();
102  double half = side / 2;
103 
104 
105  QgsTerrainTileEntity *entity = new QgsTerrainTileEntity( nodeId );
106 
107  // create geometry renderer
108 
109  Qt3DRender::QGeometryRenderer *mesh = new Qt3DRender::QGeometryRenderer;
110  mesh->setGeometry( new DemTerrainTileGeometry( mResolution, side, map.terrainVerticalScale(), mSkirtHeight, mHeightMap, mesh ) );
111  entity->addComponent( mesh ); // takes ownership if the component has no parent
112 
113  // create material
114 
115  createTextureComponent( entity, map.isTerrainShadingEnabled(), map.terrainShadingMaterial(), !map.layers().empty() );
116 
117  // create transform
118 
119  Qt3DCore::QTransform *transform = nullptr;
120  transform = new Qt3DCore::QTransform();
121  entity->addComponent( transform );
122 
123  transform->setScale( side );
124  transform->setTranslation( QVector3D( x0 + half, 0, - ( y0 + half ) ) );
125 
126  mNode->setExactBbox( QgsAABB( x0, zMin * map.terrainVerticalScale(), -y0, x0 + side, zMax * map.terrainVerticalScale(), -( y0 + side ) ) );
127 
128  entity->setEnabled( false );
129  entity->setParent( parent );
130  return entity;
131 }
132 
133 void QgsDemTerrainTileLoader::onHeightMapReady( int jobId, const QByteArray &heightMap )
134 {
135  if ( mHeightMapJobId == jobId )
136  {
137  this->mHeightMap = heightMap;
138  mHeightMapJobId = -1;
139 
140  // continue loading - texture
141  loadTexture();
142  }
143 }
144 
145 
146 // ---------------------
147 
148 #include <qgsrasterlayer.h>
149 #include <qgsrasterprojector.h>
150 #include <QtConcurrent/QtConcurrentRun>
151 #include <QFutureWatcher>
152 #include "qgsterraindownloader.h"
153 
154 QgsDemHeightMapGenerator::QgsDemHeightMapGenerator( QgsRasterLayer *dtm, const QgsTilingScheme &tilingScheme, int resolution, const QgsCoordinateTransformContext &transformContext )
155  : mDtm( dtm )
156  , mClonedProvider( dtm ? ( QgsRasterDataProvider * )dtm->dataProvider()->clone() : nullptr )
157  , mTilingScheme( tilingScheme )
158  , mResolution( resolution )
159  , mLastJobId( 0 )
160  , mDownloader( dtm ? nullptr : new QgsTerrainDownloader( transformContext ) )
161  , mTransformContext( transformContext )
162 {
163 }
164 
165 QgsDemHeightMapGenerator::~QgsDemHeightMapGenerator()
166 {
167  delete mClonedProvider;
168 }
169 
170 
171 static QByteArray _readDtmData( QgsRasterDataProvider *provider, const QgsRectangle &extent, int res, const QgsCoordinateReferenceSystem &destCrs )
172 {
173  QgsEventTracing::ScopedEvent e( QStringLiteral( "3D" ), QStringLiteral( "DEM" ) );
174 
175  // TODO: use feedback object? (but GDAL currently does not support cancellation anyway)
176  QgsRasterInterface *input = provider;
177  std::unique_ptr<QgsRasterProjector> projector;
178  if ( provider->crs() != destCrs )
179  {
180  projector.reset( new QgsRasterProjector );
181  projector->setCrs( provider->crs(), destCrs, provider->transformContext() );
182  projector->setInput( provider );
183  input = projector.get();
184  }
185  std::unique_ptr< QgsRasterBlock > block( input->block( 1, extent, res, res ) );
186 
187  QByteArray data;
188  if ( block )
189  {
190  block->convert( Qgis::DataType::Float32 ); // currently we expect just floats
191  data = block->data();
192  data.detach(); // this should make a deep copy
193 
194  if ( block->hasNoData() )
195  {
196  // turn all no-data values into NaN in the output array
197  float *floatData = reinterpret_cast<float *>( data.data() );
198  Q_ASSERT( data.count() % sizeof( float ) == 0 );
199  int count = data.count() / sizeof( float );
200  for ( int i = 0; i < count; ++i )
201  {
202  if ( block->isNoData( i ) )
203  floatData[i] = std::numeric_limits<float>::quiet_NaN();
204  }
205  }
206  }
207  return data;
208 }
209 
210 static QByteArray _readOnlineDtm( QgsTerrainDownloader *downloader, const QgsRectangle &extent, int res, const QgsCoordinateReferenceSystem &destCrs, const QgsCoordinateTransformContext &context )
211 {
212  return downloader->getHeightMap( extent, res, destCrs, context );
213 }
214 
215 int QgsDemHeightMapGenerator::render( const QgsChunkNodeId &nodeId )
216 {
217  QgsEventTracing::addEvent( QgsEventTracing::AsyncBegin, QStringLiteral( "3D" ), QStringLiteral( "DEM" ), nodeId.text() );
218 
219  // extend the rect by half-pixel on each side? to get the values in "corners"
220  QgsRectangle extent = mTilingScheme.tileToExtent( nodeId );
221  float mapUnitsPerPixel = extent.width() / mResolution;
222  extent.grow( mapUnitsPerPixel / 2 );
223  // but make sure not to go beyond the full extent (returns invalid values)
224  QgsRectangle fullExtent = mTilingScheme.tileToExtent( 0, 0, 0 );
225  extent = extent.intersect( fullExtent );
226 
227  JobData jd;
228  jd.jobId = ++mLastJobId;
229  jd.tileId = nodeId;
230  jd.extent = extent;
231  jd.timer.start();
232  QFutureWatcher<QByteArray> *fw = new QFutureWatcher<QByteArray>( nullptr );
233  connect( fw, &QFutureWatcher<QByteArray>::finished, this, &QgsDemHeightMapGenerator::onFutureFinished );
234  connect( fw, &QFutureWatcher<QByteArray>::finished, fw, &QObject::deleteLater );
235  // make a clone of the data provider so it is safe to use in worker thread
236  if ( mDtm )
237  jd.future = QtConcurrent::run( _readDtmData, mClonedProvider, extent, mResolution, mTilingScheme.crs() );
238  else
239  jd.future = QtConcurrent::run( _readOnlineDtm, mDownloader.get(), extent, mResolution, mTilingScheme.crs(), mTransformContext );
240 
241  fw->setFuture( jd.future );
242 
243  mJobs.insert( fw, jd );
244 
245  return jd.jobId;
246 }
247 
248 void QgsDemHeightMapGenerator::waitForFinished()
249 {
250  for ( QFutureWatcher<QByteArray> *fw : mJobs.keys() )
251  {
252  disconnect( fw, &QFutureWatcher<QByteArray>::finished, this, &QgsDemHeightMapGenerator::onFutureFinished );
253  disconnect( fw, &QFutureWatcher<QByteArray>::finished, fw, &QObject::deleteLater );
254  }
255  QVector<QFutureWatcher<QByteArray>*> toBeDeleted;
256  for ( QFutureWatcher<QByteArray> *fw : mJobs.keys() )
257  {
258  fw->waitForFinished();
259  JobData jobData = mJobs.value( fw );
260  toBeDeleted.push_back( fw );
261 
262  QByteArray data = jobData.future.result();
263  emit heightMapReady( jobData.jobId, data );
264  }
265 
266  for ( QFutureWatcher<QByteArray> *fw : toBeDeleted )
267  {
268  mJobs.remove( fw );
269  fw->deleteLater();
270  }
271 }
272 
273 void QgsDemHeightMapGenerator::lazyLoadDtmCoarseData( int res, const QgsRectangle &rect )
274 {
275  QMutexLocker locker( &mLazyLoadDtmCoarseDataMutex );
276  if ( mDtmCoarseData.isEmpty() )
277  {
278  std::unique_ptr< QgsRasterBlock > block( mDtm->dataProvider()->block( 1, rect, res, res ) );
279  block->convert( Qgis::DataType::Float32 );
280  mDtmCoarseData = block->data();
281  mDtmCoarseData.detach(); // make a deep copy
282  }
283 }
284 
285 float QgsDemHeightMapGenerator::heightAt( double x, double y )
286 {
287  if ( !mDtm )
288  return 0; // TODO: calculate heights for online DTM
289 
290  // TODO: this is quite a primitive implementation: better to use heightmaps currently in use
291  int res = 1024;
292  QgsRectangle rect = mDtm->extent();
293  lazyLoadDtmCoarseData( res, rect );
294 
295  int cellX = ( int )( ( x - rect.xMinimum() ) / rect.width() * res + .5f );
296  int cellY = ( int )( ( rect.yMaximum() - y ) / rect.height() * res + .5f );
297  cellX = std::clamp( cellX, 0, res - 1 );
298  cellY = std::clamp( cellY, 0, res - 1 );
299 
300  const float *data = ( const float * ) mDtmCoarseData.constData();
301  return data[cellX + cellY * res];
302 }
303 
304 void QgsDemHeightMapGenerator::onFutureFinished()
305 {
306  QFutureWatcher<QByteArray> *fw = static_cast<QFutureWatcher<QByteArray>*>( sender() );
307  Q_ASSERT( fw );
308  Q_ASSERT( mJobs.contains( fw ) );
309  JobData jobData = mJobs.value( fw );
310 
311  mJobs.remove( fw );
312  fw->deleteLater();
313 
314  QgsEventTracing::addEvent( QgsEventTracing::AsyncEnd, QStringLiteral( "3D" ), QStringLiteral( "DEM" ), jobData.tileId.text() );
315 
316  QByteArray data = jobData.future.result();
317  emit heightMapReady( jobData.jobId, data );
318 }
319 
@ Float32
Thirty two bit floating point (float)
QgsTerrainGenerator * terrainGenerator() const
Returns terrain generator. It takes care of producing terrain tiles from the input data.
double terrainVerticalScale() const
Returns vertical scale (exaggeration) of terrain.
bool isTerrainShadingEnabled() const
Returns whether terrain shading is enabled.
QgsPhongMaterialSettings terrainShadingMaterial() const
Returns terrain shading material.
QList< QgsMapLayer * > layers() const
Returns the list of 3D map layers to be rendered in the scene.
QgsVector3D origin() const
Returns coordinates in map CRS at which 3D scene has origin (0,0,0)
3
Definition: qgsaabb.h:34
This class represents a coordinate reference system (CRS).
Contains information about the context in which a coordinate transform is executed.
QgsCoordinateTransformContext transformContext() const
Returns data provider coordinate transform context.
virtual QgsCoordinateReferenceSystem crs() const =0
Returns the coordinate system for the data source.
QgsDemHeightMapGenerator * heightMapGenerator()
Returns height map generator object - takes care of extraction of elevations from the layer)
float skirtHeight() const
Returns skirt height (in world units). Skirts at the edges of terrain tiles help hide cracks between ...
QgsDemHeightMapGenerator * heightMapGenerator()
Returns height map generator object - takes care of extraction of elevations from the layer)
float skirtHeight() const
Returns skirt height (in world units). Skirts at the edges of terrain tiles help hide cracks between ...
Base class for raster data providers.
bool setInput(QgsRasterInterface *input) override
Set input.
Base class for processing filters like renderers, reprojector, resampler etc.
virtual QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr)=0
Read block of data using given extent and size.
Represents a raster layer.
QgsRasterProjector implements approximate projection support for it calculates grid of points in sour...
A rectangle specified with double values.
Definition: qgsrectangle.h:42
double yMaximum() const SIP_HOLDGIL
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:193
double xMinimum() const SIP_HOLDGIL
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:188
double yMinimum() const SIP_HOLDGIL
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:198
double height() const SIP_HOLDGIL
Returns the height of the rectangle.
Definition: qgsrectangle.h:230
void grow(double delta)
Grows the rectangle in place by the specified amount.
Definition: qgsrectangle.h:296
double width() const SIP_HOLDGIL
Returns the width of the rectangle.
Definition: qgsrectangle.h:223
QgsRectangle intersect(const QgsRectangle &rect) const
Returns the intersection with the given rectangle.
Definition: qgsrectangle.h:333
QByteArray getHeightMap(const QgsRectangle &extentOrig, int res, const QgsCoordinateReferenceSystem &destCrs, const QgsCoordinateTransformContext &context=QgsCoordinateTransformContext(), QString tmpFilenameImg=QString(), QString tmpFilenameTif=QString())
For given extent and resolution (number of pixels for width/height) in specified CRS,...
@ Dem
Terrain is built from raster layer with digital elevation model.
@ Online
Terrain is built from downloaded tiles with digital elevation model.
const QgsTilingScheme & tilingScheme() const
Returns tiling scheme of the terrain.
virtual Type type() const =0
What texture generator implementation is this.
QgsRectangle tileToExtent(int x, int y, int z) const
Returns map coordinates of the extent of a tile.
double y() const
Returns Y coordinate.
Definition: qgsvector3d.h:51
double x() const
Returns X coordinate.
Definition: qgsvector3d.h:49