QGIS API Documentation 3.27.0-Master (95e00c50d2)
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
37static 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
49static 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
165QgsTessellator::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
180QgsTessellator::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
196void 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
205static bool _isRingCounterClockWise( const QgsCurve &ring )
206{
207 double a = 0;
208 const int count = ring.numPoints();
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
221static 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
245static 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
305static 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
337static 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
370inline 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
377static 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
425static 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
470void 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
753std::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:1518
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:164
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:166
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:45
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.
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.
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.
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:2527
double _round_coord(double x)
double _minimum_distance_between_coordinates(const QgsPolygon &polygon)
std::size_t operator()(const std::pair< float, float > pair) const