[3d] Respect no-data values from DEMs in terrain generator (fixes #17558)

Author: wonder-sk

src/3d/chunks/qgschunkedentity_p.cpp

@@ -304,6 +304,9 @@ void QgsChunkedEntity::requestResidency( QgsChunkNode *node )
   }
   else if ( node->state() == QgsChunkNode::Skeleton )
   {
+    if ( !node->hasData() )
+      return;   // no need to load (we already tried but got nothing back)
+
     // add to the loading queue
     QgsChunkListEntry *entry = new QgsChunkListEntry( node );
     node->setQueuedForLoad( entry );
@@ -332,10 +335,18 @@ void QgsChunkedEntity::onActiveJobFinished()
     // mark as loaded + create entity
     Qt3DCore::QEntity *entity = node->loader()->createEntity( this );
 
-    // load into node (should be in main thread again)
-    node->setLoaded( entity );
+    if ( entity )
+    {
+      // load into node (should be in main thread again)
+      node->setLoaded( entity );
 
-    mReplacementQueue->insertFirst( node->replacementQueueEntry() );
+      mReplacementQueue->insertFirst( node->replacementQueueEntry() );
+    }
+    else
+    {
+      node->setHasData( false );
+      node->cancelLoading();
+    }
 
     // now we need an update!
     mNeedsUpdate = true;

src/3d/chunks/qgschunknode_p.cpp

@@ -60,7 +60,7 @@ bool QgsChunkNode::allChildChunksResident( const QTime &currentTime ) const
   {
     if ( !mChildren[i] )
       return false;  // not even a skeleton
-    if ( !mChildren[i]->mEntity )
+    if ( mChildren[i]->mHasData && !mChildren[i]->mEntity )
       return false;  // no there yet
     Q_UNUSED( currentTime ); // seems we do not need this extra time (it just brings extra problems)
     //if (children[i]->entityCreatedTime.msecsTo(currentTime) < 100)

src/3d/chunks/qgschunknode_p.h

@@ -172,6 +172,11 @@ class QgsChunkNode
     //! called when bounding box
     void setExactBbox( const QgsAABB &box );
 
+    //! Sets whether the node has any data to be displayed. Can be used to set to false after load returned no data
+    void setHasData( bool hasData ) { mHasData = hasData; }
+    //! Returns whether the node has any data to be displayed. If not, it will be kept as a skeleton node and will not get loaded anymore
+    bool hasData() const { return mHasData; }
+
   private:
     QgsAABB mBbox;      //!< Bounding box in world coordinates
     float mError;    //!< Error of the node in world coordinates
@@ -193,6 +198,7 @@ class QgsChunkNode
     QgsChunkQueueJob *mUpdater;                //!< Object that does update of the chunk (not null <=> Updating state)
 
     QTime mEntityCreatedTime;
+    bool mHasData = true;   //!< Whether there are (will be) any data in this node (or any descentants) and so whether it makes sense to load this node
 };
 
 /// @endcond

src/3d/terrain/qgsdemterraintilegeometry_p.cpp

@@ -18,6 +18,7 @@
 #include <Qt3DRender/qbuffer.h>
 #include <Qt3DRender/qbufferdatagenerator.h>
 #include <limits>
+#include <cmath>
 
 ///@cond PRIVATE
 
@@ -51,6 +52,10 @@ static QByteArray createPlaneVertexData( int res, float skirtHeight, const QByte
   const float du = 1.0 / ( resolution.width() - 1 );
   const float dv = 1.0 / ( resolution.height() - 1 );
 
+  // the height of vertices with no-data value... the value should not really matter
+  // as we do not create valid triangles that would use such vertices
+  const float noDataHeight = 0;
+
   // Iterate over z
   for ( int j = -1; j <= resolution.height(); ++j )
   {
@@ -71,6 +76,9 @@ static QByteArray createPlaneVertexData( int res, float skirtHeight, const QByte
       else
         height = zData[ jBound * resolution.width() + iBound ] - skirtHeight;
 
+      if ( std::isnan( height ) )
+        height = noDataHeight;
+
       // position
       *fptr++ = x;
       *fptr++ = height;
@@ -91,8 +99,23 @@ static QByteArray createPlaneVertexData( int res, float skirtHeight, const QByte
   return bufferBytes;
 }
 
+inline int ijToHeightMapIndex( int i, int j, int numVerticesX, int numVerticesZ )
+{
+  i = qBound( 1, i, numVerticesX - 1 ) - 1;
+  j = qBound( 1, j, numVerticesZ - 1 ) - 1;
+  return j * ( numVerticesX - 2 ) + i;
+}
+
+
+static bool hasNoData( int i, int j, const float *heightMap, int numVerticesX, int numVerticesZ )
+{
+  return std::isnan( heightMap[ ijToHeightMapIndex( i, j, numVerticesX, numVerticesZ ) ] ) ||
+         std::isnan( heightMap[ ijToHeightMapIndex( i + 1, j, numVerticesX, numVerticesZ ) ] ) ||
+         std::isnan( heightMap[ ijToHeightMapIndex( i, j + 1, numVerticesX, numVerticesZ ) ] ) ||
+         std::isnan( heightMap[ ijToHeightMapIndex( i + 1, j + 1, numVerticesX, numVerticesZ ) ] );
+}
 
-static QByteArray createPlaneIndexData( int res )
+static QByteArray createPlaneIndexData( int res, const QByteArray &heightMap )
 {
   QSize resolution( res, res );
   int numVerticesX = resolution.width() + 2;
@@ -106,6 +129,8 @@ static QByteArray createPlaneIndexData( int res )
   indexBytes.resize( indices * sizeof( quint32 ) );
   quint32 *indexPtr = reinterpret_cast<quint32 *>( indexBytes.data() );
 
+  const float *heightMapFloat = reinterpret_cast<const float *>( heightMap.constData() );
+
   // Iterate over z
   for ( int j = 0; j < numVerticesZ - 1; ++j )
   {
@@ -115,6 +140,20 @@ static QByteArray createPlaneIndexData( int res )
     // Iterate over x
     for ( int i = 0; i < numVerticesX - 1; ++i )
     {
+      if ( hasNoData( i, j, heightMapFloat, numVerticesX, numVerticesZ ) )
+      {
+        // at least one corner of the quad has no-data value
+        // so let's make two invalid triangles
+        *indexPtr++ = rowStartIndex + i;
+        *indexPtr++ = rowStartIndex + i;
+        *indexPtr++ = rowStartIndex + i;
+
+        *indexPtr++ = rowStartIndex + i;
+        *indexPtr++ = rowStartIndex + i;
+        *indexPtr++ = rowStartIndex + i;
+        continue;
+      }
+
       // Split quad into two triangles
       *indexPtr++ = rowStartIndex + i;
       *indexPtr++ = nextRowStartIndex + i;
@@ -169,13 +208,14 @@ class PlaneVertexBufferFunctor : public QBufferDataGenerator
 class PlaneIndexBufferFunctor : public QBufferDataGenerator
 {
   public:
-    explicit PlaneIndexBufferFunctor( int resolution )
+    explicit PlaneIndexBufferFunctor( int resolution, const QByteArray &heightMap )
       : mResolution( resolution )
+      , mHeightMap( heightMap )
     {}
 
     QByteArray operator()() final
     {
-      return createPlaneIndexData( mResolution );
+      return createPlaneIndexData( mResolution, mHeightMap );
     }
 
     bool operator ==( const QBufferDataGenerator &other ) const final
@@ -190,6 +230,7 @@ class PlaneIndexBufferFunctor : public QBufferDataGenerator
 
   private:
     int mResolution;
+    QByteArray mHeightMap;
 };
 
 
@@ -254,7 +295,7 @@ void DemTerrainTileGeometry::init()
   mIndexAttribute->setCount( faces * 3 );
 
   mVertexBuffer->setDataGenerator( QSharedPointer<PlaneVertexBufferFunctor>::create( mResolution, mSkirtHeight, mHeightMap ) );
-  mIndexBuffer->setDataGenerator( QSharedPointer<PlaneIndexBufferFunctor>::create( mResolution ) );
+  mIndexBuffer->setDataGenerator( QSharedPointer<PlaneIndexBufferFunctor>::create( mResolution, mHeightMap ) );
 
   addAttribute( mPositionAttribute );
   addAttribute( mTexCoordAttribute );

src/3d/terrain/qgsdemterraintileloader_p.cpp

@@ -31,14 +31,23 @@ static void _heightMapMinMax( const QByteArray &heightMap, float &zMin, float &z
 {
   const float *zBits = ( const float * ) heightMap.constData();
   int zCount = heightMap.count() / sizeof( float );
-  zMin = zBits[0];
-  zMax = zBits[0];
+  bool first = true;
   for ( int i = 0; i < zCount; ++i )
   {
     float z = zBits[i];
+    if ( std::isnan( z ) )
+      continue;
+    if ( first )
+    {
+      zMin = zMax = z;
+      first = false;
+    }
     zMin = qMin( zMin, z );
     zMax = qMax( zMax, z );
   }
+
+  if ( first )
+    zMin = zMax = std::numeric_limits<float>::quiet_NaN();
 }
 
 
@@ -59,6 +68,15 @@ QgsDemTerrainTileLoader::QgsDemTerrainTileLoader( QgsTerrainEntity *terrain, Qgs
 
 Qt3DCore::QEntity *QgsDemTerrainTileLoader::createEntity( Qt3DCore::QEntity *parent )
 {
+  float zMin, zMax;
+  _heightMapMinMax( mHeightMap, zMin, zMax );
+
+  if ( std::isnan( zMin ) || std::isnan( zMax ) )
+  {
+    // no data available for this tile
+    return nullptr;
+  }
+
   QgsTerrainTileEntity *entity = new QgsTerrainTileEntity;
 
   // create geometry renderer
@@ -77,9 +95,6 @@ Qt3DCore::QEntity *QgsDemTerrainTileLoader::createEntity( Qt3DCore::QEntity *par
   transform = new Qt3DCore::QTransform();
   entity->addComponent( transform );
 
-  float zMin, zMax;
-  _heightMapMinMax( mHeightMap, zMin, zMax );
-
   const Qgs3DMapSettings &map = terrain()->map3D();
   QgsRectangle extent = map.terrainGenerator()->tilingScheme().tileToExtent( mNode->tileX(), mNode->tileY(), mNode->tileZ() ); //node->extent;
   double x0 = extent.xMinimum() - map.origin().x();
@@ -156,6 +171,20 @@ static QByteArray _readDtmData( QgsRasterDataProvider *provider, const QgsRectan
     block->convert( Qgis::Float32 ); // currently we expect just floats
     data = block->data();
     data.detach();  // this should make a deep copy
+
+    if ( block->hasNoData() )
+    {
+      // turn all no-data values into NaN in the output array
+      float *floatData = reinterpret_cast<float *>( data.data() );
+      Q_ASSERT( data.count() % sizeof( float ) == 0 );
+      int count = data.count() / sizeof( float );
+      for ( int i = 0; i < count; ++i )
+      {
+        if ( block->isNoData( i ) )
+          floatData[i] = std::numeric_limits<float>::quiet_NaN();
+      }
+    }
+
     delete block;
   }
   return data;