Also port processing densify to distance to c++

- Add QgsGeometry method to densify by distance
- Fix bug in processing algorithm which resulted in duplicate
vertices and incorrectly spaced extra vertices

Author: nyalldawson

python/core/geometry/qgsgeometry.sip

@@ -543,6 +543,7 @@ class QgsGeometry
     QgsGeometry extendLine( double startDistance, double endDistance ) const;
     QgsGeometry simplify( double tolerance ) const;
     QgsGeometry densifyByCount( int extraNodesPerSegment ) const;
+    QgsGeometry densifyByDistance( double distance ) const;
     QgsGeometry centroid() const;
 
     /**

python/plugins/processing/algs/help/qgis.yaml

@@ -139,9 +139,13 @@ qgis:densifygeometries:
   The number of new vertices to add to each feature geometry is specified as an input parameter.
 
 qgis:densifygeometriesgivenaninterval:
-  This algorithm takes a polygon or line layer and generates a new one in which the geometries have a larger number of vertices than the original one.
+  This algorithm takes a polygon or line layer and generates a new one in which the geometries have a larger number of vertices than the original one. The geometries are densified by adding regularly placed extra nodes inside each segment so that the maximum distance between any two nodes does not exceed the specified distance.
+
+  E.g. specifying a distance 3 would cause the segment [0 0] -> [10 0] to be converted to [0 0] -> [2.5 0] -> [5 0] -> [7.5 0] -> [10 0], since 3 extra nodes are required on the segment and spacing these at 2.5 increments allows them to be evenly spaced over the segment.
+
+  If the geometries have z or m values present then these will be linearly interpolated at the added nodes.
 
-  The number of new vertices depends on the length of the geometry, and is specified as a distance between them. The distance is expressed in the same units used by the layer CRS.
+  The distance is expressed in the same units used by the layer CRS.
 
 qgis:difference: >
   This algorithm extracts features from the Input layer that fall outside, or partially overlap, features in the Difference layer. Input layer features that partially overlap the difference layer feature(s) are split along the boundary of the difference layer feature(s) and only the portions outside the difference layer features are retained.

python/plugins/processing/algs/qgis/DensifyGeometriesInterval.py

@@ -71,61 +71,10 @@ def processAlgorithm(self, feedback):
         for current, f in enumerate(features):
             feature = f
             if feature.hasGeometry():
-                new_geometry = self.densifyGeometry(feature.geometry(), interval,
-                                                    isPolygon)
+                new_geometry = feature.geometry().densifyByCount(float(interval))
                 feature.setGeometry(new_geometry)
             writer.addFeature(feature)
 
             feedback.setProgress(int(current * total))
 
         del writer
-
-    def densifyGeometry(self, geometry, interval, isPolygon):
-        output = []
-        if isPolygon:
-            if geometry.isMultipart():
-                polygons = geometry.asMultiPolygon()
-                for poly in polygons:
-                    p = []
-                    for ring in poly:
-                        p.append(self.densify(ring, interval))
-                    output.append(p)
-                return QgsGeometry.fromMultiPolygon(output)
-            else:
-                rings = geometry.asPolygon()
-                for ring in rings:
-                    output.append(self.densify(ring, interval))
-                return QgsGeometry.fromPolygon(output)
-        else:
-            if geometry.isMultipart():
-                lines = geometry.asMultiPolyline()
-                for points in lines:
-                    output.append(self.densify(points, interval))
-                return QgsGeometry.fromMultiPolyline(output)
-            else:
-                points = geometry.asPolyline()
-                output = self.densify(points, interval)
-                return QgsGeometry.fromPolyline(output)
-
-    def densify(self, polyline, interval):
-        output = []
-        for i in range(len(polyline) - 1):
-            p1 = polyline[i]
-            p2 = polyline[i + 1]
-            output.append(p1)
-
-            # Calculate necessary number of points between p1 and p2
-            pointsNumber = sqrt(p1.sqrDist(p2)) / interval
-            if pointsNumber > 1:
-                multiplier = 1.0 / float(pointsNumber)
-            else:
-                multiplier = 1
-            for j in range(int(pointsNumber)):
-                delta = multiplier * (j + 1)
-                x = p1.x() + delta * (p2.x() - p1.x())
-                y = p1.y() + delta * (p2.y() - p1.y())
-                output.append(QgsPoint(x, y))
-                if j + 1 == pointsNumber:
-                    break
-        output.append(polyline[len(polyline) - 1])
-        return output

src/core/geometry/qgsgeometry.cpp

@@ -1583,6 +1583,13 @@ QgsGeometry QgsGeometry::densifyByCount( int extraNodesPerSegment ) const
   return engine.densifyByCount( extraNodesPerSegment );
 }
 
+QgsGeometry QgsGeometry::densifyByDistance( double distance ) const
+{
+  QgsInternalGeometryEngine engine( *this );
+
+  return engine.densifyByDistance( distance );
+}
+
 QgsGeometry QgsGeometry::centroid() const
 {
   if ( !d->geometry )

src/core/geometry/qgsgeometry.h

@@ -618,9 +618,26 @@ class CORE_EXPORT QgsGeometry
      * at the added nodes.
      * Curved geometry types are automatically segmentized by this routine.
      * @node added in QGIS 3.0
+     * @see densifyByDistance()
      */
     QgsGeometry densifyByCount( int extraNodesPerSegment ) const;
 
+    /**
+     * Densifies the geometry by adding regularly placed extra nodes inside each segment
+     * so that the maximum distance between any two nodes does not exceed the
+     * specified \a distance.
+     * E.g. specifying a distance 3 would cause the segment [0 0] -> [10 0]
+     * to be converted to [0 0] -> [2.5 0] -> [5 0] -> [7.5 0] -> [10 0], since
+     * 3 extra nodes are required on the segment and spacing these at 2.5 increments
+     * allows them to be evenly spaced over the segment.
+     * If the geometry has z or m values present then these will be linearly interpolated
+     * at the added nodes.
+     * Curved geometry types are automatically segmentized by this routine.
+     * @node added in QGIS 3.0
+     * @see densifyByCount()
+     */
+    QgsGeometry densifyByDistance( double distance ) const;
+
     /**
      * Returns the center of mass of a geometry.
      * @note for line based geometries, the center point of the line is returned,

src/core/geometry/qgsinternalgeometryengine.cpp

@@ -477,7 +477,8 @@ QgsGeometry QgsInternalGeometryEngine::orthogonalize( double tolerance, int maxI
   }
 }
 
-QgsLineString *doDensifyByCount( QgsLineString *ring, int extraNodesPerSegment )
+// if extraNodesPerSegment < 0, then use distance based mode
+QgsLineString *doDensify( QgsLineString *ring, int extraNodesPerSegment = -1, double distance = 1 )
 {
   QgsPointSequence out;
   double multiplier = 1.0 / double( extraNodesPerSegment + 1 );
@@ -503,6 +504,7 @@ QgsLineString *doDensifyByCount( QgsLineString *ring, int extraNodesPerSegment )
   double yOut = 0;
   double zOut = 0;
   double mOut = 0;
+  int extraNodesThisSegment = extraNodesPerSegment;
   for ( int i = 0; i < nPoints - 1; ++i )
   {
     x1 = ring->xAt( i );
@@ -521,7 +523,16 @@ QgsLineString *doDensifyByCount( QgsLineString *ring, int extraNodesPerSegment )
     }
 
     out << QgsPointV2( outType, x1, y1, z1, m1 );
-    for ( int j = 0; j < extraNodesPerSegment; ++j )
+
+    if ( extraNodesPerSegment < 0 )
+    {
+      // distance mode
+      extraNodesThisSegment = floor( sqrt( ( x2 - x1 ) * ( x2 - x1 ) + ( y2 - y1 ) * ( y2 - y1 ) ) / distance );
+      if ( extraNodesThisSegment >= 1 )
+        multiplier = 1.0 / ( extraNodesThisSegment + 1 );
+    }
+
+    for ( int j = 0; j < extraNodesThisSegment; ++j )
     {
       double delta = multiplier * ( j + 1 );
       xOut = x1 + delta * ( x2 - x1 );
@@ -542,7 +553,7 @@ QgsLineString *doDensifyByCount( QgsLineString *ring, int extraNodesPerSegment )
   return result;
 }
 
-QgsAbstractGeometry *densifyGeometryByCount( const QgsAbstractGeometry *geom, int extraNodesPerSegment )
+QgsAbstractGeometry *densifyGeometry( const QgsAbstractGeometry *geom, int extraNodesPerSegment = 1, double distance = 1 )
 {
   std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
   if ( QgsWkbTypes::isCurvedType( geom->wkbType() ) )
@@ -553,20 +564,20 @@ QgsAbstractGeometry *densifyGeometryByCount( const QgsAbstractGeometry *geom, in
 
   if ( QgsWkbTypes::geometryType( geom->wkbType() ) == QgsWkbTypes::LineGeometry )
   {
-    return doDensifyByCount( static_cast< QgsLineString * >( geom->clone() ), extraNodesPerSegment );
+    return doDensify( static_cast< QgsLineString * >( geom->clone() ), extraNodesPerSegment, distance );
   }
   else
   {
     // polygon
     const QgsPolygonV2 *polygon = static_cast< const QgsPolygonV2 * >( geom );
     QgsPolygonV2 *result = new QgsPolygonV2();
 
-    result->setExteriorRing( doDensifyByCount( static_cast< QgsLineString * >( polygon->exteriorRing()->clone() ),
-                             extraNodesPerSegment ) );
+    result->setExteriorRing( doDensify( static_cast< QgsLineString * >( polygon->exteriorRing()->clone() ),
+                                        extraNodesPerSegment, distance ) );
     for ( int i = 0; i < polygon->numInteriorRings(); ++i )
     {
-      result->addInteriorRing( doDensifyByCount( static_cast< QgsLineString * >( polygon->interiorRing( i )->clone() ),
-                               extraNodesPerSegment ) );
+      result->addInteriorRing( doDensify( static_cast< QgsLineString * >( polygon->interiorRing( i )->clone() ),
+                                          extraNodesPerSegment, distance ) );
     }
 
     return result;
@@ -592,7 +603,42 @@ QgsGeometry QgsInternalGeometryEngine::densifyByCount( int extraNodesPerSegment
     geometryList.reserve( numGeom );
     for ( int i = 0; i < numGeom; ++i )
     {
-      geometryList << densifyGeometryByCount( gc->geometryN( i ), extraNodesPerSegment );
+      geometryList << densifyGeometry( gc->geometryN( i ), extraNodesPerSegment );
+    }
+
+    QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
+    Q_FOREACH ( QgsAbstractGeometry *g, geometryList )
+    {
+      first.addPart( g );
+    }
+    return first;
+  }
+  else
+  {
+    return QgsGeometry( densifyGeometry( mGeometry, extraNodesPerSegment ) );
+  }
+}
+
+QgsGeometry QgsInternalGeometryEngine::densifyByDistance( double distance ) const
+{
+  if ( !mGeometry )
+  {
+    return QgsGeometry();
+  }
+
+  if ( QgsWkbTypes::geometryType( mGeometry->wkbType() ) == QgsWkbTypes::PointGeometry )
+  {
+    return QgsGeometry( mGeometry->clone() ); // point geometry, nothing to do
+  }
+
+  if ( const QgsGeometryCollection *gc = dynamic_cast< const QgsGeometryCollection *>( mGeometry ) )
+  {
+    int numGeom = gc->numGeometries();
+    QList< QgsAbstractGeometry * > geometryList;
+    geometryList.reserve( numGeom );
+    for ( int i = 0; i < numGeom; ++i )
+    {
+      geometryList << densifyGeometry( gc->geometryN( i ), -1, distance );
     }
 
     QgsGeometry first = QgsGeometry( geometryList.takeAt( 0 ) );
@@ -604,6 +650,6 @@ QgsGeometry QgsInternalGeometryEngine::densifyByCount( int extraNodesPerSegment
   }
   else
   {
-    return QgsGeometry( densifyGeometryByCount( mGeometry, extraNodesPerSegment ) );
+    return QgsGeometry( densifyGeometry( mGeometry, -1, distance ) );
   }
 }

src/core/geometry/qgsinternalgeometryengine.h

@@ -81,6 +81,21 @@ class QgsInternalGeometryEngine
      */
     QgsGeometry densifyByCount( int extraNodesPerSegment ) const;
 
+    /**
+     * Densifies the geometry by adding regularly placed extra nodes inside each segment
+     * so that the maximum distance between any two nodes does not exceed the
+     * specified \a distance.
+     * E.g. specifying a distance 3 would cause the segment [0 0] -> [10 0]
+     * to be converted to [0 0] -> [2.5 0] -> [5 0] -> [7.5 0] -> [10 0], since
+     * 3 extra nodes are required on the segment and spacing these at 2.5 increments
+     * allows them to be evenly spaced over the segment.
+     * If the geometry has z or m values present then these will be linearly interpolated
+     * at the added nodes.
+     * Curved geometry types are automatically segmentized by this routine.
+     * @node added in QGIS 3.0
+     */
+    QgsGeometry densifyByDistance( double distance ) const;
+
   private:
     const QgsAbstractGeometry *mGeometry = nullptr;
 };

tests/src/python/test_qgsgeometry.py

@@ -3996,5 +3996,73 @@ def testDensifyByCount(self):
         self.assertTrue(compareWkt(result, exp, 0.00001),
                         "densify by count: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result))
 
+    def testDensifyByDistance(self):
+        empty = QgsGeometry()
+        o = empty.densifyByDistance(4)
+        self.assertFalse(o)
+
+        # point
+        input = QgsGeometry.fromWkt("PointZ( 1 2 3 )")
+        o = input.densifyByDistance(0.1)
+        exp = "PointZ( 1 2 3 )"
+        result = o.exportToWkt()
+        self.assertTrue(compareWkt(result, exp, 0.00001),
+                        "densify by count: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result))
+        input = QgsGeometry.fromWkt(
+            "MULTIPOINT ((155 271), (150 360), (260 360), (271 265), (280 260), (270 370), (154 354), (150 260))")
+        o = input.densifyByDistance(0.1)
+        exp = "MULTIPOINT ((155 271), (150 360), (260 360), (271 265), (280 260), (270 370), (154 354), (150 260))"
+        result = o.exportToWkt()
+        self.assertTrue(compareWkt(result, exp, 0.00001),
+                        "densify by count: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result))
+
+        # line
+        input = QgsGeometry.fromWkt("LineString( 0 0, 10 0, 10 10 )")
+        o = input.densifyByDistance(100)
+        exp = "LineString( 0 0, 10 0, 10 10 )"
+        result = o.exportToWkt()
+        self.assertTrue(compareWkt(result, exp, 0.00001),
+                        "densify by count: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result))
+        o = input.densifyByDistance(3)
+        exp = "LineString (0 0, 2.5 0, 5 0, 7.5 0, 10 0, 10 2.5, 10 5, 10 7.5, 10 10)"
+
+        result = o.exportToWkt()
+        self.assertTrue(compareWkt(result, exp, 0.00001),
+                        "densify by count: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result))
+        input = QgsGeometry.fromWkt("LineStringZ( 0 0 1, 10 0 2, 10 10 0)")
+        o = input.densifyByDistance(6)
+        exp = "LineStringZ (0 0 1, 5 0 1.5, 10 0 2, 10 5 1, 10 10 0)"
+        result = o.exportToWkt()
+        self.assertTrue(compareWkt(result, exp, 0.00001),
+                        "densify by count: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result))
+        input = QgsGeometry.fromWkt("LineStringM( 0 0 0, 10 0 2, 10 10 0)")
+        o = input.densifyByDistance(3)
+        exp = "LineStringM (0 0 0, 2.5 0 0.5, 5 0 1, 7.5 0 1.5, 10 0 2, 10 2.5 1.5, 10 5 1, 10 7.5 0.5, 10 10 0)"
+        result = o.exportToWkt()
+        self.assertTrue(compareWkt(result, exp, 0.00001),
+                        "densify by count: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result))
+        input = QgsGeometry.fromWkt("LineStringZM( 0 0 1 10, 10 0 2 8, 10 10 0 4)")
+        o = input.densifyByDistance(6)
+        exp = "LineStringZM (0 0 1 10, 5 0 1.5 9, 10 0 2 8, 10 5 1 6, 10 10 0 4)"
+        result = o.exportToWkt()
+        self.assertTrue(compareWkt(result, exp, 0.00001),
+                        "densify by count: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result))
+
+        # polygon
+        input = QgsGeometry.fromWkt("Polygon(( 0 0, 20 0, 20 20, 0 0 ))")
+        o = input.densifyByDistance(110)
+        exp = "Polygon(( 0 0, 20 0, 20 20, 0 0 ))"
+        result = o.exportToWkt()
+        self.assertTrue(compareWkt(result, exp, 0.00001),
+                        "densify by count: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result))
+
+        input = QgsGeometry.fromWkt("PolygonZ(( 0 0 1, 20 0 2, 20 20 0, 0 0 1 ))")
+        o = input.densifyByDistance(6)
+        exp = "PolygonZ ((0 0 1, 5 0 1.25, 10 0 1.5, 15 0 1.75, 20 0 2, 20 5 1.5, 20 10 1, 20 15 0.5, 20 20 0, 16 16 0.2, 12 12 0.4, 8 8 0.6, 4 4 0.8, 0 0 1))"
+        result = o.exportToWkt()
+        self.assertTrue(compareWkt(result, exp, 0.00001),
+                        "densify by count: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result))
+
+
 if __name__ == '__main__':
     unittest.main()