QGIS API Documentation  3.0.2-Girona (307d082)
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 "qgsterrainentity_p.h"
24 #include "qgsterraintileentity_p.h"
25 
26 #include <Qt3DRender/QGeometryRenderer>
27 
29 
30 static void _heightMapMinMax( const QByteArray &heightMap, float &zMin, float &zMax )
31 {
32  const float *zBits = ( const float * ) heightMap.constData();
33  int zCount = heightMap.count() / sizeof( float );
34  bool first = true;
35 
36  zMin = zMax = std::numeric_limits<float>::quiet_NaN();
37  for ( int i = 0; i < zCount; ++i )
38  {
39  float z = zBits[i];
40  if ( std::isnan( z ) )
41  continue;
42  if ( first )
43  {
44  zMin = zMax = z;
45  first = false;
46  }
47  zMin = qMin( zMin, z );
48  zMax = qMax( zMax, z );
49  }
50 }
51 
52 
53 QgsDemTerrainTileLoader::QgsDemTerrainTileLoader( QgsTerrainEntity *terrain, QgsChunkNode *node )
54  : QgsTerrainTileLoader( terrain, node )
55  , mResolution( 0 )
56 {
57 
58  const Qgs3DMapSettings &map = terrain->map3D();
59  QgsDemTerrainGenerator *generator = static_cast<QgsDemTerrainGenerator *>( map.terrainGenerator() );
60 
61  // get heightmap asynchronously
62  connect( generator->heightMapGenerator(), &QgsDemHeightMapGenerator::heightMapReady, this, &QgsDemTerrainTileLoader::onHeightMapReady );
63  mHeightMapJobId = generator->heightMapGenerator()->render( node->tileX(), node->tileY(), node->tileZ() );
64  mResolution = generator->heightMapGenerator()->resolution();
65  mSkirtHeight = generator->skirtHeight();
66 }
67 
68 Qt3DCore::QEntity *QgsDemTerrainTileLoader::createEntity( Qt3DCore::QEntity *parent )
69 {
70  float zMin, zMax;
71  _heightMapMinMax( mHeightMap, zMin, zMax );
72 
73  if ( std::isnan( zMin ) || std::isnan( zMax ) )
74  {
75  // no data available for this tile
76  return nullptr;
77  }
78 
79  QgsTerrainTileEntity *entity = new QgsTerrainTileEntity;
80 
81  // create geometry renderer
82 
83  Qt3DRender::QGeometryRenderer *mesh = new Qt3DRender::QGeometryRenderer;
84  mesh->setGeometry( new DemTerrainTileGeometry( mResolution, mSkirtHeight, mHeightMap, mesh ) );
85  entity->addComponent( mesh ); // takes ownership if the component has no parent
86 
87  // create material
88 
89  createTextureComponent( entity );
90 
91  // create transform
92 
93  Qt3DCore::QTransform *transform = nullptr;
94  transform = new Qt3DCore::QTransform();
95  entity->addComponent( transform );
96 
97  const Qgs3DMapSettings &map = terrain()->map3D();
98  QgsRectangle extent = map.terrainGenerator()->tilingScheme().tileToExtent( mNode->tileX(), mNode->tileY(), mNode->tileZ() ); //node->extent;
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  transform->setScale3D( QVector3D( side, map.terrainVerticalScale(), side ) );
105  transform->setTranslation( QVector3D( x0 + half, 0, - ( y0 + half ) ) );
106 
107  mNode->setExactBbox( QgsAABB( x0, zMin * map.terrainVerticalScale(), -y0, x0 + side, zMax * map.terrainVerticalScale(), -( y0 + side ) ) );
108 
109  entity->setEnabled( false );
110  entity->setParent( parent );
111  return entity;
112 }
113 
114 void QgsDemTerrainTileLoader::onHeightMapReady( int jobId, const QByteArray &heightMap )
115 {
116  if ( mHeightMapJobId == jobId )
117  {
118  this->mHeightMap = heightMap;
119  mHeightMapJobId = -1;
120 
121  // continue loading - texture
122  loadTexture();
123  }
124 }
125 
126 
127 // ---------------------
128 
129 #include <qgsrasterlayer.h>
130 #include <qgsrasterprojector.h>
131 #include <QtConcurrent/QtConcurrentRun>
132 #include <QFutureWatcher>
133 
134 QgsDemHeightMapGenerator::QgsDemHeightMapGenerator( QgsRasterLayer *dtm, const QgsTilingScheme &tilingScheme, int resolution )
135  : mDtm( dtm )
136  , mClonedProvider( ( QgsRasterDataProvider * )dtm->dataProvider()->clone() )
137  , mTilingScheme( tilingScheme )
138  , mResolution( resolution )
139  , mLastJobId( 0 )
140 {
141 }
142 
143 QgsDemHeightMapGenerator::~QgsDemHeightMapGenerator()
144 {
145  delete mClonedProvider;
146 }
147 
148 #include <QElapsedTimer>
149 
150 static QByteArray _readDtmData( QgsRasterDataProvider *provider, const QgsRectangle &extent, int res, const QgsCoordinateReferenceSystem &destCrs )
151 {
152  QElapsedTimer t;
153  t.start();
154 
155  // TODO: use feedback object? (but GDAL currently does not support cancelation anyway)
156  QgsRasterInterface *input = provider;
157  std::unique_ptr<QgsRasterProjector> projector;
158  if ( provider->crs() != destCrs )
159  {
160  projector.reset( new QgsRasterProjector );
161  projector->setCrs( provider->crs(), destCrs );
162  projector->setInput( provider );
163  input = projector.get();
164  }
165  QgsRasterBlock *block = input->block( 1, extent, res, res );
166 
167  QByteArray data;
168  if ( block )
169  {
170  block->convert( Qgis::Float32 ); // currently we expect just floats
171  data = block->data();
172  data.detach(); // this should make a deep copy
173 
174  if ( block->hasNoData() )
175  {
176  // turn all no-data values into NaN in the output array
177  float *floatData = reinterpret_cast<float *>( data.data() );
178  Q_ASSERT( data.count() % sizeof( float ) == 0 );
179  int count = data.count() / sizeof( float );
180  for ( int i = 0; i < count; ++i )
181  {
182  if ( block->isNoData( i ) )
183  floatData[i] = std::numeric_limits<float>::quiet_NaN();
184  }
185  }
186 
187  delete block;
188  }
189  return data;
190 }
191 
192 int QgsDemHeightMapGenerator::render( int x, int y, int z )
193 {
194  Q_ASSERT( mJobs.isEmpty() ); // should be always just one active job...
195 
196  // extend the rect by half-pixel on each side? to get the values in "corners"
197  QgsRectangle extent = mTilingScheme.tileToExtent( x, y, z );
198  float mapUnitsPerPixel = extent.width() / mResolution;
199  extent.grow( mapUnitsPerPixel / 2 );
200  // but make sure not to go beyond the full extent (returns invalid values)
201  QgsRectangle fullExtent = mTilingScheme.tileToExtent( 0, 0, 0 );
202  extent = extent.intersect( &fullExtent );
203 
204  JobData jd;
205  jd.jobId = ++mLastJobId;
206  jd.extent = extent;
207  jd.timer.start();
208  // make a clone of the data provider so it is safe to use in worker thread
209  jd.future = QtConcurrent::run( _readDtmData, mClonedProvider, extent, mResolution, mTilingScheme.crs() );
210 
211  QFutureWatcher<QByteArray> *fw = new QFutureWatcher<QByteArray>( nullptr );
212  fw->setFuture( jd.future );
213  connect( fw, &QFutureWatcher<QByteArray>::finished, this, &QgsDemHeightMapGenerator::onFutureFinished );
214 
215  mJobs.insert( fw, jd );
216 
217  return jd.jobId;
218 }
219 
220 QByteArray QgsDemHeightMapGenerator::renderSynchronously( int x, int y, int z )
221 {
222  // extend the rect by half-pixel on each side? to get the values in "corners"
223  QgsRectangle extent = mTilingScheme.tileToExtent( x, y, z );
224  float mapUnitsPerPixel = extent.width() / mResolution;
225  extent.grow( mapUnitsPerPixel / 2 );
226  // but make sure not to go beyond the full extent (returns invalid values)
227  QgsRectangle fullExtent = mTilingScheme.tileToExtent( 0, 0, 0 );
228  extent = extent.intersect( &fullExtent );
229 
230  QgsRasterBlock *block = mDtm->dataProvider()->block( 1, extent, mResolution, mResolution );
231 
232  QByteArray data;
233  if ( block )
234  {
235  block->convert( Qgis::Float32 ); // currently we expect just floats
236  data = block->data();
237  data.detach(); // this should make a deep copy
238  delete block;
239  }
240 
241  return data;
242 }
243 
244 float QgsDemHeightMapGenerator::heightAt( double x, double y )
245 {
246  // TODO: this is quite a primitive implementation: better to use heightmaps currently in use
247  int res = 1024;
248  QgsRectangle rect = mDtm->extent();
249  if ( mDtmCoarseData.isEmpty() )
250  {
251  QgsRasterBlock *block = mDtm->dataProvider()->block( 1, rect, res, res );
252  block->convert( Qgis::Float32 );
253  mDtmCoarseData = block->data();
254  mDtmCoarseData.detach(); // make a deep copy
255  delete block;
256  }
257 
258  int cellX = ( int )( ( x - rect.xMinimum() ) / rect.width() * res + .5f );
259  int cellY = ( int )( ( rect.yMaximum() - y ) / rect.height() * res + .5f );
260  cellX = qBound( 0, cellX, res - 1 );
261  cellY = qBound( 0, cellY, res - 1 );
262 
263  const float *data = ( const float * ) mDtmCoarseData.constData();
264  return data[cellX + cellY * res];
265 }
266 
267 void QgsDemHeightMapGenerator::onFutureFinished()
268 {
269  QFutureWatcher<QByteArray> *fw = static_cast<QFutureWatcher<QByteArray>*>( sender() );
270  Q_ASSERT( fw );
271  Q_ASSERT( mJobs.contains( fw ) );
272  JobData jobData = mJobs.value( fw );
273 
274  mJobs.remove( fw );
275  fw->deleteLater();
276 
277  QByteArray data = jobData.future.result();
278  emit heightMapReady( jobData.jobId, data );
279 }
280 
3 Axis-aligned bounding box - in world coords.
Definition: qgsaabb.h:30
A rectangle specified with double values.
Definition: qgsrectangle.h:39
bool setInput(QgsRasterInterface *input) override
Set input.
QgsTerrainGenerator * clone() const override
Makes a copy of the current instance.
This class provides qgis with the ability to render raster datasets onto the mapcanvas.
QgsRectangle intersect(const QgsRectangle *rect) const
Return the intersection with the given rectangle.
int resolution() const
Returns resolution of the generator (how many elevation samples on one side of a terrain tile) ...
Thirty two bit floating point (float)
Definition: qgis.h:99
double y() const
Returns Y coordinate.
Definition: qgsvector3d.h:44
3 Definition of the world
bool isNoData(int row, int column)
Check if value at position is no data.
QgsVector3D origin() const
Returns coordinates in map CRS at which 3D scene has origin (0,0,0)
bool convert(Qgis::DataType destDataType)
Convert data to different type.
Raster data container.
void grow(double delta)
Grows the rectangle in place by the specified amount.
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.
virtual QgsCoordinateReferenceSystem crs() const =0
Returns the coordinate system for the data source.
double width() const
Returns the width of the rectangle.
Definition: qgsrectangle.h:142
float skirtHeight() const
Returns skirt height (in world units). Skirts at the edges of terrain tiles help hide cracks between ...
QgsRectangle extent() const override
extent of the terrain in terrain&#39;s CRS
const QgsTilingScheme & tilingScheme() const
Returns tiling scheme of the terrain.
Base class for processing filters like renderers, reprojector, resampler etc.
QByteArray data() const
Get access to raw data.
double terrainVerticalScale() const
Returns vertical scale (exaggeration) of terrain.
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:130
QgsRasterProjector implements approximate projection support for it calculates grid of points in sour...
QgsDemHeightMapGenerator * heightMapGenerator()
Returns height map generator object - takes care of extraction of elevations from the layer) ...
3 Implementation of terrain generator that uses a raster layer with DEM to build terrain.
This class represents a coordinate reference system (CRS).
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:120
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:125
bool hasNoData() const
Returns true if the block may contain no data.
3 The class encapsulates tiling scheme (just like with WMTS / TMS / XYZ layers).
QgsRectangle tileToExtent(int x, int y, int z) const
Returns map coordinates of the extent of a tile.
double x() const
Returns X coordinate.
Definition: qgsvector3d.h:42
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:149
Base class for raster data providers.
QgsTerrainGenerator * terrainGenerator() const
Returns terrain generator. It takes care of producing terrain tiles from the input data...