Add QgsLineString method to interpolate nan m values along line

Fills in any nan m values by interpolating from non-nan values
in neighbouring vertices

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

@@ -863,13 +863,22 @@ Calculates the minimal 3D bounding box for the geometry.
 .. versionadded:: 3.34
 %End
 
-
     QgsLineString *measuredLine( double start, double end ) const /Factory/;
 %Docstring
 Re-write the measure ordinate (or add one, if it isn't already there) interpolating
 the measure between the supplied ``start`` and ``end`` values.
 
 .. versionadded:: 3.36
+%End
+
+    QgsLineString *interpolateM() const /Factory/;
+%Docstring
+Returns a copy of this line with all missing (nan) m values interpolated
+from m values of surrounding vertices.
+
+If the line does not contain m values, ``None`` is returned.
+
+.. versionadded:: 3.38
 %End
 
   protected:

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

@@ -863,13 +863,22 @@ Calculates the minimal 3D bounding box for the geometry.
 .. versionadded:: 3.34
 %End
 
-
     QgsLineString *measuredLine( double start, double end ) const /Factory/;
 %Docstring
 Re-write the measure ordinate (or add one, if it isn't already there) interpolating
 the measure between the supplied ``start`` and ``end`` values.
 
 .. versionadded:: 3.36
+%End
+
+    QgsLineString *interpolateM() const /Factory/;
+%Docstring
+Returns a copy of this line with all missing (nan) m values interpolated
+from m values of surrounding vertices.
+
+If the line does not contain m values, ``None`` is returned.
+
+.. versionadded:: 3.38
 %End
 
   protected:

src/core/geometry/qgslinestring.cpp

@@ -2387,3 +2387,166 @@ QgsLineString *QgsLineString::measuredLine( double start, double end ) const
 
   return cloned.release();
 }
+
+QgsLineString *QgsLineString::interpolateM() const
+{
+  if ( !isMeasure() )
+    return nullptr;
+
+  const int totalPoints = numPoints();
+  if ( totalPoints < 2 )
+    return clone();
+
+  const double *xData = mX.constData();
+  const double *yData = mY.constData();
+  const double *mData = mM.constData();
+  const double *zData = is3D() ? mZ.constData() : nullptr;
+
+  QVector< double > xOut( totalPoints );
+  QVector< double > yOut( totalPoints );
+  QVector< double > mOut( totalPoints );
+  QVector< double > zOut( is3D() ? totalPoints : 0 );
+
+  double *xOutData = xOut.data();
+  double *yOutData = yOut.data();
+  double *mOutData = mOut.data();
+  double *zOutData = is3D() ? zOut.data() : nullptr;
+
+  int i = 0;
+  double currentSegmentLength = 0;
+  double lastValidM = std::numeric_limits< double >::quiet_NaN();
+  double prevX = *xData;
+  double prevY = *yData;
+  while ( i < totalPoints )
+  {
+    double thisX = *xData++;
+    double thisY = *yData++;
+    double thisZ = zData ? *zData++ : 0;
+    double thisM = *mData++;
+
+    currentSegmentLength = QgsGeometryUtilsBase::distance2D( prevX, prevY, thisX, thisY );
+
+    if ( !std::isnan( thisM ) )
+    {
+      *xOutData++ = thisX;
+      *yOutData++ = thisY;
+      *mOutData++ = thisM;
+      if ( zOutData )
+        *zOutData++ = thisZ;
+      lastValidM = thisM;
+    }
+    else if ( i == 0 )
+    {
+      // nan m value at start of line, read ahead to find first non-nan value and backfill
+      int j = 0;
+      double scanAheadM = thisM;
+      while ( i + j + 1 < totalPoints && std::isnan( scanAheadM ) )
+      {
+        scanAheadM = mData[ j ];
+        ++j;
+      }
+      if ( std::isnan( scanAheadM ) )
+      {
+        // no valid m values in line
+        return nullptr;
+      }
+      *xOutData++ = thisX;
+      *yOutData++ = thisY;
+      *mOutData++ = scanAheadM;
+      if ( zOutData )
+        *zOutData++ = thisZ;
+      for ( ; i < j; ++i )
+      {
+        thisX = *xData++;
+        thisY = *yData++;
+        *xOutData++ = thisX;
+        *yOutData++ = thisY;
+        *mOutData++ = scanAheadM;
+        mData++;
+        if ( zOutData )
+          *zOutData++ = *zData++;
+      }
+      lastValidM = scanAheadM;
+    }
+    else
+    {
+      // nan m value in middle of line, read ahead till next non-nan value and interpolate
+      int j = 0;
+      double scanAheadX = thisX;
+      double scanAheadY = thisY;
+      double distanceToNextValidM = currentSegmentLength;
+      std::vector< double > scanAheadSegmentLengths;
+      scanAheadSegmentLengths.emplace_back( currentSegmentLength );
+
+      double nextValidM = std::numeric_limits< double >::quiet_NaN();
+      while ( i + j < totalPoints - 1 )
+      {
+        double nextScanAheadX = xData[j];
+        double nextScanAheadY = yData[j];
+        double nextScanAheadM = mData[ j ];
+        const double scanAheadSegmentLength = QgsGeometryUtilsBase::distance2D( scanAheadX, scanAheadY, nextScanAheadX, nextScanAheadY );
+        scanAheadSegmentLengths.emplace_back( scanAheadSegmentLength );
+        distanceToNextValidM += scanAheadSegmentLength;
+
+        if ( !std::isnan( nextScanAheadM ) )
+        {
+          nextValidM = nextScanAheadM;
+          break;
+        }
+
+        scanAheadX = nextScanAheadX;
+        scanAheadY = nextScanAheadY;
+        ++j;
+      }
+
+      if ( std::isnan( nextValidM ) )
+      {
+        // no more valid m values, so just fill remainder of vertices with previous valid m value
+        *xOutData++ = thisX;
+        *yOutData++ = thisY;
+        *mOutData++ = lastValidM;
+        if ( zOutData )
+          *zOutData++ = thisZ;
+        ++i;
+        for ( ; i < totalPoints; ++i )
+        {
+          *xOutData++ = *xData++;
+          *yOutData++ = *yData++;
+          *mOutData++ = lastValidM;
+          if ( zOutData )
+            *zOutData++ = *zData++;
+        }
+        break;
+      }
+      else
+      {
+        // interpolate along segments
+        const double delta = ( nextValidM - lastValidM ) / distanceToNextValidM;
+        *xOutData++ = thisX;
+        *yOutData++ = thisY;
+        *mOutData++ = lastValidM + delta * scanAheadSegmentLengths[0];
+        double totalScanAheadLength = scanAheadSegmentLengths[0];
+        if ( zOutData )
+          *zOutData++ = thisZ;
+        for ( int k = 1; k <= j; ++i, ++k )
+        {
+          thisX = *xData++;
+          thisY = *yData++;
+          *xOutData++ = thisX;
+          *yOutData++ = thisY;
+          totalScanAheadLength += scanAheadSegmentLengths[k];
+          *mOutData++ = lastValidM + delta * totalScanAheadLength;
+          mData++;
+          if ( zOutData )
+            *zOutData++ = *zData++;
+        }
+        lastValidM = nextValidM;
+      }
+    }
+
+    prevX = thisX;
+    prevY = thisY;
+    ++i;
+  }
+  return new QgsLineString( xOut, yOut, zOut, mOut );
+}

src/core/geometry/qgslinestring.h

@@ -1187,7 +1187,6 @@ class CORE_EXPORT QgsLineString: public QgsCurve
      */
     QgsBox3D calculateBoundingBox3D() const override;
 
-
     /**
      * Re-write the measure ordinate (or add one, if it isn't already there) interpolating
      * the measure between the supplied \a start and \a end values.
@@ -1196,6 +1195,16 @@ class CORE_EXPORT QgsLineString: public QgsCurve
      */
     QgsLineString *measuredLine( double start, double end ) const SIP_FACTORY;
 
+    /**
+     * Returns a copy of this line with all missing (nan) m values interpolated
+     * from m values of surrounding vertices.
+     *
+     * If the line does not contain m values, NULLPTR is returned.
+     *
+     * \since QGIS 3.38
+     */
+    QgsLineString *interpolateM() const SIP_FACTORY;
+
   protected:
 
     int compareToSameClass( const QgsAbstractGeometry *other ) const final;