Implement smoothing algorithm for geometries

Author: nyalldawson

python/core/qgsgeometry.sip

@@ -341,6 +341,17 @@ class QgsGeometry
 
     /** Returns a simplified version of this geometry using a specified tolerance value */
     QgsGeometry* simplify( double tolerance ) /Factory/;
+    
+    /**Smooths a geometry by rounding off corners using the Chaikin algorithm. This operation
+     * roughly doubles the number of vertices in a geometry.
+     * @param iterations number of smoothing iterations to run. More iterations results
+     * in a smoother geometry
+     * @param offset fraction of line to create new vertices along, between 0 and 1.0
+     * eg the default value of 0.25 will create new vertices 25% and 75% along each line segment
+     * of the geometry for each iteration. Smaller values result in "tighter" smoothing.
+     * @note added in 2.9
+     */
+    QgsGeometry* smooth( const unsigned int iterations = 1, const double offset = 0.25 ) /Factory/;
 
     /** Returns the center of mass of a geometry
     * @note for line based geometries, the center point of the line is returned,
@@ -353,7 +364,7 @@ class QgsGeometry
     /** Returns the smallest convex polygon that contains all the points in the geometry. */
     QgsGeometry* convexHull() /Factory/;
 
-    /* Return interpolated point on line at distance */
+    /** Return interpolated point on line at distance */
     QgsGeometry* interpolate( double distance ) /Factory/;
 
     /** Returns a geometry representing the points shared by this geometry and other. */
@@ -495,7 +506,17 @@ class QgsGeometry
      * @note added in QGIS 2.9
     */
     static bool compare( const QgsPolygon& p1, const QgsPolygon& p2, double epsilon = 4 * DBL_EPSILON );
-
+    
+    /** Compares two multipolygons for equality within a specified tolerance.
+     * @param p1 first multipolygon
+     * @param p2 second multipolygon
+     * @param epsilon maximum difference for coordinates between the multipolygons
+     * @returns true if multipolygons have the same number of polygons, the polygons have the same number
+     * of rings, and each ring has the same number of points and all points are equal within the specified
+     * tolerance
+     * @note added in QGIS 2.9
+    */
+    static bool compare( const QgsMultiPolygon& p1, const QgsMultiPolygon& p2, double epsilon = 4 * DBL_EPSILON );
 
 }; // class QgsGeometry
 

src/core/qgsgeometry.cpp

@@ -5797,6 +5797,116 @@ QgsGeometry* QgsGeometry::simplify( double tolerance )
   CATCH_GEOS( 0 )
 }
 
+QgsGeometry* QgsGeometry::smooth( const unsigned int iterations, const double offset )
+{
+  switch ( wkbType() )
+  {
+    case QGis::WKBPoint:
+    case QGis::WKBPoint25D:
+    case QGis::WKBMultiPoint:
+    case QGis::WKBMultiPoint25D:
+      //can't smooth a point based geometry
+      return new QgsGeometry( *this );
+
+    case QGis::WKBLineString:
+    case QGis::WKBLineString25D:
+    {
+      QgsPolyline line = asPolyline();
+      return QgsGeometry::fromPolyline( smoothLine( line, iterations, offset ) );
+    }
+
+    case QGis::WKBMultiLineString:
+    case QGis::WKBMultiLineString25D:
+    {
+      QgsMultiPolyline multiline = asMultiPolyline();
+      QgsMultiPolyline resultMultiline;
+      QgsMultiPolyline::const_iterator lineIt = multiline.constBegin();
+      for ( ; lineIt != multiline.constEnd(); ++lineIt )
+      {
+        resultMultiline << smoothLine( *lineIt, iterations, offset );
+      }
+      return QgsGeometry::fromMultiPolyline( resultMultiline );
+    }
+
+    case QGis::WKBPolygon:
+    case QGis::WKBPolygon25D:
+    {
+      QgsPolygon poly = asPolygon();
+      return QgsGeometry::fromPolygon( smoothPolygon( poly, iterations, offset ) );
+    }
+
+    case QGis::WKBMultiPolygon:
+    case QGis::WKBMultiPolygon25D:
+    {
+      QgsMultiPolygon multipoly = asMultiPolygon();
+      QgsMultiPolygon resultMultipoly;
+      QgsMultiPolygon::const_iterator polyIt = multipoly.constBegin();
+      for ( ; polyIt != multipoly.constEnd(); ++polyIt )
+      {
+        resultMultipoly << smoothPolygon( *polyIt, iterations, offset );
+      }
+      return QgsGeometry::fromMultiPolygon( resultMultipoly );
+    }
+    break;
+
+    case QGis::WKBUnknown:
+    default:
+      return new QgsGeometry( *this );
+  }
+}
+
+inline QgsPoint interpolatePointOnLine( const QgsPoint& p1, const QgsPoint& p2, const double offset )
+{
+  double deltaX = p2.x() - p1.x();
+  double deltaY = p2.y() - p1.y();
+  return QgsPoint( p1.x() + deltaX * offset, p1.y() + deltaY * offset );
+}
+
+QgsPolyline QgsGeometry::smoothLine( const QgsPolyline& polyline, const unsigned int iterations, const double offset ) const
+{
+  QgsPolyline result = polyline;
+  for ( unsigned int iteration = 0; iteration < iterations; ++iteration )
+  {
+    QgsPolyline outputLine = QgsPolyline();
+    for ( int i = 0; i < result.count() - 1; i++ )
+    {
+      const QgsPoint& p1 = result.at( i );
+      const QgsPoint& p2 = result.at( i + 1 );
+      outputLine << ( i == 0 ? result.at( i ) : interpolatePointOnLine( p1, p2, offset ) );
+      outputLine << ( i == result.count() - 2 ? result.at( i + 1 ) : interpolatePointOnLine( p1, p2, 1.0 - offset ) );
+    }
+    result = outputLine;
+  }
+  return result;
+}
+
+QgsPolygon QgsGeometry::smoothPolygon( const QgsPolygon& polygon, const unsigned int iterations, const double offset ) const
+{
+  QgsPolygon resultPoly;
+  QgsPolygon::const_iterator ringIt = polygon.constBegin();
+  for ( ; ringIt != polygon.constEnd(); ++ringIt )
+  {
+    QgsPolyline resultRing = *ringIt;
+    for ( unsigned int iteration = 0; iteration < iterations; ++iteration )
+    {
+      QgsPolyline outputRing = QgsPolyline();
+      for ( int i = 0; i < resultRing.count() - 1; ++i )
+      {
+        const QgsPoint& p1 = resultRing.at( i );
+        const QgsPoint& p2 = resultRing.at( i + 1 );
+        outputRing << interpolatePointOnLine( p1, p2, offset );
+        outputRing << interpolatePointOnLine( p1, p2, 1.0 - offset );
+      }
+      //close polygon
+      outputRing << outputRing.at( 0 );
+
+      resultRing = outputRing;
+    }
+    resultPoly << resultRing;
+  }
+  return resultPoly;
+}
+
 QgsGeometry* QgsGeometry::centroid()
 {
   if ( mDirtyGeos )
@@ -6658,3 +6768,17 @@ bool QgsGeometry::compare( const QgsPolygon &p1, const QgsPolygon &p2, double ep
   }
   return true;
 }
+
+
+bool QgsGeometry::compare( const QgsMultiPolygon &p1, const QgsMultiPolygon &p2, double epsilon )
+{
+  if ( p1.count() != p2.count() )
+    return false;
+
+  for ( int i = 0; i < p1.count(); ++i )
+  {
+    if ( !QgsGeometry::compare( p1.at( i ), p2.at( i ), epsilon ) )
+      return false;
+  }
+  return true;
+}

src/core/qgsgeometry.h

@@ -384,6 +384,17 @@ class CORE_EXPORT QgsGeometry
     /** Returns a simplified version of this geometry using a specified tolerance value */
     QgsGeometry* simplify( double tolerance );
 
+    /**Smooths a geometry by rounding off corners using the Chaikin algorithm. This operation
+     * roughly doubles the number of vertices in a geometry.
+     * @param iterations number of smoothing iterations to run. More iterations results
+     * in a smoother geometry
+     * @param offset fraction of line to create new vertices along, between 0 and 1.0
+     * eg the default value of 0.25 will create new vertices 25% and 75% along each line segment
+     * of the geometry for each iteration. Smaller values result in "tighter" smoothing.
+     * @note added in 2.9
+     */
+    QgsGeometry* smooth( const unsigned int iterations = 1, const double offset = 0.25 );
+
     /** Returns the center of mass of a geometry
     * @note for line based geometries, the center point of the line is returned,
     * and for point based geometries, the point itself is returned */
@@ -395,7 +406,7 @@ class CORE_EXPORT QgsGeometry
     /** Returns the smallest convex polygon that contains all the points in the geometry. */
     QgsGeometry* convexHull();
 
-    /* Return interpolated point on line at distance */
+    /** Return interpolated point on line at distance */
     QgsGeometry* interpolate( double distance );
 
     /** Returns a geometry representing the points shared by this geometry and other. */
@@ -541,6 +552,17 @@ class CORE_EXPORT QgsGeometry
     */
     static bool compare( const QgsPolygon& p1, const QgsPolygon& p2, double epsilon = 4 * DBL_EPSILON );
 
+    /** Compares two multipolygons for equality within a specified tolerance.
+     * @param p1 first multipolygon
+     * @param p2 second multipolygon
+     * @param epsilon maximum difference for coordinates between the multipolygons
+     * @returns true if multipolygons have the same number of polygons, the polygons have the same number
+     * of rings, and each ring has the same number of points and all points are equal within the specified
+     * tolerance
+     * @note added in QGIS 2.9
+    */
+    static bool compare( const QgsMultiPolygon& p1, const QgsMultiPolygon& p2, double epsilon = 4 * DBL_EPSILON );
+
   private:
     // Private variables
 
@@ -640,6 +662,11 @@ class CORE_EXPORT QgsGeometry
     @return the reshaped polygon or 0 in case of error*/
     static GEOSGeometry* reshapePolygon( const GEOSGeometry* polygon, const GEOSGeometry* reshapeLineGeos );
 
+    /**Smooths a polygon using the Chaikin algorithm*/
+    QgsPolygon smoothPolygon( const QgsPolygon &polygon, const unsigned int iterations = 1, const double offset = 0.25 ) const;
+    /**Smooths a polyline using the Chaikin algorithm*/
+    QgsPolyline smoothLine( const QgsPolyline &polyline, const unsigned int iterations = 1, const double offset = 0.25 ) const;
+
     /**Nodes together a split line and a (multi-) polygon geometry in a multilinestring
      @return the noded multiline geometry or 0 in case of error. The calling function takes ownership of the node geometry*/
     static GEOSGeometry* nodeGeometries( const GEOSGeometry *splitLine, const GEOSGeometry *poly );

tests/src/core/testqgsgeometry.cpp

@@ -74,6 +74,7 @@ class TestQgsGeometry : public QObject
     void differenceCheck1();
     void differenceCheck2();
     void bufferCheck();
+    void smoothCheck();
 
   private:
     /** A helper method to do a render check to see if the geometry op is as expected */
@@ -550,6 +551,71 @@ void TestQgsGeometry::bufferCheck()
   dumpPolygon( myPolygon );
   QVERIFY( renderCheck( "geometry_bufferCheck", "Checking buffer(10,10) of B", 10 ) );
 }
+
+void TestQgsGeometry::smoothCheck()
+{
+  //can't smooth a point
+  QString wkt = "POINT(40 50)";
+  QScopedPointer<QgsGeometry> geom( QgsGeometry::fromWkt( wkt ) );
+  QgsGeometry* result = geom->smooth( 1, 0.25 );
+  QString obtained = result->exportToWkt();
+  delete result;
+  QCOMPARE( obtained, wkt );
+
+  //linestring
+  wkt = "LINESTRING(0 0, 10 0, 10 10, 20 10)";
+  geom.reset( QgsGeometry::fromWkt( wkt ) );
+  result = geom->smooth( 1, 0.25 );
+  QgsPolyline line = result->asPolyline();
+  delete result;
+  QgsPolyline expectedLine;
+  expectedLine << QgsPoint( 0, 0 ) << QgsPoint( 7.5, 0 ) << QgsPoint( 10.0, 2.5 )
+  << QgsPoint( 10.0, 7.5 ) << QgsPoint( 12.5, 10.0 ) << QgsPoint( 20.0, 10.0 );
+  QVERIFY( QgsGeometry::compare( line, expectedLine ) );
+
+  wkt = "MULTILINESTRING((0 0, 10 0, 10 10, 20 10),(30 30, 40 30, 40 40, 50 40))";
+  geom.reset( QgsGeometry::fromWkt( wkt ) );
+  result = geom->smooth( 1, 0.25 );
+  QgsMultiPolyline multiLine = result->asMultiPolyline();
+  delete result;
+  QgsMultiPolyline expectedMultiline;
+  expectedMultiline << ( QgsPolyline() << QgsPoint( 0, 0 ) << QgsPoint( 7.5, 0 ) << QgsPoint( 10.0, 2.5 )
+                         <<  QgsPoint( 10.0, 7.5 ) << QgsPoint( 12.5, 10.0 ) << QgsPoint( 20.0, 10.0 ) )
+  << ( QgsPolyline() << QgsPoint( 30.0, 30.0 ) << QgsPoint( 37.5, 30.0 ) << QgsPoint( 40.0, 32.5 )
+       << QgsPoint( 40.0, 37.5 ) << QgsPoint( 42.5, 40.0 ) << QgsPoint( 50.0, 40.0 ) );
+  QVERIFY( QgsGeometry::compare( multiLine, expectedMultiline ) );
+
+  //polygon
+  wkt = "POLYGON((0 0, 10 0, 10 10, 0 10, 0 0 ),(2 2, 4 2, 4 4, 2 4, 2 2))";
+  geom.reset( QgsGeometry::fromWkt( wkt ) );
+  result = geom->smooth( 1, 0.25 );
+  QgsPolygon poly = result->asPolygon();
+  delete result;
+  QgsPolygon expectedPolygon;
+  expectedPolygon << ( QgsPolyline() << QgsPoint( 2.5, 0 ) << QgsPoint( 7.5, 0 ) << QgsPoint( 10.0, 2.5 )
+                       <<  QgsPoint( 10.0, 7.5 ) << QgsPoint( 7.5, 10.0 ) << QgsPoint( 2.5, 10.0 ) << QgsPoint( 0, 7.5 )
+                       << QgsPoint( 0, 2.5 ) << QgsPoint( 2.5, 0 ) )
+  << ( QgsPolyline() << QgsPoint( 2.5, 2.0 ) << QgsPoint( 3.5, 2.0 ) << QgsPoint( 4.0, 2.5 )
+       << QgsPoint( 4.0, 3.5 ) << QgsPoint( 3.5, 4.0 ) << QgsPoint( 2.5, 4.0 )
+       << QgsPoint( 2.0, 3.5 ) << QgsPoint( 2.0, 2.5 ) << QgsPoint( 2.5, 2.0 ) );
+  QVERIFY( QgsGeometry::compare( poly, expectedPolygon ) );
+
+  //multipolygon
+  wkt = "MULTIPOLYGON(((0 0, 10 0, 10 10, 0 10, 0 0 )),((2 2, 4 2, 4 4, 2 4, 2 2)))";
+  geom.reset( QgsGeometry::fromWkt( wkt ) );
+  result = geom->smooth( 1, 0.1 );
+  QgsMultiPolygon multipoly = result->asMultiPolygon();
+  delete result;
+  QgsMultiPolygon expectedMultiPoly;
+  expectedMultiPoly << ( QgsPolygon() << ( QgsPolyline() << QgsPoint( 1.0, 0 ) << QgsPoint( 9, 0 ) << QgsPoint( 10.0, 1 )
+                         <<  QgsPoint( 10.0, 9 ) << QgsPoint( 9, 10.0 ) << QgsPoint( 1, 10.0 ) << QgsPoint( 0, 9 )
+                         << QgsPoint( 0, 1 ) << QgsPoint( 1, 0 ) ) )
+  << ( QgsPolygon() << ( QgsPolyline() << QgsPoint( 2.2, 2.0 ) << QgsPoint( 3.8, 2.0 ) << QgsPoint( 4.0, 2.2 )
+                         <<  QgsPoint( 4.0, 3.8 ) << QgsPoint( 3.8, 4.0 ) << QgsPoint( 2.2, 4.0 ) << QgsPoint( 2.0, 3.8 )
+                         << QgsPoint( 2, 2.2 ) << QgsPoint( 2.2, 2 ) ) );
+  QVERIFY( QgsGeometry::compare( multipoly, expectedMultiPoly ) );
+}
+
 bool TestQgsGeometry::renderCheck( QString theTestName, QString theComment, int mismatchCount )
 {
   mReport += "<h2>" + theTestName + "</h2>\n";