fix incorrect brush join computation

* also add some curvature calculations (unused at the moment)

app/Main.hs

@@ -214,7 +214,7 @@ main = do
   maxHistorySizeTVar <- STM.newTVarIO @Int                                         1000
   fitParametersTVar  <- STM.newTVarIO @FitParameters
     ( FitParameters
-      { maxSubdiv  = 5
+      { maxSubdiv  = 10
       , nbSegments = 12
       , dist_tol   = 5e-3
       , t_tol      = 1e-4

src/app/MetaBrush/Render/Document.hs

@@ -440,7 +440,7 @@ drawStroke cols@( Colours {..} ) debug zoom strokeData = do
         False -> Cairo.fill
         True  -> do
           Cairo.fillPreserve
-          withRGBA brushCenter Cairo.setSourceRGBA
+          Cairo.setSourceRGBA 0 0 0 0.75
           Cairo.setLineWidth ( 2 / zoom )
           Cairo.stroke
           ( `evalStateT` 0 ) $ traverse_ ( drawFitPoint cols zoom ) fitPts
@@ -451,7 +451,7 @@ drawStroke cols@( Colours {..} ) debug zoom strokeData = do
         False -> Cairo.fill
         True  -> do
           Cairo.fillPreserve
-          withRGBA brushCenter Cairo.setSourceRGBA
+          Cairo.setSourceRGBA 0 0 0 0.75
           Cairo.setLineWidth ( 2 / zoom )
           Cairo.stroke
           ( `evalStateT` 0 ) $ traverse_ ( drawFitPoint cols zoom ) fitPts
@@ -488,7 +488,7 @@ drawFitPoint :: Colours -> Double -> FitPoint -> StateT Double Cairo.Render ()
 drawFitPoint ( Colours {..} ) zoom ( FitPoint { fitPoint = Point2D x y } ) = do
 
   hue <- get
-  put ( hue + 0.002 )
+  put ( hue + 0.01 )
   let
     r, g, b :: Double
     ( r, g, b ) = hsl2rgb hue 0.9 0.4
@@ -496,14 +496,14 @@ drawFitPoint ( Colours {..} ) zoom ( FitPoint { fitPoint = Point2D x y } ) = do
     Cairo.save
     Cairo.translate x y
     Cairo.arc 0 0 ( 2 / zoom ) 0 ( 2 * pi )
-    Cairo.setSourceRGBA r g b 1
+    Cairo.setSourceRGBA r g b 0.8
     Cairo.fill
     Cairo.restore
 
 drawFitPoint ( Colours {..} ) zoom ( FitTangent { fitPoint = Point2D x y, fitTangent = Vector2D tx ty } ) = do
 
   hue <- get
-  put ( hue + 0.002 )
+  put ( hue + 0.01 )
   let
     r, g, b :: Double
     ( r, g, b ) = hsl2rgb hue 0.9 0.4
@@ -513,7 +513,7 @@ drawFitPoint ( Colours {..} ) zoom ( FitTangent { fitPoint = Point2D x y, fitTan
     Cairo.moveTo 0 0
     Cairo.lineTo ( 0.05 * tx ) ( 0.05 * ty )
     Cairo.setLineWidth ( 1 / zoom )
-    Cairo.setSourceRGBA r g b 1
+    Cairo.setSourceRGBA r g b 0.8
     Cairo.stroke
     Cairo.arc 0 0 ( 2 / zoom ) 0 ( 2 * pi )
     Cairo.fill

src/lib/Math/Bezier/Cubic.hs

@@ -13,7 +13,8 @@
 
 module Math.Bezier.Cubic
   ( Bezier(..)
-  , bezier, bezier'
+  , bezier, bezier', bezier''
+  , curvature, squaredCurvature
   , subdivide
   , ddist, closestPoint
   )
@@ -55,6 +56,8 @@ import Data.Group.Generics
 -- MetaBrush
 import qualified Math.Bezier.Quadratic as Quadratic
   ( Bezier(Bezier), bezier )
+import Math.Epsilon
+  ( epsilon )
 import Math.Module
   ( Module (..)
   , lerp
@@ -97,6 +100,34 @@ bezier' ( Bezier {..} ) t
       ( lerp @v t ( p0 --> p1 ) ( p1 --> p2 ) )
       ( lerp @v t ( p1 --> p2 ) ( p2 --> p3 ) )
 
+-- | Second derivative of a cubic Bézier curve.
+bezier'' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
+bezier'' ( Bezier {..} ) t
+  = ( 6 *^ )
+  $ lerp @v t
+      ( p1 --> p0 ^+^ p1 --> p2 )
+      ( p2 --> p1 ^+^ p2 --> p3 )
+
+-- | Curvature of a quadratic Bézier curve.
+curvature :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> r -> r
+curvature bez t = sqrt $ squaredCurvature @v bez t
+
+-- | Square of curvature of a quadratic Bézier curve.
+squaredCurvature :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> r -> r
+squaredCurvature bez t
+  | sq_nm_g' < epsilon
+  = 1 / 0
+  | otherwise
+  = ( sq_nm_g' * squaredNorm @v g'' - ( g' ^.^ g'' ) ^ ( 2 :: Int ) )
+  / ( sq_nm_g' ^ ( 3 :: Int ) )
+  where
+    g', g'' :: v
+    g'  = bezier'  @v bez t
+    g'' = bezier'' @v bez t
+    sq_nm_g' :: r
+    sq_nm_g' = squaredNorm @v g'
+
+
 -- | Subdivide a cubic Bézier curve into two parts.
 subdivide :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> ( Bezier p, Bezier p )
 subdivide ( Bezier {..} ) t = ( Bezier p0 q1 q2 pt, Bezier pt r1 r2 p3 )

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

@@ -134,16 +134,13 @@ fitSpline ( FitParameters {..} ) = go 0
         qs = [ fst $ curve ( dt * fromIntegral j ) | j <- [ 1 .. nbSegments - 1 ] ]
       in
         case fitPiece dist_tol t_tol maxIters p tp qs r tr of
-          ( bez, Max ( Arg sq_d t_split ) )
+          ( bez, Max ( Arg max_sq_error t_split ) )
             |  subdiv >= maxSubdiv
-            || sq_d <= dist_tol ^ ( 2 :: Int )
+            || max_sq_error <= dist_tol ^ ( 2 :: Int )
             -> ( Seq.singleton bez, ( FitTangent p tp :<| Seq.fromList ( map FitPoint qs ) ) :|> FitTangent r tr )
             | let
                 t_split_eff :: Double
-                t_split_eff
-                  | t_split < 0.2 = 0.2
-                  | t_split > 0.8 = 0.8
-                  | otherwise     = t_split
+                t_split_eff = min ( 1 - dt ) $ max dt t_split
             -> go ( subdiv + 1 ) ( \ t -> curve $ t * t_split_eff )
             <> go ( subdiv + 1 ) ( \ t -> curve $ t_split_eff + t * ( 1 - t_split_eff ) )
 

src/lib/Math/Bezier/Quadratic.hs

@@ -13,7 +13,8 @@
 
 module Math.Bezier.Quadratic
   ( Bezier(..)
-  , bezier, bezier'
+  , bezier, bezier', bezier''
+  , curvature, squaredCurvature
   , subdivide
   , ddist, closestPoint
   )
@@ -53,6 +54,8 @@ import Data.Group.Generics
   ()
 
 -- MetaBrush
+import Math.Epsilon
+  ( epsilon )
 import Math.Module
   ( Module (..)
   , lerp
@@ -87,6 +90,29 @@ bezier ( Bezier {..} ) t = lerp @v t ( lerp @v t p0 p1 ) ( lerp @v t p1 p2 )
 bezier' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
 bezier' ( Bezier {..} ) t = 2 *^ lerp @v t ( p0 --> p1 ) ( p1 --> p2 )
 
+-- | Second derivative of a quadratic Bézier curve.
+bezier'' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> v
+bezier'' ( Bezier {..} ) = 2 *^ ( p1 --> p0 ^+^ p1 --> p2 )
+
+-- | Curvature of a quadratic Bézier curve.
+curvature :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> r -> r
+curvature bez t = sqrt $ squaredCurvature @v bez t
+
+-- | Square of curvature of a quadratic Bézier curve.
+squaredCurvature :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> r -> r
+squaredCurvature bez t
+  | sq_nm_g' < epsilon
+  = 1 / 0
+  | otherwise
+  = ( sq_nm_g' * squaredNorm @v g'' - ( g' ^.^ g'' ) ^ ( 2 :: Int ) )
+  / ( sq_nm_g' ^ ( 3 :: Int ) )
+  where
+    g', g'' :: v
+    g'  = bezier'  @v bez t
+    g'' = bezier'' @v bez
+    sq_nm_g' :: r
+    sq_nm_g' = squaredNorm @v g'
+
 -- | Subdivide a quadratic Bézier curve into two parts.
 subdivide :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> ( Bezier p, Bezier p )
 subdivide ( Bezier {..} ) t = ( Bezier p0 q1 pt, Bezier pt r1 p2 )

src/lib/Math/Bezier/Stroke.hs

@@ -205,14 +205,14 @@ stroke params allPts@( spt0 :<| spt1 :<| spts )
       -- Open curve.
       | otherwise
       = pure
-        ( ( endOffset   • joinWithBrush ( withTangent           tgt_end     brush_end   ) ( withTangent ( (-1) *^ tgt_end ) brush_end   ) brush_end, Empty )
+        ( ( endOffset • joinWithBrush ( withTangent tgt_end brush_end ) ( withTangent ( (-1) *^ tgt_end ) brush_end ) brush_end, Empty )
         , ( Empty, Empty ) -- handled separately: see 'startingCap' below
         )
 
     -- Final cap for an open curve. Handled separately for correct stroke order.
     startingCap :: Seq ( StrokePoint () ) 
     startingCap
-      =   startOffset • joinWithBrush ( withTangent ( (-1) *^ tgt_start ) brush_start ) ( withTangent           tgt_start brush_start ) brush_start
+      = startOffset • joinWithBrush ( withTangent ( (-1) *^ tgt_start ) brush_start ) ( withTangent tgt_start brush_start ) brush_start
 
     go :: StrokePoint d -> Seq ( StrokePoint d ) -> Par ( ( Seq ( StrokePoint () ), Seq FitPoint ), ( Seq ( StrokePoint () ), Seq FitPoint ) )
     go _ Empty = pure mempty
@@ -466,7 +466,7 @@ joinWithBrush
         start, middle, end :: Seq ( StrokePoint d )
         ( ( middle, end ), start ) = first ( Seq.splitAt i2 ) $ Seq.splitAt i1 pts
       in
-        snd ( splitFirstPiece t1 start ) <> removePointData middle <> fst ( splitFirstPiece t2 end )
+        snd ( splitFirstPiece t1 start ) <> dropFirstPiece start <> removePointData middle <> fst ( splitFirstPiece t2 end )
   where
     t1, t2 :: Double
     t1 = fromMaybe 0.5 mb_t1