Branch data Line data Source code
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 : 0 : static std::pair<float, float> rotateCoords( float x, float y, float origin_x, float origin_y, float r )
38 : : {
39 : 0 : r = qDegreesToRadians( r );
40 : 0 : 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 : 0 : float x1 = origin_x + x0 * qCos( r ) - y0 * qSin( r );
45 : 0 : float y1 = origin_y + x0 * qSin( r ) + y0 * qCos( r );
46 : 0 : return std::make_pair( x1, y1 );
47 : : }
48 : :
49 : 0 : 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 : 0 : float dx = x1 - x0;
52 : 0 : float dy = -( y1 - y0 );
53 : :
54 : : // perpendicular vector in plane to [x,y] is [-y,x]
55 : 0 : QVector3D vn( -dy, 0, dx );
56 : 0 : vn = -vn;
57 : 0 : vn.normalize();
58 : :
59 : : float u0, v0;
60 : : float u1, v1;
61 : : float u2, v2;
62 : : float u3, v3;
63 : :
64 : 0 : QVector<double> textureCoordinates;
65 : 0 : textureCoordinates.reserve( 12 );
66 : : // select which side of the coordinates to use (x, z or y, z) depending on which side is smaller
67 : 0 : if ( fabsf( dy ) <= fabsf( dx ) )
68 : : {
69 : : // consider x and z as the texture coordinates
70 : 0 : u0 = x0;
71 : 0 : v0 = z0 + height;
72 : :
73 : 0 : u1 = x1;
74 : 0 : v1 = z1 + height;
75 : :
76 : 0 : u2 = x0;
77 : 0 : v2 = z0;
78 : :
79 : 0 : u3 = x1;
80 : 0 : v3 = z1;
81 : 0 : }
82 : : else
83 : : {
84 : : // consider y and z as the texture coowallsTextureRotationrdinates
85 : 0 : u0 = -y0;
86 : 0 : v0 = z0 + height;
87 : :
88 : 0 : u1 = -y1;
89 : 0 : v1 = z1 + height;
90 : :
91 : 0 : u2 = -y0;
92 : 0 : v2 = z0;
93 : :
94 : 0 : u3 = -y1;
95 : 0 : v3 = z1;
96 : : }
97 : :
98 : 0 : textureCoordinates.push_back( u0 );
99 : 0 : textureCoordinates.push_back( v0 );
100 : :
101 : 0 : textureCoordinates.push_back( u1 );
102 : 0 : textureCoordinates.push_back( v1 );
103 : :
104 : 0 : textureCoordinates.push_back( u2 );
105 : 0 : textureCoordinates.push_back( v2 );
106 : :
107 : 7 : textureCoordinates.push_back( u2 );
108 : 7 : textureCoordinates.push_back( v2 );
109 : :
110 : 0 : textureCoordinates.push_back( u1 );
111 : 0 : textureCoordinates.push_back( v1 );
112 : :
113 : 0 : textureCoordinates.push_back( u3 );
114 : 0 : textureCoordinates.push_back( v3 );
115 : :
116 : 0 : for ( int i = 0; i < textureCoordinates.size(); i += 2 )
117 : : {
118 : 0 : std::pair<float, float> rotated = rotateCoords( textureCoordinates[i], textureCoordinates[i + 1], 0, 0, textureRotation );
119 : 0 : textureCoordinates[i] = rotated.first;
120 : 0 : textureCoordinates[i + 1] = rotated.second;
121 : 0 : }
122 : :
123 : : // triangle 1
124 : : // vertice 1
125 : 0 : data << x0 << z0 + height << -y0;
126 : 0 : if ( addNormals )
127 : 0 : data << vn.x() << vn.y() << vn.z();
128 : 0 : if ( addTextureCoords )
129 : 0 : data << textureCoordinates[0] << textureCoordinates[1];
130 : : // vertice 2
131 : 0 : data << x1 << z1 + height << -y1;
132 : 0 : if ( addNormals )
133 : 0 : data << vn.x() << vn.y() << vn.z();
134 : 0 : if ( addTextureCoords )
135 : 0 : data << textureCoordinates[2] << textureCoordinates[3];
136 : : // verice 3
137 : 0 : data << x0 << z0 << -y0;
138 : 0 : if ( addNormals )
139 : 0 : data << vn.x() << vn.y() << vn.z();
140 : 0 : if ( addTextureCoords )
141 : 0 : data << textureCoordinates[4] << textureCoordinates[5];
142 : :
143 : : // triangle 2
144 : : // vertice 1
145 : 0 : data << x0 << z0 << -y0;
146 : 0 : if ( addNormals )
147 : 0 : data << vn.x() << vn.y() << vn.z();
148 : 0 : if ( addTextureCoords )
149 : 0 : data << textureCoordinates[6] << textureCoordinates[7];
150 : : // vertice 2
151 : 0 : data << x1 << z1 + height << -y1;
152 : 0 : if ( addNormals )
153 : 0 : data << vn.x() << vn.y() << vn.z();
154 : 0 : if ( addTextureCoords )
155 : 0 : data << textureCoordinates[8] << textureCoordinates[9];
156 : : // vertice 3
157 : 0 : data << x1 << z1 << -y1;
158 : 0 : if ( addNormals )
159 : 0 : data << vn.x() << vn.y() << vn.z();
160 : 0 : if ( addTextureCoords )
161 : 0 : data << textureCoordinates[10] << textureCoordinates[11];
162 : 0 : }
163 : :
164 : :
165 : 0 : QgsTessellator::QgsTessellator( double originX, double originY, bool addNormals, bool invertNormals, bool addBackFaces, bool noZ,
166 : : bool addTextureCoords, int facade, float textureRotation )
167 : 0 : : mOriginX( originX )
168 : 0 : , mOriginY( originY )
169 : 0 : , mAddNormals( addNormals )
170 : 0 : , mInvertNormals( invertNormals )
171 : 0 : , mAddBackFaces( addBackFaces )
172 : 0 : , mAddTextureCoords( addTextureCoords )
173 : 0 : , mNoZ( noZ )
174 : 0 : , mTessellatedFacade( facade )
175 : 0 : , mTextureRotation( textureRotation )
176 : : {
177 : 0 : init();
178 : 0 : }
179 : :
180 : 7 : QgsTessellator::QgsTessellator( const QgsRectangle &bounds, bool addNormals, bool invertNormals, bool addBackFaces, bool noZ,
181 : : bool addTextureCoords, int facade, float textureRotation )
182 : 7 : : mBounds( bounds )
183 : 7 : , mOriginX( mBounds.xMinimum() )
184 : 7 : , mOriginY( mBounds.yMinimum() )
185 : 7 : , mAddNormals( addNormals )
186 : 7 : , mInvertNormals( invertNormals )
187 : 7 : , mAddBackFaces( addBackFaces )
188 : 7 : , mAddTextureCoords( addTextureCoords )
189 : 7 : , mNoZ( noZ )
190 : 7 : , mTessellatedFacade( facade )
191 : 7 : , mTextureRotation( textureRotation )
192 : : {
193 : 7 : init();
194 : 7 : }
195 : :
196 : 7 : void QgsTessellator::init()
197 : : {
198 : 7 : mStride = 3 * sizeof( float );
199 : 7 : if ( mAddNormals )
200 : 0 : mStride += 3 * sizeof( float );
201 : 7 : if ( mAddTextureCoords )
202 : 0 : mStride += 2 * sizeof( float );
203 : 7 : }
204 : :
205 : 0 : static bool _isRingCounterClockWise( const QgsCurve &ring )
206 : : {
207 : 0 : double a = 0;
208 : 0 : int count = ring.numPoints();
209 : : QgsVertexId::VertexType vt;
210 : 0 : QgsPoint pt, ptPrev;
211 : 0 : ring.pointAt( 0, ptPrev, vt );
212 : 0 : for ( int i = 1; i < count + 1; ++i )
213 : : {
214 : 0 : ring.pointAt( i % count, pt, vt );
215 : 0 : a += ptPrev.x() * pt.y() - ptPrev.y() * pt.x();
216 : 0 : ptPrev = pt;
217 : 0 : }
218 : 0 : return a > 0; // clockwise if a is negative
219 : 0 : }
220 : :
221 : 0 : 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 : 0 : bool is_counter_clockwise = _isRingCounterClockWise( ring );
228 : :
229 : 0 : QgsPoint pt;
230 : 0 : QgsPoint ptPrev = ring.pointN( is_counter_clockwise == ccw ? 0 : ring.numPoints() - 1 );
231 : 0 : for ( int i = 1; i < ring.numPoints(); ++i )
232 : : {
233 : 0 : pt = ring.pointN( is_counter_clockwise == ccw ? i : ring.numPoints() - i - 1 );
234 : 0 : float x0 = ptPrev.x() - originX, y0 = ptPrev.y() - originY;
235 : 0 : float x1 = pt.x() - originX, y1 = pt.y() - originY;
236 : 0 : float z0 = std::isnan( ptPrev.z() ) ? 0 : ptPrev.z();
237 : 0 : float z1 = std::isnan( pt.z() ) ? 0 : pt.z();
238 : :
239 : : // make a quad
240 : 0 : make_quad( x0, y0, z0, x1, y1, z1, extrusionHeight, data, addNormals, addTextureCoords, textureRotation );
241 : 0 : ptPrev = pt;
242 : 0 : }
243 : 0 : }
244 : :
245 : 0 : 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 : 0 : if ( !QgsWkbTypes::hasZ( curve->wkbType() ) || curve->isEmpty() )
250 : 0 : 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 : 0 : bool sameZ = true;
256 : 0 : QgsPoint pt1 = curve->pointN( 0 );
257 : 0 : QgsPoint pt2;
258 : 0 : for ( int i = 1; i < curve->numPoints(); i++ )
259 : : {
260 : 0 : pt2 = curve->pointN( i );
261 : 0 : if ( pt1.z() != pt2.z() )
262 : : {
263 : 0 : sameZ = false;
264 : 0 : break;
265 : : }
266 : 0 : }
267 : 0 : if ( sameZ )
268 : 0 : 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 : 0 : double nx = 0, ny = 0, nz = 0;
276 : 0 : 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 : 0 : pt1.setX( pt1.x() - originX );
280 : 0 : pt1.setY( pt1.y() - originY );
281 : 0 : for ( int i = 1; i < curve->numPoints(); i++ )
282 : : {
283 : 0 : pt2 = curve->pointN( i );
284 : 0 : pt2.setX( pt2.x() - originX );
285 : 0 : pt2.setY( pt2.y() - originY );
286 : :
287 : 0 : if ( std::isnan( pt1.z() ) || std::isnan( pt2.z() ) )
288 : 0 : continue;
289 : :
290 : 0 : nx += ( pt1.y() - pt2.y() ) * ( pt1.z() + pt2.z() );
291 : 0 : ny += ( pt1.z() - pt2.z() ) * ( pt1.x() + pt2.x() );
292 : 0 : nz += ( pt1.x() - pt2.x() ) * ( pt1.y() + pt2.y() );
293 : :
294 : 0 : pt1 = pt2;
295 : 0 : }
296 : :
297 : 0 : QVector3D normal( nx, ny, nz );
298 : 0 : if ( invertNormal )
299 : 0 : normal = -normal;
300 : 0 : normal.normalize();
301 : 0 : return normal;
302 : 0 : }
303 : :
304 : :
305 : 0 : 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 : 0 : if ( pNormal.z() > 0.001 || pNormal.z() < -0.001 )
311 : : {
312 : 0 : pXVector = QVector3D( 1, 0, -pNormal.x() / pNormal.z() );
313 : 0 : }
314 : 0 : else if ( pNormal.y() > 0.001 || pNormal.y() < -0.001 )
315 : : {
316 : 0 : pXVector = QVector3D( 1, -pNormal.x() / pNormal.y(), 0 );
317 : 0 : }
318 : : else
319 : : {
320 : 0 : pXVector = QVector3D( -pNormal.y() / pNormal.x(), 1, 0 );
321 : : }
322 : 0 : pXVector.normalize();
323 : 0 : pYVector = QVector3D::normal( pNormal, pXVector );
324 : 0 : }
325 : :
326 : : struct float_pair_hash
327 : : {
328 : 91 : std::size_t operator()( const std::pair<float, float> pair ) const
329 : : {
330 : 91 : std::size_t h1 = std::hash<float>()( pair.first );
331 : 91 : std::size_t h2 = std::hash<float>()( pair.second );
332 : :
333 : 91 : return h1 ^ h2;
334 : : }
335 : : };
336 : :
337 : 26 : static void _ringToPoly2tri( const QgsLineString *ring, std::vector<p2t::Point *> &polyline, QHash<p2t::Point *, float> *zHash )
338 : : {
339 : 26 : const int pCount = ring->numPoints();
340 : :
341 : 26 : polyline.reserve( pCount );
342 : :
343 : 26 : const double *srcXData = ring->xData();
344 : 26 : const double *srcYData = ring->yData();
345 : 26 : const double *srcZData = ring->zData();
346 : 26 : std::unordered_set<std::pair<float, float>, float_pair_hash> foundPoints;
347 : :
348 : 117 : for ( int i = 0; i < pCount - 1; ++i )
349 : : {
350 : 91 : const float x = *srcXData++;
351 : 91 : const float y = *srcYData++;
352 : :
353 : 91 : auto res = foundPoints.insert( std::make_pair( x, y ) );
354 : 91 : if ( !res.second )
355 : : {
356 : : // already used this point, don't add a second time
357 : 0 : continue;
358 : : }
359 : :
360 : 91 : p2t::Point *pt2 = new p2t::Point( x, y );
361 : 91 : polyline.push_back( pt2 );
362 : 91 : if ( zHash )
363 : : {
364 : 0 : ( *zHash )[pt2] = *srcZData++;
365 : 0 : }
366 : 91 : }
367 : 26 : }
368 : :
369 : :
370 : 351 : inline double _round_coord( double x )
371 : : {
372 : 351 : const double exp = 1e10; // round to 10 decimal digits
373 : 351 : return round( x * exp ) / exp;
374 : : }
375 : :
376 : :
377 : 26 : static QgsCurve *_transform_ring_to_new_base( const QgsLineString &curve, const QgsPoint &pt0, const QMatrix4x4 *toNewBase, const float scale )
378 : : {
379 : 26 : int count = curve.numPoints();
380 : 26 : QVector<double> x;
381 : 26 : QVector<double> y;
382 : 26 : QVector<double> z;
383 : 26 : x.resize( count );
384 : 26 : y.resize( count );
385 : 26 : z.resize( count );
386 : 26 : double *xData = x.data();
387 : 26 : double *yData = y.data();
388 : 26 : double *zData = z.data();
389 : :
390 : 26 : const double *srcXData = curve.xData();
391 : 26 : const double *srcYData = curve.yData();
392 : 26 : const double *srcZData = curve.is3D() ? curve.zData() : nullptr;
393 : :
394 : 143 : for ( int i = 0; i < count; ++i )
395 : : {
396 : 234 : QVector4D v( *srcXData++ - pt0.x(),
397 : 117 : *srcYData++ - pt0.y(),
398 : 117 : srcZData ? *srcZData++ - pt0.z() : 0,
399 : : 0 );
400 : 117 : if ( toNewBase )
401 : 0 : v = toNewBase->map( v );
402 : :
403 : : // scale coordinates
404 : 117 : v.setX( v.x() * scale );
405 : 117 : 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 : 117 : *xData++ = _round_coord( v.x() );
418 : 117 : *yData++ = _round_coord( v.y() );
419 : 117 : *zData++ = _round_coord( v.z() );
420 : 117 : }
421 : 26 : return new QgsLineString( x, y, z );
422 : 26 : }
423 : :
424 : :
425 : 13 : static QgsPolygon *_transform_polygon_to_new_base( const QgsPolygon &polygon, const QgsPoint &pt0, const QMatrix4x4 *toNewBase, const float scale )
426 : : {
427 : 13 : QgsPolygon *p = new QgsPolygon;
428 : 13 : p->setExteriorRing( _transform_ring_to_new_base( *qgsgeometry_cast< const QgsLineString * >( polygon.exteriorRing() ), pt0, toNewBase, scale ) );
429 : 26 : for ( int i = 0; i < polygon.numInteriorRings(); ++i )
430 : 13 : p->addInteriorRing( _transform_ring_to_new_base( *qgsgeometry_cast< const QgsLineString * >( polygon.interiorRing( i ) ), pt0, toNewBase, scale ) );
431 : 13 : return p;
432 : 0 : }
433 : :
434 : 13 : static bool _check_intersecting_rings( const QgsPolygon &polygon )
435 : : {
436 : 13 : std::vector< std::unique_ptr< QgsGeometryEngine > > ringEngines;
437 : 13 : ringEngines.reserve( 1 + polygon.numInteriorRings() );
438 : 13 : ringEngines.emplace_back( QgsGeometry::createGeometryEngine( polygon.exteriorRing() ) );
439 : 26 : for ( int i = 0; i < polygon.numInteriorRings(); ++i )
440 : 13 : ringEngines.emplace_back( QgsGeometry::createGeometryEngine( polygon.interiorRing( i ) ) );
441 : :
442 : : // we need to make sure that the polygon has no rings with self-intersection: that may
443 : : // crash the tessellator. The original geometry maybe have been valid and the self-intersection
444 : : // was introduced when transforming to a new base (in a rare case when all points are not in the same plane)
445 : :
446 : 39 : for ( const std::unique_ptr< QgsGeometryEngine > &ring : ringEngines )
447 : : {
448 : 26 : if ( !ring->isSimple() )
449 : 0 : return false;
450 : : }
451 : :
452 : : // At this point we assume that input polygons are valid according to the OGC definition.
453 : : // This means e.g. no duplicate points, polygons are simple (no butterfly shaped polygon with self-intersection),
454 : : // internal rings are inside exterior rings, rings do not cross each other, no dangles.
455 : :
456 : : // There is however an issue with polygons where rings touch:
457 : : // +---+
458 : : // | |
459 : : // | +-+-+
460 : : // | | | |
461 : : // | +-+ |
462 : : // | |
463 : : // +-----+
464 : : // This is a valid polygon with one exterior and one interior ring that touch at one point,
465 : : // but poly2tri library does not allow interior rings touch each other or exterior ring.
466 : : // TODO: Handle the situation better - rather than just detecting the problem, try to fix
467 : : // it by converting touching rings into one ring.
468 : :
469 : 13 : if ( ringEngines.size() > 1 )
470 : : {
471 : 39 : for ( size_t i = 0; i < ringEngines.size(); ++i )
472 : : {
473 : 26 : std::unique_ptr< QgsGeometryEngine > &first = ringEngines.at( i );
474 : 26 : if ( polygon.numInteriorRings() > 1 )
475 : 0 : first->prepareGeometry();
476 : :
477 : : // TODO this is inefficient - QgsGeometryEngine::intersects only works with QgsAbstractGeometry
478 : : // objects and accordingly we have to use those, instead of the previously build geos
479 : : // representations available in ringEngines
480 : : // This needs addressing by extending the QgsGeometryEngine relation tests to allow testing against
481 : : // another QgsGeometryEngine object.
482 : 39 : for ( int interiorRing = static_cast< int >( i ); interiorRing < polygon.numInteriorRings(); ++interiorRing )
483 : : {
484 : 13 : if ( first->intersects( polygon.interiorRing( interiorRing ) ) )
485 : 0 : return false;
486 : 13 : }
487 : 26 : }
488 : 13 : }
489 : 13 : return true;
490 : 13 : }
491 : :
492 : :
493 : 13 : double _minimum_distance_between_coordinates( const QgsPolygon &polygon )
494 : : {
495 : 13 : double min_d = 1e20;
496 : :
497 : 13 : std::vector< const QgsLineString * > rings;
498 : 13 : rings.reserve( 1 + polygon.numInteriorRings() );
499 : 13 : rings.emplace_back( qgsgeometry_cast< const QgsLineString * >( polygon.exteriorRing() ) );
500 : 26 : for ( int i = 0; i < polygon.numInteriorRings(); ++i )
501 : 13 : rings.emplace_back( qgsgeometry_cast< const QgsLineString * >( polygon.interiorRing( i ) ) );
502 : :
503 : 39 : for ( const QgsLineString *ring : rings )
504 : : {
505 : 26 : const int numPoints = ring->numPoints();
506 : 26 : if ( numPoints <= 1 )
507 : 0 : continue;
508 : :
509 : 26 : const double *srcXData = ring->xData();
510 : 26 : const double *srcYData = ring->yData();
511 : 26 : double x0 = *srcXData++;
512 : 26 : double y0 = *srcYData++;
513 : 117 : for ( int i = 1; i < numPoints; ++i )
514 : : {
515 : 91 : double x1 = *srcXData++;
516 : 91 : double y1 = *srcYData++;
517 : 91 : double d = ( x0 - x1 ) * ( x0 - x1 ) + ( y0 - y1 ) * ( y0 - y1 );
518 : 91 : if ( d < min_d )
519 : 26 : min_d = d;
520 : 91 : x0 = x1;
521 : 91 : y0 = y1;
522 : 91 : }
523 : : }
524 : :
525 : 13 : return min_d != 1e20 ? std::sqrt( min_d ) : 1e20;
526 : 13 : }
527 : :
528 : 13 : void QgsTessellator::addPolygon( const QgsPolygon &polygon, float extrusionHeight )
529 : : {
530 : 13 : const QgsLineString *exterior = qgsgeometry_cast< const QgsLineString * >( polygon.exteriorRing() );
531 : 13 : if ( !exterior )
532 : 0 : return;
533 : :
534 : 13 : const QVector3D pNormal = !mNoZ ? _calculateNormal( exterior, mOriginX, mOriginY, mInvertNormals ) : QVector3D();
535 : 13 : const int pCount = exterior->numPoints();
536 : 13 : if ( pCount == 0 )
537 : 0 : return;
538 : :
539 : 13 : float zMin = std::numeric_limits<float>::max();
540 : 13 : float zMax = std::numeric_limits<float>::min();
541 : :
542 : 13 : const float scale = mBounds.isNull() ? 1.0 : std::max( 10000.0 / mBounds.width(), 10000.0 / mBounds.height() );
543 : :
544 : 13 : std::unique_ptr<QMatrix4x4> toNewBase, toOldBase;
545 : 13 : QgsPoint ptStart, pt0;
546 : 13 : std::unique_ptr<QgsPolygon> polygonNew;
547 : 26 : auto rotatePolygonToXYPlane = [&]()
548 : : {
549 : 13 : if ( !mNoZ && pNormal != QVector3D( 0, 0, 1 ) )
550 : : {
551 : : // this is not a horizontal plane - need to reproject the polygon to a new base so that
552 : : // we can do the triangulation in a plane
553 : 0 : QVector3D pXVector, pYVector;
554 : 0 : _normalVectorToXYVectors( pNormal, pXVector, pYVector );
555 : :
556 : : // so now we have three orthogonal unit vectors defining new base
557 : : // let's build transform matrix. We actually need just a 3x3 matrix,
558 : : // but Qt does not have good support for it, so using 4x4 matrix instead.
559 : 0 : toNewBase.reset( new QMatrix4x4(
560 : 0 : pXVector.x(), pXVector.y(), pXVector.z(), 0,
561 : 0 : pYVector.x(), pYVector.y(), pYVector.z(), 0,
562 : 0 : pNormal.x(), pNormal.y(), pNormal.z(), 0,
563 : : 0, 0, 0, 0 ) );
564 : :
565 : : // our 3x3 matrix is orthogonal, so for inverse we only need to transpose it
566 : 0 : toOldBase.reset( new QMatrix4x4( toNewBase->transposed() ) );
567 : 0 : }
568 : :
569 : 13 : ptStart = QgsPoint( exterior->startPoint() );
570 : 13 : pt0 = QgsPoint( QgsWkbTypes::PointZ, ptStart.x(), ptStart.y(), std::isnan( ptStart.z() ) ? 0 : ptStart.z() );
571 : :
572 : : // subtract ptFirst from geometry for better numerical stability in triangulation
573 : : // and apply new 3D vector base if the polygon is not horizontal
574 : :
575 : 13 : polygonNew.reset( _transform_polygon_to_new_base( polygon, pt0, toNewBase.get(), scale ) );
576 : 13 : };
577 : :
578 : 13 : if ( !mNoZ && !qgsDoubleNear( pNormal.length(), 1, 0.001 ) )
579 : 0 : return; // this should not happen - pNormal should be normalized to unit length
580 : :
581 : 13 : QVector3D upVector( 0, 0, 1 );
582 : 13 : float pNormalUpVectorDotProduct = QVector3D::dotProduct( upVector, pNormal );
583 : 13 : float radsBetwwenUpNormal = qAcos( pNormalUpVectorDotProduct );
584 : :
585 : 13 : float detectionDelta = qDegreesToRadians( 10.0f );
586 : 13 : int facade = 0;
587 : 13 : if ( radsBetwwenUpNormal > M_PI_2 - detectionDelta && radsBetwwenUpNormal < M_PI_2 + detectionDelta ) facade = 1;
588 : 0 : else if ( radsBetwwenUpNormal > - M_PI_2 - detectionDelta && radsBetwwenUpNormal < -M_PI_2 + detectionDelta ) facade = 1;
589 : 0 : else facade = 2;
590 : :
591 : 13 : if ( pCount == 4 && polygon.numInteriorRings() == 0 && ( mTessellatedFacade & facade ) )
592 : : {
593 : 0 : QgsLineString *triangle = nullptr;
594 : 0 : if ( mAddTextureCoords )
595 : : {
596 : 0 : rotatePolygonToXYPlane();
597 : 0 : triangle = qgsgeometry_cast< QgsLineString * >( polygonNew->exteriorRing() );
598 : : Q_ASSERT( polygonNew->exteriorRing()->numPoints() >= 3 );
599 : 0 : }
600 : :
601 : : // polygon is a triangle - write vertices to the output data array without triangulation
602 : 0 : const double *xData = exterior->xData();
603 : 0 : const double *yData = exterior->yData();
604 : 0 : const double *zData = !mNoZ ? exterior->zData() : nullptr;
605 : 0 : for ( int i = 0; i < 3; i++ )
606 : : {
607 : 0 : float z = !zData ? 0 : *zData;
608 : 0 : if ( z < zMin )
609 : 0 : zMin = z;
610 : 0 : if ( z > zMax )
611 : 0 : zMax = z;
612 : :
613 : 0 : mData << *xData - mOriginX << z << - *yData + mOriginY;
614 : 0 : if ( mAddNormals )
615 : 0 : mData << pNormal.x() << pNormal.z() << - pNormal.y();
616 : 0 : if ( mAddTextureCoords )
617 : : {
618 : 0 : std::pair<float, float> p( triangle->xAt( i ), triangle->yAt( i ) );
619 : 0 : if ( facade & 1 )
620 : : {
621 : 0 : p = rotateCoords( p.first, p.second, 0.0f, 0.0f, mTextureRotation );
622 : 0 : }
623 : 0 : else if ( facade & 2 )
624 : : {
625 : 0 : p = rotateCoords( p.first, p.second, 0.0f, 0.0f, mTextureRotation );
626 : 0 : }
627 : 0 : mData << p.first << p.second;
628 : 0 : }
629 : 0 : xData++; yData++;
630 : : // zData can be nullptr if mNoZ or triangle is 2D
631 : 0 : if ( zData )
632 : 0 : zData++;
633 : 0 : }
634 : :
635 : 0 : if ( mAddBackFaces )
636 : : {
637 : : // the same triangle with reversed order of coordinates and inverted normal
638 : 0 : for ( int i = 2; i >= 0; i-- )
639 : : {
640 : 0 : mData << exterior->xAt( i ) - mOriginX << ( mNoZ ? 0 : exterior->zAt( i ) ) << - exterior->yAt( i ) + mOriginY;
641 : 0 : if ( mAddNormals )
642 : 0 : mData << -pNormal.x() << -pNormal.z() << pNormal.y();
643 : 0 : if ( mAddTextureCoords )
644 : : {
645 : 0 : std::pair<float, float> p( triangle->xAt( i ), triangle->yAt( i ) );
646 : 0 : if ( facade & 1 )
647 : : {
648 : 0 : p = rotateCoords( p.first, p.second, 0.0f, 0.0f, mTextureRotation );
649 : 0 : }
650 : 0 : else if ( facade & 2 )
651 : : {
652 : 0 : p = rotateCoords( p.first, p.second, 0.0f, 0.0f, mTextureRotation );
653 : 0 : }
654 : 0 : mData << p.first << p.second;
655 : 0 : }
656 : 0 : }
657 : 0 : }
658 : 0 : }
659 : 13 : else if ( mTessellatedFacade & facade )
660 : : {
661 : :
662 : 13 : rotatePolygonToXYPlane();
663 : :
664 : 13 : if ( _minimum_distance_between_coordinates( *polygonNew ) < 0.001 )
665 : : {
666 : : // when the distances between coordinates of input points are very small,
667 : : // the triangulation likes to crash on numerical errors - when the distances are ~ 1e-5
668 : : // Assuming that the coordinates should be in a projected CRS, we should be able
669 : : // to simplify geometries that may cause problems and avoid possible crashes
670 : 0 : QgsGeometry polygonSimplified = QgsGeometry( polygonNew->clone() ).simplify( 0.001 );
671 : 0 : if ( polygonSimplified.isNull() )
672 : : {
673 : 0 : QgsMessageLog::logMessage( QObject::tr( "geometry simplification failed - skipping" ), QObject::tr( "3D" ) );
674 : 0 : return;
675 : : }
676 : 0 : const QgsPolygon *polygonSimplifiedData = qgsgeometry_cast<const QgsPolygon *>( polygonSimplified.constGet() );
677 : 0 : if ( _minimum_distance_between_coordinates( *polygonSimplifiedData ) < 0.001 )
678 : : {
679 : : // Failed to fix that. It could be a really tiny geometry... or maybe they gave us
680 : : // geometry in unprojected lat/lon coordinates
681 : 0 : QgsMessageLog::logMessage( QObject::tr( "geometry's coordinates are too close to each other and simplification failed - skipping" ), QObject::tr( "3D" ) );
682 : 0 : return;
683 : : }
684 : : else
685 : : {
686 : 0 : polygonNew.reset( polygonSimplifiedData->clone() );
687 : : }
688 : 0 : }
689 : :
690 : 13 : if ( !_check_intersecting_rings( *polygonNew ) )
691 : : {
692 : : // skip the polygon - it would cause a crash inside poly2tri library
693 : 0 : QgsMessageLog::logMessage( QObject::tr( "polygon rings self-intersect or intersect each other - skipping" ), QObject::tr( "3D" ) );
694 : 0 : return;
695 : : }
696 : :
697 : 13 : QList< std::vector<p2t::Point *> > polylinesToDelete;
698 : 13 : QHash<p2t::Point *, float> z;
699 : :
700 : : // polygon exterior
701 : 13 : std::vector<p2t::Point *> polyline;
702 : 13 : _ringToPoly2tri( qgsgeometry_cast< const QgsLineString * >( polygonNew->exteriorRing() ), polyline, mNoZ ? nullptr : &z );
703 : 13 : polylinesToDelete << polyline;
704 : :
705 : 13 : std::unique_ptr<p2t::CDT> cdt( new p2t::CDT( polyline ) );
706 : :
707 : : // polygon holes
708 : 26 : for ( int i = 0; i < polygonNew->numInteriorRings(); ++i )
709 : : {
710 : 13 : std::vector<p2t::Point *> holePolyline;
711 : 13 : const QgsLineString *hole = qgsgeometry_cast< const QgsLineString *>( polygonNew->interiorRing( i ) );
712 : :
713 : 13 : _ringToPoly2tri( hole, holePolyline, mNoZ ? nullptr : &z );
714 : :
715 : 13 : cdt->AddHole( holePolyline );
716 : 13 : polylinesToDelete << holePolyline;
717 : 13 : }
718 : :
719 : : // run triangulation and write vertices to the output data array
720 : : try
721 : : {
722 : 13 : cdt->Triangulate();
723 : :
724 : 13 : std::vector<p2t::Triangle *> triangles = cdt->GetTriangles();
725 : :
726 : 13 : mData.reserve( mData.size() + 3 * triangles.size() * ( stride() / sizeof( float ) ) );
727 : 104 : for ( size_t i = 0; i < triangles.size(); ++i )
728 : : {
729 : 91 : p2t::Triangle *t = triangles[i];
730 : 364 : for ( int j = 0; j < 3; ++j )
731 : : {
732 : 273 : p2t::Point *p = t->GetPoint( j );
733 : 273 : QVector4D pt( p->x, p->y, mNoZ ? 0 : z[p], 0 );
734 : 273 : if ( toOldBase )
735 : 0 : pt = *toOldBase * pt;
736 : 273 : const double fx = ( pt.x() / scale ) - mOriginX + pt0.x();
737 : 273 : const double fy = ( pt.y() / scale ) - mOriginY + pt0.y();
738 : 273 : const double fz = mNoZ ? 0 : ( pt.z() + extrusionHeight + pt0.z() );
739 : 273 : if ( fz < zMin )
740 : 13 : zMin = fz;
741 : 273 : if ( fz > zMax )
742 : 0 : zMax = fz;
743 : :
744 : 273 : mData << fx << fz << -fy;
745 : 273 : if ( mAddNormals )
746 : 0 : mData << pNormal.x() << pNormal.z() << - pNormal.y();
747 : 273 : if ( mAddTextureCoords )
748 : : {
749 : 0 : std::pair<float, float> pr = rotateCoords( p->x, p->y, 0.0f, 0.0f, mTextureRotation );
750 : 0 : mData << pr.first << pr.second;
751 : 0 : }
752 : 273 : }
753 : :
754 : 91 : if ( mAddBackFaces )
755 : : {
756 : : // the same triangle with reversed order of coordinates and inverted normal
757 : 0 : for ( int j = 2; j >= 0; --j )
758 : : {
759 : 0 : p2t::Point *p = t->GetPoint( j );
760 : 0 : QVector4D pt( p->x, p->y, mNoZ ? 0 : z[p], 0 );
761 : 0 : if ( toOldBase )
762 : 0 : pt = *toOldBase * pt;
763 : 0 : const double fx = ( pt.x() / scale ) - mOriginX + pt0.x();
764 : 0 : const double fy = ( pt.y() / scale ) - mOriginY + pt0.y();
765 : 0 : const double fz = mNoZ ? 0 : ( pt.z() + extrusionHeight + pt0.z() );
766 : 0 : mData << fx << fz << -fy;
767 : 0 : if ( mAddNormals )
768 : 0 : mData << -pNormal.x() << -pNormal.z() << pNormal.y();
769 : 0 : if ( mAddTextureCoords )
770 : : {
771 : 0 : std::pair<float, float> pr = rotateCoords( p->x, p->y, 0.0f, 0.0f, mTextureRotation );
772 : 0 : mData << pr.first << pr.second;
773 : 0 : }
774 : 0 : }
775 : 0 : }
776 : 91 : }
777 : 13 : }
778 : : catch ( ... )
779 : : {
780 : 0 : QgsMessageLog::logMessage( QObject::tr( "Triangulation failed. Skipping polygon…" ), QObject::tr( "3D" ) );
781 : 0 : }
782 : :
783 : 39 : for ( int i = 0; i < polylinesToDelete.count(); ++i )
784 : 26 : qDeleteAll( polylinesToDelete[i] );
785 : 13 : }
786 : :
787 : : // add walls if extrusion is enabled
788 : 13 : if ( extrusionHeight != 0 && ( mTessellatedFacade & 1 ) )
789 : : {
790 : 0 : _makeWalls( *exterior, false, extrusionHeight, mData, mAddNormals, mAddTextureCoords, mOriginX, mOriginY, mTextureRotation );
791 : :
792 : 0 : for ( int i = 0; i < polygon.numInteriorRings(); ++i )
793 : 0 : _makeWalls( *qgsgeometry_cast< const QgsLineString * >( polygon.interiorRing( i ) ), true, extrusionHeight, mData, mAddNormals, mAddTextureCoords, mOriginX, mOriginY, mTextureRotation );
794 : :
795 : 0 : zMax += extrusionHeight;
796 : 0 : }
797 : :
798 : 13 : if ( zMin < mZMin )
799 : 7 : mZMin = zMin;
800 : 13 : if ( zMax > mZMax )
801 : 0 : mZMax = zMax;
802 : 13 : }
803 : :
804 : 0 : QgsPoint getPointFromData( QVector< float >::const_iterator &it )
805 : : {
806 : : // tessellator geometry is x, z, -y
807 : 0 : double x = *it;
808 : 0 : ++it;
809 : 0 : double z = *it;
810 : 0 : ++it;
811 : 0 : double y = -( *it );
812 : 0 : ++it;
813 : 0 : return QgsPoint( x, y, z );
814 : : }
815 : :
816 : 7 : int QgsTessellator::dataVerticesCount() const
817 : : {
818 : 7 : return mData.size() / ( stride() / sizeof( float ) );
819 : : }
820 : :
821 : 0 : std::unique_ptr<QgsMultiPolygon> QgsTessellator::asMultiPolygon() const
822 : : {
823 : 0 : std::unique_ptr< QgsMultiPolygon > mp = std::make_unique< QgsMultiPolygon >();
824 : 0 : const QVector<float> data = mData;
825 : 0 : mp->reserve( mData.size() );
826 : 0 : for ( auto it = data.constBegin(); it != data.constEnd(); )
827 : : {
828 : 0 : QgsPoint p1 = getPointFromData( it );
829 : 0 : QgsPoint p2 = getPointFromData( it );
830 : 0 : QgsPoint p3 = getPointFromData( it );
831 : 0 : mp->addGeometry( new QgsTriangle( p1, p2, p3 ) );
832 : 0 : }
833 : 0 : return mp;
834 : 0 : }
|
