diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs
index 82dfd88..6534688 100644
--- a/brush-strokes/src/cusps/bench/Main.hs
+++ b/brush-strokes/src/cusps/bench/Main.hs
@@ -133,8 +133,8 @@ benchTestCase ( TestCase { testName, testBrushStroke, testCuspOptions, testStart
   after <- getMonotonicTime
   let dt = after - before
   putStrLn $ unlines
-    [ "  -  #sols: " ++ show (length sols)
-    , "  - #dunno: " ++ show (length dunno)
+    [ "  -  #sols: " ++ show sols --( length sols )
+    , "  - #dunno: " ++ show ( length dunno )
     , "  - Time elapsed: " ++ show dt ++ "s"
     , "" ]
   return dt
@@ -181,7 +181,7 @@ testCases = benchCases
 
 benchCases :: [ TestCase ]
 benchCases =
-  [ ellipseTestCase opts "full" ( 0, 1 ) pi $ defaultStartBoxes [ 0 .. 3 ] ]
+  [ ellipseTestCase opts "full" ( 0, 1 ) pi $ defaultStartBoxes [ 2 ] ] -- [ 0 .. 3 ] ]
   where
     opts = defaultRootIsolationOptions
 
@@ -190,7 +190,7 @@ benchCases =
 data BrushStroke =
   forall nbParams. ParamsCt nbParams =>
   BrushStroke
-  { brush  :: !( Brush nbParams )
+  { brush  :: !( Brush ( ℝ nbParams ) )
   , stroke :: !( Point nbParams, Curve Open () ( Point nbParams ) )
   }
 
@@ -500,7 +500,7 @@ ellipseTestCase opts str k0k1 rot startBoxes =
 ellipseBrushStroke :: ( Double, Double ) -> Double -> BrushStroke
 ellipseBrushStroke ( k0, k1 ) rot =
   BrushStroke
-    { brush = ellipseBrush
+    { brush  = ellipseBrush
     , stroke = ( p0, LineTo ( NextPoint p1 ) () ) }
   where
     mkPt x y w h phi =
@@ -568,18 +568,10 @@ data Point nbParams =
   deriving stock Generic
 deriving stock instance Show ( ℝ nbParams ) => Show ( Point nbParams )
 
-type Brush nbParams
-  = forall {t} k (i :: t)
-  .  DiffInterp k i ( ℝ nbParams )
-  => Proxy# i
-  -> ( forall a. a -> I i a )
-  -> C k ( I i ( ℝ nbParams ) )
-         ( Spline Closed () ( I i ( ℝ 2 ) ) )
-
 getStrokeFunctions
   :: forall nbParams
   .  ParamsCt nbParams
-  => Brush nbParams
+  => Brush ( ℝ nbParams )
       -- ^ brush shape
   -> Point nbParams
       -- ^ start point
@@ -587,7 +579,7 @@ getStrokeFunctions
       -- ^ curve points
   -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () )
      , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-getStrokeFunctions brush sp0 crv =
+getStrokeFunctions ( Brush brushShape brushShapeI mbRot ) sp0 crv =
   let
     usedParams :: C 2 ( ℝ 1 ) ( ℝ nbParams )
     path :: C 2 ( ℝ 1 ) ( ℝ 2 )
@@ -600,11 +592,14 @@ getStrokeFunctions brush sp0 crv =
       pathAndUsedParams @3 @𝕀 coerce singleton ( getParams . pointParams )
         sp0 crv
   in ( brushStrokeData @2 @( ℝ nbParams ) coerce coerce
-        path usedParams $
-        brush @2 @() proxy# id
+        path usedParams
+        brushShape
+        mbRot
      , brushStrokeData @3 @( ℝ nbParams ) coerce coerce
-        pathI usedParamsI $
-        brush @3 @𝕀 proxy# singleton )
+        pathI usedParamsI
+        brushShapeI
+        ( fmap nonDecreasing mbRot )
+     )
 {-# INLINEABLE getStrokeFunctions #-}
 
 defaultStartBoxes :: [ Int ] -> [ ( Int, [ Box 2 ] ) ]
diff --git a/brush-strokes/src/cusps/inspect/Main.hs b/brush-strokes/src/cusps/inspect/Main.hs
index a412124..485330f 100644
--- a/brush-strokes/src/cusps/inspect/Main.hs
+++ b/brush-strokes/src/cusps/inspect/Main.hs
@@ -1,5 +1,4 @@
 {-# LANGUAGE AllowAmbiguousTypes   #-}
-{-# LANGUAGE PartialTypeSignatures #-}
 {-# LANGUAGE PolyKinds             #-}
 {-# LANGUAGE RebindableSyntax      #-}
 {-# LANGUAGE ScopedTypeVariables   #-}
diff --git a/brush-strokes/src/lib/Calligraphy/Brushes.hs b/brush-strokes/src/lib/Calligraphy/Brushes.hs
index 774c2b1..4880b8c 100644
--- a/brush-strokes/src/lib/Calligraphy/Brushes.hs
+++ b/brush-strokes/src/lib/Calligraphy/Brushes.hs
@@ -6,7 +6,8 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 
 module Calligraphy.Brushes
-  ( circleBrush
+  ( BrushFn, Brush(..)
+  , circleBrush
   , ellipseBrush
   , tearDropBrush
   ) where
@@ -15,10 +16,12 @@ module Calligraphy.Brushes
 import Prelude
   hiding
     ( Num(..), Floating(..), (^), (/), fromInteger, fromRational )
+import Data.Kind
+  ( Type )
 import GHC.Exts
-  ( Proxy# )
+  ( Proxy#, proxy# )
 import GHC.TypeNats
-  ( type (<=) )
+  ( Nat, type (<=) )
 
 -- containers
 import qualified Data.Sequence as Seq
@@ -29,11 +32,71 @@ import Math.Algebra.Dual
 import Math.Bezier.Spline
 import Math.Differentiable
   ( I )
+import Math.Interval
+  ( 𝕀, singleton )
 import Math.Linear
 import Math.Module
   ( Module((^+^), (*^)) )
 import Math.Ring
 
+--------------------------------------------------------------------------------
+
+-- | The shape of a brush (before applying any rotation).
+type BrushFn :: forall {kd}. kd -> Nat -> Type -> Type
+type BrushFn i k brushParams
+  = C k ( I i brushParams )
+        ( Spline Closed () ( I i ( ℝ 2 ) ) )
+
+-- | A brush, described as a base shape + an optional rotation.
+data Brush brushParams
+  = Brush
+  { brushShape  :: BrushFn () 2 brushParams
+  , brushShapeI :: BrushFn 𝕀 3 brushParams
+  , mbRotation  :: Maybe ( brushParams -> Double )
+  }
+
+--------------------------------------------------------------------------------
+-- Brushes
+
+-- Some convenience type synonyms for brush types... a bit horrible
+type ParamsCt rec = ( ParamsICt 2 () rec, ParamsICt 3 𝕀 rec )
+type ParamsICt k i rec =
+    ( Module
+       ( D k ( I i rec ) ( I i Double ) )
+       ( D k ( I i rec ) ( I i ( ℝ 2 ) ) )
+    , Module ( I i Double ) ( T ( I i Double ) )
+    , HasChainRule ( I i Double ) k ( I i rec )
+    , Representable ( I i Double ) ( I i rec )
+    , Applicative ( D k ( I i rec ) )
+    )
+
+{-# INLINEABLE circleBrush #-}
+circleBrush :: ( 1 <= RepDim params, ParamsCt params ) => Brush params
+circleBrush =
+  Brush
+    { brushShape  = circleBrushFn @() @2 proxy# id
+    , brushShapeI = circleBrushFn @𝕀 @3 proxy# singleton
+    , mbRotation  = Nothing
+    }
+
+{-# INLINEABLE ellipseBrush #-}
+ellipseBrush :: ( 3 <= RepDim params, ParamsCt params ) => Brush params
+ellipseBrush =
+  Brush
+    { brushShape  = ellipseBrushFn @() @2 proxy# id
+    , brushShapeI = ellipseBrushFn @𝕀 @3 proxy# singleton
+    , mbRotation  = Just ( `index` ( Fin 3 ) )
+    }
+
+{-# INLINEABLE tearDropBrush #-}
+tearDropBrush :: ( 3 <= RepDim params, ParamsCt params ) => Brush params
+tearDropBrush =
+  Brush
+    { brushShape  = tearDropBrushFn @() @2 proxy# id
+    , brushShapeI = tearDropBrushFn @𝕀 @3 proxy# singleton
+    , mbRotation  = Just ( `index` ( Fin 3 ) )
+    }
+
 --------------------------------------------------------------------------------
 -- Circle & ellipse
 
@@ -60,93 +123,55 @@ circleSpline p = sequenceA $
         Bezier3To ( p  κ -1 ) ( p  1 -κ ) BackToStart ()
 {-# INLINE circleSpline #-}
 
-circleBrush :: forall {t} (i :: t) k irec
-            . ( 1 <= RepDim irec
-              , Module
-                 ( D k irec ( I i Double ) )
-                 ( D k irec ( I i ( ℝ 2 ) ) )
-              , Module ( I i Double ) ( T ( I i Double ) )
-              , HasChainRule ( I i Double ) k irec
-              , Representable ( I i Double ) irec
-              , Applicative ( D k irec )
-              )
-            => Proxy# i
-            -> ( forall a. a -> I i a )
-            -> C k irec ( Spline 'Closed () ( I i ( ℝ 2 ) ) )
-circleBrush _ mkI =
+circleBrushFn :: forall {t} (i :: t) k rec
+              . ( 1 <= RepDim ( I i rec )
+                , ParamsICt k i rec
+                )
+              => Proxy# i
+              -> ( forall a. a -> I i a )
+              -> C k ( I i rec ) ( Spline 'Closed () ( I i ( ℝ 2 ) ) )
+circleBrushFn _ mkI =
   D \ params ->
-    let r :: D k irec ( I i Double )
+    let r :: D k ( I i rec ) ( I i Double )
         r = runD ( var @_ @k $ Fin 1 ) params
-        mkPt :: Double -> Double -> D k irec ( I i ( ℝ 2 ) )
+        mkPt :: Double -> Double -> D k ( I i rec ) ( I i ( ℝ 2 ) )
         mkPt x y
           =   ( r `scaledBy` x ) *^ e_x
           ^+^ ( r `scaledBy` y ) *^ e_y
     in circleSpline mkPt
   where
-    e_x, e_y :: D k irec ( I i ( ℝ 2 ) )
+    e_x, e_y :: D k ( I i rec ) ( I i ( ℝ 2 ) )
     e_x = pure $ mkI $ ℝ2 1 0
     e_y = pure $ mkI $ ℝ2 0 1
 
     scaledBy d x = fmap ( mkI x * ) d
-{-# INLINEABLE circleBrush #-}
+{-# INLINEABLE circleBrushFn #-}
 
-ellipseBrush :: forall {t} (i :: t) k irec
-             . ( 3 <= RepDim irec
-               , Module
-                  ( D k irec ( I i Double ) )
-                  ( D k irec ( I i ( ℝ 2 ) ) )
-               , Module ( I i Double ) ( T ( I i Double ) )
-               , HasChainRule ( I i Double ) k irec
-               , Representable ( I i Double ) irec
-               , Applicative ( D k irec )
-               , Transcendental ( D k irec ( I i Double ) )
-                 -- TODO: make a synonym for the above...
-                 -- it seems DiffInterp isn't exactly right
-               )
-              => Proxy# i
-              -> ( forall a. a -> I i a )
-              -> C k irec ( Spline 'Closed () ( I i ( ℝ 2 ) ) )
-ellipseBrush _ mkI =
+ellipseBrushFn :: forall {t} (i :: t) k rec
+               . ( 3 <= RepDim ( I i rec )
+                 , ParamsICt k i rec
+                 )
+                => Proxy# i
+                -> ( forall a. a -> I i a )
+                -> C k ( I i rec ) ( Spline 'Closed () ( I i ( ℝ 2 ) ) )
+ellipseBrushFn _ mkI =
   D \ params ->
-    let a, b, phi :: D k irec ( I i Double )
-        a   = runD ( var @_ @k $ Fin 1 ) params
-        b   = runD ( var @_ @k $ Fin 2 ) params
-        phi = runD ( var @_ @k $ Fin 3 ) params
-        mkPt :: Double -> Double -> D k irec ( I i ( ℝ 2 ) )
+    let a, b :: D k ( I i rec ) ( I i Double )
+        a = runD ( var @_ @k $ Fin 1 ) params
+        b = runD ( var @_ @k $ Fin 2 ) params
+        mkPt :: Double -> Double -> D k ( I i rec ) ( I i ( ℝ 2 ) )
         mkPt x y
           = let !x' = a `scaledBy` x
                 !y' = b `scaledBy` y
--- {-
-            in
-                ( x' * cos phi - y' * sin phi ) *^ e_x
-            ^+^ ( y' * cos phi + x' * sin phi ) *^ e_y
--- -}
-{-
-                r  = sqrt ( x' ^ 2 + y' ^ 2 )
-                arctgt = atan ( y' / x' )
-                -- a and b are always strictly positive, so we can compute
-                -- the quadrant using only x and y, which are constants.
-                !theta
-                  | x > 0
-                  = arctgt
-                  | x < 0
-                  = if y >= 0 then arctgt + pi else arctgt - pi
-                  | otherwise
-                  = if y >= 0 then 0.5 * pi else -0.5 * pi
-                !phi' = phi + theta
-            in
-                ( r * cos phi' ) *^ e_x
-            ^+^ ( r * sin phi' ) *^ e_y
--}
-
+            in x' *^ e_x ^+^ y' *^ e_y
     in circleSpline mkPt
   where
-    e_x, e_y :: D k irec ( I i ( ℝ 2 ) )
+    e_x, e_y :: D k ( I i rec ) ( I i ( ℝ 2 ) )
     e_x = pure $ mkI $ ℝ2 1 0
     e_y = pure $ mkI $ ℝ2 0 1
 
     scaledBy d x = fmap ( mkI x * ) d
-{-# INLINEABLE ellipseBrush #-}
+{-# INLINEABLE ellipseBrushFn #-}
 
 --------------------------------------------------------------------------------
 -- Tear drop
@@ -168,27 +193,19 @@ tearHeight = 3 * sqrt 3 / 8
 sqrt3_over_2 :: Double
 sqrt3_over_2 = 0.5 * sqrt 3
 
-tearDropBrush :: forall {t} (i :: t) k irec
-             . ( Module
-                  ( D k irec ( I i Double ) )
-                  ( D k irec ( I i ( ℝ 2 ) ) )
-               , Module ( I i Double ) ( T ( I i Double ) )
-               , HasChainRule ( I i Double ) k irec
-               , Representable ( I i Double ) irec
-               , Applicative ( D k irec )
-               , Transcendental ( D k irec ( I i Double ) )
-               )
-              => Proxy# i
-              -> ( forall a. a -> I i a )
-              -> C k irec ( Spline 'Closed () ( I i ( ℝ 2 ) ) )
-tearDropBrush _ mkI =
+tearDropBrushFn :: forall {t} (i :: t) k rec
+                . ( 3 <= RepDim ( I i rec )
+                 , ParamsICt k i rec
+                 )
+                 => Proxy# i
+                 -> ( forall a. a -> I i a )
+                 -> C k ( I i rec ) ( Spline 'Closed () ( I i ( ℝ 2 ) ) )
+tearDropBrushFn _ mkI =
   D \ params ->
-    let w, h, phi :: D k irec ( I i Double )
-        w   = runD ( var @_ @k ( Fin 1 ) ) params
-        h   = runD ( var @_ @k ( Fin 2 ) ) params
-        phi = runD ( var @_ @k ( Fin 3 ) ) params
-
-        mkPt :: Double -> Double -> D k irec ( I i ( ℝ 2 ) )
+    let w, h :: D k ( I i rec ) ( I i Double )
+        w = runD ( var @_ @k ( Fin 1 ) ) params
+        h = runD ( var @_ @k ( Fin 2 ) ) params
+        mkPt :: Double -> Double -> D k ( I i rec ) ( I i ( ℝ 2 ) )
         mkPt x y
           -- 1. translate the teardrop so that the center of mass is at the origin
           -- 2. scale the teardrop so that it has the requested width/height
@@ -196,8 +213,8 @@ tearDropBrush _ mkI =
           = let !x' =   w `scaledBy` (x / tearWidth)
                 !y' = ( h `scaledBy` ( ( y - tearCenter ) / tearHeight) )
             in
-                ( x' * cos phi - y' * sin phi ) *^ e_x
-            ^+^ ( y' * cos phi + x' * sin phi ) *^ e_y
+                x'  *^ e_x
+            ^+^ y'  *^ e_y
 
     in sequenceA $
          Spline { splineStart  = mkPt 0 0
@@ -207,9 +224,9 @@ tearDropBrush _ mkI =
                     ( mkPt -0.5 sqrt3_over_2 )
                     BackToStart () }
   where
-    e_x, e_y :: D k irec ( I i ( ℝ 2 ) )
+    e_x, e_y :: D k ( I i rec ) ( I i ( ℝ 2 ) )
     e_x = pure $ mkI $ ℝ2 1 0
     e_y = pure $ mkI $ ℝ2 0 1
 
     scaledBy d x = fmap ( mkI x * ) d
-{-# INLINEABLE tearDropBrush #-}
+{-# INLINEABLE tearDropBrushFn #-}
diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs
index cba086e..7c64a2e 100644
--- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs
+++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs
@@ -12,7 +12,7 @@ module Math.Bezier.Stroke
 
     -- * Brush stroking
 
-  , brushStroke, envelopeEquation
+  , envelopeEquation
   , line, bezier2, bezier3
   , brushStrokeData, pathAndUsedParams
 
@@ -39,6 +39,8 @@ import Data.Fixed
   ( divMod' )
 import Data.Foldable
   ( for_ )
+import Data.Functor
+  ( (<&>) )
 import Data.Functor.Identity
   ( Identity(..) )
 import Data.List
@@ -49,11 +51,12 @@ import qualified Data.List.NonEmpty as NE
   ( cons, singleton, unzip )
 import Data.Maybe
   ( fromMaybe, isJust, listToMaybe, mapMaybe )
+import Data.Proxy
+  ( Proxy(..) )
 import Data.Semigroup
   ( sconcat )
 import GHC.Exts
   ( newMutVar#, runRW#, inline
-  , Proxy#, proxy#
   )
 import GHC.STRef
   ( STRef(..), readSTRef, writeSTRef )
@@ -105,6 +108,8 @@ import Control.Monad.Trans.Writer.CPS
   ( WriterT, execWriterT, runWriter, tell )
 
 -- MetaBrush
+import Calligraphy.Brushes
+  ( Brush(..) )
 import Math.Algebra.Dual
 import qualified Math.Bezier.Cubic     as Cubic
 import Math.Bezier.Cubic.Fit
@@ -131,6 +136,8 @@ import Math.Module
   ( Module(..), Inner((^.^)), Cross((×)), Interpolatable
   , lerp, convexCombination, strictlyParallel
   )
+import Math.Ring
+  ( Transcendental )
 import Math.Orientation
   ( Orientation(..), splineOrientation
   , between
@@ -250,20 +257,14 @@ computeStrokeOutline ::
   -> FitParameters
   -> ( ptData -> usedParams )
   -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
-  -> ( forall {t} k (i :: t)
-       .  DiffInterp k i brushParams
-       => Proxy# i
-       -> ( forall a. a -> I i a )
-       -> C k ( I i brushParams )
-              ( Spline Closed () ( I i ( ℝ 2 ) ) )
-      )
+  -> Brush brushParams
   -> Spline clo crvData ptData
   -> ST s
         ( Either ( SplinePts Closed ) ( SplinePts Closed, SplinePts Closed )
         , Seq FitPoint
         , [ Cusp ]
         )
-computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams brushFn spline@( Spline { splineStart = spt0 } ) = case ssplineType @clo of
+computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline@( Spline { splineStart = spt0 } ) = case ssplineType @clo of
   -- Open brush path with at least one segment.
   -- Need to add caps at both ends of the path.
   SOpen
@@ -365,7 +366,7 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru
   where
 
     outlineInfo :: ptData -> Curve Open crvData ptData -> OutlineInfo
-    outlineInfo = inline ( outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams brushFn )
+    outlineInfo = inline ( outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams brush )
 
     outlineFns :: Seq OutlineInfo
     outlineFns = go spt0 ( openCurves $ splineCurves ( adjustSplineType @Open spline ) )
@@ -379,7 +380,17 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru
           outlineInfo p0 crv :<| go ( openCurveEnd crv ) crvs
 
     brushShape :: ptData -> SplinePts Closed
-    brushShape pt = fun @Double ( brushFn @2 @() proxy# id ) $ toBrushParams $ ptParams pt
+    brushShape pt =
+      let Brush { brushShape = shapeFn, mbRotation = mbRot } = brush
+          brushParams = toBrushParams $ ptParams pt
+          shape = fun @Double shapeFn brushParams
+      in case mbRot of
+           Nothing -> shape
+           Just getθ ->
+            let θ = getθ brushParams
+                cosθ = cos θ
+                sinθ = sin θ
+            in fmap ( unT . rotate cosθ sinθ . T ) shape
 
     updateSpline :: ( T ( ℝ 2 ), T ( ℝ 2 ), T ( ℝ 2 ) ) -> ST s OutlineData
     updateSpline ( lastTgt, lastTgtFwd, lastTgtBwd )
@@ -528,17 +539,12 @@ outlineFunction
   -> Maybe ( RootIsolationOptions N 3 )
   -> ( ptData -> usedParams )
   -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
-  -> ( forall {t} k (i :: t)
-       .  DiffInterp k i brushParams
-       => Proxy# i
-       -> ( forall a. a -> I i a )
-       -> C k ( I i brushParams )
-              ( Spline Closed () ( I i ( ℝ 2 ) ) )
-      )
+  -> Brush brushParams
   -> ptData
   -> Curve Open crvData ptData
   -> OutlineInfo
-outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams brushFromParams = \ sp0 crv ->
+outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams
+  ( Brush { brushShape, brushShapeI, mbRotation } ) = \ sp0 crv ->
   let
 
     usedParams :: C 2 ( ℝ 1 ) usedParams
@@ -552,22 +558,19 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams brushFromParams =
       brushStrokeData @2 @brushParams
         coerce coerce
         path
-        ( chainRule @Double @2
-            usedParams
-            ( linear toBrushParams )
-        )
-        ( brushFromParams @2 @() proxy# id )
+        ( fmap toBrushParams usedParams )
+        brushShape
+        mbRotation
 
     curvesI :: 𝕀ℝ 1 -- t
             -> Seq ( 𝕀ℝ 1 {- s -} -> StrokeDatum 3 𝕀 )
-    curvesI = brushStrokeData @3 @brushParams
-               coerce coerce
-               pathI
-               ( chainRule @( 𝕀 Double ) @3
-                   usedParamsI
-                   ( linear ( nonDecreasing toBrushParams ) )
-               )
-               ( brushFromParams @3 @𝕀 proxy# singleton )
+    curvesI =
+      brushStrokeData @3 @brushParams
+        coerce coerce
+        pathI
+        ( fmap ( nonDecreasing toBrushParams ) usedParamsI )
+        brushShapeI
+        ( fmap nonDecreasing mbRotation )
 
     usedParamsI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀 usedParams )
     pathI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 )
@@ -585,7 +588,7 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams brushFromParams =
         D21 path_t path'_t _ = runD path t
         D21 params_t _ _     = runD usedParams t
         brush_t              = value @Double @2 @brushParams
-                             $ runD ( brushFromParams @2 @() proxy# id )
+                             $ runD brushShape
                              $ toBrushParams params_t
 
     ( potentialCusps, definiteCusps ) =
@@ -981,9 +984,12 @@ brushStrokeData :: forall k brushParams i arr
                    , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 )
                    , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 )
                    , D k ( I i ( ℝ 2 ) ) ~ D k ( ℝ 2 )
+                   , Transcendental ( I i Double )
+                   , Module ( I i Double ) ( T ( I i ( ℝ 1 ) ) )
                    , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) )
                    , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) )
                    , Show brushParams
+                   , Representable ( I i Double ) ( I i ( ℝ 2 ) ), RepDim ( I i ( ℝ 2 ) ) ~ 2
                    )
                 => ( I i Double -> I i ( ℝ 1 ) )
                 -> ( I i ( ℝ 1 ) -> I i Double )
@@ -993,8 +999,10 @@ brushStrokeData :: forall k brushParams i arr
                      -- ^ brush parameters
                 -> ( I i brushParams `arr` Spline Closed () ( I i ( ℝ 2 ) ) )
                      -- ^ brush from parameters
+                -> ( Maybe ( I i brushParams -> I i Double ) )
+                     -- ^ rotation parameter
                 -> ( I i ( ℝ 1 ) -> Seq ( I i ( ℝ 1 ) -> StrokeDatum k i ) )
-brushStrokeData co1 co2 path params brush =
+brushStrokeData co1 co2 path params brush mbBrushRotation =
   \ t ->
     let
         dpath_t :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) )
@@ -1009,35 +1017,17 @@ brushStrokeData co1 co2 path params brush =
         !dbrushes_t = force $ fmap ( uncurryD @k . chain @( I i Double ) @k dparams_t ) splines
           -- This is the crucial use of the chain rule.
 
-    in fmap ( mkStrokeDatum dpath_t ) dbrushes_t
+    in fmap ( mkStrokeDatum dpath_t dparams_t ) dbrushes_t
   where
 
     mkStrokeDatum :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) )
+                  -> D k ( I i ( ℝ 1 ) ) ( I i brushParams )
                   -> ( I i ( ℝ 1 ) -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) )
                   -> ( I i ( ℝ 1 ) -> StrokeDatum k i )
-    mkStrokeDatum dpath_t dbrush_t s =
+    mkStrokeDatum dpath_t dparams_t dbrush_t s =
       let dbrush_t_s = dbrush_t s
-          dstroke = brushStroke @k dpath_t dbrush_t_s
-          ( ee, 𝛿E𝛿sdcdt ) = envelopeEquation @k @i co1 dstroke
-      in  -- trace
-          --   ( unlines
-          --     [ "envelopeEquation:"
-          --     , "      t = " ++ show t
-          --     , "      s = " ++ show s
-          --     , "      c = " ++ show _c
-          --     , "  ∂c/∂t = " ++ show _𝛿c𝛿t
-          --     , "  ∂c/∂s = " ++ show _𝛿c𝛿s
-          --     , "      E = " ++ show ee
-          --     , "  ∂E/∂t = " ++ show _𝛿E𝛿t
-          --     , "  ∂E/∂s = " ++ show ee_s
-          --     , "  dc/dt = " ++ show dcdt ] ) $
-          StrokeDatum
-            { dpath  = dpath_t
-            , dbrush = dbrush_t_s
-            , dstroke
-            , ee
-            , 𝛿E𝛿sdcdt
-            }
+          mbRotation = mbBrushRotation <&> \ getTheta -> fmap getTheta dparams_t
+      in  envelopeEquation @k @i Proxy co1 dpath_t dbrush_t_s mbRotation
 
 --------------------------------------------------------------------------------
 -- Solving the envelolpe equation: root-finding.
@@ -1105,7 +1095,7 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok
       let (i, s) = fromDomain is in
       case evalStrokeDatum fs is of -- TODO: a bit redundant to have to compute this again...
         StrokeDatum
-          { dstroke
+          { stroke
           , ee = D12 ( ℝ1 _ee ) _ ( T ( ℝ1 𝛿E𝛿s ) )
           , 𝛿E𝛿sdcdt = D0 𝛿E𝛿sdcdt
           } ->
@@ -1119,7 +1109,7 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok
                      -> recip 𝛿E𝛿s' *^ 𝛿E𝛿sdcdt'
                 | otherwise
                 = recip 𝛿E𝛿s *^ 𝛿E𝛿sdcdt
-          in ( value @Double @2 @( ℝ 2 ) dstroke, dcdt )
+          in ( stroke, dcdt )
 
     evalStrokeDatum :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> ( Double -> StrokeDatum 2 () )
     evalStrokeDatum fs is =
@@ -1154,8 +1144,7 @@ cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) )
            -> Cusp
 cuspCoords eqs ( i, box )
   | StrokeDatum
-    { dpath
-    , dstroke = D22 { _D22_v = stroke } }
+    { dpath, stroke }
       <- ( eqs ( ℝ1 t_mid ) `Seq.index` i ) ( ℝ1 s_mid )
   = Cusp
   { cuspParameters   = ℝ2 t_mid s_mid
@@ -1168,8 +1157,6 @@ cuspCoords eqs ( i, box )
     t_mid = 0.5 * ( t_lo + t_hi )
     s_mid = 0.5 * ( s_lo + s_hi )
 
-type N = 2
-
 -- | Find cusps in the envelope for values of the parameters in
 -- \( 0 \leqslant t, s \leqslant 1 \), using interval arithmetic.
 --
@@ -1233,10 +1220,12 @@ findCuspsIn opts boxStrokeData initBoxes =
         let t = 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi )
             s = 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi )
             StrokeDatum
-              { dstroke =
-                D32
-                  { _D32_dx = T ( 𝕀 ( ℝ2 cx_t_lo cy_t_lo ) ( ℝ2 cx_t_hi cy_t_hi ) )
-                  , _D32_dy = T ( 𝕀 ( ℝ2 cx_s_lo cy_s_lo ) ( ℝ2 cx_s_hi cy_s_hi ) )}
+              { du =
+                D22
+                  { _D22_v = 𝕀 ( ℝ2 ux_lo uy_lo ) ( ℝ2 ux_hi uy_hi )}
+              , dv =
+                D22
+                  { _D22_v = 𝕀 ( ℝ2 vx_lo vy_lo ) ( ℝ2 vx_hi vy_hi ) }
               , ee =
                 D22
                   { _D22_dx = T ( 𝕀 ( ℝ1 ee_t_lo ) ( ℝ1 ee_t_hi ) )
@@ -1244,9 +1233,9 @@ findCuspsIn opts boxStrokeData initBoxes =
               } = ( boxStrokeData t `Seq.index` i ) s
         -- λ = ∂E/∂t / ∂E/∂s
             λ1 = 𝕀 ee_t_lo ee_t_hi `extendedDivide` 𝕀 ee_s_lo ee_s_hi
-        -- λ = ∂c/∂t / ∂c/∂s
-            λ2 = 𝕀 cx_t_lo cx_t_hi `extendedDivide` 𝕀 cx_s_lo cx_s_hi
-            λ3 = 𝕀 cy_t_lo cy_t_hi `extendedDivide` 𝕀 cy_s_lo cy_s_hi
+        -- λ = u / v
+            λ2 = 𝕀 ux_lo ux_hi `extendedDivide` 𝕀 vx_lo vx_hi
+            λ3 = 𝕀 uy_lo uy_hi `extendedDivide` 𝕀 vy_lo vy_hi
             λ = [ 𝕀 ( recip -0 ) ( recip 0 ) ]
                 `intersectMany` λ1
                 `intersectMany` λ2
@@ -1264,9 +1253,12 @@ findCuspsIn opts boxStrokeData initBoxes =
             λ = 𝕀 λ_lo λ_hi
 
             StrokeDatum
-              { dstroke =
-                 D32 { _D32_dx = T c_t, _D32_dy = T c_s
-                     , _D32_dxdx = T c_tt, _D32_dxdy = T c_ts, _D32_dydy = T c_ss }
+              { du =
+                 D22 { _D22_v = u
+                     , _D22_dx   = T u_t , _D22_dy   = T u_s }
+              , dv =
+                 D22 { _D22_v = v
+                     , _D22_dx   = T v_t , _D22_dy   = T v_s }
               , ee =
                 D22 { _D22_dx = T ee_t, _D22_dy = T ee_s
                     , _D22_dxdx = T ee_tt, _D22_dxdy = T ee_ts, _D22_dydy = T ee_ss }
@@ -1278,11 +1270,11 @@ findCuspsIn opts boxStrokeData initBoxes =
             𝕀 ( ℝ1 f1_s_lo ) ( ℝ1 f1_s_hi ) = unT $ λ *^ T ee_ss ^-^ T ee_ts
             𝕀 ( ℝ1 f1_λ_lo ) ( ℝ1 f1_λ_hi ) = ee_s
 
-            -- λ ∂c/∂s - ∂c/∂t = 0
-            𝕀 ( ℝ2 f2_lo   f3_lo   ) ( ℝ2 f2_hi   f3_hi   ) = unT $ λ *^ T c_s  ^-^ T c_t
-            𝕀 ( ℝ2 f2_t_lo f3_t_lo ) ( ℝ2 f2_t_hi f3_t_hi ) = unT $ λ *^ T c_ts ^-^ T c_tt
-            𝕀 ( ℝ2 f2_s_lo f3_s_lo ) ( ℝ2 f2_s_hi f3_s_hi ) = unT $ λ *^ T c_ss ^-^ T c_ts
-            𝕀 ( ℝ2 f2_λ_lo f3_λ_lo ) ( ℝ2 f2_λ_hi f3_λ_hi ) = c_s
+            -- λ v - u = 0
+            𝕀 ( ℝ2 f2_lo   f3_lo   ) ( ℝ2 f2_hi   f3_hi   ) = unT $ λ *^ T v   ^-^ T u
+            𝕀 ( ℝ2 f2_t_lo f3_t_lo ) ( ℝ2 f2_t_hi f3_t_hi ) = unT $ λ *^ T v_t ^-^ T u_t
+            𝕀 ( ℝ2 f2_s_lo f3_s_lo ) ( ℝ2 f2_s_hi f3_s_hi ) = unT $ λ *^ T v_s ^-^ T u_s
+            𝕀 ( ℝ2 f2_λ_lo f3_λ_lo ) ( ℝ2 f2_λ_hi f3_λ_hi ) = v
 
         in D13 ( 𝕀 ( ℝ3 f1_lo f2_lo f3_lo ) ( ℝ3 f1_hi f2_hi f3_hi ) )
                ( T $ 𝕀 ( ℝ3 f1_t_lo f2_t_lo f3_t_lo ) ( ℝ3 f1_t_hi f2_t_hi f3_t_hi ) )
diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs
index 9c69a6c..7c5482a 100644
--- a/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs
+++ b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs
@@ -10,11 +10,13 @@ module Math.Bezier.Stroke.EnvelopeEquation
   ) where
 
 -- base
-import Prelude hiding ( Num(..), (^) )
+import Prelude hiding ( Num(..), (^), pi, sin, cos )
 import Data.Kind
   ( Type, Constraint )
 import Data.List.NonEmpty
   ( NonEmpty(..) )
+import Data.Proxy
+  ( Proxy(..) )
 import GHC.TypeNats
   ( Nat, type (-) )
 
@@ -37,20 +39,26 @@ import Math.Ring
 --------------------------------------------------------------------------------
 
 -- | The value and derivative of a brush stroke at a given coordinate
--- \( (t_0, s_0) \), together with the value of the envelope equation at that
+-- \( (t, s) \), together with the value of the envelope equation at that
 -- point.
 type StrokeDatum :: Nat -> k -> Type
 data StrokeDatum k i
   = StrokeDatum
-  { -- | Path \( p(t_0) \).
+  { -- | Path \( p(t) \).
     dpath    :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) )
-    -- | Brush shape \( b(t_0, s_0) \).
+    -- | Brush shape \( b(t, s) \).
   , dbrush   :: D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) )
+    -- | (Optional) rotation angle \( \theta(t) \).
+  , mbRotation :: Maybe ( D k ( I i ( ℝ 1 ) ) ( I i Double ) )
 
-  -- Everything below can be computed in terms of the first two fields.
+   -- Everything below is computed in terms of the first three fields.
 
-    -- | Stroke \( c(t_0,s_0) = p(t_0) + b(t_0,s_0) \).
-  , dstroke  :: D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) )
+    -- | Stroke shape \( c(t,s) = p(t) + R(\theta(t)) b(t,s) \).
+  , stroke :: I i ( ℝ 2 )
+
+    -- | \( u(t,s) = R(-\theta(t)) \frac{\partial c}{\partial t} \),
+    --   \( v(t,s) = R(-\theta(t)) \frac{\partial c}{\partial s} \)
+  , du, dv :: D ( k - 1 ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) )
 
     -- | Envelope function
     --
@@ -108,19 +116,6 @@ class HasEnvelopeEquation k where
   uncurryD :: D k a ~ D k ( ℝ 1 )
            => D k ( ℝ 1 ) ( C k a b ) -> a -> D k ( ℝ 2 ) b
 
-  -- | A brush stroke, as described by the equation
-  --
-  -- \[ c(t,s) = p(t) + b(t,s) \]
-  --
-  -- where:
-  --
-  --   - \( p(t) \) is the path that the brush follows, and
-  --   - \( b(t,s) \) is the brush shape, as it varies along the path.
-  brushStroke :: Module r ( T v )
-              => D k ( ℝ 1 ) v -- ^ stroke path \( p(t) \)
-              -> D k ( ℝ 2 ) v -- ^ brush \( b(t,s) \)
-              -> D k ( ℝ 2 ) v
-
   -- | The envelope function
   --
   -- \[ E = \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s}, \]
@@ -136,12 +131,19 @@ class HasEnvelopeEquation k where
                    .  ( D ( k - 2 ) ( I i ( ℝ 2 ) ) ~ D ( k - 2 ) ( ℝ 2 )
                       , D ( k - 1 ) ( I i ( ℝ 2 ) ) ~ D ( k - 1 ) ( ℝ 2 )
                       , D k ( I i ( ℝ 2 ) ) ~ D k ( ℝ 2 )
+                      , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 )
+                      , Module ( I i Double ) ( T ( I i ( ℝ 1 ) ) )
                       , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) )
+                      , Transcendental ( I i Double )
+                      , Representable ( I i Double ) ( I i ( ℝ 2 ) )
+                      , RepDim ( I i ( ℝ 2 ) ) ~ 2
                       )
-                   => ( I i Double -> I i ( ℝ 1 ) )
+                   => Proxy i
+                   -> ( I i Double -> I i ( ℝ 1 ) )
+                   -> D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) )
                    -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) )
-                   -> ( D ( k - 1 ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 1 ) )
-                      , D ( k - 2 ) ( I i ( ℝ 2 ) ) ( T ( I i ( ℝ 2 ) ) ) )
+                   -> Maybe ( D k ( I i ( ℝ 1 ) ) ( I i Double ) )
+                   -> StrokeDatum k i
 
 instance HasBézier 2 () where
   line co ( Segment a b :: Segment b ) =
@@ -166,34 +168,85 @@ instance HasEnvelopeEquation 2 where
 
   uncurryD = uncurryD2
 
-  brushStroke ( D21 p p_t p_tt ) ( D22 b b_t b_s b_tt b_ts b_ss ) =
-    D22 ( unT $ T p ^+^ T b )
-       -- c = p + b
+  envelopeEquation ( _ :: Proxy i ) co
+    dp@( D21 ( T -> p ) p_t p_tt )
+    db@( D22 ( T -> b ) b_t b_s
+             b_tt b_ts b_ss )
+    mbRotation =
+      StrokeDatum
+        { dpath = dp
+        , dbrush = db
+        , mbRotation
+        , stroke = c
+        , du, dv, ee, 𝛿E𝛿sdcdt
+        }
+      where
+        (ee, 𝛿E𝛿sdcdt) =
+          let
+            D12 (T -> u) u_t u_s = du
+            D12 (T -> v) v_t v_s = dv
+            ee_val, ee_t, ee_s :: I i Double
+            ee_val= u × v
+            ee_t  = u_t × v
+                  + u × v_t
+            ee_s  = u_s × v
+                  + u × v_s
 
-     ( p_t ^+^ b_t ) b_s
-       -- ∂c/∂t = dp/dt + ∂b/∂t
-       -- ∂c/∂s = ∂b/∂s
+            𝛿E𝛿sdcdt_val = ee_s *^ u ^-^ ee_t *^ v
+          in ( D12
+                 ( co ee_val )
+                 ( T $ co ee_t ) ( T $ co ee_s )
+             , D0 𝛿E𝛿sdcdt_val
+             )
+        (c, du, dv) = case mbRotation of
+          Nothing ->
+            -- c(t,s) = p(t) + b(t,s)
+            ( unT $ p ^+^ b
+            , D12 ( unT $ p_t ^+^ b_t )
+                  ( p_tt ^+^ b_tt ) b_ts
+            , D12 ( unT $ b_s )
+                  b_ts b_ss
+            )
+          Just ( D21 θ ( T θ_t ) ( T θ_tt ) ) ->
+            -- c(t,s) = p(t) + R(θ(t)) b(t,s)
+            -- E = ∂c/∂t × ∂c/ds
+            --   = ( R(-θ(t)) ∂c/∂t ) × ( R(-θ(t)) ∂c/ds )
+            --   = ( R(-θ(t)) p'(t) + θ'(t) S b(t,s) + ∂b/∂t ) × ∂b/ds
+            --
+            let rot, rot' :: T ( I i ( ℝ 2 ) ) -> T ( I i ( ℝ 2 ) )
+                cosθ = cos θ
+                sinθ = sin θ
+                -- rot = R(-θ), rot' = R'(-θ)
+                -- NB: rot' is not the derivative of f(θ) = R(-θ)
+                rot = rotate cosθ -sinθ
+                rot' = rotate sinθ cosθ
+                swap :: T ( I i ( ℝ 2 ) ) -> T ( I i ( ℝ 2 ) )
+                swap ( T xy ) =
+                  let x = xy `index` Fin 1
+                      y = xy `index` Fin 2
+                  in T $ tabulate \ case
+                        Fin 1 -> -y
+                        _     ->  x
 
-     ( p_tt ^+^ b_tt ) b_ts b_ss
-       -- ∂²c/∂t²  = d²p/dt² + ∂²b/∂t²
-       -- ∂²c/∂t∂s = ∂²b/∂t∂s
-       -- ∂²c/∂s²  = ∂²b/∂s²
+                u, v, u_t, u_s, v_t, v_s :: T ( I i ( ℝ 2 ) )
+                u = rot p_t ^+^ θ_t *^ swap b ^+^ b_t
+                v = b_s
+                u_t  = ( -θ_t *^ rot' p_t
+                       ^+^ rot p_tt
+                       )
+                     ^+^
+                       (   θ_tt *^ swap b
+                       ^+^ θ_t  *^ swap b_t
+                       )
+                     ^+^ b_tt
+                u_s  = θ_t *^ swap b_s ^+^ b_ts
+                v_t  = b_ts
+                v_s  = b_ss
 
-  envelopeEquation co ( D22 _ c_t c_s c_tt c_ts c_ss ) =
-    let ee   = c_t × c_s
-        ee_t = c_tt × c_s + c_t × c_ts
-        ee_s = c_ts × c_s + c_t × c_ss
-        𝛿E𝛿sdcdt = ee_s *^ c_t ^-^ ee_t *^ c_s
-          -- TODO: we get c_t * c_t and c_s * c_s terms...
-          -- These could be squares (better with interval arithmetic)?
-    in ( D12 ( co ee ) ( T $ co ee_t ) ( T $ co ee_s )
-       , D0 𝛿E𝛿sdcdt )
-    -- Computation of total derivative dc/dt:
-    --
-    -- dc/dt = ∂c/∂t + ∂c/∂s ∂s/∂t
-    -- ∂s/∂t = - ∂E/∂t / ∂E/∂s
-    --
-    -- ∂E/∂s dc/dt = ∂E/∂s ∂c/∂t - ∂E/∂t ∂c/∂s.
+            in ( unT $ p ^+^ rotate cosθ sinθ b
+               , D12 ( unT u ) u_t u_s
+               , D12 ( unT v ) v_t v_s
+               )
 
 instance HasBézier 3 () where
 
@@ -222,42 +275,123 @@ instance HasEnvelopeEquation 3 where
 
   uncurryD = uncurryD3
 
-  brushStroke
-    ( D31 p p_t p_tt p_ttt )
-    ( D32 b b_t b_s b_tt b_ts b_ss b_ttt b_tts b_tss b_sss ) =
-      D32
-        ( unT $ T p ^+^ T b )
-        ( p_t ^+^ b_t ) b_s
-        ( p_tt ^+^ b_tt ) b_ts b_ss
-        ( p_ttt ^+^ b_ttt ) b_tts b_tss b_sss
+  envelopeEquation ( _ :: Proxy i ) co
+    dp@( D31 ( T -> p ) p_t p_tt p_ttt )
+    db@( D32 ( T -> b ) b_t b_s
+             b_tt b_ts b_ss
+             b_ttt b_tts b_tss b_sss )
+    mbRotation =
+      StrokeDatum
+        { dpath = dp
+        , dbrush = db
+        , mbRotation
+        , stroke = c
+        , du, dv, ee, 𝛿E𝛿sdcdt
+        }
+      where
+        (ee, 𝛿E𝛿sdcdt) =
+          let
+            D22 (T -> u) u_t u_s u_tt u_ts u_ss = du
+            D22 (T -> v) v_t v_s v_tt v_ts v_ss = dv
+            ee_val, ee_t, ee_s, ee_tt, ee_ts, ee_ss :: I i Double
+            ee_val= u × v
+            ee_t  = u_t × v
+                  + u × v_t
+            ee_s  = u_s × v
+                  + u × v_s
+            ee_tt = u_tt × v
+                  + 2 * ( u_t × v_t )
+                  + u × v_tt
+            ee_ts = u_ts × v
+                  + u_t  × v_s
+                  + u_s  × v_t
+                  + u    × v_ts
+            ee_ss = u_ss × v
+                  + 2 * ( u_s × v_s )
+                  + u × v_ss
 
-  envelopeEquation co
-    ( D32 _ c_t c_s
-          c_tt c_ts c_ss
-          c_ttt c_tts c_tss c_sss )
-    = let ee    = c_t × c_s
-          ee_t  = c_tt × c_s + c_t × c_ts
-          ee_s  = c_ts × c_s + c_t × c_ss
-          ee_tt = c_ttt × c_s
-                + c_tt  × c_ts * 2
-                + c_t   × c_tts
-          ee_ts = c_tts × c_s
-                + c_tt  × c_ss
-             -- + c_ts  × c_ts -- cancels out
-                + c_t   × c_tss
-          ee_ss = c_tss × c_s
-                + c_ts  × c_ss * 2
-                + c_t   × c_sss
-          𝛿E𝛿sdcdt   = ee_s *^ c_t ^-^ ee_t *^ c_s
-          𝛿E𝛿sdcdt_t =    ee_ts *^ c_t ^+^ ee_s *^ c_tt
-                    ^-^ ( ee_tt *^ c_s ^+^ ee_t *^ c_ts )
-          𝛿E𝛿sdcdt_s =   ee_ss *^ c_t ^+^ ee_s *^ c_ts
-                   ^-^ ( ee_ts *^ c_s ^+^ ee_t *^ c_ss )
-      in ( D22
-            ( co ee )
-            ( T $ co ee_t ) ( T $ co ee_s )
-            ( T $ co ee_tt) ( T $ co ee_ts ) ( T $ co ee_ss )
-         , D12 𝛿E𝛿sdcdt ( T 𝛿E𝛿sdcdt_t ) ( T 𝛿E𝛿sdcdt_s ) )
+            𝛿E𝛿sdcdt_val = ee_s *^ u ^-^ ee_t *^ v
+            𝛿E𝛿sdcdt_t =    ee_ts *^ u ^+^ ee_s *^ u_t
+                      ^-^ ( ee_tt *^ v ^+^ ee_t *^ v_t )
+            𝛿E𝛿sdcdt_s =   ee_ss *^ u ^+^ ee_s *^ u_s
+                     ^-^ ( ee_ts *^ v ^+^ ee_t *^ v_s )
+          in ( D22
+                 ( co ee_val )
+                 ( T $ co ee_t ) ( T $ co ee_s )
+                 ( T $ co ee_tt) ( T $ co ee_ts ) ( T $ co ee_ss )
+             , D12 𝛿E𝛿sdcdt_val ( T 𝛿E𝛿sdcdt_t ) ( T 𝛿E𝛿sdcdt_s )
+             )
+        (c, du, dv) = case mbRotation of
+          Nothing ->
+            -- c(t,s) = p(t) + b(t,s)
+            ( unT $ p ^+^ b
+            , D22 ( unT $ p_t ^+^ b_t )
+                  ( p_tt ^+^ b_tt ) b_ts
+                  ( p_ttt ^+^ b_ttt ) b_tts b_tss
+            , D22 ( unT $ b_s )
+                  b_ts b_ss
+                  b_tts b_tss b_sss
+            )
+          Just ( D31 θ ( T θ_t ) ( T θ_tt ) ( T θ_ttt ) ) ->
+            -- c(t,s) = p(t) + R(θ(t)) b(t,s)
+            -- E = ∂c/∂t × ∂c/ds
+            --   = ( R(-θ(t)) ∂c/∂t ) × ( R(-θ(t)) ∂c/ds )
+            --   = ( R(-θ(t)) p'(t) + θ'(t) S b(t,s) + ∂b/∂t ) × ∂b/ds
+            --
+            let rot, rot', rot'' :: T ( I i ( ℝ 2 ) ) -> T ( I i ( ℝ 2 ) )
+                cosθ = cos θ
+                sinθ = sin θ
+                -- rot = R(-θ), rot' = R'(-θ), rot'' = R''(-θ)
+                -- NB: rot' is not the derivative of f(θ) = R(-θ)
+                rot = rotate cosθ -sinθ
+                rot' = rotate sinθ cosθ
+                rot'' z = -1 *^ rot z
+                swap :: T ( I i ( ℝ 2 ) ) -> T ( I i ( ℝ 2 ) )
+                swap ( T xy ) =
+                  let x = xy `index` Fin 1
+                      y = xy `index` Fin 2
+                  in T $ tabulate \ case
+                        Fin 1 -> -y
+                        _     ->  x
+
+                u, v, u_t, u_s, u_tt, u_ts, u_ss, v_t, v_s, v_tt, v_ts, v_ss :: T ( I i ( ℝ 2 ) )
+
+                u = rot p_t ^+^ θ_t *^ swap b ^+^ b_t
+                v = b_s
+                u_t  = ( -θ_t *^ rot' p_t
+                       ^+^ rot p_tt
+                       )
+                     ^+^
+                       (   θ_tt *^ swap b
+                       ^+^ θ_t  *^ swap b_t
+                       )
+                     ^+^ b_tt
+                u_s  = θ_t *^ swap b_s ^+^ b_ts
+                u_tt = (   ( θ_t ^ 2 ) *^ rot'' p_t
+                       ^-^ θ_tt *^ rot' p_t
+                       ^+^ ( 2 * θ_t ) *^ rot' p_tt
+                       ^+^ rot p_ttt
+                       )
+                     ^+^ ( θ_ttt *^ swap b
+                       ^+^ ( 2 * θ_tt ) *^ swap b_t
+                       ^+^ θ_t *^ swap b_tt
+                       )
+                     ^+^ b_ttt
+                u_ts =  θ_tt *^ swap b_s
+                    ^+^ θ_t *^ swap b_ts
+                    ^+^ b_tts
+                u_ss = θ_t *^ swap b_ss
+                    ^+^ b_tss
+                v_t  = b_ts
+                v_s  = b_ss
+                v_tt = b_tts
+                v_ts = b_tss
+                v_ss = b_sss
+
+            in ( unT $ p ^+^ rotate cosθ sinθ b
+               , D22 ( unT u ) u_t u_s u_tt u_ts u_ss
+               , D22 ( unT v ) v_t v_s v_tt v_ts v_ss
+               )
 
 instance HasBézier 3 𝕀 where
 
@@ -292,7 +426,7 @@ are convex combinations
 
   b_0(t) [p_0,q_0] + b_1(t) [p_1,q_1] + ... + b_n(t) [p_n,q_n]
 
-    -- Here b_1, ..., b_n are Bernstein polynomials.
+    -- Here b_0, ..., b_n are Bernstein polynomials.
 
 This means that the minimum value attained by the Bézier curve as we vary
 both the time parameter and the values of the points within their respective
@@ -334,31 +468,3 @@ evaluateQuadratic bez t =
       maxs = fmap (Quadratic.bezier @( T Double ) sup_bez)
            $ inf t :| ( sup t : filter ( `inside` t ) ( Quadratic.extrema sup_bez ) )
   in 𝕀 ( minimum mins ) ( maximum maxs )
-
-{-
-
-evaluateCubic :: Cubic.Bezier ( 𝕀 Double ) -> 𝕀 Double -> 𝕀 Double
-evaluateCubic bez t =
-  -- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t ⊆ [0,1]
-  let inf_bez = Cubic.restrict @( T Double ) ( fmap inf bez ) ( inf t, sup t )
-      sup_bez = Cubic.restrict @( T Double ) ( fmap sup bez ) ( inf t, sup t )
-      mins = fmap (Cubic.bezier @( T Double ) inf_bez)
-           $ 0 :| ( 1 : Cubic.extrema inf_bez )
-      maxs = fmap (Cubic.bezier @( T Double ) sup_bez)
-           $ 0 :| ( 1 : Cubic.extrema sup_bez )
-  in 𝕀 ( minimum mins ) ( maximum maxs )
-
--- | Evaluate a quadratic Bézier curve, when both the coefficients and the
--- parameter are intervals.
-evaluateQuadratic :: Quadratic.Bezier ( 𝕀 Double ) -> 𝕀 Double -> 𝕀 Double
-evaluateQuadratic bez t =
-  -- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t ⊆ [0,1]
-  let inf_bez = Quadratic.restrict @( T Double ) ( fmap inf bez ) ( inf t, sup t )
-      sup_bez = Quadratic.restrict @( T Double ) ( fmap sup bez ) ( inf t, sup t )
-      mins = fmap (Quadratic.bezier @( T Double ) inf_bez)
-           $ 0 :| ( 1 : Quadratic.extrema inf_bez )
-      maxs = fmap (Quadratic.bezier @( T Double ) sup_bez)
-           $ 0 :| ( 1 : Quadratic.extrema sup_bez )
-  in 𝕀 ( minimum mins ) ( maximum maxs )
-
--}
\ No newline at end of file
diff --git a/brush-strokes/src/lib/Math/Linear.hs b/brush-strokes/src/lib/Math/Linear.hs
index c5f0554..7225b35 100644
--- a/brush-strokes/src/lib/Math/Linear.hs
+++ b/brush-strokes/src/lib/Math/Linear.hs
@@ -15,6 +15,8 @@ module Math.Linear
   , Representable(..), set, injection, projection
   , Vec(..), (!), find, zipIndices
 
+  , rotate
+
   ) where
 
 -- base
@@ -52,6 +54,9 @@ import Data.Group.Generics
 
 -- brush-strokes
 import Math.Linear.Internal
+import Math.Ring
+  ( Ring )
+import qualified Math.Ring as Ring
 
 --------------------------------------------------------------------------------
 
@@ -151,3 +156,20 @@ zipIndices ( Vec v ) = zipIndices_ 1 v
     zipIndices_ ::  Word -> [ a ] -> [ ( Fin n, a ) ]
     zipIndices_ _ [] = []
     zipIndices_ i (a : as) = ( Fin i, a ) : zipIndices_ ( i + 1 ) as
+
+--------------------------------------------------------------------------------
+
+-- | Rotate a vector by the given angle (counter-clockwise),
+-- given the cosine and sine of the angle (in that order)
+rotate :: ( Representable r m, RepDim m ~ 2, Ring r )
+       => r -- \( \cos \theta \)
+       -> r -- \( \sin \theta \)
+       -> T m
+       -> T m
+rotate cosθ sinθ ( T xy ) =
+  let x = xy `index` Fin 1
+      y = xy `index` Fin 2
+  in  T $ tabulate \ case
+        Fin 1 -> x Ring.* cosθ Ring.- y Ring.* sinθ
+        _     -> y Ring.* cosθ Ring.+ x Ring.* sinθ
+{-# INLINEABLE rotate #-}
diff --git a/brush-strokes/src/lib/Math/Linear/Internal.hs b/brush-strokes/src/lib/Math/Linear/Internal.hs
index d01d724..0459a5e 100644
--- a/brush-strokes/src/lib/Math/Linear/Internal.hs
+++ b/brush-strokes/src/lib/Math/Linear/Internal.hs
@@ -45,15 +45,15 @@ data instance ℝ 0 = ℝ0
 newtype instance ℝ 1 = ℝ1 { unℝ1 :: Double }
   deriving stock   ( Generic )
   deriving newtype ( Eq, Ord, NFData )
-data    instance ℝ 2 = ℝ2 { _ℝ2_x, _ℝ2_y :: {-# UNPACK #-} !Double }
+data    instance ℝ 2 = ℝ2 { _ℝ2_x, _ℝ2_y :: Double }
   deriving stock    Generic
   deriving anyclass NFData
   deriving stock    ( Eq, Ord )
-data    instance ℝ 3 = ℝ3 { _ℝ3_x, _ℝ3_y, _ℝ3_z :: {-# UNPACK #-} !Double }
+data    instance ℝ 3 = ℝ3 { _ℝ3_x, _ℝ3_y, _ℝ3_z :: Double }
   deriving stock    Generic
   deriving anyclass NFData
   deriving stock    ( Eq, Ord )
-data    instance ℝ 4 = ℝ4 { _ℝ4_x, _ℝ4_y, _ℝ4_z, _ℝ4_w :: {-# UNPACK #-} !Double }
+data    instance ℝ 4 = ℝ4 { _ℝ4_x, _ℝ4_y, _ℝ4_z, _ℝ4_w :: Double }
   deriving stock    Generic
   deriving anyclass NFData
   deriving stock    ( Eq, Ord )
diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs
index 7613a52..8c364b2 100644
--- a/brush-strokes/src/lib/Math/Root/Isolation.hs
+++ b/brush-strokes/src/lib/Math/Root/Isolation.hs
@@ -24,6 +24,10 @@ module Math.Root.Isolation
     -- ** Trees recording search space of root isolation algorithms
   , RootIsolationTree(..), showRootIsolationTree
   , RootIsolationStep(..)
+
+    -- * Hack for changing between 2 and 3 d formulations
+    -- for my personal testing
+  , N
   )
   where
 
@@ -130,7 +134,8 @@ showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")"
 
 type Box n = 𝕀ℝ n
 type BoxHistory n = [ NE.NonEmpty ( RootIsolationStep, Box n ) ]
-type BoxCt n d = ( n ~ 2, d ~ 3 )
+type N = 2
+type BoxCt n d = ( n ~ N, d ~ 3 )
 {-
   ( Show ( 𝕀ℝ n ), Show ( ℝ n )
   , Eq ( ℝ n )
@@ -167,9 +172,9 @@ data RootIsolationAlgorithm n d
   -- | Gauss–Seidel step with the given preconditioning method.
   | GaussSeidel !( GaussSeidelOptions n d )
   -- | @box(1)@-consistency.
-  | Box1 !Box1Options
+  | Box1 !( Box1Options n d )
   -- | @box(2)@-consistency.
-  | Box2 !Box2Options
+  | Box2 !( Box2Options n d )
 
 -- | Options for the bisection method.
 type BisectionOptions :: Nat -> Nat -> Type
@@ -191,15 +196,20 @@ data GaussSeidelOptions n d =
     , gsDims :: ( 𝕀ℝ d -> 𝕀ℝ n ) }
 
 -- | Options for the @box(1)@-consistency method.
-data Box1Options =
+data Box1Options n d =
   Box1Options
-    { box1EpsEq :: !Double }
+    { box1EpsEq :: !Double
+    , box1CoordsToNarrow :: [ Fin n ]
+    , box1EqsToUse :: [ Fin d ]
+    }
 
 -- | Options for the @box(2)@-consistency method.
-data Box2Options =
+data Box2Options n d =
   Box2Options
     { box2EpsEq :: !Double
     , box2LambdaMin :: !Double
+    , box2CoordsToNarrow :: [ Fin n ]
+    , box2EqsToUse :: [ Fin d ]
     }
 
 defaultRootIsolationOptions :: BoxCt n d => RootIsolationOptions n d
@@ -236,8 +246,16 @@ defaultRootIsolationAlgorithms minWidth narrowAbs box history
       -- Otherwise, do a normal round.
       -- Currently: we try an interval Gauss–Seidel step followed by box(1)-consistency.
       _ -> GaussSeidel defaultGaussSeidelOptions
-           NE.:| [ Box1 ( Box1Options { box1EpsEq = narrowAbs } ) ]
+           NE.:| [ Box1 box1Options ]
   where
+    box1Options :: Box1Options n d
+    box1Options =
+      Box1Options
+        { box1EpsEq = narrowAbs
+        , box1CoordsToNarrow = toList $ universe @n -- [ Fin 1, Fin 2 ]
+        , box1EqsToUse = toList $ universe @d
+        }
+
     -- Did we reduce the box width by at least "narrowAbs" in at least one of the dimensions?
     sufficientlySmallerThan :: Box n -> Box n -> Bool
     b1 `sufficientlySmallerThan` b2 =
@@ -247,12 +265,12 @@ defaultRootIsolationAlgorithms minWidth narrowAbs box history
           <*> coordinates b2
 {-# INLINEABLE defaultRootIsolationAlgorithms #-}
 
-defaultGaussSeidelOptions :: GaussSeidelOptions 2 3
+defaultGaussSeidelOptions :: GaussSeidelOptions N 3
 defaultGaussSeidelOptions =
   GaussSeidelOptions
     { gsPreconditioner = InverseMidJacobian
     , gsDims = \ ( 𝕀 ( ℝ3 _a_lo b_lo c_lo ) ( ℝ3 _a_hi b_hi c_hi ) )
-             -> 𝕀 ( ℝ2 b_lo c_lo ) ( ℝ2 b_hi c_hi )
+               -> 𝕀 ( ℝ2 b_lo c_lo ) ( ℝ2 b_hi c_hi )
 }
 
 defaultBisectionOptions
@@ -431,10 +449,10 @@ doStrategy roundHistory previousRoundsHistory eqs minWidth algo box =
     GaussSeidel gsOptions -> do
       boxes <- intervalGaussSeidel gsOptions eqs box
       return ( GaussSeidelStep, boxes )
-    Box1 ( Box1Options { box1EpsEq } ) ->
-      return ( Box1Step, makeBox1Consistent eqs minWidth box1EpsEq box )
-    Box2 ( Box2Options { box2LambdaMin, box2EpsEq } ) ->
-      return ( Box2Step, [ makeBox2Consistent eqs minWidth box2EpsEq box2LambdaMin box ] )
+    Box1 box1Options ->
+      return ( Box1Step, makeBox1Consistent minWidth box1Options eqs box )
+    Box2 box2Options ->
+      return ( Box2Step, [ makeBox2Consistent minWidth box2Options eqs box ] )
     Bisection ( BisectionOptions { canHaveSols, fallbackBisectionDim } ) -> do
       let ( boxes, ( whatBis, mid ) ) = bisect ( canHaveSols eqs ) ( fallbackBisectionDim roundHistory previousRoundsHistory eqs ) box
       return ( BisectionStep whatBis mid, boxes )
@@ -625,23 +643,29 @@ data Preconditioner
 -- "Presentation of a highly tuned multithreaded interval solver for underdetermined and well-determined nonlinear systems"
 makeBox1Consistent
   :: BoxCt n d
-  => ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
-  -> Double -> Double
+  => Double -> Box1Options n d
+  -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
   -> Box n -> [ Box n ]
-makeBox1Consistent eqs minWidth epsEq x =
+makeBox1Consistent minWidth box1Options eqs x =
   ( `State.evalState` False ) $
-    pipeFunctionsWhileTrue ( allNarrowingOperators epsEq minWidth eqs ) x
+    pipeFunctionsWhileTrue ( allNarrowingOperators minWidth box1Options eqs ) x
 
 -- | An implementation of "bound-consistency" from the paper
 -- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
 makeBox2Consistent
   :: forall n d
   .  BoxCt n d
-  => ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
-  -> Double -> Double -> Double
+  => Double
+  -> Box2Options n d
+  -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
   -> Box n -> Box n
-makeBox2Consistent eqs minWidth epsEq lambdaMin x0 = ( `State.evalState` False ) $ doLoop 0.25 x0
+makeBox2Consistent minWidth (Box2Options epsEq lambdaMin coordsToNarrow eqsToUse) eqs x0
+  = ( `State.evalState` False ) $ doLoop 0.25 x0
   where
+    box1Options :: Box1Options n d
+    box1Options = Box1Options epsEq coordsToNarrow eqsToUse
+    doBox1 :: Box n -> [ Box n ]
+    doBox1 = makeBox1Consistent minWidth box1Options eqs
     doLoop :: Double -> Box n -> State Bool ( Box n )
     doLoop lambda x = do
       x'' <- forEachDim @n x $ \ i ->
@@ -660,11 +684,11 @@ makeBox2Consistent eqs minWidth epsEq lambdaMin x0 = ( `State.evalState` False )
           c1 = ( 1 - lambda ) * x_inf + lambda * x_sup
           c2 = lambda * x_inf + ( 1 - lambda ) * x_sup
           x'_inf =
-            case makeBox1Consistent eqs minWidth epsEq ( setter ( 𝕀 x_inf c1 ) box ) of
+            case doBox1 ( setter ( 𝕀 x_inf c1 ) box ) of
               []  -> c1
               x's -> minimum $ map ( inf . getter ) x's
           x'_sup =
-            case makeBox1Consistent eqs minWidth epsEq ( setter ( 𝕀 c2 x_sup ) box ) of
+            case doBox1 ( setter ( 𝕀 c2 x_sup ) box ) of
               []  -> c2
               x's -> maximum $ map ( sup . getter ) x's
           x' = 𝕀 x'_inf x'_sup
@@ -832,10 +856,10 @@ allNarrowingOperators
   :: forall n d
   .  BoxCt n d
   => Double
-  -> Double
+  -> Box1Options n d
   -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
   -> [ Box n -> State Bool [ Box n ] ]
-allNarrowingOperators eps_eq eps_bis eqs =
+allNarrowingOperators eps_bis ( Box1Options eps_eq coordsToNarrow eqsToUse ) eqs =
   [ \ cand ->
     let getter = ( `index` coordIndex )
         setter = set coordIndex
@@ -848,8 +872,8 @@ allNarrowingOperators eps_eq eps_bis eqs =
   | narrowFn <- [ leftNarrow, rightNarrow ]
   , ( coordIndex, fn ) <-
       [ ( i, ff' i d )
-      | i <- toList $ universe @n
-      , d <- toList $ universe @d
+      | i <- coordsToNarrow
+      , d <- eqsToUse
       ]
   ]
   where
diff --git a/src/app/MetaBrush/Application.hs b/src/app/MetaBrush/Application.hs
index 08d9617..946b0b1 100644
--- a/src/app/MetaBrush/Application.hs
+++ b/src/app/MetaBrush/Application.hs
@@ -73,13 +73,16 @@ import Control.Monad.Trans.Reader
   ( runReaderT )
 
 -- MetaBrush
+import Math.Root.Isolation
+  ( RootIsolationOptions(..), defaultRootIsolationOptions
+  , N
+  )
 import Math.Bezier.Cubic.Fit
   ( FitParameters(..) )
 import Math.Bezier.Spline
   ( Spline(..), Curves(..), Curve(..), NextPoint(..) )
 import Math.Bezier.Stroke
   ( RootSolvingAlgorithm(..)
-  , CuspFindingOptions(..), Preconditioner(..)
   , invalidateCache
   )
 import Math.Linear
@@ -172,9 +175,9 @@ runApplication application = do
               Spline
                 { splineStart  = mkPoint ( ℝ2 0 0 ) 10 25 0
                 , splineCurves = OpenCurves $ Seq.fromList
-                  [ LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 100 0 ) 15 40 (0.1 * pi) ), curveData = invalidateCache undefined }
-                  --, LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10  10 )  8 5 ( pi / 4 ) ), curveData = invalidateCache undefined }
-                  --, LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10 -20 ) 10 7 ( pi / 2 ) ), curveData = invalidateCache undefined }
+                  [ LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 100 0 ) 15 40 0 ), curveData = invalidateCache undefined }
+                  , LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10  10 )  8 5 ( pi / 4 ) ), curveData = invalidateCache undefined }
+                  , LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10 -20 ) 10 7 ( pi / 2 ) ), curveData = invalidateCache undefined }
                   ]
                 }
             }
@@ -215,8 +218,8 @@ runApplication application = do
                             --HalleyM2
                             NewtonRaphson
                               { maxIters = 20, precision = 8 }
-  cuspFindingOptionsTVar <- STM.newTVarIO @( Maybe CuspFindingOptions ) $
-                              Just defaultCuspFindingOptions
+  cuspFindingOptionsTVar <- STM.newTVarIO @( Maybe ( RootIsolationOptions N 3 ) ) $
+                              Just defaultRootIsolationOptions
 
   -- Put all these stateful variables in a record for conciseness.
   let
diff --git a/src/app/MetaBrush/Render/Document.hs b/src/app/MetaBrush/Render/Document.hs
index 7fc2b1f..4cb9f7f 100644
--- a/src/app/MetaBrush/Render/Document.hs
+++ b/src/app/MetaBrush/Render/Document.hs
@@ -19,8 +19,6 @@ import Data.Functor.Compose
   ( Compose(..) )
 import Data.Int
   ( Int32 )
-import GHC.Exts
-  ( proxy# )
 import GHC.Generics
   ( Generic, Generic1, Generically1(..) )
 
@@ -56,6 +54,8 @@ import Control.Monad.Trans.State.Strict
   ( StateT, evalStateT, get, put )
 
 -- MetaBrush
+import Calligraphy.Brushes
+  ( Brush(..) )
 import Math.Algebra.Dual
   ( D2𝔸1(..), fun )
 import qualified Math.Bezier.Cubic     as Cubic
@@ -74,16 +74,20 @@ import Math.Bezier.Spline
 import Math.Bezier.Stroke
   ( Cusp(..), CachedStroke(..), invalidateCache
   , computeStrokeOutline
-  , RootSolvingAlgorithm, CuspFindingOptions
+  , RootSolvingAlgorithm
   )
 import Math.Linear
-  ( ℝ(..), T(..) )
+  ( ℝ(..), T(..)
+  , rotate
+  )
 import Math.Module
   ( Module((*^)), normalise )
+import Math.Root.Isolation
+  ( RootIsolationOptions, N )
 import MetaBrush.Asset.Colours
   ( Colours, ColourRecord(..) )
 import MetaBrush.Brush
-  ( Brush(..), WithParams(..) )
+  ( NamedBrush(..), WithParams(..) )
 import MetaBrush.Context
   ( Modifier(..)
   , HoldAction(..), PartialPath(..)
@@ -149,7 +153,7 @@ blankRender _ = pure ()
 
 getDocumentRender
   :: Colours
-  -> RootSolvingAlgorithm -> Maybe CuspFindingOptions -> FitParameters
+  -> RootSolvingAlgorithm -> Maybe ( RootIsolationOptions N 3 ) -> FitParameters
   -> Mode -> Bool
   -> Set Modifier -> Maybe ( ℝ 2 ) -> Maybe HoldAction -> Maybe PartialPath
   -> Document
@@ -285,22 +289,22 @@ instance NFData StrokeRenderData where
 -- - Otherwise, this consists of the underlying spline path only.
 strokeRenderData
   :: RootSolvingAlgorithm
-  -> Maybe CuspFindingOptions
+  -> Maybe ( RootIsolationOptions N 3 )
   -> FitParameters
   -> Stroke
   -> Maybe ( ST RealWorld StrokeRenderData )
 strokeRenderData rootAlgo mbCuspOptions fitParams
   ( Stroke
     { strokeSpline = spline :: StrokeSpline clo ( Record pointFields )
-    , strokeBrush  = ( strokeBrush :: Maybe ( Brush brushFields ) )
+    , strokeBrush  = ( strokeBrush :: Maybe ( NamedBrush brushFields ) )
     , ..
     }
   ) | strokeVisible
     = Just $ case strokeBrush of
-        Just ( BrushData { brushFunction = fn } )
+        Just ( NamedBrush { brushFunction = fn } )
           | WithParams
              { defaultParams = brush_defaults
-             , withParams    = brushFn
+             , withParams    = brush@( Brush { brushShape, mbRotation = mbRot } )
              } <- fn
           -> -- This is the key place where we need to perform impedance matching
              -- between the collection of parameters supplied along a stroke and
@@ -315,15 +319,26 @@ strokeRenderData rootAlgo mbCuspOptions fitParams
                   -- Compute the outline using the brush function.
                     ( outline, fitPts, cusps ) <-
                       computeStrokeOutline @clo rootAlgo mbCuspOptions fitParams
-                        ( toUsedParams . brushParams ) embedUsedParams brushFn
-                        spline
+                        ( toUsedParams . brushParams ) embedUsedParams
+                        brush spline
                     pure $
                       StrokeWithOutlineRenderData
                         { strokeDataSpline    = spline
                         , strokeOutlineData   = ( outline, fitPts, cusps )
-                        , strokeBrushFunction = fun @Double ( brushFn @2 @() proxy# id )
-                                              . embedUsedParams
-                                              . toUsedParams
+                        , strokeBrushFunction =
+                            \ params ->
+                            let brushParams = embedUsedParams $ toUsedParams params
+                                shape = fun @Double brushShape brushParams
+                                -- TODO: remove this logic which is duplicated
+                                -- from elsewhere. The type should make it
+                                -- impossible to forget to apply the rotation.
+                            in case mbRot of
+                                 Nothing -> shape
+                                 Just getθ ->
+                                   let θ = getθ brushParams
+                                       cosθ = cos θ
+                                       sinθ = sin θ
+                                   in fmap ( unT . rotate cosθ sinθ . T ) shape
                         }
         _ -> pure $
                 StrokeRenderData
diff --git a/src/metabrushes/MetaBrush/Asset/Brushes.hs b/src/metabrushes/MetaBrush/Asset/Brushes.hs
index 7f6b160..45139af 100644
--- a/src/metabrushes/MetaBrush/Asset/Brushes.hs
+++ b/src/metabrushes/MetaBrush/Asset/Brushes.hs
@@ -15,6 +15,8 @@ module MetaBrush.Asset.Brushes
 import Prelude
   hiding
     ( Num(..), Floating(..), (^), (/), fromInteger, fromRational )
+import Data.Coerce
+  ( coerce )
 import GHC.Exts
   ( fromString )
 
@@ -30,12 +32,13 @@ import qualified Data.HashMap.Strict as HashMap
 
 -- brush-strokes
 import Calligraphy.Brushes
+  ( circleBrush, ellipseBrush, tearDropBrush )
 import Math.Linear
 import Math.Ring
 
 -- MetaBrush
 import MetaBrush.Brush
-  ( Brush(..), SomeBrush(..), WithParams(..) )
+  ( NamedBrush(..), SomeBrush(..), WithParams(..) )
 import MetaBrush.Records
   ( Record(MkR) )
 
@@ -48,7 +51,7 @@ lookupBrush nm = HashMap.lookup nm brushes
 brushes :: HashMap Text SomeBrush
 brushes = HashMap.fromList
   [ ( nm, b )
-  | b@( SomeBrush ( BrushData { brushName = nm } ) )
+  | b@( SomeBrush ( NamedBrush { brushName = nm } ) )
     <- [ SomeBrush circle, SomeBrush ellipse, SomeBrush tearDrop ]
   ]
 
@@ -56,8 +59,8 @@ brushes = HashMap.fromList
 
 type CircleBrushFields = '[ "r" ]
 -- | A circular brush with the given radius.
-circle :: Brush CircleBrushFields
-circle = BrushData "circle" ( WithParams deflts circleBrush )
+circle :: NamedBrush CircleBrushFields
+circle = NamedBrush "circle" ( WithParams deflts $ coerce circleBrush )
   where
     deflts :: Record CircleBrushFields
     deflts = MkR ( ℝ1 1 )
@@ -66,8 +69,8 @@ circle = BrushData "circle" ( WithParams deflts circleBrush )
 type EllipseBrushFields = '[ "a", "b", "phi" ]
 -- | An elliptical brush with the given semi-major and semi-minor axes and
 -- angle of rotation.
-ellipse :: Brush EllipseBrushFields
-ellipse = BrushData "ellipse" ( WithParams deflts ellipseBrush )
+ellipse :: NamedBrush EllipseBrushFields
+ellipse = NamedBrush "ellipse" ( WithParams deflts $ coerce ellipseBrush )
   where
     deflts :: Record EllipseBrushFields
     deflts = MkR ( ℝ3 1 1 0 )
@@ -75,8 +78,8 @@ ellipse = BrushData "ellipse" ( WithParams deflts ellipseBrush )
 
 type TearDropBrushFields = '[ "w", "h", "phi" ]
 -- | A tear-drop shape with the given width, height and angle of rotation.
-tearDrop :: Brush TearDropBrushFields
-tearDrop = BrushData "tear-drop" ( WithParams deflts tearDropBrush )
+tearDrop :: NamedBrush TearDropBrushFields
+tearDrop = NamedBrush "tear-drop" ( WithParams deflts $ coerce tearDropBrush )
   where
     deflts :: Record TearDropBrushFields
     deflts = MkR ( ℝ3 1 2.25 0 )
diff --git a/src/metabrushes/MetaBrush/Brush.hs b/src/metabrushes/MetaBrush/Brush.hs
index d9618e9..2cffed5 100644
--- a/src/metabrushes/MetaBrush/Brush.hs
+++ b/src/metabrushes/MetaBrush/Brush.hs
@@ -6,7 +6,7 @@
 
 module MetaBrush.Brush
   ( WithParams(..)
-  , Brush(..), SomeBrush(..), BrushFunction
+  , NamedBrush(..), SomeBrush(..), BrushFunction
   , PointFields, provePointFields, duplicates
   )
   where
@@ -39,16 +39,18 @@ import Data.Text
 import qualified Data.Text as Text
   ( unpack )
 
--- MetaBrush
-import Math.Algebra.Dual
-  ( C )
+-- brush-strokes
+import Calligraphy.Brushes
+  ( Brush(..) )
 import Math.Bezier.Spline
   ( SplineType(Closed), Spline )
 import Math.Differentiable
-  ( DiffInterp, I )
+  ( DiffInterp )
 import Math.Interval
   ( 𝕀 )
 import Math.Linear
+
+-- MetaBrush
 import MetaBrush.Records
   ( KnownSymbols, Length, Record )
 import MetaBrush.Serialisable
@@ -56,18 +58,12 @@ import MetaBrush.Serialisable
 
 --------------------------------------------------------------------------------
 
--- | A differentiable function from a given record type,
--- with provided default values that can be overridden.
-type WithParams :: Type -> ( Type -> Type ) -> Type
-data WithParams params f =
+-- | A brush, with default parameter values.
+type WithParams :: Type -> Type
+data WithParams params =
   WithParams
     { defaultParams :: params
-    , withParams
-        :: forall {t} k (i :: t)
-        .  ( DiffInterp k i params )
-        => Proxy# i
-        -> ( forall a. a -> I i a )
-        -> C k ( I i params ) ( f ( I i ( ℝ 2 ) ) )
+    , withParams :: Brush params
     }
 
 --------------------------------------------------------------------------------
@@ -75,11 +71,11 @@ data WithParams params f =
 -- | A brush function: a function from a record of parameters to a closed spline.
 type BrushFunction :: [ Symbol ] -> Type
 type BrushFunction brushFields =
-  WithParams ( Record brushFields ) ( Spline Closed () )
+  WithParams ( Record brushFields )
 
-type Brush :: [ Symbol ] -> Type
-data Brush brushFields where
-  BrushData
+type NamedBrush :: [ Symbol ] -> Type
+data NamedBrush brushFields where
+  NamedBrush
     :: forall brushFields
     .  ( KnownSymbols brushFields, Typeable brushFields
        , Representable Double ( ℝ ( Length brushFields ) )
@@ -89,28 +85,28 @@ data Brush brushFields where
     => { brushName     :: !Text
        , brushFunction :: BrushFunction brushFields
        }
-    -> Brush brushFields
+    -> NamedBrush brushFields
 
 data SomeBrush where
   SomeBrush
-    :: { someBrush :: !( Brush brushFields ) }
+    :: { someBrush :: !( NamedBrush brushFields ) }
     -> SomeBrush
 
-instance Show ( Brush brushFields ) where
-  show ( BrushData { brushName } ) =
-    "BrushData\n\
+instance Show ( NamedBrush brushFields ) where
+  show ( NamedBrush { brushName } ) =
+    "NamedBrush\n\
     \  { brushName = "  <> Text.unpack brushName <> "\n\
     \  }"
 
-instance NFData ( Brush brushFields ) where
-  rnf ( BrushData { brushName } )
+instance NFData ( NamedBrush brushFields ) where
+  rnf ( NamedBrush { brushName } )
     = rnf brushName
-instance Eq ( Brush brushFields ) where
-  BrushData name1 _ == BrushData name2 _ = name1 == name2
-instance Ord ( Brush brushFields ) where
-  compare ( BrushData name1 _ ) ( BrushData name2 _ ) = compare name1 name2
-instance Hashable ( Brush brushFields ) where
-  hashWithSalt salt ( BrushData { brushName } ) =
+instance Eq ( NamedBrush brushFields ) where
+  NamedBrush name1 _ == NamedBrush name2 _ = name1 == name2
+instance Ord ( NamedBrush brushFields ) where
+  compare ( NamedBrush name1 _ ) ( NamedBrush name2 _ ) = compare name1 name2
+instance Hashable ( NamedBrush brushFields ) where
+  hashWithSalt salt ( NamedBrush { brushName } ) =
     hashWithSalt salt brushName
 
 type PointFields :: [ Symbol ] -> Constraint
@@ -153,6 +149,12 @@ provePointFields fieldNames k =
       , SomeSymbol @f2 _ <- someSymbolVal ( Text.unpack f2 )
       , SomeSymbol @f3 _ <- someSymbolVal ( Text.unpack f3 )
       -> k ( proxy# @'[ f1, f2, f3 ] )
+    [ f1, f2, f3, f4 ]
+      | SomeSymbol @f1 _ <- someSymbolVal ( Text.unpack f1 )
+      , SomeSymbol @f2 _ <- someSymbolVal ( Text.unpack f2 )
+      , SomeSymbol @f3 _ <- someSymbolVal ( Text.unpack f3 )
+      , SomeSymbol @f4 _ <- someSymbolVal ( Text.unpack f4 )
+      -> k ( proxy# @'[ f1, f2, f3, f4 ] )
     _ -> error $ "I haven't defined ℝ " ++ show ( length fieldNames )
 {-# INLINE provePointFields #-}
 
diff --git a/src/metabrushes/MetaBrush/Document.hs b/src/metabrushes/MetaBrush/Document.hs
index bcb76d5..a33b4c1 100644
--- a/src/metabrushes/MetaBrush/Document.hs
+++ b/src/metabrushes/MetaBrush/Document.hs
@@ -89,7 +89,7 @@ import Math.Module
 import Math.Linear
   ( ℝ(..), T(..) )
 import MetaBrush.Brush
-  ( Brush, PointFields )
+  ( NamedBrush, PointFields )
 import MetaBrush.Records
 import MetaBrush.Unique
   ( UniqueSupply, Unique, freshUnique )
@@ -174,7 +174,7 @@ data Stroke where
       { strokeName     :: !Text
       , strokeVisible  :: !Bool
       , strokeUnique   :: Unique
-      , strokeBrush    :: !( Maybe ( Brush brushFields ) )
+      , strokeBrush    :: !( Maybe ( NamedBrush brushFields ) )
       , strokeSpline   :: !( StrokeSpline clo pointParams )
       }
     -> Stroke
diff --git a/src/metabrushes/MetaBrush/Document/Draw.hs b/src/metabrushes/MetaBrush/Document/Draw.hs
index 6ae752d..1f577e8 100644
--- a/src/metabrushes/MetaBrush/Document/Draw.hs
+++ b/src/metabrushes/MetaBrush/Document/Draw.hs
@@ -63,7 +63,7 @@ import Math.Linear
 import MetaBrush.Assert
   ( assert )
 import MetaBrush.Brush
-  ( Brush(..), PointFields )
+  ( NamedBrush(..), PointFields )
 import MetaBrush.Document
   ( Document(..), DocumentContent(..)
   , Stroke(..), StrokeHierarchy(..), StrokeSpline
@@ -123,7 +123,7 @@ getOrCreateDrawAnchor uniqueSupply c doc@( Document { zoomFactor } ) =
             , strokeVisible = True
             , strokeUnique  = uniq
             , strokeSpline  = newSpline
-            , strokeBrush   = Nothing :: Maybe ( Brush ( '[] :: [ Symbol ] ) )
+            , strokeBrush   = Nothing :: Maybe ( NamedBrush ( '[] :: [ Symbol ] ) )
             }
         newDoc' :: Document
         newDoc'
@@ -234,7 +234,7 @@ withAnchorBrushData
      . ( pointParams ~ Record pointFields
        , PointFields pointFields
        )
-     => Maybe ( Brush brushFields )
+     => Maybe ( NamedBrush brushFields )
      -> pointParams
      -> r
      )
diff --git a/src/metabrushes/MetaBrush/Document/Serialise.hs b/src/metabrushes/MetaBrush/Document/Serialise.hs
index 750cac0..993a386 100644
--- a/src/metabrushes/MetaBrush/Document/Serialise.hs
+++ b/src/metabrushes/MetaBrush/Document/Serialise.hs
@@ -110,7 +110,7 @@ import Math.Linear
 import MetaBrush.Asset.Brushes
   ( lookupBrush )
 import MetaBrush.Brush
-  ( Brush(..), SomeBrush(..), provePointFields, duplicates )
+  ( NamedBrush(..), SomeBrush(..), provePointFields, duplicates )
 import MetaBrush.Document
   ( Document(..), DocumentContent(..), Guide(..)
   , Stroke(..), StrokeHierarchy(..), StrokeSpline
@@ -223,9 +223,9 @@ decodeFields = do
     dups  -> throwError ( JSON.ParseFailed $ "Duplicate field names in brush record type:\n" <> Text.unwords dups )
 
 
-encodeBrush :: Applicative f => JSON.Encoder f (Brush brushFields)
+encodeBrush :: Applicative f => JSON.Encoder f (NamedBrush brushFields)
 encodeBrush = JSON.Encoder.mapLikeObj
-  \ ( BrushData { brushName } ) ->
+  \ ( NamedBrush { brushName } ) ->
     JSON.Encoder.atKey' "name" JSON.Encoder.text brushName
 
 decodeBrush :: MonadIO m => JSON.Decoder m SomeBrush
@@ -277,14 +277,14 @@ decodeStroke uniqueSupply = do
       strokeSpline <- JSON.Decoder.atKey "spline" ( decodeSpline @Closed @( PointData ( Record pointFields ) ) decodePointData )
       pure $ case mbSomeBrush of
         Nothing ->
-          Stroke { strokeName, strokeVisible, strokeUnique, strokeSpline, strokeBrush = Nothing :: Maybe ( Brush '[] ) }
+          Stroke { strokeName, strokeVisible, strokeUnique, strokeSpline, strokeBrush = Nothing :: Maybe ( NamedBrush '[] ) }
         Just (SomeBrush brush) ->
           Stroke { strokeName, strokeVisible, strokeUnique, strokeSpline, strokeBrush = Just brush }
     else do
       strokeSpline <- JSON.Decoder.atKey "spline" ( decodeSpline @Open   @( PointData ( Record pointFields ) ) decodePointData )
       pure $ case mbSomeBrush of
         Nothing ->
-          Stroke { strokeName, strokeVisible, strokeUnique, strokeSpline, strokeBrush = Nothing :: Maybe ( Brush '[] ) }
+          Stroke { strokeName, strokeVisible, strokeUnique, strokeSpline, strokeBrush = Nothing :: Maybe ( NamedBrush '[] ) }
         Just (SomeBrush brush) ->
           Stroke { strokeName, strokeVisible, strokeUnique, strokeSpline, strokeBrush = Just brush }
 
diff --git a/src/metabrushes/MetaBrush/Records.hs b/src/metabrushes/MetaBrush/Records.hs
index d154ac2..07fe0a6 100644
--- a/src/metabrushes/MetaBrush/Records.hs
+++ b/src/metabrushes/MetaBrush/Records.hs
@@ -25,7 +25,7 @@ import GHC.TypeLits
   , SomeSymbol(..), someSymbolVal
   )
 import GHC.TypeNats
-  ( Nat, type (+) )
+  ( Nat, KnownNat, type (+) )
 import Unsafe.Coerce
   ( unsafeCoerce )
 
@@ -139,7 +139,8 @@ instance ( Torsor ( T ( 𝕀ℝ ( Length ks ) ) ) ( 𝕀ℝ ( Length ks ) )
             T ( 𝕀 ( MkR c_lo ) ( MkR c_hi ) )
 
 type instance RepDim ( Record ks ) = Length ks
-deriving newtype instance Representable r ( ℝ ( Length ks ) )
+deriving newtype instance ( KnownNat (Length ks)
+                          , Representable r ( ℝ ( Length ks ) ) )
                        => Representable r ( Record ks )
 
 type instance D k ( Record ks ) = D k ( ℝ ( Length ks ) )
@@ -238,16 +239,16 @@ doIntersection k =
       | ( _ :: Proxy# r1r2 ) <- proxy# @'[ ]
       , Refl <- ( unsafeCoerce Refl :: r1r2 :~: Intersect r1 r2 )
       -> k @r1r2 proxy#
-           VZ
-           VZ
+           ( Vec [] )
+           ( Vec [] )
 
     [   ( f1, r1_i1, r2_i1 ) ]
       | SomeSymbol @f1 _     <- someSymbolVal ( Text.unpack f1 )
       , ( _ :: Proxy# r1r2 ) <- proxy# @'[ f1 ]
       , Refl <- ( unsafeCoerce Refl :: r1r2 :~: Intersect r1 r2 )
       -> k @r1r2 proxy#
-           ( VS ( Fin r1_i1 ) VZ )
-           ( VS ( Fin r2_i1 ) VZ )
+           ( Vec [ Fin r1_i1 ] )
+           ( Vec [ Fin r2_i1 ] )
 
     [   ( f1, r1_i1, r2_i1 )
       , ( f2, r1_i2, r2_i2 ) ]
@@ -256,8 +257,8 @@ doIntersection k =
       , ( _ :: Proxy# r1r2 ) <- proxy# @'[ f1, f2 ]
       , Refl <- ( unsafeCoerce Refl :: r1r2 :~: Intersect r1 r2 )
       -> k @r1r2 proxy#
-           ( VS ( Fin r1_i1 ) $ VS ( Fin r1_i2 ) VZ )
-           ( VS ( Fin r2_i1 ) $ VS ( Fin r2_i2 ) VZ )
+           ( Vec [ Fin r1_i1, Fin r1_i2 ] )
+           ( Vec [ Fin r2_i1, Fin r2_i2 ] )
 
     [   ( f1, r1_i1, r2_i1 )
       , ( f2, r1_i2, r2_i2 )
@@ -268,8 +269,8 @@ doIntersection k =
       , ( _ :: Proxy# r1r2 ) <- proxy# @'[ f1, f2, f3 ]
       , Refl <- ( unsafeCoerce Refl :: r1r2 :~: Intersect r1 r2 )
       -> k @r1r2 proxy#
-           ( VS ( Fin r1_i1 ) $ VS ( Fin r1_i2 ) $ VS ( Fin r1_i3 ) VZ )
-           ( VS ( Fin r2_i1 ) $ VS ( Fin r2_i2 ) $ VS ( Fin r2_i3 ) VZ )
+           ( Vec [ Fin r1_i1, Fin r1_i2, Fin r1_i3 ] )
+           ( Vec [ Fin r2_i1, Fin r2_i2, Fin r2_i3 ] )
 
     [   ( f1, r1_i1, r2_i1 )
       , ( f2, r1_i2, r2_i2 )
@@ -282,8 +283,8 @@ doIntersection k =
       , ( _ :: Proxy# r1r2 ) <- proxy# @'[ f1, f2, f3, f4 ]
       , Refl <- ( unsafeCoerce Refl :: r1r2 :~: Intersect r1 r2 )
       -> k @r1r2 proxy#
-           ( VS ( Fin r1_i1 ) $ VS ( Fin r1_i2 ) $ VS ( Fin r1_i3 ) $ VS ( Fin r1_i4 ) VZ )
-           ( VS ( Fin r2_i1 ) $ VS ( Fin r2_i2 ) $ VS ( Fin r2_i3 ) $ VS ( Fin r2_i4 ) VZ )
+           ( Vec [ Fin r1_i1, Fin r1_i2, Fin r1_i3, Fin r1_i4 ] )
+           ( Vec [ Fin r2_i1, Fin r2_i2, Fin r2_i3, Fin r2_i4 ] )
 
     other -> error $ "Intersection not defined in dimension " ++ show ( length other )