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