Merge pull request #5295 from liberostelios/3d-improve-triangulation

Improve triangulation for QGIS 3D
qgis · Oct 10, 2017 · 59fa541 · 59fa541
2 parents 904a877 + 8eaa8c2
commit 59fa541
Show file tree

Hide file tree

Showing 2 changed files with 158 additions and 53 deletions.
diff --git a/src/3d/qgstessellator.cpp b/src/3d/qgstessellator.cpp
@@ -19,12 +19,13 @@
 #include "qgsgeometry.h"
 #include "qgspoint.h"
 #include "qgspolygon.h"
+#include "qgis_sip.h"
 
 #include "poly2tri/poly2tri.h"
 
 #include <QtDebug>
-
 #include <QVector3D>
+#include <algorithm>
 
 static void make_quad( float x0, float y0, float x1, float y1, float zLow, float zHigh, QVector<float> &data, bool addNormals )
 {
@@ -110,6 +111,29 @@ static void _makeWalls( const QgsCurve &ring, bool ccw, float extrusionHeight, Q
   }
 }
 
+static QVector3D _calculateNormal( const QgsCurve *curve )
+{
+  // Calculate the polygon's normal vector, based on Newell's method
+  // https://www.khronos.org/opengl/wiki/Calculating_a_Surface_Normal
+  QgsVertexId::VertexType vt;
+  QgsPoint pt1, pt2;
+  QVector3D normal( 0, 0, 0 );
+
+  for ( int i = 0; i < curve->numPoints() - 1; i++ )
+  {
+    curve->pointAt( i, pt1, vt );
+    curve->pointAt( i + 1, pt2, vt );
+
+    normal.setX( normal.x() + ( pt1.y() - pt2.y() ) * ( pt1.z() + pt2.z() ) );
+    normal.setY( normal.y() + ( pt1.z() - pt2.z() ) * ( pt1.x() + pt2.x() ) );
+    normal.setZ( normal.z() + ( pt1.x() - pt2.x() ) * ( pt1.y() + pt2.y() ) );
+  }
+
+  normal.normalize();
+
+  return normal;
+}
+
 void QgsTessellator::addPolygon( const QgsPolygonV2 &polygon, float extrusionHeight )
 {
   // At this point we assume that input polygons are valid according to the OGC definition.
@@ -157,68 +181,156 @@ void QgsTessellator::addPolygon( const QgsPolygonV2 &polygon, float extrusionHei
   polyline.reserve( exterior->numPoints() );
 
   QgsVertexId::VertexType vt;
-  QgsPoint pt, ptPrev;
+  QgsPoint pt;
+
+  const QVector3D pNormal = _calculateNormal( exterior );
+  const int pCount = exterior->numPoints();
 
-  for ( int i = 0; i < exterior->numPoints() - 1; ++i )
+  // Polygon is a triangle
+  if ( pCount == 4 )
   {
-    exterior->pointAt( i, pt, vt );
-    if ( i == 0 || pt != ptPrev )
+    QgsPoint pt;
+    for ( int i = 0; i < 3; i++ )
     {
-      p2t::Point *pt2 = new p2t::Point( pt.x() - mOriginX, pt.y() - mOriginY );
-      polyline.push_back( pt2 );
-      float zPt = qIsNaN( pt.z() ) ? 0 : pt.z();
-      z[pt2] = zPt;
+      exterior->pointAt( i, pt, vt );
+      mData << pt.x() - mOriginX << pt.z() << - pt.y() + mOriginY;
+      if ( mAddNormals )
+        mData << pNormal.x() << pNormal.z() << - pNormal.y();
     }
-    ptPrev = pt;
   }
-  polylinesToDelete << polyline;
+  else
+  {
 
-  p2t::CDT *cdt = new p2t::CDT( polyline );
+    if ( !qgsDoubleNear( pNormal.length(), 1, 0.001 ) )
+    {
+      return;
+    }
 
-  // polygon holes
-  for ( int i = 0; i < polygon.numInteriorRings(); ++i )
-  {
-    std::vector<p2t::Point *> holePolyline;
-    holePolyline.reserve( exterior->numPoints() );
-    const QgsCurve *hole = polygon.interiorRing( i );
-    for ( int j = 0; j < hole->numPoints() - 1; ++j )
+    const QgsPoint ptFirst( exterior->startPoint() );
+    QVector3D pOrigin( ptFirst.x(), ptFirst.y(), ptFirst.z() ), pXVector;
+    // Here we define the two perpendicular vectors that define the local
+    // 2D space on the plane. They will act as axis for which we will
+    // calculate the projection coordinates of a 3D point to the plane.
+    if ( pNormal.z() > 0.001 || pNormal.z() < -0.001 )
+    {
+      pXVector = QVector3D( 1, 0, -pNormal.x() / pNormal.z() );
+    }
+    else if ( pNormal.y() > 0.001 || pNormal.y() < -0.001 )
     {
-      hole->pointAt( j, pt, vt );
-      if ( j == 0 || pt != ptPrev )
+      pXVector = QVector3D( 1, -pNormal.x() / pNormal.y(), 0 );
+    }
+    else
+    {
+      pXVector = QVector3D( -pNormal.y() / pNormal.x(), 1, 0 );
+    }
+    pXVector.normalize();
+    const QVector3D pYVector = QVector3D::normal( pNormal, pXVector );
+
+    for ( int i = 0; i < pCount - 1; ++i )
+    {
+      exterior->pointAt( i, pt, vt );
+      QVector3D tempPt( pt.x(), pt.y(), ( qIsNaN( pt.z() ) ? 0 : pt.z() ) );
+      const float x = QVector3D::dotProduct( tempPt - pOrigin, pXVector );
+      const float y = QVector3D::dotProduct( tempPt - pOrigin, pYVector );
+
+      const bool found = std::find_if( polyline.begin(), polyline.end(), [x, y]( p2t::Point *&p ) { return *p == p2t::Point( x, y ); } ) != polyline.end();
+
+      if ( found )
       {
-        p2t::Point *pt2 = new p2t::Point( pt.x() - mOriginX, pt.y() - mOriginY );
-        holePolyline.push_back( pt2 );
-        float zPt = qIsNaN( pt.z() ) ? 0 : pt.z();
-        z[pt2] = zPt;
+        continue;
       }
-      ptPrev = pt;
-    }
-    cdt->AddHole( holePolyline );
-    polylinesToDelete << holePolyline;
-  }
 
-  // TODO: robustness (no nearly duplicate points, invalid geometries ...)
+      p2t::Point *pt2 = new p2t::Point( x, y );
+      polyline.push_back( pt2 );
 
-  cdt->Triangulate();
+      z[pt2] = qIsNaN( pt.z() ) ? 0 : pt.z();
+    }
+    polylinesToDelete << polyline;
 
-  std::vector<p2t::Triangle *> triangles = cdt->GetTriangles();
+    // TODO: robustness (no nearly duplicate points, invalid geometries ...)
 
-  for ( size_t i = 0; i < triangles.size(); ++i )
-  {
-    p2t::Triangle *t = triangles[i];
-    for ( int j = 0; j < 3; ++j )
+    if ( polyline.size() == 3 && polygon.numInteriorRings() == 0 )
     {
-      p2t::Point *p = t->GetPoint( j );
-      float zPt = z[p];
-      mData << p->x << extrusionHeight + zPt << -p->y;
-      if ( mAddNormals )
-        mData << 0.f << 1.f << 0.f;
+      for ( std::vector<p2t::Point *>::iterator it = polyline.begin(); it != polyline.end(); it++ )
+      {
+        p2t::Point *p = *it;
+        const double zPt = z[p];
+        QVector3D nPoint = pOrigin + pXVector * p->x + pYVector * p->y;
+        const double fx = nPoint.x() - mOriginX;
+        const double fy = nPoint.y() - mOriginY;
+        const double fz = extrusionHeight + ( qIsNaN( zPt ) ? 0 : zPt );
+        mData << fx << fz << -fy;
+        if ( mAddNormals )
+          mData << pNormal.x() << pNormal.z() << - pNormal.y();
+      }
     }
-  }
+    else if ( polyline.size() >= 3 )
+    {
+      p2t::CDT *cdt = new p2t::CDT( polyline );
+
+      // polygon holes
+      for ( int i = 0; i < polygon.numInteriorRings(); ++i )
+      {
+        std::vector<p2t::Point *> holePolyline;
+        holePolyline.reserve( exterior->numPoints() );
+        const QgsCurve *hole = polygon.interiorRing( i );
+        for ( int j = 0; j < hole->numPoints() - 1; ++j )
+        {
+          hole->pointAt( j, pt, vt );
+          QVector3D tempPt( pt.x(), pt.y(), ( qIsNaN( pt.z() ) ? 0 : pt.z() ) );
+
+          const float x = QVector3D::dotProduct( tempPt - pOrigin, pXVector );
+          const float y = QVector3D::dotProduct( tempPt - pOrigin, pYVector );
+
+          const bool found = std::find_if( polyline.begin(), polyline.end(), [x, y]( p2t::Point *&p ) { return *p == p2t::Point( x, y ); } ) != polyline.end();
+
+          if ( found )
+          {
+            continue;
+          }
 
-  delete cdt;
-  for ( int i = 0; i < polylinesToDelete.count(); ++i )
-    qDeleteAll( polylinesToDelete[i] );
+          p2t::Point *pt2 = new p2t::Point( x, y );
+          holePolyline.push_back( pt2 );
+
+          z[pt2] = qIsNaN( pt.z() ) ? 0 : pt.z();
+        }
+        cdt->AddHole( holePolyline );
+        polylinesToDelete << holePolyline;
+      }
+      try
+      {
+        cdt->Triangulate();
+
+        std::vector<p2t::Triangle *> triangles = cdt->GetTriangles();
+
+        for ( size_t i = 0; i < triangles.size(); ++i )
+        {
+          p2t::Triangle *t = triangles[i];
+          for ( int j = 0; j < 3; ++j )
+          {
+            p2t::Point *p = t->GetPoint( j );
+            float zPt = z[p];
+            QVector3D nPoint = pOrigin + pXVector * p->x + pYVector * p->y;
+            float fx = nPoint.x() - mOriginX;
+            float fy = nPoint.y() - mOriginY;
+            float fz = extrusionHeight + ( qIsNaN( zPt ) ? 0 : zPt );
+            mData << fx << fz << -fy;
+            if ( mAddNormals )
+              mData << pNormal.x() << pNormal.z() << - pNormal.y();
+          }
+        }
+      }
+      catch ( ... )
+      {
+        qDebug() << "Triangulation failed. Skipping polygon...";
+      }
+
+      delete cdt;
+    }
+
+    for ( int i = 0; i < polylinesToDelete.count(); ++i )
+      qDeleteAll( polylinesToDelete[i] );
+  }
 
   // add walls if extrusion is enabled
   if ( extrusionHeight != 0 )

diff --git a/src/3d/symbols/qgspolygon3dsymbol_p.cpp b/src/3d/symbols/qgspolygon3dsymbol_p.cpp
@@ -119,13 +119,6 @@ Qt3DRender::QGeometryRenderer *QgsPolygon3DSymbolEntityNode::renderer( const Qgs
     if ( QgsWkbTypes::isCurvedType( geom.geometry()->wkbType() ) )
       geom = QgsGeometry( geom.geometry()->segmentize() );
 
-    if ( !geom.isGeosValid() )
-    {
-      // invalid geometries break tessellation
-      qDebug() << "skipping invalid geometry" << f.id();
-      continue;
-    }
-
     const QgsAbstractGeometry *g = geom.geometry();
 
     if ( const QgsPolygonV2 *poly = qgsgeometry_cast< const QgsPolygonV2 *>( g ) )