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