hashirama/metabrush
1
0
Fork 0
mirror of https://gitlab.com/sheaf/metabrush.git synced 2024-11-05 06:43:37 +00:00
Code Issues Actions Projects Releases Wiki Activity

handle brush corners

Browse source
This commit is contained in:
sheaf 2024-11-02 10:57:04 +01:00
parent 7a82470db9
commit 197adec8d0
13 changed files with 852 additions and 552 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

176
brush-strokes/src/lib/Calligraphy/Brushes.hs
View file

@ -4,9 +4,10 @@
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE RebindableSyntax #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE UndecidableInstances #-}
module Calligraphy.Brushes
( BrushFn, Brush(..)
( BrushFn, Brush(..), Corner(..)
, circleBrush
, ellipseBrush
, tearDropBrush
@ -21,19 +22,25 @@ import Data.Kind
import GHC.Exts
( Proxy#, proxy# )
import GHC.TypeNats
( Nat, type (<=) )
( Nat )
-- acts
import Data.Act
( Torsor((-->)) )
-- containers
import Data.Sequence
( Seq )
import qualified Data.Sequence as Seq
( empty, fromList )
( empty, fromList, singleton )
-- brush-strokes
import Math.Algebra.Dual
import Math.Bezier.Spline
import Math.Differentiable
( I )
( I, IVness(..), SingIVness (..) )
import Math.Interval
( 𝕀, singleton, point )
( 𝕀ℝ, singleton, point )
import Math.Linear
import Math.Module
( Module((^+^), (*^)) )
@ -42,18 +49,33 @@ import Math.Ring
--------------------------------------------------------------------------------
-- | The shape of a brush (before applying any rotation).
type BrushFn :: forall {kd}. kd -> Nat -> Nat -> Type
type BrushFn :: IVness -> Nat -> Nat -> Type
type BrushFn i k nbBrushParams
= C k ( I i nbBrushParams )
( Spline Closed () ( I i 2 ) )
-- | A corner of a brush.
type Corner :: Type -> Type
data Corner a
= Corner
{ cornerPoint :: a
, cornerStartTangent, cornerEndTangent :: !( T a )
}
deriving stock ( Eq, Show, Functor, Foldable, Traversable )
-- | A brush, described as a base shape + an optional rotation.
data Brush nbBrushParams
= Brush
{ -- | Base brush shape, before applying any rotation (if any).
brushBaseShape :: BrushFn ℝ 2 nbBrushParams
-- | Base brush shape, before applying any rotation (if any).
, brushBaseShapeI :: BrushFn 𝕀 3 nbBrushParams
{ -- | Base brush shape, before applying rotation (if any).
brushBaseShape :: BrushFn NonIV 2 nbBrushParams
-- | Base brush shape, before applying rotation (if any).
, brushBaseShapeI :: BrushFn IsIV 3 nbBrushParams
-- | Brush corners, before applying rotation (if any).
, brushCorners :: Seq ( C 2 ( ℝ nbBrushParams ) ( Corner ( I NonIV 2 ) ) )
-- | Brush corners, before applying rotation (if any).
, brushCornersI :: Seq ( C 3 ( 𝕀ℝ nbBrushParams ) ( Corner ( I IsIV 2 ) ) )
-- | Optional rotation angle function
-- (assumed to be a linear function).
, mbRotation :: Maybe ( ℝ nbBrushParams -> Double )
@ -63,7 +85,7 @@ data Brush nbBrushParams
-- Brushes
-- Some convenience type synonyms for brush types... a bit horrible
type ParamsICt :: Nat -> k -> Nat -> Constraint
type ParamsICt :: Nat -> IVness -> Nat -> Constraint
type ParamsICt k i nbParams =
( Module
( D k nbParams ( I i Double ) )
@ -79,8 +101,10 @@ type ParamsICt k i nbParams =
circleBrush :: Brush 1
circleBrush =
Brush
{ brushBaseShape = circleBrushFn @ℝ @2 @1 proxy# id id
, brushBaseShapeI = circleBrushFn @𝕀 @3 @1 proxy# singleton point
{ brushBaseShape = circleBrushFn @2 SNonIV
, brushBaseShapeI = circleBrushFn @3 SIsIV
, brushCorners = Seq.empty
, brushCornersI = Seq.empty
, mbRotation = Nothing
}
@ -88,8 +112,10 @@ circleBrush =
ellipseBrush :: Brush 3
ellipseBrush =
Brush
{ brushBaseShape = ellipseBrushFn @ℝ @2 @3 proxy# id id
, brushBaseShapeI = ellipseBrushFn @𝕀 @3 @3 proxy# singleton point
{ brushBaseShape = ellipseBrushFn @2 SNonIV
, brushBaseShapeI = ellipseBrushFn @3 SIsIV
, brushCorners = Seq.empty
, brushCornersI = Seq.empty
, mbRotation = Just ( `index` ( Fin 3 ) )
}
@ -97,11 +123,24 @@ ellipseBrush =
tearDropBrush :: Brush 3
tearDropBrush =
Brush
{ brushBaseShape = tearDropBrushFn @ℝ @2 @3 proxy# id id
, brushBaseShapeI = tearDropBrushFn @𝕀 @3 @3 proxy# singleton point
{ brushBaseShape = tearDropBrushFn @2 SNonIV
, brushBaseShapeI = tearDropBrushFn @3 SIsIV
, brushCorners = Seq.singleton $ tearDropCorner @2 proxy# SNonIV
, brushCornersI = Seq.singleton $ tearDropCorner @3 proxy# SIsIV
, mbRotation = Just ( `index` ( Fin 3 ) )
}
--------------------------------------------------------------------------------
type SplineFn k i nbParams
= SingIVness i
-> C k ( I i nbParams ) ( Spline 'Closed () ( I i 2 ) )
type CornerFn k i nbParams
= Proxy# k
-> SingIVness i
-> C k ( I i nbParams ) ( Corner ( I i 2 ) )
--------------------------------------------------------------------------------
-- Circle & ellipse
@ -128,15 +167,24 @@ circleSpline p = sequenceA $
Bezier3To ( p κ -1 ) ( p 1 -κ ) BackToStart ()
{-# INLINE circleSpline #-}
circleBrushFn :: forall {t} (i :: t) k nbParams
. ( 1 <= nbParams
mkI1 :: SingIVness i -> Double -> I i Double
mkI1 = \case
SNonIV -> id
SIsIV -> singleton
{-# INLINE mkI1 #-}
mkI2 :: SingIVness i -> ℝ 2 -> I i 2
mkI2 = \case
SNonIV -> id
SIsIV -> point
{-# INLINE mkI2 #-}
circleBrushFn :: forall k i nbParams
. ( nbParams ~ 1
, ParamsICt k i nbParams
)
=> Proxy# i
-> ( Double -> I i Double )
-> ( I ℝ 2 -> I i 2 )
-> C k ( I i nbParams ) ( Spline 'Closed () ( I i 2 ) )
circleBrushFn _ mkI1 mkI2 =
=> SplineFn k i nbParams
circleBrushFn ivness =
D \ params ->
let r :: D k nbParams ( I i Double )
r = runD ( var @_ @k $ Fin 1 ) params
@ -147,22 +195,20 @@ circleBrushFn _ mkI1 mkI2 =
in circleSpline mkPt
where
e_x, e_y :: D k nbParams ( I i 2 )
e_x = pure $ mkI2 $ ℝ2 1 0
e_y = pure $ mkI2 $ ℝ2 0 1
e_x = pure $ mkI2 ivness $ ℝ2 1 0
e_y = pure $ mkI2 ivness $ ℝ2 0 1
scaledBy :: D k nbParams ( I i Double ) -> Double -> D k nbParams ( I i Double )
scaledBy d x = fmap ( mkI1 x * ) d
{-# INLINEABLE circleBrushFn #-}
scaledBy d x = fmap ( mkI1 ivness x * ) d
{-# SPECIALISE circleBrushFn :: SplineFn 2 NonIV 1 #-}
{-# SPECIALISE circleBrushFn :: SplineFn 3 IsIV 1 #-}
ellipseBrushFn :: forall {t} (i :: t) k nbParams
. ( 3 <= nbParams
ellipseBrushFn :: forall k i nbParams
. ( nbParams ~ 3
, ParamsICt k i nbParams
)
=> Proxy# i
-> ( Double -> I i Double )
-> ( I ℝ 2 -> I i 2 )
-> C k ( I i nbParams ) ( Spline 'Closed () ( I i 2 ) )
ellipseBrushFn _ mkI1 mkI2 =
=> SplineFn k i nbParams
ellipseBrushFn ivness =
D \ params ->
let a, b :: D k nbParams ( I i Double )
a = runD ( var @_ @k $ Fin 1 ) params
@ -175,12 +221,13 @@ ellipseBrushFn _ mkI1 mkI2 =
in circleSpline mkPt
where
e_x, e_y :: D k nbParams ( I i 2 )
e_x = pure $ mkI2 $ ℝ2 1 0
e_y = pure $ mkI2 $ ℝ2 0 1
e_x = pure $ mkI2 ivness $ ℝ2 1 0
e_y = pure $ mkI2 ivness $ ℝ2 0 1
scaledBy :: D k nbParams ( I i Double ) -> Double -> D k nbParams ( I i Double )
scaledBy d x = fmap ( mkI1 x * ) d
{-# INLINEABLE ellipseBrushFn #-}
scaledBy d x = fmap ( mkI1 ivness x * ) d
{-# SPECIALISE ellipseBrushFn :: SplineFn 2 NonIV 3 #-}
{-# SPECIALISE ellipseBrushFn :: SplineFn 3 IsIV 3 #-}
--------------------------------------------------------------------------------
-- Tear drop
@ -202,15 +249,12 @@ tearHeight = 3 * sqrt 3 / 8
sqrt3_over_2 :: Double
sqrt3_over_2 = 0.5 * sqrt 3
tearDropBrushFn :: forall {t} (i :: t) k nbParams
. ( 3 <= nbParams
tearDropBrushFn :: forall k i nbParams
. ( nbParams ~ 3
, ParamsICt k i nbParams
)
=> Proxy# i
-> ( Double -> I i Double )
-> ( I ℝ 2 -> I i 2 )
-> C k ( I i nbParams ) ( Spline 'Closed () ( I i 2 ) )
tearDropBrushFn _ mkI1 mkI2 =
=> SplineFn k i nbParams
tearDropBrushFn ivness =
D \ params ->
let w, h :: D k nbParams ( I i Double )
w = 2 * runD ( var @_ @k ( Fin 1 ) ) params
@ -233,9 +277,41 @@ tearDropBrushFn _ mkI1 mkI2 =
BackToStart () }
where
e_x, e_y :: D k nbParams ( I i 2 )
e_x = pure $ mkI2 $ ℝ2 1 0
e_y = pure $ mkI2 $ ℝ2 0 1
e_x = pure $ mkI2 ivness $ ℝ2 1 0
e_y = pure $ mkI2 ivness $ ℝ2 0 1
scaledBy :: D k nbParams ( I i Double ) -> Double -> D k nbParams ( I i Double )
scaledBy d x = fmap ( mkI1 x * ) d
{-# INLINEABLE tearDropBrushFn #-}
scaledBy d x = fmap ( mkI1 ivness x * ) d
{-# SPECIALISE tearDropBrushFn :: SplineFn 2 NonIV 3 #-}
{-# SPECIALISE tearDropBrushFn :: SplineFn 3 IsIV 3 #-}
tearDropCorner :: forall k i nbParams. ( nbParams ~ 3, ParamsICt k i nbParams, Torsor ( T ( I i 2 ) ) ( I i 2 ) ) => CornerFn k i nbParams
tearDropCorner _ ivness =
D \ params ->
let w, h :: D k nbParams ( I i Double )
w = 2 * runD ( var @_ @k ( Fin 1 ) ) params
h = 2 * runD ( var @_ @k ( Fin 2 ) ) params
mkPt :: Double -> Double -> D k nbParams ( I i 2 )
mkPt x y
-- 1. translate the teardrop so that the centre of mass is at the origin
-- 2. scale the teardrop so that it has the requested width/height
-- 3. rotate (done separately)
= let !x' = w `scaledBy` ( x / tearWidth )
!y' = h `scaledBy` ( ( y - tearCenter ) / tearHeight )
in x' *^ e_x ^+^ y' *^ e_y
!cornerPoint = mkPt 0 0
in sequenceA $
Corner
{ cornerPoint
, cornerStartTangent = T $ liftA2 ( \ x y -> unT $ x --> y ) ( mkPt -0.5 sqrt3_over_2 ) cornerPoint
, cornerEndTangent = T $ liftA2 ( \ x y -> unT $ x --> y ) cornerPoint ( mkPt 0.5 sqrt3_over_2 )
}
where
e_x, e_y :: D k nbParams ( I i 2 )
e_x = pure $ mkI2 ivness $ ℝ2 1 0
e_y = pure $ mkI2 ivness $ ℝ2 0 1
scaledBy :: D k nbParams ( I i Double ) -> Double -> D k nbParams ( I i Double )
scaledBy d x = fmap ( mkI1 ivness x * ) d
{-# SPECIALISE tearDropCorner :: CornerFn 2 NonIV 3 #-}
{-# SPECIALISE tearDropCorner :: CornerFn 3 IsIV 3 #-}

206
brush-strokes/src/lib/Math/Algebra/Dual.hs
View file

@ -9,12 +9,15 @@
module Math.Algebra.Dual
( C(..), D, Dim
, Differential(..)
, HasChainRule(..), chainRule
, uncurryD2, uncurryD3
, linear, fun, var
, chainRuleN1Q, chainRule1NQ
, DiffModule(..)
, D𝔸0(..)
, D1𝔸1(..), D2𝔸1(..), D3𝔸1(..)
, D1𝔸2(..), D2𝔸2(..), D3𝔸2(..)
@ -30,7 +33,7 @@ import Prelude
import Data.Coerce
( coerce )
import Data.Kind
( Type )
( Constraint, Type )
import GHC.TypeNats
( Nat )
@ -94,18 +97,27 @@ instance ( Applicative ( D k ( Dim u ) ), Module r ( T v ) )
T ( D f ) ^-^ T ( D g ) = T $ D \ t -> liftA2 ( coerce $ (^-^) @r @( T v ) ) ( f t ) ( g t )
a *^ T ( D f ) = T $ D \ t -> fmap ( coerce $ (*^) @r @( T v ) a ) $ f t
newtype DiffModule k n v = DiffModule { getDiffModule :: D k n v }
instance ( Applicative ( D k n ), Module r ( T v ) )
=> Module r ( T ( DiffModule k n v ) ) where
origin = T $ DiffModule $ pure $ coerce $ origin @r @( T v )
T ( DiffModule f ) ^+^ T ( DiffModule g ) = T $ DiffModule $ liftA2 ( coerce $ (^+^) @r @( T v ) ) f g
T ( DiffModule f ) ^-^ T ( DiffModule g ) = T $ DiffModule $ liftA2 ( coerce $ (^-^) @r @( T v ) ) f g
a *^ T ( DiffModule f ) = T $ DiffModule $ fmap ( coerce $ (*^) @r @( T v ) a ) f
--------------------------------------------------------------------------------
-- TODO: split up this class into the chain rule operation
-- and all the other operations.
type Differential :: Nat -> Nat -> Constraint
class Differential k n where
konst :: AbelianGroup w => w -> D k n w
value :: D k n w -> w
-- | @HasChainRule r k v@ means we have a chain rule
-- with @D k v w@ in the middle, for any @r@-module @w@.
class HasChainRule r k v where
class Differential k ( Dim v ) => HasChainRule r k v where
chain :: Module r ( T w )
=> D k 1 v -> D k ( Dim v ) w -> D k 1 w
konst :: AbelianGroup w => w -> D k ( Dim v ) w
value :: D k ( Dim v ) w -> w
linearD :: Module r ( T w ) => ( v -> w ) -> v -> D k ( Dim v ) w
linear :: forall k r v w
@ -122,27 +134,24 @@ chainRule ( D df ) ( D dg ) =
D \ x ->
case df x of
df_x ->
chain @r @k @v df_x ( dg $ value @r @k @u df_x )
chain @r @k @v df_x ( dg $ value @k @( Dim u ) df_x )
uncurryD2 :: Dim a ~ 1
=> D 2 1 ( C 2 a b ) -> a -> D 2 2 b
uncurryD2 ( D21 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ) s0 =
let D21 b_t0s0 dbds_t0s0 d2bds2_t0s0 = b_t0 s0
D21 dbdt_t0s0 d2bdtds_t0s0 _ = dbdt_t0 s0
D21 d2bdt2_t0s0 _ _ = d2bdt2_t0 s0
uncurryD2 :: D 2 1 ( D 2 1 b ) -> D 2 2 b
uncurryD2 ( D21 ( b_t0 ) ( T ( dbdt_t0 ) ) ( T ( d2bdt2_t0 ) ) ) =
let D21 b_t0s0 dbds_t0s0 d2bds2_t0s0 = b_t0
D21 dbdt_t0s0 d2bdtds_t0s0 _ = dbdt_t0
D21 d2bdt2_t0s0 _ _ = d2bdt2_t0
in D22
b_t0s0
( T dbdt_t0s0 ) dbds_t0s0
( T d2bdt2_t0s0 ) d2bdtds_t0s0 d2bds2_t0s0
uncurryD3 :: Dim a ~ 1
=> D 3 1 ( C 3 a b ) -> a -> D 3 2 b
uncurryD3 ( D31 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ( T ( D d3bdt3_t0 ) ) ) s0 =
let D31 b_t0s0 dbds_t0s0 d2bds2_t0s0 d3bds3_t0s0 = b_t0 s0
D31 dbdt_t0s0 d2bdtds_t0s0 d3bdtds2_t0s0 _ = dbdt_t0 s0
D31 d2bdt2_t0s0 d3bdt2ds_t0s0 _ _ = d2bdt2_t0 s0
D31 d3bdt3_t0s0 _ _ _ = d3bdt3_t0 s0
uncurryD3 :: D 3 1 ( D 3 1 b ) -> D 3 2 b
uncurryD3 ( D31 b_t0 ( T dbdt_t0 ) ( T d2bdt2_t0 ) ( T d3bdt3_t0 ) ) =
let D31 b_t0s0 dbds_t0s0 d2bds2_t0s0 d3bds3_t0s0 = b_t0
D31 dbdt_t0s0 d2bdtds_t0s0 d3bdtds2_t0s0 _ = dbdt_t0
D31 d2bdt2_t0s0 d3bdt2ds_t0s0 _ _ = d2bdt2_t0
D31 d3bdt3_t0s0 _ _ _ = d3bdt3_t0
in D32
b_t0s0
( T dbdt_t0s0 ) dbds_t0s0
@ -150,8 +159,8 @@ uncurryD3 ( D31 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ( T ( D d3b
( T d3bdt3_t0s0 ) d3bdt2ds_t0s0 d3bdtds2_t0s0 d3bds3_t0s0
-- | Recover the underlying function, discarding all infinitesimal information.
fun :: forall r k v w. HasChainRule r k v => C k v w -> ( v -> w )
fun ( D df ) = value @r @k @v . df
fun :: forall k v w. Differential k ( Dim v ) => C k v w -> ( v -> w )
fun ( D df ) = value @k @( Dim v ) . df
{-# INLINE fun #-}
-- | The differentiable germ of a coordinate variable.
@ -216,7 +225,7 @@ instance ( Applicative ( D k ( Dim u ) )
negate ( ApAp2 !x ) = ApAp2 $ fmap ( negate @r ) x
-- DO NOT USE PURE!!
fromInteger !i = ApAp2 $ konst @Double @k @u ( fromInteger @r i )
fromInteger !i = ApAp2 $ konst @k @( Dim u ) ( fromInteger @r i )
deriving newtype instance AbelianGroup r => AbelianGroup ( D𝔸0 r )
@ -498,7 +507,7 @@ deriving newtype instance Field r => Field ( D𝔸0 r )
--instance Field r => Field ( D1𝔸3 r ) where
--instance Field r => Field ( D1𝔸4 r ) where
instance Field r => Field ( D2𝔸1 r ) where
fromRational q = konst @Double @2 @( ℝ 1 ) ( fromRational q )
fromRational q = konst @2 @1 ( fromRational q )
recip df =
let
fromInt = fromInteger @r
@ -511,7 +520,7 @@ instance Field r => Field ( D2𝔸1 r ) where
{-# SPECIALISE instance Field ( D2𝔸1 Double ) #-}
{-# SPECIALISE instance Field ( D2𝔸1 𝕀 ) #-}
instance Field r => Field ( D2𝔸2 r ) where
fromRational q = konst @Double @2 @( ℝ 2 ) ( fromRational q )
fromRational q = konst @2 @2 ( fromRational q )
recip df =
let
fromInt = fromInteger @r
@ -524,7 +533,7 @@ instance Field r => Field ( D2𝔸2 r ) where
{-# SPECIALISE instance Field ( D2𝔸2 Double ) #-}
{-# SPECIALISE instance Field ( D2𝔸2 𝕀 ) #-}
instance Field r => Field ( D2𝔸3 r ) where
fromRational q = konst @Double @2 @( ℝ 3 ) ( fromRational q )
fromRational q = konst @2 @3 ( fromRational q )
recip df =
let
fromInt = fromInteger @r
@ -537,7 +546,7 @@ instance Field r => Field ( D2𝔸3 r ) where
{-# SPECIALISE instance Field ( D2𝔸3 Double ) #-}
{-# SPECIALISE instance Field ( D2𝔸3 𝕀 ) #-}
instance Field r => Field ( D2𝔸4 r ) where
fromRational q = konst @Double @2 @( ℝ 4 ) ( fromRational q )
fromRational q = konst @2 @4 ( fromRational q )
recip df =
let
fromInt = fromInteger @r
@ -550,7 +559,7 @@ instance Field r => Field ( D2𝔸4 r ) where
{-# SPECIALISE instance Field ( D2𝔸4 Double ) #-}
{-# SPECIALISE instance Field ( D2𝔸4 𝕀 ) #-}
instance Field r => Field ( D3𝔸1 r ) where
fromRational q = konst @Double @3 @( ℝ 1 ) ( fromRational q )
fromRational q = konst @3 @1 ( fromRational q )
recip df =
let
fromInt = fromInteger @r
@ -563,7 +572,7 @@ instance Field r => Field ( D3𝔸1 r ) where
{-# SPECIALISE instance Field ( D3𝔸1 Double ) #-}
{-# SPECIALISE instance Field ( D3𝔸1 𝕀 ) #-}
instance Field r => Field ( D3𝔸2 r ) where
fromRational q = konst @Double @3 @( ℝ 2 ) ( fromRational q )
fromRational q = konst @3 @2 ( fromRational q )
recip df =
let
fromInt = fromInteger @r
@ -576,7 +585,7 @@ instance Field r => Field ( D3𝔸2 r ) where
{-# SPECIALISE instance Field ( D3𝔸2 Double ) #-}
{-# SPECIALISE instance Field ( D3𝔸2 𝕀 ) #-}
instance Field r => Field ( D3𝔸3 r ) where
fromRational q = konst @Double @3 @( ℝ 3 ) ( fromRational q )
fromRational q = konst @3 @3 ( fromRational q )
recip df =
let
fromInt = fromInteger @r
@ -589,7 +598,7 @@ instance Field r => Field ( D3𝔸3 r ) where
{-# SPECIALISE instance Field ( D3𝔸3 Double ) #-}
{-# SPECIALISE instance Field ( D3𝔸3 𝕀 ) #-}
instance Field r => Field ( D3𝔸4 r ) where
fromRational q = konst @Double @3 @( ℝ 4 ) ( fromRational q )
fromRational q = konst @3 @4 ( fromRational q )
recip df =
let
fromInt = fromInteger @r
@ -786,7 +795,7 @@ deriving newtype instance Transcendental r => Transcendental ( D𝔸0 r )
--instance Transcendental r => Transcendental ( D1𝔸3 r ) where
--instance Transcendental r => Transcendental ( D1𝔸4 r ) where
instance Transcendental r => Transcendental ( D2𝔸1 r ) where
pi = konst @Double @2 @( ℝ 1 ) pi
pi = konst @2 @1 pi
sin df =
let
fromInt = fromInteger @r
@ -817,7 +826,7 @@ instance Transcendental r => Transcendental ( D2𝔸1 r ) where
{-# SPECIALISE instance Transcendental ( D2𝔸1 Double ) #-}
{-# SPECIALISE instance Transcendental ( D2𝔸1 𝕀 ) #-}
instance Transcendental r => Transcendental ( D2𝔸2 r ) where
pi = konst @Double @2 @( ℝ 2 ) pi
pi = konst @2 @2 pi
sin df =
let
fromInt = fromInteger @r
@ -848,7 +857,7 @@ instance Transcendental r => Transcendental ( D2𝔸2 r ) where
{-# SPECIALISE instance Transcendental ( D2𝔸2 Double ) #-}
{-# SPECIALISE instance Transcendental ( D2𝔸2 𝕀 ) #-}
instance Transcendental r => Transcendental ( D2𝔸3 r ) where
pi = konst @Double @2 @( ℝ 3 ) pi
pi = konst @2 @3 pi
sin df =
let
fromInt = fromInteger @r
@ -880,7 +889,7 @@ instance Transcendental r => Transcendental ( D2𝔸3 r ) where
{-# SPECIALISE instance Transcendental ( D2𝔸3 𝕀 ) #-}
instance Transcendental r => Transcendental ( D2𝔸4 r ) where
pi = konst @Double @2 @( ℝ 4 ) pi
pi = konst @2 @4 pi
sin df =
let
fromInt = fromInteger @r
@ -912,7 +921,7 @@ instance Transcendental r => Transcendental ( D2𝔸4 r ) where
{-# SPECIALISE instance Transcendental ( D2𝔸4 𝕀 ) #-}
instance Transcendental r => Transcendental ( D3𝔸1 r ) where
pi = konst @Double @3 @( ℝ 1 ) pi
pi = konst @3 @1 pi
sin df =
let
fromInt = fromInteger @r
@ -944,7 +953,7 @@ instance Transcendental r => Transcendental ( D3𝔸1 r ) where
{-# SPECIALISE instance Transcendental ( D3𝔸1 𝕀 ) #-}
instance Transcendental r => Transcendental ( D3𝔸2 r ) where
pi = konst @Double @3 @( ℝ 2 ) pi
pi = konst @3 @2 pi
sin df =
let
fromInt = fromInteger @r
@ -976,7 +985,7 @@ instance Transcendental r => Transcendental ( D3𝔸2 r ) where
{-# SPECIALISE instance Transcendental ( D3𝔸2 𝕀 ) #-}
instance Transcendental r => Transcendental ( D3𝔸3 r ) where
pi = konst @Double @3 @( ℝ 3 ) pi
pi = konst @3 @3 pi
sin df =
let
fromInt = fromInteger @r
@ -1008,7 +1017,7 @@ instance Transcendental r => Transcendental ( D3𝔸3 r ) where
{-# SPECIALISE instance Transcendental ( D3𝔸3 𝕀 ) #-}
instance Transcendental r => Transcendental ( D3𝔸4 r ) where
pi = konst @Double @3 @( ℝ 4 ) pi
pi = konst @3 @4 pi
sin df =
let
fromInt = fromInteger @r
@ -1042,35 +1051,41 @@ instance Transcendental r => Transcendental ( D3𝔸4 r ) where
--------------------------------------------------------------------------------
-- HasChainRule instances.
instance Differential 2 0 where
konst w = D0 w
value ( D0 v ) = v
{-# INLINEABLE konst #-}
{-# INLINEABLE value #-}
instance HasChainRule Double 2 ( ℝ 0 ) where
konst w = D0 w
linearD f v = D0 ( f v )
value ( D0 v ) = v
chain _ ( D0 gfx ) = D21 gfx origin origin
{-# INLINEABLE konst #-}
{-# INLINEABLE linearD #-}
{-# INLINEABLE value #-}
{-# INLINEABLE chain #-}
instance HasChainRule Double 3 ( ℝ 0 ) where
instance Differential 3 0 where
konst w = D0 w
linearD f v = D0 ( f v )
value ( D0 v ) = v
chain _ ( D0 gfx ) = D31 gfx origin origin origin
{-# INLINEABLE konst #-}
{-# INLINEABLE linearD #-}
{-# INLINEABLE value #-}
instance HasChainRule Double 3 ( ℝ 0 ) where
linearD f v = D0 ( f v )
chain _ ( D0 gfx ) = D31 gfx origin origin origin
{-# INLINEABLE linearD #-}
{-# INLINEABLE chain #-}
instance HasChainRule Double 2 ( ℝ 1 ) where
instance Differential 2 1 where
konst :: forall w. AbelianGroup w => w -> D2𝔸1 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
{-# INLINEABLE konst #-}
{-# INLINEABLE value #-}
instance HasChainRule Double 2 ( ℝ 1 ) where
linearD :: forall w. Module Double ( T w ) => ( ℝ 1 -> w ) -> ℝ 1 -> D2𝔸1 w
linearD f v =
let !o = origin @Double @( T w )
@ -1087,7 +1102,6 @@ instance HasChainRule Double 2 ( ℝ 1 ) where
| otherwise
-> [|| unT o ||]
)
chain :: forall w. Module Double ( T w ) => D2𝔸1 ( ℝ 1 ) -> D2𝔸1 w -> D2𝔸1 w
chain !df !dg =
let !o = origin @Double @( T w )
@ -1096,20 +1110,28 @@ instance HasChainRule Double 2 ( ℝ 1 ) where
in $$( chainRule1NQ
[|| o ||] [|| p ||] [|| s ||]
[|| df ||] [|| dg ||] )
{-# INLINEABLE konst #-}
{-# INLINEABLE linearD #-}
{-# INLINEABLE value #-}
{-# INLINEABLE chain #-}
instance HasChainRule Double 3 ( ℝ 1 ) where
instance Differential 3 1 where
konst :: forall w. AbelianGroup w => w -> D3𝔸1 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
{-# INLINEABLE konst #-}
{-# INLINEABLE value #-}
instance HasChainRule Double 3 ( ℝ 1 ) where
chain :: forall w. Module Double ( T w ) => D3𝔸1 ( ℝ 1 ) -> D3𝔸1 w -> D3𝔸1 w
chain !df !dg =
let !o = origin @Double @( T w )
!p = (^+^) @Double @( T w )
!s = (^*) @Double @( T w )
in $$( chainRule1NQ
[|| o ||] [|| p ||] [|| s ||]
[|| df ||] [|| dg ||] )
linearD :: forall w. Module Double ( T w ) => ( ℝ 1 -> w ) -> ℝ 1 -> D3𝔸1 w
linearD f v =
let !o = origin @Double @( T w )
@ -1126,29 +1148,19 @@ instance HasChainRule Double 3 ( ℝ 1 ) where
| otherwise
-> [|| unT o ||]
)
chain :: forall w. Module Double ( T w ) => D3𝔸1 ( ℝ 1 ) -> D3𝔸1 w -> D3𝔸1 w
chain !df !dg =
let !o = origin @Double @( T w )
!p = (^+^) @Double @( T w )
!s = (^*) @Double @( T w )
in $$( chainRule1NQ
[|| o ||] [|| p ||] [|| s ||]
[|| df ||] [|| dg ||] )
{-# INLINEABLE konst #-}
{-# INLINEABLE linearD #-}
{-# INLINEABLE value #-}
{-# INLINEABLE chain #-}
instance HasChainRule Double 2 ( ℝ 2 ) where
instance Differential 2 2 where
konst :: forall w. AbelianGroup w => w -> D2𝔸2 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
{-# INLINEABLE konst #-}
{-# INLINEABLE value #-}
instance HasChainRule Double 2 ( ℝ 2 ) where
linearD :: forall w. Module Double ( T w ) => ( ℝ 2 -> w ) -> ℝ 2 -> D2𝔸2 w
linearD f v =
let !o = origin @Double @( T w )
@ -1165,7 +1177,6 @@ instance HasChainRule Double 2 ( ℝ 2 ) where
| otherwise
-> [|| unT o ||]
)
chain :: forall w. Module Double ( T w ) => D2𝔸1 ( ℝ 2 ) -> D2𝔸2 w -> D2𝔸1 w
chain !df !dg =
let !o = origin @Double @( T w )
@ -1174,20 +1185,19 @@ instance HasChainRule Double 2 ( ℝ 2 ) where
in $$( chainRule1NQ
[|| o ||] [|| p ||] [|| s ||]
[|| df ||] [|| dg ||] )
{-# INLINEABLE konst #-}
{-# INLINEABLE linearD #-}
{-# INLINEABLE value #-}
{-# INLINEABLE chain #-}
instance HasChainRule Double 3 ( ℝ 2 ) where
instance Differential 3 2 where
konst :: forall w. AbelianGroup w => w -> D3𝔸2 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
{-# INLINEABLE konst #-}
{-# INLINEABLE value #-}
instance HasChainRule Double 3 ( ℝ 2 ) where
linearD :: forall w. Module Double ( T w ) => ( ℝ 2 -> w ) -> ℝ 2 -> D3𝔸2 w
linearD f v =
let !o = origin @Double @( T w )
@ -1204,7 +1214,6 @@ instance HasChainRule Double 3 ( ℝ 2 ) where
| otherwise
-> [|| unT o ||]
)
chain :: forall w. Module Double ( T w ) => D3𝔸1 ( ℝ 2 ) -> D3𝔸2 w -> D3𝔸1 w
chain !df !dg =
let !o = origin @Double @( T w )
@ -1213,20 +1222,19 @@ instance HasChainRule Double 3 ( ℝ 2 ) where
in $$( chainRule1NQ
[|| o ||] [|| p ||] [|| s ||]
[|| df ||] [|| dg ||] )
{-# INLINEABLE konst #-}
{-# INLINEABLE linearD #-}
{-# INLINEABLE value #-}
{-# INLINEABLE chain #-}
instance HasChainRule Double 2 ( ℝ 3 ) where
instance Differential 2 3 where
konst :: forall w. AbelianGroup w => w -> D2𝔸3 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
{-# INLINEABLE konst #-}
{-# INLINEABLE value #-}
instance HasChainRule Double 2 ( ℝ 3 ) where
linearD :: forall w. Module Double ( T w ) => ( ℝ 3 -> w ) -> ℝ 3 -> D2𝔸3 w
linearD f v =
let !o = origin @Double @( T w )
@ -1243,7 +1251,6 @@ instance HasChainRule Double 2 ( ℝ 3 ) where
| otherwise
-> [|| unT o ||]
)
chain :: forall w. Module Double ( T w ) => D2𝔸1 ( ℝ 3 ) -> D2𝔸3 w -> D2𝔸1 w
chain !df !dg =
let !o = origin @Double @( T w )
@ -1252,20 +1259,19 @@ instance HasChainRule Double 2 ( ℝ 3 ) where
in $$( chainRule1NQ
[|| o ||] [|| p ||] [|| s ||]
[|| df ||] [|| dg ||] )
{-# INLINEABLE konst #-}
{-# INLINEABLE linearD #-}
{-# INLINEABLE value #-}
{-# INLINEABLE chain #-}
instance HasChainRule Double 3 ( ℝ 3 ) where
instance Differential 3 3 where
konst :: forall w. AbelianGroup w => w -> D3𝔸3 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
{-# INLINEABLE konst #-}
{-# INLINEABLE value #-}
instance HasChainRule Double 3 ( ℝ 3 ) where
linearD :: forall w. Module Double ( T w ) => ( ℝ 3 -> w ) -> ℝ 3 -> D3𝔸3 w
linearD f v =
let !o = origin @Double @( T w )
@ -1282,7 +1288,6 @@ instance HasChainRule Double 3 ( ℝ 3 ) where
| otherwise
-> [|| unT o ||]
)
chain :: forall w. Module Double ( T w ) => D3𝔸1 ( ℝ 3 ) -> D3𝔸3 w -> D3𝔸1 w
chain !df !dg =
let !o = origin @Double @( T w )
@ -1291,20 +1296,20 @@ instance HasChainRule Double 3 ( ℝ 3 ) where
in $$( chainRule1NQ
[|| o ||] [|| p ||] [|| s ||]
[|| df ||] [|| dg ||] )
{-# INLINEABLE konst #-}
{-# INLINEABLE linearD #-}
{-# INLINEABLE value #-}
{-# INLINEABLE chain #-}
{-# INLINEABLE linearD #-}
instance HasChainRule Double 2 ( ℝ 4 ) where
instance Differential 2 4 where
konst :: forall w. AbelianGroup w => w -> D2𝔸4 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
{-# INLINEABLE konst #-}
{-# INLINEABLE value #-}
instance HasChainRule Double 2 ( ℝ 4 ) where
linearD :: forall w. Module Double ( T w ) => ( ℝ 4 -> w ) -> ℝ 4 -> D2𝔸4 w
linearD f v =
let !o = origin @Double @( T w )
@ -1321,7 +1326,6 @@ instance HasChainRule Double 2 ( ℝ 4 ) where
| otherwise
-> [|| unT o ||]
)
chain :: forall w. Module Double ( T w ) => D2𝔸1 ( ℝ 4 ) -> D2𝔸4 w -> D2𝔸1 w
chain !df !dg =
let !o = origin @Double @( T w )
@ -1330,20 +1334,19 @@ instance HasChainRule Double 2 ( ℝ 4 ) where
in $$( chainRule1NQ
[|| o ||] [|| p ||] [|| s ||]
[|| df ||] [|| dg ||] )
{-# INLINEABLE konst #-}
{-# INLINEABLE linearD #-}
{-# INLINEABLE value #-}
{-# INLINEABLE chain #-}
instance HasChainRule Double 3 ( ℝ 4 ) where
instance Differential 3 4 where
konst :: forall w. AbelianGroup w => w -> D3𝔸4 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
{-# INLINEABLE konst #-}
{-# INLINEABLE value #-}
instance HasChainRule Double 3 ( ℝ 4 ) where
linearD :: forall w. Module Double ( T w ) => ( ℝ 4 -> w ) -> ℝ 4 -> D3𝔸4 w
linearD f v =
let !o = origin @Double @( T w )
@ -1360,7 +1363,6 @@ instance HasChainRule Double 3 ( ℝ 4 ) where
| otherwise
-> [|| unT o ||]
)
chain :: forall w. Module Double ( T w ) => D3𝔸1 ( ℝ 4 ) -> D3𝔸4 w -> D3𝔸1 w
chain !df !dg =
let !o = origin @Double @( T w )
@ -1369,7 +1371,5 @@ instance HasChainRule Double 3 ( ℝ 4 ) where
in $$( chainRule1NQ
[|| o ||] [|| p ||] [|| s ||]
[|| df ||] [|| dg ||] )
{-# INLINEABLE konst #-}
{-# INLINEABLE linearD #-}
{-# INLINEABLE value #-}
{-# INLINEABLE chain #-}

13
brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs
View file

@ -312,6 +312,19 @@ instance MonomialBasisQ D1𝔸1 where
Mon ( Vec [ 1 ] ) -> [|| unT $ _D11_dx $$d ||]
_ -> [|| _D11_v $$d ||]
instance MonomialBasis D1𝔸1 where
monTabulate f =
let
_D11_v = f $ Mon ( Vec [ 0 ] )
_D11_dx = T $ f $ Mon ( Vec [ 1 ] )
in D11 { .. }
{-# INLINE monTabulate #-}
monIndex d = \ case
Mon ( Vec [ 1 ] ) -> unT $ _D11_dx d
_ -> _D11_v d
{-# INLINE monIndex #-}
type instance Deg D2𝔸1 = 2
type instance Vars D2𝔸1 = 1
instance MonomialBasisQ D2𝔸1 where

2
brush-strokes/src/lib/Math/Bezier/Cubic/Fit.hs
View file

@ -99,6 +99,8 @@ data FitPoint
{ fitPoint :: !( ℝ 2 )
, fitTangent :: !( T ( ℝ 2 ) )
}
| JoinPoint
{ joinPoint :: !( ℝ 2 ) }
deriving stock ( Show, Generic )
deriving anyclass NFData

586
brush-strokes/src/lib/Math/Bezier/Stroke.hs
View file

@ -1,4 +1,5 @@
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE PartialTypeSignatures #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
@ -51,9 +52,9 @@ import Data.List.NonEmpty
( NonEmpty )
import qualified Data.List.NonEmpty as NE
import Data.Maybe
( catMaybes, fromMaybe, isJust, listToMaybe, mapMaybe )
import Data.Proxy
( Proxy(..) )
( catMaybes, fromMaybe, isJust, isNothing, listToMaybe, mapMaybe )
import Data.Monoid
( First(..) )
import Data.Semigroup
( sconcat )
import GHC.Exts
@ -108,35 +109,26 @@ import Control.Monad.Trans.Writer.CPS
-- MetaBrush
import Calligraphy.Brushes
( Brush(..) )
( Brush(..), Corner(..) )
import Math.Algebra.Dual
import qualified Math.Bezier.Cubic as Cubic
import Math.Bezier.Cubic.Fit
( FitPoint, FitParameters, fitSpline )
( FitPoint (..), FitParameters, fitSpline )
import Math.Bezier.Spline
( SplineType(..), SSplineType(..), SplineTypeI
, ssplineType, adjustSplineType, reverseSpline
, NextPoint(..), fromNextPoint
, KnownSplineType
( bifoldSpline, ibifoldSpline )
, Spline(..), SplinePts, Curves(..), Curve(..)
, openCurveStart, openCurveEnd, splitSplineAt, dropCurves
, showSplinePoints
)
import qualified Math.Bezier.Quadratic as Quadratic
import Math.Bezier.Stroke.EnvelopeEquation
import Math.Differentiable
( Differentiable, DiffInterp, I )
( I, IVness(..), SingIVness (..)
, Differentiable, DiffInterp
)
import Math.Epsilon
( epsilon )
import Math.Interval
import Math.Linear
import Math.Module
( Module(..), Inner((^.^)), Cross((×)), Interpolatable
, lerp, convexCombination, strictlyParallel
)
import Math.Ring
( Transcendental )
import qualified Math.Ring as Ring
import Math.Orientation
( Orientation(..), splineOrientation
, between
@ -240,9 +232,10 @@ computeStrokeOutline ::
-- Differentiability.
, Interpolatable Double ( ℝ nbUsedParams )
, DiffInterp 2 ℝ nbBrushParams
, DiffInterp 3 𝕀 nbBrushParams
, DiffInterp 2 NonIV nbBrushParams
, DiffInterp 3 IsIV nbBrushParams
, HasChainRule Double 2 ( ℝ nbUsedParams )
, HasChainRule 𝕀 2 ( 𝕀ℝ nbBrushParams )
, HasChainRule 𝕀 3 ( 𝕀ℝ nbUsedParams )
, HasChainRule Double 2 ( ℝ nbBrushParams )
, HasChainRule 𝕀 3 ( 𝕀ℝ nbBrushParams )
@ -254,13 +247,13 @@ computeStrokeOutline ::
-- Debugging.
, Show ptData, Show crvData, Show ( ℝ nbBrushParams )
)
=> RootSolvingAlgorithm
-> Maybe ( RootIsolationOptions N 3 )
-> Maybe ( RootIsolationOptions 1 2, RootIsolationOptions 2 3 )
-> FitParameters
-> ( ptData -> ℝ nbUsedParams )
-> ( ℝ nbUsedParams -> ℝ nbBrushParams ) -- ^ assumed to be linear and non-decreasing
-> ( ℝ nbUsedParams -> ℝ nbBrushParams )
-- ^ assumed to be linear and non-decreasing
-> Brush nbBrushParams
-> Spline clo crvData ptData
-> ST s
@ -302,21 +295,29 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru
etx, ety :: Double
V2 etx ety = endTgt
startCap, endCap :: SplinePts Open
startCap
mbStartCap, mbEndCap :: Maybe ( SplinePts Open )
mbStartCap
| startTgtBwd `strictlyParallel` startTgtFwd
= Nothing
| isJust $ between CCW startTgtBwd startTgtFwd startTestTgt
= fmap ( T ( coords spt0 ) • )
= Just
. fmap ( T ( coords spt0 ) • )
$ joinWithBrush startBrush startTgtBwd startTgtFwd
| otherwise
= fmap ( T ( coords spt0 ) • )
= Just
. fmap ( T ( coords spt0 ) • )
. reverseSpline
$ joinWithBrush startBrush startTgtFwd startTgtBwd
endCap
mbEndCap
| endTgtBwd `strictlyParallel` endTgtFwd
= Nothing
| isJust $ between CCW endTgtBwd endTgtFwd endTestTgt
= fmap ( T ( coords endPt ) • )
= Just
. fmap ( T ( coords endPt ) • )
$ joinWithBrush endBrush endTgtFwd endTgtBwd
| otherwise
= fmap ( T ( coords endPt ) • )
= Just
. fmap ( T ( coords endPt ) • )
. reverseSpline
$ joinWithBrush endBrush endTgtBwd endTgtFwd
@ -326,8 +327,14 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru
cusps
<- updateSpline ( startTgt, startTgtFwd, startTgtBwd )
pure
( Left ( adjustSplineType @Closed $ startCap <> fwdPts <> endCap <> bwdPts )
, fwdFits <> bwdFits
( Left ( adjustSplineType @Closed $ sconcat $ NE.fromList $ catMaybes [ mbStartCap, Just fwdPts, mbEndCap, Just bwdPts ] )
, ( Seq.fromList $ fmap JoinPoint $ catMaybes
[ fmap splineEnd mbStartCap
, fmap splineStart mbEndCap
, fmap splineEnd mbEndCap
, fmap splineStart mbStartCap ]
)
<> fwdFits <> bwdFits
, cusps
)
-- Closed brush path with at least one segment.
@ -387,7 +394,7 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru
brushShape pt =
let Brush { brushBaseShape = shapeFn, mbRotation = mbRot } = brush
brushParams = toBrushParams $ ptParams pt
shape = fun @Double shapeFn brushParams
shape = fun shapeFn brushParams
in case mbRot of
Nothing -> shape
Just getθ ->
@ -463,7 +470,10 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru
-> WriterT OutlineData m ()
tellBrushJoin ( prevTgt, prevTgtFwd, prevTgtBwd ) sp0 ( tgt, tgtFwd, tgtBwd ) =
tell $ OutlineData
( TwoSided ( fwdJoin, Empty ) ( bwdJoin, Empty ) )
( TwoSided
( fwdJoin, Seq.fromList $ map JoinPoint [ splineStart fwdJoin, splineEnd fwdJoin ] )
( bwdJoin, Seq.fromList $ map JoinPoint [ splineStart bwdJoin, splineEnd bwdJoin ] )
)
[]
where
ptOffset :: T ( ℝ 2 )
@ -524,9 +534,10 @@ outlineFunction
-- Differentiability.
, Interpolatable Double ( ℝ nbUsedParams )
, DiffInterp 2 ℝ nbBrushParams
, DiffInterp 3 𝕀 nbBrushParams
, DiffInterp 2 NonIV nbBrushParams
, DiffInterp 3 IsIV nbBrushParams
, HasChainRule Double 2 ( ℝ nbUsedParams )
, HasChainRule 𝕀 2 ( 𝕀ℝ nbBrushParams )
, HasChainRule 𝕀 3 ( 𝕀ℝ nbUsedParams )
, HasChainRule Double 2 ( ℝ nbBrushParams )
, HasChainRule 𝕀 3 ( 𝕀ℝ nbBrushParams )
@ -542,7 +553,7 @@ outlineFunction
, Show ptData, Show crvData, Show ( ℝ nbBrushParams )
)
=> RootSolvingAlgorithm
-> Maybe ( RootIsolationOptions N 3 )
-> Maybe ( RootIsolationOptions 1 2, RootIsolationOptions 2 3 )
-> ( ptData -> ℝ nbUsedParams )
-> ( ℝ nbUsedParams -> ℝ nbBrushParams ) -- ^ assumed to be linear and non-decreasing
-> Brush nbBrushParams
@ -550,37 +561,42 @@ outlineFunction
-> Curve Open crvData ptData
-> OutlineInfo
outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams
( Brush { brushBaseShape, brushBaseShapeI, mbRotation } ) = \ sp0 crv ->
( Brush { brushBaseShape, brushBaseShapeI
, brushCorners, brushCornersI
, mbRotation } ) = \ sp0 crv ->
let
usedParams :: C 2 ( ℝ 1 ) ( ℝ nbUsedParams )
path :: C 2 ( ℝ 1 ) ( ℝ 2 )
( path, usedParams )
= pathAndUsedParams @2 @ℝ @nbUsedParams coerce id id ptParams sp0 crv
= pathAndUsedParams @2 @NonIV @nbUsedParams coerce id id ptParams sp0 crv
curves :: ℝ 1 -- t
-> Seq ( ℝ 1 {- s -} -> StrokeDatum 2 ℝ )
-> ( Seq ( CornerStrokeDatum 2 NonIV )
, Seq ( ℝ 1 {- s -} -> StrokeDatum 2 NonIV ) )
curves =
brushStrokeData @2 @nbBrushParams
coerce coerce
brushStrokeData @2 @nbBrushParams SNonIV
path
( fmap toBrushParams usedParams )
brushBaseShape
brushCorners
mbRotation
curvesI :: 𝕀ℝ 1 -- t
-> Seq ( 𝕀ℝ 1 {- s -} -> StrokeDatum 3 𝕀 )
-> ( Seq ( CornerStrokeDatum 3 IsIV )
, Seq ( 𝕀ℝ 1 {- s -} -> StrokeDatum 3 IsIV )
)
curvesI =
brushStrokeData @3 @nbBrushParams
coerce coerce
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 @𝕀 @nbUsedParams coerce point point ptParams sp0 crv
( pathI, usedParamsI ) = pathAndUsedParams @3 @IsIV @nbUsedParams coerce point point ptParams sp0 crv
fwdBwd :: OutlineFn
fwdBwd t
@ -604,12 +620,12 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams
in fmap ( unT . rotate cosθ sinθ . T )
in
applyRot $
value @Double @2 @( ℝ nbBrushParams ) $
value @2 @nbBrushParams $
runD brushBaseShape brushParams
( potentialCusps, definiteCusps ) =
( potentialCusps :: [Cusp], definiteCusps :: [Cusp] ) =
case mbCuspOptions of
Just cuspOptions
Just ( cusp1Opts, cusp2Opts )
-- 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.
@ -617,19 +633,28 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams
-- 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, normalCusps ) = findCusps ( cusp1Opts, cusp2Opts ) curvesI
->
foldMap
( \ ( i, ( _trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = potCusps } ) ) ->
( map ( ( i , ) . fst ) potCusps, map ( i , ) defCusps )
( map ( \ ( box, _ ) -> cornerCuspCoords ( fst . curves ) ( i, box ) ) potCusps
, map ( \ box -> cornerCuspCoords ( fst . curves ) ( i, box ) ) defCusps )
)
( IntMap.toList $ findCusps cuspOptions curvesI )
( IntMap.toList cornerCusps )
<>
foldMap
( \ ( i, ( _trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = potCusps } ) ) ->
( map ( \ ( box, _ ) -> cuspCoords ( snd . curves ) ( i, box ) ) potCusps
, map ( \ box -> cuspCoords ( snd . curves ) ( i, box ) ) defCusps )
)
( IntMap.toList normalCusps )
_ ->
( [], [] )
in OutlineInfo
{ outlineFn = fwdBwd
, outlineDefiniteCusps = map ( cuspCoords curves ) definiteCusps
, outlinePotentialCusps = map ( cuspCoords curves ) potentialCusps
, outlineDefiniteCusps = definiteCusps
, outlinePotentialCusps = potentialCusps
}
{-# INLINEABLE outlineFunction #-}
@ -648,15 +673,15 @@ pathAndUsedParams :: forall k i (nbUsedParams :: Nat) arr crvData ptData
, 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 ( I ℝ nbUsedParams ) )
, Torsor ( T ( ( I ℝ nbUsedParams ) ) ) ( I ℝ nbUsedParams )
, Module Double ( T ( ℝ nbUsedParams ) )
, Torsor ( T ( ( ℝ nbUsedParams ) ) ) ( ℝ nbUsedParams )
, Representable Double ( ℝ nbUsedParams )
, Representable 𝕀 ( 𝕀ℝ nbUsedParams )
)
=> ( I i 1 -> I i Double )
-> ( I ℝ 2 -> I i 2 )
-> ( I ℝ nbUsedParams -> I i nbUsedParams )
-> ( ptData -> I ℝ nbUsedParams )
-> ( ℝ 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 )
@ -1004,15 +1029,19 @@ instance Applicative ZipSeq where
liftA2 f ( ZipSeq xs ) ( ZipSeq ys ) = ZipSeq ( Seq.zipWith f xs ys )
{-# INLINE liftA2 #-}
brushStrokeData :: forall {kd} (k :: Nat) (nbBrushParams :: Nat) (i :: kd) arr
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 𝕀 k ( 𝕀ℝ 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 ) )
@ -1020,88 +1049,161 @@ brushStrokeData :: forall {kd} (k :: Nat) (nbBrushParams :: Nat) (i :: kd) arr
, Torsor ( T ( I i 2 ) ) ( I i 2 )
, Show ( ℝ nbBrushParams )
, Show ( StrokeDatum k i )
, Show ( StrokeDatum k IsIV )
, Representable ( I i Double ) ( I i 2 ), RepDim ( I i 2 ) ~ 2
)
=> ( I i Double -> I i 1 )
-> ( I i 1 -> I i Double )
=> SingIVness i
-> ( I i 1 `arr` I i 2 )
-- ^ path
-> ( I i 1 `arr` I i nbBrushParams )
-- ^ brush parameters
-> ( I i nbBrushParams `arr` Spline Closed () ( I i 2 ) )
-- ^ brush from parameters
-> ( Seq ( 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 ( I i 1 -> StrokeDatum k i ) )
brushStrokeData co1 co2 path params brush mbBrushRotation =
-> ( 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 brush $ value @( I i Double ) @k @( I i 1 ) dparams_t
!dbrush_params = runD brush brushParams
splines :: Seq ( D k nbBrushParams ( I i 1 `arr` I i 2 ) )
!splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @k @i co2 ) dbrush_params
!splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @k @i co ) dbrush_params
dbrushes_t :: Seq ( I i 1 -> D k 2 ( I i 2 ) )
!dbrushes_t = force $ fmap ( uncurryD @k . chain @( I i Double ) @k dparams_t ) splines
!dbrushes_t =
force $
fmap
( fmap ( 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 cornerFn brushParams
!corner = value @k @nbBrushParams dcorner
!dcornerPt = chain @_ @k dparams_t $ fmap cornerPoint dcorner
in mkCornerStrokeDatum dpath_t dparams_t corner dcornerPt
, fmap ( mkStrokeDatum dpath_t dparams_t . ) dbrushes_t )
in fmap ( mkStrokeDatum dpath_t dparams_t ) dbrushes_t
where
mkStrokeDatum :: D k 1 ( I i 2 )
-> D k 1 ( I i nbBrushParams )
-> ( I i 1 -> D k 2 ( I i 2 ) )
-> ( I i 1 -> StrokeDatum k i )
mkStrokeDatum dpath_t dparams_t dbrush_t s =
let dbrush_t_s = dbrush_t s
mbRotation = mbBrushRotation <&> \ getTheta -> fmap getTheta dparams_t
in envelopeEquation @k @i Proxy co1 dpath_t dbrush_t_s mbRotation
{-# INLINEABLE brushStrokeData #-}
-> 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
{-
{-# SPECIALISE brushStrokeData
:: ( 𝕀 -> ( 𝕀ℝ 1 ) )
-> ( 𝕀ℝ 1 -> 𝕀 )
mkCornerStrokeDatum :: D k 1 ( I i 2 )
-> D k 1 ( I i nbBrushParams )
-> Corner ( I i 2 )
-> D k 1 ( I i 2 )
-> CornerStrokeDatum k i
mkCornerStrokeDatum dpath_t dparams_t corner dbrush_t =
let mbRotation = mbBrushRotation <&> \ getTheta -> fmap getTheta dparams_t
in CornerStrokeDatum
{ dpath = dpath_t
, dbrush = dbrush_t
, corner
, stroke = unT $
T ( value @k @1 dpath_t )
^+^ ( let b = T $ value @k @1 dbrush_t in case mbRotation of
Nothing -> b
Just dθ -> let { θ = value @k @1 dθ; cosθ = Ring.cos θ; sinθ = Ring.sin θ } in rotate cosθ sinθ b
)
, mbRotation }
{-# SPECIALISE INLINE brushStrokeData
:: SingIVness NonIV
-> ( C 2 ( ℝ 1 ) ( ℝ 2 ) )
-> ( C 2 ( ℝ 1 ) ( ℝ 1 ) )
-> ( C 2 ( ℝ 1 ) ( Spline Closed () ( ℝ 2 ) ) )
-> ( Seq ( C 2 ( ℝ 1 ) ( Corner ( ℝ 2 ) ) ) )
-> ( Maybe ( ℝ 1 -> Double ) )
-> ( ℝ 1 -> ( Seq ( CornerStrokeDatum 2 NonIV ), Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) ) )
#-}
{-# SPECIALISE INLINE brushStrokeData
:: SingIVness NonIV
-> ( C 2 ( ℝ 1 ) ( ℝ 2 ) )
-> ( C 2 ( ℝ 1 ) ( ℝ 2 ) )
-> ( C 2 ( ℝ 2 ) ( Spline Closed () ( ℝ 2 ) ) )
-> ( Seq ( C 2 ( ℝ 2 ) ( Corner ( ℝ 2 ) ) ) )
-> ( Maybe ( ℝ 2 -> Double ) )
-> ( ℝ 1 -> ( Seq ( CornerStrokeDatum 2 NonIV ), Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) ) )
#-}
{-# SPECIALISE INLINE brushStrokeData
:: SingIVness NonIV
-> ( C 2 ( ℝ 1 ) ( ℝ 2 ) )
-> ( C 2 ( ℝ 1 ) ( ℝ 3 ) )
-> ( C 2 ( ℝ 3 ) ( Spline Closed () ( ℝ 2 ) ) )
-> ( Seq ( C 2 ( ℝ 3 ) ( Corner ( ℝ 2 ) ) ) )
-> ( Maybe ( ℝ 3 -> Double ) )
-> ( ℝ 1 -> ( Seq ( CornerStrokeDatum 2 NonIV ), Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) ) )
#-}
{-# SPECIALISE INLINE brushStrokeData
:: SingIVness NonIV
-> ( C 2 ( ℝ 1 ) ( ℝ 2 ) )
-> ( C 2 ( ℝ 1 ) ( ℝ 4 ) )
-> ( C 2 ( ℝ 4 ) ( Spline Closed () ( ℝ 2 ) ) )
-> ( Seq ( C 2 ( ℝ 4 ) ( Corner ( ℝ 2 ) ) ) )
-> ( Maybe ( ℝ 4 -> Double ) )
-> ( ℝ 1 -> ( Seq ( CornerStrokeDatum 2 NonIV ), Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) ) )
#-}
{-# SPECIALISE INLINE brushStrokeData
:: SingIVness IsIV
-> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) )
-> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 1 ) )
-> ( C 3 ( 𝕀ℝ 1 ) ( Spline Closed () ( 𝕀ℝ 2 ) ) )
-> ( Seq ( C 3 ( 𝕀ℝ 1 ) ( Corner ( 𝕀ℝ 2 ) ) ) )
-> ( Maybe ( 𝕀ℝ 1 -> 𝕀 ) )
-> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-> ( 𝕀ℝ 1 -> ( Seq ( CornerStrokeDatum 3 IsIV ), Seq ( 𝕀ℝ 1 -> StrokeDatum 3 IsIV ) ) )
#-}
{-# SPECIALISE brushStrokeData
:: ( 𝕀 -> ( 𝕀ℝ 1 ) )
-> ( 𝕀ℝ 1 -> 𝕀 )
{-# SPECIALISE INLINE brushStrokeData
:: SingIVness IsIV
-> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) )
-> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) )
-> ( C 3 ( 𝕀ℝ 2 ) ( Spline Closed () ( 𝕀ℝ 2 ) ) )
-> ( Seq ( C 3 ( 𝕀ℝ 2 ) ( Corner ( 𝕀ℝ 2 ) ) ) )
-> ( Maybe ( 𝕀ℝ 2 -> 𝕀 ) )
-> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-> ( 𝕀ℝ 1 -> ( Seq ( CornerStrokeDatum 3 IsIV ), Seq ( 𝕀ℝ 1 -> StrokeDatum 3 IsIV ) ) )
#-}
{-# SPECIALISE brushStrokeData
:: ( 𝕀 -> ( 𝕀ℝ 1 ) )
-> ( 𝕀ℝ 1 -> 𝕀 )
{-# SPECIALISE INLINE brushStrokeData
:: SingIVness IsIV
-> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) )
-> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 3 ) )
-> ( C 3 ( 𝕀ℝ 3 ) ( Spline Closed () ( 𝕀ℝ 2 ) ) )
-> ( Seq ( C 3 ( 𝕀ℝ 3 ) ( Corner ( 𝕀ℝ 2 ) ) ) )
-> ( Maybe ( 𝕀ℝ 3 -> 𝕀 ) )
-> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-> ( 𝕀ℝ 1 -> ( Seq ( CornerStrokeDatum 3 IsIV ), Seq ( 𝕀ℝ 1 -> StrokeDatum 3 IsIV ) ) )
#-}
{-# SPECIALISE brushStrokeData
:: ( 𝕀 -> ( 𝕀ℝ 1 ) )
-> ( 𝕀ℝ 1 -> 𝕀 )
{-# SPECIALISE INLINE brushStrokeData
:: SingIVness IsIV
-> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) )
-> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 4 ) )
-> ( C 3 ( 𝕀ℝ 4 ) ( Spline Closed () ( 𝕀ℝ 2 ) ) )
-> ( Seq ( C 3 ( 𝕀ℝ 4 ) ( Corner ( 𝕀ℝ 2 ) ) ) )
-> ( Maybe ( 𝕀ℝ 4 -> 𝕀 ) )
-> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-> ( 𝕀ℝ 1 -> ( Seq ( CornerStrokeDatum 3 IsIV ), Seq ( 𝕀ℝ 1 -> StrokeDatum 3 IsIV ) ) )
#-}
-}
-- TODO: these specialisations fire in the benchmarking code because
-- we instantiate brushParams with ( ℝ nbBrushParams ), but they won't fire
-- in the main app code because we are using types such as "Record EllipseBrushFields".
--------------------------------------------------------------------------------
-- Solving the envelolpe equation: root-finding.
@ -1113,6 +1215,11 @@ data RootSolvingAlgorithm
-- | Use the modified Halley M2 method.
| HalleyM2 { maxIters :: Word, precision :: Int }
data FwdBwd
= Fwd
| Bwd
deriving stock ( Eq, Show )
-- | 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
@ -1120,39 +1227,44 @@ solveEnvelopeEquations :: RootSolvingAlgorithm
-> ℝ 2
-> T ( ℝ 2 )
-> ( Offset, Offset )
-> Seq ( ℝ 1 -> StrokeDatum 2 ℝ )
-> ( Seq ( CornerStrokeDatum 2 NonIV ), Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) )
-> ( ( ℝ 2, T ( ℝ 2 ) ), ( ℝ 2, T ( ℝ 2 ) ) )
solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffset ) strokeData
= ( fwdSol, ( bwdPt, -1 *^ bwdTgt ) )
solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffset ) ( corners, strokeData )
= ( fromMaybe fwdSol mbFwdCornerSol
, maybe ( bwdPt, -1 *^ bwdTgt ) ( \ ( pt, tgt ) -> ( pt, -1 *^ tgt ) ) mbBwdCornerSol
)
where
fwdSol = findSolFrom "fwd" fwdOffset
( bwdPt, bwdTgt ) = findSolFrom "bwd" bwdOffset
fwdSol = findSolFrom Fwd fwdOffset
( bwdPt, bwdTgt ) = findSolFrom Bwd bwdOffset
findSolFrom :: String -> Offset -> ( ℝ 2, T ( ℝ 2 ) )
findSolFrom desc ( Offset { offsetIndex = i00, offsetParameter = s00, offset = off } )
mbFwdCornerSol = getFirst $ foldMap ( First . cornerSol Fwd ) corners
mbBwdCornerSol = getFirst $ foldMap ( First . cornerSol Bwd ) corners
findSolFrom :: FwdBwd -> Offset -> ( ℝ 2, T ( ℝ 2 ) )
findSolFrom fwdOrBwd ( Offset { offsetIndex = i00, offsetParameter = s00, offset = off } )
= go ( fromIntegral i00 + fromMaybe 0.5 s00 )
where
go :: Double -> ( ℝ 2, T ( ℝ 2 ) )
go is0 =
case sol desc strokeData is0 of
case sol fwdOrBwd strokeData is0 of
( goodSoln, pt, tgt )
| goodSoln
-> ( pt, tgt )
| otherwise
-> ( off • path_t, path'_t )
sol :: String -> Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) -> Double -> ( Bool, ℝ 2, T ( ℝ 2 ) )
sol desc f is0 =
sol :: FwdBwd -> Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) -> Double -> ( Bool, ℝ 2, T ( ℝ 2 ) )
sol fwdOrBwd f is0 =
let ( solRes, _solSteps ) = runSolveMethod ( eqn f ) is0
( good, is ) =
case solRes of
Nothing -> ( False, is0 )
Just is1 -> ( if desc == "fwd"
then sgn >= 0
else sgn <= 0
Just is1 -> ( case fwdOrBwd of
Fwd -> sgn >= 0
Bwd -> sgn <= 0
, is1
)
( sgn, ds, dcdt ) = finish f is
@ -1164,7 +1276,7 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse
NewtonRaphson { maxIters, precision } ->
newtonRaphson maxIters precision domain
finish :: Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) -> Double -> ( Double, ℝ 2, T ( ℝ 2 ) )
finish :: Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) -> Double -> ( Double, ℝ 2, T ( ℝ 2 ) )
finish fs is =
let (i, s) = fromDomain is in
case evalStrokeDatum fs is of -- TODO: a bit redundant to have to compute this again...
@ -1192,21 +1304,67 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse
let cosθ = cos θ
sinθ = sin θ
in rotate cosθ sinθ $ unrot_dcdt
-- TODO: reduce duplication with 'cornerSol'.
in ( T u ^.^ T v, stroke, dcdt )
-- 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 ℝ ) -> ( Double -> StrokeDatum 2 ℝ )
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 1e-6 $ min (1 - 1e-6) $ s )
eqn :: Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) -> ( Double -> (# Double, Double #) )
eqn :: Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) -> ( Double -> (# Double, Double #) )
eqn fs is =
case evalStrokeDatum fs is of
StrokeDatum { ee = D12 ee _ ee_s } -> coerce (# ee, ee_s #)
cornerSol :: FwdBwd -> CornerStrokeDatum 2 NonIV -> Maybe ( ℝ 2, T ( ℝ 2 ) )
cornerSol fwdOrBwd
cornerDatum@(
CornerStrokeDatum
{ dpath = _dp@( D21 p _ _ )
, dbrush = D21 b _b_t _b_tt
, mbRotation
, corner = Corner _ startTgt endTgt }
)
=
let D11 u _u_t = cornerStrokeDatum cornerDatum
flipWhenBwd =
case fwdOrBwd of
Fwd -> id
Bwd -> ( -1 *^ )
b_s = flipWhenBwd ( T u )
ori = CCW --case fwdOrBwd of { Fwd -> CCW ; Bwd -> CW }
res@( _pt, _tgt ) = case mbRotation of
Nothing -> ( unT $ T p ^+^ T b
, T u )
Just ( D21 θ _ _ ) ->
let cosθ = cos θ
sinθ = sin θ
in ( unT $ T p ^+^ rotate cosθ sinθ ( T b )
, rotate cosθ sinθ ( T u ) )
in
{-
trace (unlines
[ "cornerSol"
, "fwdOrBwd: " ++ show fwdOrBwd
, "ori: " ++ show ori
, "dp: " ++ show _dp
, "startTgt: " ++ show startTgt
, "endTgt: " ++ show endTgt
, "u: " ++ show u
, "b_s: " ++ show b_s
, "ok: " ++ show ( between ori startTgt endTgt b_s )
, "pt: " ++ show _pt
, "tgt: " ++ show _tgt
]) $
-}
if isNothing $ between ori startTgt endTgt b_s
then Nothing
else Just res
n :: Int
n = Seq.length strokeData
domain :: ( Double, Double )
@ -1222,42 +1380,69 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse
-- | Computes the brush stroke coordinates of a cusp from
-- the @(t,s)@ parameter values.
--
-- TODO: use Newton's method starting at the midpoint of the box,
-- instead of just taking the midpoint of the box.
cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) )
cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) )
-> ( Int, Box N )
-> Cusp
-- TODO: use Newton's method starting at the midpoint of the box,
-- instead of just taking the midpoint of the box.
cuspCoords eqs ( i, box )
| StrokeDatum
{ dpath, stroke }
<- ( eqs ( ℝ1 t_mid ) `Seq.index` i ) ( ℝ1 s_mid )
= Cusp
{ cuspParameters = ℝ2 t_mid s_mid
, cuspPathCoords = coerce dpath
, cuspStrokeCoords = coerce stroke
, cuspPathCoords = dpath
, cuspStrokeCoords = stroke
}
where
𝕀 t_lo t_hi = box `index` Fin 1
𝕀 s_lo s_hi = box `index` Fin 2
t_mid = 0.5 * ( t_lo + t_hi )
s_mid = 0.5 * ( s_lo + s_hi )
t = box `index` Fin 1
s = box `index` Fin 2
t_mid = midpoint t
s_mid = midpoint s
cornerCuspCoords :: ( ℝ 1 -> Seq ( CornerStrokeDatum 2 NonIV ) )
-> ( Int, Box 1 )
-> Cusp
-- TODO: use Newton's method.
cornerCuspCoords eqs ( i, box )
| CornerStrokeDatum { dpath, stroke } <- eqs ( ℝ1 t_mid ) `Seq.index` i
= Cusp
{ cuspParameters = ℝ2 t_mid 0
, cuspPathCoords = dpath
, cuspStrokeCoords = stroke
}
where
t = box `index` Fin 1
t_mid = midpoint t
-- | 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 N 3
-> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], DoneBoxes N )
findCusps opts boxStrokeData =
findCuspsIn opts boxStrokeData $
IntMap.fromList
[ ( i, [ 𝕀ℝ2 ( 𝕀 zero one ) ( 𝕀 zero one ) ] )
| i <- [ 0 .. length ( boxStrokeData ( 𝕀ℝ1 ( 𝕀 zero one ) ) ) - 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 2, RootIsolationTree ( Box 2 ) ) ], DoneBoxes 2 )
)
findCusps ( opts1, opts2 ) boxStrokeData =
( findCornerCuspsIn opts1 ( fst . boxStrokeData )
( IntMap.fromList
[ ( i, [ 𝕀ℝ1 unit ] )
| i <- [ 0 .. length ( fst $ boxStrokeData ( 𝕀ℝ1 unit ) ) - 1 ]
]
)
, findCuspsIn opts2 ( snd . boxStrokeData )
( IntMap.fromList
[ ( i, [ 𝕀ℝ2 unit unit ] )
| i <- [ 0 .. length ( snd $ boxStrokeData ( 𝕀ℝ1 unit ) ) - 1 ]
]
)
)
where
unit :: 𝕀
unit = 𝕀 zero one
{-# INLINE unit #-}
zero, one :: Double
zero = 1e-6
one = 1 - zero
@ -1268,7 +1453,7 @@ findCusps opts boxStrokeData =
-- root isolation method.
findCuspsIn
:: RootIsolationOptions N 3
-> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 IsIV ) )
-> IntMap [ Box 2 ]
-> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], DoneBoxes N )
findCuspsIn opts boxStrokeData initBoxes =
@ -1289,95 +1474,58 @@ findCuspsIn opts boxStrokeData initBoxes =
( T $ 𝕀ℝ3 ee_t f_t g_t )
( T $ 𝕀ℝ3 ee_s f_s g_s )
{-
findCuspsIn
:: RootIsolationOptions N 3
-> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-> IntMap [ Box 2 ]
-> IntMap ( [ ( Box 3, RootIsolationTree ( Box 3 ) ) ], ( [ Box 3 ], [ Box 3 ] ) )
findCuspsIn opts boxStrokeData initBoxes =
IntMap.mapWithKey ( \ i boxes -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) $ concatMap ( mkInitBox i ) boxes ) initBoxes
-- | Like 'findCuspsIn' but in the case that the envelope is traced out
-- by a brush corner.
findCornerCuspsIn
:: RootIsolationOptions 1 2
-> ( 𝕀ℝ 1 -> Seq ( CornerStrokeDatum 3 IsIV ) )
-> IntMap [ Box 1 ]
-> IntMap ( [ ( Box 1, RootIsolationTree ( Box 1 ) ) ], DoneBoxes 1 )
findCornerCuspsIn opts boxStrokeData initBoxes =
IntMap.mapWithKey ( \ i -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) ) initBoxes
where
mkInitBox :: Int -> Box 2 -> [ Box 3 ]
mkInitBox i ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) =
let t = 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi )
s = 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi )
StrokeDatum
{ du =
D22
{ _D22_v = 𝕀 ( ℝ2 ux_lo uy_lo ) ( ℝ2 ux_hi uy_hi )}
, dv =
D22
{ _D22_v = 𝕀 ( ℝ2 vx_lo vy_lo ) ( ℝ2 vx_hi vy_hi ) }
, ee =
D22
{ _D22_dx = T ( 𝕀 ( ℝ1 ee_t_lo ) ( ℝ1 ee_t_hi ) )
, _D22_dy = T ( 𝕀 ( ℝ1 ee_s_lo ) ( ℝ1 ee_s_hi ) ) }
} = ( boxStrokeData t `Seq.index` i ) s
-- λ = ∂E/∂t / ∂E/∂s
λ1 = 𝕀 ee_t_lo ee_t_hi ⊘ 𝕀 ee_s_lo ee_s_hi
-- λ = u / v
λ2 = 𝕀 ux_lo ux_hi ⊘ 𝕀 vx_lo vx_hi
λ3 = 𝕀 uy_lo uy_hi ⊘ 𝕀 vy_lo vy_hi
λ = [ 𝕀 ( recip -0 ) ( recip 0 ) ]
`intersectMany` λ1
`intersectMany` λ2
`intersectMany` λ3
eqnPiece i ( 𝕀ℝ1 t ) =
let CornerStrokeDatum
{ dpath =
D31 { _D31_dx = p'
, _D31_dxdx = p''
}
, dbrush =
D31 { _D31_v = b
, _D31_dx = b'
, _D31_dxdx = b''
}
, mbRotation
} = ( boxStrokeData ( 𝕀ℝ1 t ) `Seq.index` i )
in case mbRotation of
Nothing ->
D11 ( unT $ p' ^+^ b' )
( p'' ^+^ b'' )
Just ( D31 { _D31_v = θ, _D31_dx = T θ', _D31_dxdx = T θ'' } ) ->
let cosθ = Ring.cos θ
sinθ = Ring.sin θ
swap ( T ( 𝕀ℝ2 x y ) ) = T ( 𝕀ℝ2 -y x )
rot = rotate cosθ sinθ
rot' = rot . swap
rot'' = rot . ( -1 *^ )
in
let boxes = [ 𝕀 ( ℝ3 t_lo s_lo λ_lo ) ( ℝ3 t_hi s_hi λ_hi )
| 𝕀 λ_lo' λ_hi' <- λ
, let λ_lo = max λ_lo' ( min ( λ_hi' - 10 ) -100 )
λ_hi = min λ_hi' ( max ( λ_lo' + 10 ) 100 ) ]
in boxes
eqnPiece :: Int -> 𝕀 ( ℝ 3 ) -> D1𝔸3 ( 𝕀 ( ℝ 3 ) )
eqnPiece i ( 𝕀 ( ℝ3 t_lo s_lo λ_lo ) ( ℝ3 t_hi s_hi λ_hi ) ) =
let t = 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi )
s = 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi )
λ = 𝕀 λ_lo λ_hi
-- TODO: I shouldn't have to do this by hand.
--
-- c(t) = p(t) + rot(θ(t)) b(t)
--
-- abbreviate; let R = rot(θ(t)), R'=rot'(θ(t)), ...
--
-- c = p + R b
-- c' = p' + R b' + θ' R' b
-- c'' = p'' + R b'' + 2 θ' R' b' + θ'² R'' b + θ'' R' b
--
D11 ( unT $ p' ^+^ rot b' ^+^ ( θ' *^ rot' ( T b ) ) )
( p'' ^+^ rot b''
^+^ ( 2 * θ' ) *^ rot' b'
^+^ ( θ' Ring.^ 2 ) *^ rot'' ( T b )
^+^ θ'' *^ rot' ( T b )
)
StrokeDatum
{ du =
D22 { _D22_v = u
, _D22_dx = T u_t , _D22_dy = T u_s }
, dv =
D22 { _D22_v = v
, _D22_dx = T v_t , _D22_dy = T v_s }
, ee =
D22 { _D22_dx = T ee_t, _D22_dy = T ee_s
, _D22_dxdx = T ee_tt, _D22_dxdy = T ee_ts, _D22_dydy = T ee_ss }
} = ( boxStrokeData t `Seq.index` i ) s
-- λ ∂E/∂s - ∂E/∂t = 0
𝕀 ( ℝ1 f1_lo ) ( ℝ1 f1_hi ) = unT $ λ *^ T ee_s ^-^ T ee_t
𝕀 ( ℝ1 f1_t_lo ) ( ℝ1 f1_t_hi ) = unT $ λ *^ T ee_ts ^-^ T ee_tt
𝕀 ( ℝ1 f1_s_lo ) ( ℝ1 f1_s_hi ) = unT $ λ *^ T ee_ss ^-^ T ee_ts
𝕀 ( ℝ1 f1_λ_lo ) ( ℝ1 f1_λ_hi ) = ee_s
-- λ v - u = 0
𝕀 ( ℝ2 f2_lo f3_lo ) ( ℝ2 f2_hi f3_hi ) = unT $ λ *^ T v ^-^ T u
𝕀 ( ℝ2 f2_t_lo f3_t_lo ) ( ℝ2 f2_t_hi f3_t_hi ) = unT $ λ *^ T v_t ^-^ T u_t
𝕀 ( ℝ2 f2_s_lo f3_s_lo ) ( ℝ2 f2_s_hi f3_s_hi ) = unT $ λ *^ T v_s ^-^ T u_s
𝕀 ( ℝ2 f2_λ_lo f3_λ_lo ) ( ℝ2 f2_λ_hi f3_λ_hi ) = v
in D13 ( 𝕀 ( ℝ3 f1_lo f2_lo f3_lo ) ( ℝ3 f1_hi f2_hi f3_hi ) )
( T $ 𝕀 ( ℝ3 f1_t_lo f2_t_lo f3_t_lo ) ( ℝ3 f1_t_hi f2_t_hi f3_t_hi ) )
( T $ 𝕀 ( ℝ3 f1_s_lo f2_s_lo f3_s_lo ) ( ℝ3 f1_s_hi f2_s_hi f3_s_hi ) )
( T $ 𝕀 ( ℝ3 f1_λ_lo f2_λ_lo f3_λ_lo ) ( ℝ3 f1_λ_hi f2_λ_hi f3_λ_hi ) )
intersectMany :: [𝕀] -> [𝕀] -> [𝕀]
intersectMany _ [] = []
intersectMany is (j : js) = intersectOne is j ++ intersectMany is js
intersectOne :: [ 𝕀 ] -> 𝕀 -> [ 𝕀 ]
intersectOne is i = concatMap ( intersect i ) is
intersect :: 𝕀 -> 𝕀 -> [ 𝕀 ]
intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 )
| lo > hi
= [ ]
| otherwise
= [ 𝕀 lo hi ]
where
lo = max lo1 lo2
hi = min hi1 hi2
-}
-- TODO: instead of c'(t) and c''(t), it might be simpler
-- to compute g(t) = rot(-θ(t)) c'(t) and g'(t)?

121
brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs
View file

@ -1,22 +1,25 @@
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE DuplicateRecordFields #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE RebindableSyntax #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE UndecidableInstances #-}
module Math.Bezier.Stroke.EnvelopeEquation
( StrokeDatum(..)
( StrokeDatum(..), CornerStrokeDatum(..)
, HasBézier(..)
, HasEnvelopeEquation(..)
, cornerStrokeDatum
) where
-- base
import Prelude hiding ( Num(..), (^), pi, sin, cos )
import Data.Coerce
( coerce )
import Data.Kind
( Type, Constraint )
import Data.List.NonEmpty
( NonEmpty(..) )
import Data.Proxy
( Proxy(..) )
import GHC.TypeNats
( Nat, type (-) )
@ -25,11 +28,13 @@ import Data.Act
( Torsor ((-->)) )
-- brush-strokes
import Calligraphy.Brushes
( Corner )
import Math.Algebra.Dual
import qualified Math.Bezier.Cubic as Cubic
import qualified Math.Bezier.Quadratic as Quadratic
import Math.Differentiable
( type I )
( I, IVness(..), SingIVness (..) )
import Math.Interval
import Math.Linear
import Math.Module
@ -41,7 +46,7 @@ import Math.Ring
-- | The value and derivative of a brush stroke at a given coordinate
-- \( (t, s) \), together with the value of the envelope equation at that
-- point.
type StrokeDatum :: Nat -> k -> Type
type StrokeDatum :: Nat -> IVness -> Type
data StrokeDatum k i
= StrokeDatum
{ -- | Path \( p(t) \).
@ -62,7 +67,11 @@ data StrokeDatum k i
-- | Envelope function
--
-- \[ E(t_0,s_0) = \left ( \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \right )_{(t_0,s_0)}. \]
-- \[ \begin{align} E(t,s) &= \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \\
-- &= \left ( R(-\theta(t)) \frac{\partial c}{\partial t} \right ) \times \left ( R(-\theta(t)) \frac{\partial c}{\partial s} \right ) \\
-- &= \left ( R(-\theta(t)) p'(t) + \theta'(t) S b(t,s) + \frac{\partial b}{\partial t} \right ) \times \frac{\partial b}{\partial s} \]
--
-- where \( S \) denotes the 2D swap function \( S(x,y) = (y,x) \).
, ee :: D ( k - 1 ) 2 ( I i 1 )
-- \[ \frac{\partial E}{\partial s} \frac{\mathrm{d} c}{\mathrm{d} t}, \]
@ -73,12 +82,43 @@ data StrokeDatum k i
, 𝛿E𝛿sdcdt :: D ( k - 2 ) 2 ( T ( I i 2 ) )
}
deriving stock instance Show ( StrokeDatum 2 ℝ )
deriving stock instance Show ( StrokeDatum 3 𝕀 )
deriving stock instance
( Show ( I i 2 )
, Show ( D k 1 ( I i Double ) )
, Show ( D k 1 ( I i 2 ) )
, Show ( D k 2 ( I i 2 ) )
, Show ( D ( k - 1 ) 2 ( I i 1 ) )
, Show ( D ( k - 1 ) 2 ( I i 2 ) )
, Show ( D ( k - 2 ) 2 ( T ( I i 2 ) ) )
)
=> Show ( StrokeDatum k i )
-- | A slimmed down version of 'StrokeDatum' for use at brush corners.
type CornerStrokeDatum :: Nat -> IVness -> Type
data CornerStrokeDatum k i
= CornerStrokeDatum
{ -- | Path \( p(t) \).
dpath :: D k 1 ( I i 2 )
-- | Brush shape \( b(t) \).
, dbrush :: D k 1 ( I i 2 )
-- | (Optional) rotation angle \( \theta(t) \).
, mbRotation :: Maybe ( D k 1 ( I i Double ) )
, corner :: Corner ( I i 2 )
, stroke :: I i 2
}
deriving stock instance
( Show ( I i Double )
, Show ( I i 2 )
, Show ( D k 1 ( I i Double ) )
, Show ( D k 1 ( I i 2 ) )
)
=> Show ( CornerStrokeDatum k i )
--------------------------------------------------------------------------------
type HasBézier :: forall {t}. Nat -> t -> Constraint
type HasBézier :: Nat -> IVness -> Constraint
class HasBézier k i where
-- | Linear interpolation, as a differentiable function.
@ -112,8 +152,7 @@ class HasBézier k i where
type HasEnvelopeEquation :: Nat -> Constraint
class HasEnvelopeEquation k where
uncurryD :: Dim a ~ 1
=> D k 1 ( C k a b ) -> a -> D k 2 b
uncurryD :: D k 1 ( D k 1 b ) -> D k 2 b
-- | The envelope function
--
@ -135,14 +174,13 @@ class HasEnvelopeEquation k where
, Show ( StrokeDatum k i )
)
=> Proxy i
-> ( I i Double -> I i 1 )
=> SingIVness i
-> D k 1 ( I i 2 )
-> D k 2 ( I i 2 )
-> Maybe ( D k 1 ( I i Double ) )
-> StrokeDatum k i
instance HasBézier 2 ℝ where
instance HasBézier 2 NonIV where
line co ( Segment a b :: Segment b ) =
D \ ( co -> t ) ->
D21 ( lerp @( T b ) t a b )
@ -169,7 +207,7 @@ instance HasEnvelopeEquation 2 where
uncurryD = uncurryD2
{-# INLINEABLE uncurryD #-}
envelopeEquation ( _ :: Proxy i ) co
envelopeEquation ( ivness :: SingIVness i )
dp@( D21 ( T -> p ) p_t p_tt )
db@( D22 ( T -> b ) b_t b_s
b_tt b_ts b_ss )
@ -182,6 +220,10 @@ instance HasEnvelopeEquation 2 where
, du, dv, ee, 𝛿E𝛿sdcdt
}
where
co :: I i Double -> I i 1
co = case ivness of
SNonIV -> coerce
SIsIV -> coerce
(ee, 𝛿E𝛿sdcdt) =
let
D12 (T -> u) u_t u_s = du
@ -250,7 +292,46 @@ instance HasEnvelopeEquation 2 where
)
{-# INLINEABLE envelopeEquation #-}
instance HasBézier 3 ℝ where
-- | Computes \( u \) and \( u_t \), where \( u \) is as in 'StrokeDatum'.
cornerStrokeDatum
:: CornerStrokeDatum 2 NonIV
-> D 1 1 ( ℝ 2 )
cornerStrokeDatum
CornerStrokeDatum
{ dpath = D21 _ p_t p_tt
, dbrush = D21 ( T -> b ) b_t b_tt
, mbRotation } =
case mbRotation of
Nothing -> D11 ( unT $ p_t ^+^ b_t ) ( p_tt ^+^ b_tt )
Just ( D21 θ ( T θ_t ) ( T θ_tt ) ) ->
let
rot, rot' :: T ( ℝ 2 ) -> T ( ℝ 2 )
cosθ = cos θ
sinθ = sin θ
rot = rotate cosθ -sinθ
rot' = rotate sinθ cosθ
swap :: T ( ℝ 2 ) -> T ( ℝ 2 )
swap ( T xy ) =
let x = xy `index` Fin 1
y = xy `index` Fin 2
in T $ tabulate \ case
Fin 1 -> -y
_ -> x
u, u_t :: T ( ℝ 2 )
u = rot p_t ^+^ θ_t *^ swap b ^+^ b_t
u_t = ( -θ_t *^ rot' p_t
^+^ rot p_tt
)
^+^
( θ_tt *^ swap b
^+^ θ_t *^ swap b_t
)
^+^ b_tt
in D11 ( unT u ) u_t
instance HasBézier 3 NonIV where
line co ( Segment a b :: Segment b ) =
D \ ( co -> t ) ->
@ -281,7 +362,7 @@ instance HasEnvelopeEquation 3 where
uncurryD = uncurryD3
{-# INLINEABLE uncurryD #-}
envelopeEquation ( _ :: Proxy i ) co
envelopeEquation ( ivness :: SingIVness i )
dp@( D31 ( T -> p ) p_t p_tt p_ttt )
db@( D32 ( T -> b ) b_t b_s
b_tt b_ts b_ss
@ -295,6 +376,10 @@ instance HasEnvelopeEquation 3 where
, du, dv, ee, 𝛿E𝛿sdcdt
}
where
co :: I i Double -> I i 1
co = case ivness of
SNonIV -> coerce
SIsIV -> coerce
(ee, 𝛿E𝛿sdcdt) =
let
D22 (T -> u) u_t u_s u_tt u_ts u_ss = du
@ -400,7 +485,7 @@ instance HasEnvelopeEquation 3 where
)
{-# INLINEABLE envelopeEquation #-}
instance HasBézier 3 𝕀 where
instance HasBézier 3 IsIV where
line co ( Segment a b :: Segment b ) =
D \ ( co -> t ) ->

29
brush-strokes/src/lib/Math/Differentiable.hs
View file

@ -2,7 +2,9 @@
{-# LANGUAGE UndecidableInstances #-}
module Math.Differentiable
( I, Differentiable, DiffInterp )
( IVness(..), SingIVness(..)
, I, Differentiable, DiffInterp
)
where
-- base
@ -23,17 +25,24 @@ import Math.Ring
--------------------------------------------------------------------------------
data IVness
= NonIV
| IsIV
type SingIVness :: IVness -> Type
data SingIVness iv where
SNonIV :: SingIVness NonIV
SIsIV :: SingIVness IsIV
-- | Type family to parametrise over both pointwise and interval computations.
--
-- Use 'ℝ' parameter for points, and '𝕀' parameter for intervals.
type I :: k -> l -> Type
type I :: forall l. IVness -> l -> Type
type family I i t
type instance I @(Nat -> Type) @Type ℝ Double = Double
type instance I @Type @Type 𝕀 Double = 𝕀
type instance I @(Nat -> Type) @Nat ℝ n = ℝ n
type instance I @Type @Nat 𝕀 n = 𝕀ℝ n
type instance I @Type NonIV Double = Double
type instance I @Type IsIV Double = 𝕀
type instance I @Nat NonIV n = ℝ n
type instance I @Nat IsIV n = 𝕀ℝ n
type Differentiable :: Nat -> k -> Nat -> Constraint
type Differentiable :: Nat -> IVness -> Nat -> Constraint
class
( Module ( I i Double ) ( T ( I i u ) )
, HasChainRule ( I i Double ) k ( I i u )
@ -47,7 +56,7 @@ instance
, Representable ( I i Double ) ( I i u )
) => Differentiable k i u
type DiffInterp :: Nat -> k -> Nat -> Constraint
type DiffInterp :: Nat -> IVness -> Nat -> Constraint
class
( Differentiable k i u
, Interpolatable ( I i Double ) ( I i u )

60
brush-strokes/src/lib/Math/Interval.hs
View file

@ -242,26 +242,15 @@ deriving via ( ViaModule 𝕀 ( 𝕀ℝ n ) )
-- HasChainRule instances.
instance HasChainRule 𝕀 2 ( 𝕀ℝ 0 ) where
konst w = D0 w
linearD f v = D0 ( f v )
value ( D0 v ) = v
chain _ ( D0 gfx ) = D21 gfx origin origin
instance HasChainRule 𝕀 3 ( 𝕀ℝ 0 ) where
konst w = D0 w
linearD f v = D0 ( f v )
value ( D0 v ) = v
chain _ ( D0 gfx ) = D31 gfx origin origin origin
instance HasChainRule 𝕀 2 ( 𝕀ℝ 1 ) where
konst :: forall w. AbelianGroup w => w -> D2𝔸1 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 1 -> w ) -> 𝕀ℝ 1 -> D2𝔸1 w
linearD f v =
let !o = origin @𝕀 @( T w )
@ -290,13 +279,6 @@ instance HasChainRule 𝕀 2 ( 𝕀ℝ 1 ) where
instance HasChainRule 𝕀 3 ( 𝕀ℝ 1 ) where
konst :: forall w. AbelianGroup w => w -> D3𝔸1 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 1 -> w ) -> 𝕀ℝ 1 -> D3𝔸1 w
linearD f v =
let !o = origin @𝕀 @( T w )
@ -325,13 +307,6 @@ instance HasChainRule 𝕀 3 ( 𝕀ℝ 1 ) where
instance HasChainRule 𝕀 2 ( 𝕀ℝ 2 ) where
konst :: forall w. AbelianGroup w => w -> D2𝔸2 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 2 -> w ) -> 𝕀ℝ 2 -> D2𝔸2 w
linearD f v =
let !o = origin @𝕀 @( T w )
@ -360,13 +335,6 @@ instance HasChainRule 𝕀 2 ( 𝕀ℝ 2 ) where
instance HasChainRule 𝕀 3 ( 𝕀ℝ 2 ) where
konst :: forall w. AbelianGroup w => w -> D3𝔸2 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 2 -> w ) -> 𝕀ℝ 2 -> D3𝔸2 w
linearD f v =
let !o = origin @𝕀 @( T w )
@ -395,13 +363,6 @@ instance HasChainRule 𝕀 3 ( 𝕀ℝ 2 ) where
instance HasChainRule 𝕀 2 ( 𝕀ℝ 3 ) where
konst :: forall w. AbelianGroup w => w -> D2𝔸3 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 3 -> w ) -> 𝕀ℝ 3 -> D2𝔸3 w
linearD f v =
let !o = origin @𝕀 @( T w )
@ -430,13 +391,6 @@ instance HasChainRule 𝕀 2 ( 𝕀ℝ 3 ) where
instance HasChainRule 𝕀 3 ( 𝕀ℝ 3 ) where
konst :: forall w. AbelianGroup w => w -> D3𝔸3 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 3 -> w ) -> 𝕀ℝ 3 -> D3𝔸3 w
linearD f v =
let !o = origin @𝕀 @( T w )
@ -465,13 +419,6 @@ instance HasChainRule 𝕀 3 ( 𝕀ℝ 3 ) where
instance HasChainRule 𝕀 2 ( 𝕀ℝ 4 ) where
konst :: forall w. AbelianGroup w => w -> D2𝔸4 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 4 -> w ) -> 𝕀ℝ 4 -> D2𝔸4 w
linearD f v =
let !o = origin @𝕀 @( T w )
@ -500,13 +447,6 @@ instance HasChainRule 𝕀 2 ( 𝕀ℝ 4 ) where
instance HasChainRule 𝕀 3 ( 𝕀ℝ 4 ) where
konst :: forall w. AbelianGroup w => w -> D3𝔸4 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 4 -> w ) -> 𝕀ℝ 4 -> D3𝔸4 w
linearD f v =
let !o = origin @𝕀 @( T w )

5
src/app/MetaBrush/Application.hs
View file

@ -232,7 +232,10 @@ runApplication application = do
NewtonRaphson
{ maxIters = 20, precision = 8 }
cuspFindingOptionsTVar <- STM.newTVarIO $
Just defaultRootIsolationOptions
Just
( defaultRootIsolationOptions @1 @2
, defaultRootIsolationOptions @2 @3
)
-- Put all these stateful variables in a record for conciseness.
let

2
src/app/MetaBrush/Application/Context.hs
View file

@ -127,7 +127,7 @@ data Variables
, maxHistorySizeTVar :: !( STM.TVar Int )
, fitParametersTVar :: !( STM.TVar FitParameters )
, rootsAlgoTVar :: !( STM.TVar RootSolvingAlgorithm )
, cuspFindingOptionsTVar :: !( STM.TVar ( Maybe ( RootIsolationOptions 2 3 ) ) )
, cuspFindingOptionsTVar :: !( STM.TVar ( Maybe ( RootIsolationOptions 1 2, RootIsolationOptions 2 3 ) ) )
}
--------------------------------------------------------------------------------

20
src/app/MetaBrush/Render/Document.hs
View file

@ -148,7 +148,9 @@ blankRender _ = pure ()
getDocumentRender
:: Colours
-> RootSolvingAlgorithm -> Maybe ( RootIsolationOptions 2 3 ) -> FitParameters
-> RootSolvingAlgorithm
-> Maybe ( RootIsolationOptions 1 2, RootIsolationOptions 2 3 )
-> FitParameters
-> Mode -> Bool
-> Set Modifier -> Maybe ( ℝ 2 ) -> Maybe HoldAction -> Maybe PartialPath
-> Document
@ -300,7 +302,7 @@ instance NFData StrokeRenderData where
-- - Otherwise, this consists of the underlying spline path only.
strokeRenderData
:: RootSolvingAlgorithm
-> Maybe ( RootIsolationOptions 2 3 )
-> Maybe ( RootIsolationOptions 1 2, RootIsolationOptions 2 3 )
-> FitParameters
-> Stroke
-> ST RealWorld StrokeRenderData
@ -325,6 +327,7 @@ strokeRenderData rootAlgo mbCuspOptions fitParams
-> do
let embedUsedParams = inject2 $ MkR brush_defaults
-- Compute the outline using the brush function.
( outline, fitPts, cusps ) <-
computeStrokeOutline @clo rootAlgo mbCuspOptions fitParams
@ -337,7 +340,7 @@ strokeRenderData rootAlgo mbCuspOptions fitParams
, strokeBrushFunction =
\ params ->
let MkR brushParams = embedUsedParams $ toUsedParams params
shape = fun @Double brushBaseShape brushParams
shape = fun brushBaseShape brushParams
-- TODO: remove this logic which is duplicated
-- from elsewhere. The type should make it
-- impossible to forget to apply the rotation.
@ -703,6 +706,17 @@ drawFitPoint _ ( Zoom { zoomFactor } ) ( FitPoint { fitPoint = ℝ2 x y } ) = do
Cairo.fill
Cairo.restore
drawFitPoint _ ( Zoom { zoomFactor } ) ( JoinPoint { joinPoint = ℝ2 x y } ) = lift do
-- Draw a little white circle.
Cairo.save
Cairo.translate x y
Cairo.arc 0 0 ( 3 / zoomFactor ) 0 ( 2 * pi )
Cairo.setSourceRGBA 255 255 255 0.9
Cairo.setLineWidth ( 1 / zoomFactor )
Cairo.stroke
Cairo.restore
drawFitPoint _ ( Zoom { zoomFactor } ) ( FitTangent { fitPoint = ℝ2 x y, fitTangent = V2 tx ty } ) = do
hue <- get

34
src/metabrushes/MetaBrush/Brush.hs
View file

@ -44,10 +44,12 @@ import qualified Data.Text as Text
-- brush-strokes
import Calligraphy.Brushes
( Brush(..) )
import Math.Algebra.Dual
( HasChainRule )
import Math.Differentiable
( DiffInterp )
( DiffInterp, IVness(..) )
import Math.Interval
( 𝕀 )
( 𝕀, 𝕀ℝ )
import Math.Linear
-- MetaBrush
@ -78,13 +80,19 @@ type BrushFunction brushFields =
type NamedBrush :: [ Symbol ] -> Type
data NamedBrush brushFields where
NamedBrush
:: forall brushFields
. ( KnownSymbols brushFields
, Representable Double ( ℝ ( Length brushFields ) )
, DiffInterp 2 ℝ ( Length brushFields )
, DiffInterp 3 𝕀 ( Length brushFields )
, Show ( ℝ ( Length brushFields ) )
, NFData ( ℝ ( Length brushFields ) )
:: forall brushFields nbBrushFields
. ( nbBrushFields ~ Length brushFields
, KnownSymbols brushFields
, Representable Double ( ℝ nbBrushFields )
, DiffInterp 2 NonIV nbBrushFields
, DiffInterp 3 IsIV nbBrushFields
-- TODO: clean up this constraint?
, HasChainRule 𝕀 2 ( 𝕀ℝ nbBrushFields )
, Show ( ℝ nbBrushFields )
, NFData ( ℝ nbBrushFields )
)
=> { brushName :: !Text
, brushFunction :: !( BrushFunction brushFields )
@ -122,8 +130,8 @@ class ( KnownSymbols pointFields, Typeable pointFields
, NFData ( Record pointFields )
, Representable Double ( ℝ ( Length pointFields ) )
, RepDim ( ℝ ( Length pointFields ) ) ~ Length pointFields
, DiffInterp 2 ℝ ( Length pointFields )
, DiffInterp 3 𝕀 ( Length pointFields )
, DiffInterp 2 NonIV ( Length pointFields )
, DiffInterp 3 IsIV ( Length pointFields )
)
=> PointFields pointFields where { }
instance ( KnownSymbols pointFields, Typeable pointFields
@ -133,8 +141,8 @@ instance ( KnownSymbols pointFields, Typeable pointFields
, NFData ( Record pointFields )
, Representable Double ( ℝ ( Length pointFields ) )
, RepDim ( ℝ ( Length pointFields ) ) ~ Length pointFields
, DiffInterp 2 ℝ ( Length pointFields )
, DiffInterp 3 𝕀 ( Length pointFields )
, DiffInterp 2 NonIV ( Length pointFields )
, DiffInterp 3 IsIV ( Length pointFields )
)
=> PointFields pointFields where { }

42
src/metabrushes/MetaBrush/Records.hs
View file

@ -55,7 +55,6 @@ import qualified Data.Text as Text
-- MetaBrush
import Math.Differentiable
import Math.Interval
import Math.Linear
import Math.Module
( Module )
@ -146,8 +145,9 @@ intersect :: forall r1 r2 l1 l2
, Representable Double ( ℝ l2 )
, Show ( ℝ l1 )
, Show ( ℝ l2 )
, Differentiable 2 ℝ l2
, Differentiable 3 𝕀 l2
, Differentiable 2 NonIV l2
, Differentiable 2 IsIV l2
, Differentiable 3 IsIV l2
)
=> Intersection r1 r2
intersect
@ -178,8 +178,9 @@ data Intersection r1 r2 where
, KnownSymbols r1r2
, Representable Double ( ℝ l12 )
, Show ( ℝ l12 )
, Differentiable 2 ℝ l12
, Differentiable 3 𝕀 l12
, Differentiable 2 NonIV l12
, Differentiable 2 IsIV l12
, Differentiable 3 IsIV l12
)
=> { project1 :: Record r1 -> Record r1r2
-- ^ project out fields present in both rows
@ -204,8 +205,9 @@ doIntersection
=> ( forall r1r2 l12.
( r1r2 ~ Intersect r1 r2
, KnownSymbols r1r2, l12 ~ Length r1r2
, Differentiable 2 ℝ l12
, Differentiable 3 𝕀 l12
, Differentiable 2 NonIV l12
, Differentiable 2 IsIV l12
, Differentiable 3 IsIV l12
, Representable Double ( ℝ l12 )
, Show ( ℝ l12 )
)
@ -304,17 +306,17 @@ type family Elem k ks where
data Union r1 r2 where
Union
:: forall r1r2_rocky r1 r2 l12_rocky
. ( l12_rocky ~ Length r1r2_rocky
, KnownSymbols r1r2_rocky
, Representable Double ( ℝ l12_rocky )
, Show ( ℝ l12_rocky )
, NFData ( ℝ l12_rocky )
, DiffInterp 2 ℝ l12_rocky
, DiffInterp 3 𝕀 l12_rocky
:: forall r1r2 r1 r2 l12
. ( l12 ~ Length r1r2
, KnownSymbols r1r2
, Representable Double ( ℝ l12 )
, Show ( ℝ l12 )
, NFData ( ℝ l12 )
, DiffInterp 2 NonIV l12
, DiffInterp 3 IsIV l12
)
=> { unionWith :: ( Double -> Double -> Double )
-> Record r1 -> Record r2 -> Record r1r2_rocky
-> Record r1 -> Record r2 -> Record r1r2
-- ^ union of two records
} -> Union r1 r2
@ -327,8 +329,8 @@ union :: forall r1 r2 l1 l2
, Show ( ℝ l1 )
, Show ( ℝ l2 )
, NFData ( ℝ l2 )
, DiffInterp 2 ℝ l2
, DiffInterp 3 𝕀 l2
, DiffInterp 2 NonIV l2
, DiffInterp 3 IsIV l2
)
=> Union r1 r2
union
@ -380,8 +382,8 @@ doUnion
)
=> ( forall r1r2 l12.
( KnownSymbols r1r2, l12 ~ Length r1r2
, DiffInterp 2 ℝ l12
, DiffInterp 3 𝕀 l12
, DiffInterp 2 NonIV l12
, DiffInterp 3 IsIV l12
, Representable Double ( ℝ l12 )
, Show ( ℝ l12 )
, NFData ( ℝ l12 )