QGIS API Documentation  3.8.0-Zanzibar (11aff65)
qgsdemterraintilegeometry_p.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsdemterraintilegeometry_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 #include <QMatrix4x4>
18 #include <Qt3DRender/qattribute.h>
19 #include <Qt3DRender/qbuffer.h>
20 #include <Qt3DRender/qbufferdatagenerator.h>
21 #include <limits>
22 #include <cmath>
23 #include "qgsraycastingutils_p.h"
24 
26 
27 using namespace Qt3DRender;
28 
29 
30 static QByteArray createPlaneVertexData( int res, float side, float vertScale, float skirtHeight, const QByteArray &heights )
31 {
32  Q_ASSERT( res >= 2 );
33  Q_ASSERT( heights.count() == res * res * static_cast<int>( sizeof( float ) ) );
34 
35  const float *zData = ( const float * ) heights.constData();
36  const float *zBits = zData;
37 
38  const int nVerts = ( res + 2 ) * ( res + 2 );
39 
40  // Populate a buffer with the interleaved per-vertex data with
41  // vec3 pos, vec2 texCoord, vec3 normal, vec4 tangent
42  const quint32 elementSize = 3 + 2 + 3;
43  const quint32 stride = elementSize * sizeof( float );
44  QByteArray bufferBytes;
45  bufferBytes.resize( stride * nVerts );
46  float *fptr = reinterpret_cast<float *>( bufferBytes.data() );
47 
48  float w = 1, h = 1;
49  QSize resolution( res, res );
50  const float x0 = -w / 2.0f;
51  const float z0 = -h / 2.0f;
52  const float dx = w / ( resolution.width() - 1 );
53  const float dz = h / ( resolution.height() - 1 );
54  const float du = 1.0 / ( resolution.width() - 1 );
55  const float dv = 1.0 / ( resolution.height() - 1 );
56 
57  // the height of vertices with no-data value... the value should not really matter
58  // as we do not create valid triangles that would use such vertices
59  const float noDataHeight = 0;
60 
61  const int iMax = resolution.width() - 1;
62  const int jMax = resolution.height() - 1;
63 
64  // Iterate over z
65  for ( int j = -1; j <= resolution.height(); ++j )
66  {
67  int jBound = qBound( 0, j, jMax );
68  const float z = z0 + static_cast<float>( jBound ) * dz;
69  const float v = static_cast<float>( jBound ) * dv;
70 
71  // Iterate over x
72  for ( int i = -1; i <= resolution.width(); ++i )
73  {
74  int iBound = qBound( 0, i, iMax );
75  const float x = x0 + static_cast<float>( iBound ) * dx;
76  const float u = static_cast<float>( iBound ) * du;
77 
78  float height;
79  if ( i == iBound && j == jBound )
80  height = *zBits++;
81  else
82  height = zData[ jBound * resolution.width() + iBound ] - skirtHeight;
83 
84  if ( std::isnan( height ) )
85  height = noDataHeight;
86 
87  // position
88  *fptr++ = x;
89  *fptr++ = height / side * vertScale;
90  *fptr++ = z;
91 
92  // texture coordinates
93  *fptr++ = u;
94  *fptr++ = v;
95 
96  // calculate normal coordinates
97 #define zAt( ii, jj ) zData[ jj * resolution.width() + ii ] * vertScale
98  float zi0 = zAt( qBound( 0, i - 1, iMax ), jBound );
99  float zi1 = zAt( qBound( 0, i + 1, iMax ), jBound );
100  float zj0 = zAt( iBound, qBound( 0, j - 1, jMax ) );
101  float zj1 = zAt( iBound, qBound( 0, j + 1, jMax ) );
102 
103  QVector3D n;
104  if ( std::isnan( zi0 ) || std::isnan( zi1 ) || std::isnan( zj0 ) || std::isnan( zj1 ) )
105  n = QVector3D( 0, 1, 0 );
106  else
107  {
108  float di, dj;
109  float zij = height * vertScale;
110 
111  if ( i == 0 )
112  di = 2 * ( zij - zi1 );
113  else if ( i == iMax )
114  di = 2 * ( zi0 - zij );
115  else
116  di = zi0 - zi1;
117 
118  if ( j == 0 )
119  dj = 2 * ( zij - zj1 );
120  else if ( j == jMax )
121  dj = 2 * ( zj0 - zij );
122  else
123  dj = zj0 - zj1;
124 
125  n = QVector3D( di, 2 * side / res, dj );
126  n.normalize();
127  }
128 
129  *fptr++ = n.x();
130  *fptr++ = n.y();
131  *fptr++ = n.z();
132  }
133  }
134 
135  return bufferBytes;
136 }
137 
138 inline int ijToHeightMapIndex( int i, int j, int resX, int resZ )
139 {
140  i = qBound( 1, i, resX ) - 1;
141  j = qBound( 1, j, resZ ) - 1;
142  return j * resX + i;
143 }
144 
145 
146 static bool hasNoData( int i, int j, const float *heightMap, int resX, int resZ )
147 {
148  return std::isnan( heightMap[ ijToHeightMapIndex( i, j, resX, resZ ) ] ) ||
149  std::isnan( heightMap[ ijToHeightMapIndex( i + 1, j, resX, resZ ) ] ) ||
150  std::isnan( heightMap[ ijToHeightMapIndex( i, j + 1, resX, resZ ) ] ) ||
151  std::isnan( heightMap[ ijToHeightMapIndex( i + 1, j + 1, resX, resZ ) ] );
152 }
153 
154 static QByteArray createPlaneIndexData( int res, const QByteArray &heightMap )
155 {
156  QSize resolution( res, res );
157  int numVerticesX = resolution.width() + 2;
158  int numVerticesZ = resolution.height() + 2;
159 
160  // Create the index data. 2 triangles per rectangular face
161  const int faces = 2 * ( numVerticesX - 1 ) * ( numVerticesZ - 1 );
162  const quint32 indices = 3 * faces;
163  Q_ASSERT( indices < std::numeric_limits<quint32>::max() );
164  QByteArray indexBytes;
165  indexBytes.resize( indices * sizeof( quint32 ) );
166  quint32 *indexPtr = reinterpret_cast<quint32 *>( indexBytes.data() );
167 
168  const float *heightMapFloat = reinterpret_cast<const float *>( heightMap.constData() );
169 
170  // Iterate over z
171  for ( int j = 0; j < numVerticesZ - 1; ++j )
172  {
173  const int rowStartIndex = j * numVerticesX;
174  const int nextRowStartIndex = ( j + 1 ) * numVerticesX;
175 
176  // Iterate over x
177  for ( int i = 0; i < numVerticesX - 1; ++i )
178  {
179  if ( hasNoData( i, j, heightMapFloat, res, res ) )
180  {
181  // at least one corner of the quad has no-data value
182  // so let's make two invalid triangles
183  *indexPtr++ = rowStartIndex + i;
184  *indexPtr++ = rowStartIndex + i;
185  *indexPtr++ = rowStartIndex + i;
186 
187  *indexPtr++ = rowStartIndex + i;
188  *indexPtr++ = rowStartIndex + i;
189  *indexPtr++ = rowStartIndex + i;
190  continue;
191  }
192 
193  // Split quad into two triangles
194  *indexPtr++ = rowStartIndex + i;
195  *indexPtr++ = nextRowStartIndex + i;
196  *indexPtr++ = rowStartIndex + i + 1;
197 
198  *indexPtr++ = nextRowStartIndex + i;
199  *indexPtr++ = nextRowStartIndex + i + 1;
200  *indexPtr++ = rowStartIndex + i + 1;
201  }
202  }
203 
204  return indexBytes;
205 }
206 
207 
208 
210 class PlaneVertexBufferFunctor : public QBufferDataGenerator
211 {
212  public:
213  explicit PlaneVertexBufferFunctor( int resolution, float side, float vertScale, float skirtHeight, const QByteArray &heightMap )
214  : mResolution( resolution )
215  , mSide( side )
216  , mVertScale( vertScale )
217  , mSkirtHeight( skirtHeight )
218  , mHeightMap( heightMap )
219  {}
220 
221  QByteArray operator()() final
222  {
223  return createPlaneVertexData( mResolution, mSide, mVertScale, mSkirtHeight, mHeightMap );
224  }
225 
226  bool operator ==( const QBufferDataGenerator &other ) const final
227  {
228  const PlaneVertexBufferFunctor *otherFunctor = functor_cast<PlaneVertexBufferFunctor>( &other );
229  if ( otherFunctor != nullptr )
230  return ( otherFunctor->mResolution == mResolution &&
231  otherFunctor->mSide == mSide &&
232  otherFunctor->mVertScale == mVertScale &&
233  otherFunctor->mSkirtHeight == mSkirtHeight &&
234  otherFunctor->mHeightMap == mHeightMap );
235  return false;
236  }
237 
238  QT3D_FUNCTOR( PlaneVertexBufferFunctor )
239 
240  private:
241  int mResolution;
242  float mSide;
243  float mVertScale;
244  float mSkirtHeight;
245  QByteArray mHeightMap;
246 };
247 
248 
250 class PlaneIndexBufferFunctor : public QBufferDataGenerator
251 {
252  public:
253  explicit PlaneIndexBufferFunctor( int resolution, const QByteArray &heightMap )
254  : mResolution( resolution )
255  , mHeightMap( heightMap )
256  {}
257 
258  QByteArray operator()() final
259  {
260  return createPlaneIndexData( mResolution, mHeightMap );
261  }
262 
263  bool operator ==( const QBufferDataGenerator &other ) const final
264  {
265  const PlaneIndexBufferFunctor *otherFunctor = functor_cast<PlaneIndexBufferFunctor>( &other );
266  if ( otherFunctor != nullptr )
267  return ( otherFunctor->mResolution == mResolution );
268  return false;
269  }
270 
271  QT3D_FUNCTOR( PlaneIndexBufferFunctor )
272 
273  private:
274  int mResolution;
275  QByteArray mHeightMap;
276 };
277 
278 
279 // ------------
280 
281 
282 DemTerrainTileGeometry::DemTerrainTileGeometry( int resolution, float side, float vertScale, float skirtHeight, const QByteArray &heightMap, DemTerrainTileGeometry::QNode *parent )
283  : QGeometry( parent )
284  , mResolution( resolution )
285  , mSide( side )
286  , mVertScale( vertScale )
287  , mSkirtHeight( skirtHeight )
288  , mHeightMap( heightMap )
289 {
290  init();
291 }
292 
293 static bool intersectionDemTriangles( const QByteArray &vertexBuf, const QByteArray &indexBuf, const QgsRayCastingUtils::Ray3D &r, const QMatrix4x4 &worldTransform, QVector3D &intPt )
294 {
295  // WARNING! this code is specific to how vertex buffers are built for DEM tiles,
296  // it is not usable for any mesh...
297 
298  const float *vertices = reinterpret_cast<const float *>( vertexBuf.constData() );
299  const uint *indices = reinterpret_cast<const uint *>( indexBuf.constData() );
300 #ifdef QGISDEBUG
301  int vertexCnt = vertexBuf.count() / sizeof( float );
302  Q_ASSERT( vertexCnt % 8 == 0 );
303 #endif
304  int indexCnt = indexBuf.count() / sizeof( uint );
305  Q_ASSERT( indexCnt % 3 == 0 );
306  int triangleCount = indexCnt / 3;
307 
308  QVector3D intersectionPt, minIntersectionPt;
309  float distance;
310  float minDistance = -1;
311 
312  for ( int i = 0; i < triangleCount; ++i )
313  {
314  int v0 = indices[i * 3], v1 = indices[i * 3 + 1], v2 = indices[i * 3 + 2];
315  QVector3D a( vertices[v0 * 8], vertices[v0 * 8 + 1], vertices[v0 * 8 + 2] );
316  QVector3D b( vertices[v1 * 8], vertices[v1 * 8 + 1], vertices[v1 * 8 + 2] );
317  QVector3D c( vertices[v2 * 8], vertices[v2 * 8 + 1], vertices[v2 * 8 + 2] );
318 
319  const QVector3D tA = worldTransform * a;
320  const QVector3D tB = worldTransform * b;
321  const QVector3D tC = worldTransform * c;
322 
323  QVector3D uvw;
324  float t = 0;
325  if ( QgsRayCastingUtils::rayTriangleIntersection( r, tA, tB, tC, uvw, t ) )
326  {
327  intersectionPt = r.point( t * r.distance() );
328  distance = r.projectedDistance( intersectionPt );
329 
330  // we only want the first intersection of the ray with the mesh (closest to the ray origin)
331  if ( minDistance == -1 || distance < minDistance )
332  {
333  minDistance = distance;
334  minIntersectionPt = intersectionPt;
335  }
336  }
337  }
338 
339  if ( minDistance != -1 )
340  {
341  intPt = minIntersectionPt;
342  return true;
343  }
344  else
345  return false;
346 }
347 
348 bool DemTerrainTileGeometry::rayIntersection( const QgsRayCastingUtils::Ray3D &ray, const QMatrix4x4 &worldTransform, QVector3D &intersectionPoint )
349 {
350  return intersectionDemTriangles( mVertexBuffer->data(), mIndexBuffer->data(), ray, worldTransform, intersectionPoint );
351 }
352 
353 void DemTerrainTileGeometry::init()
354 {
355  mPositionAttribute = new QAttribute( this );
356  mNormalAttribute = new QAttribute( this );
357  mTexCoordAttribute = new QAttribute( this );
358  mIndexAttribute = new QAttribute( this );
359  mVertexBuffer = new Qt3DRender::QBuffer( Qt3DRender::QBuffer::VertexBuffer, this );
360  mIndexBuffer = new Qt3DRender::QBuffer( Qt3DRender::QBuffer::IndexBuffer, this );
361 
362  int nVertsX = mResolution + 2;
363  int nVertsZ = mResolution + 2;
364  const int nVerts = nVertsX * nVertsZ;
365  const int stride = ( 3 + 2 + 3 ) * sizeof( float );
366  const int faces = 2 * ( nVertsX - 1 ) * ( nVertsZ - 1 );
367 
368  mPositionAttribute->setName( QAttribute::defaultPositionAttributeName() );
369  mPositionAttribute->setVertexBaseType( QAttribute::Float );
370  mPositionAttribute->setVertexSize( 3 );
371  mPositionAttribute->setAttributeType( QAttribute::VertexAttribute );
372  mPositionAttribute->setBuffer( mVertexBuffer );
373  mPositionAttribute->setByteStride( stride );
374  mPositionAttribute->setCount( nVerts );
375 
376  mTexCoordAttribute->setName( QAttribute::defaultTextureCoordinateAttributeName() );
377  mTexCoordAttribute->setVertexBaseType( QAttribute::Float );
378  mTexCoordAttribute->setVertexSize( 2 );
379  mTexCoordAttribute->setAttributeType( QAttribute::VertexAttribute );
380  mTexCoordAttribute->setBuffer( mVertexBuffer );
381  mTexCoordAttribute->setByteStride( stride );
382  mTexCoordAttribute->setByteOffset( 3 * sizeof( float ) );
383  mTexCoordAttribute->setCount( nVerts );
384 
385  mNormalAttribute->setName( QAttribute::defaultNormalAttributeName() );
386  mNormalAttribute->setVertexBaseType( QAttribute::Float );
387  mNormalAttribute->setVertexSize( 3 );
388  mNormalAttribute->setAttributeType( QAttribute::VertexAttribute );
389  mNormalAttribute->setBuffer( mVertexBuffer );
390  mNormalAttribute->setByteStride( stride );
391  mNormalAttribute->setByteOffset( 5 * sizeof( float ) );
392  mNormalAttribute->setCount( nVerts );
393 
394  mIndexAttribute->setAttributeType( QAttribute::IndexAttribute );
395  mIndexAttribute->setVertexBaseType( QAttribute::UnsignedInt );
396  mIndexAttribute->setBuffer( mIndexBuffer );
397 
398  // Each primitive has 3 vertives
399  mIndexAttribute->setCount( faces * 3 );
400 
401  // switched to setting data instead of just setting data generators because we also need the buffers
402  // available for ray-mesh intersections and we can't access the private copy of data in Qt (if there is any)
403  mVertexBuffer->setData( PlaneVertexBufferFunctor( mResolution, mSide, mVertScale, mSkirtHeight, mHeightMap )() );
404  mIndexBuffer->setData( PlaneIndexBufferFunctor( mResolution, mHeightMap )() );
405 
406  addAttribute( mPositionAttribute );
407  addAttribute( mTexCoordAttribute );
408  addAttribute( mNormalAttribute );
409  addAttribute( mIndexAttribute );
410 }
411 
bool operator==(const QgsFeatureIterator &fi1, const QgsFeatureIterator &fi2)
As part of the API refactoring and improvements which landed in the Processing API was substantially reworked from the x version This was done in order to allow much of the underlying Processing framework to be ported into c