diff --git a/src/core/geometry/qgslinestring.cpp b/src/core/geometry/qgslinestring.cpp
index 711314d8a5a44..a290a13650323 100644
--- a/src/core/geometry/qgslinestring.cpp
+++ b/src/core/geometry/qgslinestring.cpp
@@ -1086,6 +1086,37 @@ std::tuple<std::unique_ptr<QgsCurve>, std::unique_ptr<QgsCurve> > QgsLineString:
     return std::make_tuple( std::make_unique< QgsLineString >( x1, y1, z1, m1 ), std::make_unique< QgsLineString >( x2, y2, z2, m2 ) );
 }
 
+QVector<QgsLineString> QgsLineString::splitToDisjointParts() const
+{
+  const double *allPointsX = xData();
+  const double *allPointsY = yData();
+  size_t allPointsCount = numPoints();
+  QgsPointSequence partPoints;
+  QSet<QgsPoint> partPointSet;
+
+  QVector<QgsLineString> disjointParts;
+  for ( size_t i = 0; i < allPointsCount; i++ )
+  {
+    QgsPoint point( *allPointsX++, *allPointsY++ );
+    if ( partPointSet.contains( point ) )
+    {
+      // This point is used multiple times, cut the curve and add the
+      // current part
+      disjointParts.append( QgsLineString( partPoints ) );
+      // Now start a new part containing the last line
+      partPoints = { partPoints.last() };
+      partPointSet = { partPoints[0] };
+    }
+    partPoints.push_back( point );
+    partPointSet.insert( point );
+  }
+  // Add the last part (if we didn't stop by closing the loop)
+  if ( partPoints.size() > 1 || disjointParts.size() == 0 )
+    disjointParts.append( QgsLineString( partPoints ) );
+
+  return disjointParts;
+}
+
 double QgsLineString::length3D() const
 {
   if ( is3D() )
diff --git a/src/core/geometry/qgslinestring.h b/src/core/geometry/qgslinestring.h
index 1fe960eb924eb..54bd1e5c55f86 100644
--- a/src/core/geometry/qgslinestring.h
+++ b/src/core/geometry/qgslinestring.h
@@ -991,6 +991,12 @@ class CORE_EXPORT QgsLineString: public QgsCurve
     std::tuple< std::unique_ptr< QgsCurve >, std::unique_ptr< QgsCurve > > splitCurveAtVertex( int index ) const final;
 #endif
 
+    /**
+     * Divides the linestring into parts that don't share any points or lines.
+     * \since QGIS 3.40
+     */
+    QVector<QgsLineString> splitToDisjointParts() const;
+
     /**
      * Returns the length in 3D world of the line string.
      * If it is not a 3D line string, return its 2D length.
diff --git a/src/core/vector/qgsvectorlayerprofilegenerator.cpp b/src/core/vector/qgsvectorlayerprofilegenerator.cpp
index a77d368089bf1..a2a911ce86171 100644
--- a/src/core/vector/qgsvectorlayerprofilegenerator.cpp
+++ b/src/core/vector/qgsvectorlayerprofilegenerator.cpp
@@ -15,6 +15,7 @@
  *                                                                         *
  ***************************************************************************/
 #include "qgsvectorlayerprofilegenerator.h"
+#include "qgsabstractgeometry.h"
 #include "qgspolyhedralsurface.h"
 #include "qgsprofilerequest.h"
 #include "qgscurve.h"
@@ -721,6 +722,65 @@ bool QgsVectorLayerProfileGenerator::generateProfile( const QgsProfileGeneration
   if ( !mProfileCurve || mFeedback->isCanceled() )
     return false;
 
+  auto profileLine = qgsgeometry_cast<QgsLineString *>( mProfileCurve.get() );
+  if ( profileLine )
+  {
+    // The profile generation code can't deal with curves that enter a single
+    // point multiple times. We handle this for line strings by splitting them
+    // into multiple parts, each with no repeated points, and computing the
+    // profile for each by itself.
+    QgsCurve *origCurve = mProfileCurve.release();
+    std::unique_ptr< QgsVectorLayerProfileResults > totalResults;
+    double distanceProcessed = 0;
+
+    for ( QgsLineString &part : profileLine->splitToDisjointParts() )
+    {
+      mProfileCurve.reset( &part );
+      if ( !generateProfileInner() ) return false;
+
+      if ( !totalResults )
+        // Use the first result set as a base
+        totalResults.reset( mResults.release() );
+      else
+      {
+        // Merge the results, shifting them by distanceProcessed
+        totalResults->mRawPoints.append( mResults->mRawPoints );
+        totalResults->minZ = std::min( totalResults->minZ, mResults->minZ );
+        totalResults->maxZ = std::max( totalResults->maxZ, mResults->maxZ );
+        for ( auto it = mResults->mDistanceToHeightMap.constKeyValueBegin();
+              it != mResults->mDistanceToHeightMap.constKeyValueEnd();
+              ++it )
+        {
+          totalResults->mDistanceToHeightMap[it->first + distanceProcessed] = it->second;
+        }
+        for ( auto it = mResults->features.constKeyValueBegin();
+              it != mResults->features.constKeyValueEnd();
+              ++it )
+        {
+          for ( QgsVectorLayerProfileResults::Feature feature : it->second )
+          {
+            feature.crossSectionGeometry.translate( distanceProcessed, 0 );
+            totalResults->features[it->first].push_back( feature );
+          }
+        }
+      }
+
+      distanceProcessed += mProfileCurve->length();
+      // We filled mProfileCurve with a reference to part of a vector. We now
+      // need to get it out without free()ing the memory, since we don't own it.
+      Q_UNUSED( mProfileCurve.release() );
+    }
+
+    mProfileCurve.reset( origCurve );
+    mResults.reset( totalResults.release() );
+    return true;
+  }
+
+  return generateProfileInner();
+}
+
+bool QgsVectorLayerProfileGenerator::generateProfileInner( const QgsProfileGenerationContext & )
+{
   // we need to transform the profile curve to the vector's CRS
   mTransformedCurve.reset( mProfileCurve->clone() );
   mLayerToTargetTransform = QgsCoordinateTransform( mSourceCrs, mTargetCrs, mTransformContext );
diff --git a/src/core/vector/qgsvectorlayerprofilegenerator.h b/src/core/vector/qgsvectorlayerprofilegenerator.h
index 17ae6130298cf..8bc6478011075 100644
--- a/src/core/vector/qgsvectorlayerprofilegenerator.h
+++ b/src/core/vector/qgsvectorlayerprofilegenerator.h
@@ -120,6 +120,10 @@ class CORE_EXPORT QgsVectorLayerProfileGenerator : public QgsAbstractProfileSurf
 
   private:
 
+    // We may need to split mProfileCurve into multiple parts, this will be
+    // called for each part.
+    bool generateProfileInner( const QgsProfileGenerationContext &context = QgsProfileGenerationContext() );
+
     bool generateProfileForPoints();
     bool generateProfileForLines();
     bool generateProfileForPolygons();
diff --git a/tests/src/core/geometry/testqgsgeometry.cpp b/tests/src/core/geometry/testqgsgeometry.cpp
index 8818388ddca04..f31c5d4386681 100644
--- a/tests/src/core/geometry/testqgsgeometry.cpp
+++ b/tests/src/core/geometry/testqgsgeometry.cpp
@@ -89,6 +89,7 @@ class TestQgsGeometry : public QgsTest
     void curveIndexOf();
     void splitCurve_data();
     void splitCurve();
+    void splitToDisjointParts();
 
     void fromBox3d();
     void fromPoint();
@@ -751,6 +752,20 @@ void TestQgsGeometry::splitCurve()
   QCOMPARE( p2->asWkt(), curve2 );
 }
 
+void TestQgsGeometry::splitToDisjointParts()
+{
+  QgsLineString onePartLine( QgsPoint( 1.0, 1.0 ), QgsPoint( 2.0, 2.0 ) );
+  QVector<QgsLineString> onePartParts = onePartLine.splitToDisjointParts();
+  QCOMPARE( onePartParts.size(), 1 );
+  QCOMPARE( onePartParts[0].asWkt(), onePartLine.asWkt() );
+
+  QgsLineString triangle( QVector {QgsPoint( 0.0, 0.0 ), QgsPoint( 1.0, 0.0 ), QgsPoint( 1.0, 1.0 ), QgsPoint( 0.0, 0.0 )} );
+  QVector<QgsLineString> triangleParts = triangle.splitToDisjointParts();
+  QCOMPARE( triangleParts.size(), 2 );
+  QCOMPARE( triangleParts[0].asWkt(), "" );
+  QCOMPARE( triangleParts[1].asWkt(), "" );
+}
+
 void TestQgsGeometry::fromBox3d()
 {
   QgsBox3D box3d( 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 );