Fix crash using Visvalingam simplification, add test

Author: nyalldawson

src/core/CMakeLists.txt

@@ -44,6 +44,8 @@ SET(QGIS_CORE_SRCS
   symbology-ng/qgssymbol.cpp
   symbology-ng/qgsvectorfieldsymbollayer.cpp
 
+  simplify/effectivearea.cpp
+
   diagram/qgsdiagram.cpp
   diagram/qgshistogramdiagram.cpp
   diagram/qgspiediagram.cpp

src/core/qgsmaptopixelgeometrysimplifier.cpp

@@ -70,17 +70,7 @@ bool QgsMapToPixelSimplifier::equalSnapToGrid( double x1, double y1, double x2,
 // https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/effectivearea.h
 // https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/effectivearea.c
 
-#define LWDEBUG //
-#define LWDEBUGF //
-#define FP_MAX qMax
-#define FLAGS_GET_Z( flags ) ( ( flags ) & 0x01 )
-#define LW_MSG_MAXLEN 256
-#define lwalloc qgsMalloc
-#define lwfree qgsFree
-#define lwerror qWarning
-
 #include "simplify/effectivearea.h"
-#include "simplify/effectivearea.c"
 
 //////////////////////////////////////////////////////////////////////////////////////////////
 
@@ -175,75 +165,80 @@ QgsGeometry QgsMapToPixelSimplifier::simplifyGeometry(
     }
 
     // Process each vertex...
-    if ( simplifyAlgorithm == SnapToGrid )
+    switch ( simplifyAlgorithm )
     {
-      double gridOriginX = envelope.xMinimum();
-      double gridOriginY = envelope.yMinimum();
-
-      // Use a factor for the maximum displacement distance for simplification, similar as GeoServer does
-      float gridInverseSizeXY = map2pixelTol != 0 ? ( float )( 1.0f / ( 0.8 * map2pixelTol ) ) : 0.0f;
-
-      for ( int i = 0; i < numPoints; ++i )
+      case SnapToGrid:
       {
-        x = srcCurve.xAt( i );
-        y = srcCurve.yAt( i );
+        double gridOriginX = envelope.xMinimum();
+        double gridOriginY = envelope.yMinimum();
 
-        if ( i == 0 ||
-             !isGeneralizable ||
-             !equalSnapToGrid( x, y, lastX, lastY, gridOriginX, gridOriginY, gridInverseSizeXY ) ||
-             ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
+        // Use a factor for the maximum displacement distance for simplification, similar as GeoServer does
+        float gridInverseSizeXY = map2pixelTol != 0 ? ( float )( 1.0f / ( 0.8 * map2pixelTol ) ) : 0.0f;
+
+        for ( int i = 0; i < numPoints; ++i )
         {
-          output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPointV2( x, y ) );
-          lastX = x;
-          lastY = y;
+          x = srcCurve.xAt( i );
+          y = srcCurve.yAt( i );
+
+          if ( i == 0 ||
+               !isGeneralizable ||
+               !equalSnapToGrid( x, y, lastX, lastY, gridOriginX, gridOriginY, gridInverseSizeXY ) ||
+               ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
+          {
+            output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPointV2( x, y ) );
+            lastX = x;
+            lastY = y;
+          }
+
+          r.combineExtentWith( x, y );
         }
-
-        r.combineExtentWith( x, y );
+        break;
       }
-    }
-    else if ( simplifyAlgorithm == Visvalingam )
-    {
-      map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'Area' calculations.
 
-      EFFECTIVE_AREAS* ea;
-      ea = initiate_effectivearea( srcCurve );
+      case Visvalingam:
+      {
+        map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'Area' calculations.
 
-      int set_area = 0;
-      ptarray_calc_areas( ea, isaLinearRing ? 4 : 2, set_area, map2pixelTol );
+        EFFECTIVE_AREAS ea( srcCurve );
 
-      for ( int i = 0; i < numPoints; ++i )
-      {
-        if ( ea->res_arealist[ i ] > map2pixelTol )
+        int set_area = 0;
+        ptarray_calc_areas( &ea, isaLinearRing ? 4 : 2, set_area, map2pixelTol );
+
+        for ( int i = 0; i < numPoints; ++i )
         {
-          output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), ea->inpts.at( i ) );
+          if ( ea.res_arealist[ i ] > map2pixelTol )
+          {
+            output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), ea.inpts.at( i ) );
+          }
         }
+        break;
       }
-      destroy_effectivearea( ea );
-    }
-    else
-    {
-      map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'LengthSquare' calculations.
 
-      for ( int i = 0; i < numPoints; ++i )
+      case Distance:
       {
-        x = srcCurve.xAt( i );
-        y = srcCurve.yAt( i );
-
-        isLongSegment = false;
+        map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'LengthSquare' calculations.
 
-        if ( i == 0 ||
-             !isGeneralizable ||
-             ( isLongSegment = ( calculateLengthSquared2D( x, y, lastX, lastY ) > map2pixelTol ) ) ||
-             ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
+        for ( int i = 0; i < numPoints; ++i )
         {
-          output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPointV2( x, y ) );
-          lastX = x;
-          lastY = y;
+          x = srcCurve.xAt( i );
+          y = srcCurve.yAt( i );
 
-          hasLongSegments |= isLongSegment;
-        }
+          isLongSegment = false;
+
+          if ( i == 0 ||
+               !isGeneralizable ||
+               ( isLongSegment = ( calculateLengthSquared2D( x, y, lastX, lastY ) > map2pixelTol ) ) ||
+               ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
+          {
+            output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPointV2( x, y ) );
+            lastX = x;
+            lastY = y;
 
-        r.combineExtentWith( x, y );
+            hasLongSegments |= isLongSegment;
+          }
+
+          r.combineExtentWith( x, y );
+        }
       }
     }
 

src/core/simplify/effectivearea.c src/core/simplify/effectivearea.cpp

@@ -24,25 +24,6 @@
 
 #include "effectivearea.h"
 
-EFFECTIVE_AREAS* initiate_effectivearea( const QgsCurve& inpts )
-{
-  //LWDEBUG( 2, "Entered  initiate_effectivearea" );
-  EFFECTIVE_AREAS *ea;
-  ea = ( EFFECTIVE_AREAS* )lwalloc( sizeof( EFFECTIVE_AREAS ) );
-  inpts.points( ea->inpts );
-  ea->is3d = inpts.is3D();
-  ea->initial_arealist = ( areanode* )lwalloc( ea->inpts.size() * sizeof( areanode ) );
-  ea->res_arealist = ( double* )lwalloc( ea->inpts.size() * sizeof( double ) );
-  return ea;
-}
-
-void destroy_effectivearea( EFFECTIVE_AREAS *ea )
-{
-  lwfree( ea->initial_arealist );
-  lwfree( ea->res_arealist );
-  lwfree( ea );
-}
-
 static MINHEAP initiate_minheap( int npoints )
 {
   MINHEAP tree;

src/core/simplify/effectivearea.h

@@ -22,48 +22,75 @@
  *
  **********************************************************************/
 
+#include "qgsabstractgeometry.h"
+#include "qgscurve.h"
+#include "qgspointv2.h"
+
 #ifndef _EFFECTIVEAREA_H
 #define _EFFECTIVEAREA_H 1
 
+
+#define LWDEBUG //
+#define LWDEBUGF //
+#define FP_MAX qMax
+#define FLAGS_GET_Z( flags ) ( ( flags ) & 0x01 )
+#define LW_MSG_MAXLEN 256
+#define lwalloc qgsMalloc
+#define lwfree qgsFree
+#define lwerror qWarning
+
+
 /**
  * This structure is placed in an array with one member per point.
  * It has links into the minheap rtee and kepps track of eliminated points.
  */
-typedef struct
+struct areanode
 {
   double area;
   int treeindex;
   int prev;
   int next;
-} areanode;
+};
 
 /**
  * This structure holds a minheap tree that is used to keep track of what points
  * that has the smallest effective area.
  * When elliminating points the neighbor points has its effective area affected
  * and the minheap helps to resort efficient.
  */
-typedef struct
+struct MINHEAP
 {
   int maxSize;
   int usedSize;
   areanode **key_array;
-} MINHEAP;
+};
 
 /**
  * Structure to hold pointarray and it's arealist.
  */
-typedef struct
+struct EFFECTIVE_AREAS
 {
+  EFFECTIVE_AREAS( const QgsCurve& curve )
+      : is3d( curve.is3D() )
+      , initial_arealist( nullptr )
+      , res_arealist( nullptr )
+  {
+    curve.points( inpts );
+    initial_arealist = new areanode[ inpts.size()];
+    res_arealist = new double[ inpts.size()];
+  }
+
+  ~EFFECTIVE_AREAS()
+  {
+    delete [] initial_arealist;
+    delete [] res_arealist;
+  }
+
   bool is3d;
   QgsPointSequence inpts;
   areanode *initial_arealist;
   double *res_arealist;
-} EFFECTIVE_AREAS;
-
-EFFECTIVE_AREAS* initiate_effectivearea( const QgsCurve &inpts );
-
-void destroy_effectivearea( EFFECTIVE_AREAS *ea );
+};
 
 void ptarray_calc_areas( EFFECTIVE_AREAS *ea, int avoid_collaps, int set_area, double trshld );
 

tests/src/core/testqgsmaptopixelgeometrysimplifier.cpp

@@ -74,6 +74,7 @@ class TestQgsMapToPixelGeometrySimplifier : public QObject
     void testIsGeneralizableByMapBoundingBox();
     void testWkbDimensionMismatch();
     void testCircularString();
+    void testVisvalingam();
 
 };
 
@@ -194,5 +195,16 @@ void TestQgsMapToPixelGeometrySimplifier::testCircularString()
   QCOMPARE( simplifier.simplify( g ).exportToWkt(), WKT );
 }
 
+void TestQgsMapToPixelGeometrySimplifier::testVisvalingam()
+{
+  QString wkt( "LineString (0 0, 30 0, 31 30, 32 0, 40 0, 41 100, 42 0, 50 0)" );
+  QgsGeometry g = QgsGeometry::fromWkt( wkt );
+
+  const QgsMapToPixelSimplifier simplifier( QgsMapToPixelSimplifier::SimplifyGeometry, 7, QgsMapToPixelSimplifier::Visvalingam );
+  QString expectedWkt( "LineString (0 0, 40 0, 41 100, 42 0, 50 0)" );
+
+  QCOMPARE( simplifier.simplify( g ).exportToWkt(), expectedWkt );
+}
+
 QTEST_MAIN( TestQgsMapToPixelGeometrySimplifier )
 #include "testqgsmaptopixelgeometrysimplifier.moc"