QGIS API Documentation  3.25.0-Master (dec16ba68b)
qgstessellator.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgstessellator.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 
16 #include "qgstessellator.h"
17 
18 #include "qgscurve.h"
19 #include "qgsgeometry.h"
20 #include "qgsmessagelog.h"
21 #include "qgsmultipolygon.h"
22 #include "qgspoint.h"
23 #include "qgspolygon.h"
24 #include "qgstriangle.h"
25 #include "qgis_sip.h"
26 #include "qgsgeometryengine.h"
27 
28 #include "poly2tri.h"
29 
30 #include <QtDebug>
31 #include <QMatrix4x4>
32 #include <QVector3D>
33 #include <QtMath>
34 #include <algorithm>
35 #include <unordered_set>
36 
37 static std::pair<float, float> rotateCoords( float x, float y, float origin_x, float origin_y, float r )
38 {
39  r = qDegreesToRadians( r );
40  float x0 = x - origin_x, y0 = y - origin_y;
41  // p0 = x0 + i * y0
42  // rot = cos(r) + i * sin(r)
43  // p0 * rot = x0 * cos(r) - y0 * sin(r) + i * [ x0 * sin(r) + y0 * cos(r) ]
44  const float x1 = origin_x + x0 * qCos( r ) - y0 * qSin( r );
45  const float y1 = origin_y + x0 * qSin( r ) + y0 * qCos( r );
46  return std::make_pair( x1, y1 );
47 }
48 
49 static void make_quad( float x0, float y0, float z0, float x1, float y1, float z1, float height, QVector<float> &data, bool addNormals, bool addTextureCoords, float textureRotation )
50 {
51  const float dx = x1 - x0;
52  const float dy = -( y1 - y0 );
53 
54  // perpendicular vector in plane to [x,y] is [-y,x]
55  QVector3D vn( -dy, 0, dx );
56  vn = -vn;
57  vn.normalize();
58 
59  float u0, v0;
60  float u1, v1;
61  float u2, v2;
62  float u3, v3;
63 
64  QVector<double> textureCoordinates;
65  textureCoordinates.reserve( 12 );
66  // select which side of the coordinates to use (x, z or y, z) depending on which side is smaller
67  if ( fabsf( dy ) <= fabsf( dx ) )
68  {
69  // consider x and z as the texture coordinates
70  u0 = x0;
71  v0 = z0 + height;
72 
73  u1 = x1;
74  v1 = z1 + height;
75 
76  u2 = x0;
77  v2 = z0;
78 
79  u3 = x1;
80  v3 = z1;
81  }
82  else
83  {
84  // consider y and z as the texture coowallsTextureRotationrdinates
85  u0 = -y0;
86  v0 = z0 + height;
87 
88  u1 = -y1;
89  v1 = z1 + height;
90 
91  u2 = -y0;
92  v2 = z0;
93 
94  u3 = -y1;
95  v3 = z1;
96  }
97 
98  textureCoordinates.push_back( u0 );
99  textureCoordinates.push_back( v0 );
100 
101  textureCoordinates.push_back( u1 );
102  textureCoordinates.push_back( v1 );
103 
104  textureCoordinates.push_back( u2 );
105  textureCoordinates.push_back( v2 );
106 
107  textureCoordinates.push_back( u2 );
108  textureCoordinates.push_back( v2 );
109 
110  textureCoordinates.push_back( u1 );
111  textureCoordinates.push_back( v1 );
112 
113  textureCoordinates.push_back( u3 );
114  textureCoordinates.push_back( v3 );
115 
116  for ( int i = 0; i < textureCoordinates.size(); i += 2 )
117  {
118  const std::pair<float, float> rotated = rotateCoords( textureCoordinates[i], textureCoordinates[i + 1], 0, 0, textureRotation );
119  textureCoordinates[i] = rotated.first;
120  textureCoordinates[i + 1] = rotated.second;
121  }
122 
123  // triangle 1
124  // vertice 1
125  data << x0 << z0 + height << -y0;
126  if ( addNormals )
127  data << vn.x() << vn.y() << vn.z();
128  if ( addTextureCoords )
129  data << textureCoordinates[0] << textureCoordinates[1];
130  // vertice 2
131  data << x1 << z1 + height << -y1;
132  if ( addNormals )
133  data << vn.x() << vn.y() << vn.z();
134  if ( addTextureCoords )
135  data << textureCoordinates[2] << textureCoordinates[3];
136  // verice 3
137  data << x0 << z0 << -y0;
138  if ( addNormals )
139  data << vn.x() << vn.y() << vn.z();
140  if ( addTextureCoords )
141  data << textureCoordinates[4] << textureCoordinates[5];
142 
143  // triangle 2
144  // vertice 1
145  data << x0 << z0 << -y0;
146  if ( addNormals )
147  data << vn.x() << vn.y() << vn.z();
148  if ( addTextureCoords )
149  data << textureCoordinates[6] << textureCoordinates[7];
150  // vertice 2
151  data << x1 << z1 + height << -y1;
152  if ( addNormals )
153  data << vn.x() << vn.y() << vn.z();
154  if ( addTextureCoords )
155  data << textureCoordinates[8] << textureCoordinates[9];
156  // vertice 3
157  data << x1 << z1 << -y1;
158  if ( addNormals )
159  data << vn.x() << vn.y() << vn.z();
160  if ( addTextureCoords )
161  data << textureCoordinates[10] << textureCoordinates[11];
162 }
163 
164 
165 QgsTessellator::QgsTessellator( double originX, double originY, bool addNormals, bool invertNormals, bool addBackFaces, bool noZ,
166  bool addTextureCoords, int facade, float textureRotation )
167  : mOriginX( originX )
168  , mOriginY( originY )
169  , mAddNormals( addNormals )
170  , mInvertNormals( invertNormals )
171  , mAddBackFaces( addBackFaces )
172  , mAddTextureCoords( addTextureCoords )
173  , mNoZ( noZ )
174  , mTessellatedFacade( facade )
175  , mTextureRotation( textureRotation )
176 {
177  init();
178 }
179 
180 QgsTessellator::QgsTessellator( const QgsRectangle &bounds, bool addNormals, bool invertNormals, bool addBackFaces, bool noZ,
181  bool addTextureCoords, int facade, float textureRotation )
182  : mBounds( bounds )
183  , mOriginX( mBounds.xMinimum() )
184  , mOriginY( mBounds.yMinimum() )
185  , mAddNormals( addNormals )
186  , mInvertNormals( invertNormals )
187  , mAddBackFaces( addBackFaces )
188  , mAddTextureCoords( addTextureCoords )
189  , mNoZ( noZ )
190  , mTessellatedFacade( facade )
191  , mTextureRotation( textureRotation )
192 {
193  init();
194 }
195 
196 void QgsTessellator::init()
197 {
198  mStride = 3 * sizeof( float );
199  if ( mAddNormals )
200  mStride += 3 * sizeof( float );
201  if ( mAddTextureCoords )
202  mStride += 2 * sizeof( float );
203 }
204 
205 static bool _isRingCounterClockWise( const QgsCurve &ring )
206 {
207  double a = 0;
208  const int count = ring.numPoints();
209  Qgis::VertexType vt;
210  QgsPoint pt, ptPrev;
211  ring.pointAt( 0, ptPrev, vt );
212  for ( int i = 1; i < count + 1; ++i )
213  {
214  ring.pointAt( i % count, pt, vt );
215  a += ptPrev.x() * pt.y() - ptPrev.y() * pt.x();
216  ptPrev = pt;
217  }
218  return a > 0; // clockwise if a is negative
219 }
220 
221 static void _makeWalls( const QgsLineString &ring, bool ccw, float extrusionHeight, QVector<float> &data,
222  bool addNormals, bool addTextureCoords, double originX, double originY, float textureRotation )
223 {
224  // we need to find out orientation of the ring so that the triangles we generate
225  // face the right direction
226  // (for exterior we want clockwise order, for holes we want counter-clockwise order)
227  const bool is_counter_clockwise = _isRingCounterClockWise( ring );
228 
229  QgsPoint pt;
230  QgsPoint ptPrev = ring.pointN( is_counter_clockwise == ccw ? 0 : ring.numPoints() - 1 );
231  for ( int i = 1; i < ring.numPoints(); ++i )
232  {
233  pt = ring.pointN( is_counter_clockwise == ccw ? i : ring.numPoints() - i - 1 );
234  float x0 = ptPrev.x() - originX, y0 = ptPrev.y() - originY;
235  float x1 = pt.x() - originX, y1 = pt.y() - originY;
236  const float z0 = std::isnan( ptPrev.z() ) ? 0 : ptPrev.z();
237  const float z1 = std::isnan( pt.z() ) ? 0 : pt.z();
238 
239  // make a quad
240  make_quad( x0, y0, z0, x1, y1, z1, extrusionHeight, data, addNormals, addTextureCoords, textureRotation );
241  ptPrev = pt;
242  }
243 }
244 
245 static QVector3D _calculateNormal( const QgsLineString *curve, double originX, double originY, bool invertNormal )
246 {
247  // if it is just plain 2D curve there is no need to calculate anything
248  // because it will be a flat horizontally oriented patch
249  if ( !QgsWkbTypes::hasZ( curve->wkbType() ) || curve->isEmpty() )
250  return QVector3D( 0, 0, 1 );
251 
252  // often we have 3D coordinates, but Z is the same for all vertices
253  // so in order to save calculation and avoid possible issues with order of vertices
254  // (the calculation below may decide that a polygon faces downwards)
255  bool sameZ = true;
256  QgsPoint pt1 = curve->pointN( 0 );
257  QgsPoint pt2;
258  for ( int i = 1; i < curve->numPoints(); i++ )
259  {
260  pt2 = curve->pointN( i );
261  if ( pt1.z() != pt2.z() )
262  {
263  sameZ = false;
264  break;
265  }
266  }
267  if ( sameZ )
268  return QVector3D( 0, 0, 1 );
269 
270  // Calculate the polygon's normal vector, based on Newell's method
271  // https://www.khronos.org/opengl/wiki/Calculating_a_Surface_Normal
272  //
273  // Order of vertices is important here as it determines the front/back face of the polygon
274 
275  double nx = 0, ny = 0, nz = 0;
276  pt1 = curve->pointN( 0 );
277 
278  // shift points by the tessellator's origin - this does not affect normal calculation and it may save us from losing some precision
279  pt1.setX( pt1.x() - originX );
280  pt1.setY( pt1.y() - originY );
281  for ( int i = 1; i < curve->numPoints(); i++ )
282  {
283  pt2 = curve->pointN( i );
284  pt2.setX( pt2.x() - originX );
285  pt2.setY( pt2.y() - originY );
286 
287  if ( std::isnan( pt1.z() ) || std::isnan( pt2.z() ) )
288  continue;
289 
290  nx += ( pt1.y() - pt2.y() ) * ( pt1.z() + pt2.z() );
291  ny += ( pt1.z() - pt2.z() ) * ( pt1.x() + pt2.x() );
292  nz += ( pt1.x() - pt2.x() ) * ( pt1.y() + pt2.y() );
293 
294  pt1 = pt2;
295  }
296 
297  QVector3D normal( nx, ny, nz );
298  if ( invertNormal )
299  normal = -normal;
300  normal.normalize();
301  return normal;
302 }
303 
304 
305 static void _normalVectorToXYVectors( const QVector3D &pNormal, QVector3D &pXVector, QVector3D &pYVector )
306 {
307  // Here we define the two perpendicular vectors that define the local
308  // 2D space on the plane. They will act as axis for which we will
309  // calculate the projection coordinates of a 3D point to the plane.
310  if ( pNormal.z() > 0.001 || pNormal.z() < -0.001 )
311  {
312  pXVector = QVector3D( 1, 0, -pNormal.x() / pNormal.z() );
313  }
314  else if ( pNormal.y() > 0.001 || pNormal.y() < -0.001 )
315  {
316  pXVector = QVector3D( 1, -pNormal.x() / pNormal.y(), 0 );
317  }
318  else
319  {
320  pXVector = QVector3D( -pNormal.y() / pNormal.x(), 1, 0 );
321  }
322  pXVector.normalize();
323  pYVector = QVector3D::normal( pNormal, pXVector );
324 }
325 
327 {
328  std::size_t operator()( const std::pair<float, float> pair ) const
329  {
330  const std::size_t h1 = std::hash<float>()( pair.first );
331  const std::size_t h2 = std::hash<float>()( pair.second );
332 
333  return h1 ^ h2;
334  }
335 };
336 
337 static void _ringToPoly2tri( const QgsLineString *ring, std::vector<p2t::Point *> &polyline, QHash<p2t::Point *, float> *zHash )
338 {
339  const int pCount = ring->numPoints();
340 
341  polyline.reserve( pCount );
342 
343  const double *srcXData = ring->xData();
344  const double *srcYData = ring->yData();
345  const double *srcZData = ring->zData();
346  std::unordered_set<std::pair<float, float>, float_pair_hash> foundPoints;
347 
348  for ( int i = 0; i < pCount - 1; ++i )
349  {
350  const float x = *srcXData++;
351  const float y = *srcYData++;
352 
353  const auto res = foundPoints.insert( std::make_pair( x, y ) );
354  if ( !res.second )
355  {
356  // already used this point, don't add a second time
357  continue;
358  }
359 
360  p2t::Point *pt2 = new p2t::Point( x, y );
361  polyline.push_back( pt2 );
362  if ( zHash )
363  {
364  ( *zHash )[pt2] = *srcZData++;
365  }
366  }
367 }
368 
369 
370 inline double _round_coord( double x )
371 {
372  const double exp = 1e10; // round to 10 decimal digits
373  return round( x * exp ) / exp;
374 }
375 
376 
377 static QgsCurve *_transform_ring_to_new_base( const QgsLineString &curve, const QgsPoint &pt0, const QMatrix4x4 *toNewBase, const float scale )
378 {
379  const int count = curve.numPoints();
380  QVector<double> x;
381  QVector<double> y;
382  QVector<double> z;
383  x.resize( count );
384  y.resize( count );
385  z.resize( count );
386  double *xData = x.data();
387  double *yData = y.data();
388  double *zData = z.data();
389 
390  const double *srcXData = curve.xData();
391  const double *srcYData = curve.yData();
392  const double *srcZData = curve.is3D() ? curve.zData() : nullptr;
393 
394  for ( int i = 0; i < count; ++i )
395  {
396  QVector4D v( *srcXData++ - pt0.x(),
397  *srcYData++ - pt0.y(),
398  srcZData ? *srcZData++ - pt0.z() : 0,
399  0 );
400  if ( toNewBase )
401  v = toNewBase->map( v );
402 
403  // scale coordinates
404  v.setX( v.x() * scale );
405  v.setY( v.y() * scale );
406 
407  // we also round coordinates before passing them to poly2tri triangulation in order to fix possible numerical
408  // stability issues. We had crashes with nearly collinear points where one of the points was off by a tiny bit (e.g. by 1e-20).
409  // See TestQgsTessellator::testIssue17745().
410  //
411  // A hint for a similar issue: https://github.com/greenm01/poly2tri/issues/99
412  //
413  // The collinear tests uses epsilon 1e-12. Seems rounding to 12 places you still
414  // can get problems with this test when points are pretty much on a straight line.
415  // I suggest you round to 10 decimals for stability and you can live with that
416  // precision.
417  *xData++ = _round_coord( v.x() );
418  *yData++ = _round_coord( v.y() );
419  *zData++ = _round_coord( v.z() );
420  }
421  return new QgsLineString( x, y, z );
422 }
423 
424 
425 static QgsPolygon *_transform_polygon_to_new_base( const QgsPolygon &polygon, const QgsPoint &pt0, const QMatrix4x4 *toNewBase, const float scale )
426 {
427  QgsPolygon *p = new QgsPolygon;
428  p->setExteriorRing( _transform_ring_to_new_base( *qgsgeometry_cast< const QgsLineString * >( polygon.exteriorRing() ), pt0, toNewBase, scale ) );
429  for ( int i = 0; i < polygon.numInteriorRings(); ++i )
430  p->addInteriorRing( _transform_ring_to_new_base( *qgsgeometry_cast< const QgsLineString * >( polygon.interiorRing( i ) ), pt0, toNewBase, scale ) );
431  return p;
432 }
433 
434 
436 {
437  double min_d = 1e20;
438 
439  std::vector< const QgsLineString * > rings;
440  rings.reserve( 1 + polygon.numInteriorRings() );
441  rings.emplace_back( qgsgeometry_cast< const QgsLineString * >( polygon.exteriorRing() ) );
442  for ( int i = 0; i < polygon.numInteriorRings(); ++i )
443  rings.emplace_back( qgsgeometry_cast< const QgsLineString * >( polygon.interiorRing( i ) ) );
444 
445  for ( const QgsLineString *ring : rings )
446  {
447  const int numPoints = ring->numPoints();
448  if ( numPoints <= 1 )
449  continue;
450 
451  const double *srcXData = ring->xData();
452  const double *srcYData = ring->yData();
453  double x0 = *srcXData++;
454  double y0 = *srcYData++;
455  for ( int i = 1; i < numPoints; ++i )
456  {
457  const double x1 = *srcXData++;
458  const double y1 = *srcYData++;
459  const double d = ( x0 - x1 ) * ( x0 - x1 ) + ( y0 - y1 ) * ( y0 - y1 );
460  if ( d < min_d )
461  min_d = d;
462  x0 = x1;
463  y0 = y1;
464  }
465  }
466 
467  return min_d != 1e20 ? std::sqrt( min_d ) : 1e20;
468 }
469 
470 void QgsTessellator::addPolygon( const QgsPolygon &polygon, float extrusionHeight )
471 {
472  const QgsLineString *exterior = qgsgeometry_cast< const QgsLineString * >( polygon.exteriorRing() );
473  if ( !exterior )
474  return;
475 
476  const QVector3D pNormal = !mNoZ ? _calculateNormal( exterior, mOriginX, mOriginY, mInvertNormals ) : QVector3D();
477  const int pCount = exterior->numPoints();
478  if ( pCount == 0 )
479  return;
480 
481  float zMin = std::numeric_limits<float>::max();
482  float zMaxBase = -std::numeric_limits<float>::max();
483  float zMaxExtruded = -std::numeric_limits<float>::max();
484 
485  const float scale = mBounds.isNull() ? 1.0 : std::max( 10000.0 / mBounds.width(), 10000.0 / mBounds.height() );
486 
487  std::unique_ptr<QMatrix4x4> toNewBase, toOldBase;
488  QgsPoint ptStart, pt0;
489  std::unique_ptr<QgsPolygon> polygonNew;
490  auto rotatePolygonToXYPlane = [&]()
491  {
492  if ( !mNoZ && pNormal != QVector3D( 0, 0, 1 ) )
493  {
494  // this is not a horizontal plane - need to reproject the polygon to a new base so that
495  // we can do the triangulation in a plane
496  QVector3D pXVector, pYVector;
497  _normalVectorToXYVectors( pNormal, pXVector, pYVector );
498 
499  // so now we have three orthogonal unit vectors defining new base
500  // let's build transform matrix. We actually need just a 3x3 matrix,
501  // but Qt does not have good support for it, so using 4x4 matrix instead.
502  toNewBase.reset( new QMatrix4x4(
503  pXVector.x(), pXVector.y(), pXVector.z(), 0,
504  pYVector.x(), pYVector.y(), pYVector.z(), 0,
505  pNormal.x(), pNormal.y(), pNormal.z(), 0,
506  0, 0, 0, 0 ) );
507 
508  // our 3x3 matrix is orthogonal, so for inverse we only need to transpose it
509  toOldBase.reset( new QMatrix4x4( toNewBase->transposed() ) );
510  }
511 
512  ptStart = QgsPoint( exterior->startPoint() );
513  pt0 = QgsPoint( QgsWkbTypes::PointZ, ptStart.x(), ptStart.y(), std::isnan( ptStart.z() ) ? 0 : ptStart.z() );
514 
515  // subtract ptFirst from geometry for better numerical stability in triangulation
516  // and apply new 3D vector base if the polygon is not horizontal
517 
518  polygonNew.reset( _transform_polygon_to_new_base( polygon, pt0, toNewBase.get(), scale ) );
519  };
520 
521  if ( !mNoZ && !qgsDoubleNear( pNormal.length(), 1, 0.001 ) )
522  return; // this should not happen - pNormal should be normalized to unit length
523 
524  const QVector3D upVector( 0, 0, 1 );
525  const float pNormalUpVectorDotProduct = QVector3D::dotProduct( upVector, pNormal );
526  const float radsBetwwenUpNormal = qAcos( pNormalUpVectorDotProduct );
527 
528  const float detectionDelta = qDegreesToRadians( 10.0f );
529  int facade = 0;
530  if ( radsBetwwenUpNormal > M_PI_2 - detectionDelta && radsBetwwenUpNormal < M_PI_2 + detectionDelta ) facade = 1;
531  else if ( radsBetwwenUpNormal > - M_PI_2 - detectionDelta && radsBetwwenUpNormal < -M_PI_2 + detectionDelta ) facade = 1;
532  else facade = 2;
533 
534  if ( pCount == 4 && polygon.numInteriorRings() == 0 && ( mTessellatedFacade & facade ) )
535  {
536  QgsLineString *triangle = nullptr;
537  if ( mAddTextureCoords )
538  {
539  rotatePolygonToXYPlane();
540  triangle = qgsgeometry_cast< QgsLineString * >( polygonNew->exteriorRing() );
541  Q_ASSERT( polygonNew->exteriorRing()->numPoints() >= 3 );
542  }
543 
544  // polygon is a triangle - write vertices to the output data array without triangulation
545  const double *xData = exterior->xData();
546  const double *yData = exterior->yData();
547  const double *zData = !mNoZ ? exterior->zData() : nullptr;
548  for ( int i = 0; i < 3; i++ )
549  {
550  const float z = !zData ? 0 : *zData;
551  if ( z < zMin )
552  zMin = z;
553  if ( z > zMaxBase )
554  zMaxBase = z;
555  if ( z > zMaxExtruded )
556  zMaxExtruded = z;
557 
558  mData << *xData - mOriginX << z << - *yData + mOriginY;
559  if ( mAddNormals )
560  mData << pNormal.x() << pNormal.z() << - pNormal.y();
561  if ( mAddTextureCoords )
562  {
563  std::pair<float, float> p( triangle->xAt( i ), triangle->yAt( i ) );
564  if ( facade & 1 )
565  {
566  p = rotateCoords( p.first, p.second, 0.0f, 0.0f, mTextureRotation );
567  }
568  else if ( facade & 2 )
569  {
570  p = rotateCoords( p.first, p.second, 0.0f, 0.0f, mTextureRotation );
571  }
572  mData << p.first << p.second;
573  }
574  xData++; yData++;
575  // zData can be nullptr if mNoZ or triangle is 2D
576  if ( zData )
577  zData++;
578  }
579 
580  if ( mAddBackFaces )
581  {
582  // the same triangle with reversed order of coordinates and inverted normal
583  for ( int i = 2; i >= 0; i-- )
584  {
585  mData << exterior->xAt( i ) - mOriginX << ( mNoZ ? 0 : exterior->zAt( i ) ) << - exterior->yAt( i ) + mOriginY;
586  if ( mAddNormals )
587  mData << -pNormal.x() << -pNormal.z() << pNormal.y();
588  if ( mAddTextureCoords )
589  {
590  std::pair<float, float> p( triangle->xAt( i ), triangle->yAt( i ) );
591  if ( facade & 1 )
592  {
593  p = rotateCoords( p.first, p.second, 0.0f, 0.0f, mTextureRotation );
594  }
595  else if ( facade & 2 )
596  {
597  p = rotateCoords( p.first, p.second, 0.0f, 0.0f, mTextureRotation );
598  }
599  mData << p.first << p.second;
600  }
601  }
602  }
603  }
604  else if ( mTessellatedFacade & facade )
605  {
606 
607  rotatePolygonToXYPlane();
608 
609  if ( _minimum_distance_between_coordinates( *polygonNew ) < 0.001 )
610  {
611  // when the distances between coordinates of input points are very small,
612  // the triangulation likes to crash on numerical errors - when the distances are ~ 1e-5
613  // Assuming that the coordinates should be in a projected CRS, we should be able
614  // to simplify geometries that may cause problems and avoid possible crashes
615  const QgsGeometry polygonSimplified = QgsGeometry( polygonNew->clone() ).simplify( 0.001 );
616  if ( polygonSimplified.isNull() )
617  {
618  QgsMessageLog::logMessage( QObject::tr( "geometry simplification failed - skipping" ), QObject::tr( "3D" ) );
619  return;
620  }
621  const QgsPolygon *polygonSimplifiedData = qgsgeometry_cast<const QgsPolygon *>( polygonSimplified.constGet() );
622  if ( polygonSimplifiedData == nullptr || _minimum_distance_between_coordinates( *polygonSimplifiedData ) < 0.001 )
623  {
624  // Failed to fix that. It could be a really tiny geometry... or maybe they gave us
625  // geometry in unprojected lat/lon coordinates
626  QgsMessageLog::logMessage( QObject::tr( "geometry's coordinates are too close to each other and simplification failed - skipping" ), QObject::tr( "3D" ) );
627  return;
628  }
629  else
630  {
631  polygonNew.reset( polygonSimplifiedData->clone() );
632  }
633  }
634 
635  QList< std::vector<p2t::Point *> > polylinesToDelete;
636  QHash<p2t::Point *, float> z;
637 
638  // polygon exterior
639  std::vector<p2t::Point *> polyline;
640  _ringToPoly2tri( qgsgeometry_cast< const QgsLineString * >( polygonNew->exteriorRing() ), polyline, mNoZ ? nullptr : &z );
641  polylinesToDelete << polyline;
642 
643  std::unique_ptr<p2t::CDT> cdt( new p2t::CDT( polyline ) );
644 
645  // polygon holes
646  for ( int i = 0; i < polygonNew->numInteriorRings(); ++i )
647  {
648  std::vector<p2t::Point *> holePolyline;
649  const QgsLineString *hole = qgsgeometry_cast< const QgsLineString *>( polygonNew->interiorRing( i ) );
650 
651  _ringToPoly2tri( hole, holePolyline, mNoZ ? nullptr : &z );
652 
653  cdt->AddHole( holePolyline );
654  polylinesToDelete << holePolyline;
655  }
656 
657  // run triangulation and write vertices to the output data array
658  try
659  {
660  cdt->Triangulate();
661 
662  std::vector<p2t::Triangle *> triangles = cdt->GetTriangles();
663 
664  mData.reserve( mData.size() + 3 * triangles.size() * ( stride() / sizeof( float ) ) );
665  for ( size_t i = 0; i < triangles.size(); ++i )
666  {
667  p2t::Triangle *t = triangles[i];
668  for ( int j = 0; j < 3; ++j )
669  {
670  p2t::Point *p = t->GetPoint( j );
671  QVector4D pt( p->x, p->y, mNoZ ? 0 : z[p], 0 );
672  if ( toOldBase )
673  pt = *toOldBase * pt;
674  const double fx = ( pt.x() / scale ) - mOriginX + pt0.x();
675  const double fy = ( pt.y() / scale ) - mOriginY + pt0.y();
676  const double baseHeight = mNoZ ? 0 : ( pt.z() + pt0.z() );
677  const double fz = mNoZ ? 0 : ( pt.z() + extrusionHeight + pt0.z() );
678  if ( baseHeight < zMin )
679  zMin = baseHeight;
680  if ( baseHeight > zMaxBase )
681  zMaxBase = baseHeight;
682  if ( fz > zMaxExtruded )
683  zMaxExtruded = fz;
684 
685  mData << fx << fz << -fy;
686  if ( mAddNormals )
687  mData << pNormal.x() << pNormal.z() << - pNormal.y();
688  if ( mAddTextureCoords )
689  {
690  const std::pair<float, float> pr = rotateCoords( p->x, p->y, 0.0f, 0.0f, mTextureRotation );
691  mData << pr.first << pr.second;
692  }
693  }
694 
695  if ( mAddBackFaces )
696  {
697  // the same triangle with reversed order of coordinates and inverted normal
698  for ( int j = 2; j >= 0; --j )
699  {
700  p2t::Point *p = t->GetPoint( j );
701  QVector4D pt( p->x, p->y, mNoZ ? 0 : z[p], 0 );
702  if ( toOldBase )
703  pt = *toOldBase * pt;
704  const double fx = ( pt.x() / scale ) - mOriginX + pt0.x();
705  const double fy = ( pt.y() / scale ) - mOriginY + pt0.y();
706  const double fz = mNoZ ? 0 : ( pt.z() + extrusionHeight + pt0.z() );
707  mData << fx << fz << -fy;
708  if ( mAddNormals )
709  mData << -pNormal.x() << -pNormal.z() << pNormal.y();
710  if ( mAddTextureCoords )
711  {
712  const std::pair<float, float> pr = rotateCoords( p->x, p->y, 0.0f, 0.0f, mTextureRotation );
713  mData << pr.first << pr.second;
714  }
715  }
716  }
717  }
718  }
719  catch ( ... )
720  {
721  QgsMessageLog::logMessage( QObject::tr( "Triangulation failed. Skipping polygon…" ), QObject::tr( "3D" ) );
722  }
723 
724  for ( int i = 0; i < polylinesToDelete.count(); ++i )
725  qDeleteAll( polylinesToDelete[i] );
726  }
727 
728  // add walls if extrusion is enabled
729  if ( extrusionHeight != 0 && ( mTessellatedFacade & 1 ) )
730  {
731  _makeWalls( *exterior, false, extrusionHeight, mData, mAddNormals, mAddTextureCoords, mOriginX, mOriginY, mTextureRotation );
732 
733  for ( int i = 0; i < polygon.numInteriorRings(); ++i )
734  _makeWalls( *qgsgeometry_cast< const QgsLineString * >( polygon.interiorRing( i ) ), true, extrusionHeight, mData, mAddNormals, mAddTextureCoords, mOriginX, mOriginY, mTextureRotation );
735 
736  if ( zMaxBase + extrusionHeight > zMaxExtruded )
737  zMaxExtruded = zMaxBase + extrusionHeight;
738  }
739 
740  if ( zMin < mZMin )
741  mZMin = zMin;
742  if ( zMaxExtruded > mZMax )
743  mZMax = zMaxExtruded;
744  if ( zMaxBase > mZMax )
745  mZMax = zMaxBase;
746 }
747 
749 {
750  return mData.size() / ( stride() / sizeof( float ) );
751 }
752 
753 std::unique_ptr<QgsMultiPolygon> QgsTessellator::asMultiPolygon() const
754 {
755  std::unique_ptr< QgsMultiPolygon > mp = std::make_unique< QgsMultiPolygon >();
756  const auto nVals = mData.size();
757  mp->reserve( nVals / 9 );
758  for ( auto i = decltype( nVals ) {0}; i + 8 < nVals; i += 9 )
759  {
760  // tessellator geometry is x, z, -y
761  const QgsPoint p1( mData[i + 0], -mData[i + 2], mData[i + 1] );
762  const QgsPoint p2( mData[i + 3], -mData[i + 5], mData[i + 4] );
763  const QgsPoint p3( mData[i + 6], -mData[i + 8], mData[i + 7] );
764  mp->addGeometry( new QgsTriangle( p1, p2, p3 ) );
765  }
766  return mp;
767 }
VertexType
Types of vertex.
Definition: qgis.h:1192
bool is3D() const SIP_HOLDGIL
Returns true if the geometry is 3D and contains a z-value.
QgsWkbTypes::Type wkbType() const SIP_HOLDGIL
Returns the WKB type of the geometry.
const QgsCurve * interiorRing(int i) const SIP_HOLDGIL
Retrieves an interior ring from the curve polygon.
const QgsCurve * exteriorRing() const SIP_HOLDGIL
Returns the curve polygon's exterior ring.
int numInteriorRings() const SIP_HOLDGIL
Returns the number of interior rings contained with the curve polygon.
Abstract base class for curved geometry type.
Definition: qgscurve.h:36
virtual int numPoints() const =0
Returns the number of points in the curve.
virtual bool pointAt(int node, QgsPoint &point, Qgis::VertexType &type) const =0
Returns the point and vertex id of a point within the curve.
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:125
const QgsAbstractGeometry * constGet() const SIP_HOLDGIL
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
Q_GADGET bool isNull
Definition: qgsgeometry.h:127
QgsGeometry simplify(double tolerance) const
Returns a simplified version of this geometry using a specified tolerance value.
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:44
QgsPoint startPoint() const override SIP_HOLDGIL
Returns the starting point of the curve.
const double * yData() const
Returns a const pointer to the y vertex data.
int numPoints() const override SIP_HOLDGIL
Returns the number of points in the curve.
QgsPoint pointN(int i) const
Returns the specified point from inside the line string.
bool isEmpty() const override SIP_HOLDGIL
Returns true if the geometry is empty.
double yAt(int index) const override
Returns the y-coordinate of the specified node in the line string.
double zAt(int index) const
Returns the z-coordinate of the specified node in the line string.
const double * xData() const
Returns a const pointer to the x vertex data.
const double * zData() const
Returns a const pointer to the z vertex data, or nullptr if the linestring does not have z values.
double xAt(int index) const override
Returns the x-coordinate of the specified node in the line string.
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::MessageLevel::Warning, bool notifyUser=true)
Adds a message to the log instance (and creates it if necessary).
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
void setX(double x) SIP_HOLDGIL
Sets the point's x-coordinate.
Definition: qgspoint.h:280
Q_GADGET double x
Definition: qgspoint.h:52
void setY(double y) SIP_HOLDGIL
Sets the point's y-coordinate.
Definition: qgspoint.h:291
double z
Definition: qgspoint.h:54
double y
Definition: qgspoint.h:53
Polygon geometry type.
Definition: qgspolygon.h:34
void setExteriorRing(QgsCurve *ring) override
Sets the exterior ring of the polygon.
Definition: qgspolygon.cpp:219
void addInteriorRing(QgsCurve *ring) override
Adds an interior ring to the geometry (takes ownership)
Definition: qgspolygon.cpp:188
QgsPolygon * clone() const override
Clones the geometry by performing a deep copy.
Definition: qgspolygon.cpp:54
A rectangle specified with double values.
Definition: qgsrectangle.h:42
bool isNull() const
Test if the rectangle is null (all coordinates zero or after call to setMinimal()).
Definition: qgsrectangle.h:479
double height() const SIP_HOLDGIL
Returns the height of the rectangle.
Definition: qgsrectangle.h:230
double width() const SIP_HOLDGIL
Returns the width of the rectangle.
Definition: qgsrectangle.h:223
std::unique_ptr< QgsMultiPolygon > asMultiPolygon() const
Returns the triangulation as a multipolygon geometry.
QgsTessellator(double originX, double originY, bool addNormals, bool invertNormals=false, bool addBackFaces=false, bool noZ=false, bool addTextureCoords=false, int facade=3, float textureRotation=0.0f)
Creates tessellator with a specified origin point of the world (in map coordinates)
int stride() const
Returns size of one vertex entry in bytes.
void addPolygon(const QgsPolygon &polygon, float extrusionHeight)
Tessellates a triangle and adds its vertex entries to the output data array.
int dataVerticesCount() const
Returns the number of vertices stored in the output data array.
Triangle geometry type.
Definition: qgstriangle.h:34
static bool hasZ(Type type) SIP_HOLDGIL
Tests whether a WKB type contains the z-dimension.
Definition: qgswkbtypes.h:1080
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:1990
double _round_coord(double x)
double _minimum_distance_between_coordinates(const QgsPolygon &polygon)
std::size_t operator()(const std::pair< float, float > pair) const