diff --git a/python/core/auto_generated/geometry/qgslinestring.sip.in b/python/core/auto_generated/geometry/qgslinestring.sip.in
index d5934d372686..16a816b5b381 100644
--- a/python/core/auto_generated/geometry/qgslinestring.sip.in
+++ b/python/core/auto_generated/geometry/qgslinestring.sip.in
@@ -76,6 +76,18 @@ or repeatedly calling addVertex()
 Construct a linestring from a single 2d line segment.
 
 .. versionadded:: 3.2
+%End
+
+    static QgsLineString *fromBezierCurve( const QgsPoint &start, const QgsPoint &controlPoint1, const QgsPoint &controlPoint2, const QgsPoint &end, int segments = 30 ) /Factory/;
+%Docstring
+Returns a new linestring created by segmentizing the bezier curve between ``start`` and ``end``, with
+the specified control points.
+
+The ``segments`` parameter controls how many line segments will be present in the returned linestring.
+
+Any z or m values present in the input coordinates will be interpolated along with the x and y values.
+
+.. versionadded:: 3.10
 %End
 
     virtual bool equals( const QgsCurve &other ) const;
diff --git a/src/core/geometry/qgslinestring.cpp b/src/core/geometry/qgslinestring.cpp
index 65c9c8d62f07..dd9a8deeb0d3 100644
--- a/src/core/geometry/qgslinestring.cpp
+++ b/src/core/geometry/qgslinestring.cpp
@@ -175,6 +175,73 @@ QgsLineString::QgsLineString( const QgsLineSegment2D &segment )
   mY[1] = segment.endY();
 }
 
+static double cubicInterpolate( double a, double b,
+                                double A, double B, double C, double D )
+{
+  return A * b * b * b + 3 * B * b * b * a + 3 * C * b * a * a + D * a * a * a;
+}
+
+QgsLineString *QgsLineString::fromBezierCurve( const QgsPoint &start, const QgsPoint &controlPoint1, const QgsPoint &controlPoint2, const QgsPoint &end, int segments )
+{
+  if ( segments == 0 )
+    return new QgsLineString();
+
+  QVector<double> x;
+  x.resize( segments + 1 );
+  QVector<double> y;
+  y.resize( segments + 1 );
+  QVector<double> z;
+  double *zData = nullptr;
+  if ( start.is3D() && end.is3D() && controlPoint1.is3D() && controlPoint2.is3D() )
+  {
+    z.resize( segments + 1 );
+    zData = z.data();
+  }
+  QVector<double> m;
+  double *mData = nullptr;
+  if ( start.isMeasure() && end.isMeasure() && controlPoint1.isMeasure() && controlPoint2.isMeasure() )
+  {
+    m.resize( segments + 1 );
+    mData = m.data();
+  }
+
+  double *xData = x.data();
+  double *yData = y.data();
+  const double step = 1.0 / segments;
+  double a = 0;
+  double b = 1.0;
+  for ( int i = 0; i < segments; i++, a += step, b -= step )
+  {
+    if ( i == 0 )
+    {
+      *xData++ = start.x();
+      *yData++ = start.y();
+      if ( zData )
+        *zData++ = start.z();
+      if ( mData )
+        *mData++ = start.m();
+    }
+    else
+    {
+      *xData++ = cubicInterpolate( a, b, start.x(), controlPoint1.x(), controlPoint2.x(), end.x() );
+      *yData++ = cubicInterpolate( a, b, start.y(), controlPoint1.y(), controlPoint2.y(), end.y() );
+      if ( zData )
+        *zData++ = cubicInterpolate( a, b, start.z(), controlPoint1.z(), controlPoint2.z(), end.z() );
+      if ( mData )
+        *mData++ = cubicInterpolate( a, b, start.m(), controlPoint1.m(), controlPoint2.m(), end.m() );
+    }
+  }
+
+  *xData = end.x();
+  *yData = end.y();
+  if ( zData )
+    *zData = end.z();
+  if ( mData )
+    *mData = end.m();
+
+  return new QgsLineString( x, y, z, m );
+}
+
 bool QgsLineString::equals( const QgsCurve &other ) const
 {
   const QgsLineString *otherLine = qgsgeometry_cast< const QgsLineString * >( &other );
diff --git a/src/core/geometry/qgslinestring.h b/src/core/geometry/qgslinestring.h
index b9ea127616e7..518634b59bb4 100644
--- a/src/core/geometry/qgslinestring.h
+++ b/src/core/geometry/qgslinestring.h
@@ -92,6 +92,18 @@ class CORE_EXPORT QgsLineString: public QgsCurve
      */
     explicit QgsLineString( const QgsLineSegment2D &segment );
 
+    /**
+     * Returns a new linestring created by segmentizing the bezier curve between \a start and \a end, with
+     * the specified control points.
+     *
+     * The \a segments parameter controls how many line segments will be present in the returned linestring.
+     *
+     * Any z or m values present in the input coordinates will be interpolated along with the x and y values.
+     *
+     * \since QGIS 3.10
+     */
+    static QgsLineString *fromBezierCurve( const QgsPoint &start, const QgsPoint &controlPoint1, const QgsPoint &controlPoint2, const QgsPoint &end, int segments = 30 ) SIP_FACTORY;
+
     bool equals( const QgsCurve &other ) const override;
 
 #ifndef SIP_RUN
diff --git a/tests/src/python/test_qgsgeometry.py b/tests/src/python/test_qgsgeometry.py
index 49e8cfb1b380..1dbf7baf0156 100644
--- a/tests/src/python/test_qgsgeometry.py
+++ b/tests/src/python/test_qgsgeometry.py
@@ -5151,6 +5151,35 @@ def testForceRHR(self):
             self.assertEqual(res.asWkt(1), t[1],
                              "mismatch for {}, expected:\n{}\nGot:\n{}\n".format(t[0], t[1], res.asWkt(1)))
 
+    def testLineStringFromBezier(self):
+        tests = [
+            [QgsPoint(1, 1), QgsPoint(10, 1), QgsPoint(10, 10), QgsPoint(20, 10), 5, 'LineString (1 1, 5.5 1.9, 8.7 4.2, 11.6 6.8, 15 9.1, 20 10)'],
+            [QgsPoint(1, 1), QgsPoint(10, 1), QgsPoint(10, 10), QgsPoint(1, 1), 10,
+             'LineString (1 1, 3.4 1.2, 5.3 1.9, 6.7 2.7, 7.5 3.6, 7.8 4.4, 7.5 4.9, 6.7 5, 5.3 4.5, 3.4 3.2, 1 1)'],
+            [QgsPoint(1, 1), QgsPoint(10, 1), QgsPoint(10, 10), QgsPoint(20, 10), 10,
+             'LineString (1 1, 3.4 1.3, 5.5 1.9, 7.2 2.9, 8.7 4.2, 10.1 5.5, 11.6 6.8, 13.2 8.1, 15 9.1, 17.3 9.7, 20 10)'],
+            [QgsPoint(1, 1), QgsPoint(10, 1), QgsPoint(10, 10), QgsPoint(20, 10), 1,
+             'LineString (1 1, 20 10)'],
+            [QgsPoint(1, 1), QgsPoint(10, 1), QgsPoint(10, 10), QgsPoint(20, 10), 0,
+             'LineString EMPTY'],
+            [QgsPoint(1, 1, 2), QgsPoint(10, 1), QgsPoint(10, 10), QgsPoint(20, 10), 5,
+             'LineString (1 1, 5.5 1.9, 8.7 4.2, 11.6 6.8, 15 9.1, 20 10)'],
+            [QgsPoint(1, 1), QgsPoint(10, 1, 2), QgsPoint(10, 10), QgsPoint(20, 10), 5,
+             'LineString (1 1, 5.5 1.9, 8.7 4.2, 11.6 6.8, 15 9.1, 20 10)'],
+            [QgsPoint(1, 1, 2), QgsPoint(10, 1), QgsPoint(10, 10, 2), QgsPoint(20, 10), 5,
+             'LineString (1 1, 5.5 1.9, 8.7 4.2, 11.6 6.8, 15 9.1, 20 10)'],
+            [QgsPoint(1, 1, 2), QgsPoint(10, 1), QgsPoint(10, 10), QgsPoint(20, 10, 2), 5,
+             'LineString (1 1, 5.5 1.9, 8.7 4.2, 11.6 6.8, 15 9.1, 20 10)'],
+            [QgsPoint(1, 1, 1), QgsPoint(10, 1, 2), QgsPoint(10, 10, 3), QgsPoint(20, 10, 4), 5,
+             'LineStringZ (1 1 1, 5.5 1.9 1.6, 8.7 4.2 2.2, 11.6 6.8 2.8, 15 9.1 3.4, 20 10 4)'],
+            [QgsPoint(1, 1, 1, 10), QgsPoint(10, 1, 2, 9), QgsPoint(10, 10, 3, 2), QgsPoint(20, 10, 4, 1), 5,
+             'LineStringZM (1 1 1 10, 5.5 1.9 1.6 8.8, 8.7 4.2 2.2 6.7, 11.6 6.8 2.8 4.3, 15 9.1 3.4 2.2, 20 10 4 1)']
+        ]
+        for t in tests:
+            res = QgsLineString.fromBezierCurve(t[0], t[1], t[2], t[3], t[4])
+            self.assertEqual(res.asWkt(1), t[5],
+                             "mismatch for {}, expected:\n{}\nGot:\n{}\n".format(t[0], t[5], res.asWkt(1)))
+
     def testIsGeosValid(self):
         tests = [
             ["", False, False, ''],