do the interval brush stroking at degree 3

This commit is contained in:
sheaf 2023-01-21 15:24:08 +01:00
parent 6945aac704
commit 236055b4ca
12 changed files with 205 additions and 279 deletions

View file

@ -179,6 +179,7 @@ library splines
, Math.Bezier.Quadratic
, Math.Bezier.Spline
, Math.Bezier.Stroke
, Math.Bezier.Stroke.EnvelopeEquation
, Math.Differentiable
, Math.Epsilon
, Math.Interval

View file

@ -29,8 +29,6 @@ import qualified Data.HashMap.Strict as HashMap
-- MetaBrush
import Math.Algebra.Dual
import Math.Bezier.Spline
import Math.Differentiable
( DiffInterp )
import Math.Interval
( type I )
import Math.Linear
@ -94,46 +92,65 @@ circleSpline p = sequenceA $
lastCrv =
Bezier3To ( p κ -1 ) ( p 1 -κ ) BackToStart ()
circleBrush :: forall i k
. ( k ~ 2, DiffInterp i ( Record CircleBrushFields ) )
circleBrush :: forall i k irec
. ( irec ~ I i ( Record CircleBrushFields )
, Module
( D k irec ( I i Double ) )
( D k irec ( I i ( 2 ) ) )
, Module ( I i Double ) ( T ( I i Double ) )
, HasChainRule ( I i Double ) k irec
, Representable ( I i Double ) irec
, Applicative ( D k irec )
)
=> Proxy# i
-> ( forall a. a -> I i a )
-> C k ( I i ( Record CircleBrushFields ) ) ( Spline 'Closed () ( I i ( 2 ) ) )
-> C k irec ( Spline 'Closed () ( I i ( 2 ) ) )
circleBrush _ mkI =
D \ params ->
let r :: D k ( I i ( Record CircleBrushFields ) ) ( I i Double )
let r :: D k irec( I i Double )
r = runD ( var @_ @k ( Fin 1 ) ) params
mkPt :: Double -> Double -> D k ( I i ( Record CircleBrushFields ) ) ( I i ( 2 ) )
mkPt :: Double -> Double -> D k irec ( I i ( 2 ) )
mkPt ( kon -> x ) ( kon -> y )
= ( x * r ) *^ e_x
^+^ ( y * r ) *^ e_y
in circleSpline @i @k @( Record CircleBrushFields ) @( 2 ) mkPt
where
e_x, e_y :: D k ( I i ( Record CircleBrushFields ) ) ( I i ( 2 ) )
e_x, e_y :: D k irec ( I i ( 2 ) )
e_x = pure $ mkI $ 2 1 0
e_y = pure $ mkI $ 2 0 1
kon = konst @( I i Double ) @k @( I i ( Record CircleBrushFields ) ) . mkI
kon = konst @( I i Double ) @k @irec . mkI
ellipseBrush :: forall i k
. ( k ~ 2, DiffInterp i ( Record EllipseBrushFields ) )
ellipseBrush :: forall i k irec
. ( irec ~ I i ( Record EllipseBrushFields )
, Module
( D k irec ( I i Double ) )
( D k irec ( I i ( 2 ) ) )
, Module ( I i Double ) ( T ( I i Double ) )
, HasChainRule ( I i Double ) k irec
, Representable ( I i Double ) irec
, Applicative ( D k irec )
, Transcendental ( D k irec ( I i Double ) )
-- TODO: make a synonym for the above...
-- it seems DiffInterp isn't exactly right
)
=> Proxy# i
-> ( forall a. a -> I i a )
-> C k ( I i ( Record EllipseBrushFields ) ) ( Spline 'Closed () ( I i ( 2 ) ) )
-> C k irec ( Spline 'Closed () ( I i ( 2 ) ) )
ellipseBrush _ mkI =
D \ params ->
let a, b, phi :: D k ( I i ( Record EllipseBrushFields ) ) ( I i Double )
let a, b, phi :: D k irec ( I i Double )
a = runD ( var @_ @k ( Fin 1 ) ) params
b = runD ( var @_ @k ( Fin 2 ) ) params
phi = runD ( var @_ @k ( Fin 3 ) ) params
mkPt :: Double -> Double -> D k ( I i ( Record EllipseBrushFields ) ) ( I i ( 2 ) )
mkPt :: Double -> Double -> D k irec ( I i ( 2 ) )
mkPt ( kon -> x ) ( kon -> y )
= ( x * a * cos phi - y * b * sin phi ) *^ e_x
^+^ ( y * b * cos phi + x * a * sin phi ) *^ e_y
in circleSpline @i @k @( Record EllipseBrushFields ) @( 2 ) mkPt
where
e_x, e_y :: D k ( I i ( Record EllipseBrushFields ) ) ( I i ( 2 ) )
e_x, e_y :: D k irec ( I i ( 2 ) )
e_x = pure $ mkI $ 2 1 0
e_y = pure $ mkI $ 2 0 1
kon = konst @( I i Double ) @k @( I i ( Record EllipseBrushFields ) ) . mkI
kon = konst @( I i Double ) @k @irec . mkI

View file

@ -40,15 +40,14 @@ import qualified Data.Text as Text
-- MetaBrush
import Math.Algebra.Dual
( type (~>) )
import Math.Linear
import Math.Interval
( type I, Extent(Point, Interval) )
( C )
import Math.Bezier.Spline
( SplineType(Closed), Spline )
import Math.Differentiable
( DiffInterp )
( DiffInterp, ExtentOrder )
import Math.Interval
( type I, Extent(Point, Interval) )
import Math.Linear
import MetaBrush.Records
( KnownSymbols, Length, Record )
import MetaBrush.Serialisable
@ -67,7 +66,7 @@ data WithParams params f =
. ( DiffInterp i params )
=> Proxy# i
-> ( forall a. a -> I i a )
-> I i params ~> f ( I i ( 2 ) )
-> C ( ExtentOrder i ) ( I i params ) ( f ( I i ( 2 ) ) )
}
--------------------------------------------------------------------------------

View file

@ -150,6 +150,9 @@ deriving newtype instance HasChainRule Double 2 ( ( Length ks ) )
deriving via 𝕀 ( Length ks )
instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀 ( Length ks ) )
=> HasChainRule ( 𝕀 Double ) 2 ( 𝕀 ( Record ks ) )
deriving via 𝕀 ( Length ks )
instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀 ( Length ks ) )
=> HasChainRule ( 𝕀 Double ) 3 ( 𝕀 ( Record ks ) )
--------------------------------------------------------------------------------

View file

@ -10,7 +10,7 @@
-dsuppress-unfoldings -dsuppress-coercions #-}
module Math.Algebra.Dual
( C(..), D, type (~>), type (~~>)
( C(..), D
, HasChainRule(..), chainRule
, uncurryD2, uncurryD3
, linear, fun, var
@ -49,11 +49,6 @@ type C :: Nat -> Type -> Type -> Type
newtype C k u v = D { runD :: u -> D k u v }
deriving stock instance Functor ( D k u ) => Functor ( C k u )
-- | \( C^2 \)-differentiable mappings.
type (~>) = C 2
-- | \( C^3 \)-differentiable mappings.
type (~~>) = C 3
-- | @D k u v@ is the space of @k@-th order germs of functions from @u@ to @v@,
-- represented by the algebra:
--
@ -64,6 +59,10 @@ type D :: Nat -> Type -> Type -> Type
type family D k u
type instance D k ( 0 ) = D𝔸0
type instance D 0 ( 1 ) = D𝔸0
type instance D 0 ( 2 ) = D𝔸0
type instance D 0 ( 3 ) = D𝔸0
type instance D 0 ( 4 ) = D𝔸0
type instance D 1 ( 1 ) = D1𝔸1
type instance D 1 ( 2 ) = D1𝔸2
@ -93,6 +92,9 @@ instance ( Applicative ( D k u ), Module r ( T v ) )
--------------------------------------------------------------------------------
-- TODO: split up this class into the chain rule operation
-- and all the other operations.
-- | @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

View file

@ -5,7 +5,7 @@
module Math.Bezier.Cubic
( Bezier(..)
, fromQuadratic
, bezier, bezier', bezier''
, bezier, bezier', bezier'', bezier'''
, curvature, squaredCurvature, signedCurvature
, subdivide
, ddist, closestPoint
@ -56,8 +56,8 @@ import Math.Epsilon
import Math.Module
( Module (..)
, lerp
, Inner(..), norm, squaredNorm
, cross
, Inner((^.^)), norm, squaredNorm
, Cross((×))
)
import Math.Roots
( realRoots, solveQuadratic )
@ -119,6 +119,12 @@ bezier'' ( Bezier {..} ) t
( p1 --> p0 ^+^ p1 --> p2 )
( p2 --> p1 ^+^ p2 --> p3 )
-- | Third derivative of a cubic Bézier curve.
bezier''' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> v
bezier''' ( Bezier {..} )
= ( Ring.fromInteger 6 *^ )
$ ( ( p0 --> p3 ) ^+^ Ring.fromInteger 3 *^ ( p2 --> p1 ) )
-- | Curvature of a cubic Bézier curve.
curvature :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> r -> r
curvature bez t = sqrt $ squaredCurvature @v bez t
@ -140,7 +146,7 @@ squaredCurvature bez t
-- | Signed curvature of a planar cubic Bézier curve.
signedCurvature :: Bezier ( 2 ) -> Double -> Double
signedCurvature bez t = ( g' `cross` g'' ) / norm g' ^ ( 3 :: Int )
signedCurvature bez t = ( g' × g'' ) / norm g' ^ ( 3 :: Int )
where
g', g'' :: T ( 2 )
g' = bezier' @( T ( 2 ) ) bez t

View file

@ -54,7 +54,7 @@ import Math.Module
( Module (..)
, lerp
, Inner(..), norm, squaredNorm
, cross
, Cross((×))
)
import Math.Roots
( realRoots )
@ -119,7 +119,7 @@ squaredCurvature bez t
-- | Signed curvature of a planar quadratic Bézier curve.
signedCurvature :: Bezier ( 2 ) -> Double -> Double
signedCurvature bez t = ( g' `cross` g'' ) / norm g' ^ ( 3 :: Int )
signedCurvature bez t = ( g' × g'' ) / norm g' ^ ( 3 :: Int )
where
g', g'' :: T ( 2 )
g' = bezier' @( T ( 2 ) ) bez t

View file

@ -1,5 +1,4 @@
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE PartialTypeSignatures #-}
{-# LANGUAGE QuantifiedConstraints #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE UndecidableInstances #-}
@ -32,7 +31,7 @@ import Control.Monad.ST
import Data.Bifunctor
( Bifunctor(bimap) )
import Data.Coerce
( Coercible, coerce )
( Coercible )
import Data.Foldable
( for_, toList )
import Data.Functor.Identity
@ -49,6 +48,8 @@ import GHC.STRef
( STRef(..), readSTRef, writeSTRef )
import GHC.Generics
( Generic, Generic1 )
import GHC.TypeNats
( type (-) )
-- acts
import Data.Act
@ -113,14 +114,17 @@ import Math.Bezier.Spline
, showSplinePoints
)
import qualified Math.Bezier.Quadratic as Quadratic
import Math.Bezier.Stroke.EnvelopeEquation
import Math.Differentiable
( Differentiable, DiffInterp )
( Differentiable, DiffInterp
, type ExtentOrder
)
import Math.Epsilon
( epsilon )
import Math.Interval
import Math.Linear
import Math.Module
( Module(..), Inner((^.^)), Cross(cross), Interpolatable
( Module(..), Inner((^.^)), Cross((×)), Interpolatable
, lerp, convexCombination, strictlyParallel
)
import Math.Orientation
@ -197,15 +201,15 @@ computeStrokeOutline ::
, NFData ptData, NFData crvData
-- Differentiability.
, DiffInterp 'Point brushParams
, DiffInterp 'Interval brushParams
, Interpolatable Double usedParams
, Interpolatable ( 𝕀 Double ) ( 𝕀 usedParams )
, HasChainRule Double 2 brushParams
, HasChainRule ( 𝕀 Double ) 2 ( 𝕀 brushParams )
, HasChainRule Double 2 usedParams
, HasChainRule ( 𝕀 Double ) 2 ( 𝕀 usedParams )
, Traversable ( D 2 brushParams )
, DiffInterp 'Point brushParams
, DiffInterp 'Interval brushParams
, HasChainRule Double ( ExtentOrder 'Point ) usedParams
, HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 usedParams )
, HasChainRule Double ( ExtentOrder 'Point ) brushParams
, HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 brushParams )
, Traversable ( D ( ExtentOrder 'Point ) brushParams )
-- Debugging.
, Show ptData, Show brushParams
@ -218,7 +222,9 @@ computeStrokeOutline ::
. DiffInterp i brushParams
=> Proxy# i
-> ( forall a. a -> I i a )
-> I i brushParams ~> Spline Closed () ( I i ( 2 ) )
-> C ( ExtentOrder i )
( I i brushParams )
( Spline Closed () ( I i ( 2 ) ) )
)
-> Spline clo crvData ptData
-> ST s
@ -398,7 +404,7 @@ computeStrokeOutline fitParams ptParams toBrushParams brushFn spline@( Spline {
ori = splineOrientation brush0
fwdCond, bwdCond :: Bool
( fwdCond, bwdCond )
| prevTgt `cross` tgt < 0 && prevTgt ^.^ tgt < 0
| prevTgt × tgt < 0 && prevTgt ^.^ tgt < 0
= ( isJust $ between ori prevTgtFwd tgtFwd testTgt1
, isJust $ between ori prevTgtBwd tgtBwd ( -1 *^ testTgt1 )
)
@ -448,14 +454,11 @@ outlineFunction
, Interpolatable ( 𝕀 Double ) ( 𝕀 usedParams )
, DiffInterp 'Point brushParams
, DiffInterp 'Interval brushParams
, HasChainRule Double 2 usedParams
, HasChainRule ( 𝕀 Double ) 2 ( 𝕀 usedParams )
, HasChainRule Double 2 brushParams
, HasChainRule ( 𝕀 Double ) 2 ( 𝕀 brushParams )
, Traversable ( D 2 brushParams )
-- , Diffy Double usedParams
-- , Diffy ( 𝕀 Double ) ( 𝕀 usedParams )
, HasChainRule Double ( ExtentOrder 'Point ) usedParams
, HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 usedParams )
, HasChainRule Double ( ExtentOrder 'Point ) brushParams
, HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 brushParams )
, Traversable ( D ( ExtentOrder 'Point ) brushParams )
-- Debugging.
, Show ptData, Show brushParams
@ -466,15 +469,19 @@ outlineFunction
. DiffInterp i brushParams
=> Proxy# i
-> ( forall a. a -> I i a )
-> I i brushParams ~> Spline Closed () ( I i ( 2 ) )
-> C ( ExtentOrder i )
( I i brushParams )
( Spline Closed () ( I i ( 2 ) ) )
)
-> ptData
-> Curve Open crvData ptData
-> OutlineFn
outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
let
pathAndUsedParams :: forall i
. ( D 2 ( I i ( 1 ) ) ~ D 2 ( 1 )
pathAndUsedParams :: forall i k arr
. ( k ~ ExtentOrder i, CurveOrder k
, arr ~ C k
, D k ( I i ( 1 ) ) ~ D k ( 1 )
, Coercible ( I i ( 1 ) ) ( I i Double )
, Module ( I i Double ) ( T ( I i ( 2 ) ) )
, Torsor ( T ( I i ( 2 ) ) ) ( I i ( 2 ) )
@ -482,34 +489,37 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
, Torsor ( T ( I i usedParams ) ) ( I i usedParams )
)
=> ( forall a. a -> I i a )
-> ( I i ( 1 ) ~> I i ( 2 ), I i ( 1 ) ~> I i usedParams )
-> ( I i ( 1 ) `arr` I i ( 2 ), I i ( 1 ) `arr` I i usedParams )
pathAndUsedParams toI =
case crv of
LineTo { curveEnd = NextPoint sp1 }
| let seg = Segment sp0 sp1
-> ( line @i ( fmap ( toI . coords ) seg )
, line @i ( fmap ( toI . ptParams ) seg ) )
-> ( line @k @i ( fmap ( toI . coords ) seg )
, line @k @i ( fmap ( toI . ptParams ) seg ) )
Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 }
| let bez2 = Quadratic.Bezier sp0 sp1 sp2
-> ( bezier2 @i ( fmap ( toI . coords ) bez2 )
, bezier2 @i ( fmap ( toI . ptParams ) bez2 ) )
-> ( bezier2 @k @i ( fmap ( toI . coords ) bez2 )
, bezier2 @k @i ( fmap ( toI . ptParams ) bez2 ) )
Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 }
| let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3
-> ( bezier3 @i ( fmap ( toI . coords ) bez3 )
, bezier3 @i ( fmap ( toI . ptParams ) bez3 ) )
-> ( bezier3 @k @i ( fmap ( toI . coords ) bez3 )
, bezier3 @k @i ( fmap ( toI . ptParams ) bez3 ) )
usedParams :: 1 ~> usedParams
path :: 1 ~> 2
usedParams :: C ( ExtentOrder 'Point ) ( 1 ) usedParams
path :: C ( ExtentOrder 'Point ) ( 1 ) ( 2 )
( path, usedParams ) = pathAndUsedParams @Point id
curvesI :: 𝕀 1 -> Seq ( 𝕀 1 -> StrokeDatum 'Interval )
curvesI = brushStrokeData @'Interval @brushParams
curvesI = brushStrokeData @'Interval @( ExtentOrder 'Interval ) @brushParams
pathI
( chainRule @( 𝕀 Double ) @2 usedParamsI $ linear ( nonDecreasing toBrushParams ) )
( chainRule @( 𝕀 Double ) @( ExtentOrder 'Interval )
usedParamsI
( linear ( nonDecreasing toBrushParams ) )
)
( brushFromParams @'Interval proxy# singleton )
usedParamsI :: 𝕀 1 ~> 𝕀 usedParams
pathI :: 𝕀 1 ~> 𝕀 2
usedParamsI :: C ( ExtentOrder 'Interval ) ( 𝕀 1 ) ( 𝕀 usedParams )
pathI :: C ( ExtentOrder 'Interval ) ( 𝕀 1 ) ( 𝕀 2 )
( pathI, usedParamsI ) = pathAndUsedParams @'Interval singleton
fwdBwd :: OutlineFn
@ -520,9 +530,12 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
where
curves :: Seq ( 1 -> StrokeDatum Point )
curves = brushStrokeData @Point @brushParams
curves = brushStrokeData @Point @( ExtentOrder 'Point ) @brushParams
path
( chainRule @Double usedParams $ linear toBrushParams )
( chainRule @Double @( ExtentOrder 'Point )
usedParams
( linear toBrushParams )
)
( brushFromParams @Point proxy# id )
t
@ -537,11 +550,11 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
bisSols = bisection 0.0001 curvesI
in trace
( unlines $
( "bisectionMethod: #(possible zeroes) = " ++ show ( length bisSols ) ) :
"" :
map show bisSols )
in --trace
-- ( unlines $
-- ( "bisectionMethod: #(possible zeroes) = " ++ show ( length bisSols ) ) :
-- "" :
-- map show bisSols )
fwdBwd
-----------------------------------
@ -799,9 +812,9 @@ withTangent tgt_wanted spline@( Spline { splineStart } )
bez :: Cubic.Bezier ( 2 )
bez = Cubic.Bezier {..}
c01, c12, c23 :: Double
c01 = tgt_wanted `cross` tgt0
c12 = tgt_wanted `cross` tgt1
c23 = tgt_wanted `cross` tgt2
c01 = tgt_wanted × tgt0
c12 = tgt_wanted × tgt1
c23 = tgt_wanted × tgt2
correctTangentParam :: Double -> Maybe Double
correctTangentParam t
| t > -epsilon && t < 1 + epsilon
@ -823,109 +836,15 @@ withTangent tgt_wanted spline@( Spline { splineStart } )
, offset = T $ Cubic.bezier @( T ( 2 ) ) bez t
}
--------------------------------------------------------------------------------
-- | 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 2 ( 1 ) v -- ^ stroke path \( p(t) \)
-> D 2 ( 2 ) v -- ^ brush \( b(t,s) \)
-> D 2 ( 2 ) v
brushStroke ( D21 p dpdt d2pdt2 ) ( D22 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) =
D22 ( unT $ T p ^+^ T b )
-- c = p + b
( dpdt ^+^ dbdt ) dbds
-- ∂c/∂t = dp/dt + ∂b/∂t
-- ∂c/∂s = ∂b/∂s
( d2pdt2 ^+^ d2bdt2 ) d2bdtds d2bds2
-- ∂²c/∂t² = d²p/dt² + ∂²b/∂t²
-- ∂²c/∂t∂s = ∂²b/∂t∂s
-- ∂²c/∂s² = ∂²b/∂s²
-- | The envelope equation
--
-- \[ E = \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} = 0, \]
--
-- together with the total derivative
--
-- \[ \frac{\mathrm{d} c}{\mathrm{d} t}, \]
--
-- and the partial derivatives
--
-- \[ \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 :: forall i
. ( D 2 ( I i ( 2 ) ) ~ D 2 ( 2 )
, Cross ( I i Double ) ( T ( I i ( 2 ) ) )
, Fractional ( I i Double )
)
=> D 2 ( I i ( 2 ) ) ( I i ( 2 ) )
-> ( I i Double, T ( I i ( 2 ) ), T ( I i ( 2 ) ), I i Double, I i Double )
envelopeEquation ( D22 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ) =
let ee = dcdt `cross` dcds
dEdt = d2cdt2 `cross` dcds + dcdt `cross` d2cdtds
dEds = d2cdtds `cross` dcds + dcdt `cross` d2cds2
tot = dcdt -- ^-^ ( dEdt / dEds ) *^ dcds
dEdsTot = dEds *^ dcdt ^-^ dEdt *^ dcds
in ( ee, tot, dEdsTot, dEdt, dEds )
-- 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.
-- | Linear interpolation, as a differentiable function.
line :: forall i b
. ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
, D 2 ( I i ( 1 ) ) ~ D 2 ( 1 )
, Coercible ( I i ( 1 ) ) ( I i Double )
)
=> Segment b -> I i ( 1 ) ~> b
line ( Segment a b ) = D \ ( coerce -> t ) ->
D21 ( lerp @( T b ) t a b )
( a --> b )
origin
-- | A quadratic Bézier curve, as a differentiable function.
bezier2 :: forall i b
. ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
, D 2 ( I i ( 1 ) ) ~ D 2 ( 1 )
, Coercible ( I i ( 1 ) ) ( I i Double )
)
=> Quadratic.Bezier b -> I i ( 1 ) ~> b
bezier2 bez = D \ ( coerce -> t ) ->
D21 ( Quadratic.bezier @( T b ) bez t )
( Quadratic.bezier' bez t )
( Quadratic.bezier'' bez )
-- | A cubic Bézier curve, as a differentiable function.
bezier3 :: forall i b
. ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
, D 2 ( I i ( 1 ) ) ~ D 2 ( 1 )
, Coercible ( I i ( 1 ) ) ( I i Double )
)
=> Cubic.Bezier b -> I i ( 1 ) ~> b
bezier3 bez = D \ ( coerce -> t ) ->
D21 ( Cubic.bezier @( T b ) bez t )
( Cubic.bezier' bez t )
( Cubic.bezier'' bez t )
splineCurveFns :: forall i
. ( D 2 ( I i ( 1 ) ) ~ D 2 ( 1 )
splineCurveFns :: forall i k
. ( k ~ ExtentOrder i, CurveOrder k
, D k ( I i ( 1 ) ) ~ D k ( 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 ) )
=> Spline Closed () ( I i ( 2 ) ) -> Seq ( C k ( I i ( 1 ) ) ( I i ( 2 ) ) )
splineCurveFns spls
= runIdentity
. bifoldSpline
@ -936,14 +855,14 @@ splineCurveFns spls
where
curveFn :: I i ( 2 )
-> Curve Open () ( I i ( 2 ) )
-> ( I i ( 1 ) ~> I i ( 2 ) )
-> ( C k ( I i ( 1 ) ) ( I i ( 2 ) ) )
curveFn p0 = \case
LineTo { curveEnd = NextPoint p1 }
-> line @i $ Segment p0 p1
-> line @k @i $ Segment p0 p1
Bezier2To { controlPoint = p1, curveEnd = NextPoint p2 }
-> bezier2 @i $ Quadratic.Bezier p0 p1 p2
-> bezier2 @k @i $ Quadratic.Bezier p0 p1 p2
Bezier3To { controlPoint1 = p1, controlPoint2 = p2, curveEnd = NextPoint p3 }
-> bezier3 @i $ Cubic.Bezier p0 p1 p2 p3
-> bezier3 @k @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.
@ -1019,15 +938,21 @@ solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData
Nothing -> ( False, initialGuess )
Just s0 -> ( True , s0 )
in case f ( 1 s ) of -- TODO: a bit redundant to have to compute this again...
StrokeDatum { ee = _ee, dstroke, 𝛿E𝛿t = _𝛿E𝛿t, 𝛿E𝛿s, dcdt } ->
StrokeDatum
{ dstroke
, ee = D12 ( 1 _ee ) ( T ( 1 _𝛿E𝛿t ) ) ( T ( 1 ee_s ) )
, 𝛿E𝛿sdcdt = D0 𝛿E𝛿sdcdt
} ->
-- The total derivative dc/dt is computed by dividing by ∂E/∂s,
-- so check it isn't zero first. This corresponds to cusps in the envelope.
let totDeriv
| abs 𝛿E𝛿s < epsilon
let dcdt
| abs ee_s < epsilon
, let s' = if s >= 0.5 then s - 1e-9 else s + 1e-9
= case f ( 1 s' ) of { StrokeDatum { dcdt = dcdt_s' } -> dcdt_s' }
= case f ( 1 s' ) of
StrokeDatum { ee = D12 _ _ ( T ( 1 ee_s' ) ), 𝛿E𝛿sdcdt = D0 𝛿E𝛿sdcdt' }
-> recip ee_s' *^ 𝛿E𝛿sdcdt'
| otherwise
= dcdt
= recip ee_s *^ 𝛿E𝛿sdcdt
in --trace
-- ( unlines
-- [ "solveEnvelopeEquations"
@ -1036,17 +961,16 @@ solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData
-- , " c = " ++ show dstroke
-- , " E = " ++ show _ee
-- , " ∂E/∂t = " ++ show _𝛿E𝛿t
-- , " ∂E/∂s = " ++ show 𝛿E𝛿s
-- , " dc/dt = " ++ show totDeriv
-- , " ∂E/∂s = " ++ show ee_s
-- , " dc/dt = " ++ show dcdt
-- ] )
( good, 1 s, value @Double @2 @( 2 ) dstroke
, totDeriv )
( good, 1 s, value @Double @2 @( 2 ) dstroke, dcdt )
eqn :: ( 1 -> StrokeDatum Point ) -> ( Double -> ( Double, Double ) )
eqn f s =
case f ( 1 s ) of
StrokeDatum { ee, 𝛿E𝛿s } ->
( ee, 𝛿E𝛿s )
StrokeDatum { ee = D12 ( 1 ee ) _ ( T ( 1 ee_s ) ) } ->
( ee, ee_s )
maxIters :: Word
maxIters = 5 --30
@ -1061,48 +985,55 @@ 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 i brushParams
. ( Differentiable i brushParams
brushStrokeData :: forall i k brushParams arr
. ( k ~ ExtentOrder i, CurveOrder k, arr ~ C k
, Differentiable i brushParams
, Fractional ( I i Double )
, D 2 ( I i ( 1 ) ) ~ D 2 ( 1 )
, D 2 ( I i ( 2 ) ) ~ D 2 ( 2 )
, HasChainRule ( I i Double ) k ( I i ( 1 ) )
, Applicative ( D k ( 1 ) )
, D ( k - 2 ) ( I i ( 2 ) ) ~ D ( k - 2 ) ( 2 )
, D ( k - 1 ) ( I i ( 2 ) ) ~ D ( k - 1 ) ( 2 )
, D k ( I i ( 1 ) ) ~ D k ( 1 )
, D k ( I i ( 1 ) ) ~ D k ( 1 )
, D k ( I i ( 2 ) ) ~ D k ( 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 ) )
=> ( I i ( 1 ) `arr` I i ( 2 ) )
-- ^ path
-> ( I i ( 1 ) ~> I i brushParams )
-> ( I i ( 1 ) `arr` I i brushParams )
-- ^ brush parameters
-> ( I i brushParams ~> Spline Closed () ( I i ( 2 ) ) )
-> ( I i brushParams `arr` 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 2 ( I i ( 1 ) ) ( I i ( 2 ) )
dpath_t :: D k ( I i ( 1 ) ) ( I i ( 2 ) )
!dpath_t = runD path t
dparams_t :: D 2 ( I i ( 1 ) ) ( I i brushParams )
!dparams_t@( D21 { _D21_v = params_t } ) = runD params t
dbrush_params :: D 2 ( I i brushParams ) ( Spline Closed () ( I i ( 2 ) ) )
!dbrush_params = runD brush params_t
splines :: Seq ( D 2 ( I i brushParams ) ( I i ( 1 ) ~> I i ( 2 ) ) )
!splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @i ) dbrush_params
dbrushes_t :: Seq ( I i ( 1 ) -> D 2 ( I i ( 2 ) ) ( I i ( 2 ) ) )
!dbrushes_t = force $ fmap ( uncurryD2 . ( chain @(I i Double) @2 dparams_t ) ) splines
dparams_t :: D k ( I i ( 1 ) ) ( I i brushParams )
!dparams_t = runD params t
dbrush_params :: D k ( I i brushParams ) ( Spline Closed () ( I i ( 2 ) ) )
!dbrush_params = runD brush $ value @( I i Double ) @k @( I i ( 1 ) ) dparams_t
splines :: Seq ( D k ( I i brushParams ) ( I i ( 1 ) `arr` I i ( 2 ) ) )
!splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @i @k ) dbrush_params
dbrushes_t :: Seq ( I i ( 1 ) -> D k ( I i ( 2 ) ) ( I i ( 2 ) ) )
!dbrushes_t = force $ fmap ( uncurryD @k . ( chain @(I i Double) @k dparams_t ) ) splines
-- This is the crucial use of the chain rule.
in fmap ( mkStrokeDatum dpath_t ) dbrushes_t
where
mkStrokeDatum :: D 2 ( I i ( 1 ) ) ( I i ( 2 ) )
-> ( I i ( 1 ) -> D 2 ( I i ( 2 ) ) ( I i ( 2 ) ) )
mkStrokeDatum :: D k ( I i ( 1 ) ) ( I i ( 2 ) )
-> ( I i ( 1 ) -> D k ( 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@( D22 _c _𝛿c𝛿t _𝛿c𝛿s _ _ _ ) = brushStroke dpath_t dbrush_t_s
( ee, dcdt, 𝛿E𝛿sdcdt, 𝛿E𝛿t, 𝛿E𝛿s ) = envelopeEquation @i dstroke
dstroke = brushStroke @k dpath_t dbrush_t_s
( ee, 𝛿E𝛿sdcdt ) = envelopeEquation @k @i dstroke
in -- trace
-- ( unlines
-- [ "envelopeEquation:"
@ -1113,62 +1044,29 @@ brushStrokeData path params brush =
-- , " ∂c/∂s = " ++ show _𝛿c𝛿s
-- , " E = " ++ show ee
-- , " ∂E/∂t = " ++ show _𝛿E𝛿t
-- , " ∂E/∂s = " ++ show 𝛿E𝛿s
-- , " ∂E/∂s = " ++ show ee_s
-- , " dc/dt = " ++ show dcdt ] ) $
StrokeDatum
{ dpath = dpath_t
, dbrush = dbrush_t_s
, dstroke
, ee, dcdt, 𝛿E𝛿sdcdt, 𝛿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 i
= StrokeDatum
{ -- | Path \( p(t_0) \) (with its 1st and 2nd derivatives).
dpath :: D 2 ( I i ( 1 ) ) ( I i ( 2 ) )
-- | Brush shape \( b(t_0, s_0) \) (with its 1st and 2nd derivatives).
, dbrush :: D 2 ( 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 ( 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 :: I i Double
-- | \( \left ( \frac{\partial E}{\partial s} \right )_{(t_0,s_0)}. \)
, 𝛿E𝛿s :: I i Double
-- | \( \left ( \frac{\partial E}{\partial t} \right )_{(t_0,s_0)}. \)
, 𝛿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, 𝛿E𝛿sdcdt :: T ( I i ( 2 ) )
}
deriving stock instance Show ( StrokeDatum 'Point )
deriving stock instance Show ( StrokeDatum 'Interval )
, ee
, 𝛿E𝛿sdcdt
}
--------------------------------------------------------------------------------
bisection :: Double
-> ( 𝕀 1 -> Seq ( 𝕀 1 -> StrokeDatum 'Interval ) )
-> [ ( 𝕀 1, Int, 𝕀 1, 𝕀 Double, 𝕀 2 ) ]
bisection minWidth eqs = bisect initialCands [] []
-> [ ( 𝕀 1, Int, 𝕀 1, 𝕀 1, 𝕀 2 ) ]
bisection minWidth eqs =
bisect initialCands [] []
where
bisect :: [ ( 𝕀 1, Int, 𝕀 1, 𝕀 Double, 𝕀 2 ) ] -- have solutions, need bisection to refine
bisect :: [ ( 𝕀 1, Int, 𝕀 1, 𝕀 1, 𝕀 2 ) ] -- have solutions, need bisection to refine
-> [ ( 𝕀 1, Int, 𝕀 1 ) ] -- have been bisected, don't know if they contain solutions yet
-> [ ( 𝕀 1, Int, 𝕀 1, 𝕀 Double, 𝕀 2 ) ] -- have solutions, don't bisect further
-> [ ( 𝕀 1, Int, 𝕀 1, 𝕀 Double, 𝕀 2 ) ]
-> [ ( 𝕀 1, Int, 𝕀 1, 𝕀 1, 𝕀 2 ) ] -- have solutions, don't bisect further
-> [ ( 𝕀 1, Int, 𝕀 1, 𝕀 1, 𝕀 2 ) ]
bisect [] [] sols = sols
bisect cands ( ( t, i, s ) : toTry ) sols
| Just ( ee, 𝛿E𝛿sdcdt ) <- isCand t i s
@ -1211,17 +1109,19 @@ bisection minWidth eqs = bisect initialCands [] []
[ (t, i, s, ee, 𝛿E𝛿sdcdt )
| let !eqs_t = eqs t
, ( eq_t, i ) <- zip ( toList eqs_t ) ( [0,1..] :: [Int] )
, let !( StrokeDatum { ee, 𝛿E𝛿sdcdt = T 𝛿E𝛿sdcdt } ) = eq_t s
, Interval.inf ( ival ee ) < 0 && Interval.sup ( ival ee ) > 0
, cmp2 (<) ( getRounded ( Interval.inf $ ival 𝛿E𝛿sdcdt ) ) ( 2 0 0 )
&& cmp2 (>) ( getRounded ( Interval.sup $ ival 𝛿E𝛿sdcdt ) ) ( 2 0 0 )
, let !( StrokeDatum { ee = D22 ee _ _ _ _ _, 𝛿E𝛿sdcdt = D12 ( T 𝛿E𝛿sdcdt ) _ _ } ) = eq_t s
, Interval.inf ( ival ee ) < Rounded ( 1 0 )
, Interval.sup ( ival ee ) > Rounded ( 1 0 )
, cmp2 (<) ( getRounded ( Interval.inf $ ival 𝛿E𝛿sdcdt ) ) ( 2 0 0 )
, cmp2 (>) ( getRounded ( Interval.sup $ ival 𝛿E𝛿sdcdt ) ) ( 2 0 0 )
]
isCand :: 𝕀 1 -> Int -> 𝕀 1 -> Maybe ( 𝕀 Double, 𝕀 2 )
isCand :: 𝕀 1 -> Int -> 𝕀 1 -> Maybe ( 𝕀 1, 𝕀 2 )
isCand t i s = case ( ( eqs t ) `Seq.index` i ) s of
StrokeDatum { ee, 𝛿E𝛿sdcdt = T 𝛿E𝛿sdcdt } ->
StrokeDatum { ee = D22 ee _ _ _ _ _, 𝛿E𝛿sdcdt = D12 ( T 𝛿E𝛿sdcdt ) _ _ } ->
do guard $
Interval.inf ( ival ee ) < 0 && Interval.sup ( ival ee ) > 0
Interval.inf ( ival ee ) < Rounded ( 1 0 )
&& Interval.sup ( ival ee ) > Rounded ( 1 0 )
&& cmp2 (<) ( getRounded ( Interval.inf $ ival 𝛿E𝛿sdcdt ) ) ( 2 0 0 )
&& cmp2 (>) ( getRounded ( Interval.sup $ ival 𝛿E𝛿sdcdt ) ) ( 2 0 0 )
return ( ee, 𝛿E𝛿sdcdt )

View file

@ -1,7 +1,7 @@
{-# LANGUAGE UndecidableInstances #-}
module Math.Differentiable
( Differentiable, DiffInterp )
( ExtentOrder, Differentiable, DiffInterp )
where
-- base
@ -24,12 +24,10 @@ import Math.Ring
type ExtentOrder :: Extent -> Nat
type family ExtentOrder e where
ExtentOrder i = 2
--ExtentOrder 'Point = 2
--ExtentOrder 'Interval = 2
ExtentOrder 'Point = 2
ExtentOrder 'Interval = 3
-- Currently we're doing order 2 derivatives for the brush stroke fitting,
-- but order 3 derivatives for the interval Newton method to find cusps.
-- TODO: using 2 for both until migration finishes.
type Differentiable :: Extent -> Type -> Constraint
class

View file

@ -120,7 +120,7 @@ instance Inner ( 𝕀 Double ) ( T ( 𝕀 2 ) ) where
in x1x2 + y1y2
instance Cross ( 𝕀 Double ) ( T ( 𝕀 2 ) ) where
T ( 𝕀 ( 2 x1_lo y1_lo ) ( 2 x1_hi y1_hi ) ) `cross`
T ( 𝕀 ( 2 x1_lo y1_lo ) ( 2 x1_hi y1_hi ) ) ×
T ( 𝕀 ( 2 x2_lo y2_lo ) ( 2 x2_hi y2_hi ) )
= let !x1y2 = 𝕀 x1_lo x1_hi * 𝕀 y2_lo y2_hi
!y2x1 = 𝕀 x2_lo x2_hi * 𝕀 y1_lo y1_hi

View file

@ -63,7 +63,7 @@ lerp :: forall v r p. ( Module r v, Torsor v p ) => r -> p -> p -> p
lerp t p0 p1 = ( t *^ ( p0 --> p1 :: v ) ) p0
class Module r m => Cross r m where
cross :: m -> m -> r
(×) :: m -> m -> r
-- | Norm of a vector, computed using the inner product.
norm :: forall m r. ( Floating r, Inner r m ) => m -> r
@ -122,7 +122,7 @@ instance Inner Double ( T ( 2 ) ) where
V2 x1 y1 ^.^ V2 x2 y2 = x1 Ring.* x2 + y1 Ring.* y2
instance Cross Double ( T ( 2 ) ) where
cross ( V2 x1 y1 ) ( V2 x2 y2 ) = x1 Ring.* y2 Ring.- x2 Ring.* y1
V2 x1 y1 × V2 x2 y2 = x1 Ring.* y2 Ring.- x2 Ring.* y1
-- | Compute whether two vectors point in the same direction,
-- that is, whether each vector is a (strictly) positive multiple of the other.
@ -130,8 +130,8 @@ instance Cross Double ( T ( 2 ) ) where
-- Returns @False@ if either of the vectors is zero.
strictlyParallel :: T ( 2 ) -> T ( 2 ) -> Bool
strictlyParallel u v
= abs ( u `cross` v ) < epsilon -- vectors are collinear
&& u ^.^ v > epsilon -- vectors point in the same direction (parallel and not anti-parallel)
= abs ( u × v ) < epsilon -- vectors are collinear
&& u ^.^ v > epsilon -- vectors point in the same direction (parallel and not anti-parallel)
-- | Finds whether the query vector @ u @ is a convex combination of the two provided vectors @ v0 @, @ v1 @.
--
@ -159,8 +159,8 @@ convexCombination v0 v1 u
where
c0, c10 :: Double
c0 = v0 `cross` u
c10 = ( v0 ^-^ v1 ) `cross` u
c0 = v0 × u
c10 = ( v0 ^-^ v1 ) × u
--------------------------------------------------------------------------------
-- Not sure how to set things up to automate the following...

View file

@ -32,7 +32,7 @@ import Data.Generics.Internal.VL
import Math.Epsilon
( nearZero )
import Math.Module
( cross )
( (×) )
import Math.Bezier.Spline
( Spline(..), Curves(..), Curve(..), NextPoint(..)
, SplineType(..), KnownSplineType(..), SSplineType(..)
@ -63,7 +63,7 @@ convexOrientation ( v1 : v2 : vs )
= CW
where
crossProduct :: Double
crossProduct = v1 `cross` v2
crossProduct = v1 × v2
convexOrientation _ = CCW -- default
-- | Compute the orientation of a spline, assuming tangent vectors have a monotone angle.