QGIS API Documentation  3.26.3-Buenos Aires (65e4edfdad)
quantizedmeshgeometry.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  quantizedmeshgeometry.cpp
3  ---------------------
4  begin : 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 #include "quantizedmeshgeometry.h"
16 
17 #include <zlib.h>
18 
19 #include <QFile>
20 #include <QByteArray>
21 
22 
23 // gzip decompression snipped from https://stackoverflow.com/questions/2690328/qt-quncompress-gzip-data
24 
25 #define GZIP_WINDOWS_BIT 15 + 16
26 #define GZIP_CHUNK_SIZE 32 * 1024
27 
34 bool gzipDecompress( QByteArray input, QByteArray &output )
35 {
36  // Prepare output
37  output.clear();
38 
39  // Is there something to do?
40  if ( input.isEmpty() )
41  return true;
42 
43  // Prepare inflater status
44  z_stream strm;
45  strm.zalloc = Z_NULL;
46  strm.zfree = Z_NULL;
47  strm.opaque = Z_NULL;
48  strm.avail_in = 0;
49  strm.next_in = Z_NULL;
50 
51  // Initialize inflater
52  int ret = inflateInit2( &strm, GZIP_WINDOWS_BIT );
53 
54  if ( ret != Z_OK )
55  return false;
56 
57  // Extract pointer to input data
58  const char *input_data = input.constData();
59  int input_data_left = input.length();
60 
61  // Decompress data until available
62  do
63  {
64  // Determine current chunk size
65  int chunk_size = std::min( GZIP_CHUNK_SIZE, input_data_left );
66 
67  // Check for termination
68  if ( chunk_size <= 0 )
69  break;
70 
71  // Set inflater references
72  strm.next_in = ( unsigned char * )input_data;
73  strm.avail_in = chunk_size;
74 
75  // Update interval variables
76  input_data += chunk_size;
77  input_data_left -= chunk_size;
78 
79  // Inflate chunk and cumulate output
80  do
81  {
82 
83  // Declare vars
84  char out[GZIP_CHUNK_SIZE];
85 
86  // Set inflater references
87  strm.next_out = ( unsigned char * )out;
88  strm.avail_out = GZIP_CHUNK_SIZE;
89 
90  // Try to inflate chunk
91  ret = inflate( &strm, Z_NO_FLUSH );
92 
93  switch ( ret )
94  {
95  case Z_NEED_DICT:
96  ret = Z_DATA_ERROR;
97  case Z_DATA_ERROR:
98  case Z_MEM_ERROR:
99  case Z_STREAM_ERROR:
100  // Clean-up
101  inflateEnd( &strm );
102 
103  // Return
104  return ( false );
105  }
106 
107  // Determine decompressed size
108  int have = ( GZIP_CHUNK_SIZE - strm.avail_out );
109 
110  // Cumulate result
111  if ( have > 0 )
112  output.append( ( char * )out, have );
113 
114  }
115  while ( strm.avail_out == 0 );
116 
117  }
118  while ( ret != Z_STREAM_END );
119 
120  // Clean-up
121  inflateEnd( &strm );
122 
123  // Return
124  return ( ret == Z_STREAM_END );
125 }
126 
127 
128 const char *read_zigzag_encoded_int16_array( const char *dataPtr, int count, qint16 *out )
129 {
130  for ( int i = 0; i < count; ++i )
131  {
132  quint16 encoded = *( quint16 * ) dataPtr;
133  dataPtr += 2;
134  qint16 decoded = ( encoded >> 1 ) ^ ( -( encoded & 1 ) );
135  *out++ = decoded;
136  }
137  return dataPtr;
138 }
139 
140 
141 static QString _tileFilename( int tx, int ty, int tz )
142 {
143  return QString( "/tmp/terrain-%1-%2-%3" ).arg( tz ).arg( tx ).arg( ty );
144 }
145 
146 QuantizedMeshTile *QuantizedMeshGeometry::readTile( int tx, int ty, int tz, const QgsRectangle &extent )
147 {
148  QString filename = _tileFilename( tx, ty, tz );
149  QFile f( filename );
150  if ( !f.open( QIODevice::ReadOnly ) )
151  return nullptr;
152 
153  QByteArray data;
154  if ( !gzipDecompress( f.readAll(), data ) )
155  return nullptr;
156 
157  if ( data.isEmpty() )
158  return nullptr;
159 
161  t->extent = extent;
162 
163  const char *dataPtr = data.constData();
164  memcpy( &t->header, dataPtr, sizeof( QuantizedMeshHeader ) );
165  dataPtr += sizeof( QuantizedMeshHeader );
166 
167  // vertex data - immediately after header
168  // with zig-zag encoding
169 
170  //struct VertexData
171  //{
172  // unsigned int vertexCount;
173  // unsigned short u[vertexCount];
174  // unsigned short v[vertexCount];
175  // unsigned short height[vertexCount];
176  //};
177 
178  quint32 vertexCount = *( quint32 * ) dataPtr;
179  dataPtr += 4;
180  t->uvh.resize( 3 * vertexCount );
181  qint16 *vptr = t->uvh.data();
182  dataPtr = read_zigzag_encoded_int16_array( dataPtr, vertexCount * 3, vptr );
183  // the individual values are just deltas of previous values!
184  qint16 u = 0, v = 0, h = 0;
185  for ( uint i = 0; i < vertexCount; ++i )
186  {
187  qint16 du = vptr[i], dv = vptr[vertexCount + i], dh = vptr[vertexCount * 2 + i];
188  u += du;
189  v += dv;
190  h += dh;
191  vptr[i] = u;
192  vptr[vertexCount + i] = v;
193  vptr[vertexCount * 2 + i] = h;
194  }
195 
196  Q_ASSERT( vertexCount < 65537 ); // supporting currently only 2-byte vertex indices
197 
198  // index data - if less than 65537 vertices (otherwise indices would be 4-byte)
199  // with "high watermark" encoding
200 
201  //struct IndexData16
202  //{
203  // unsigned int triangleCount;
204  // unsigned short indices[triangleCount * 3];
205  //}
206 
207  quint32 triangleCount = *( quint32 * ) dataPtr;
208  dataPtr += 4;
209  t->indices.resize( 3 * triangleCount );
210  quint16 *indicesPtr = t->indices.data();
211  quint16 *srcIdxPtr = ( quint16 * )dataPtr;
212  int highest = 0;
213  for ( uint i = 0; i < triangleCount * 3; ++i )
214  {
215  quint16 code = *srcIdxPtr++;
216  *indicesPtr++ = highest - code;
217  if ( code == 0 )
218  ++highest;
219  }
220 
221  // TODO: edge indices
222 
223  // TODO: extensions
224 
225  //qDebug() << hdr.CenterX << " " << hdr.CenterY << " " << hdr.CenterZ;
226  //qDebug() << "VC " << vertexCount;
227  //qDebug() << "TC " << triangleCount;
228 
229  return t;
230 }
231 
232 
233 #include <QGuiApplication>
234 #include <QNetworkRequest>
235 #include <QNetworkReply>
236 #include "qgsnetworkaccessmanager.h"
237 
238 void QuantizedMeshGeometry::downloadTileIfMissing( int tx, int ty, int tz )
239 {
240  // i am not proud of this bit of code... quick&dirty!
241 
242  QString tileFilename = _tileFilename( tx, ty, tz );
243  if ( !QFile::exists( tileFilename ) )
244  {
245  qDebug() << "downloading tile " << tx << " " << ty << " " << tz;
246  bool downloaded = false;
247  QString url = QString( "http://assets.agi.com/stk-terrain/tilesets/world/tiles/%1/%2/%3.terrain" ).arg( tz ).arg( tx ).arg( ty );
248  QNetworkRequest request( url );
249  request.setRawHeader( QByteArray( "Accept-Encoding" ), QByteArray( "gzip" ) );
250  request.setRawHeader( QByteArray( "Accept" ), QByteArray( "application/vnd.quantized-mesh,application/octet-stream;q=0.9" ) );
251  request.setRawHeader( QByteArray( "User-Agent" ), QByteArray( "Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) Ubuntu Chromium/58.0.3029.110 Chrome/58.0.3029.110 Safari/537.36" ) );
252  QNetworkReply *reply = QgsNetworkAccessManager::instance()->get( request );
253  connect( reply, &QNetworkReply::finished, [reply, tileFilename, &downloaded]
254  {
255  QFile fOut( tileFilename );
256  fOut.open( QIODevice::WriteOnly );
257  fOut.write( reply->readAll() );
258  fOut.close();
259  reply->deleteLater();
260  downloaded = true;
261  } );
262 
263  while ( !downloaded )
264  {
265  qApp->processEvents();
266  }
267  }
268 }
269 
270 // --------------
271 
272 #include "qgscoordinatetransform.h"
273 #include "qgsmaptopixel.h"
274 #include "map3d.h"
275 
276 QuantizedMeshGeometry::QuantizedMeshGeometry( QuantizedMeshTile *t, const Map3D &map, const QgsMapToPixel &mapToPixel, const QgsCoordinateTransform &terrainToMap, QNode *parent )
277  : QGeometry( parent )
278 {
279  int vertexCount = t->uvh.count() / 3;
280  int indexCount = t->indices.count();
281 
282  int vertexEntrySize = sizeof( float ) * ( 3 + 2 );
283 
284  double xMinWgs = t->extent.xMinimum();
285  double yMinWgs = t->extent.yMinimum();
286  double widthWgs = t->extent.width();
287  double heightWgs = t->extent.height();
288  QgsPointXY ptMinProjected = QgsPointXY( map.originX, map.originY );
289 
290  QByteArray vb;
291  const qint16 *uvh = t->uvh.constData();
292  vb.resize( vertexCount * vertexEntrySize );
293  float *vbptr = ( float * ) vb.data();
294  for ( int i = 0; i < vertexCount; ++i )
295  {
296  qint16 u = uvh[i], v = uvh[vertexCount + i], h = uvh[vertexCount * 2 + i];
297  float uNorm = u / 32767.f; // 0...1
298  float vNorm = v / 32767.f; // 0...1
299  float hNorm = h / 32767.f; // 0...1
300  float xWgs = xMinWgs + widthWgs * uNorm;
301  float yWgs = yMinWgs + heightWgs * vNorm;
302  float hWgs = t->header.MinimumHeight + hNorm * ( t->header.MaximumHeight - t->header.MinimumHeight );
303 
304  QgsPointXY ptProjected = terrainToMap.transform( xWgs, yWgs );
305  QgsPointXY ptFinal( ptProjected.x() - ptMinProjected.x(), ptProjected.y() - ptMinProjected.y() );
306 
307  // our plane is (x,-z) with y growing towards camera
308  *vbptr++ = ptFinal.x();
309  *vbptr++ = hWgs;
310  *vbptr++ = -ptFinal.y();
311 
312  QgsPointXY uv = mapToPixel.transform( ptProjected ) / map.tileTextureSize;
313  // texture coords
314  *vbptr++ = uv.x();
315  *vbptr++ = uv.y();
316  }
317 
318  QByteArray ib;
319  ib.resize( indexCount * 2 );
320  memcpy( ib.data(), t->indices.constData(), ib.count() );
321  /*
322  quint16* ibptr = (quint16*)ib.data();
323  const quint16* srcptr = t->indices.constData();
324  // reverse order of indices to triangles
325  for (int i = 0; i < indexCount/3; ++i)
326  {
327  ibptr[i*3] = srcptr[i*3+2];
328  ibptr[i*3+1] = srcptr[i*3+1];
329  ibptr[i*3+2] = srcptr[i*3];
330  }
331  */
332  m_vertexBuffer = new Qt3DRender::QBuffer( this );
333  m_indexBuffer = new Qt3DRender::QBuffer( this );
334 
335  m_vertexBuffer->setData( vb );
336  m_indexBuffer->setData( ib );
337 
338  m_positionAttribute = new Qt3DRender::QAttribute( this );
339  m_positionAttribute->setName( Qt3DRender::QAttribute::defaultPositionAttributeName() );
340  m_positionAttribute->setVertexBaseType( Qt3DRender::QAttribute::Float );
341  m_positionAttribute->setVertexSize( 3 );
342  m_positionAttribute->setAttributeType( Qt3DRender::QAttribute::VertexAttribute );
343  m_positionAttribute->setBuffer( m_vertexBuffer );
344  m_positionAttribute->setByteStride( vertexEntrySize );
345  m_positionAttribute->setCount( vertexCount );
346 
347  m_texCoordAttribute = new Qt3DRender::QAttribute( this );
348  m_texCoordAttribute->setName( Qt3DRender::QAttribute::defaultTextureCoordinateAttributeName() );
349  m_texCoordAttribute->setVertexBaseType( Qt3DRender::QAttribute::Float );
350  m_texCoordAttribute->setVertexSize( 2 );
351  m_texCoordAttribute->setAttributeType( Qt3DRender::QAttribute::VertexAttribute );
352  m_texCoordAttribute->setBuffer( m_vertexBuffer );
353  m_texCoordAttribute->setByteStride( vertexEntrySize );
354  m_texCoordAttribute->setByteOffset( 3 * sizeof( float ) );
355  m_texCoordAttribute->setCount( vertexCount );
356 
357  m_indexAttribute = new Qt3DRender::QAttribute( this );
358  m_indexAttribute->setAttributeType( Qt3DRender::QAttribute::IndexAttribute );
359  m_indexAttribute->setVertexBaseType( Qt3DRender::QAttribute::UnsignedShort );
360  m_indexAttribute->setBuffer( m_indexBuffer );
361  m_indexAttribute->setCount( indexCount );
362 
363  addAttribute( m_positionAttribute );
364  addAttribute( m_texCoordAttribute );
365  addAttribute( m_indexAttribute );
366 }
QgsRectangle::height
double height() const SIP_HOLDGIL
Returns the height of the rectangle.
Definition: qgsrectangle.h:230
QuantizedMeshGeometry::downloadTileIfMissing
static void downloadTileIfMissing(int tx, int ty, int tz)
Downloads a tile to to a file in the local disk cache.
Definition: quantizedmeshgeometry.cpp:238
QgsPointXY::y
double y
Definition: qgspointxy.h:63
quantizedmeshgeometry.h
qgsmaptopixel.h
QgsRectangle::yMinimum
double yMinimum() const SIP_HOLDGIL
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:198
QuantizedMeshGeometry::QuantizedMeshGeometry
QuantizedMeshGeometry(QuantizedMeshTile *t, const Map3D &map, const QgsMapToPixel &mapToPixel, const QgsCoordinateTransform &terrainToMap, QNode *parent=nullptr)
Constructs geometry based on the loaded tile data.
Definition: quantizedmeshgeometry.cpp:276
read_zigzag_encoded_int16_array
const char * read_zigzag_encoded_int16_array(const char *dataPtr, int count, qint16 *out)
Definition: quantizedmeshgeometry.cpp:128
QgsRectangle
A rectangle specified with double values.
Definition: qgsrectangle.h:41
QuantizedMeshTile::uvh
QVector< qint16 > uvh
Definition: quantizedmeshgeometry.h:62
GZIP_WINDOWS_BIT
#define GZIP_WINDOWS_BIT
Definition: quantizedmeshgeometry.cpp:25
qgsnetworkaccessmanager.h
QuantizedMeshTile::extent
QgsRectangle extent
Definition: quantizedmeshgeometry.h:60
qgscoordinatetransform.h
QgsRectangle::xMinimum
double xMinimum() const SIP_HOLDGIL
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:188
QgsNetworkAccessManager::instance
static QgsNetworkAccessManager * instance(Qt::ConnectionType connectionType=Qt::BlockingQueuedConnection)
Returns a pointer to the active QgsNetworkAccessManager for the current thread.
Definition: qgsnetworkaccessmanager.cpp:202
QgsMapToPixel::transform
QgsPointXY transform(const QgsPointXY &p) const
Transforms a point p from map (world) coordinates to device coordinates.
Definition: qgsmaptopixel.h:90
QuantizedMeshHeader
Definition: quantizedmeshgeometry.h:26
QuantizedMeshTile::indices
QVector< quint16 > indices
Definition: quantizedmeshgeometry.h:63
QuantizedMeshHeader::MaximumHeight
float MaximumHeight
Definition: quantizedmeshgeometry.h:52
QgsPointXY
A class to represent a 2D point.
Definition: qgspointxy.h:58
QuantizedMeshHeader::MinimumHeight
float MinimumHeight
Definition: quantizedmeshgeometry.h:51
QuantizedMeshTile::header
QuantizedMeshHeader header
Definition: quantizedmeshgeometry.h:61
QgsMapToPixel
Perform transforms between map coordinates and device coordinates.
Definition: qgsmaptopixel.h:38
QgsRectangle::width
double width() const SIP_HOLDGIL
Returns the width of the rectangle.
Definition: qgsrectangle.h:223
QgsPointXY::x
double x
Definition: qgspointxy.h:62
QuantizedMeshTile
Definition: quantizedmeshgeometry.h:58
GZIP_CHUNK_SIZE
#define GZIP_CHUNK_SIZE
Definition: quantizedmeshgeometry.cpp:26
gzipDecompress
bool gzipDecompress(QByteArray input, QByteArray &output)
Decompresses the given buffer using the standard GZIP algorithm.
Definition: quantizedmeshgeometry.cpp:34
QgsCoordinateTransform::transform
QgsPointXY transform(const QgsPointXY &point, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward) const SIP_THROW(QgsCsException)
Transform the point from the source CRS to the destination CRS.
Definition: qgscoordinatetransform.cpp:272
QgsCoordinateTransform
Class for doing transforms between two map coordinate systems.
Definition: qgscoordinatetransform.h:57
QuantizedMeshGeometry::readTile
static QuantizedMeshTile * readTile(int tx, int ty, int tz, const QgsRectangle &extent)
Reads a tile from a file in the local disk cache.
Definition: quantizedmeshgeometry.cpp:146