add separate EnvelopeEquation module (forgot to add earlier)

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 ) )