Refactor curveToLine to emit equidistant segments and fix some issues

Fixes #16717
Fixes #16722

Include tests

src/core/geometry/qgscircularstring.cpp

@@ -23,6 +23,7 @@
 #include "qgsmaptopixel.h"
 #include "qgspoint.h"
 #include "qgswkbptr.h"
+#include "qgslogger.h"
 #include <QPainter>
 #include <QPainterPath>
 
@@ -353,6 +354,7 @@ QgsLineString *QgsCircularString::curveToLine( double tolerance, SegmentationTol
   QgsPointSequence points;
   int nPoints = numPoints();
 
+  QgsDebugMsg( QString( "curveToLine input: %1" ) .arg( asWkt( 2 ) ) );
   for ( int i = 0; i < ( nPoints - 2 ) ; i += 2 )
   {
     QgsGeometryUtils::segmentizeArc( pointN( i ), pointN( i + 1 ), pointN( i + 2 ), points, tolerance, toleranceType, is3D(), isMeasure() );

src/core/geometry/qgsgeometryutils.cpp

@@ -20,6 +20,7 @@ email                : marco.hugentobler at sourcepole dot com
 #include "qgsgeometrycollection.h"
 #include "qgslinestring.h"
 #include "qgswkbptr.h"
+#include "qgslogger.h"
 
 #include <memory>
 #include <QStringList>
@@ -631,25 +632,44 @@ double QgsGeometryUtils::circleTangentDirection( const QgsPoint &tangentPoint, c
 
 void QgsGeometryUtils::segmentizeArc( const QgsPoint &p1, const QgsPoint &p2, const QgsPoint &p3, QgsPointSequence &points, double tolerance, QgsAbstractGeometry::SegmentationToleranceType toleranceType, bool hasZ, bool hasM )
 {
+  bool reversed = false;
   bool clockwise = false;
   int segSide = segmentSide( p1, p3, p2 );
   if ( segSide == -1 )
   {
     clockwise = true;
   }
 
-  QgsPoint circlePoint1 = clockwise ? p3 : p1;
-  QgsPoint circlePoint2 = p2;
-  QgsPoint circlePoint3 = clockwise ? p1 : p3 ;
+  QgsPoint circlePoint1;
+  const QgsPoint circlePoint2 = p2;
+  QgsPoint circlePoint3;
+
+  if ( clockwise )
+  {
+    // Reverse !
+    circlePoint1 = p3;
+    circlePoint3 = p1;
+    clockwise = false;
+    reversed = true;
+  }
+  else
+  {
+    circlePoint1 = p1;
+    circlePoint3 = p3;
+  }
 
   //adapted code from PostGIS
   double radius = 0;
   double centerX = 0;
   double centerY = 0;
   circleCenterRadius( circlePoint1, circlePoint2, circlePoint3, radius, centerX, centerY );
 
+  QgsDebugMsg( QString( "Center: POINT(%1 %2) - Radius: %3 - Clockwise: %4" )
+               .arg( centerX ) .arg( centerY ) .arg( radius ) .arg( clockwise ) );
+
   if ( circlePoint1 != circlePoint3 && ( radius < 0 || qgsDoubleNear( segSide, 0.0 ) ) ) //points are colinear
   {
+    QgsDebugMsg( QString( "Collinear curve" ) );
     points.append( p1 );
     points.append( p2 );
     points.append( p3 );
@@ -662,24 +682,65 @@ void QgsGeometryUtils::segmentizeArc( const QgsPoint &p1, const QgsPoint &p2, co
     double halfAngle = std::acos( -tolerance / radius + 1 );
     increment = 2 * halfAngle;
   }
+  QgsDebugMsg( QString( "Increment: %1" ).arg( increment ) );
 
   //angles of pt1, pt2, pt3
   double a1 = std::atan2( circlePoint1.y() - centerY, circlePoint1.x() - centerX );
   double a2 = std::atan2( circlePoint2.y() - centerY, circlePoint2.x() - centerX );
   double a3 = std::atan2( circlePoint3.y() - centerY, circlePoint3.x() - centerX );
 
-  /* Adjust a3 up so we can increment from a1 to a3 cleanly */
-  if ( a3 <= a1 )
-    a3 += 2.0 * M_PI;
-  if ( a2 < a1 )
-    a2 += 2.0 * M_PI;
+  QgsDebugMsg( QString( "a1:%1 (%4) a2:%2 (%5) a3:%3 (%6)" )
+               .arg( a1 ).arg( a2 ).arg( a3 )
+               .arg( a1 * 180 / M_PI ).arg( a2 * 180 / M_PI ).arg( a3 * 180 / M_PI )
+             );
+
+  // Make segmentation symmetric
+  const bool symmetric = true;
+  if ( symmetric )
+  {
+    double angle = clockwise ? a1 - a3 : a3 - a1;
+    if ( angle < 0 ) angle += M_PI * 2;
+    QgsDebugMsg( QString( "total angle: %1 (%2)" )
+                 .arg( angle ) .arg( angle * 180 / M_PI )
+               );
+
+    /* Number of segments in output */
+    int segs = ceil( angle / increment );
+    /* Tweak increment to be regular for all the arc */
+    increment = angle / segs;
+
+    QgsDebugMsg( QString( "symmetric adjusted increment:%1" ) .arg( increment ) );
+  }
+
+  if ( clockwise )
+  {
+    increment *= -1;
+    /* Adjust a3 down so we can increment from a1 to a3 cleanly */
+    if ( a3 > a1 )
+      a3 -= 2.0 * M_PI;
+    if ( a2 > a1 )
+      a2 -= 2.0 * M_PI;
+  }
+  else
+  {
+    /* Adjust a3 up so we can increment from a1 to a3 cleanly */
+    if ( a3 < a1 )
+      a3 += 2.0 * M_PI;
+    if ( a2 < a1 )
+      a2 += 2.0 * M_PI;
+  }
+
+  QgsDebugMsg( QString( "ADJUSTED - a1:%1 (%4) a2:%2 (%5) a3:%3 (%6)" )
+               .arg( a1 ).arg( a2 ).arg( a3 )
+               .arg( a1 * 180 / M_PI ).arg( a2 * 180 / M_PI ).arg( a3 * 180 / M_PI )
+             );
 
   double x, y;
   double z = 0;
   double m = 0;
 
   QList<QgsPoint> stringPoints;
-  stringPoints.insert( clockwise ? 0 : stringPoints.size(), circlePoint1 );
+  stringPoints.insert( 0, circlePoint1 );
   if ( circlePoint2 != circlePoint3 && circlePoint1 != circlePoint2 ) //draw straight line segment if two points have the same position
   {
     QgsWkbTypes::Type pointWkbType = QgsWkbTypes::Point;
@@ -688,30 +749,58 @@ void QgsGeometryUtils::segmentizeArc( const QgsPoint &p1, const QgsPoint &p2, co
     if ( hasM )
       pointWkbType = QgsWkbTypes::addM( pointWkbType );
 
+    QgsDebugMsg( QString( "a1:%1 (%2), a3:%3 (%4), inc:%5, shi:?, cw:%6" )
+                 . arg( a1 ).  arg( a1 * 180 / M_PI )
+                 . arg( a3 ).  arg( a3 * 180 / M_PI )
+                 . arg( increment ).  arg( clockwise )
+               );
+
     //make sure the curve point p2 is part of the segmentized vertices. But only if p1 != p3
+    // TODO: make this a parameter
     bool addP2 = true;
     if ( qgsDoubleNear( circlePoint1.x(), circlePoint3.x() ) && qgsDoubleNear( circlePoint1.y(), circlePoint3.y() ) )
     {
       addP2 = false;
     }
-
-    for ( double angle = a1 + increment; angle < a3; angle += increment )
+    addP2 = false;
+
+    // As we're adding the last point in any case, we'll avoid
+    // including a point which is at less than 1% increment distance
+    // from it (may happen to find them due to numbers approximation).
+    // NOTE that this effectively allows in output some segments which
+    //      are more distant than requested. This is at most 1% off
+    //      from requested MaxAngle and less for MaxError.
+    double tolError = increment / 100;
+    double stopAngle = clockwise ? a3 - tolError : a3 - tolError;
+    QgsDebugMsg( QString( "stopAngle: %1 (%2)" ) . arg( stopAngle ) .arg( stopAngle * 180 / M_PI ) );
+    for ( double angle = a1 + increment; clockwise ? angle > stopAngle : angle < stopAngle; angle += increment )
     {
-      if ( ( addP2 && angle > a2 ) )
+      if ( addP2 && angle > a2 )
       {
-        stringPoints.insert( clockwise ? 0 : stringPoints.size(), circlePoint2 );
+        if ( clockwise )
+        {
+          if ( *stringPoints.begin() != circlePoint2 )
+          {
+            QgsDebugMsg( QString( "Adding control point, with angle %1 (%2)" ) . arg( a2 ) .arg( a2 * 180 / M_PI ) );
+            stringPoints.insert( 0, circlePoint2 );
+          }
+        }
+        else
+        {
+          if ( *stringPoints.rbegin() != circlePoint2 )
+          {
+            QgsDebugMsg( QString( "Adding control point, with angle %1 (%2)" ) . arg( a2 ) .arg( a2 * 180 / M_PI ) );
+            stringPoints.insert( stringPoints.size(), circlePoint2 );
+          }
+        }
         addP2 = false;
       }
 
+      QgsDebugMsg( QString( "SA - %1 (%2)" ) . arg( angle ) .arg( angle * 180 / M_PI ) );
+
       x = centerX + radius * std::cos( angle );
       y = centerY + radius * std::sin( angle );
 
-      if ( !hasZ && !hasM )
-      {
-        stringPoints.insert( clockwise ? 0 : stringPoints.size(), QgsPoint( x, y ) );
-        continue;
-      }
-
       if ( hasZ )
       {
         z = interpolateArcValue( angle, a1, a2, a3, circlePoint1.z(), circlePoint2.z(), circlePoint3.z() );
@@ -721,10 +810,18 @@ void QgsGeometryUtils::segmentizeArc( const QgsPoint &p1, const QgsPoint &p2, co
         m = interpolateArcValue( angle, a1, a2, a3, circlePoint1.m(), circlePoint2.m(), circlePoint3.m() );
       }
 
-      stringPoints.insert( clockwise ? 0 : stringPoints.size(), QgsPoint( pointWkbType, x, y, z, m ) );
+      QgsDebugMsg( QString( " -> POINT(%1 %2)" ) . arg( x ) .arg( y ) );
+      stringPoints.insert( stringPoints.size(), QgsPoint( pointWkbType, x, y, z, m ) );
     }
   }
-  stringPoints.insert( clockwise ? 0 : stringPoints.size(), circlePoint3 );
+  QgsDebugMsg( QString( " appending last point -> POINT(%1 %2)" ) .  arg( circlePoint3.x() ) .arg( circlePoint3.y() ) );
+  stringPoints.insert( stringPoints.size(), circlePoint3 );
+
+  // TODO: check if or implement QgsPointSequence directly taking an iterator to append
+  if ( reversed )
+  {
+    std::reverse( stringPoints.begin(), stringPoints.end() );
+  }
   points.append( stringPoints );
 }
 

tests/src/core/CMakeLists.txt

@@ -100,6 +100,7 @@ SET(TESTS
  testcontrastenhancements.cpp
  testqgscoordinatereferencesystem.cpp
  testqgscoordinatetransform.cpp
+ testqgscurve.cpp
  testqgsdatadefinedsizelegend.cpp
  testqgsdataitem.cpp
  testqgsdatasourceuri.cpp

tests/src/core/testqgscurve.cpp

@@ -0,0 +1,119 @@
+/***************************************************************************
+     testqgscurve.cpp
+     --------------------------------------
+    Date                 : 21 July 2017
+    Copyright            : (C) 2017 by Sandro Santilli
+    Email                : strk @ kbt.io
+ ***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+#include <QObject>
+#include <QString>
+#include <QApplication>
+#include <memory> // for unique_ptr
+
+//qgis includes...
+#include "qgsabstractgeometry.h"
+#include "qgscircularstring.h"
+#include "qgsgeometry.h"
+#include "qgsgeometryfactory.h"
+#include "qgslinestring.h"
+#include "qgspoint.h"
+#include "qgstest.h"
+#include "qgstestutils.h"
+
+/** \ingroup UnitTests
+ * This is a unit test for the operations on curve geometries
+ */
+class TestQgsCurve : public QObject
+{
+    Q_OBJECT
+
+  public:
+    TestQgsCurve() {};
+
+  private slots:
+    //void initTestCase();// will be called before the first testfunction is executed.
+    //void cleanupTestCase();// will be called after the last testfunction was executed.
+    //void init();// will be called before each testfunction is executed.
+    //void cleanup();// will be called after every testfunction.
+
+    void curveToLine();
+};
+
+
+#define TEST_C2L(circularString, tol, toltype, exp, prec) { \
+    std::unique_ptr< QgsLineString > lineString( \
+        circularString->curveToLine(tol, toltype) \
+                                               ); \
+    QVERIFY( lineString.get() ); \
+    QString wkt_out = lineString->asWkt(prec); \
+    QCOMPARE( wkt_out, QString( exp ) ); \
+    /* Test reverse */ \
+    std::unique_ptr< QgsCircularString > reversed( \
+        circularString->reversed() \
+                                                 ); \
+    lineString.reset( \
+                      reversed->curveToLine(tol, toltype) \
+                    ); \
+    wkt_out = lineString->asWkt(prec); \
+    lineString.reset( \
+                      reversed->curveToLine(tol, toltype) \
+                    ); \
+    std::unique_ptr< QgsLineString > expgeom( \
+        dynamic_cast<QgsLineString *>( \
+                                       QgsGeometryFactory::geomFromWkt( exp ) \
+                                     ) \
+                                            ); \
+    expgeom.reset( expgeom->reversed() ); \
+    QString exp_reversed = expgeom->asWkt(prec); \
+    QCOMPARE( wkt_out, exp_reversed ); \
+  }
+
+void TestQgsCurve::curveToLine()
+{
+  std::unique_ptr< QgsCircularString > circularString;
+
+  /* input: 2 quadrants arc (180 degrees, PI radians) */
+  circularString.reset( dynamic_cast<QgsCircularString *>(
+                          QgsGeometryFactory::geomFromWkt( QString(
+                                "CIRCULARSTRING(0 0,100 100,200 0)"
+                              ) )
+                        ) );
+  QVERIFY( circularString.get() );
+
+  /* op: Maximum of 10 units of difference, symmetric */
+  TEST_C2L( circularString, 10, QgsAbstractGeometry::MaximumDifference,
+            "LineString (0 0, 29.29 70.71, 100 100, 170.71 70.71, 200 0)", 2 );
+
+  /* op: Maximum of M_PI / 8 degrees of angle, (a)symmetric */
+  /* See https://issues.qgis.org/issues/16717 */
+  TEST_C2L( circularString, M_PI / 8, QgsAbstractGeometry::MaximumAngle,
+            "LineString (0 0, 7.61 38.27, 29.29 70.71, 61.73 92.39, 100 100, 138.27 92.39, 170.71 70.71, 192.39 38.27, 200 0)", 2 );
+
+  /* op: Maximum of 70 degrees of angle, symmetric */
+  /* See https://issues.qgis.org/issues/16722 */
+  TEST_C2L( circularString, 70 * M_PI / 180, QgsAbstractGeometry::MaximumAngle,
+            "LineString (0 0, 50 86.6, 150 86.6, 200 0)", 2 );
+
+  /* input: 2 arcs of 2 quadrants each (180 degrees + 180 degrees other direction) */
+  circularString.reset( dynamic_cast<QgsCircularString *>(
+                          QgsGeometryFactory::geomFromWkt( QString(
+                                "CIRCULARSTRING(0 0,100 100,200 0,300 -100,400 0)"
+                              ) )
+                        ) );
+  QVERIFY( circularString.get() );
+
+  /* op: Maximum of M_PI / 3 degrees of angle */
+  TEST_C2L( circularString, M_PI / 3, QgsAbstractGeometry::MaximumAngle,
+            "LineString (0 0, 50 86.6, 150 86.6, 200 0, 200 0, 250 -86.6, 350 -86.6, 400 0)", 2 );
+}
+
+
+QGSTEST_MAIN( TestQgsCurve )
+#include "testqgscurve.moc"

tests/src/core/testqgsgeometry.cpp

@@ -5319,6 +5319,7 @@ void TestQgsGeometry::directionNeutralSegmentation()
   QgsLineString *CCWLineString = CCWCircularString->curveToLine();
   QgsLineString *reversedCCWLineString = CCWLineString->reversed();
 
+  QCOMPARE( CWLineString->asWkt(), reversedCCWLineString->asWkt() );
   bool equal = ( *CWLineString == *reversedCCWLineString );
 
   delete CWCircularString;