[tesselator] Allow internal scaling of coordinates by input coordinate

bounds to avoid numerical instability with close coordinates (e.g. calculating
tesselation of geometries in geographic CRS)

python/core/auto_generated/qgstessellator.sip.in

@@ -30,6 +30,15 @@ Optionally provides extrusion by adding triangles that serve as walls when extru
     QgsTessellator( double originX, double originY, bool addNormals, bool invertNormals = false, bool addBackFaces = false );
 %Docstring
 Creates tessellator with a specified origin point of the world (in map coordinates)
+%End
+
+    QgsTessellator( const QgsRectangle &bounds, bool addNormals, bool invertNormals = false, bool addBackFaces = false );
+%Docstring
+Creates tessellator with a specified ``bounds`` of input geometry coordinates.
+This constructor allows the tesselator to map input coordinates to a desirable range for numerically
+stability during calculations.
+
+.. versionadded:: 3.10
 %End
 
     void addPolygon( const QgsPolygon &polygon, float extrusionHeight );

src/3d/processing/qgsalgorithmtessellate.cpp

@@ -89,7 +89,7 @@ QgsFeatureList QgsTessellateAlgorithm::processFeature( const QgsFeature &feature
     else
     {
       QgsRectangle bounds = f.geometry().boundingBox();
-      QgsTessellator t( bounds.xMinimum(), bounds.yMinimum(), false );
+      QgsTessellator t( bounds, false );
 
       if ( f.geometry().isMultipart() )
       {

src/core/qgstessellator.cpp

@@ -72,13 +72,28 @@ QgsTessellator::QgsTessellator( double originX, double originY, bool addNormals,
   , mAddNormals( addNormals )
   , mInvertNormals( invertNormals )
   , mAddBackFaces( addBackFaces )
+{
+  init();
+}
+
+QgsTessellator::QgsTessellator( const QgsRectangle &bounds, bool addNormals, bool invertNormals, bool addBackFaces )
+  : mBounds( bounds )
+  , mOriginX( mBounds.xMinimum() )
+  , mOriginY( mBounds.yMinimum() )
+  , mAddNormals( addNormals )
+  , mInvertNormals( invertNormals )
+  , mAddBackFaces( addBackFaces )
+{
+  init();
+}
+
+void QgsTessellator::init()
 {
   mStride = 3 * sizeof( float );
-  if ( addNormals )
+  if ( mAddNormals )
     mStride += 3 * sizeof( float );
 }
 
-
 static bool _isRingCounterClockWise( const QgsCurve &ring )
 {
   double a = 0;
@@ -240,7 +255,7 @@ inline double _round_coord( double x )
 }
 
 
-static QgsCurve *_transform_ring_to_new_base( const QgsCurve &curve, const QgsPoint &pt0, const QMatrix4x4 *toNewBase )
+static QgsCurve *_transform_ring_to_new_base( const QgsCurve &curve, const QgsPoint &pt0, const QMatrix4x4 *toNewBase, float scaleX, float scaleY )
 {
   int count = curve.numPoints();
   QVector<QgsPoint> pts;
@@ -255,6 +270,10 @@ static QgsCurve *_transform_ring_to_new_base( const QgsCurve &curve, const QgsPo
     if ( toNewBase )
       v = toNewBase->map( v );
 
+    // scale coordinates
+    v.setX( v.x() * scaleX );
+    v.setY( v.y() * scaleY );
+
     // we also round coordinates before passing them to poly2tri triangulation in order to fix possible numerical
     // 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).
     // See TestQgsTessellator::testIssue17745().
@@ -272,12 +291,12 @@ static QgsCurve *_transform_ring_to_new_base( const QgsCurve &curve, const QgsPo
 }
 
 
-static QgsPolygon *_transform_polygon_to_new_base( const QgsPolygon &polygon, const QgsPoint &pt0, const QMatrix4x4 *toNewBase )
+static QgsPolygon *_transform_polygon_to_new_base( const QgsPolygon &polygon, const QgsPoint &pt0, const QMatrix4x4 *toNewBase, float scaleX, float scaleY )
 {
   QgsPolygon *p = new QgsPolygon;
-  p->setExteriorRing( _transform_ring_to_new_base( *polygon.exteriorRing(), pt0, toNewBase ) );
+  p->setExteriorRing( _transform_ring_to_new_base( *polygon.exteriorRing(), pt0, toNewBase, scaleX, scaleY ) );
   for ( int i = 0; i < polygon.numInteriorRings(); ++i )
-    p->addInteriorRing( _transform_ring_to_new_base( *polygon.interiorRing( i ), pt0, toNewBase ) );
+    p->addInteriorRing( _transform_ring_to_new_base( *polygon.interiorRing( i ), pt0, toNewBase, scaleX, scaleY ) );
   return p;
 }
 
@@ -412,9 +431,13 @@ void QgsTessellator::addPolygon( const QgsPolygon &polygon, float extrusionHeigh
     const QgsPoint ptStart( exterior->startPoint() );
     const QgsPoint pt0( QgsWkbTypes::PointZ, ptStart.x(), ptStart.y(), std::isnan( ptStart.z() ) ? 0 : ptStart.z() );
 
+    const float scaleX = !mBounds.isNull() ? 10000.0 / mBounds.width() : 1.0;
+    const float scaleY = !mBounds.isNull() ? 10000.0 / mBounds.height() : 1.0;
+
     // subtract ptFirst from geometry for better numerical stability in triangulation
     // and apply new 3D vector base if the polygon is not horizontal
-    std::unique_ptr<QgsPolygon> polygonNew( _transform_polygon_to_new_base( polygon, pt0, toNewBase.get() ) );
+
+    std::unique_ptr<QgsPolygon> polygonNew( _transform_polygon_to_new_base( polygon, pt0, toNewBase.get(), scaleX, scaleY ) );
 
     if ( _minimum_distance_between_coordinates( *polygonNew ) < 0.001 )
     {
@@ -487,8 +510,8 @@ void QgsTessellator::addPolygon( const QgsPolygon &polygon, float extrusionHeigh
           QVector4D pt( p->x, p->y, z[p], 0 );
           if ( toOldBase )
             pt = *toOldBase * pt;
-          const double fx = pt.x() - mOriginX + pt0.x();
-          const double fy = pt.y() - mOriginY + pt0.y();
+          const double fx = ( pt.x() / scaleX ) - mOriginX + pt0.x();
+          const double fy = ( pt.y() / scaleY ) - mOriginY + pt0.y();
           const double fz = pt.z() + extrusionHeight + pt0.z();
           mData << fx << fz << -fy;
           if ( mAddNormals )
@@ -504,8 +527,8 @@ void QgsTessellator::addPolygon( const QgsPolygon &polygon, float extrusionHeigh
             QVector4D pt( p->x, p->y, z[p], 0 );
             if ( toOldBase )
               pt = *toOldBase * pt;
-            const double fx = pt.x() - mOriginX + pt0.x();
-            const double fy = pt.y() - mOriginY + pt0.y();
+            const double fx = ( pt.x() / scaleX ) - mOriginX + pt0.x();
+            const double fy = ( pt.y() / scaleY ) - mOriginY + pt0.y();
             const double fz = pt.z() + extrusionHeight + pt0.z();
             mData << fx << fz << -fy;
             if ( mAddNormals )

src/core/qgstessellator.h

@@ -18,6 +18,7 @@
 
 #include "qgis_core.h"
 #include "qgis_sip.h"
+#include "qgsrectangle.h"
 
 class QgsPolygon;
 class QgsMultiPolygon;
@@ -42,6 +43,14 @@ class CORE_EXPORT QgsTessellator
     //! Creates tessellator with a specified origin point of the world (in map coordinates)
     QgsTessellator( double originX, double originY, bool addNormals, bool invertNormals = false, bool addBackFaces = false );
 
+    /**
+     * Creates tessellator with a specified \a bounds of input geometry coordinates.
+     * This constructor allows the tesselator to map input coordinates to a desirable range for numerically
+     * stability during calculations.
+     * \since QGIS 3.10
+     */
+    QgsTessellator( const QgsRectangle &bounds, bool addNormals, bool invertNormals = false, bool addBackFaces = false );
+
     //! Tessellates a triangle and adds its vertex entries to the output data array
     void addPolygon( const QgsPolygon &polygon, float extrusionHeight );
 
@@ -64,6 +73,9 @@ class CORE_EXPORT QgsTessellator
     std::unique_ptr< QgsMultiPolygon > asMultiPolygon() const SIP_SKIP;
 
   private:
+    void init();
+
+    QgsRectangle mBounds;
     double mOriginX = 0, mOriginY = 0;
     bool mAddNormals = false;
     bool mInvertNormals = false;

tests/src/3d/testqgstessellator.cpp

@@ -135,6 +135,7 @@ class TestQgsTessellator : public QObject
     void testIssue17745();
     void testCrashSelfIntersection();
     void testCrashEmptyPolygon();
+    void testBoundsScaling();
 
   private:
 };
@@ -348,6 +349,26 @@ void TestQgsTessellator::testCrashEmptyPolygon()
   t.addPolygon( p, 0 );  // must not crash - that's all we test here
 }
 
+void TestQgsTessellator::testBoundsScaling()
+{
+  QgsPolygon polygon;
+  polygon.fromWkt( "POLYGON((1 1, 1.00000001 1, 1.00000001 1.00000001, 1 1.0000000001, 1 1))" );
+
+  QList<TriangleCoords> tc;
+  tc << TriangleCoords( QVector3D( 0, 1e-10, 0 ), QVector3D( 1e-08, 0, 0 ), QVector3D( 1e-08, 1e-08, 0 ), QVector3D( 0, 0, 1 ), QVector3D( 0, 0, 1 ), QVector3D( 0, 0, 1 ) );
+  tc << TriangleCoords( QVector3D( 0, 1e-10, 0 ), QVector3D( 0, 0, 0 ), QVector3D( 1e-08, 0, 0 ), QVector3D( 0, 0, 1 ), QVector3D( 0, 0, 1 ), QVector3D( 0, 0, 1 ) );
+
+  // without using bounds -- numerically unstable, expect no result
+  QgsTessellator t( 0, 0, true );
+  t.addPolygon( polygon, 0 );
+  QCOMPARE( t.data().size(), 0 );
+
+  // using bounds scaling, expect good result
+  QgsTessellator t2( polygon.boundingBox(), true );
+  t2.addPolygon( polygon, 0 );
+  QVERIFY( checkTriangleOutput( t2.data(), true, tc ) );
+}
+
 
 QGSTEST_MAIN( TestQgsTessellator )
 #include "testqgstessellator.moc"