continue interval arithmetic integration

This commit is contained in:
sheaf 2023-01-13 15:23:33 +01:00
parent a4e9c1cf32
commit 09c1bdd948
5 changed files with 115 additions and 77 deletions

View file

@ -76,6 +76,7 @@ common common
PatternSynonyms
RankNTypes
RecordWildCards
RoleAnnotations
StandaloneDeriving
StandaloneKindSignatures
TupleSections

View file

@ -18,7 +18,7 @@ import Data.Functor.Const
import Data.Functor.Identity
( Identity(..) )
import Data.Kind
( Constraint )
( Type, Constraint )
import Data.Monoid
( Ap(..) )
import Data.Semigroup
@ -87,7 +87,8 @@ instance SplineTypeI Open where
instance SplineTypeI Closed where
ssplineType = SClosed
data family NextPoint ( clo :: SplineType ) ptData
type NextPoint :: SplineType -> Type -> Type
data family NextPoint clo
newtype instance NextPoint Open ptData = NextPoint { nextPoint :: ptData }
deriving stock ( Show, Generic, Generic1, Functor, Foldable, Traversable )
deriving anyclass ( NFData, NFData1 )
@ -111,7 +112,8 @@ toNextPoint pt = case ssplineType @clo of
SOpen -> NextPoint pt
SClosed -> BackToStart
data Curve ( clo :: SplineType ) crvData ptData
type Curve :: SplineType -> Type -> Type -> Type
data Curve clo crvData ptData
= LineTo
{ curveEnd :: !( NextPoint clo ptData )
, curveData :: !crvData
@ -157,7 +159,8 @@ openCurveStart ( LineTo ( NextPoint p ) _ ) = p
openCurveStart ( Bezier2To p _ _ ) = p
openCurveStart ( Bezier3To p _ _ _ ) = p
data family Curves ( clo :: SplineType ) crvData ptData
type Curves :: SplineType -> Type -> Type -> Type
data family Curves clo
newtype instance Curves Open crvData ptData
= OpenCurves { openCurves :: Seq ( Curve Open crvData ptData ) }

View file

@ -1,8 +1,7 @@
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE QuantifiedConstraints #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE DuplicateRecordFields #-}
{-# LANGUAGE UndecidableInstances #-}
module Math.Bezier.Stroke
( Offset(..)
@ -32,11 +31,13 @@ import Control.Monad.ST
import Data.Bifunctor
( Bifunctor(bimap) )
import Data.Coerce
( coerce )
( Coercible, coerce )
import Data.Foldable
( for_ )
import Data.Functor.Identity
( Identity(..) )
import Data.Kind
( Type )
import Data.List.NonEmpty
( unzip )
import Data.Maybe
@ -434,13 +435,16 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
case crv of
LineTo { curveEnd = NextPoint sp1 }
| let seg = Segment sp0 sp1
-> ( line ( fmap coords seg ), line ( fmap ptParams seg ) )
-> ( line @Point ( fmap coords seg )
, line @Point ( fmap ptParams seg ) )
Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 }
| let bez2 = Quadratic.Bezier sp0 sp1 sp2
-> ( bezier2 ( fmap coords bez2 ), bezier2 ( fmap ptParams bez2 ) )
-> ( bezier2 @Point ( fmap coords bez2 )
, bezier2 @Point ( fmap ptParams bez2 ) )
Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 }
| let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3
-> ( bezier3 ( fmap coords bez3 ), bezier3 ( fmap ptParams bez3 ) )
-> ( bezier3 @Point ( fmap coords bez3 )
, bezier3 @Point ( fmap ptParams bez3 ) )
fwdBwd :: OutlineFn
fwdBwd t
@ -449,8 +453,12 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
-- , ( offset bwdOffset • path_t, -1 *^ path'_t ) )
where
curves :: Seq ( 1 -> StrokeDatum )
curves = brushStrokeData path ( usedParams `chainRule` toBrushParams ) brushFromParams t
curves :: Seq ( 1 -> StrokeDatum Point )
curves = brushStrokeData @Point @brushParams
path
( usedParams `chainRule` toBrushParams )
brushFromParams
t
fwdOffset = withTangent path'_t brush_t
bwdOffset = withTangent ( -1 *^ path'_t ) brush_t
@ -771,7 +779,7 @@ brushStroke ( D1 p dpdt d2pdt2 ) ( D2 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) =
-- | The envelope equation
--
-- \[ E = \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} = 0, ]
-- \[ E = \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} = 0, \]
--
-- together with the total derivative
--
@ -782,12 +790,13 @@ brushStroke ( D1 p dpdt d2pdt2 ) ( D2 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) =
-- \[ \frac{\partial E}{\partial s}, \qquad \frac{\partial E}{\partial s}. \]
--
-- NB: if \( \frac{\partial E}{\partial s} \) is zero, the total derivative is ill-defined.
envelopeEquation :: ( D ( i ( 2 ) ) ~ D ( 2 )
, Cross ( i Double ) ( T ( i ( 2 ) ) )
, Fractional ( i Double )
envelopeEquation :: forall i
. ( D ( I i ( 2 ) ) ~ D ( 2 )
, Cross ( I i Double ) ( T ( I i ( 2 ) ) )
, Fractional ( I i Double )
)
=> D ( i ( 2 ) ) ( i ( 2 ) )
-> ( i Double, T ( i ( 2 ) ), i Double, i Double )
=> D ( I i ( 2 ) ) ( I i ( 2 ) )
-> ( I i Double, T ( I i ( 2 ) ), I i Double, I i Double )
envelopeEquation ( D2 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ) =
let ee = dcdt `cross` dcds
dEdt = d2cdt2 `cross` dcds + dcdt `cross` d2cdtds
@ -802,30 +811,47 @@ envelopeEquation ( D2 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ) =
-- ∂E/∂s dc/dt = ∂E/∂s ∂c/∂t - ∂E/∂t ∂c/∂s.
-- | Linear interpolation, as a differentiable function.
line :: forall b. ( Module Double ( T b ), Torsor ( T b ) b )
=> Segment b -> 1 ~> b
line ( Segment a b ) = D \ ( 1 t ) ->
line :: forall i b
. ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
, D ( I i ( 1 ) ) ~ D ( 1 )
, Coercible ( I i ( 1 ) ) ( I i Double )
)
=> Segment b -> I i ( 1 ) ~> b
line ( Segment a b ) = D \ ( coerce -> t ) ->
D1 ( lerp @( T b ) t a b )
( a --> b )
origin
-- | A quadratic Bézier curve, as a differentiable function.
bezier2 :: forall b. ( Module Double ( T b ), Torsor ( T b ) b )
=> Quadratic.Bezier b -> 1 ~> b
bezier2 bez = D \ ( 1 t ) ->
bezier2 :: forall i b
. ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
, D ( I i ( 1 ) ) ~ D ( 1 )
, Coercible ( I i ( 1 ) ) ( I i Double )
)
=> Quadratic.Bezier b -> I i ( 1 ) ~> b
bezier2 bez = D \ ( coerce -> t ) ->
D1 ( Quadratic.bezier @( T b ) bez t )
( Quadratic.bezier' bez t )
( Quadratic.bezier'' bez )
-- | A cubic Bézier curve, as a differentiable function.
bezier3 :: forall b. ( Module Double ( T b ), Torsor ( T b ) b )
=> Cubic.Bezier b -> 1 ~> b
bezier3 bez = D \ ( 1 t ) ->
bezier3 :: forall i b
. ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
, D ( I i ( 1 ) ) ~ D ( 1 )
, Coercible ( I i ( 1 ) ) ( I i Double )
)
=> Cubic.Bezier b -> I i ( 1 ) ~> b
bezier3 bez = D \ ( coerce -> t ) ->
D1 ( Cubic.bezier @( T b ) bez t )
( Cubic.bezier' bez t )
( Cubic.bezier'' bez t )
splineCurveFns :: SplinePts Closed -> Seq ( 1 ~> 2 )
splineCurveFns :: forall i
. ( D ( I i ( 1 ) ) ~ D ( 1 )
, Module ( I i Double ) ( T ( I i ( 2 ) ) )
, Torsor ( T ( I i ( 2 ) ) ) ( I i ( 2 ) )
, Coercible ( I i ( 1 ) ) ( I i Double ) )
=> Spline Closed () ( I i ( 2 ) ) -> Seq ( I i ( 1 ) ~> I i ( 2 ) )
splineCurveFns spls
= runIdentity
. bifoldSpline
@ -834,21 +860,23 @@ splineCurveFns spls
. adjustSplineType @Open
$ spls
where
curveFn :: 2 -> Curve Open () ( 2 ) -> ( 1 ~> 2 )
curveFn :: I i ( 2 )
-> Curve Open () ( I i ( 2 ) )
-> ( I i ( 1 ) ~> I i ( 2 ) )
curveFn p0 = \case
LineTo { curveEnd = NextPoint p1 }
-> line $ Segment p0 p1
-> line @i $ Segment p0 p1
Bezier2To { controlPoint = p1, curveEnd = NextPoint p2 }
-> bezier2 $ Quadratic.Bezier p0 p1 p2
-> bezier2 @i $ Quadratic.Bezier p0 p1 p2
Bezier3To { controlPoint1 = p1, controlPoint2 = p2, curveEnd = NextPoint p3 }
-> bezier3 $ Cubic.Bezier p0 p1 p2 p3
-> bezier3 @i $ Cubic.Bezier p0 p1 p2 p3
-- | 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 :: 2
-> T ( 2 )
-> ( Offset, Offset )
-> Seq ( 1 -> StrokeDatum )
-> Seq ( 1 -> StrokeDatum Point )
-> ( ( 2, T ( 2 ) ), ( 2, T ( 2 ) ) )
solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData
= ( fwdSol, ( bwdPt, -1 *^ bwdTgt ) )
@ -909,7 +937,7 @@ solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData
| otherwise
= i + 1
sol :: Double -> ( 1 -> StrokeDatum ) -> ( Bool, 1, 2, T ( 2 ) )
sol :: Double -> ( 1 -> StrokeDatum Point ) -> ( Bool, 1, 2, T ( 2 ) )
sol initialGuess f =
let (good, s) = case newtonRaphson maxIters precision domain ( eqn f ) initialGuess of
Nothing -> ( False, initialGuess )
@ -937,7 +965,7 @@ solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData
-- ] )
( good, 1 s, value @Double @( 2 ) dstroke, totDeriv )
eqn :: ( 1 -> StrokeDatum ) -> ( Double -> ( Double, Double ) )
eqn :: ( 1 -> StrokeDatum Point ) -> ( Double -> ( Double, Double ) )
eqn f s =
case f ( 1 s ) of
StrokeDatum { ee, 𝛿E𝛿s } ->
@ -956,36 +984,47 @@ 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 )
brushStrokeData :: forall brushParams
. ( Diffy Double brushParams, Show brushParams )
=> ( 1 ~> 2 ) -- ^ path
-> ( 1 ~> brushParams ) -- ^ brush parameters
-> ( brushParams ~> SplinePts Closed ) -- ^ brush from parameters
-> ( 1 -> Seq ( 1 -> StrokeDatum ) )
brushStrokeData :: forall i brushParams
. ( Diffy ( I i Double ) ( I i brushParams )
, Fractional ( I i Double )
, D ( I i ( 1 ) ) ~ D ( 1 )
, D ( I i ( 2 ) ) ~ D ( 2 )
, Cross ( I i Double ) ( T ( I i ( 2 ) ) )
, Torsor ( T ( I i ( 2 ) ) ) ( I i ( 2 ) )
, Coercible ( I i ( 1 ) ) ( I i Double )
, Show brushParams
)
=> ( I i ( 1 ) ~> I i ( 2 ) )
-- ^ path
-> ( I i ( 1 ) ~> I i brushParams )
-- ^ brush parameters
-> ( I i brushParams ~> Spline Closed () ( I i ( 2 ) ) )
-- ^ brush from parameters
-> ( I i ( 1 ) -> Seq ( I i ( 1 ) -> StrokeDatum i ) )
brushStrokeData path params brush =
\ t ->
let
dpath_t :: D ( 1 ) ( 2 )
dpath_t :: D ( I i ( 1 ) ) ( I i ( 2 ) )
!dpath_t = runD path t
dparams_t :: D ( 1 ) brushParams
dparams_t :: D ( I i ( 1 ) ) ( I i brushParams )
!dparams_t@( D1 { v = params_t } ) = runD params t
dbrush_params :: D brushParams ( SplinePts Closed )
dbrush_params :: D ( I i brushParams ) ( Spline Closed () ( I i ( 2 ) ) )
!dbrush_params = runD brush params_t
splines :: Seq ( D brushParams ( 1 ~> 2 ) )
!splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns ) dbrush_params
dbrushes_t :: Seq ( 1 -> D ( 2 ) ( 2 ) )
splines :: Seq ( D ( I i brushParams ) ( I i ( 1 ) ~> I i ( 2 ) ) )
!splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @i ) dbrush_params
dbrushes_t :: Seq ( I i ( 1 ) -> D ( I i ( 2 ) ) ( I i ( 2 ) ) )
!dbrushes_t = force $ fmap ( uncurryD . ( dparams_t `chain` ) ) splines
in fmap ( mkStrokeDatum dpath_t ) dbrushes_t
where
mkStrokeDatum :: D ( 1 ) ( 2 )
-> ( 1 -> D ( 2 ) ( 2 ) )
-> ( 1 -> StrokeDatum )
mkStrokeDatum :: D ( I i ( 1 ) ) ( I i ( 2 ) )
-> ( I i ( 1 ) -> D ( I i ( 2 ) ) ( I i ( 2 ) ) )
-> ( I i ( 1 ) -> StrokeDatum i )
mkStrokeDatum dpath_t dbrush_t s =
let dbrush_t_s = dbrush_t s
dstroke@( D2 _c _𝛿c𝛿t _𝛿c𝛿s _ _ _ ) = brushStroke dpath_t dbrush_t_s
( ee, dcdt, 𝛿E𝛿t, 𝛿E𝛿s ) = coerce $ envelopeEquation @Identity $ coerce dstroke
( ee, dcdt, 𝛿E𝛿t, 𝛿E𝛿s ) = envelopeEquation @i dstroke
in -- trace
-- ( unlines
-- [ "envelopeEquation:"
@ -1004,34 +1043,46 @@ brushStrokeData path params brush =
, dstroke
, ee, dcdt, 𝛿E𝛿t, 𝛿E𝛿s }
-- | The value and derivative of a brush stroke at a given coordinate
-- \( (t_0, s_0) \), together with the value of the envelope equation at that
-- point.
data StrokeDatum
data StrokeDatum i
= StrokeDatum
{ -- | Path \( p(t_0) \) (with its 1st and 2nd derivatives).
dpath :: D ( 1 ) ( 2 )
dpath :: D ( I i ( 1 ) ) ( I i ( 2 ) )
-- | Brush shape \( b(t_0, s_0) \) (with its 1st and 2nd derivatives).
, dbrush :: D ( 2 ) ( 2 )
, dbrush :: D ( I i ( 2 ) ) ( I i ( 2 ) )
-- Everything below can be computed in terms of the first two fields.
-- | Stroke \( c(t_0,s_0) = p(t_0) + b(t_0,s_0) \) (with its 1st and 2nd derivatives).
, dstroke :: D ( 2 ) ( 2 )
, dstroke :: D ( I i ( 2 ) ) ( I i ( 2 ) )
-- | Envelope
--
-- \[ E(t_0,s_0) = \left ( \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \right )_{(t_0,s_0)}. \]
, ee :: Double
, ee :: I i Double
-- | \( \left ( \frac{\partial E}{\partial s} \right )_{(t_0,s_0)}. \)
, 𝛿E𝛿s :: Double
, 𝛿E𝛿s :: I i Double
-- | \( \left ( \frac{\partial E}{\partial t} \right )_{(t_0,s_0)}. \)
, 𝛿E𝛿t :: Double
, 𝛿E𝛿t :: I i Double
-- | Total derivative
--
-- \[ \left ( \frac{\mathrm{d} c}{\mathrm{d} t} \right )_{(t_0,s_0)}. \]
--
-- This is ill-defined when \( \frac{\partial E}{\partial s} = 0 \).
, dcdt :: T ( 2 )
, dcdt :: T ( I i ( 2 ) )
}
deriving stock Show
deriving stock instance Show ( StrokeDatum Point )
deriving stock instance Show ( StrokeDatum Interval )
-- Handling points and intervals uniformly.
data Extent = Point | Interval
type I :: Extent -> Type -> Type
type family I i a where
I Point a = a
I Interval a = 𝕀 a

View file

@ -11,8 +11,6 @@ import Control.Applicative
( liftA2 )
import Data.Coerce
( coerce )
import Data.Functor.Identity
( Identity(..) )
import Data.Kind
( Type, Constraint )
import GHC.Generics
@ -45,7 +43,6 @@ type instance D ( 1 ) = D1
type instance D ( 2 ) = D2
type instance D ( 3 ) = D3
type instance D ( Identity a ) = D a
type instance D ( 𝕀 u ) = D u
newtype D0 v = D0 { v :: v }

View file

@ -18,10 +18,6 @@ import Control.Applicative
( liftA2 )
import Control.Monad
( guard )
import Data.Coerce
( coerce )
import Data.Functor.Identity
( Identity(..) )
import Data.Kind
( Type, Constraint )
import Data.Monoid
@ -176,16 +172,6 @@ instance Inner Double ( T ( 2 ) ) where
instance Cross Double ( T ( 2 ) ) where
cross ( V2 x1 y1 ) ( V2 x2 y2 ) = x1 * y2 - x2 * y1
instance Module r ( T m ) => Module ( Identity r ) ( T ( Identity m ) ) where
origin = coerce $ origin @r @( T m )
(^+^) = coerce $ (^+^) @r @( T m )
(^-^) = coerce $ (^-^) @r @( T m )
(*^) = coerce $ (*^) @r @( T m )
instance Inner r ( T m ) => Inner ( Identity r ) ( T ( Identity m ) ) where
(^.^) = coerce $ (^.^) @r @( T m )
instance Cross r ( T m ) => Cross ( Identity r ) ( T ( Identity m ) ) where
cross = coerce $ cross @r @( T m )
-- | Compute whether two vectors point in the same direction,
-- that is, whether each vector is a (strictly) positive multiple of the other.
--