diff --git a/MetaBrush.cabal b/MetaBrush.cabal
index 5797aec..626c74c 100644
--- a/MetaBrush.cabal
+++ b/MetaBrush.cabal
@@ -83,6 +83,7 @@ library
     exposed-modules:
         Math.Bezier.Cubic
       , Math.Bezier.Cubic.Fit
+      , Math.Bezier.Envelope
       , Math.Bezier.Quadratic
       , Math.Bezier.Spline
       , Math.Bezier.Stroke
diff --git a/src/lib/Math/Bezier/Envelope.hs b/src/lib/Math/Bezier/Envelope.hs
new file mode 100644
index 0000000..823f247
--- /dev/null
+++ b/src/lib/Math/Bezier/Envelope.hs
@@ -0,0 +1,455 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications    #-}
+
+module Math.Bezier.Envelope where
+
+-- acts
+import Data.Act
+  ( Torsor((-->)) )
+
+-- deepseq
+import Control.DeepSeq
+  ( NFData )
+
+-- primitive
+import Data.Primitive.Types
+  ( Prim )
+
+-- MetaBrush
+import Math.Roots
+  ( realRoots )
+import qualified Math.Bezier.Cubic as Cubic
+  ( Bezier(..), bezier, bezier' )
+import qualified Math.Bezier.Quadratic as Quadratic
+  ( Bezier(..), bezier, bezier' )
+import Math.Module
+  ( Module((^+^),(*^)), lerp, cross )
+import Math.Vector2D
+  ( Point2D(..), Vector2D(..), Segment(..) )
+
+--------------------------------------------------------------------------------
+
+-- | Find the roots of the envelope equation for a family of cubic Bézier curves
+-- varying along a cubic Bézier path.
+--
+-- \[ c(t,u) = p(t) + b(t,u), \]
+--
+-- where \( t \mapsto p(t) \) describes the underlying path,
+-- and \( u \mapsto b(t_0,u) \) describes the brush shape at point \( t = t_0 \).
+--
+-- The envelope equation is then:
+--
+-- \[ \frac{\partial c}{\partial t} \cross \frac{\partial c}{\partial u} = 0. \]
+--
+-- Given \( t_0 \), this function returns a (possibly empty) list of values \( u_i \)
+-- satisfying the envelope equation at \( (t_0, u_i) \).
+--
+-- The points \( c(t_0,u_i) \) are thus potential outline points on the contour stroked
+-- by the brush as it moves along the path.
+envelope33
+  :: forall r
+  .  ( RealFloat r, Prim r, NFData r )
+  => Cubic.Bezier ( Point2D r ) -> Cubic.Bezier ( Cubic.Bezier ( Point2D r ) ) -> r -> [ r ]
+envelope33 path
+  ( Cubic.Bezier
+    ( Cubic.Bezier b00 b01 b02 b03 )
+    ( Cubic.Bezier b10 b11 b12 b13 )
+    ( Cubic.Bezier b20 b21 b22 b23 )
+    ( Cubic.Bezier b30 b31 b32 b33 )
+  ) t0 = realRoots 50 [ a5, a4, a3, a2, a1, a0 ]
+
+  where
+
+    -- Compute ∂p/∂t(t0).
+    dpdt :: Vector2D r
+    dpdt = Cubic.bezier' @( Vector2D r ) path t0
+
+    -- Compute ∂b/∂t(t0,u) using the Bernstein basis:
+    --
+    -- ∂b/∂t(t0,u) = Cubic.bezier ( Cubic.Bezier ct0 ct1 ct2 ct3 ) u.
+    ct0, ct1, ct2, ct3, dt0, dt1, dt2, dt3 :: Vector2D r
+    ct0 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b00 b10 b20 b30 ) t0
+    ct1 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b01 b11 b21 b31 ) t0
+    ct2 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b02 b12 b22 b32 ) t0
+    ct3 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b03 b13 b23 b33 ) t0
+
+    -- Add ∂p/∂t and convert the Bernstein representation to the monomial basis to obtain
+    --
+    -- ∂c/∂t(t0,u) = dt0 + u dt1 + u² dt2 + u³ dt3.
+    dt0 = ct0 ^+^ dpdt
+    dt1 = 3 *^ ( ct0 --> ct1 )
+    dt2 = 3 *^ ( ct1 --> ct0 ^+^ ct1 --> ct2 )
+    dt3 = ct0 --> ct3 ^+^ 3 *^ ( ct2 --> ct1 )
+
+    -- Compute ∂c/∂u(t0,u) using the Bernstein basis:
+    --
+    -- ∂c/∂u(t0,u) = Cubic.bezier' ( Cubic.Bezier cu0 cu1 cu2 cu3 ) u.
+    cu0, cu1, cu2, cu3 :: Point2D r
+    cu0 = Cubic.bezier  @( Vector2D r ) ( Cubic.Bezier b00 b10 b20 b30 ) t0
+    cu1 = Cubic.bezier  @( Vector2D r ) ( Cubic.Bezier b01 b11 b21 b31 ) t0
+    cu2 = Cubic.bezier  @( Vector2D r ) ( Cubic.Bezier b02 b12 b22 b32 ) t0
+    cu3 = Cubic.bezier  @( Vector2D r ) ( Cubic.Bezier b03 b13 b23 b33 ) t0
+
+    -- Convert the Bernstein representation to the monomial basis to obtain
+    --
+    -- ∂c/∂u(t0,u) = du0 + u du1 + u² du2.
+    du0, du1, du2 :: Vector2D r
+    du0 = cu0 --> cu1
+    du1 = 2 *^ ( cu1 --> cu0 ^+^ cu1 --> cu2 )
+    du2 = cu0 --> cu3 ^+^ 3 *^ ( cu2 --> cu1 )
+
+    -- Expand out the cross-product ∂c/∂t × ∂c/∂u to obtain the envelope equation:
+    --
+    -- a0 + a1 u + a2 u² + a3 u³ + a4 u⁴ + a5 u⁵.
+    a0, a1, a2, a3, a4, a5 :: r
+    a0 = dt0 `cross` du0
+    a1 = dt1 `cross` du0 + dt0 `cross` du1
+    a2 = dt2 `cross` du0 + dt1 `cross` du1 + dt0 `cross` du2
+    a3 = dt3 `cross` du0 + dt2 `cross` du1 + dt1 `cross` du2
+    a4 = dt3 `cross` du1 + dt2 `cross` du2
+    a5 = dt3 `cross` du2
+
+-- | Find the roots of the envelope equation for a family of cubic Bézier curves
+-- varying along a quadratic Bézier path.
+--
+-- See 'envelope33' for more information.
+envelope23
+  :: forall r
+  .  ( RealFloat r, Prim r, NFData r )
+  => Quadratic.Bezier ( Point2D r ) -> Quadratic.Bezier ( Cubic.Bezier ( Point2D r ) ) -> r -> [ r ]
+envelope23 path
+  ( Quadratic.Bezier
+    ( Cubic.Bezier b00 b01 b02 b03 )
+    ( Cubic.Bezier b10 b11 b12 b13 )
+    ( Cubic.Bezier b20 b21 b22 b23 )
+  ) t0 = realRoots 50 [ a5, a4, a3, a2, a1, a0 ]
+
+  where
+
+    dpdt :: Vector2D r
+    dpdt = Quadratic.bezier' @( Vector2D r ) path t0
+
+    ct0, ct1, ct2, ct3, dt0, dt1, dt2, dt3 :: Vector2D r
+    ct0 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b00 b10 b20 ) t0
+    ct1 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b01 b11 b21 ) t0
+    ct2 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b02 b12 b22 ) t0
+    ct3 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b03 b13 b23 ) t0
+
+    dt0 = ct0 ^+^ dpdt
+    dt1 = 3 *^ ( ct0 --> ct1 )
+    dt2 = 3 *^ ( ct1 --> ct0 ^+^ ct1 --> ct2 )
+    dt3 = ct0 --> ct3 ^+^ 3 *^ ( ct2 --> ct1 )
+
+    cu0, cu1, cu2, cu3 :: Point2D r
+    cu0 = Quadratic.bezier  @( Vector2D r ) ( Quadratic.Bezier b00 b10 b20 ) t0
+    cu1 = Quadratic.bezier  @( Vector2D r ) ( Quadratic.Bezier b01 b11 b21 ) t0
+    cu2 = Quadratic.bezier  @( Vector2D r ) ( Quadratic.Bezier b02 b12 b22 ) t0
+    cu3 = Quadratic.bezier  @( Vector2D r ) ( Quadratic.Bezier b03 b13 b23 ) t0
+
+    du0, du1, du2 :: Vector2D r
+    du0 = cu0 --> cu1
+    du1 = 2 *^ ( cu1 --> cu0 ^+^ cu1 --> cu2 )
+    du2 = cu0 --> cu3 ^+^ 3 *^ ( cu2 --> cu1 )
+
+    a0, a1, a2, a3, a4, a5 :: r
+    a0 = dt0 `cross` du0
+    a1 = dt1 `cross` du0 + dt0 `cross` du1
+    a2 = dt2 `cross` du0 + dt1 `cross` du1 + dt0 `cross` du2
+    a3 = dt3 `cross` du0 + dt2 `cross` du1 + dt1 `cross` du2
+    a4 = dt3 `cross` du1 + dt2 `cross` du2
+    a5 = dt3 `cross` du2
+
+-- | Find the roots of the envelope equation for a family of cubic Bézier curves
+-- varying along a straight line path.
+--
+-- See 'envelope33' for more information.
+envelope13
+  :: forall r
+  .  ( RealFloat r, Prim r, NFData r )
+  => Segment ( Point2D r ) -> Segment ( Cubic.Bezier ( Point2D r ) ) -> r -> [ r ]
+envelope13 ( Segment p0 p1 )
+  ( Segment
+    ( Cubic.Bezier b00 b01 b02 b03 )
+    ( Cubic.Bezier b10 b11 b12 b13 )
+  ) t0 = realRoots 50 [ a5, a4, a3, a2, a1, a0 ]
+
+  where
+
+    dpdt :: Vector2D r
+    dpdt = p0 --> p1
+
+    ct0, ct1, ct2, ct3, dt0, dt1, dt2, dt3 :: Vector2D r
+    ct0 = b00 --> b10
+    ct1 = b01 --> b11
+    ct2 = b02 --> b12
+    ct3 = b03 --> b13
+
+    dt0 = ct0 ^+^ dpdt
+    dt1 = 3 *^ ( ct0 --> ct1 )
+    dt2 = 3 *^ ( ct1 --> ct0 ^+^ ct1 --> ct2 )
+    dt3 = ct0 --> ct3 ^+^ 3 *^ ( ct2 --> ct1 )
+
+    cu0, cu1, cu2, cu3 :: Point2D r
+    cu0 = lerp @( Vector2D r ) t0 b00 b10
+    cu1 = lerp @( Vector2D r ) t0 b01 b11
+    cu2 = lerp @( Vector2D r ) t0 b02 b12
+    cu3 = lerp @( Vector2D r ) t0 b03 b13
+
+    du0, du1, du2 :: Vector2D r
+    du0 = cu0 --> cu1
+    du1 = 2 *^ ( cu1 --> cu0 ^+^ cu1 --> cu2 )
+    du2 = cu0 --> cu3 ^+^ 3 *^ ( cu2 --> cu1 )
+
+    a0, a1, a2, a3, a4, a5 :: r
+    a0 = dt0 `cross` du0
+    a1 = dt1 `cross` du0 + dt0 `cross` du1
+    a2 = dt2 `cross` du0 + dt1 `cross` du1 + dt0 `cross` du2
+    a3 = dt3 `cross` du0 + dt2 `cross` du1 + dt1 `cross` du2
+    a4 = dt3 `cross` du1 + dt2 `cross` du2
+    a5 = dt3 `cross` du2
+
+-- | Find the roots of the envelope equation for a family of quadratic Bézier curves
+-- varying along a cubic Bézier path.
+--
+-- See 'envelope33' for more information.
+envelope32
+  :: forall r
+  .  ( RealFloat r, Prim r, NFData r )
+  => Cubic.Bezier ( Point2D r ) -> Cubic.Bezier ( Quadratic.Bezier ( Point2D r ) ) -> r -> [ r ]
+envelope32 path
+  ( Cubic.Bezier
+    ( Quadratic.Bezier b00 b01 b02 )
+    ( Quadratic.Bezier b10 b11 b12 )
+    ( Quadratic.Bezier b20 b21 b22 )
+    ( Quadratic.Bezier b30 b31 b32 )
+  ) t0 = realRoots 50 [ a3, a2, a1, a0 ]
+
+  where
+
+    dpdt :: Vector2D r
+    dpdt = Cubic.bezier' @( Vector2D r ) path t0
+
+    ct0, ct1, ct2, dt0, dt1, dt2 :: Vector2D r
+    ct0 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b00 b10 b20 b30 ) t0
+    ct1 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b01 b11 b21 b31 ) t0
+    ct2 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b02 b12 b22 b32 ) t0
+
+    dt0 = ct0 ^+^ dpdt
+    dt1 = 2 *^ ( ct0 --> ct1 )
+    dt2 = ct1 --> ct0 ^+^ ct1 --> ct2
+
+    cu0, cu1, cu2 :: Point2D r
+    cu0 = Cubic.bezier  @( Vector2D r ) ( Cubic.Bezier b00 b10 b20 b30 ) t0
+    cu1 = Cubic.bezier  @( Vector2D r ) ( Cubic.Bezier b01 b11 b21 b31 ) t0
+    cu2 = Cubic.bezier  @( Vector2D r ) ( Cubic.Bezier b02 b12 b22 b32 ) t0
+
+    du0, du1 :: Vector2D r
+    du0 = cu0 --> cu1
+    du1 = cu1 --> cu0 ^+^ cu1 --> cu2
+
+    a0, a1, a2, a3 :: r
+    a0 = dt0 `cross` du0
+    a1 = dt1 `cross` du0 + dt0 `cross` du1
+    a2 = dt2 `cross` du0 + dt1 `cross` du1
+    a3 = dt2 `cross` du1
+
+-- | Find the roots of the envelope equation for a family of quadratic Bézier curves
+-- varying along a quadratic Bézier path.
+--
+-- See 'envelope33' for more information.
+envelope22
+  :: forall r
+  .  ( RealFloat r, Prim r, NFData r )
+  => Quadratic.Bezier ( Point2D r ) -> Quadratic.Bezier ( Quadratic.Bezier ( Point2D r ) ) -> r -> [ r ]
+envelope22 path
+  ( Quadratic.Bezier
+    ( Quadratic.Bezier b00 b01 b02 )
+    ( Quadratic.Bezier b10 b11 b12 )
+    ( Quadratic.Bezier b20 b21 b22 )
+  ) t0 = realRoots 50 [ a3, a2, a1, a0 ]
+
+  where
+
+    dpdt :: Vector2D r
+    dpdt = Quadratic.bezier' @( Vector2D r ) path t0
+
+    ct0, ct1, ct2, dt0, dt1, dt2 :: Vector2D r
+    ct0 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b00 b10 b20 ) t0
+    ct1 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b01 b11 b21 ) t0
+    ct2 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b02 b12 b22 ) t0
+
+    dt0 = ct0 ^+^ dpdt
+    dt1 = 2 *^ ( ct0 --> ct1 )
+    dt2 = ct1 --> ct0 ^+^ ct1 --> ct2
+
+    cu0, cu1, cu2 :: Point2D r
+    cu0 = Quadratic.bezier  @( Vector2D r ) ( Quadratic.Bezier b00 b10 b20 ) t0
+    cu1 = Quadratic.bezier  @( Vector2D r ) ( Quadratic.Bezier b01 b11 b21 ) t0
+    cu2 = Quadratic.bezier  @( Vector2D r ) ( Quadratic.Bezier b02 b12 b22 ) t0
+
+    du0, du1 :: Vector2D r
+    du0 = cu0 --> cu1
+    du1 = cu1 --> cu0 ^+^ cu1 --> cu2
+
+    a0, a1, a2, a3 :: r
+    a0 = dt0 `cross` du0
+    a1 = dt1 `cross` du0 + dt0 `cross` du1
+    a2 = dt2 `cross` du0 + dt1 `cross` du1
+    a3 = dt2 `cross` du1
+
+-- | Find the roots of the envelope equation for a family of quadratic Bézier curves
+-- varying along a straight line.
+--
+-- See 'envelope33' for more information.
+envelope12
+  :: forall r
+  .  ( RealFloat r, Prim r, NFData r )
+  => Segment ( Point2D r ) -> Segment ( Quadratic.Bezier ( Point2D r ) ) -> r -> [ r ]
+envelope12 ( Segment p0 p1 )
+  ( Segment
+    ( Quadratic.Bezier b00 b01 b02 )
+    ( Quadratic.Bezier b10 b11 b12 )
+  ) t0 = realRoots 50 [ a3, a2, a1, a0 ]
+
+  where
+
+    dpdt :: Vector2D r
+    dpdt = p0 --> p1
+
+    ct0, ct1, ct2, dt0, dt1, dt2 :: Vector2D r
+    ct0 = b00 --> b10
+    ct1 = b01 --> b11
+    ct2 = b02 --> b12
+
+    dt0 = ct0 ^+^ dpdt
+    dt1 = 2 *^ ( ct0 --> ct1 )
+    dt2 = ct1 --> ct0 ^+^ ct1 --> ct2
+
+    cu0, cu1, cu2 :: Point2D r
+    cu0 = lerp @( Vector2D r ) t0 b00 b10
+    cu1 = lerp @( Vector2D r ) t0 b01 b11
+    cu2 = lerp @( Vector2D r ) t0 b02 b12
+
+    du0, du1 :: Vector2D r
+    du0 = cu0 --> cu1
+    du1 = cu1 --> cu0 ^+^ cu1 --> cu2
+
+    a0, a1, a2, a3 :: r
+    a0 = dt0 `cross` du0
+    a1 = dt1 `cross` du0 + dt0 `cross` du1
+    a2 = dt2 `cross` du0 + dt1 `cross` du1
+    a3 = dt2 `cross` du1
+
+-- | Find the roots of the envelope equation for a family of line segments
+-- varying along a cubic Bézier curve.
+--
+-- See 'envelope33' for more information.
+envelope31
+  :: forall r
+  .  ( RealFloat r, Prim r, NFData r )
+  => Cubic.Bezier ( Point2D r ) -> Cubic.Bezier ( Segment ( Point2D r ) ) -> r -> [ r ]
+envelope31 path
+  ( Cubic.Bezier
+    ( Segment b00 b01 )
+    ( Segment b10 b11 )
+    ( Segment b20 b21 )
+    ( Segment b30 b31 )
+  ) t0 = [ - a1 / a0 ]
+
+  where
+
+    dpdt :: Vector2D r
+    dpdt = Cubic.bezier' @( Vector2D r ) path t0
+
+    ct0, ct1, dt0, dt1 :: Vector2D r
+    ct0 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b00 b10 b20 b30 ) t0
+    ct1 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b01 b11 b21 b31 ) t0
+
+    dt0 = ct0 ^+^ dpdt
+    dt1 = ct0 --> ct1
+
+    cu0, cu1 :: Point2D r
+    cu0 = Cubic.bezier  @( Vector2D r ) ( Cubic.Bezier b00 b10 b20 b30 ) t0
+    cu1 = Cubic.bezier  @( Vector2D r ) ( Cubic.Bezier b01 b11 b21 b31 ) t0
+
+    du0 :: Vector2D r
+    du0 = cu0 --> cu1
+
+    a0, a1 :: r
+    a0 = dt0 `cross` du0
+    a1 = dt1 `cross` du0
+
+-- | Find the roots of the envelope equation for a family of line segments
+-- varying along a quadratic Bézier curve.
+--
+-- See 'envelope33' for more information.
+envelope21
+  :: forall r
+  .  ( RealFloat r, Prim r, NFData r )
+  => Quadratic.Bezier ( Point2D r ) -> Quadratic.Bezier ( Segment ( Point2D r ) ) -> r -> [ r ]
+envelope21 path
+  ( Quadratic.Bezier
+    ( Segment b00 b01 )
+    ( Segment b10 b11 )
+    ( Segment b20 b21 )
+  ) t0 = [ - a1 / a0 ]
+
+  where
+
+    dpdt :: Vector2D r
+    dpdt = Quadratic.bezier' @( Vector2D r ) path t0
+
+    ct0, ct1, dt0, dt1 :: Vector2D r
+    ct0 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b00 b10 b20 ) t0
+    ct1 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b01 b11 b21 ) t0
+
+    dt0 = ct0 ^+^ dpdt
+    dt1 = ct0 --> ct1
+
+    cu0, cu1 :: Point2D r
+    cu0 = Quadratic.bezier  @( Vector2D r ) ( Quadratic.Bezier b00 b10 b20 ) t0
+    cu1 = Quadratic.bezier  @( Vector2D r ) ( Quadratic.Bezier b01 b11 b21 ) t0
+
+    du0 :: Vector2D r
+    du0 = cu0 --> cu1
+
+    a0, a1 :: r
+    a0 = dt0 `cross` du0
+    a1 = dt1 `cross` du0
+
+-- | Find the roots of the envelope equation for a family of line segments
+-- varying along a straight line path.
+--
+-- See 'envelope33' for more information.
+envelope11
+  :: forall r
+  .  ( RealFloat r, Prim r, NFData r )
+  => Segment ( Point2D r ) -> Segment ( Segment ( Point2D r ) ) -> r -> [ r ]
+envelope11 ( Segment p0 p1 )
+  ( Segment
+    ( Segment b00 b01 )
+    ( Segment b10 b11 )
+  ) t0 = [ - a1 / a0 ]
+
+  where
+
+    dpdt :: Vector2D r
+    dpdt = p0 --> p1
+
+    ct0, ct1, dt0, dt1 :: Vector2D r
+    ct0 = b00 --> b10
+    ct1 = b01 --> b11
+
+    dt0 = ct0 ^+^ dpdt
+    dt1 = ct0 --> ct1
+
+    cu0, cu1 :: Point2D r
+    cu0 = lerp @( Vector2D r ) t0 b00 b10
+    cu1 = lerp @( Vector2D r ) t0 b01 b11
+
+    du0 :: Vector2D r
+    du0 = cu0 --> cu1
+
+    a0, a1 :: r
+    a0 = dt0 `cross` du0
+    a1 = dt1 `cross` du0