diff --git a/src/lib/Math/Bezier/Cubic.hs b/src/lib/Math/Bezier/Cubic.hs
index 4307b80..930fe21 100644
--- a/src/lib/Math/Bezier/Cubic.hs
+++ b/src/lib/Math/Bezier/Cubic.hs
@@ -21,7 +21,7 @@ module Math.Bezier.Cubic
   , curvature, squaredCurvature
   , subdivide
   , ddist, closestPoint
-  , drag
+  , drag, selfIntersectionParameters
   )
   where
 
@@ -73,7 +73,9 @@ import Math.Module
   , Inner(..), squaredNorm
   )
 import Math.Roots
-  ( realRoots )
+  ( realRoots, solveQuadratic )
+import Math.Vector2D
+  ( Point2D(..) )
 
 --------------------------------------------------------------------------------
 
@@ -111,11 +113,9 @@ bezier ( Bezier {..} ) t =
 
 -- | Derivative of cubic Bézier curve.
 bezier' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
-bezier' ( Bezier {..} ) t
+bezier' ( Bezier {..} )
   = ( 3 *^ )
-  $ lerp @v t
-      ( lerp @v t ( p0 --> p1 ) ( p1 --> p2 ) )
-      ( lerp @v t ( p1 --> p2 ) ( p2 --> p3 ) )
+  . Quadratic.bezier @v ( Quadratic.Bezier ( p0 --> p1 ) ( p1 --> p2 ) ( p2 --> p3 ) )
 
 -- | Second derivative of a cubic Bézier curve.
 bezier'' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
@@ -179,7 +179,7 @@ ddist ( Bezier {..} ) c = [ a5, a4, a3, a2, a1, a0 ]
 closestPoint
   :: forall v r p. ( Torsor v p, Inner r v, RealFloat r, Prim r, NFData r )
   => Bezier p -> p -> ArgMin r ( r, p )
-closestPoint pts c = pickClosest ( 0 :| 1 : roots ) -- todo: also include the self-intersection point if one exists
+closestPoint pts c = pickClosest ( 0 :| 1 : roots )
   where
     roots :: [ r ]
     roots = filter ( \ r -> r > 0 && r < 1 ) ( realRoots 50 $ ddist @v pts c )
@@ -222,3 +222,29 @@ drag ( Bezier {..} ) t q = Bezier { p0, p1 = p1', p2 = p2', p3 }
     p1', p2' :: p
     p1' = ( ( 1 - t ) *^ delta ) • p1
     p2' = (     t     *^ delta ) • p2
+
+-- | Compute parameter values for the self-intersection of a planar cubic Bézier curve, if such exist.
+--
+-- The parameter values might lie outside the interval [0,1], 
+-- indicating a self-intersection of the extended curve.
+--
+-- Formula taken from:
+--  "A Basis for the Implicit Representation of Planar Rational Cubic Bézier Curves"
+--      – Oliver J. D. Barrowclough, 2016
+selfIntersectionParameters :: forall r. RealFloat r => Bezier ( Point2D r ) -> [ r ]
+selfIntersectionParameters ( Bezier {..} ) = solveQuadratic c0 c1 c2
+  where
+    areaConstant :: Point2D r -> Point2D r -> Point2D r -> r
+    areaConstant ( Point2D x1 y1 ) ( Point2D x2 y2 ) ( Point2D x3 y3 ) =
+      x1 * ( y2 - y3 ) + x2 * ( y3 - y1 ) + x3 * ( y1 - y2 )
+    l0, l1, l2, l3, f1, f2, f3, c0, c1, c2 :: r
+    l0 = areaConstant p3 p2 p1
+    l1 = areaConstant p2 p3 p0
+    l2 = areaConstant p1 p0 p3
+    l3 = areaConstant p0 p1 p2
+    f1 = 3 * ( l1 * l1 - 3 * l0 * l2 )
+    f2 = 3 * ( l2 * l2 - 3 * l1 * l3 )
+    f3 = 3 * ( 9 * l0 * l3 - l1 * l2 )
+    c0 = f2
+    c1 = f3 - 2 * f2
+    c2 = f1 + f2 - f3