From 09c1bdd9488ead6e8ca76fa808debf8ec44905be Mon Sep 17 00:00:00 2001 From: sheaf Date: Fri, 13 Jan 2023 15:23:33 +0100 Subject: [PATCH] continue interval arithmetic integration --- MetaBrush.cabal | 1 + src/splines/Math/Bezier/Spline.hs | 11 +- src/splines/Math/Bezier/Stroke.hs | 163 ++++++++++++++++++++---------- src/splines/Math/Linear/Dual.hs | 3 - src/splines/Math/Module.hs | 14 --- 5 files changed, 115 insertions(+), 77 deletions(-) diff --git a/MetaBrush.cabal b/MetaBrush.cabal index e5f9374..8a46e24 100644 --- a/MetaBrush.cabal +++ b/MetaBrush.cabal @@ -76,6 +76,7 @@ common common PatternSynonyms RankNTypes RecordWildCards + RoleAnnotations StandaloneDeriving StandaloneKindSignatures TupleSections diff --git a/src/splines/Math/Bezier/Spline.hs b/src/splines/Math/Bezier/Spline.hs index 90b7706..369dbcc 100644 --- a/src/splines/Math/Bezier/Spline.hs +++ b/src/splines/Math/Bezier/Spline.hs @@ -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 ) } diff --git a/src/splines/Math/Bezier/Stroke.hs b/src/splines/Math/Bezier/Stroke.hs index 411264b..c4f75c2 100644 --- a/src/splines/Math/Bezier/Stroke.hs +++ b/src/splines/Math/Bezier/Stroke.hs @@ -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 diff --git a/src/splines/Math/Linear/Dual.hs b/src/splines/Math/Linear/Dual.hs index ee66200..fdf8ab7 100644 --- a/src/splines/Math/Linear/Dual.hs +++ b/src/splines/Math/Linear/Dual.hs @@ -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 ) = Dℝ1 type instance D ( ℝ 2 ) = Dℝ2 type instance D ( ℝ 3 ) = Dℝ3 -type instance D ( Identity a ) = D a type instance D ( 𝕀 u ) = D u newtype Dℝ0 v = D0 { v :: v } diff --git a/src/splines/Math/Module.hs b/src/splines/Math/Module.hs index 92d44f0..61a1937 100644 --- a/src/splines/Math/Module.hs +++ b/src/splines/Math/Module.hs @@ -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. --