QGIS API Documentation  3.22.4-Białowieża (ce8e65e95e)
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 zMax = std::numeric_limits<float>::min();
483 
484  const float scale = mBounds.isNull() ? 1.0 : std::max( 10000.0 / mBounds.width(), 10000.0 / mBounds.height() );
485 
486  std::unique_ptr<QMatrix4x4> toNewBase, toOldBase;
487  QgsPoint ptStart, pt0;
488  std::unique_ptr<QgsPolygon> polygonNew;
489  auto rotatePolygonToXYPlane = [&]()
490  {
491  if ( !mNoZ && pNormal != QVector3D( 0, 0, 1 ) )
492  {
493  // this is not a horizontal plane - need to reproject the polygon to a new base so that
494  // we can do the triangulation in a plane
495  QVector3D pXVector, pYVector;
496  _normalVectorToXYVectors( pNormal, pXVector, pYVector );
497 
498  // so now we have three orthogonal unit vectors defining new base
499  // let's build transform matrix. We actually need just a 3x3 matrix,
500  // but Qt does not have good support for it, so using 4x4 matrix instead.
501  toNewBase.reset( new QMatrix4x4(
502  pXVector.x(), pXVector.y(), pXVector.z(), 0,
503  pYVector.x(), pYVector.y(), pYVector.z(), 0,
504  pNormal.x(), pNormal.y(), pNormal.z(), 0,
505  0, 0, 0, 0 ) );
506 
507  // our 3x3 matrix is orthogonal, so for inverse we only need to transpose it
508  toOldBase.reset( new QMatrix4x4( toNewBase->transposed() ) );
509  }
510 
511  ptStart = QgsPoint( exterior->startPoint() );
512  pt0 = QgsPoint( QgsWkbTypes::PointZ, ptStart.x(), ptStart.y(), std::isnan( ptStart.z() ) ? 0 : ptStart.z() );
513 
514  // subtract ptFirst from geometry for better numerical stability in triangulation
515  // and apply new 3D vector base if the polygon is not horizontal
516 
517  polygonNew.reset( _transform_polygon_to_new_base( polygon, pt0, toNewBase.get(), scale ) );
518  };
519 
520  if ( !mNoZ && !qgsDoubleNear( pNormal.length(), 1, 0.001 ) )
521  return; // this should not happen - pNormal should be normalized to unit length
522 
523  const QVector3D upVector( 0, 0, 1 );
524  const float pNormalUpVectorDotProduct = QVector3D::dotProduct( upVector, pNormal );
525  const float radsBetwwenUpNormal = qAcos( pNormalUpVectorDotProduct );
526 
527  const float detectionDelta = qDegreesToRadians( 10.0f );
528  int facade = 0;
529  if ( radsBetwwenUpNormal > M_PI_2 - detectionDelta && radsBetwwenUpNormal < M_PI_2 + detectionDelta ) facade = 1;
530  else if ( radsBetwwenUpNormal > - M_PI_2 - detectionDelta && radsBetwwenUpNormal < -M_PI_2 + detectionDelta ) facade = 1;
531  else facade = 2;
532 
533  if ( pCount == 4 && polygon.numInteriorRings() == 0 && ( mTessellatedFacade & facade ) )
534  {
535  QgsLineString *triangle = nullptr;
536  if ( mAddTextureCoords )
537  {
538  rotatePolygonToXYPlane();
539  triangle = qgsgeometry_cast< QgsLineString * >( polygonNew->exteriorRing() );
540  Q_ASSERT( polygonNew->exteriorRing()->numPoints() >= 3 );
541  }
542 
543  // polygon is a triangle - write vertices to the output data array without triangulation
544  const double *xData = exterior->xData();
545  const double *yData = exterior->yData();
546  const double *zData = !mNoZ ? exterior->zData() : nullptr;
547  for ( int i = 0; i < 3; i++ )
548  {
549  const float z = !zData ? 0 : *zData;
550  if ( z < zMin )
551  zMin = z;
552  if ( z > zMax )
553  zMax = z;
554 
555  mData << *xData - mOriginX << z << - *yData + mOriginY;
556  if ( mAddNormals )
557  mData << pNormal.x() << pNormal.z() << - pNormal.y();
558  if ( mAddTextureCoords )
559  {
560  std::pair<float, float> p( triangle->xAt( i ), triangle->yAt( i ) );
561  if ( facade & 1 )
562  {
563  p = rotateCoords( p.first, p.second, 0.0f, 0.0f, mTextureRotation );
564  }
565  else if ( facade & 2 )
566  {
567  p = rotateCoords( p.first, p.second, 0.0f, 0.0f, mTextureRotation );
568  }
569  mData << p.first << p.second;
570  }
571  xData++; yData++;
572  // zData can be nullptr if mNoZ or triangle is 2D
573  if ( zData )
574  zData++;
575  }
576 
577  if ( mAddBackFaces )
578  {
579  // the same triangle with reversed order of coordinates and inverted normal
580  for ( int i = 2; i >= 0; i-- )
581  {
582  mData << exterior->xAt( i ) - mOriginX << ( mNoZ ? 0 : exterior->zAt( i ) ) << - exterior->yAt( i ) + mOriginY;
583  if ( mAddNormals )
584  mData << -pNormal.x() << -pNormal.z() << pNormal.y();
585  if ( mAddTextureCoords )
586  {
587  std::pair<float, float> p( triangle->xAt( i ), triangle->yAt( i ) );
588  if ( facade & 1 )
589  {
590  p = rotateCoords( p.first, p.second, 0.0f, 0.0f, mTextureRotation );
591  }
592  else if ( facade & 2 )
593  {
594  p = rotateCoords( p.first, p.second, 0.0f, 0.0f, mTextureRotation );
595  }
596  mData << p.first << p.second;
597  }
598  }
599  }
600  }
601  else if ( mTessellatedFacade & facade )
602  {
603 
604  rotatePolygonToXYPlane();
605 
606  if ( _minimum_distance_between_coordinates( *polygonNew ) < 0.001 )
607  {
608  // when the distances between coordinates of input points are very small,
609  // the triangulation likes to crash on numerical errors - when the distances are ~ 1e-5
610  // Assuming that the coordinates should be in a projected CRS, we should be able
611  // to simplify geometries that may cause problems and avoid possible crashes
612  const QgsGeometry polygonSimplified = QgsGeometry( polygonNew->clone() ).simplify( 0.001 );
613  if ( polygonSimplified.isNull() )
614  {
615  QgsMessageLog::logMessage( QObject::tr( "geometry simplification failed - skipping" ), QObject::tr( "3D" ) );
616  return;
617  }
618  const QgsPolygon *polygonSimplifiedData = qgsgeometry_cast<const QgsPolygon *>( polygonSimplified.constGet() );
619  if ( _minimum_distance_between_coordinates( *polygonSimplifiedData ) < 0.001 )
620  {
621  // Failed to fix that. It could be a really tiny geometry... or maybe they gave us
622  // geometry in unprojected lat/lon coordinates
623  QgsMessageLog::logMessage( QObject::tr( "geometry's coordinates are too close to each other and simplification failed - skipping" ), QObject::tr( "3D" ) );
624  return;
625  }
626  else
627  {
628  polygonNew.reset( polygonSimplifiedData->clone() );
629  }
630  }
631 
632  QList< std::vector<p2t::Point *> > polylinesToDelete;
633  QHash<p2t::Point *, float> z;
634 
635  // polygon exterior
636  std::vector<p2t::Point *> polyline;
637  _ringToPoly2tri( qgsgeometry_cast< const QgsLineString * >( polygonNew->exteriorRing() ), polyline, mNoZ ? nullptr : &z );
638  polylinesToDelete << polyline;
639 
640  std::unique_ptr<p2t::CDT> cdt( new p2t::CDT( polyline ) );
641 
642  // polygon holes
643  for ( int i = 0; i < polygonNew->numInteriorRings(); ++i )
644  {
645  std::vector<p2t::Point *> holePolyline;
646  const QgsLineString *hole = qgsgeometry_cast< const QgsLineString *>( polygonNew->interiorRing( i ) );
647 
648  _ringToPoly2tri( hole, holePolyline, mNoZ ? nullptr : &z );
649 
650  cdt->AddHole( holePolyline );
651  polylinesToDelete << holePolyline;
652  }
653 
654  // run triangulation and write vertices to the output data array
655  try
656  {
657  cdt->Triangulate();
658 
659  std::vector<p2t::Triangle *> triangles = cdt->GetTriangles();
660 
661  mData.reserve( mData.size() + 3 * triangles.size() * ( stride() / sizeof( float ) ) );
662  for ( size_t i = 0; i < triangles.size(); ++i )
663  {
664  p2t::Triangle *t = triangles[i];
665  for ( int j = 0; j < 3; ++j )
666  {
667  p2t::Point *p = t->GetPoint( j );
668  QVector4D pt( p->x, p->y, mNoZ ? 0 : z[p], 0 );
669  if ( toOldBase )
670  pt = *toOldBase * pt;
671  const double fx = ( pt.x() / scale ) - mOriginX + pt0.x();
672  const double fy = ( pt.y() / scale ) - mOriginY + pt0.y();
673  const double fz = mNoZ ? 0 : ( pt.z() + extrusionHeight + pt0.z() );
674  if ( fz < zMin )
675  zMin = fz;
676  if ( fz > zMax )
677  zMax = fz;
678 
679  mData << fx << fz << -fy;
680  if ( mAddNormals )
681  mData << pNormal.x() << pNormal.z() << - pNormal.y();
682  if ( mAddTextureCoords )
683  {
684  const std::pair<float, float> pr = rotateCoords( p->x, p->y, 0.0f, 0.0f, mTextureRotation );
685  mData << pr.first << pr.second;
686  }
687  }
688 
689  if ( mAddBackFaces )
690  {
691  // the same triangle with reversed order of coordinates and inverted normal
692  for ( int j = 2; j >= 0; --j )
693  {
694  p2t::Point *p = t->GetPoint( j );
695  QVector4D pt( p->x, p->y, mNoZ ? 0 : z[p], 0 );
696  if ( toOldBase )
697  pt = *toOldBase * pt;
698  const double fx = ( pt.x() / scale ) - mOriginX + pt0.x();
699  const double fy = ( pt.y() / scale ) - mOriginY + pt0.y();
700  const double fz = mNoZ ? 0 : ( pt.z() + extrusionHeight + pt0.z() );
701  mData << fx << fz << -fy;
702  if ( mAddNormals )
703  mData << -pNormal.x() << -pNormal.z() << pNormal.y();
704  if ( mAddTextureCoords )
705  {
706  const std::pair<float, float> pr = rotateCoords( p->x, p->y, 0.0f, 0.0f, mTextureRotation );
707  mData << pr.first << pr.second;
708  }
709  }
710  }
711  }
712  }
713  catch ( ... )
714  {
715  QgsMessageLog::logMessage( QObject::tr( "Triangulation failed. Skipping polygon…" ), QObject::tr( "3D" ) );
716  }
717 
718  for ( int i = 0; i < polylinesToDelete.count(); ++i )
719  qDeleteAll( polylinesToDelete[i] );
720  }
721 
722  // add walls if extrusion is enabled
723  if ( extrusionHeight != 0 && ( mTessellatedFacade & 1 ) )
724  {
725  _makeWalls( *exterior, false, extrusionHeight, mData, mAddNormals, mAddTextureCoords, mOriginX, mOriginY, mTextureRotation );
726 
727  for ( int i = 0; i < polygon.numInteriorRings(); ++i )
728  _makeWalls( *qgsgeometry_cast< const QgsLineString * >( polygon.interiorRing( i ) ), true, extrusionHeight, mData, mAddNormals, mAddTextureCoords, mOriginX, mOriginY, mTextureRotation );
729 
730  zMax += extrusionHeight;
731  }
732 
733  if ( zMin < mZMin )
734  mZMin = zMin;
735  if ( zMax > mZMax )
736  mZMax = zMax;
737 }
738 
740 {
741  return mData.size() / ( stride() / sizeof( float ) );
742 }
743 
744 std::unique_ptr<QgsMultiPolygon> QgsTessellator::asMultiPolygon() const
745 {
746  std::unique_ptr< QgsMultiPolygon > mp = std::make_unique< QgsMultiPolygon >();
747  const auto nVals = mData.size();
748  mp->reserve( nVals / 9 );
749  for ( auto i = decltype( nVals ) {0}; i + 8 < nVals; i += 9 )
750  {
751  // tessellator geometry is x, z, -y
752  const QgsPoint p1( mData[i + 0], -mData[i + 2], mData[i + 1] );
753  const QgsPoint p2( mData[i + 3], -mData[i + 5], mData[i + 4] );
754  const QgsPoint p3( mData[i + 6], -mData[i + 8], mData[i + 7] );
755  mp->addGeometry( new QgsTriangle( p1, p2, p3 ) );
756  }
757  return mp;
758 }
VertexType
Types of vertex.
Definition: qgis.h:990
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:1246
double _round_coord(double x)
double _minimum_distance_between_coordinates(const QgsPolygon &polygon)
std::size_t operator()(const std::pair< float, float > pair) const