diff --git a/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs b/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs new file mode 100644 index 0000000..9e9c19a --- /dev/null +++ b/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs @@ -0,0 +1,246 @@ +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE RebindableSyntax #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module Math.Bezier.Stroke.EnvelopeEquation + ( StrokeDatum(..) + , CurveOrder(..) + ) where + +-- base +import Prelude hiding ( Num(..), (^) ) +import Data.Coerce + ( Coercible, coerce ) +import Data.Kind + ( Constraint ) +import GHC.TypeNats + ( Nat, type (-) ) + +-- acts +import Data.Act + ( Torsor ((-->)) ) + +-- MetaBrush +import Math.Algebra.Dual +import qualified Math.Bezier.Cubic as Cubic +import qualified Math.Bezier.Quadratic as Quadratic +import Math.Differentiable + ( type ExtentOrder ) +import Math.Interval +import Math.Linear +import Math.Module + ( Module(..), Cross((×)), lerp ) +import Math.Ring + +-------------------------------------------------------------------------------- + +-- | 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 i + = StrokeDatum + { -- | Path \( p(t_0) \). + dpath :: D ( ExtentOrder i ) ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) + -- | Brush shape \( b(t_0, s_0) \). + , dbrush :: D ( ExtentOrder i ) ( 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) \). + , dstroke :: D ( ExtentOrder i ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) + + -- | Envelope function + -- + -- \[ E(t_0,s_0) = \left ( \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \right )_{(t_0,s_0)}. \] + , ee :: D ( ExtentOrder i - 1 ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 1 ) ) + + -- \[ \frac{\partial E}{\partial s} \frac{\mathrm{d} c}{\mathrm{d} t}, \] + -- + -- where \( \frac{\mathrm{d} c}{\mathrm{d} t} \) + -- + -- denotes a total derivative. + , 𝛿E𝛿sdcdt :: D ( ExtentOrder i - 2 ) ( I i ( ℝ 2 ) ) ( T ( I i ( ℝ 2 ) ) ) + + } + +deriving stock instance Show ( StrokeDatum 'Point ) +deriving stock instance Show ( StrokeDatum 'Interval ) + +-------------------------------------------------------------------------------- + +type CurveOrder :: Nat -> Constraint +class CurveOrder k where + + -- | Linear interpolation, as a differentiable function. + line :: forall i b + . ( Module ( I i Double ) ( T b ), Torsor ( T b ) b + , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) + , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) + ) + => Segment b -> C k ( I i ( ℝ 1 ) ) b + + -- | A quadratic Bézier curve, as a differentiable function. + bezier2 :: forall i b + . ( Module ( I i Double ) ( T b ), Torsor ( T b ) b + , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) + , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) + ) + => Quadratic.Bezier b -> C k ( I i ( ℝ 1 ) ) b + + -- | A cubic Bézier curve, as a differentiable function. + bezier3 :: forall i b + . ( Module ( I i Double ) ( T b ), Torsor ( T b ) b + , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) + , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) + ) + => Cubic.Bezier b -> C k ( I i ( ℝ 1 ) ) b + + uncurryD :: D k a ~ D k ( ℝ 1 ) + => D k ( ℝ 1 ) ( C k a b ) -> a -> D k ( ℝ 2 ) b + + -- | A brush stroke, as described by the equation + -- + -- \[ c(t,s) = p(t) + b(t,s) \] + -- + -- where: + -- + -- - \( p(t) \) is the path that the brush follows, and + -- - \( b(t,s) \) is the brush shape, as it varies along the path. + brushStroke :: Module r ( T v ) + => D k ( ℝ 1 ) v -- ^ stroke path \( p(t) \) + -> D k ( ℝ 2 ) v -- ^ brush \( b(t,s) \) + -> D k ( ℝ 2 ) v + + -- | The envelope function + -- + -- \[ E = \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s}, \] + -- + -- together with the 2-vector + -- + -- \[ \frac{\partial E}{\partial s} \frac{\mathrm{d} c}{\mathrm{d} t}, \] + -- + -- where \( \frac{\mathrm{d} c}{\mathrm{d} t} \) + -- + -- denotes a total derivative. + envelopeEquation :: forall i + . ( k ~ ExtentOrder i + , D ( k - 2 ) ( I i ( ℝ 2 ) ) ~ D ( k - 2 ) ( ℝ 2 ) + , D ( k - 1 ) ( I i ( ℝ 2 ) ) ~ D ( k - 1 ) ( ℝ 2 ) + , D k ( I i ( ℝ 2 ) ) ~ D k ( ℝ 2 ) + , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) + , Coercible ( I i Double ) ( I i ( ℝ 1 ) ) + ) + => D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) + -> ( D ( k - 1 ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 1 ) ) + , D ( k - 2 ) ( I i ( ℝ 2 ) ) ( T ( I i ( ℝ 2 ) ) ) ) + + +instance CurveOrder 2 where + line ( Segment a b :: Segment b ) = + D \ ( coerce -> t ) -> + D21 ( lerp @( T b ) t a b ) + ( a --> b ) + origin + + bezier2 ( bez :: Quadratic.Bezier b ) = + D \ ( coerce -> t ) -> + D21 ( Quadratic.bezier @( T b ) bez t ) + ( Quadratic.bezier' bez t ) + ( Quadratic.bezier'' bez ) + + bezier3 ( bez :: Cubic.Bezier b ) = + D \ ( coerce -> t ) -> + D21 ( Cubic.bezier @( T b ) bez t ) + ( Cubic.bezier' bez t ) + ( Cubic.bezier'' bez t ) + + uncurryD = uncurryD2 + + brushStroke ( D21 p p_t p_tt ) ( D22 b b_t b_s b_tt b_ts b_ss ) = + D22 ( unT $ T p ^+^ T b ) + -- c = p + b + + ( p_t ^+^ b_t ) b_s + -- ∂c/∂t = dp/dt + ∂b/∂t + -- ∂c/∂s = ∂b/∂s + + ( p_tt ^+^ b_tt ) b_ts b_ss + -- ∂²c/∂t² = d²p/dt² + ∂²b/∂t² + -- ∂²c/∂t∂s = ∂²b/∂t∂s + -- ∂²c/∂s² = ∂²b/∂s² + + envelopeEquation ( D22 _ c_t c_s c_tt c_ts c_ss ) = + let !ee = c_t × c_s + !ee_t = c_tt × c_s + c_t × c_ts + !ee_s = c_ts × c_s + c_t × c_ss + !𝛿E𝛿sdcdt = ee_s *^ c_t ^-^ ee_t *^ c_s + in ( D12 ( coerce ee ) ( T $ coerce ee_t ) ( T $ coerce ee_s ) + , D0 𝛿E𝛿sdcdt ) + -- Computation of total derivative dc/dt: + -- + -- dc/dt = ∂c/∂t + ∂c/∂s ∂s/∂t + -- ∂s/∂t = - ∂E/∂t / ∂E/∂s + -- + -- ∂E/∂s dc/dt = ∂E/∂s ∂c/∂t - ∂E/∂t ∂c/∂s. + +instance CurveOrder 3 where + + line ( Segment a b :: Segment b ) = + D \ ( coerce -> t ) -> + D31 ( lerp @( T b ) t a b ) + ( a --> b ) + origin + origin + + bezier2 ( bez :: Quadratic.Bezier b ) = + D \ ( coerce -> t ) -> + D31 ( Quadratic.bezier @( T b ) bez t ) + ( Quadratic.bezier' bez t ) + ( Quadratic.bezier'' bez ) + origin + + bezier3 ( bez :: Cubic.Bezier b ) = + D \ ( coerce -> t ) -> + D31 ( Cubic.bezier @( T b ) bez t ) + ( Cubic.bezier' bez t ) + ( Cubic.bezier'' bez t ) + ( Cubic.bezier''' bez ) + + uncurryD = uncurryD3 + + brushStroke + ( D31 p p_t p_tt p_ttt ) + ( D32 b b_t b_s b_tt b_ts b_ss b_ttt b_tts b_tss b_sss ) = + D32 + ( unT $ T p ^+^ T b ) + ( p_t ^+^ b_t ) b_s + ( p_tt ^+^ b_tt ) b_ts b_ss + ( p_ttt ^+^ b_ttt ) b_tts b_tss b_sss + + envelopeEquation + ( D32 _ c_t c_s + c_tt c_ts c_ss + c_ttt c_tts c_tss c_sss ) + = let !ee = c_t × c_s + !ee_t = c_tt × c_s + c_t × c_ts + !ee_s = c_ts × c_s + c_t × c_ss + !ee_tt = c_ttt × c_s + + c_tt × c_ts * 2 + + c_t × c_tts + !ee_ts = c_tts × c_s + + c_tt × c_ss + -- + c_ts × c_ts -- cancels out + + c_t × c_tss + !ee_ss = c_tss × c_s + + c_ts × c_ss * 2 + + c_t × c_sss + !𝛿E𝛿sdcdt = ee_s *^ c_t ^-^ ee_t *^ c_s + !𝛿E𝛿sdcdt_t = ee_ts *^ c_t ^+^ ee_s *^ c_tt + ^-^ ( ee_tt *^ c_s ^+^ ee_t *^ c_ts ) + !𝛿E𝛿sdcdt_s = ee_ss *^ c_t ^+^ ee_s *^ c_ts + ^-^ ( ee_ts *^ c_s ^+^ ee_t *^ c_ss ) + in ( D22 + ( coerce ee ) + ( T $ coerce ee_t ) ( T $ coerce ee_s ) + ( T $ coerce ee_tt) ( T $ coerce ee_ts ) ( T $ coerce ee_ss ) + , D12 𝛿E𝛿sdcdt ( T 𝛿E𝛿sdcdt_t ) ( T 𝛿E𝛿sdcdt_s ) )