QGIS API Documentation  3.0.2-Girona (307d082)
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 
27 #include "poly2tri.h"
28 
29 #include <QtDebug>
30 #include <QMatrix4x4>
31 #include <QVector3D>
32 #include <algorithm>
33 
34 
35 static void make_quad( float x0, float y0, float z0, float x1, float y1, float z1, float height, QVector<float> &data, bool addNormals )
36 {
37  float dx = x1 - x0;
38  float dy = -( y1 - y0 );
39 
40  // perpendicular vector in plane to [x,y] is [-y,x]
41  QVector3D vn( -dy, 0, dx );
42  vn = -vn;
43  vn.normalize();
44 
45  // triangle 1
46  data << x0 << z0 + height << -y0;
47  if ( addNormals )
48  data << vn.x() << vn.y() << vn.z();
49  data << x1 << z1 + height << -y1;
50  if ( addNormals )
51  data << vn.x() << vn.y() << vn.z();
52  data << x0 << z0 << -y0;
53  if ( addNormals )
54  data << vn.x() << vn.y() << vn.z();
55 
56  // triangle 2
57  data << x0 << z0 << -y0;
58  if ( addNormals )
59  data << vn.x() << vn.y() << vn.z();
60  data << x1 << z1 + height << -y1;
61  if ( addNormals )
62  data << vn.x() << vn.y() << vn.z();
63  data << x1 << z1 << -y1;
64  if ( addNormals )
65  data << vn.x() << vn.y() << vn.z();
66 }
67 
68 
69 QgsTessellator::QgsTessellator( double originX, double originY, bool addNormals, bool invertNormals )
70  : mOriginX( originX )
71  , mOriginY( originY )
72  , mAddNormals( addNormals )
73  , mInvertNormals( invertNormals )
74 {
75  mStride = 3 * sizeof( float );
76  if ( addNormals )
77  mStride += 3 * sizeof( float );
78 }
79 
80 
81 static bool _isRingCounterClockWise( const QgsCurve &ring )
82 {
83  double a = 0;
84  int count = ring.numPoints();
86  QgsPoint pt, ptPrev;
87  ring.pointAt( 0, ptPrev, vt );
88  for ( int i = 1; i < count + 1; ++i )
89  {
90  ring.pointAt( i % count, pt, vt );
91  a += ptPrev.x() * pt.y() - ptPrev.y() * pt.x();
92  ptPrev = pt;
93  }
94  return a > 0; // clockwise if a is negative
95 }
96 
97 static void _makeWalls( const QgsCurve &ring, bool ccw, float extrusionHeight, QVector<float> &data, bool addNormals, double originX, double originY )
98 {
99  // we need to find out orientation of the ring so that the triangles we generate
100  // face the right direction
101  // (for exterior we want clockwise order, for holes we want counter-clockwise order)
102  bool is_counter_clockwise = _isRingCounterClockWise( ring );
103 
105  QgsPoint pt;
106 
107  QgsPoint ptPrev;
108  ring.pointAt( is_counter_clockwise == ccw ? 0 : ring.numPoints() - 1, ptPrev, vt );
109  for ( int i = 1; i < ring.numPoints(); ++i )
110  {
111  ring.pointAt( is_counter_clockwise == ccw ? i : ring.numPoints() - i - 1, pt, vt );
112  float x0 = ptPrev.x() - originX, y0 = ptPrev.y() - originY;
113  float x1 = pt.x() - originX, y1 = pt.y() - originY;
114  float z0 = std::isnan( ptPrev.z() ) ? 0 : ptPrev.z();
115  float z1 = std::isnan( pt.z() ) ? 0 : pt.z();
116 
117  // make a quad
118  make_quad( x0, y0, z0, x1, y1, z1, extrusionHeight, data, addNormals );
119  ptPrev = pt;
120  }
121 }
122 
123 static QVector3D _calculateNormal( const QgsCurve *curve, double originX, double originY, bool invertNormal )
124 {
126  QgsPoint pt1, pt2;
127 
128  // if it is just plain 2D curve there is no need to calculate anything
129  // because it will be a flat horizontally oriented patch
130  if ( !QgsWkbTypes::hasZ( curve->wkbType() ) )
131  return QVector3D( 0, 0, 1 );
132 
133  // often we have 3D coordinates, but Z is the same for all vertices
134  // so in order to save calculation and avoid possible issues with order of vertices
135  // (the calculation below may decide that a polygon faces downwards)
136  bool sameZ = true;
137  curve->pointAt( 0, pt1, vt );
138  for ( int i = 1; i < curve->numPoints(); i++ )
139  {
140  curve->pointAt( i, pt2, vt );
141  if ( pt1.z() != pt2.z() )
142  {
143  sameZ = false;
144  break;
145  }
146  }
147  if ( sameZ )
148  return QVector3D( 0, 0, 1 );
149 
150  // Calculate the polygon's normal vector, based on Newell's method
151  // https://www.khronos.org/opengl/wiki/Calculating_a_Surface_Normal
152  //
153  // Order of vertices is important here as it determines the front/back face of the polygon
154 
155  double nx = 0, ny = 0, nz = 0;
156  for ( int i = 0; i < curve->numPoints() - 1; i++ )
157  {
158  curve->pointAt( i, pt1, vt );
159  curve->pointAt( i + 1, pt2, vt );
160 
161  // shift points by the tessellator's origin - this does not affect normal calculation and it may save us from losing some precision
162  pt1.setX( pt1.x() - originX );
163  pt1.setY( pt1.y() - originY );
164  pt2.setX( pt2.x() - originX );
165  pt2.setY( pt2.y() - originY );
166 
167  if ( std::isnan( pt1.z() ) || std::isnan( pt2.z() ) )
168  continue;
169 
170  nx += ( pt1.y() - pt2.y() ) * ( pt1.z() + pt2.z() );
171  ny += ( pt1.z() - pt2.z() ) * ( pt1.x() + pt2.x() );
172  nz += ( pt1.x() - pt2.x() ) * ( pt1.y() + pt2.y() );
173  }
174 
175  QVector3D normal( nx, ny, nz );
176  if ( invertNormal )
177  normal = -normal;
178  normal.normalize();
179  return normal;
180 }
181 
182 
183 static void _normalVectorToXYVectors( const QVector3D &pNormal, QVector3D &pXVector, QVector3D &pYVector )
184 {
185  // Here we define the two perpendicular vectors that define the local
186  // 2D space on the plane. They will act as axis for which we will
187  // calculate the projection coordinates of a 3D point to the plane.
188  if ( pNormal.z() > 0.001 || pNormal.z() < -0.001 )
189  {
190  pXVector = QVector3D( 1, 0, -pNormal.x() / pNormal.z() );
191  }
192  else if ( pNormal.y() > 0.001 || pNormal.y() < -0.001 )
193  {
194  pXVector = QVector3D( 1, -pNormal.x() / pNormal.y(), 0 );
195  }
196  else
197  {
198  pXVector = QVector3D( -pNormal.y() / pNormal.x(), 1, 0 );
199  }
200  pXVector.normalize();
201  pYVector = QVector3D::normal( pNormal, pXVector );
202 }
203 
204 
205 static void _ringToPoly2tri( const QgsCurve *ring, std::vector<p2t::Point *> &polyline, QHash<p2t::Point *, float> &zHash )
206 {
208  QgsPoint pt;
209 
210  const int pCount = ring->numPoints();
211 
212  polyline.reserve( pCount );
213 
214  for ( int i = 0; i < pCount - 1; ++i )
215  {
216  ring->pointAt( i, pt, vt );
217  const float x = pt.x();
218  const float y = pt.y();
219  const float z = pt.z();
220 
221  const bool found = std::find_if( polyline.begin(), polyline.end(), [x, y]( p2t::Point *&p ) { return *p == p2t::Point( x, y ); } ) != polyline.end();
222 
223  if ( found )
224  {
225  continue;
226  }
227 
228  p2t::Point *pt2 = new p2t::Point( x, y );
229  polyline.push_back( pt2 );
230  zHash[pt2] = z;
231  }
232 }
233 
234 
235 inline double _round_coord( double x )
236 {
237  const double exp = 1e10; // round to 10 decimal digits
238  return round( x * exp ) / exp;
239 }
240 
241 
242 static QgsCurve *_transform_ring_to_new_base( const QgsCurve &curve, const QgsPoint &pt0, const QMatrix4x4 *toNewBase )
243 {
244  int count = curve.numPoints();
245  QVector<QgsPoint> pts;
246  pts.reserve( count );
248  for ( int i = 0; i < count; ++i )
249  {
250  QgsPoint pt;
251  curve.pointAt( i, pt, vt );
252  QgsPoint pt2( QgsWkbTypes::PointZ, pt.x() - pt0.x(), pt.y() - pt0.y(), std::isnan( pt.z() ) ? 0 : pt.z() - pt0.z() );
253  QVector4D v( pt2.x(), pt2.y(), pt2.z(), 0 );
254  if ( toNewBase )
255  v = toNewBase->map( v );
256 
257  // we also round coordinates before passing them to poly2tri triangulation in order to fix possible numerical
258  // 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).
259  // See TestQgsTessellator::testIssue17745().
260  //
261  // A hint for a similar issue: https://github.com/greenm01/poly2tri/issues/99
262  //
263  // The collinear tests uses epsilon 1e-12. Seems rounding to 12 places you still
264  // can get problems with this test when points are pretty much on a straight line.
265  // I suggest you round to 10 decimals for stability and you can live with that
266  // precision.
267 
268  pts << QgsPoint( QgsWkbTypes::PointZ, _round_coord( v.x() ), _round_coord( v.y() ), _round_coord( v.z() ) );
269  }
270  return new QgsLineString( pts );
271 }
272 
273 
274 static QgsPolygon *_transform_polygon_to_new_base( const QgsPolygon &polygon, const QgsPoint &pt0, const QMatrix4x4 *toNewBase )
275 {
276  QgsPolygon *p = new QgsPolygon;
277  p->setExteriorRing( _transform_ring_to_new_base( *polygon.exteriorRing(), pt0, toNewBase ) );
278  for ( int i = 0; i < polygon.numInteriorRings(); ++i )
279  p->addInteriorRing( _transform_ring_to_new_base( *polygon.interiorRing( i ), pt0, toNewBase ) );
280  return p;
281 }
282 
283 static bool _check_intersecting_rings( const QgsPolygon &polygon )
284 {
285  // At this point we assume that input polygons are valid according to the OGC definition.
286  // This means e.g. no duplicate points, polygons are simple (no butterfly shaped polygon with self-intersection),
287  // internal rings are inside exterior rings, rings do not cross each other, no dangles.
288 
289  // There is however an issue with polygons where rings touch:
290  // +---+
291  // | |
292  // | +-+-+
293  // | | | |
294  // | +-+ |
295  // | |
296  // +-----+
297  // This is a valid polygon with one exterior and one interior ring that touch at one point,
298  // but poly2tri library does not allow interior rings touch each other or exterior ring.
299  // TODO: Handle the situation better - rather than just detecting the problem, try to fix
300  // it by converting touching rings into one ring.
301 
302  if ( polygon.numInteriorRings() > 0 )
303  {
304  QList<QgsGeometry> geomRings;
305  geomRings << QgsGeometry( polygon.exteriorRing()->clone() );
306  for ( int i = 0; i < polygon.numInteriorRings(); ++i )
307  geomRings << QgsGeometry( polygon.interiorRing( i )->clone() );
308 
309  for ( int i = 0; i < geomRings.count(); ++i )
310  for ( int j = i + 1; j < geomRings.count(); ++j )
311  {
312  if ( geomRings[i].intersects( geomRings[j] ) )
313  return false;
314  }
315  }
316  return true;
317 }
318 
319 
321 {
322  double min_d = 1e20;
323  auto it = polygon.vertices_begin();
324 
325  if ( it == polygon.vertices_end() )
326  return min_d;
327 
328  QgsPoint p0 = *it;
329  ++it;
330  for ( ; it != polygon.vertices_end(); ++it )
331  {
332  QgsPoint p1 = *it;
333  double d = p0.distance( p1 );
334  if ( d < min_d )
335  min_d = d;
336  p0 = p1;
337  }
338  return min_d;
339 }
340 
341 
342 void QgsTessellator::addPolygon( const QgsPolygon &polygon, float extrusionHeight )
343 {
344  const QgsCurve *exterior = polygon.exteriorRing();
345 
346  const QVector3D pNormal = _calculateNormal( exterior, mOriginX, mOriginY, mInvertNormals );
347  const int pCount = exterior->numPoints();
348 
349  if ( pCount == 4 && polygon.numInteriorRings() == 0 )
350  {
351  // polygon is a triangle - write vertices to the output data array without triangulation
352  QgsPoint pt;
354  for ( int i = 0; i < 3; i++ )
355  {
356  exterior->pointAt( i, pt, vt );
357  mData << pt.x() - mOriginX << pt.z() << - pt.y() + mOriginY;
358  if ( mAddNormals )
359  mData << pNormal.x() << pNormal.z() << - pNormal.y();
360  }
361  }
362  else
363  {
364  if ( !qgsDoubleNear( pNormal.length(), 1, 0.001 ) )
365  return; // this should not happen - pNormal should be normalized to unit length
366 
367  std::unique_ptr<QMatrix4x4> toNewBase, toOldBase;
368  if ( pNormal != QVector3D( 0, 0, 1 ) )
369  {
370  // this is not a horizontal plane - need to reproject the polygon to a new base so that
371  // we can do the triangulation in a plane
372 
373  QVector3D pXVector, pYVector;
374  _normalVectorToXYVectors( pNormal, pXVector, pYVector );
375 
376  // so now we have three orthogonal unit vectors defining new base
377  // let's build transform matrix. We actually need just a 3x3 matrix,
378  // but Qt does not have good support for it, so using 4x4 matrix instead.
379  toNewBase.reset( new QMatrix4x4(
380  pXVector.x(), pXVector.y(), pXVector.z(), 0,
381  pYVector.x(), pYVector.y(), pYVector.z(), 0,
382  pNormal.x(), pNormal.y(), pNormal.z(), 0,
383  0, 0, 0, 0 ) );
384 
385  // our 3x3 matrix is orthogonal, so for inverse we only need to transpose it
386  toOldBase.reset( new QMatrix4x4( toNewBase->transposed() ) );
387  }
388 
389  const QgsPoint ptStart( exterior->startPoint() );
390  const QgsPoint pt0( QgsWkbTypes::PointZ, ptStart.x(), ptStart.y(), std::isnan( ptStart.z() ) ? 0 : ptStart.z() );
391 
392  // subtract ptFirst from geometry for better numerical stability in triangulation
393  // and apply new 3D vector base if the polygon is not horizontal
394  std::unique_ptr<QgsPolygon> polygonNew( _transform_polygon_to_new_base( polygon, pt0, toNewBase.get() ) );
395 
396  if ( _minimum_distance_between_coordinates( *polygonNew ) < 0.001 )
397  {
398  // when the distances between coordinates of input points are very small,
399  // the triangulation likes to crash on numerical errors - when the distances are ~ 1e-5
400  // Assuming that the coordinates should be in a projected CRS, we should be able
401  // to simplify geometries that may cause problems and avoid possible crashes
402  QgsGeometry polygonSimplified = QgsGeometry( polygonNew->clone() ).simplify( 0.001 );
403  const QgsPolygon *polygonSimplifiedData = qgsgeometry_cast<const QgsPolygon *>( polygonSimplified.constGet() );
404  if ( _minimum_distance_between_coordinates( *polygonSimplifiedData ) < 0.001 )
405  {
406  // Failed to fix that. It could be a really tiny geometry... or maybe they gave us
407  // geometry in unprojected lat/lon coordinates
408  QgsMessageLog::logMessage( QObject::tr( "geometry's coordinates are too close to each other and simplification failed - skipping" ), QObject::tr( "3D" ) );
409  return;
410  }
411  else
412  {
413  polygonNew.reset( polygonSimplifiedData->clone() );
414  }
415  }
416 
417  if ( !_check_intersecting_rings( *polygonNew.get() ) )
418  {
419  // skip the polygon - it would cause a crash inside poly2tri library
420  QgsMessageLog::logMessage( QObject::tr( "polygon rings intersect each other - skipping" ), QObject::tr( "3D" ) );
421  return;
422  }
423 
424  QList< std::vector<p2t::Point *> > polylinesToDelete;
425  QHash<p2t::Point *, float> z;
426 
427  // polygon exterior
428  std::vector<p2t::Point *> polyline;
429  _ringToPoly2tri( polygonNew->exteriorRing(), polyline, z );
430  polylinesToDelete << polyline;
431 
432  std::unique_ptr<p2t::CDT> cdt( new p2t::CDT( polyline ) );
433 
434  // polygon holes
435  for ( int i = 0; i < polygonNew->numInteriorRings(); ++i )
436  {
437  std::vector<p2t::Point *> holePolyline;
438  const QgsCurve *hole = polygonNew->interiorRing( i );
439 
440  _ringToPoly2tri( hole, holePolyline, z );
441 
442  cdt->AddHole( holePolyline );
443  polylinesToDelete << holePolyline;
444  }
445 
446  // run triangulation and write vertices to the output data array
447  try
448  {
449  cdt->Triangulate();
450 
451  std::vector<p2t::Triangle *> triangles = cdt->GetTriangles();
452 
453  for ( size_t i = 0; i < triangles.size(); ++i )
454  {
455  p2t::Triangle *t = triangles[i];
456  for ( int j = 0; j < 3; ++j )
457  {
458  p2t::Point *p = t->GetPoint( j );
459  QVector4D pt( p->x, p->y, z[p], 0 );
460  if ( toOldBase )
461  pt = *toOldBase * pt;
462  const double fx = pt.x() - mOriginX + pt0.x();
463  const double fy = pt.y() - mOriginY + pt0.y();
464  const double fz = pt.z() + extrusionHeight + pt0.z();
465  mData << fx << fz << -fy;
466  if ( mAddNormals )
467  mData << pNormal.x() << pNormal.z() << - pNormal.y();
468  }
469  }
470  }
471  catch ( ... )
472  {
473  QgsMessageLog::logMessage( QObject::tr( "Triangulation failed. Skipping polygon…" ), QObject::tr( "3D" ) );
474  }
475 
476  for ( int i = 0; i < polylinesToDelete.count(); ++i )
477  qDeleteAll( polylinesToDelete[i] );
478  }
479 
480  // add walls if extrusion is enabled
481  if ( extrusionHeight != 0 )
482  {
483  _makeWalls( *exterior, false, extrusionHeight, mData, mAddNormals, mOriginX, mOriginY );
484 
485  for ( int i = 0; i < polygon.numInteriorRings(); ++i )
486  _makeWalls( *polygon.interiorRing( i ), true, extrusionHeight, mData, mAddNormals, mOriginX, mOriginY );
487  }
488 }
489 
490 QgsPoint getPointFromData( QVector< float >::const_iterator &it )
491 {
492  // tessellator geometry is x, z, -y
493  double x = *it;
494  ++it;
495  double z = *it;
496  ++it;
497  double y = -( *it );
498  ++it;
499  return QgsPoint( x, y, z );
500 }
501 
502 std::unique_ptr<QgsMultiPolygon> QgsTessellator::asMultiPolygon() const
503 {
504  std::unique_ptr< QgsMultiPolygon > mp = qgis::make_unique< QgsMultiPolygon >();
505  const QVector<float> data = mData;
506  for ( auto it = data.constBegin(); it != data.constEnd(); )
507  {
508  QgsPoint p1 = getPointFromData( it );
509  QgsPoint p2 = getPointFromData( it );
510  QgsPoint p3 = getPointFromData( it );
511  mp->addGeometry( new QgsTriangle( p1, p2, p3 ) );
512  }
513  return mp;
514 }
const QgsCurve * interiorRing(int i) const
QgsTessellator(double originX, double originY, bool addNormals, bool invertNormals=false)
Creates tessellator with a specified origin point of the world (in map coordinates) ...
double y
Definition: qgspoint.h:42
double distance(double x, double y) const
Returns the distance between this point and a specified x, y coordinate.
Definition: qgspoint.cpp:587
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:111
virtual bool pointAt(int node, QgsPoint &point, QgsVertexId::VertexType &type) const =0
Returns the point and vertex id of a point within the curve.
std::unique_ptr< QgsMultiPolygon > asMultiPolygon() const
Returns the triangulation as a multipolygon geometry.
QgsPoint getPointFromData(QVector< float >::const_iterator &it)
Triangle geometry type.
Definition: qgstriangle.h:33
static bool hasZ(Type type)
Tests whether a WKB type contains the z-dimension.
Definition: qgswkbtypes.h:768
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Definition: qgis.h:251
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::Warning)
add a message to the instance (and create it if necessary)
void addInteriorRing(QgsCurve *ring) override
Adds an interior ring to the geometry (takes ownership)
Definition: qgspolygon.cpp:148
int numInteriorRings() const
double _round_coord(double x)
void addPolygon(const QgsPolygon &polygon, float extrusionHeight)
Tessellates a triangle and adds its vertex entries to the output data array.
vertex_iterator vertices_end() const
Returns STL-style iterator pointing to the imaginary vertex after the last vertex of the geometry...
T qgsgeometry_cast(const QgsAbstractGeometry *geom)
Abstract base class for curved geometry type.
Definition: qgscurve.h:35
const QgsCurve * exteriorRing() const
QgsWkbTypes::Type wkbType() const
Returns the WKB type of the geometry.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
double _minimum_distance_between_coordinates(const QgsPolygon &polygon)
void setX(double x)
Sets the point&#39;s x-coordinate.
Definition: qgspoint.h:192
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
void setY(double y)
Sets the point&#39;s y-coordinate.
Definition: qgspoint.h:203
void setExteriorRing(QgsCurve *ring) override
Sets the exterior ring of the polygon.
Definition: qgspolygon.cpp:179
QgsCurve * clone() const override=0
Clones the geometry by performing a deep copy.
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:41
QgsPolygon * clone() const override
Clones the geometry by performing a deep copy.
Definition: qgspolygon.cpp:42
double z
Definition: qgspoint.h:43
vertex_iterator vertices_begin() const
Returns STL-style iterator pointing to the first vertex of the geometry.
virtual QgsPoint startPoint() const =0
Returns the starting point of the curve.
Polygon geometry type.
Definition: qgspolygon.h:31
virtual int numPoints() const =0
Returns the number of points in the curve.
double x
Definition: qgspoint.h:41