[FEATURE] Use native interpolate point method instead of GEOS method

Because:
- Exactly follows curves and doesn't require segmentizing input geometry
- Also interpolates z/m values if they are present in input geometry
- Is faster

Author: nyalldawson

python/core/auto_generated/geometry/qgscircularstring.sip.in

@@ -143,6 +143,8 @@ Sets the circular string's points
 
     virtual QgsCircularString *reversed() const  /Factory/;
 
+    virtual QgsPoint *interpolatePoint( double distance ) const /Factory/;
+
     virtual QgsCircularString *curveSubstring( double startDistance, double endDistance ) const /Factory/;
 
     virtual bool addZValue( double zValue = 0 );

python/core/auto_generated/geometry/qgscompoundcurve.sip.in

@@ -144,6 +144,8 @@ Appends first point if not already closed.
 
     virtual QgsCompoundCurve *reversed() const /Factory/;
 
+    virtual QgsPoint *interpolatePoint( double distance ) const /Factory/;
+
     virtual QgsCompoundCurve *curveSubstring( double startDistance, double endDistance ) const /Factory/;
 
 

python/core/auto_generated/geometry/qgscurve.sip.in

@@ -180,6 +180,19 @@ Returns the y-coordinate of the specified node in the line string.
     virtual QPolygonF asQPolygonF() const;
 %Docstring
 Returns a QPolygonF representing the points.
+%End
+
+    virtual QgsPoint *interpolatePoint( double distance ) const = 0 /Factory/;
+%Docstring
+Returns an interpolated point on the curve at the specified ``distance``.
+
+If z or m values are present, the output z and m will be interpolated using
+the existing vertices' z or m values.
+
+If distance is negative, or is greater than the length of the curve, a None
+will be returned.
+
+.. versionadded:: 3.4
 %End
 
     virtual QgsCurve *curveSubstring( double startDistance, double endDistance ) const = 0 /Factory/;

python/core/auto_generated/geometry/qgsgeometry.sip.in

@@ -1154,12 +1154,16 @@ by calling `error()` on the returned geometry.
 
     QgsGeometry interpolate( double distance ) const;
 %Docstring
-Returns interpolated point on line at distance.
+Returns an interpolated point on the geometry at the specified ``distance``.
 
-If the input is a NULL geometry, the output will also be a NULL geometry.
+If the original geometry is a polygon type, the boundary of the polygon
+will be used during interpolation. If the original geometry is a point
+type, a null geometry will be returned.
 
-If an error was encountered while creating the result, more information can be retrieved
-by calling `error()` on the returned geometry.
+If z or m values are present, the output z and m will be interpolated using
+the existing vertices' z or m values.
+
+If the input is a NULL geometry, the output will also be a NULL geometry.
 
 .. seealso:: :py:func:`lineLocatePoint`
 

python/core/auto_generated/geometry/qgslinestring.sip.in

@@ -283,6 +283,8 @@ of the curve.
 
     virtual QgsLineString *reversed() const /Factory/;
 
+    virtual QgsPoint *interpolatePoint( double distance ) const /Factory/;
+
     virtual QgsLineString *curveSubstring( double startDistance, double endDistance ) const /Factory/;
 
 

python/plugins/processing/tests/testdata/custom/circular_strings.gpkg

@@ -71,9 +71,14 @@ QgsProcessing::SourceType QgsInterpolatePointAlgorithm::outputLayerType() const
   return QgsProcessing::TypeVectorPoint;
 }
 
-QgsWkbTypes::Type QgsInterpolatePointAlgorithm::outputWkbType( QgsWkbTypes::Type ) const
-{
-  return QgsWkbTypes::Point;
+QgsWkbTypes::Type QgsInterpolatePointAlgorithm::outputWkbType( QgsWkbTypes::Type inputType ) const
+{
+  QgsWkbTypes::Type out = QgsWkbTypes::Point;
+  if ( QgsWkbTypes::hasZ( inputType ) )
+    out = QgsWkbTypes::addZ( out );
+  if ( QgsWkbTypes::hasM( inputType ) )
+    out = QgsWkbTypes::addM( out );
+  return out;
 }
 
 QgsInterpolatePointAlgorithm *QgsInterpolatePointAlgorithm::createInstance() const

src/analysis/processing/qgsalgorithminterpolatepoint.cpp

@@ -1178,6 +1178,75 @@ QgsCircularString *QgsCircularString::reversed() const
   return copy;
 }
 
+QgsPoint *QgsCircularString::interpolatePoint( const double distance ) const
+{
+  if ( distance < 0 )
+    return nullptr;
+
+  double distanceTraversed = 0;
+  const int totalPoints = numPoints();
+  if ( totalPoints == 0 )
+    return nullptr;
+
+  QgsWkbTypes::Type pointType = QgsWkbTypes::Point;
+  if ( is3D() )
+    pointType = QgsWkbTypes::PointZ;
+  if ( isMeasure() )
+    pointType = QgsWkbTypes::addM( pointType );
+
+  const double *x = mX.constData();
+  const double *y = mY.constData();
+  const double *z = is3D() ? mZ.constData() : nullptr;
+  const double *m = isMeasure() ? mM.constData() : nullptr;
+
+  double prevX = *x++;
+  double prevY = *y++;
+  double prevZ = z ? *z++ : 0.0;
+  double prevM = m ? *m++ : 0.0;
+
+  if ( qgsDoubleNear( distance, 0.0 ) )
+  {
+    return new QgsPoint( pointType, prevX, prevY, prevZ, prevM );
+  }
+
+  for ( int i = 0; i < ( totalPoints - 2 ) ; i += 2 )
+  {
+    double x1 = prevX;
+    double y1 = prevY;
+    double z1 = prevZ;
+    double m1 = prevM;
+
+    double x2 = *x++;
+    double y2 = *y++;
+    double z2 = z ? *z++ : 0.0;
+    double m2 = m ? *m++ : 0.0;
+
+    double x3 = *x++;
+    double y3 = *y++;
+    double z3 = z ? *z++ : 0.0;
+    double m3 = m ? *m++ : 0.0;
+
+    const double segmentLength = QgsGeometryUtils::circleLength( x1, y1, x2, y2, x3, y3 );
+    if ( distance < distanceTraversed + segmentLength || qgsDoubleNear( distance, distanceTraversed + segmentLength ) )
+    {
+      // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
+      const double distanceToPoint = std::min( distance - distanceTraversed, segmentLength );
+      return new QgsPoint( QgsGeometryUtils::interpolatePointOnArc( QgsPoint( pointType, x1, y1, z1, m1 ),
+                           QgsPoint( pointType, x2, y2, z2, m2 ),
+                           QgsPoint( pointType, x3, y3, z3, m3 ), distanceToPoint ) );
+    }
+
+    distanceTraversed += segmentLength;
+
+    prevX = x3;
+    prevY = y3;
+    prevZ = z3;
+    prevM = m3;
+  }
+
+  return nullptr;
+}
+
 QgsCircularString *QgsCircularString::curveSubstring( double startDistance, double endDistance ) const
 {
   if ( startDistance < 0 && endDistance < 0 )

src/core/geometry/qgscircularstring.cpp

@@ -118,6 +118,7 @@ class CORE_EXPORT QgsCircularString: public QgsCurve
     double vertexAngle( QgsVertexId vertex ) const override;
     double segmentLength( QgsVertexId startVertex ) const override;
     QgsCircularString *reversed() const override  SIP_FACTORY;
+    QgsPoint *interpolatePoint( double distance ) const override SIP_FACTORY;
     QgsCircularString *curveSubstring( double startDistance, double endDistance ) const override SIP_FACTORY;
     bool addZValue( double zValue = 0 ) override;
     bool addMValue( double mValue = 0 ) override;

src/core/geometry/qgscircularstring.h

@@ -864,6 +864,30 @@ QgsCompoundCurve *QgsCompoundCurve::reversed() const
   return clone;
 }
 
+QgsPoint *QgsCompoundCurve::interpolatePoint( const double distance ) const
+{
+  if ( distance < 0 )
+    return nullptr;
+
+  double distanceTraversed = 0;
+  for ( const QgsCurve *curve : mCurves )
+  {
+    const double thisCurveLength = curve->length();
+    if ( distanceTraversed + thisCurveLength > distance || qgsDoubleNear( distanceTraversed + thisCurveLength, distance ) )
+    {
+      // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
+      const double distanceToPoint = std::min( distance - distanceTraversed, thisCurveLength );
+
+      // point falls on this curve
+      return curve->interpolatePoint( distanceToPoint );
+    }
+
+    distanceTraversed += thisCurveLength;
+  }
+
+  return nullptr;
+}
+
 QgsCompoundCurve *QgsCompoundCurve::curveSubstring( double startDistance, double endDistance ) const
 {
   if ( startDistance < 0 && endDistance < 0 )

src/core/geometry/qgscompoundcurve.cpp

@@ -115,6 +115,7 @@ class CORE_EXPORT QgsCompoundCurve: public QgsCurve
     double vertexAngle( QgsVertexId vertex ) const override;
     double segmentLength( QgsVertexId startVertex ) const override;
     QgsCompoundCurve *reversed() const override SIP_FACTORY;
+    QgsPoint *interpolatePoint( double distance ) const override SIP_FACTORY;
     QgsCompoundCurve *curveSubstring( double startDistance, double endDistance ) const override SIP_FACTORY;
 
     bool addZValue( double zValue = 0 ) override;

src/core/geometry/qgscompoundcurve.h

@@ -166,6 +166,19 @@ class CORE_EXPORT QgsCurve: public QgsAbstractGeometry
      */
     virtual QPolygonF asQPolygonF() const;
 
+    /**
+     * Returns an interpolated point on the curve at the specified \a distance.
+     *
+     * If z or m values are present, the output z and m will be interpolated using
+     * the existing vertices' z or m values.
+     *
+     * If distance is negative, or is greater than the length of the curve, a nullptr
+     * will be returned.
+     *
+     * \since QGIS 3.4
+     */
+    virtual QgsPoint *interpolatePoint( double distance ) const = 0 SIP_FACTORY;
+
     /**
      * Returns a new curve representing a substring of this curve.
      *

src/core/geometry/qgscurve.h

@@ -45,6 +45,7 @@ email                : morb at ozemail dot com dot au
 #include "qgspolygon.h"
 #include "qgslinestring.h"
 #include "qgscircle.h"
+#include "qgscurve.h"
 
 struct QgsGeometryPrivate
 {
@@ -2005,25 +2006,42 @@ QgsGeometry QgsGeometry::subdivide( int maxNodes ) const
   return QgsGeometry( std::move( result ) );
 }
 
-QgsGeometry QgsGeometry::interpolate( double distance ) const
+QgsGeometry QgsGeometry::interpolate( const double distance ) const
 {
   if ( !d->geometry )
   {
     return QgsGeometry();
   }
 
   QgsGeometry line = *this;
-  if ( type() == QgsWkbTypes::PolygonGeometry )
+  if ( type() == QgsWkbTypes::PointGeometry )
+    return QgsGeometry();
+  else if ( type() == QgsWkbTypes::PolygonGeometry )
+  {
     line = QgsGeometry( d->geometry->boundary() );
+  }
 
-  QgsGeos geos( line.constGet() );
-  mLastError.clear();
-  std::unique_ptr< QgsAbstractGeometry > result( geos.interpolate( distance, &mLastError ) );
+  const QgsCurve *curve = nullptr;
+  if ( line.isMultipart() )
+  {
+    // if multi part, just use first part
+    const QgsGeometryCollection *collection = qgsgeometry_cast< const QgsGeometryCollection * >( line.constGet() );
+    if ( collection && collection->numGeometries() > 0 )
+    {
+      curve = qgsgeometry_cast< const QgsCurve * >( collection->geometryN( 0 ) );
+    }
+  }
+  else
+  {
+    curve = qgsgeometry_cast< const QgsCurve * >( line.constGet() );
+  }
+  if ( !curve )
+    return QgsGeometry();
+
+  std::unique_ptr< QgsPoint > result( curve->interpolatePoint( distance ) );
   if ( !result )
   {
-    QgsGeometry geom;
-    geom.mLastError = mLastError;
-    return geom;
+    return QgsGeometry();
   }
   return QgsGeometry( std::move( result ) );
 }

src/core/geometry/qgsgeometry.cpp

@@ -1134,12 +1134,16 @@ class CORE_EXPORT QgsGeometry
     QgsGeometry subdivide( int maxNodes = 256 ) const;
 
     /**
-     * Returns interpolated point on line at distance.
+     * Returns an interpolated point on the geometry at the specified \a distance.
      *
-     * If the input is a NULL geometry, the output will also be a NULL geometry.
+     * If the original geometry is a polygon type, the boundary of the polygon
+     * will be used during interpolation. If the original geometry is a point
+     * type, a null geometry will be returned.
      *
-     * If an error was encountered while creating the result, more information can be retrieved
-     * by calling `error()` on the returned geometry.
+     * If z or m values are present, the output z and m will be interpolated using
+     * the existing vertices' z or m values.
+     *
+     * If the input is a NULL geometry, the output will also be a NULL geometry.
      *
      * \see lineLocatePoint()
      * \since QGIS 2.0

src/core/geometry/qgsgeometry.h

@@ -711,6 +711,68 @@ QgsLineString *QgsLineString::reversed() const
   return copy;
 }
 
+QgsPoint *QgsLineString::interpolatePoint( const double distance ) const
+{
+  if ( distance < 0 )
+    return nullptr;
+
+  double distanceTraversed = 0;
+  const int totalPoints = numPoints();
+  if ( totalPoints == 0 )
+    return nullptr;
+
+  const double *x = mX.constData();
+  const double *y = mY.constData();
+  const double *z = is3D() ? mZ.constData() : nullptr;
+  const double *m = isMeasure() ? mM.constData() : nullptr;
+
+  QgsWkbTypes::Type pointType = QgsWkbTypes::Point;
+  if ( is3D() )
+    pointType = QgsWkbTypes::PointZ;
+  if ( isMeasure() )
+    pointType = QgsWkbTypes::addM( pointType );
+
+  double prevX = *x++;
+  double prevY = *y++;
+  double prevZ = z ? *z++ : 0.0;
+  double prevM = m ? *m++ : 0.0;
+
+  if ( qgsDoubleNear( distance, 0.0 ) )
+  {
+    return new QgsPoint( pointType, prevX, prevY, prevZ, prevM );
+  }
+
+  for ( int i = 1; i < totalPoints; ++i )
+  {
+    double thisX = *x++;
+    double thisY = *y++;
+    double thisZ = z ? *z++ : 0.0;
+    double thisM = m ? *m++ : 0.0;
+
+    const double segmentLength = std::sqrt( ( thisX - prevX ) * ( thisX - prevX ) + ( thisY - prevY ) * ( thisY - prevY ) );
+    if ( distance < distanceTraversed + segmentLength || qgsDoubleNear( distance, distanceTraversed + segmentLength ) )
+    {
+      // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
+      const double distanceToPoint = std::min( distance - distanceTraversed, segmentLength );
+      double pX, pY;
+      double pZ = 0;
+      double pM = 0;
+      QgsGeometryUtils::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY,
+          z ? &prevZ : nullptr, z ? &thisZ : nullptr, z ? &pZ : nullptr,
+          m ? &prevM : nullptr, m ? &thisM : nullptr, m ? &pM : nullptr );
+      return new QgsPoint( pointType, pX, pY, pZ, pM );
+    }
+
+    distanceTraversed += segmentLength;
+    prevX = thisX;
+    prevY = thisY;
+    prevZ = thisZ;
+    prevM = thisM;
+  }
+
+  return nullptr;
+}
+
 QgsLineString *QgsLineString::curveSubstring( double startDistance, double endDistance ) const
 {
   if ( startDistance < 0 && endDistance < 0 )

src/core/geometry/qgslinestring.cpp

@@ -308,6 +308,7 @@ class CORE_EXPORT QgsLineString: public QgsCurve
     bool deleteVertex( QgsVertexId position ) override;
 
     QgsLineString *reversed() const override SIP_FACTORY;
+    QgsPoint *interpolatePoint( double distance ) const override SIP_FACTORY;
     QgsLineString *curveSubstring( double startDistance, double endDistance ) const override SIP_FACTORY;
 
     double closestSegment( const QgsPoint &pt, QgsPoint &segmentPt SIP_OUT, QgsVertexId &vertexAfter SIP_OUT, int *leftOf SIP_OUT = nullptr, double epsilon = 4 * std::numeric_limits<double>::epsilon() ) const override;