some optimisations

MetaBrush.cabal

@@ -90,10 +90,12 @@ common common
         ViewPatterns
 
     ghc-options:
-        -O1
+        -O2
         -fexpose-all-unfoldings
-        -funfolding-use-threshold=16
+      --  -funfolding-use-threshold=1000
         -fspecialise-aggressively
+        -flate-dmd-anal
+        -fmax-worker-args=200
         -optc-O3
         -Wall
         -Wcompat
@@ -102,7 +104,6 @@ common common
         -fwarn-incomplete-uni-patterns
         -fwarn-missing-deriving-strategies
         -fno-warn-unticked-promoted-constructors
-        -rtsopts
 
     if flag(asserts)
       cpp-options:

src/app/MetaBrush/Render/Document.hs

@@ -57,7 +57,7 @@ import Control.Monad.Trans.State.Strict
 
 -- MetaBrush
 import Math.Algebra.Dual
-  ( fun )
+  ( D2𝔸1(..), fun )
 import qualified Math.Bezier.Cubic     as Cubic
   ( Bezier(..), fromQuadratic )
 import Math.Bezier.Cubic.Fit
@@ -79,6 +79,8 @@ import Math.Interval
   ( Extent(Point) )
 import Math.Linear
   ( ℝ(..), T(..) )
+import Math.Module
+  ( Module((*^)), normalise )
 import MetaBrush.Asset.Colours
   ( Colours, ColourRecord(..) )
 import MetaBrush.Brush
@@ -648,22 +650,28 @@ drawFitPoint _ zoom ( FitTangent { fitPoint = ℝ2 x y, fitTangent = V2 tx ty }
 
 drawCusp :: Colours -> Double -> Cusp -> Cairo.Render ()
 drawCusp _ zoom
-  ( Cusp { cuspPathCoords   = ℝ2 px py
+  ( Cusp { cuspPathCoords   = D21 { _D21_v  = ℝ2 px py
+                                  , _D21_dx = tgt }
          , cuspStrokeCoords = ℝ2 cx cy } ) = do
 
-    -- Draw a white circle on the path at the cusp point.
+    -- Draw a line perpendicular to the underlying path at the cusp.
+    let
+      !( V2 tx ty ) = ( 6 / zoom ) *^ normalise tgt
     Cairo.save
     Cairo.translate px py
-    Cairo.arc 0 0 ( 4 / zoom ) 0 ( 2 * pi )
+    Cairo.moveTo -ty  tx
-    Cairo.setSourceRGBA 1.0 1.0 1.0 1.0
+    Cairo.lineTo  ty -tx
+    --withRGBA path Cairo.setSourceRGBA
+    Cairo.setSourceRGBA 0 0 0 0.75
+    Cairo.setLineWidth ( 2 / zoom )
     Cairo.stroke
     Cairo.restore
 
-    -- Draw a black circle on the outline at the cusp point.
+    -- Draw a circle around the outline cusp point.
     Cairo.save
     Cairo.translate cx cy
     Cairo.arc 0 0 ( 4 / zoom ) 0 ( 2 * pi )
-    Cairo.setSourceRGBA 0 0 0 1.0
+    Cairo.setSourceRGBA 0 0 0 0.75
     Cairo.stroke
     Cairo.restore
 

src/metabrushes/MetaBrush/Asset/Brushes.hs

@@ -1,9 +1,12 @@
 {-# LANGUAGE AllowAmbiguousTypes #-}
 {-# LANGUAGE OverloadedStrings   #-}
 {-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE TemplateHaskell     #-}
 
-module MetaBrush.Asset.Brushes where
+module MetaBrush.Asset.Brushes
+  ( lookupBrush, brushes
+  ,  CircleBrushFields,  circle
+  , EllipseBrushFields, ellipse
+  ) where
 
 -- base
 import Prelude
@@ -30,7 +33,7 @@ import qualified Data.HashMap.Strict as HashMap
 import Math.Algebra.Dual
 import Math.Bezier.Spline
 import Math.Interval
-  ( type I )
+
 import Math.Linear
 import Math.Module
   ( Module((^+^), (*^)) )
@@ -65,6 +68,7 @@ circle = BrushData "circle" ( WithParams deflts circleBrush )
   where
     deflts :: Record CircleBrushFields
     deflts = MkR ( ℝ1 1 )
+{-# INLINE circle #-}
 
 type EllipseBrushFields = '[ "a", "b", "phi" ]
 ellipse :: Brush EllipseBrushFields
@@ -72,6 +76,7 @@ ellipse = BrushData "ellipse" ( WithParams deflts ellipseBrush )
   where
     deflts :: Record EllipseBrushFields
     deflts = MkR ( ℝ3 1 1 0 )
+{-# INLINE ellipse #-}
 
 --------------------------------------------------------------------------------
 -- Differentiable brushes.
@@ -91,6 +96,7 @@ circleSpline p = sequenceA $
       ]
     lastCrv =
         Bezier3To ( p  κ -1 ) ( p  1 -κ ) BackToStart ()
+{-# INLINE circleSpline #-}
 
 circleBrush :: forall i k irec
             . ( irec ~ I i ( Record CircleBrushFields )
@@ -120,6 +126,7 @@ circleBrush _ mkI =
     e_y = pure $ mkI $ ℝ2 0 1
 
     kon = konst @( I i Double ) @k @irec . mkI
+{-# INLINEABLE circleBrush #-}
 
 ellipseBrush :: forall i k irec
              . ( irec ~ I i ( Record EllipseBrushFields )
@@ -154,3 +161,4 @@ ellipseBrush _ mkI =
     e_y = pure $ mkI $ ℝ2 0 1
 
     kon = konst @( I i Double ) @k @irec . mkI
+{-# INLINEABLE ellipseBrush #-}

src/splines/Math/Algebra/Dual.hs

@@ -37,6 +37,7 @@ import GHC.TypeNats
 
 -- MetaBrush
 import Math.Algebra.Dual.Internal
+import Math.Interval.Internal
 import Math.Linear
 import Math.Module
 import Math.Monomial
@@ -258,29 +259,33 @@ instance Ring r => Ring ( D2𝔸1 r ) where
   !dr1 * !dr2 =
     let
       o :: r
-      o = fromInteger 0
+      !o = fromInteger 0
       p :: r -> r -> r
-      p = (+) @r
+      !p = (+) @r
       m :: r -> r -> r
-      m = (*) @r
+      !m = (*) @r
    in
        $$( prodRuleQ
              [|| o ||] [|| p ||] [|| m ||]
              [|| dr1 ||] [|| dr2 ||] )
+  {-# SPECIALISE instance Ring ( D2𝔸1 Double ) #-}
+  {-# SPECIALISE instance Ring ( D2𝔸1 ( 𝕀 Double ) ) #-}
 
 instance Ring r => Ring ( D3𝔸1 r ) where
   !dr1 * !dr2 =
     let
       o :: r
-      o = fromInteger 0
+      !o = fromInteger 0
       p :: r -> r -> r
-      p = (+) @r
+      !p = (+) @r
       m :: r -> r -> r
-      m = (*) @r
+      !m = (*) @r
    in
        $$( prodRuleQ
              [|| o ||] [|| p ||] [|| m ||]
              [|| dr1 ||] [|| dr2 ||] )
+  {-# SPECIALISE instance Ring ( D3𝔸1 Double ) #-}
+  {-# SPECIALISE instance Ring ( D3𝔸1 ( 𝕀 Double ) ) #-}
 
 --instance Ring r => Ring ( D1𝔸2 r ) where
 --  !dr1 * !dr2 =
@@ -300,29 +305,33 @@ instance Ring r => Ring ( D2𝔸2 r ) where
   !dr1 * !dr2 =
     let
       o :: r
-      o = fromInteger 0
+      !o = fromInteger 0
       p :: r -> r -> r
-      p = (+) @r
+      !p = (+) @r
       m :: r -> r -> r
-      m = (*) @r
+      !m = (*) @r
    in
        $$( prodRuleQ
              [|| o ||] [|| p ||] [|| m ||]
              [|| dr1 ||] [|| dr2 ||] )
+  {-# SPECIALISE instance Ring ( D2𝔸2 Double ) #-}
+  {-# SPECIALISE instance Ring ( D2𝔸2 ( 𝕀 Double ) ) #-}
 
 instance Ring r => Ring ( D3𝔸2 r ) where
   !dr1 * !dr2 =
     let
       o :: r
-      o = fromInteger 0
+      !o = fromInteger 0
       p :: r -> r -> r
-      p = (+) @r
+      !p = (+) @r
       m :: r -> r -> r
-      m = (*) @r
+      !m = (*) @r
    in
        $$( prodRuleQ
              [|| o ||] [|| p ||] [|| m ||]
              [|| dr1 ||] [|| dr2 ||] )
+  {-# SPECIALISE instance Ring ( D3𝔸2 Double ) #-}
+  {-# SPECIALISE instance Ring ( D3𝔸2 ( 𝕀 Double ) ) #-}
 
 --instance Ring r => Ring ( D1𝔸3 r ) where
 --  !dr1 * !dr2 =
@@ -342,29 +351,33 @@ instance Ring r => Ring ( D2𝔸3 r ) where
   !dr1 * !dr2 =
     let
       o :: r
-      o = fromInteger 0
+      !o = fromInteger 0
       p :: r -> r -> r
-      p = (+) @r
+      !p = (+) @r
       m :: r -> r -> r
-      m = (*) @r
+      !m = (*) @r
    in
        $$( prodRuleQ
              [|| o ||] [|| p ||] [|| m ||]
              [|| dr1 ||] [|| dr2 ||] )
+  {-# SPECIALISE instance Ring ( D2𝔸3 Double ) #-}
+  {-# SPECIALISE instance Ring ( D2𝔸3 ( 𝕀 Double ) ) #-}
 
 instance Ring r => Ring ( D3𝔸3 r ) where
   !dr1 * !dr2 =
     let
       o :: r
-      o = fromInteger 0
+      !o = fromInteger 0
       p :: r -> r -> r
-      p = (+) @r
+      !p = (+) @r
       m :: r -> r -> r
-      m = (*) @r
+      !m = (*) @r
    in
        $$( prodRuleQ
              [|| o ||] [|| p ||] [|| m ||]
              [|| dr1 ||] [|| dr2 ||] )
+  {-# SPECIALISE instance Ring ( D3𝔸3 Double ) #-}
+  {-# SPECIALISE instance Ring ( D3𝔸3 ( 𝕀 Double ) ) #-}
 
 --instance Ring r => Ring ( D1𝔸4 r ) where
 --  !dr1 * !dr2 =
@@ -384,29 +397,33 @@ instance Ring r => Ring ( D2𝔸4 r ) where
   !dr1 * !dr2 =
     let
       o :: r
-      o = fromInteger 0
+      !o = fromInteger 0
       p :: r -> r -> r
-      p = (+) @r
+      !p = (+) @r
       m :: r -> r -> r
-      m = (*) @r
+      !m = (*) @r
    in
        $$( prodRuleQ
              [|| o ||] [|| p ||] [|| m ||]
              [|| dr1 ||] [|| dr2 ||] )
+  {-# SPECIALISE instance Ring ( D2𝔸4 Double ) #-}
+  {-# SPECIALISE instance Ring ( D2𝔸4 ( 𝕀 Double ) ) #-}
 
 instance Ring r => Ring ( D3𝔸4 r ) where
   !dr1 * !dr2 =
     let
       o :: r
-      o = fromInteger 0
+      !o = fromInteger 0
       p :: r -> r -> r
-      p = (+) @r
+      !p = (+) @r
       m :: r -> r -> r
-      m = (*) @r
+      !m = (*) @r
    in
        $$( prodRuleQ
              [|| o ||] [|| p ||] [|| m ||]
              [|| dr1 ||] [|| dr2 ||] )
+  {-# SPECIALISE instance Ring ( D3𝔸4 Double ) #-}
+  {-# SPECIALISE instance Ring ( D3𝔸4 ( 𝕀 Double ) ) #-}
 
 --------------------------------------------------------------------------------
 -- Field & transcendental instances

src/splines/Math/Algebra/Dual/Internal.hs

@@ -21,6 +21,10 @@ import GHC.Generics
 import GHC.TypeNats
   ( KnownNat )
 
+-- deepseq
+import Control.DeepSeq
+  ( NFData )
+
 -- template-haskell
 import Language.Haskell.TH
   ( CodeQ )
@@ -50,6 +54,7 @@ import TH.Utils
 -- | \( \mathbb{Z} \).
 newtype D𝔸0 v = D0 { _D0_v :: v }
   deriving stock   ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving anyclass NFData
   deriving Applicative
     via Generically1 D𝔸0
 
@@ -59,6 +64,7 @@ data D1𝔸1 v =
       , _D11_dx :: !( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving anyclass NFData
   deriving Applicative
     via Generically1 D1𝔸1
 
@@ -69,6 +75,7 @@ data D2𝔸1 v =
       , _D21_dxdx :: !( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving anyclass NFData
   deriving Applicative
     via Generically1 D2𝔸1
 
@@ -89,6 +96,7 @@ data D1𝔸2 v =
       , _D12_dx, _D12_dy :: !( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving anyclass NFData
   deriving Applicative
     via Generically1 D1𝔸2
 
@@ -110,6 +118,7 @@ data D3𝔸2 v =
       , _D32_dxdxdx, _D32_dxdxdy, _D32_dxdydy, _D32_dydydy :: !( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving anyclass NFData
   deriving Applicative
     via Generically1 D3𝔸2
 
@@ -119,6 +128,7 @@ data D1𝔸3 v =
       , _D13_dx, _D13_dy, _D13_dz :: !( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving anyclass NFData
   deriving Applicative
     via Generically1 D1𝔸3
 
@@ -129,6 +139,7 @@ data D2𝔸3 v =
       , _D23_dxdx, _D23_dxdy, _D23_dydy, _D23_dxdz, _D23_dydz, _D23_dzdz :: !( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving anyclass NFData
   deriving Applicative
     via Generically1 D2𝔸3
 
@@ -141,6 +152,7 @@ data D3𝔸3 v =
          , _D33_dxdxdz, _D33_dxdydz, _D33_dxdzdz, _D33_dydydz, _D33_dydzdz, _D33_dzdzdz :: !( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving anyclass NFData
   deriving Applicative
     via Generically1 D3𝔸3
 
@@ -150,6 +162,7 @@ data D1𝔸4 v =
       , _D14_dx, _D14_dy, _D14_dz, _D14_dw :: !( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving anyclass NFData
   deriving Applicative
     via Generically1 D1𝔸4
 
@@ -161,6 +174,7 @@ data D2𝔸4 v =
          , _D24_dydz, _D24_dzdz, _D24_dxdw, _D24_dydw, _D24_dzdw, _D24_dwdw :: !( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving anyclass NFData
   deriving Applicative
     via Generically1 D2𝔸4
 
@@ -176,6 +190,7 @@ data D3𝔸4 v =
          _D34_dxdwdw, _D34_dydwdw, _D34_dzdwdw, _D34_dwdwdw :: !( T v )
       }
   deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
+  deriving anyclass NFData
   deriving Applicative
     via Generically1 D3𝔸4
 

src/splines/Math/Bezier/Cubic/Fit.hs

@@ -74,7 +74,7 @@ import Math.Epsilon
 import Math.Linear.Solve
   ( linearSolve )
 import Math.Module
-  ( Module((*^), (^-^))
+  ( Module((*^), (^-^)), lerp
   , Inner((^.^)), quadrance
   )
 import Math.Roots
@@ -82,6 +82,8 @@ import Math.Roots
 import Math.Linear
   ( Mat22(..), ℝ(..), T(..) )
 
+import Debug.Utils
+
 --------------------------------------------------------------------------------
 
 -- | Parameters to the curve fitting algorithm.
@@ -121,37 +123,40 @@ data FitPoint
 fitSpline
   :: FitParameters
   -> ( ℝ 1 -> ( ℝ 2, T ( ℝ 2 ) ) ) -- ^ curve \( t \mapsto ( C(t), C'(t) ) \) to fit
+  -> ( ℝ 1, ℝ 1 )                  -- ^ interval \( [t_min, t_max] \)
   -> ( SplinePts Open, Seq FitPoint )
-fitSpline ( FitParameters {..} ) = go 0
+fitSpline ( FitParameters {..} ) curveFn = go 0
   where
     dt :: Double
     dt = recip ( fromIntegral nbSegments )
     go
       :: Int
-      -> ( ℝ 1 -> ( ℝ 2, T ( ℝ 2 ) ) )
+      -> ( ℝ 1, ℝ 1 )
       -> ( SplinePts Open, Seq FitPoint )
-    go subdiv curveFn =
+    go subdiv (t_min, t_max) =
       let
         p, r :: ℝ 2
         tp, tr :: T ( ℝ 2 )
         qs :: [ ℝ 2 ]
-        (p, tp) = curveFn $ ℝ1 0
+        (p, tp) = curveFn $ T ( ℝ1  0.0001 ) • t_min
-        (r, tr) = curveFn $ ℝ1 1
+        (r, tr) = curveFn $ T ( ℝ1 -0.0001 ) • t_max
-        qs = [ fst $ curveFn ( ℝ1 $ dt * fromIntegral j ) | j <- [ 1 .. nbSegments - 1 ] ]
+        qs = [ fst $ curveFn ( lerp @( T ( ℝ 1 ) ) ( dt * fromIntegral j ) t_min t_max )
+             | j <- [ 1 .. nbSegments - 1 ] ]
       in
         case fitPiece dist_tol t_tol maxIters p tp qs r tr of
-          ( bez, Max ( Arg max_sq_error t_split ) )
+          ( bez, Max ( Arg max_sq_error t_split_0 ) )
             |  subdiv >= maxSubdiv
             || max_sq_error <= dist_tol ^ ( 2 :: Int )
-            -> ( openCubicBezierCurveSpline () bez, ( FitTangent p tp :<| Seq.fromList ( map FitPoint qs ) ) :|> FitTangent r tr )
+            -> trace ( unlines [ "fitSpline: piece is OK", "t_min = " ++ show t_min, "start = " ++ show p, "start tgt = " ++ show tp, "t_max = " ++ show t_max, "end = " ++ show r, "end tgt = " ++ show tr ] )
+             $ ( openCubicBezierCurveSpline () bez, ( FitTangent p tp :<| Seq.fromList ( map FitPoint qs ) ) :|> FitTangent r tr )
             | let
-                t_split_eff :: Double
+                t_split :: ℝ 1
-                t_split_eff = min ( 1 - dt ) $ max dt t_split
+                t_split = ℝ1 $ min ( 1 - dt ) $ max dt t_split_0
                 c1, c2 :: SplinePts Open
                 ps1, ps2 :: Seq FitPoint
                 ( ( c1, ps1 ), ( c2, ps2 ) )
-                  = ( go ( subdiv + 1 ) ( \ ( ℝ1 t ) -> curveFn $ ℝ1 $ t * t_split_eff )
+                  = ( go ( subdiv + 1 ) (t_min  , t_split)
-                    , go ( subdiv + 1 ) ( \ ( ℝ1 t ) -> curveFn $ ℝ1 $ t_split_eff + t * ( 1 - t_split_eff ) )
+                    , go ( subdiv + 1 ) (t_split, t_max  )
                     ) `Parallel.Strategy.using`
                        ( Parallel.Strategy.parTuple2 Parallel.Strategy.rdeepseq Parallel.Strategy.rdeepseq )
             -> ( c1 <> c2, ps1 <> ps2 )

src/splines/Math/Bezier/Stroke.hs

@@ -18,8 +18,6 @@ module Math.Bezier.Stroke
   where
 
 -- base
-import Prelude
-  hiding ( unzip )
 import Control.Arrow
   ( first, (***) )
 import Control.Applicative
@@ -37,11 +35,15 @@ import Data.Foldable
 import Data.Functor.Identity
   ( Identity(..) )
 import Data.List
-  ( nub, partition )
+  ( nub, partition, sort )
 import Data.List.NonEmpty
-  ( unzip )
+  ( NonEmpty )
+import qualified Data.List.NonEmpty as NE
+  ( cons, singleton, unzip )
 import Data.Maybe
   ( fromMaybe, isJust, listToMaybe, mapMaybe )
+import Data.Semigroup
+  ( sconcat )
 import GHC.Exts
   ( newMutVar#, runRW#
   , Proxy#, proxy#
@@ -211,8 +213,8 @@ data Cusp
   = Cusp
   { cuspParameters   :: !( ℝ 2 )
     -- ^ @(t,s)@ parameter values of the cusp
-  , cuspPathCoords   :: !( ℝ 2 )
+  , cuspPathCoords   :: !( D 2 ( ℝ 1 ) ( ℝ 2 ) )
-    -- ^ path point coordinates
+    -- ^ path point coordinates and tangent
   , cuspStrokeCoords :: !( ℝ 2 )
     -- ^ brush stroke point coordinates
   }
@@ -417,14 +419,19 @@ computeStrokeOutline fitParams ptParams toBrushParams brushFn spline@( Spline {
         -- No cached fit: compute the fit anew.
         Nothing -> do
           let
+            -- Split up the path at the cusps
             cusps :: [ Cusp ]
             cusps = outlineDefiniteCusps  fwdBwd
                  ++ outlinePotentialCusps fwdBwd
-            -- SLD TODO: here!!!
+            intervals :: NonEmpty ( ℝ 1, ℝ 1 )
+            intervals = splitInterval ( ℝ1 0, ℝ1 1 )
+                      $ sort [ ℝ1 t | Cusp { cuspParameters = ℝ2 t _ } <- cusps ]
+
             fwdData, bwdData :: ( SplinePts Open, Seq FitPoint )
             ( fwdData, bwdData ) =
-              ( fitSpline fitParams ( fst . outlineFn fwdBwd )
+              ( sconcat $ fmap ( fitSpline fitParams ( fst . outlineFn fwdBwd ) ) intervals
-              , fitSpline fitParams ( snd . outlineFn fwdBwd ) )
+              , sconcat $ fmap ( fitSpline fitParams ( snd . outlineFn fwdBwd ) ) intervals )
+                -- TODO: use foldMap1 once that's in base.
                 `Strats.using`
                   ( Strats.parTuple2 Strats.rdeepseq Strats.rdeepseq )
             outlineData :: OutlineData
@@ -564,7 +571,8 @@ outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv ->
     path :: C ( ExtentOrder 'Point ) ( ℝ 1 ) ( ℝ 2 )
     ( path, usedParams ) = pathAndUsedParams @Point id
 
-    curves :: ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum Point )
+    curves :: ℝ 1 -- t
+           -> Seq ( ℝ 1 {- s -} -> StrokeDatum Point )
     curves =
       brushStrokeData @Point @( ExtentOrder 'Point ) @brushParams
         path
@@ -574,7 +582,8 @@ outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv ->
         )
         ( brushFromParams @Point proxy# id )
 
-    curvesI :: 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 'Interval )
+    curvesI :: 𝕀ℝ 1 -- t
+            -> Seq ( 𝕀ℝ 1 {- s -} -> StrokeDatum 'Interval )
     curvesI = brushStrokeData @'Interval @( ExtentOrder 'Interval ) @brushParams
                pathI
                ( chainRule @( 𝕀 Double ) @( ExtentOrder 'Interval )
@@ -602,7 +611,7 @@ outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv ->
                              $ runD ( brushFromParams @Point proxy# id )
                              $ toBrushParams params_t
 
-    --( newtDunno, newtSols ) = intervalNewtonGS InverseMidJacobian 0.0001 curvesI
+    ( newtDunno, newtSols ) = intervalNewtonGS InverseMidJacobian 0.0001 curvesI
 
   in --trace
      --   ( unlines $
@@ -617,8 +626,8 @@ outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv ->
      --   ) $
      OutlineInfo
        { outlineFn = fwdBwd
-       , outlineDefiniteCusps  = [] -- map ( cuspCoords curves ) newtSols
+       , outlineDefiniteCusps  = map ( cuspCoords curves ) newtSols
-       , outlinePotentialCusps = [] -- map ( cuspCoords curves ) newtDunno
+       , outlinePotentialCusps = map ( cuspCoords curves ) newtDunno
        }
 {-# INLINEABLE outlineFunction #-}
 
@@ -640,6 +649,20 @@ lastTangent ( Spline { splineCurves = NoCurves } ) = Nothing
 lastTangent ( Spline { splineStart, splineCurves = ClosedCurves Empty          lst } ) = Just $ endTangent splineStart splineStart lst
 lastTangent ( Spline { splineStart, splineCurves = ClosedCurves ( _ :|> prev ) lst } ) = Just $ endTangent splineStart ( openCurveEnd prev ) lst
 
+-- | Split up the given interval at the list of points provided.
+--
+-- Assumes the list is sorted.
+splitInterval :: Ord a => ( a, a ) -> [ a ] -> NonEmpty ( a, a )
+splitInterval ( start, end ) []
+  = NE.singleton ( start, end )
+splitInterval ( start, end ) ( split : splits )
+  | split >= end
+  = NE.singleton ( start, end )
+  | split <= start
+  = splitInterval ( start, end ) splits
+  | otherwise
+  = NE.cons ( start, split ) $ splitInterval ( split, end ) splits
+
 --------------------------------------------------------------------------------
 
 -- | Compute the join at a point of discontinuity of the tangent vector direction (G1 discontinuity).
@@ -668,7 +691,7 @@ joinBetweenOffsets
     = let
         pcs, lastAndRest :: Maybe ( SplinePts Open )
         ( pcs, lastAndRest )
-          =  unzip
+          =  NE.unzip
           $ ( discardCurveData *** discardCurveData )
           .  splitSplineAt ( i2 - i1 )
          <$> dropCurves i1 openSpline
@@ -1049,6 +1072,7 @@ newtype ZipSeq a = ZipSeq { getZipSeq :: Seq a }
 instance Applicative ZipSeq where
   pure _ = error "only use Applicative ZipSeq with non-empty Traversable functors"
   liftA2 f ( ZipSeq xs ) ( ZipSeq ys ) = ZipSeq ( Seq.zipWith f xs ys )
+  {-# INLINE liftA2 #-}
 
 brushStrokeData :: forall i k brushParams arr
                 .  ( k ~ ExtentOrder i, CurveOrder k, arr ~ C k
@@ -1181,12 +1205,12 @@ cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 'Point ) )
            -> Cusp
 cuspCoords eqs ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), i, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
   | StrokeDatum
-    { dpath   = D21 { _D21_v = path  }
+    { dpath
     , dstroke = D22 { _D22_v = stroke } }
       <- ( eqs ( ℝ1 t_mid ) `Seq.index` i ) ( ℝ1 s_mid )
   = Cusp
   { cuspParameters   = ℝ2 t_mid s_mid
-  , cuspPathCoords   = path
+  , cuspPathCoords   = dpath
   , cuspStrokeCoords = stroke
   }
   where

src/splines/Math/Module.hs

@@ -9,7 +9,8 @@ module Math.Module
   ( Module(..), lerp
   , Inner(..), Cross(..)
   , Interpolatable
-  , norm, squaredNorm, quadrance, distance
+  , norm, squaredNorm, normalise
+  , quadrance, distance
   , proj, projC, closestPointOnSegment
   , strictlyParallel, convexCombination
 
@@ -73,6 +74,9 @@ norm = sqrt . squaredNorm
 squaredNorm :: forall m r. Inner r m => m -> r
 squaredNorm v = v ^.^ v
 
+normalise :: ( Floating r, Inner r m ) => m -> m
+normalise v = recip ( norm v ) *^ v
+
 -- | Squared distance between two points.
 quadrance :: forall v r p. ( Inner r v, Torsor v p ) => p -> p -> r
 quadrance p1 p2 = squaredNorm ( p1 --> p2 :: v )