compute correct open stroke caps

src/lib/Math/Bezier/Stroke.hs

@@ -22,7 +22,6 @@ module Math.Bezier.Stroke
   , CachedStroke(..), discardCache, invalidateCache
   , computeStrokeOutline, joinWithBrush
   , withTangent
-  , between, parallel
   )
   where
 
@@ -32,7 +31,7 @@ import Prelude
 import Control.Arrow
   ( first, (***) )
 import Control.Monad
-  ( guard, unless )
+  ( unless )
 import Control.Monad.ST
   ( RealWorld, ST )
 import Data.Bifunctor
@@ -42,7 +41,7 @@ import Data.Foldable
 import Data.List.NonEmpty
   ( unzip )
 import Data.Maybe
-  ( fromMaybe, mapMaybe )
+  ( fromMaybe, isJust, mapMaybe )
 import GHC.Exts
   ( newMutVar#, runRW# )
 import GHC.STRef
@@ -103,17 +102,20 @@ import Math.Bezier.Spline
   , KnownSplineType
     ( bifoldSpline, ibifoldSpline )
   , Spline(..), SplinePts, Curves(..), Curve(..)
-  , openCurveEnd, splitSplineAt, dropCurves
+  , openCurveStart, openCurveEnd, splitSplineAt, dropCurves
   )
 import qualified Math.Bezier.Quadratic as Quadratic
 import Math.Epsilon
   ( epsilon )
 import Math.Module
-  ( Module((^-^), (*^)), Inner((^.^))
+  ( Module((*^)), Inner((^.^))
   , lerp, squaredNorm, cross
+  , convexCombination, strictlyParallel
   )
 import Math.Orientation
-  ( Orientation(..), convexOrientation, splineOrientation )
+  ( Orientation(..), convexOrientation, splineOrientation
+  , between
+  )
 import Math.Roots
   ( solveQuadratic )
 import Math.Vector2D
@@ -152,7 +154,7 @@ instance Monoid OutlineData where
 
 newtype CachedStroke s = CachedStroke { cachedStrokeRef :: STRef s ( Maybe OutlineData ) }
 instance Show ( CachedStroke s ) where
-  show _ = "CachedStroke..."
+  show _ = "<<CachedStroke>>"
 instance NFData ( CachedStroke s ) where
   rnf _ = ()
 
@@ -198,7 +200,8 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart =
   SOpen
     | ( firstOutlineFwd, firstOutlineBwd ) :<| _ <- outlineFns
     , _ :|> ( lastOutlineFwd, lastOutlineBwd )   <- outlineFns
-    , _ :|> lastCurve <- openCurves $ splineCurves spline
+    , firstCurve :<| _ <- openCurves $ splineCurves spline
+    , prevCurves :|> lastCurve <- openCurves $ splineCurves spline
     , let
         endPt :: ptData
         endPt = openCurveEnd lastCurve
@@ -210,9 +213,26 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart =
         startBrush, endBrush :: SplinePts Closed
         startBrush = brushShape spt0
         endBrush = brushShape endPt
+
+        -- Computation of which brush segment to use for the end caps.
+        startTgt, endTgt :: Vector2D Double
+        startTgt = coords spt0 --> coords ( openCurveStart firstCurve )
+        endTgt = case prevCurves of
+          Empty          -> endTangent spt0 spt0                      lastCurve
+          _ :|> lastPrev -> endTangent spt0 ( openCurveEnd lastPrev ) lastCurve
+        startTestTgt, endTestTgt :: Vector2D Double
+        ( startTestTgt, endTestTgt ) =
+          case brushOrientation of
+            CCW -> ( Vector2D   sty (-stx), Vector2D   ety (-etx) )
+            CW  -> ( Vector2D (-sty)  stx , Vector2D (-ety)  etx  )
+          where
+            stx, sty, etx, ety :: Double
+            Vector2D stx sty = startTgt
+            Vector2D etx ety = endTgt
+        
         startCap, endCap :: SplinePts Open
         startCap
-          | brushOrientation == convexOrientation [ startTgtBwd, startTgtFwd ]
+          | isJust $ between brushOrientation startTgtBwd startTgtFwd startTestTgt
           = fmap ( MkVector2D ( coords spt0 ) • )
           $ joinWithBrush ( withTangent startTgtBwd startBrush ) ( withTangent startTgtFwd startBrush ) startBrush
           | otherwise
@@ -220,7 +240,7 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart =
           . reverseSpline
           $ joinWithBrush ( withTangent startTgtFwd startBrush ) ( withTangent startTgtBwd startBrush )startBrush
         endCap
-          | brushOrientation == convexOrientation [ endTgtFwd, endTgtBwd ]
+          | not . isJust $ between brushOrientation endTgtBwd endTgtFwd endTestTgt
           = fmap ( MkVector2D ( coords endPt ) • )
           . reverseSpline
           $ joinWithBrush ( withTangent endTgtBwd endBrush ) ( withTangent endTgtFwd endBrush ) endBrush
@@ -263,7 +283,7 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart =
   where
 
     brushOrientation :: Orientation
-    brushOrientation = splineOrientation ( brushShape spt0 )
+    brushOrientation = splineOrientation @Double ( brushShape spt0 )
 
     outlineFns
       :: Seq
@@ -349,7 +369,7 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart =
         brush0 = brushShape sp0
         fwdJoin, bwdJoin :: SplinePts Open
         fwdJoin
-          | tgtFwd `parallel` prevTgtFwd
+          | tgtFwd `strictlyParallel` prevTgtFwd
           = Spline { splineStart = coords sp0, splineCurves = OpenCurves Empty }
           | brushOrientation == convexOrientation [ prevTgtFwd, tgtFwd ]
           = fmap ( ptOffset • )
@@ -359,7 +379,7 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart =
           . reverseSpline
           $ joinWithBrush ( withTangent     tgtFwd brush0 ) ( withTangent prevTgtFwd brush0 ) brush0
         bwdJoin
-          | tgtBwd `parallel` prevTgtBwd
+          | tgtBwd `strictlyParallel` prevTgtBwd
           = Spline { splineStart = coords sp0, splineCurves = OpenCurves Empty }
           | brushOrientation == convexOrientation [ tgtBwd, prevTgtBwd ]
           = fmap ( ptOffset • )
@@ -433,33 +453,45 @@ outlineFunctions ptParams brushFn sp0 crv =
         -> ( brush3, Cubic.bezier @( Vector2D Double ) bez, Cubic.bezier' bez )
     fwd, bwd :: Double -> ( Point2D Double, Vector2D Double )
     fwd t
-      = ( off t
+      = ( offset ( withTangent ( fwd' t ) ( brush t ) ) • f t
-        , if squaredNorm offTgt < epsilon then f' t else offTgt
+        , fwd' t
         )
       where
         off :: Double -> Point2D Double
         off u = offset ( withTangent ( f' u ) ( brush u ) ) • f u
-        offTgt :: Vector2D Double
+        offTgt :: Double -> Vector2D Double
-        offTgt
+        offTgt u
-          | t < 0.5
+          | u < 0.5
-          = 1e9 *^ ( off t --> off (t + 1e-9) )
+          = 1e9 *^ ( off u --> off (u + 1e-9) )
           | otherwise
-          = 1e9 *^ ( off (t - 1e-9) --> off t )
+          = 1e9 *^ ( off (u - 1e-9) --> off u )
+        fwd' :: Double -> Vector2D Double
+        fwd' u
+          | squaredNorm ( offTgt u ) < epsilon
+          = f' u
+          | otherwise
+          = offTgt u
     bwd t
-      = ( off s
+      = ( offset ( withTangent ( (-1) *^ bwd' s ) ( brush s ) ) • f s
-        , if squaredNorm offTgt < epsilon then (-1) *^ f' s else offTgt
+        , bwd' s
         )
       where
         s :: Double
         s = 1 - t
         off :: Double -> Point2D Double
         off u = offset ( withTangent ( (-1) *^ f' u ) ( brush u ) ) • f u
-        offTgt :: Vector2D Double
+        offTgt :: Double -> Vector2D Double
-        offTgt
+        offTgt u
-          | s < 0.5
+          | u < 0.5
-          = 1e9 *^ ( off s --> off (s + 1e-9) )
+          = 1e9 *^ ( off u --> off (u + 1e-9) )
           | otherwise
-          = 1e9 *^ ( off (s - 1e-9) --> off s )
+          = 1e9 *^ ( off (u - 1e-9) --> off u )
+        bwd' :: Double -> Vector2D Double
+        bwd' u
+          | squaredNorm ( offTgt u ) < epsilon
+          = (-1) *^ f' u
+          | otherwise
+          = offTgt u
   in ( fwd, bwd )
 
 -----------------------------------
@@ -647,6 +679,8 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) = case lastTangent spli
             \tangent vector: " <> show tgt_wanted <> "\n\
             \spline: " <> show spline <> "\n"
   where
+    ori :: Orientation
+    ori = splineOrientation @Double spline
     go :: Int -> ptData -> Curve Open crvData ptData -> StateT ( Vector2D Double ) ( Except Offset ) ()
     go i cp cseg = do
       tgt_prev <- get
@@ -659,8 +693,8 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) = case lastTangent spli
         tgt_start = startTangent splineStart cp cseg
         tgt_end   = endTangent   splineStart cp cseg
       -- Handle corner.
-      unless ( parallel tgt_prev tgt_start ) do
+      unless ( tgt_prev `strictlyParallel` tgt_start ) do
-        for_ ( between tgt_wanted tgt_prev tgt_start ) \ _ ->
+        for_ ( between ori tgt_prev tgt_start tgt_wanted ) \ _ ->
           lift . throwE $
             Offset
               { offsetIndex     = i
@@ -673,7 +707,7 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) = case lastTangent spli
 
     handleSegment :: Int -> Point2D Double -> Curve Open crvData ( Point2D Double ) -> Vector2D Double -> Except Offset ()
     handleSegment i p0 ( LineTo ( NextPoint p1 ) _ ) tgt0
-      | parallel tgt_wanted tgt0
+      | tgt_wanted `strictlyParallel` tgt0
       , let
           offset :: Vector2D Double
           offset = MkVector2D $ lerp @( Vector2D Double ) 0.5 p0 p1
@@ -684,7 +718,7 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) = case lastTangent spli
       let
         tgt1 :: Vector2D Double
         tgt1 = p1 --> p2
-      in for_ ( between tgt_wanted tgt0 tgt1 ) \ t ->
+      in for_ ( convexCombination tgt0 tgt1 tgt_wanted ) \ t ->
           throwE $
             Offset
               { offsetIndex     = i
@@ -715,7 +749,7 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) = case lastTangent spli
           mbParam =
             case mapMaybe correctTangentParam $ solveQuadratic c01 ( 2 * ( c12 - c01 ) ) ( c01 + c23 - 2 * c12 ) of
               ( t : _ ) -> Just t
-              _         -> between tgt_wanted tgt0 tgt2 -- fallback in case we couldn't solve the quadratic for some reason
+              _         -> between ori tgt0 tgt2 tgt_wanted -- fallback in case we couldn't solve the quadratic for some reason
         in for_ mbParam \ t ->
             throwE $
               Offset
@@ -723,41 +757,3 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) = case lastTangent spli
                 , offsetParameter = Just t
                 , offset          = MkVector2D $ Cubic.bezier @( Vector2D Double ) bez t
                 }
-
--- | Finds whether the query vector @ u @ lies between the two provided vectors @ v0 @, @ v1 @.
---
--- If so, returns @ t @ in @ [ 0, 1 ] @ such that @ ( 1 - t ) v0 + t v1 @ is a positive multiple of @ u @.
-between
-  :: Vector2D Double -- ^ query vector
-  -> Vector2D Double -- ^ first vector
-  -> Vector2D Double -- ^ second vector
-  -> Maybe Double
-between u v0 v1
-  | abs c10 < epsilon
-  = if parallel u v0
-    then Just 0
-    else if parallel u v1
-      then Just 1
-      else Nothing
-  | otherwise
-  = do
-      let
-        t :: Double
-        t = c0 / c10
-      guard ( t > - epsilon && t < 1 + epsilon )
-      guard ( epsilon < u ^.^ ( lerp @( Vector2D Double ) t v0 v1 ) )
-      Just $ min 1 ( max 0 t )
-
-  where
-    c0, c10 :: Double
-    c0  = v0 `cross` u
-    c10 = ( v0 ^-^ v1 ) `cross` u
-
--- | Compute whether two vectors point in the same direction,
--- that is, whether each vector is a (strictly) positive multiple of the other.
---
--- Returns @False@ if either of the vectors is zero.
-parallel :: Vector2D Double -> Vector2D Double -> Bool
-parallel u v
-  =  abs ( u `cross` v ) < epsilon -- vectors are collinear
-  &&       u   ^.^   v   > epsilon -- vectors point in the same direction (parallel and not anti-parallel)

src/lib/Math/Module.hs

@@ -12,12 +12,15 @@ module Math.Module
   , squaredNorm, quadrance, distance
   , proj, projC, closestPointOnSegment
   , cross
+  , strictlyParallel, convexCombination
   )
   where
 
 -- base
 import Control.Applicative
   ( liftA2 )
+import Control.Monad
+  ( guard )
 import Data.Monoid
   ( Ap(..), Sum(..) )
 
@@ -34,6 +37,8 @@ import Data.Group
   ( invert )
 
 -- MetaBrush
+import Math.Epsilon
+  ( epsilon )
 import Math.Vector2D
   ( Vector2D(..), Segment(..) )
 
@@ -107,6 +112,7 @@ closestPointOnSegment c ( Segment p0 p1 )
     t :: r
     t = projC ( p0 --> c ) v01
 
+--------------------------------------------------------------------------------
 
 instance Num a => Module a ( Sum a ) where
 
@@ -136,6 +142,47 @@ instance Num a => Inner a ( Vector2D a ) where
   ( Vector2D x1 y1 ) ^.^ ( Vector2D x2 y2 )
     = x1 * x2 + y1 * y2
 
+-- | Cross-product of two 2D vectors.
 cross :: Num a => Vector2D a -> Vector2D a -> a
 cross ( Vector2D x1 y1 ) ( Vector2D x2 y2 )
-  = x1 * y2 - x2 * y1
+  = x1 * y2 - x2 * y1
+
+-- | Compute whether two vectors point in the same direction,
+-- that is, whether each vector is a (strictly) positive multiple of the other.
+--
+-- Returns @False@ if either of the vectors is zero.
+strictlyParallel :: RealFloat r => Vector2D r -> Vector2D r -> Bool
+strictlyParallel u v
+  =  abs ( u `cross` v ) < epsilon -- vectors are collinear
+  &&       u   ^.^   v   > epsilon -- vectors point in the same direction (parallel and not anti-parallel)
+
+-- | Finds whether the query vector @ u @ is a convex combination of the two provided vectors @ v0 @, @ v1 @.
+--
+-- If so, returns @ t @ in @ [ 0, 1 ] @ such that @ ( 1 - t ) v0 + t v1 @ is a positive multiple of @ u @.
+convexCombination
+  :: forall r
+  .  RealFloat r
+  => Vector2D r -- ^ first vector
+  -> Vector2D r -- ^ second vector
+  -> Vector2D r -- ^ query vector
+  -> Maybe r
+convexCombination v0 v1 u
+  | abs c10 < epsilon
+  = if strictlyParallel u v0
+    then Just 0
+    else if strictlyParallel u v1
+      then Just 1
+      else Nothing
+  | otherwise
+  = do
+      let
+        t :: r
+        t = c0 / c10
+      guard ( t > - epsilon && t < 1 + epsilon )
+      guard ( epsilon < u ^.^ ( lerp @( Vector2D r ) t v0 v1 ) )
+      Just $ min 1 ( max 0 t )
+
+  where
+    c0, c10 :: r
+    c0  = v0 `cross` u
+    c10 = ( v0 ^-^ v1 ) `cross` u

src/lib/Math/Orientation.hs

@@ -1,5 +1,7 @@
+{-# LANGUAGE AllowAmbiguousTypes #-}
 {-# LANGUAGE DataKinds           #-}
 {-# LANGUAGE DerivingStrategies  #-}
+{-# LANGUAGE FlexibleContexts    #-}
 {-# LANGUAGE GADTs               #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeApplications    #-}
@@ -7,9 +9,16 @@
 module Math.Orientation
   ( Orientation(..), reverseOrientation
   , convexOrientation, splineOrientation, splineTangents
+  , between
   )
   where
 
+-- base
+import Control.Monad
+  ( guard )
+import Data.Fixed
+  ( mod' )
+
 -- acts
 import Data.Act
   ( Torsor((-->)) )
@@ -18,6 +27,12 @@ import Data.Act
 import Data.Sequence
   ( Seq(..) )
 
+-- generic-lens
+import Data.Generics.Product.Typed
+  ( HasType(typed) )
+import Data.GenericLens.Internal
+  ( view )
+
 -- MetaBrush
 import Math.Epsilon
   ( nearZero )
@@ -29,7 +44,7 @@ import Math.Bezier.Spline
   , ssplineType
   )
 import Math.Vector2D
-  ( Point2D, Vector2D )
+  ( Point2D, Vector2D(..) )
 
 --------------------------------------------------------------------------------
 
@@ -58,34 +73,78 @@ convexOrientation _ = CCW -- default
 
 -- | Compute the orientation of a spline, assuming tangent vectors have a monotone angle.
 splineOrientation
-  :: ( KnownSplineType clo, RealFloat r )
+  :: forall r clo crvData ptData
-  => Spline clo crvData ( Point2D r )
+  .  ( KnownSplineType clo, RealFloat r, HasType ( Point2D r ) ptData )
+  => Spline clo crvData ptData
   -> Orientation
-splineOrientation = convexOrientation . splineTangents
+splineOrientation = convexOrientation . splineTangents @r
 
 -- | Compute the sequence of tangent vectors given by the control points of a Bézier spline.
 splineTangents
-  :: forall clo crvData r
+  :: forall r clo crvData ptData
-  .  ( Num r, KnownSplineType clo )
+  .  ( Num r, KnownSplineType clo, HasType ( Point2D r ) ptData )
-  => Spline clo crvData ( Point2D r )
+  => Spline clo crvData ptData
   -> [ Vector2D r ]
-splineTangents spline@( Spline { splineStart = p0, splineCurves = curves } ) =
+splineTangents spline@( Spline { splineStart = sp0, splineCurves = curves } )
-  case ssplineType @clo of
+  | let
-    SOpen
+      p0 :: Point2D r
-      | OpenCurves { openCurves = cs } <- curves
+      p0 = view typed sp0
-      -> go p0 cs
+  = case ssplineType @clo of
-    SClosed
+      SOpen
-      | OpenCurves cs@( c :<| _ ) <- splineCurves $ adjustSplineType @Open spline
+        | OpenCurves { openCurves = cs } <- curves
-      -> go p0 ( cs :|> c )
+        -> go p0 cs
-    _ -> []
+      SClosed
+        | OpenCurves cs@( c :<| _ ) <- splineCurves $ adjustSplineType @Open spline
+        -> go p0 ( cs :|> c )
+      _ -> []
   where
-    go :: Point2D r -> Seq ( Curve Open crvData ( Point2D r ) ) -> [ Vector2D r ]
+    go :: Point2D r -> Seq ( Curve Open crvData ptData ) -> [ Vector2D r ]
     go _ Empty = []
     go p ( crv :<| crvs ) =
       case crv of
-        LineTo { curveEnd = NextPoint q } ->
+        LineTo { curveEnd = NextPoint sq }
-          ( p --> q ) : go q crvs
+          | let
-        Bezier2To { controlPoint = cp, curveEnd = NextPoint q } ->
+              q :: Point2D r
-          ( p --> cp ) : ( cp --> q ) : go q crvs
+              q = view typed sq
-        Bezier3To { controlPoint1 = cp1, controlPoint2 = cp2, curveEnd = NextPoint q } ->
+          -> ( p --> q ) : go q crvs
-          ( p --> cp1 ) : ( cp2 --> q ) : go q crvs
+        Bezier2To { controlPoint = scp, curveEnd = NextPoint sq }
+          | let
+              cp, q :: Point2D r
+              cp = view typed scp
+              q  = view typed sq
+          -> ( p --> cp ) : ( cp --> q ) : go q crvs
+        Bezier3To { controlPoint1 = scp1, controlPoint2 = scp2, curveEnd = NextPoint sq }
+          | let
+              cp1, cp2, q :: Point2D r
+              cp1 = view typed scp1
+              cp2 = view typed scp2
+              q   = view typed sq
+          -> ( p --> cp1 ) : ( cp2 --> q ) : go q crvs
+
+-- | Checks whether a 2D vector lies "in between" two other vectors according to a given orientation,
+-- i.e. whether the angle of the query vector lies in between the angles of the start and end vectors.
+--
+-- Returns the proportion of the angle the vector is in between, or @Nothing@ if the query vector
+-- is not in between.
+--
+-- >>> between CCW ( Vector2D 1 0 ) ( Vector2D (-1) 1 ) ( Vector2D 1 1 )
+-- Just 0.3333333333333333
+between
+  :: forall r
+  .  RealFloat r
+  => Orientation
+  -> Vector2D r -- ^ start vector
+  -> Vector2D r -- ^ end vector
+  -> Vector2D r -- ^ query vector: is in between the start and end vectors w.r.t. the provided orientation?
+  -> Maybe r
+between CCW ( Vector2D x1 y1 ) ( Vector2D x2 y2 ) ( Vector2D a b ) =
+  let
+    τ, η, φ, θ :: r
+    τ = 2 * pi
+    η = atan2 y1 x1
+    φ = ( atan2 y2 x2 - η ) `mod'` τ
+    θ = ( atan2 b  a  - η ) `mod'` τ
+  in do
+    guard ( θ < φ )
+    pure  ( θ / φ )
+between CW v1 v2 u = ( 1 - ) <$> between CCW v2 v1 u