brush-strokes-0.1.0.0: src/lib/Math/Bezier/Stroke.hs
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE PartialTypeSignatures #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators #-}
{-# OPTIONS_GHC -fspecialise-aggressively #-}
module Math.Bezier.Stroke
( Offset(..), Cusp(..)
, CachedStroke(..), discardCache, newCache, invalidateCache
, computeStrokeOutline_specialise
, computeStrokeOutline, joinWithBrush
, withTangent
, RootSolvingAlgorithm(..) -- TODO: move this?
-- * Brush stroking
, envelopeEquation
, line, bezier2, bezier3
, brushStrokeData, pathAndUsedParams
-- ** Cusp finding
, findCusps, findCuspsIn
-- TODO: hack for switching between 2 and 3 dim formulations of cusp-finding
, N
)
where
-- base
import Prelude
hiding ( unzip )
import Control.Arrow
( first )
import Control.Monad
( unless )
import Control.Monad.ST
( RealWorld, ST )
import Data.Bifunctor
( Bifunctor(bimap) )
import Data.Coerce
( coerce )
import Data.Fixed
( divMod' )
import Data.Foldable
( for_ )
import Data.Functor
( (<&>), unzip )
import Data.Functor.Identity
( Identity(..) )
import Data.List
( sort )
import Data.List.NonEmpty
( NonEmpty )
import qualified Data.List.NonEmpty as NE
import Data.Maybe
( catMaybes, fromMaybe, isNothing, listToMaybe, mapMaybe )
import Data.Monoid
( First(..) )
import Data.Proxy
( Proxy(..))
import Data.Semigroup
( sconcat )
import Data.Type.Equality
( (:~:)(Refl) )
import GHC.Exts
( newMutVar#, runRW#, proxy# )
import GHC.STRef
( STRef(..), readSTRef, newSTRef, writeSTRef )
import GHC.Generics
( Generic, Generic1, Generically(..) )
import GHC.TypeNats
( Nat, KnownNat
, type (<=)
, sameNat, natVal'
)
-- acts
import Data.Act
( Act
( (•) )
, Torsor
( (-->) )
)
-- containers
import Data.IntMap.Strict
( IntMap )
import qualified Data.IntMap.Strict as IntMap
( fromList, mapWithKey, toList )
import Data.Sequence
( Seq(..) )
import qualified Data.Sequence as Seq
-- deepseq
import Control.DeepSeq
( NFData(..), NFData1, deepseq, force )
-- generic-lens
import Data.Generics.Product.Fields
( field' )
import Data.Generics.Product.Typed
( HasType(typed) )
import qualified Data.Generics.Internal.VL as Lens
( set, view, over )
-- parallel
import qualified Control.Parallel.Strategies as Strats
( rdeepseq, parTuple2, using )
-- transformers
import Control.Monad.Trans.Class
( lift )
import Control.Monad.Trans.Except
( Except, runExcept, throwE )
import Control.Monad.Trans.State.Strict as State
( StateT, runStateT, evalStateT, get, put )
import Control.Monad.Trans.Writer.CPS
( WriterT, execWriterT, runWriter, tell )
-- brush-strokes
import Calligraphy.Brushes
( Brush(..), Corner(..), Unrotated(..), BrushFn )
import Math.Algebra.Dual
import qualified Math.Bezier.Cubic as Cubic
import Math.Bezier.Cubic.Fit
( FwdBwd(..), FitPoint(..), OutlinePoint(..)
, FitParameters, fitSpline
)
import Math.Bezier.Spline
import qualified Math.Bezier.Quadratic as Quadratic
import Math.Bezier.Stroke.EnvelopeEquation
import Math.Differentiable
( I, IVness(..), SingIVness (..)
, Differentiable, DiffInterp
)
import Math.Epsilon
( epsilon )
import Math.Interval
import Math.Linear
import Math.Module
import Math.Ring
( Transcendental )
import qualified Math.Ring as Ring
import Math.Orientation
( Orientation(..), splineOrientation
, between
)
import Math.Roots
import Math.Root.Isolation
import Math.Root.Isolation.Utils
( boxMidpoint )
import Debug.Utils
--------------------------------------------------------------------------------
data Offset
= Offset
{ offsetIndex :: !Int
, offsetParameter :: !( Maybe Double )
, offset :: !( T ( ℝ 2 ) )
}
deriving stock ( Show, Generic )
deriving anyclass NFData
data TwoSided a
= TwoSided
{ fwd :: !a
, bwd :: !a
}
deriving stock ( Show, Generic, Generic1, Functor, Foldable, Traversable )
deriving anyclass ( NFData, NFData1 )
data OutlineData =
OutlineData
{ outline :: !( TwoSided ( SplinePts Open, Seq FitPoint ) )
, cusps :: ![ Cusp ] }
deriving stock Generic
deriving ( Semigroup, Monoid )
via Generically OutlineData
deriving anyclass NFData
instance Semigroup ( TwoSided ( SplinePts Open, Seq FitPoint ) ) where
TwoSided ( fwdSpline1, fwdPts1 ) ( bwdSpline1, bwdPts1 ) <> TwoSided ( fwdSpline2, fwdPts2 ) ( bwdSpline2, bwdPts2 ) =
TwoSided
( fwdSpline1 <> fwdSpline2, fwdPts1 <> fwdPts2 )
( bwdSpline2 <> bwdSpline1, bwdPts2 <> bwdPts1 )
instance Monoid ( TwoSided ( SplinePts Open, Seq FitPoint ) ) where
mempty = TwoSided empt empt
where
empt :: ( SplinePts Open, Seq FitPoint )
empt = ( Spline { splineStart = ℝ2 0 0, splineCurves = OpenCurves Empty }, Empty )
data CachedStroke s = NoCache | CachedStroke { cachedStrokeRef :: STRef s ( Maybe OutlineData ) }
instance Show ( CachedStroke s ) where
show _ = "<<CachedStroke>>"
instance NFData ( CachedStroke s ) where
rnf _ = ()
discardCache :: forall crvData s. HasType ( CachedStroke s ) crvData => crvData -> ST s ()
discardCache ( Lens.view ( typed @( CachedStroke s ) ) -> cache ) =
case cache of
NoCache -> return ()
CachedStroke { cachedStrokeRef } ->
writeSTRef cachedStrokeRef Nothing
{-# INLINE invalidateCache #-}
invalidateCache :: forall crvData. HasType ( CachedStroke RealWorld ) crvData => crvData -> crvData
invalidateCache = runRW# \ s ->
case newMutVar# Nothing s of
(# _, mutVar #) ->
Lens.set ( typed @( CachedStroke RealWorld ) )
( CachedStroke $ STRef mutVar )
newCache :: forall s. ST s ( CachedStroke s )
newCache = CachedStroke <$> newSTRef Nothing
coords :: forall ptData. HasType ( ℝ 2 ) ptData => ptData -> ℝ 2
coords = Lens.view typed
--------------------------------------------------------------------------------
-- | Forward and backward outlines.
type OutlineFn = ℝ 1 -> TwoSided OutlinePoint
data Cusp
= Cusp
{ cuspParameters :: !( ℝ 2 )
-- ^ @(t,s)@ parameter values of the cusp
, cuspPathCoords :: !( D 2 1 ( ℝ 2 ) )
-- ^ path point coordinates and tangent
, cuspStrokeCoords :: !( ℝ 2 )
-- ^ brush stroke point coordinates
, cornerCusp :: !Bool
}
deriving stock ( Generic, Show )
deriving anyclass NFData
data OutlineInfo =
OutlineInfo
{ outlineFn :: OutlineFn
, outlineDefiniteCusps, outlinePotentialCusps :: [ Cusp ] }
type N = 2
-- | A version of 'computeStrokeOutline' which dispatches to the correctly
-- specialised version at runtime.
computeStrokeOutline_specialise ::
forall ( clo :: SplineType ) ( nbUsedParams :: Nat ) ( nbBrushParams :: Nat ) crvData ptData s
. ( KnownSplineType clo
, HasType ( ℝ 2 ) ptData
, HasType ( CachedStroke s ) crvData
, NFData ptData, NFData crvData
-- Debugging.
, Show ptData, Show crvData
-- Constraints for runtime dispatch to specialisations
, KnownNat nbUsedParams, KnownNat nbBrushParams
, nbUsedParams <= 4, nbBrushParams <= 4, 1 <= nbBrushParams
, nbUsedParams <= nbBrushParams
)
=> RootSolvingAlgorithm
-> Maybe ( RootIsolationOptions 1 1, RootIsolationOptions 1 2, RootIsolationOptions 2 3 )
-> FitParameters
-> ( ptData -> ℝ nbUsedParams )
-> ( ℝ nbUsedParams -> ℝ nbBrushParams )
-> Brush nbBrushParams
-> Spline clo crvData ptData
-> ST s
( Either ( SplinePts Closed ) ( SplinePts Closed, SplinePts Closed )
, Seq FitPoint
, [ Cusp ]
)
computeStrokeOutline_specialise rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| Just Refl <- sameNat @nbUsedParams @0 Proxy Proxy
, Just Refl <- sameNat @nbBrushParams @1 Proxy Proxy
= computeStrokeOutline @0 @1 rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| Just Refl <- sameNat @nbUsedParams @1 Proxy Proxy
, Just Refl <- sameNat @nbBrushParams @1 Proxy Proxy
= computeStrokeOutline @1 @1 rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| Just Refl <- sameNat @nbUsedParams @0 Proxy Proxy
, Just Refl <- sameNat @nbBrushParams @2 Proxy Proxy
= computeStrokeOutline @0 @2 rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| Just Refl <- sameNat @nbUsedParams @1 Proxy Proxy
, Just Refl <- sameNat @nbBrushParams @2 Proxy Proxy
= computeStrokeOutline @1 @2 rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| Just Refl <- sameNat @nbUsedParams @2 Proxy Proxy
, Just Refl <- sameNat @nbBrushParams @2 Proxy Proxy
= computeStrokeOutline @2 @2 rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| Just Refl <- sameNat @nbUsedParams @0 Proxy Proxy
, Just Refl <- sameNat @nbBrushParams @3 Proxy Proxy
= computeStrokeOutline @0 @3 rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| Just Refl <- sameNat @nbUsedParams @1 Proxy Proxy
, Just Refl <- sameNat @nbBrushParams @3 Proxy Proxy
= computeStrokeOutline @1 @3 rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| Just Refl <- sameNat @nbUsedParams @2 Proxy Proxy
, Just Refl <- sameNat @nbBrushParams @3 Proxy Proxy
= computeStrokeOutline @2 @3 rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| Just Refl <- sameNat @nbUsedParams @3 Proxy Proxy
, Just Refl <- sameNat @nbBrushParams @3 Proxy Proxy
= computeStrokeOutline @3 @3 rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| Just Refl <- sameNat @nbUsedParams @0 Proxy Proxy
, Just Refl <- sameNat @nbBrushParams @4 Proxy Proxy
= computeStrokeOutline @0 @4 rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| Just Refl <- sameNat @nbUsedParams @1 Proxy Proxy
, Just Refl <- sameNat @nbBrushParams @4 Proxy Proxy
= computeStrokeOutline @1 @4 rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| Just Refl <- sameNat @nbUsedParams @2 Proxy Proxy
, Just Refl <- sameNat @nbBrushParams @4 Proxy Proxy
= computeStrokeOutline @2 @4 rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| Just Refl <- sameNat @nbUsedParams @3 Proxy Proxy
, Just Refl <- sameNat @nbBrushParams @4 Proxy Proxy
= computeStrokeOutline @3 @4 rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| Just Refl <- sameNat @nbUsedParams @4 Proxy Proxy
, Just Refl <- sameNat @nbBrushParams @4 Proxy Proxy
= computeStrokeOutline @4 @4 rootAlgo mbCuspOptions fitParams ptParams toBrushParams brush spline
| otherwise
= error $ unlines
[ "computeStrokeOutline_specialise: expected 'nbUsedParams <= 4', 'nbBrushParams <= 4'"
, " nbUsedParams: " ++ show (natVal' @nbUsedParams proxy#)
, "nbBrushParams: " ++ show (natVal' @nbBrushParams proxy#)
]
computeStrokeOutline ::
forall { clo :: SplineType } ( nbUsedParams :: Nat ) ( nbBrushParams :: Nat ) crvData ptData s
. ( KnownSplineType clo
, HasType ( ℝ 2 ) ptData
, HasType ( CachedStroke s ) crvData
, NFData ptData, NFData crvData
-- Differentiability.
, Interpolatable Double ( ℝ nbUsedParams )
, DiffInterp 2 NonIV nbBrushParams
, DiffInterp 3 IsIV nbBrushParams
, HasChainRule Double 3 ( ℝ nbBrushParams )
, Module 𝕀 ( T ( 𝕀ℝ nbUsedParams ) )
-- AABB computations
, Representable Double ( ℝ nbUsedParams )
, Representable 𝕀 ( 𝕀ℝ nbUsedParams )
-- Debugging.
, Show ptData, Show crvData, Show ( ℝ nbBrushParams )
)
=> RootSolvingAlgorithm
-> Maybe ( RootIsolationOptions 1 1, RootIsolationOptions 1 2, RootIsolationOptions 2 3 )
-> FitParameters
-> ( ptData -> ℝ nbUsedParams )
-> ( ℝ nbUsedParams -> ℝ nbBrushParams )
-- ^ assumed to be linear and non-decreasing
-> Brush nbBrushParams
-> Spline clo crvData ptData
-> ST s
( Either ( SplinePts Closed ) ( SplinePts Closed, SplinePts Closed )
, Seq FitPoint
, [ Cusp ]
)
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
| firstCurve :<| _ <- openCurves $ splineCurves spline
, prevCurves :|> lastCurve <- openCurves $ splineCurves spline
, firstOutlineFn :<| _ <- outlineFns
, _ :|> lastOutlineFn <- outlineFns
, let
endPt :: ptData
endPt = openCurveEnd lastCurve
startTgtFwd, startTgtBwd, endTgtFwd, endTgtBwd :: T ( ℝ 2 )
TwoSided
{ fwd = OutlinePoint { outlineTangent = startTgtFwd }
, bwd = OutlinePoint { outlineTangent = startTgtBwd }
} = outlineFn firstOutlineFn $ ℝ1 0
TwoSided
{ fwd = OutlinePoint { outlineTangent = endTgtFwd }
, bwd = OutlinePoint { outlineTangent = endTgtBwd }
} = outlineFn lastOutlineFn $ ℝ1 1
startBrush, endBrush :: SplinePts Closed
startBrush = brushShape spt0
endBrush = brushShape endPt
-- Computation of which brush segment to use for the end caps.
startTgt, endTgt :: T ( ℝ 2 )
startTgt = coords spt0 --> coords ( openCurveStart firstCurve )
endTgt = case prevCurves of
Empty -> endTangent spt0 spt0 lastCurve
_ :|> lastPrev -> endTangent spt0 ( openCurveEnd lastPrev ) lastCurve
startCap =
brushJoin ( Just CCW ) Bwd startTgtBwd ( coords spt0 , ( -1 *^ startTgt, startTgt ), startBrush ) startTgtFwd
endCap =
brushJoin ( Just CCW ) Fwd endTgtFwd ( coords endPt, ( endTgt , -1 *^ endTgt ), endBrush ) endTgtBwd
-> do
OutlineData
( TwoSided ( fwdPts, fwdFits ) ( bwdPts, bwdFits ) )
cusps
<- updateSpline ( startTgt, startTgtFwd, startTgtBwd )
pure
-- Add on start/end caps.
( Left ( adjustSplineType @Closed $ sconcat $ NE.fromList $ [ startCap, fwdPts, endCap, bwdPts ] )
, ( Seq.fromList $
[ JoinPoint Fwd $ splineStart endCap
, JoinPoint Fwd $ splineEnd startCap
, JoinPoint Bwd $ splineEnd endCap
, JoinPoint Bwd $ splineStart startCap
]
)
<> fwdFits <> bwdFits
, cusps
)
-- Closed brush path with at least one segment.
SClosed
| ClosedCurves prevCurves lastCurve <- splineCurves spline
, firstOutlineFn :<| _ <- outlineFns
, _ :|> lastOutlineFn <- outlineFns
, let
startTgt, endTgt, startTgtFwd, startTgtBwd, endTgtFwd, endTgtBwd :: T ( ℝ 2 )
startTgt = case prevCurves of
Empty -> startTangent spt0 spt0 lastCurve
firstCrv :<| _ -> startTangent spt0 spt0 firstCrv
endTgt = case prevCurves of
Empty -> endTangent spt0 spt0 lastCurve
_ :|> lastPrev -> endTangent spt0 ( openCurveEnd lastPrev ) lastCurve
TwoSided
{ fwd = OutlinePoint { outlineTangent = startTgtFwd }
, bwd = OutlinePoint { outlineTangent = startTgtBwd }
} = outlineFn firstOutlineFn $ ℝ1 0
TwoSided
{ fwd = OutlinePoint { outlineTangent = endTgtFwd }
, bwd = OutlinePoint { outlineTangent = endTgtBwd }
} = outlineFn lastOutlineFn $ ℝ1 1
fwdStartCap, bwdStartCap :: SplinePts Open
OutlineData ( fmap fst -> TwoSided fwdStartCap bwdStartCap ) _
= snd . runWriter
$ tellBrushJoin ( endTgt, endTgtFwd, endTgtBwd ) spt0 ( startTgt, startTgtFwd, startTgtBwd )
-> do
OutlineData
( TwoSided ( fwdPts, fwdFits ) ( bwdPts, bwdFits ) )
cusps
<- updateSpline ( endTgt, endTgtFwd, endTgtBwd )
pure
( Right ( adjustSplineType @Closed ( fwdStartCap <> fwdPts ), adjustSplineType @Closed ( bwdPts <> bwdStartCap ) )
, fwdFits <> bwdFits
, cusps
)
-- Single point.
_ ->
pure
( Left $ fmap ( T ( coords spt0 ) • ) ( brushShape spt0 )
, Empty
, []
)
where
outlineInfo :: ptData -> Curve Open crvData ptData -> OutlineInfo
outlineInfo = outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams brush
outlineFns :: Seq OutlineInfo
outlineFns = go spt0 ( openCurves $ splineCurves ( adjustSplineType @Open spline ) )
where
go
:: ptData
-> Seq ( Curve Open crvData ptData )
-> Seq OutlineInfo
go _ Empty = Empty
go p0 ( crv :<| crvs ) =
outlineInfo p0 crv :<| go ( openCurveEnd crv ) crvs
brushShape :: ptData -> SplinePts Closed
brushShape pt =
let Brush { brushBaseShape = brushFn, mbRotation = mbRot } = brush
brushParams = toBrushParams $ ptParams pt
in applyRotation @NonIV brushParams brushFn mbRot
updateSpline :: ( T ( ℝ 2 ), T ( ℝ 2 ), T ( ℝ 2 ) ) -> ST s OutlineData
updateSpline ( lastTgt, lastTgtFwd, lastTgtBwd )
= execWriterT
. ( `evalStateT` ( lastTgt, lastTgtFwd, lastTgtBwd ) )
$ bifoldSpline
( \ ptData curve -> do
( prevTgt, prev_tgtFwd, prev_tgtBwd ) <- get
let
fwdBwd :: OutlineInfo
fwdBwd = outlineInfo ptData curve
tgt, next_tgt, tgtFwd, next_tgtFwd, tgtBwd, next_tgtBwd :: T ( ℝ 2 )
tgt = startTangent spt0 ptData curve
next_tgt = endTangent spt0 ptData curve
TwoSided
{ fwd = OutlinePoint { outlineTangent = tgtFwd }
, bwd = OutlinePoint { outlineTangent = tgtBwd }
} = outlineFn fwdBwd $ ℝ1 0
TwoSided
{ fwd = OutlinePoint { outlineTangent = next_tgtFwd }
, bwd = OutlinePoint { outlineTangent = next_tgtBwd }
} = outlineFn fwdBwd $ ℝ1 1
lift $ tellBrushJoin ( prevTgt, prev_tgtFwd, tgtBwd ) ptData ( tgt, tgtFwd, prev_tgtBwd )
lift $ updateCurveData ( curveData curve ) fwdBwd
put ( next_tgt, next_tgtFwd, next_tgtBwd )
)
( const ( pure () ) )
( adjustSplineType @Open spline )
updateCurveData
:: crvData
-> OutlineInfo
-> WriterT OutlineData ( ST s ) ()
updateCurveData ( Lens.view ( typed @( CachedStroke s ) ) -> cache ) fwdBwd = do
let
mbCacheRef =
case cache of
NoCache -> Nothing
CachedStroke { cachedStrokeRef } -> Just cachedStrokeRef
mbOutline <- case mbCacheRef of
Nothing -> return Nothing
Just cacheRef -> lift ( readSTRef cacheRef )
case mbOutline of
-- Cached fit data is available: use it.
Just outline -> tell outline
-- No cached fit: compute the fit anew.
Nothing -> do
let
-- Split up the path at the cusps
cusps :: [ Cusp ]
cusps = outlineDefiniteCusps fwdBwd
++ outlinePotentialCusps fwdBwd
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 ) =
( sconcat $ fmap ( fitSpline Fwd fitParams ( fwd . outlineFn fwdBwd ) ) intervals
, sconcat $ fmap ( fitSpline Bwd fitParams ( bwd . outlineFn fwdBwd ) ) intervals
)
-- TODO: use foldMap1 once that's in base.
`Strats.using`
( Strats.parTuple2 Strats.rdeepseq Strats.rdeepseq )
outlineData :: OutlineData
outlineData =
OutlineData
( TwoSided fwdData ( bimap reverseSpline Seq.reverse bwdData ) )
cusps
outlineData `deepseq` tell outlineData
-- Write the computed outline to the cache.
case mbCacheRef of
Nothing -> return ()
Just cacheRef ->
lift $ writeSTRef cacheRef ( Just outlineData )
-- Connecting paths at a point of discontinuity of the tangent vector direction (G1 discontinuity).
-- This happens at corners of the brush path (including endpoints of an open brush path, where the tangent flips direction).
tellBrushJoin
:: Monad m
=> ( T ( ℝ 2 ), T ( ℝ 2 ), T ( ℝ 2 ) )
-> ptData
-> ( T ( ℝ 2 ), T ( ℝ 2 ), T ( ℝ 2 ) )
-> WriterT OutlineData m ()
tellBrushJoin ( prevTgt, prevTgtFwd, prevTgtBwd ) sp0 ( tgt, tgtFwd, tgtBwd ) =
unless noJoin $
tell $ OutlineData
( TwoSided
( fwdJoin, Seq.fromList $ map ( JoinPoint Fwd ) [ splineStart fwdJoin, splineEnd fwdJoin ] )
( bwdJoin, Seq.fromList $ map ( JoinPoint Bwd ) [ splineStart bwdJoin, splineEnd bwdJoin ] )
)
[]
where
noJoin = prevTgt `strictlyParallel` tgt
joinPointData = ( coords sp0, ( prevTgt, tgt ), brushShape sp0 )
fwdJoin = brushJoin Nothing Fwd prevTgtFwd joinPointData tgtFwd
bwdJoin = brushJoin Nothing Bwd prevTgtBwd joinPointData tgtBwd
-- | Compute a brush join, i.e. a section of a brush that connects up
-- at a corner or cap of the base path.
brushJoin :: Maybe Orientation
-> FwdBwd
-> T ( ℝ 2 )
-> ( ℝ 2, ( T ( ℝ 2 ), T ( ℝ 2 ) ), SplinePts Closed )
-> T ( ℝ 2 )
-> SplinePts Open
brushJoin mbJoinOri fwdBwd incTgt ( pt, ( incPathTgt, outPathTgt ), brushShape ) outTgt =
let
-- Flip tangents to account for regions where the outline moves
-- backwards relative to the base path (e.g. in between two cusps).
--
-- This is admittedly rather ad-hoc, but some logic is necessary to handle
-- such situations.
( incTgt', outTgt' ) =
( if case fwdBwd of { Fwd -> incTgt ^.^ incPathTgt < 0; Bwd -> incTgt ^.^ outPathTgt > 0 }
then -1 *^ incTgt
else incTgt
, if case fwdBwd of { Fwd -> outTgt ^.^ outPathTgt < 0; Bwd -> outTgt ^.^ incPathTgt > 0 }
then -1 *^ outTgt
else outTgt
)
brushOri, joinOri :: Orientation
brushOri = splineOrientation brushShape
joinOri =
case mbJoinOri of
Just ori -> ori
Nothing ->
if incTgt' × outTgt' >= 0 then CCW else CW
in
if joinOri == brushOri
then
fmap ( T pt • )
$ joinWithBrush brushShape incTgt' outTgt'
else
fmap ( T pt • )
. reverseSpline
$ joinWithBrush brushShape outTgt' incTgt'
{-# SPECIALISE computeStrokeOutline @0 @1 #-}
{-# SPECIALISE computeStrokeOutline @1 @1 #-}
{-# SPECIALISE computeStrokeOutline @0 @2 #-}
{-# SPECIALISE computeStrokeOutline @1 @2 #-}
{-# SPECIALISE computeStrokeOutline @2 @2 #-}
{-# SPECIALISE computeStrokeOutline @0 @3 #-}
{-# SPECIALISE computeStrokeOutline @1 @3 #-}
{-# SPECIALISE computeStrokeOutline @2 @3 #-}
{-# SPECIALISE computeStrokeOutline @3 @3 #-}
{-# SPECIALISE computeStrokeOutline @0 @4 #-}
{-# SPECIALISE computeStrokeOutline @1 @4 #-}
{-# SPECIALISE computeStrokeOutline @2 @4 #-}
{-# SPECIALISE computeStrokeOutline @3 @4 #-}
{-# SPECIALISE computeStrokeOutline @4 @4 #-}
-- | Computes the forward and backward stroke outline functions for a single curve.
outlineFunction
:: forall (nbUsedParams :: Nat) (nbBrushParams :: Nat) crvData ptData
. ( HasType ( ℝ 2 ) ptData
-- Differentiability.
, Interpolatable Double ( ℝ nbUsedParams )
, DiffInterp 2 NonIV nbBrushParams
, DiffInterp 3 IsIV nbBrushParams
, HasChainRule Double 3 ( ℝ nbBrushParams )
, Module 𝕀 ( T ( 𝕀ℝ nbUsedParams ) )
-- Computing AABBs
, Representable Double ( ℝ nbUsedParams )
, Representable 𝕀 ( 𝕀ℝ nbUsedParams )
-- Debugging.
, Show ptData, Show crvData, Show ( ℝ nbBrushParams )
)
=> RootSolvingAlgorithm
-> Maybe ( RootIsolationOptions 1 1, RootIsolationOptions 1 2, RootIsolationOptions 2 3 )
-> ( ptData -> ℝ nbUsedParams )
-> ( ℝ nbUsedParams -> ℝ nbBrushParams ) -- ^ assumed to be linear and non-decreasing
-> Brush nbBrushParams
-> ptData
-> Curve Open crvData ptData
-> OutlineInfo
outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams
( Brush { brushBaseShape, brushBaseShape3, brushBaseShapeI
, brushCorners, brushCorners3, brushCornersI
, mbRotation } ) = \ sp0 crv ->
let
usedParams :: C 2 ( ℝ 1 ) ( ℝ nbUsedParams )
path :: C 2 ( ℝ 1 ) ( ℝ 2 )
( path, usedParams )
= pathAndUsedParams @2 @NonIV @nbUsedParams coerce id id ptParams sp0 crv
curves :: ℝ 1 -- t
-> ( Seq ( CornerStrokeDatum 2 NonIV )
, Seq ( ℝ 1 {- s -} -> StrokeDatum 2 NonIV ) )
curves =
brushStrokeData @2 @nbBrushParams SNonIV
path
( fmap toBrushParams usedParams )
brushBaseShape
brushCorners
mbRotation
usedParams3 :: C 3 ( ℝ 1 ) ( ℝ nbUsedParams )
path3 :: C 3 ( ℝ 1 ) ( ℝ 2 )
( path3, usedParams3 )
= pathAndUsedParams @3 @NonIV @nbUsedParams coerce id id ptParams sp0 crv
curves3 :: ℝ 1 -- t
-> ( Seq ( CornerStrokeDatum 3 NonIV )
, Seq ( ℝ 1 {- s -} -> StrokeDatum 3 NonIV ) )
curves3 =
brushStrokeData @3 @nbBrushParams SNonIV
path3
( fmap toBrushParams usedParams3 )
brushBaseShape3
brushCorners3
mbRotation
curvesI :: 𝕀ℝ 1 -- t
-> ( Seq ( CornerStrokeDatum 3 IsIV )
, Seq ( 𝕀ℝ 1 {- s -} -> StrokeDatum 3 IsIV )
)
curvesI =
brushStrokeData @3 @nbBrushParams SIsIV
pathI
( fmap ( nonDecreasing toBrushParams ) usedParamsI )
brushBaseShapeI
brushCornersI
( fmap ( \ rot -> un𝕀ℝ1 . nonDecreasing ( ℝ1 . rot ) ) mbRotation )
usedParamsI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ nbUsedParams )
pathI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 )
( pathI, usedParamsI ) = pathAndUsedParams @3 @IsIV @nbUsedParams coerce point point ptParams sp0 crv
fwdBwd :: OutlineFn
fwdBwd t
= solveEnvelopeEquations rootAlgo t path_t path'_t ( fwdOffset, bwdOffset )
( curves t )
where
fwdOffset = withTangent path'_t brush_t
bwdOffset = withTangent ( -1 *^ path'_t ) brush_t
D21 path_t path'_t _ = runD path t
D21 params_t _ _ = runD usedParams t
brush_t =
applyRotation @NonIV
( toBrushParams params_t )
brushBaseShape
mbRotation
( potentialCusps :: [Cusp], definiteCusps :: [Cusp] ) =
case mbCuspOptions of
Just cuspOpts
-- Don't try to compute cusps for a trivial curve
-- (e.g. a line segment with identical start- and end-points),
-- as the root isolation code chokes on this.
| not ( trivialCurve sp0 crv )
-- TODO: introduce a maximum time budget for the cusp computation,
-- and bail out if the computation exceeds the budget.
-- (Record such bailing out and warn the user if this happens often.)
, let ( ( cornerCusps, startCornerCusps, endCornerCusps ), normalCusps ) = findCusps cuspOpts curvesI
->
-- TODO: we are applying Newton's method to the returned boxes.
--
-- This is correct for "definite cusp" boxes: Newton's method will converge.
--
-- For "potential cusp" boxes:
--
-- 1. There might not be a cusp.
-- 2. There might be multiple cusps.
foldMap
( \ ( i, ( _trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = potCusps } ) ) ->
( map ( \ ( box, _ ) -> fst $ cornerCuspCoords cornerCuspEqnPiece' ( fst . curves3 ) ( i, box ) ) potCusps
, map ( \ box -> fst $ cornerCuspCoords cornerCuspEqnPiece' ( fst . curves3 ) ( i, box ) ) defCusps )
)
( IntMap.toList cornerCusps )
<>
foldMap
( \ ( i, ( _trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = potCusps } ) ) ->
( map ( \ ( box, _ ) -> fst $ cornerCuspCoords startCornerEqnPiece ( fst . curves3 ) ( i, box ) ) potCusps
, map ( \ box -> fst $ cornerCuspCoords startCornerEqnPiece ( fst . curves3 ) ( i, box ) ) defCusps )
)
( IntMap.toList startCornerCusps )
<>
foldMap
( \ ( i, ( _trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = potCusps } ) ) ->
( map ( \ ( box, _ ) -> fst $ cornerCuspCoords endCornerEqnPiece ( fst . curves3 ) ( i, box ) ) potCusps
, map ( \ box -> fst $ cornerCuspCoords endCornerEqnPiece ( fst . curves3 ) ( i, box ) ) defCusps )
)
( IntMap.toList endCornerCusps )
<>
foldMap
( \ ( i, ( _trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = potCusps } ) ) ->
( map ( \ ( box, _ ) -> fst $ cuspCoords ( snd . curves3 ) ( i, box ) ) potCusps
, map ( \ box -> fst $ cuspCoords ( snd . curves3 ) ( i, box ) ) defCusps )
)
( IntMap.toList normalCusps )
_ ->
( [], [] )
{- DEBUG code
!_ = trace ( unlines [ "cusp finding results:"
, "definite cusps: " ++ show ( map cuspParameters definiteCusps )
, "potential cusps: " ++ show ( map cuspParameters potentialCusps )
]
) ()
!_ = trace
( unlines
[ "guessed cusp: " ++ show ( cuspParameters $ fst $ cuspCoords ( snd . curves3 ) ( 0, point $ ℝ2 0.33 0.35 ) )
]
)
()
-}
in
OutlineInfo
{ outlineFn = fwdBwd
, outlineDefiniteCusps = definiteCusps
, outlinePotentialCusps = potentialCusps
}
{-# INLINEABLE outlineFunction #-}
-- TODO: move this out
trivialCurve :: HasType ( ℝ 2 ) ptData => ptData -> Curve Open crvData ptData -> Bool
trivialCurve p0 = \case
LineTo ( NextPoint p1 ) _ -> coords p0 == coords p1
Bezier2To cp1 ( NextPoint p2 ) _ -> all ( ( == coords p0 ) . coords ) [ cp1, p2 ]
Bezier3To cp1 cp2 ( NextPoint p3 ) _ -> all ( ( == coords p0 ) . coords ) [ cp1, cp2, p3 ]
pathAndUsedParams :: forall k i (nbUsedParams :: Nat) arr crvData ptData
. ( HasType ( ℝ 2 ) ptData
, HasBézier k i
, arr ~ C k
, Module ( I i Double ) ( T ( I i 2 ) )
, Torsor ( T ( I i 2 ) ) ( I i 2 )
, Module ( I i Double ) ( T ( I i nbUsedParams ) )
, Torsor ( T ( I i nbUsedParams ) ) ( I i nbUsedParams )
, Module Double ( T ( ℝ nbUsedParams ) )
, Representable Double ( ℝ nbUsedParams )
, Representable 𝕀 ( 𝕀ℝ nbUsedParams )
)
=> ( I i 1 -> I i Double )
-> ( ℝ 2 -> I i 2 )
-> ( ℝ nbUsedParams -> I i nbUsedParams )
-> ( ptData -> ℝ nbUsedParams )
-> ptData
-> Curve Open crvData ptData
-> ( I i 1 `arr` I i 2, I i 1 `arr` I i nbUsedParams )
pathAndUsedParams co toI1 toI2 ptParams sp0 crv =
case crv of
LineTo { curveEnd = NextPoint sp1 }
| let seg = Segment sp0 sp1
-> ( line @k @i @2 co ( fmap ( toI1 . coords ) seg )
, line @k @i @nbUsedParams co ( fmap ( toI2 . ptParams ) seg ) )
Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 }
| let bez2 = Quadratic.Bezier sp0 sp1 sp2
-> ( bezier2 @k @i @2 co ( fmap ( toI1 . coords ) bez2 )
, bezier2 @k @i @nbUsedParams co ( fmap ( toI2 . ptParams ) bez2 ) )
Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 }
| let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3
-> ( bezier3 @k @i @2 co ( fmap ( toI1 . coords ) bez3 )
, bezier3 @k @i @nbUsedParams co ( fmap ( toI2 . ptParams ) bez3 ) )
{-# INLINEABLE pathAndUsedParams #-}
-----------------------------------
-- Various utility functions
-- used in the "stroke" function.
-----
startTangent, endTangent
:: ( SplineTypeI clo, HasType ( ℝ 2 ) ptData )
=> ptData -- ^ start point of the spline (to handle 'NextPoint' correctly)
-> ptData -- ^ start of the curve segment
-> Curve clo crvData ptData
-> T ( ℝ 2 )
startTangent sp p0 ( LineTo mp1 _ ) = coords p0 --> coords ( fromNextPoint sp mp1 )
startTangent _ p0 ( Bezier2To p1 _ _ ) = coords p0 --> coords p1
startTangent _ p0 ( Bezier3To p1 _ _ _ ) = coords p0 --> coords p1
endTangent sp p0 ( LineTo mp1 _ ) = coords p0 --> coords ( fromNextPoint sp mp1 )
endTangent sp _ ( Bezier2To p0 mp1 _ ) = coords p0 --> coords ( fromNextPoint sp mp1 )
endTangent sp _ ( Bezier3To _ p0 mp1 _ ) = coords p0 --> coords ( fromNextPoint sp mp1 )
lastTangent :: HasType ( ℝ 2 ) ptData => Spline Closed crvData ptData -> Maybe ( T ( ℝ 2 ) )
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).
joinWithBrush
:: SplinePts Closed -> T ( ℝ 2 ) -> T ( ℝ 2 ) -> SplinePts Open
joinWithBrush brush startTgt endTgt =
joinBetweenOffsets brush startOffset endOffset
where
startOffset, endOffset :: Offset
startOffset = withTangent startTgt brush
endOffset = withTangent endTgt brush
-- | Select the section of a spline in between two offsets.
joinBetweenOffsets :: SplinePts Closed -> Offset -> Offset -> SplinePts Open
joinBetweenOffsets
spline
( Offset { offsetIndex = i1, offsetParameter = fromMaybe 1 -> t1 } )
( Offset { offsetIndex = i2, offsetParameter = fromMaybe 0 -> t2 } )
| i2 > i1
= let
pcs, lastAndRest :: Maybe ( SplinePts Open )
( pcs, lastAndRest ) =
unzip $
fmap ( splitSplineAt ( i2 - i1 ) ) ( dropCurves i1 openSpline )
in
fromMaybe empty $
mconcat
[ snd <$> ( splitFirstPiece t1 =<< pcs )
, dropFirstPiece =<< pcs
, fst <$> ( splitFirstPiece t2 =<< lastAndRest )
]
| i2 == i1 && t2 >= t1
= let
pcs :: Maybe ( SplinePts Open )
pcs = dropCurves i1 openSpline
-- We want to split: 0 -- t1 -- t2 -- 1
-- First split at t1, then split again but after rescaling
-- using t |-> ( t - t1 ) / ( 1 - t1 ).
t2'
| t1 >= 1
= 1
| otherwise
= min 1 $ ( t2 - t1 ) / ( 1 - t1 )
in
maybe empty fst
( splitFirstPiece t2' . snd =<< ( splitFirstPiece t1 =<< pcs ) )
| otherwise
= let
start, middle, end :: SplinePts Open
( ( middle, end ), start )
= first ( splitSplineAt i2 )
$ splitSplineAt i1 openSpline
in
sconcat $ NE.fromList $ catMaybes $
[ snd <$> splitFirstPiece t1 start
, dropFirstPiece start
, Just middle
, fst <$> splitFirstPiece t2 ( if i1 == i2 then start else end )
]
where
empty :: SplinePts Open
empty = Spline { splineStart = ℝ2 0 0, splineCurves = OpenCurves Empty }
openSpline :: Spline Open () ( ℝ 2 )
openSpline = adjustSplineType spline
-- | Drop the first curve in a Bézier spline.
dropFirstPiece :: SplinePts Open -> Maybe ( SplinePts Open )
dropFirstPiece ( Spline { splineCurves = OpenCurves curves } ) = case curves of
Empty -> Nothing
fstPiece :<| laterPieces ->
Just $ Spline
{ splineStart = openCurveEnd fstPiece
, splineCurves = OpenCurves laterPieces
}
-- | Subdivide the first piece at the given parameter, discarding the subsequent pieces.
splitFirstPiece :: Double -> SplinePts Open -> Maybe ( SplinePts Open, SplinePts Open )
splitFirstPiece t ( Spline { splineStart = p0, splineCurves = OpenCurves curves } ) =
case curves of
Empty -> Nothing
fstPiece :<| _
-- Handle t = 0 and t = 1 first.
| t <= 0
-- t = 0: single point, then first piece.
-> Just ( Spline { splineStart = p0, splineCurves = OpenCurves Empty }
, Spline { splineStart = p0, splineCurves = OpenCurves ( Seq.singleton fstPiece ) } )
| t >= 1
-- t = 1: first piece, then single point.
-> Just ( Spline { splineStart = p0, splineCurves = OpenCurves ( Seq.singleton fstPiece ) }
, Spline { splineStart = openCurveEnd fstPiece, splineCurves = OpenCurves Empty }
)
| otherwise
-- Normal code path: actual splitting.
-> let ( leftCurve, splitPt, rightCurve ) = subdivideCurveAt @( T ( ℝ 2 ) ) p0 fstPiece t
in Just
( Spline { splineStart = p0 , splineCurves = OpenCurves ( Seq.singleton leftCurve ) }
, Spline { splineStart = splitPt, splineCurves = OpenCurves ( Seq.singleton rightCurve ) }
)
--------------------------------------------------------------------------------
-- | Finds the point at which a convex nib (given by a piecewise Bézier curve) has the given tangent vector.
--
-- Does /not/ check that the provided nib shape is convex.
withTangent
:: forall crvData ptData
. ( HasType ( ℝ 2 ) ptData, Show crvData, Show ptData )
=> T ( ℝ 2 ) -> Spline Closed crvData ptData -> Offset
withTangent tgt_wanted spline@( Spline { splineStart } )
-- only allow non-empty splines
| Just tgt_last <- lastTangent spline
-- only allow well-defined query tangent vectors
, not $ badTangent tgt_wanted
= case runExcept . ( `runStateT` tgt_last ) $ ibifoldSpline go ( \ _ -> pure () ) $ adjustSplineType @Open spline of
Left off ->
off
Right _ ->
trace
( unlines
[ "withTangent: could not find any point with given tangent vector"
, "tangent vector: " <> show tgt_wanted
, "spline:"
, showSplinePoints spline
] ) $
Offset { offsetIndex = 0, offsetParameter = Just 0.5, offset = T ( coords splineStart ) }
| otherwise
= Offset { offsetIndex = 0, offsetParameter = Just 0.5, offset = T ( coords splineStart ) }
where
badTangent :: T ( ℝ 2 ) -> Bool
badTangent ( V2 tx ty ) =
isNaN tx || isNaN ty || isInfinite tx || isInfinite ty
|| ( abs tx < epsilon && abs ty < epsilon )
ori :: Orientation
ori = splineOrientation spline
go :: Int -> ptData -> Curve Open crvData ptData -> StateT ( T ( ℝ 2 ) ) ( Except Offset ) ()
go i cp cseg = do
tgt_prev <- get
let
p :: ℝ 2
p = coords cp
seg :: Curve Open crvData ( ℝ 2 )
seg = fmap coords cseg
tgt_start, tgt_end :: T ( ℝ 2 )
tgt_start = startTangent splineStart cp cseg
tgt_end = endTangent splineStart cp cseg
-- Handle corner.
unless ( tgt_prev `strictlyParallel` tgt_start ) do
for_ ( between ori tgt_prev tgt_start tgt_wanted ) \ a -> do
let
-- Decide whether the offset should be at the start of this
-- segment or at the end of the previous segment, depending
-- on which tangent is closer.
--
-- This is important to get the best initial guess possible
-- for Newton's method when solving the envelope equation.
off
| a < 0.5
= Offset
{ offsetIndex =
if i == 0
then nbCurves ( splineCurves spline ) - 1
else i - 1
, offsetParameter = Just 1
, offset = T p
}
| otherwise
= Offset
{ offsetIndex = i
, offsetParameter = Just 0
, offset = T p
}
lift $ throwE $ off
-- Handle segment.
lift $ handleSegment i p seg tgt_start
put tgt_end
handleSegment :: Int -> ℝ 2 -> Curve Open crvData ( ℝ 2 ) -> T ( ℝ 2 ) -> Except Offset ()
handleSegment i p0 ( LineTo ( NextPoint p1 ) _ ) tgt0
| tgt_wanted `strictlyParallel` tgt0
, let
offset :: T ( ℝ 2 )
offset = T $ lerp @( T ( ℝ 2 ) ) 0.5 p0 p1
= throwE ( Offset { offsetIndex = i, offsetParameter = Nothing, offset } )
| otherwise
= pure ()
handleSegment i p0 ( Bezier2To p1 ( NextPoint p2 ) _ ) tgt0 =
let
tgt1 :: T ( ℝ 2 )
tgt1 = p1 --> p2
in for_ ( convexCombination tgt0 tgt1 tgt_wanted ) \ s ->
throwE $
Offset
{ offsetIndex = i
, offsetParameter = Just s
, offset = T $ Quadratic.bezier @( T ( ℝ 2 ) ) ( Quadratic.Bezier {..} ) s
}
handleSegment i p0 ( Bezier3To p1 p2 ( NextPoint p3 ) _ ) tgt0 =
let
tgt1, tgt2 :: T ( ℝ 2 )
tgt1 = p1 --> p2
tgt2 = p2 --> p3
bez :: Cubic.Bezier ( ℝ 2 )
bez = Cubic.Bezier {..}
c01, c12, c23 :: Double
c01 = tgt_wanted × tgt0
c12 = tgt_wanted × tgt1
c23 = tgt_wanted × tgt2
correctTangentParam :: Double -> Maybe Double
correctTangentParam s
| s > -epsilon && s < 1 + epsilon
, tgt_wanted ^.^ Cubic.bezier' bez s > epsilon
= Just ( max 0 ( min 1 s ) )
| otherwise
= Nothing
in
let
mbParam :: Maybe Double
mbParam = listToMaybe
. mapMaybe correctTangentParam
$ solveQuadratic c01 ( 2 * ( c12 - c01 ) ) ( c01 + c23 - 2 * c12 )
in for_ mbParam \ s ->
throwE $
Offset
{ offsetIndex = i
, offsetParameter = Just s
, offset = T $ Cubic.bezier @( T ( ℝ 2 ) ) bez s
}
--------------------------------------------------------------------------------
splineCurveFns :: forall k i
. ( HasBézier k i
, Module ( I i Double ) ( T ( I i 2 ) )
, Torsor ( T ( I i 2 ) ) ( I i 2 ) )
=> ( I i 1 -> I i Double )
-> Spline Closed () ( I i 2 )
-> Seq ( C k ( I i 1 ) ( I i 2 ) )
splineCurveFns co spls
= runIdentity
. bifoldSpline
( \ pt crv -> Identity . Seq.singleton $ curveFn pt crv )
( const $ Identity Seq.empty )
. adjustSplineType @Open
$ spls
where
curveFn :: I i 2
-> Curve Open () ( I i 2 )
-> C k ( I i 1 ) ( I i 2 )
curveFn p0 = \case
LineTo { curveEnd = NextPoint p1 }
-> line @k @i @2 co $ Segment p0 p1
Bezier2To { controlPoint = p1, curveEnd = NextPoint p2 }
-> bezier2 @k @i @2 co $ Quadratic.Bezier p0 p1 p2
Bezier3To { controlPoint1 = p1, controlPoint2 = p2, curveEnd = NextPoint p3 }
-> bezier3 @k @i @2 co $ Cubic.Bezier p0 p1 p2 p3
newtype ZipSeq a = ZipSeq { getZipSeq :: Seq a }
deriving stock Functor
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 (k :: Nat) (nbBrushParams :: Nat) (i :: IVness) arr
. ( arr ~ C k
, HasBézier k i, HasEnvelopeEquation k
, Differentiable k i nbBrushParams
, HasChainRule ( I i Double ) k ( I i 1 )
, HasChainRule ( I i Double ) k ( I i nbBrushParams )
, Applicative ( D k 1 )
, Dim ( I i 1 ) ~ 1
, Dim ( I i nbBrushParams ) ~ nbBrushParams
, Traversable ( D k nbBrushParams )
, Traversable ( D k 1 )
, Module ( D k 1 𝕀 ) ( D k 1 ( 𝕀ℝ 2 ) )
, Transcendental ( D k 1 𝕀 )
, 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 ( ℝ nbBrushParams )
, Show ( StrokeDatum k i )
, Show ( StrokeDatum k IsIV )
, Show ( CornerStrokeDatum k i )
, Show ( CornerStrokeDatum k IsIV )
, Representable ( I i Double ) ( I i 2 ), RepDim ( I i 2 ) ~ 2
)
=> SingIVness i
-> ( I i 1 `arr` I i 2 )
-- ^ path
-> ( I i 1 `arr` I i nbBrushParams )
-- ^ brush parameters
-> ( Unrotated ( I i nbBrushParams `arr` Spline Closed () ( I i 2 ) ) )
-- ^ brush from parameters
-> ( Seq ( Unrotated ( I i nbBrushParams `arr` Corner ( I i 2 ) ) ) )
-- ^ brush corners from parameters
-> ( Maybe ( I i nbBrushParams -> I i Double ) )
-- ^ rotation parameter
-> ( I i 1 -> ( Seq ( CornerStrokeDatum k i ), Seq ( I i 1 -> StrokeDatum k i ) ) )
brushStrokeData ivness path params brush brushCorners mbBrushRotation =
\ t ->
let
co :: I i 1 -> I i Double
!co = case ivness of
SNonIV -> coerce
SIsIV -> coerce
dpath_t :: D k 1 ( I i 2 )
!dpath_t = runD path t
dparams_t :: D k 1 ( I i nbBrushParams )
!dparams_t = runD params t
brushParams :: I i nbBrushParams
!brushParams = value @k @1 dparams_t
dbrush_params :: D k nbBrushParams ( Spline Closed () ( I i 2 ) )
!dbrush_params = runD ( getUnrotated brush ) brushParams
splines :: Seq ( D k nbBrushParams ( I i 1 `arr` I i 2 ) )
!splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @k @i co ) dbrush_params
dbrushes_t :: Seq ( I i 1 -> Unrotated ( D k 2 ( I i 2 ) ) )
!dbrushes_t =
force $
fmap
( fmap ( Unrotated . uncurryD @k )
. ( \ db s -> fmap ( `runD` s ) db )
. chain @( I i Double ) @k dparams_t
-- This is the crucial use of the chain rule.
)
splines
in ( brushCorners <&> \ cornerFn ->
let !dcorner = runD ( getUnrotated cornerFn ) brushParams
!dcorner_t = Unrotated $ chain @( I i Double ) @k dparams_t dcorner
in mkCornerStrokeDatum dpath_t dparams_t dcorner_t
, fmap ( mkStrokeDatum dpath_t dparams_t . ) dbrushes_t )
where
mkStrokeDatum :: D k 1 ( I i 2 )
-> D k 1 ( I i nbBrushParams )
-> Unrotated ( D k 2 ( I i 2 ) )
-> StrokeDatum k i
mkStrokeDatum dpath_t dparams_t dbrush_t =
let mbRotation = mbBrushRotation <&> \ getTheta -> fmap getTheta dparams_t
in envelopeEquation @k ivness dpath_t dbrush_t mbRotation
mkCornerStrokeDatum :: D k 1 ( I i 2 )
-> D k 1 ( I i nbBrushParams )
-> Unrotated ( D k 1 ( Corner ( I i 2 ) ) )
-> CornerStrokeDatum k i
mkCornerStrokeDatum dpath_t dparams_t dcorner_t =
let mbRotation = mbBrushRotation <&> \ getTheta -> fmap getTheta dparams_t
in cornerEnvelopeEquation @k dpath_t dcorner_t mbRotation
{-# INLINEABLE brushStrokeData #-}
{-# SPECIALISE brushStrokeData @2 @1 @NonIV #-}
{-# SPECIALISE brushStrokeData @2 @2 @NonIV #-}
{-# SPECIALISE brushStrokeData @2 @3 @NonIV #-}
{-# SPECIALISE brushStrokeData @2 @4 @NonIV #-}
{-# SPECIALISE brushStrokeData @3 @1 @IsIV #-}
{-# SPECIALISE brushStrokeData @3 @2 @IsIV #-}
{-# SPECIALISE brushStrokeData @3 @3 @IsIV #-}
{-# SPECIALISE brushStrokeData @3 @4 @IsIV #-}
--------------------------------------------------------------------------------
-- Solving the envelolpe equation: root-finding.
-- | Which method to use to solve the envelope equation?
data RootSolvingAlgorithm
-- | Use the Newton–Raphson method.
= NewtonRaphson { maxIters :: Word, precision :: Int }
-- | Use the modified Halley M2 method.
| HalleyM2 { maxIters :: Word, precision :: Int }
-- | Solve the envelope equations at a given point \( t = t_0 \), to find
-- \( s_0 \) such that \( c(t_0, s_0) \) is on the envelope of the brush stroke.
solveEnvelopeEquations :: RootSolvingAlgorithm
-> ℝ 1 -- ^ @t@ (for debugging only)
-> ℝ 2
-> T ( ℝ 2 )
-> ( Offset, Offset )
-> ( Seq ( CornerStrokeDatum 2 NonIV ), Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) )
-> TwoSided OutlinePoint
solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffset ) ( corners, strokeData )
= TwoSided
{ fwd = fromMaybe fwdSol mbFwdCornerSol
, bwd = Lens.over ( field' @"outlineTangent" ) ( -1 *^ )
$ fromMaybe bwdSol mbBwdCornerSol
}
where
fwdSol = findSolFrom Fwd fwdOffset
bwdSol = findSolFrom Bwd bwdOffset
mbFwdCornerSol = getFirst $ foldMap ( First . cornerSol Fwd ) corners
mbBwdCornerSol = getFirst $ foldMap ( First . cornerSol Bwd ) corners
findSolFrom :: FwdBwd -> Offset -> OutlinePoint
findSolFrom fwdOrBwd ( Offset { offsetIndex = i00, offsetParameter = s00 } ) =
go ( fromIntegral i00 + maybe 0.5 ( min 0.99999999999999989 ) s00 )
-- NB: the 'fromDomain' function requires values in the
-- half-open interval [0,1[, hence the @min 0.99999999999999989@ above.
where
go :: Double -> OutlinePoint
go is0 =
case sol fwdOrBwd strokeData is0 of
( goodSoln, pt )
| goodSoln
-> pt
| otherwise
-> trace
( unlines
[ "solveEnvelopeEquations: bad solution"
, "falling back to naive guess; this will lead to inaccurate fitting"
, " t: " ++ show _t
, " dir: " ++ show fwdOrBwd
, " p: " ++ show path_t
, " p': " ++ show path'_t
, " is0: " ++ show is0
, " BAD pt: " ++ show pt
]
) finish fwdOrBwd strokeData is0
sol :: FwdBwd -> Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) -> Double -> ( Bool, OutlinePoint )
sol fwdOrBwd f is0 =
let solRes = runSolveMethod ( eqn fwdOrBwd f ) is0
( good, is ) =
case solRes of
Nothing -> ( False, is0 )
Just is1 -> ( case fwdOrBwd of
Fwd -> outlineDotProduct outlinePt >= 0
Bwd -> outlineDotProduct outlinePt <= 0
, is1
)
outlinePt = finish fwdOrBwd f is
in ( good, outlinePt )
runSolveMethod = case rootAlgo of
HalleyM2 { maxIters, precision } ->
halleyM2 maxIters precision
NewtonRaphson { maxIters, precision } ->
newtonRaphson maxIters precision
finish :: FwdBwd -> Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) -> Double -> OutlinePoint
finish _fwdOrBwd fs is =
let (_i, s) = fromDomain is in
case evalStrokeDatum fs is of -- TODO: a bit redundant to have to compute this again...
StrokeDatum
{ stroke
, mbRotation
, ee = D12 ( ℝ1 _ee ) _ ( T ( ℝ1 𝛿E𝛿s ) )
, du = D12 u _ _
, dv = D12 v _ _
, 𝛿E𝛿s_unrot_dcdt = Unrotated ( D0 𝛿E𝛿s_unrot_dcdt )
} ->
let
-- The total derivative dc/dt is computed by dividing by ∂E/∂s.
-- However:
--
-- - the outline fitting algorithm only cares about the direction
-- of the tangent and not its magnitude,
-- - ∂E/∂s is a scalar
--
-- Thus we only care about the sign of ∂E/∂s.
unrot_dcdt = ( if 𝛿E𝛿s < 0 then (-1 *^) else id ) 𝛿E𝛿s_unrot_dcdt
dcdt = case mbRotation of
Nothing -> unrot_dcdt
Just ( D21 θ _ _ ) ->
let cosθ = cos θ
sinθ = sin θ
in rotate cosθ sinθ $ unrot_dcdt
-- TODO: reduce duplication with 'cornerSol'.
in OutlinePoint
{ outlinePoint = stroke
, outlineParam = s
, outlineTangent = dcdt
, outlineDotProduct = T u ^.^ T v
}
-- Compute the dot product of u and v (which are rotated versions of ∂c/∂t and ∂c/∂s).
-- The sign of this quantity determines which side of the envelope
-- we are on.
evalStrokeDatum :: Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) -> ( Double -> StrokeDatum 2 NonIV )
evalStrokeDatum fs is =
let (i, s) = fromDomain is
in ( fs `Seq.index` i ) ( ℝ1 $ max 0 $ min 1 $ s )
eqn :: FwdBwd -> Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) -> ( Double -> (# Double, Double #) )
eqn fwdOrBwd fs is =
case evalStrokeDatum fs is of
StrokeDatum
{ du = D12 u _ u_s
, dv = D12 v _ v_s
, ee = D12 ( ℝ1 ee ) _ ( T ( ℝ1 ee_s ) ) } ->
let dp = T u ^.^ T v
dp_s = u_s ^.^ T v + T u ^.^ v_s
sgn = case fwdOrBwd of
Fwd -> 1
Bwd -> -1
in
if sgn * dp >= 0
then (# ee, ee_s #)
else
-- Perturb the envelope equation to avoid solutions
-- on the "wrong" side (e.g. with the wrong sign of the dot product).
--
-- NB: use a square root to avoid introducing double roots,
-- as roots with multiplicity >= 2 cause the Newton-Raphson algorithm
-- to only converge linearly instead of quadratically.
let ee_perturbed = sqrt $ ee * ee + dp * dp
in (# ee_perturbed, ( ee * ee_s + dp * dp_s ) / ee_perturbed #)
cornerSol :: FwdBwd -> CornerStrokeDatum 2 NonIV -> Maybe OutlinePoint
cornerSol fwdOrBwd
( CornerStrokeDatum
{ dstroke =
D21
( Corner pt startTgt endTgt )
( T ( Corner ( T -> tgt ) _ _ ) )
_
}
) =
let flipWhenBwd =
case fwdOrBwd of
Fwd -> id
Bwd -> ( -1 *^ )
tgt_flipped = flipWhenBwd tgt
in
if isNothing $ between CCW startTgt endTgt tgt_flipped
then Nothing
else
let midTgt = startTgt ^+^ endTgt
in
Just $
OutlinePoint
{ outlinePoint = pt
, outlineParam = 0
, outlineTangent = tgt
, outlineDotProduct = tgt ^.^ midTgt
}
n :: Int
n = Seq.length strokeData
fromDomain :: Double -> ( Int, Double )
fromDomain is =
let ( i0, s ) = is `divMod'` 1
in ( i0 `mod` n, s )
--------------------------------------------------------------------------------
-- Finding cusps in the envelope equation using root isolation.
-- | Computes the brush stroke coordinates of a cusp from
-- the @(t,s)@ parameter values.
cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 3 NonIV ) )
-> ( Int, 𝕀ℝ 2 )
-> ( Cusp, Bool )
cuspCoords eqs ( i, box )
= let p0 = boxMidpoint box
( ℝ2 t_final s_final )
= newton ( cuspEqnPiece eqs i ) p0
StrokeDatum
{ ee = D22 ( ℝ1 ee ) _ _ _ _ _
, 𝛿E𝛿s_unrot_dcdt = Unrotated ( D12 { _D12_v = T ( ℝ2 f g ) } )
, dpath
, stroke
}
= ( eqs ( ℝ1 t_final ) `Seq.index` i ) ( ℝ1 s_final )
in
( Cusp
{ cuspParameters = ℝ2 t_final s_final
, cuspPathCoords =
case dpath of
D31 p p' p'' _ -> D21 p p' p''
, cuspStrokeCoords = stroke
, cornerCusp = False
}
, abs ee < 1e-3 && abs f < 1e-2 && abs g < 1e-2
)
cornerCuspCoords :: ( ( I NonIV 1 -> Seq ( CornerStrokeDatum 3 NonIV ) ) -> Int -> I NonIV 1 -> D1𝔸1 Double )
-> ( ℝ 1 -> Seq ( CornerStrokeDatum 3 NonIV ) )
-> ( Int, 𝕀ℝ 1 )
-> ( Cusp, Double )
cornerCuspCoords mk_eq eqs ( i, box )
= let p0@( ℝ1 t0 ) = boxMidpoint box
f t =
case mk_eq eqs i ( ℝ1 t ) of
D11 x ( T x' ) -> (# x, x' #)
mb_t_final
= newtonRaphson 40 8 f t0
ℝ1 t_final =
case mb_t_final of { Nothing -> p0; Just t -> ℝ1 t }
(# ok, _ #) = f t_final
CornerStrokeDatum
{ dpath, dstroke = D31 ( Corner stroke _ _ ) _ _ _ }
= eqs ( ℝ1 t_final ) `Seq.index` i
in
( Cusp
{ cuspParameters = ℝ2 t_final 0
, cuspPathCoords =
case dpath of
D31 p p' p'' _ -> D21 p p' p''
, cuspStrokeCoords = stroke
, cornerCusp = True
}
, ok
)
-- | Find cusps in the envelope for values of the parameters in
-- \( 0 \leqslant t, s \leqslant 1 \), using interval arithmetic.
--
-- See 'isolateRootsIn' for details of the algorithms used.
findCusps
:: ( RootIsolationOptions 1 1, RootIsolationOptions 1 2, RootIsolationOptions 2 3 )
-> ( 𝕀ℝ 1 -> ( Seq ( CornerStrokeDatum 3 IsIV ), Seq ( 𝕀ℝ 1 -> StrokeDatum 3 IsIV ) ) )
-> ( ( IntMap ( [ ( Box 1, RootIsolationTree ( Box 1 ) ) ], DoneBoxes 1 )
, IntMap ( [ ( Box 1, RootIsolationTree ( Box 1 ) ) ], DoneBoxes 1 )
, IntMap ( [ ( Box 1, RootIsolationTree ( Box 1 ) ) ], DoneBoxes 1 )
)
, IntMap ( [ ( Box 2, RootIsolationTree ( Box 2 ) ) ], DoneBoxes 2 )
)
findCusps ( opts11, opts12, opts23 ) boxStrokeData =
( findCornerCuspsIn opts11 opts12 ( fst . boxStrokeData )
( IntMap.fromList
[ ( i, [ 𝕀ℝ1 unit ] )
| i <- [ 0 .. length ( fst $ boxStrokeData ( 𝕀ℝ1 unit ) ) - 1 ]
]
)
, findCuspsIn opts23 ( snd . boxStrokeData )
( IntMap.fromList
[ ( i, [ 𝕀ℝ2 unit unit ] )
| i <- [ 0 .. length ( snd $ boxStrokeData ( 𝕀ℝ1 unit ) ) - 1 ]
]
)
)
where
unit :: 𝕀
unit = 𝕀 0 1
{-# INLINE unit #-}
-- | Like 'findCusps', but parametrised over the initial boxes for the
-- root isolation method.
findCuspsIn
:: RootIsolationOptions N 3
-> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 IsIV ) )
-> IntMap [ Box 2 ]
-> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], DoneBoxes N )
findCuspsIn opts eqs =
IntMap.mapWithKey
( \ i -> foldMap ( isolateRootsIn opts ( cuspEqnPieceI eqs i ) ) )
-- TODO: code duplication between 'cuspEqnPiece' and 'cuspEqnPieceI'.
cuspEqnPiece :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 3 NonIV ) )
-> Int
-> ℝ 2
-> D1𝔸2 (ℝ 3)
cuspEqnPiece eqs i ( ℝ2 t s ) =
let StrokeDatum
{ ee =
D22 { _D22_v = ℝ1 ee
, _D22_dx = T ( ℝ1 ee_t )
, _D22_dy = T ( ℝ1 ee_s ) }
, 𝛿E𝛿s_unrot_dcdt =
Unrotated
D12 { _D12_v = T ( ℝ2 f g )
, _D12_dx = T ( T ( ℝ2 f_t g_t ) )
, _D12_dy = T ( T ( ℝ2 f_s g_s ) ) }
} = ( eqs ( ℝ1 t ) `Seq.index` i ) ( ℝ1 s )
in D12 ( ℝ3 ee f g )
( T $ ℝ3 ee_t f_t g_t )
( T $ ℝ3 ee_s f_s g_s )
cuspEqnPieceI :: ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 IsIV ) )
-> Int -> 𝕀ℝ 2 -> D1𝔸2 ( 𝕀ℝ 3 )
cuspEqnPieceI eqs i ( 𝕀ℝ2 t s ) =
let StrokeDatum
{ ee =
D22 { _D22_v = 𝕀ℝ1 ee
, _D22_dx = T ( 𝕀ℝ1 ee_t )
, _D22_dy = T ( 𝕀ℝ1 ee_s ) }
, 𝛿E𝛿s_unrot_dcdt =
Unrotated
D12 { _D12_v = T ( 𝕀ℝ2 f g )
, _D12_dx = T ( T ( 𝕀ℝ2 f_t g_t ) )
, _D12_dy = T ( T ( 𝕀ℝ2 f_s g_s ) ) }
} = ( eqs ( 𝕀ℝ1 t ) `Seq.index` i ) ( 𝕀ℝ1 s )
in D12 ( 𝕀ℝ3 ee f g )
( T $ 𝕀ℝ3 ee_t f_t g_t )
( T $ 𝕀ℝ3 ee_s f_s g_s )
-- | Like 'findCuspsIn' but in the case that the envelope is traced out
-- by a brush corner.
findCornerCuspsIn
:: RootIsolationOptions 1 1
-> RootIsolationOptions 1 2
-> ( 𝕀ℝ 1 -> Seq ( CornerStrokeDatum 3 IsIV ) )
-> IntMap [ Box 1 ]
-> ( IntMap ( [ ( Box 1, RootIsolationTree ( Box 1 ) ) ], DoneBoxes 1 )
, IntMap ( [ ( Box 1, RootIsolationTree ( Box 1 ) ) ], DoneBoxes 1 )
, IntMap ( [ ( Box 1, RootIsolationTree ( Box 1 ) ) ], DoneBoxes 1 )
)
findCornerCuspsIn opts1 opts2 eqs initBoxes =
( IntMap.mapWithKey
( \ i -> foldMap ( isolateRootsIn opts2 ( cornerCuspEqnPiece eqs i ) ) )
initBoxes
, IntMap.mapWithKey
( \ i -> foldMap ( isolateRootsIn opts1 ( fmap 𝕀ℝ1 . startCornerEqnPiece eqs i ) ) )
initBoxes
, IntMap.mapWithKey
( \ i -> foldMap ( isolateRootsIn opts1 ( fmap 𝕀ℝ1 . endCornerEqnPiece eqs i ) ) )
initBoxes
)
cornerCuspEqnPiece
:: ( I i 1 -> Seq ( CornerStrokeDatum 3 i ) )
-> Int -> I i 1 -> D1𝔸1 ( I i 2 )
cornerCuspEqnPiece eqs i t =
let CornerStrokeDatum
{ dstroke =
D31 _
( T ( Corner ( T -> c' ) _ _ ) )
( T ( Corner ( T -> c'' ) _ _ ) )
_
} = ( eqs t `Seq.index` i )
in D11 ( unT c' ) c''
cornerCuspEqnPiece'
:: ( ℝ 1 -> Seq ( CornerStrokeDatum 3 NonIV ) )
-> Int -> ℝ 1 -> D1𝔸1 Double
cornerCuspEqnPiece' eqs i t =
case cornerCuspEqnPiece eqs i t of
D11 ( ℝ2 e1 e2 ) ( T ( ℝ2 de1 de2 ) ) ->
let v = Ring.sqrt ( e1 * e1 + e2 * e2 )
in D11 v ( T $ if v < Ring.fromRational 1e-12
then Ring.fromInteger 0
else ( e1 * de1 + e2 * de2 ) / v
)
startCornerEqnPiece
:: forall i. Cross ( I i Double ) ( T ( I i 2 ) )
=> ( I i 1 -> Seq ( CornerStrokeDatum 3 i ) )
-> Int -> I i 1 -> D1𝔸1 ( I i Double )
startCornerEqnPiece eqs i t =
let CornerStrokeDatum
{ dstroke =
D31 ( Corner _ t_s _ )
( T ( Corner c' t'_s _ ) )
( T ( Corner ( T -> c'' ) _ _ ) )
_
} = ( eqs t `Seq.index` i )
in parallelEqn @i ( D11 c' c'' ) ( D11 ( unT t_s ) t'_s )
endCornerEqnPiece
:: forall i. Cross ( I i Double ) ( T ( I i 2 ) )
=> ( I i 1 -> Seq ( CornerStrokeDatum 3 i ) )
-> Int -> I i 1 -> D1𝔸1 ( I i Double )
endCornerEqnPiece eqs i t =
let CornerStrokeDatum
{ dstroke =
D31 ( Corner _ _ t_e )
( T ( Corner c' _ t'_e ) )
( T ( Corner ( T -> c'' ) _ _ ) )
_
} = ( eqs t `Seq.index` i )
in parallelEqn @i ( D11 c' c'' ) ( D11 ( unT t_e ) t'_e )
-- | An equation whose roots are when the input vectors are parallel.
parallelEqn
:: Cross ( I i Double ) ( T ( I i 2 ) )
=> D1𝔸1 ( I i 2 ) -> D1𝔸1 ( I i 2 ) -> D1𝔸1 ( I i Double )
parallelEqn ( D11 ( T -> u ) u' ) ( D11 ( T -> v ) v' ) =
D11 ( u × v ) ( T $ u' × v Ring.+ u × v' )
-- TODO: this allows anti-parallel vectors, but the equations below
-- (which require strict parallelism) seem to perform worse
-- with interval arithmetic. Perhaps due to the square roots?
{-
D11 ( c Ring.- n ) ( T $ r Ring.- ( c Ring.* r Ring.+ s Ring.* t ) Ring./ n )
where
c = u ^.^ v
s = u × v
n = Ring.sqrt ( c Ring.^ 2 Ring.+ s Ring.^ 2 )
p = c Ring.- n
r = u ^.^ v' Ring.+ u' ^.^ v
t = u × v' Ring.+ u' × v
-}
--------------------------------------------------------------------------------
-- Dealing with rotation
applyRotation :: forall i nbBrushParams
. ( Representable ( I i Double ) ( I i 2 )
, RepDim ( I i 2 ) ~ 2
, Ring.Transcendental ( I i Double )
, Differential 2 ( Dim ( I i nbBrushParams ) )
)
=> I i nbBrushParams
-> BrushFn i 2 nbBrushParams
-> Maybe ( I i nbBrushParams -> I i Double )
-> Spline Closed () ( I i 2 )
applyRotation brushParams ( Unrotated unrotatedShapeFn ) mbRot =
let unrotatedShape = fun unrotatedShapeFn brushParams
in case mbRot of
Nothing -> unrotatedShape
Just getθ ->
let θ = getθ brushParams
cosθ = Ring.cos θ
sinθ = Ring.sin θ
in fmap ( unT . rotate cosθ sinθ . T ) unrotatedShape
{-# INLINEABLE applyRotation #-}