brush-strokes-0.1.0.0: src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE DuplicateRecordFields #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE RebindableSyntax #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE UndecidableInstances #-}
module Math.Bezier.Stroke.EnvelopeEquation
( StrokeDatum(..), CornerStrokeDatum(..)
, HasBézier(..)
, HasEnvelopeEquation(..)
) where
-- base
import Prelude hiding ( Num(..), (^), pi, sin, cos )
import Data.Coerce
( coerce )
import Data.Kind
( Type, Constraint )
import Data.List.NonEmpty
( NonEmpty(..) )
import GHC.TypeNats
( Nat, type (-) )
-- acts
import Data.Act
( Torsor ((-->)) )
-- brush-strokes
import Calligraphy.Brushes
( Corner (..), Unrotated (..) )
import Math.Algebra.Dual
import qualified Math.Bezier.Cubic as Cubic
import qualified Math.Bezier.Quadratic as Quadratic
import Math.Differentiable
( I, IVness(..), SingIVness (..) )
import Math.Interval
import Math.Linear
import Math.Module
( Module(..), Cross((×)), lerp )
import Math.Ring
import qualified Math.Ring as Ring
--------------------------------------------------------------------------------
-- | Quantities related to the envelope equation for a brush stroke.
--
-- These are used to formulate:
--
-- - the envelope equation, that determines whether a point is on the outline,
-- - the cusp equation, that determines the presence of cusps in the outline.
type StrokeDatum :: Nat -> IVness -> Type
data StrokeDatum k i
= StrokeDatum
{ -- | Path \( p(t) \).
dpath :: D k 1 ( I i 2 )
-- | Brush shape \( b(t, s) \), not including rotation.
, dbrush :: Unrotated ( D k 2 ( I i 2 ) )
-- | (Optional) rotation angle \( \theta(t) \).
, mbRotation :: Maybe ( D k 1 ( I i Double ) )
-- Everything below is computed in terms of the first three fields.
-- | Stroke shape, including rotation (if any):
--
-- \( c(t,s) = p(t) + R(\theta(t)) b(t,s) \).
, stroke :: I i 2
-- TODO: return this as D 2 1 ( I i 2 ), i.e. include
-- the total derivative with respect to t.
-- | \( u(t,s) = R(-\theta(t)) \frac{\partial c}{\partial t} \),
-- \( v(t,s) = R(-\theta(t)) \frac{\partial c}{\partial s} \)
, du, dv :: D ( k - 1 ) 2 ( I i 2 )
-- | Envelope function
--
-- \[ \begin{align} E(t,s) &= \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \\
-- &= \left ( R(-\theta(t)) \frac{\partial c}{\partial t} \right ) \times \left ( R(-\theta(t)) \frac{\partial c}{\partial s} \right ) \\
-- &= \left ( R(-\theta(t)) p'(t) + \theta'(t) S b(t,s) + \frac{\partial b}{\partial t} \right ) \times \frac{\partial b}{\partial s} \]
--
-- where \( S = R(\pi/2) \) is a constant rotation.
, ee :: D ( k - 1 ) 2 ( I i 1 )
-- | \[ \frac{\partial E}{\partial s} R(-\theta(t)) \frac{\mathrm{d} c}{\mathrm{d} t}, \]
--
-- where \( \frac{\mathrm{d} c}{\mathrm{d} t} \)
--
-- denotes a total derivative, which is computed by implicit differentiation
-- of the envelope equation \( E(t,s) = 0 \):
--
-- \[ \frac{\mathrm{d} c}{\mathrm{d} t} = \frac{\partial c}{\partial t} - \frac{\partial c}{\partial s} \frac{\partial E / \partial t }{ \partial E / \partial s } \]
--
-- \[ R(-\theta(t)) \frac{\mathrm{d} c}{\mathrm{d} t} = u - v \frac{\partial E / \partial t }{ \partial E / \partial s } \]
, 𝛿E𝛿s_unrot_dcdt :: Unrotated ( D ( k - 2 ) 2 ( T ( I i 2 ) ) )
}
deriving stock instance
( Show ( I i 2 )
, Show ( D k 1 ( I i Double ) )
, Show ( D k 1 ( I i 2 ) )
, Show ( D k 2 ( I i 2 ) )
, Show ( D ( k - 1 ) 2 ( I i 1 ) )
, Show ( D ( k - 1 ) 2 ( I i 2 ) )
, Show ( D ( k - 2 ) 2 ( T ( I i 2 ) ) )
)
=> Show ( StrokeDatum k i )
-- | Like 'StrokeDatum', except that it applies for the parts of the stroke
-- which are drawn out by a corner of a brush.
type CornerStrokeDatum :: Nat -> IVness -> Type
data CornerStrokeDatum k i
= CornerStrokeDatum
{ -- | Path \( p(t) \).
dpath :: D k 1 ( I i 2 )
-- | The coordinate of the brush corner \( b(t) \) (not including rotation),
-- together with the start/end angles of the corner (not including rotation).
, dcorner :: Unrotated ( D k 1 ( Corner ( I i 2 ) ) )
-- | (Optional) rotation angle \( \theta(t) \).
, mbRotation :: Maybe ( D k 1 ( I i Double ) )
-- Everything below is computed in terms of the first three fields.
-- | The path being traced out by the corner of the brush,
-- including the rotation (if any):
--
-- \( c(t) = p(t) + R(\theta(t)) b(t) \)
, dstroke :: D k 1 ( Corner ( I i 2 ) )
}
deriving stock instance
( Show ( I i Double )
, Show ( I i 2 )
, Show ( D k 1 ( I i Double ) )
, Show ( D k 1 ( I i 2 ) )
, Show ( D k 1 ( Corner ( I i 2 ) ) )
)
=> Show ( CornerStrokeDatum k i )
--------------------------------------------------------------------------------
type HasBézier :: Nat -> IVness -> Constraint
class HasBézier k i where
-- | Linear interpolation, as a differentiable function.
line :: forall ( n :: Nat )
. ( Module Double ( T ( ℝ n ) ), Torsor ( T ( ℝ n ) ) ( ℝ n )
, Module ( I i Double ) ( T ( I i n ) ), Torsor ( T ( I i n ) ) ( I i n )
)
=> ( I i 1 -> I i Double )
-> Segment ( I i n ) -> C k ( I i 1 ) ( I i n )
-- | A quadratic Bézier curve, as a differentiable function.
bezier2 :: forall ( n :: Nat )
. ( Module Double ( T ( ℝ n ) ), Torsor ( T ( ℝ n ) ) ( ℝ n )
, Module ( I i Double ) ( T ( I i n ) ), Torsor ( T ( I i n ) ) ( I i n )
, Representable Double ( ℝ n )
, Representable 𝕀 ( 𝕀ℝ n )
)
=> ( I i 1 -> I i Double )
-> Quadratic.Bezier ( I i n ) -> C k ( I i 1 ) ( I i n )
-- | A cubic Bézier curve, as a differentiable function.
bezier3 :: forall ( n :: Nat )
. ( Module Double ( T ( ℝ n ) ), Torsor ( T ( ℝ n ) ) ( ℝ n )
, Module ( I i Double ) ( T ( I i n ) ), Torsor ( T ( I i n ) ) ( I i n )
, Representable Double ( ℝ n )
, Representable 𝕀 ( 𝕀ℝ n )
)
=> ( I i 1 -> I i Double )
-> Cubic.Bezier ( I i n ) -> C k ( I i 1 ) ( I i n )
type HasEnvelopeEquation :: Nat -> Constraint
class HasEnvelopeEquation k where
uncurryD :: D k 1 ( D k 1 b ) -> D k 2 b
-- | 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
. ( Module ( I i Double ) ( T ( I i 1 ) )
, Cross ( I i Double ) ( T ( I i 2 ) )
, Transcendental ( I i Double )
, Representable ( I i Double ) ( I i 2 )
, RepDim ( I i 2 ) ~ 2
, Show ( StrokeDatum k i )
)
=> SingIVness i
-> D k 1 ( I i 2 )
-> Unrotated ( D k 2 ( I i 2 ) )
-> Maybe ( D k 1 ( I i Double ) )
-> StrokeDatum k i
cornerEnvelopeEquation
:: ( Module ( I i Double ) ( T ( I i 1 ) )
, Cross ( I i Double ) ( T ( I i 2 ) )
, Transcendental ( I i Double )
, Representable ( I i Double ) ( I i 2 )
, RepDim ( I i 2 ) ~ 2
, Show ( CornerStrokeDatum k i )
)
=> D k 1 ( I i 2 )
-> Unrotated ( D k 1 ( Corner ( I i 2 ) ) )
-> Maybe ( D k 1 ( I i Double ) )
-> CornerStrokeDatum k i
instance HasBézier 2 NonIV where
line co ( Segment a b :: Segment b ) =
D \ ( co -> t ) ->
D21 ( lerp @( T b ) t a b )
( a --> b )
origin
{-# INLINEABLE line #-}
bezier2 co ( bez :: Quadratic.Bezier b ) =
D \ ( co -> t ) ->
D21 ( Quadratic.bezier @( T b ) bez t )
( Quadratic.bezier' bez t )
( Quadratic.bezier'' bez )
{-# INLINEABLE bezier2 #-}
bezier3 co ( bez :: Cubic.Bezier b ) =
D \ ( co -> t ) ->
D21 ( Cubic.bezier @( T b ) bez t )
( Cubic.bezier' bez t )
( Cubic.bezier'' bez t )
{-# INLINEABLE bezier3 #-}
instance HasEnvelopeEquation 2 where
uncurryD = uncurryD2
{-# INLINEABLE uncurryD #-}
envelopeEquation ( ivness :: SingIVness i )
dp@( D21 ( T -> p ) p_t p_tt )
db@( Unrotated ( D22 ( T -> b ) b_t b_s b_tt b_ts b_ss ) )
mbRotation =
StrokeDatum
{ dpath = dp
, dbrush = db
, mbRotation
, stroke = c
, du, dv, ee, 𝛿E𝛿s_unrot_dcdt
}
where
co :: I i Double -> I i 1
co = case ivness of
SNonIV -> coerce
SIsIV -> coerce
(ee, 𝛿E𝛿s_unrot_dcdt) =
let
D12 (T -> u) u_t u_s = du
D12 (T -> v) v_t v_s = dv
ee_val, ee_t, ee_s :: I i Double
ee_val= u × v
ee_t = u_t × v
+ u × v_t
ee_s = u_s × v
+ u × v_s
𝛿E𝛿s_unrot_dcdt_val = ee_s *^ u ^-^ ee_t *^ v
in ( D12
( co ee_val )
( T $ co ee_t ) ( T $ co ee_s )
, Unrotated $ D0 𝛿E𝛿s_unrot_dcdt_val
)
(c, du, dv) = case mbRotation of
Nothing ->
-- c(t,s) = p(t) + b(t,s)
( unT $ p ^+^ b
, D12 ( unT $ p_t ^+^ b_t )
( p_tt ^+^ b_tt ) b_ts
, D12 ( unT $ b_s )
b_ts b_ss
)
Just ( D21 θ ( T θ_t ) ( T θ_tt ) ) ->
-- c(t,s) = p(t) + R(θ(t)) b(t,s)
-- E = ∂c/∂t × ∂c/ds
-- = ( R(-θ(t)) ∂c/∂t ) × ( R(-θ(t)) ∂c/ds )
-- = ( R(-θ(t)) p'(t) + θ'(t) S b(t,s) + ∂b/∂t ) × ∂b/ds
--
let rot, rot' :: T ( I i 2 ) -> T ( I i 2 )
cosθ = cos θ
sinθ = sin θ
-- rot = R(-θ), rot' = R'(-θ)
-- NB: rot' is not the derivative of f(θ) = R(-θ)
rot = rotate cosθ -sinθ
rot' = rotate sinθ cosθ
rot90 = rotate 0 1
u, v, u_t, u_s, v_t, v_s :: T ( I i 2 )
u = rot p_t ^+^ θ_t *^ rot90 b ^+^ b_t
v = b_s
u_t = ( -θ_t *^ rot' p_t
^+^ rot p_tt
)
^+^
( θ_tt *^ rot90 b
^+^ θ_t *^ rot90 b_t
)
^+^ b_tt
u_s = θ_t *^ rot90 b_s ^+^ b_ts
v_t = b_ts
v_s = b_ss
in ( unT $ p ^+^ rotate cosθ sinθ b
, D12 ( unT u ) u_t u_s
, D12 ( unT v ) v_t v_s
)
{-# INLINEABLE envelopeEquation #-}
cornerEnvelopeEquation
dpath@( D21 ( T -> p ) p' p'' )
dcorner@(
Unrotated
( D21
( Corner ( T -> b ) t_s t_e )
( T ( Corner ( T -> b' ) t'_s t'_e ) )
( T ( Corner ( T -> b'' ) t''_s t''_e ) )
) )
mbRotation
=
let
dstroke = case mbRotation of
Nothing ->
D21 ( Corner ( unT $ p ^+^ b ) t_s t_e )
( T $ Corner ( unT $ p' ^+^ b' ) t'_s t'_e )
( T $ Corner ( unT $ p'' ^+^ b'' ) t''_s t''_e )
Just ( D21 θ ( T θ' ) ( T θ'' ) ) ->
let cosθ = Ring.cos θ
sinθ = Ring.sin θ
rot = rotate cosθ sinθ
rot90 = rotate 0 1
rot' = rot . rot90
-- TODO: I shouldn't have to do this by hand.
--
-- c(t) = p(t) + rot(θ(t)) b(t)
--
-- abbreviate; let R = rot(θ(t)), R'=rot'(θ(t)), ...
--
-- c = p + R b
-- c' = p' + R b' + θ' R' b
-- c'' = p'' + R b'' + 2 θ' R' b' + θ'² R'' b + θ'' R' b
d0 v = rot v
d1 v v' = rot v' ^+^ θ' *^ rot' v
d2 v v' v'' = rot v''
^+^ ( 2 * θ' ) *^ rot' v'
^-^ ( θ' Ring.^ 2 ) *^ rot v
^+^ θ'' *^ rot' v
in
D21
( Corner
( unT $ p ^+^ d0 b )
( d0 t_s )
( d0 t_e )
)
( T $ Corner
( unT $ p' ^+^ d1 b b' )
( d1 t_s t'_s )
( d1 t_e t'_e )
)
( T $ Corner
( unT $ p'' ^+^ d2 b b' b'' )
( d2 t_s t'_s t''_s )
( d2 t_e t'_e t''_e )
)
in
CornerStrokeDatum
{ dpath, dcorner, dstroke, mbRotation }
{-# INLINEABLE cornerEnvelopeEquation #-}
instance HasBézier 3 NonIV where
line co ( Segment a b :: Segment b ) =
D \ ( co -> t ) ->
D31 ( lerp @( T b ) t a b )
( a --> b )
origin
origin
{-# INLINEABLE line #-}
bezier2 co ( bez :: Quadratic.Bezier b ) =
D \ ( co -> t ) ->
D31 ( Quadratic.bezier @( T b ) bez t )
( Quadratic.bezier' bez t )
( Quadratic.bezier'' bez )
origin
{-# INLINEABLE bezier2 #-}
bezier3 co ( bez :: Cubic.Bezier b ) =
D \ ( co -> t ) ->
D31 ( Cubic.bezier @( T b ) bez t )
( Cubic.bezier' bez t )
( Cubic.bezier'' bez t )
( Cubic.bezier''' bez )
{-# INLINEABLE bezier3 #-}
instance HasEnvelopeEquation 3 where
uncurryD = uncurryD3
{-# INLINEABLE uncurryD #-}
envelopeEquation ( ivness :: SingIVness i )
dp@( D31 ( T -> p ) p_t p_tt p_ttt )
db@( Unrotated ( D32 ( T -> b ) b_t b_s b_tt b_ts b_ss b_ttt b_tts b_tss b_sss ) )
mbRotation =
StrokeDatum
{ dpath = dp
, dbrush = db
, mbRotation
, stroke = c
, du, dv, ee, 𝛿E𝛿s_unrot_dcdt
}
where
co :: I i Double -> I i 1
co = case ivness of
SNonIV -> coerce
SIsIV -> coerce
(ee, 𝛿E𝛿s_unrot_dcdt) =
let
D22 (T -> u) u_t u_s u_tt u_ts u_ss = du
D22 (T -> v) v_t v_s v_tt v_ts v_ss = dv
ee_val, ee_t, ee_s, ee_tt, ee_ts, ee_ss :: I i Double
ee_val= u × v
ee_t = u_t × v
+ u × v_t
ee_s = u_s × v
+ u × v_s
ee_tt = u_tt × v
+ 2 * ( u_t × v_t )
+ u × v_tt
ee_ts = u_ts × v
+ u_t × v_s
+ u_s × v_t
+ u × v_ts
ee_ss = u_ss × v
+ 2 * ( u_s × v_s )
+ u × v_ss
𝛿E𝛿s_unrot_dcdt_val =
ee_s *^ u ^-^ ee_t *^ v
𝛿E𝛿s_unrot_dcdt_t =
ee_ts *^ u ^+^ ee_s *^ u_t
^-^ ( ee_tt *^ v ^+^ ee_t *^ v_t )
𝛿E𝛿s_unrot_dcdt_s =
ee_ss *^ u ^+^ ee_s *^ u_s
^-^ ( ee_ts *^ v ^+^ ee_t *^ v_s )
in ( D22
( co ee_val )
( T $ co ee_t ) ( T $ co ee_s )
( T $ co ee_tt) ( T $ co ee_ts ) ( T $ co ee_ss )
, Unrotated $ D12 𝛿E𝛿s_unrot_dcdt_val ( T 𝛿E𝛿s_unrot_dcdt_t ) ( T 𝛿E𝛿s_unrot_dcdt_s )
)
(c, du, dv) = case mbRotation of
Nothing ->
-- c(t,s) = p(t) + b(t,s)
( unT $ p ^+^ b
, D22 ( unT $ p_t ^+^ b_t )
( p_tt ^+^ b_tt ) b_ts
( p_ttt ^+^ b_ttt ) b_tts b_tss
, D22 ( unT $ b_s )
b_ts b_ss
b_tts b_tss b_sss
)
Just ( D31 θ ( T θ_t ) ( T θ_tt ) ( T θ_ttt ) ) ->
-- c(t,s) = p(t) + R(θ(t)) b(t,s)
-- E = ∂c/∂t × ∂c/ds
-- = ( R(-θ(t)) ∂c/∂t ) × ( R(-θ(t)) ∂c/ds )
-- = ( R(-θ(t)) p'(t) + θ'(t) S b(t,s) + ∂b/∂t ) × ∂b/ds
--
let rot, rot' :: T ( I i 2 ) -> T ( I i 2 )
cosθ = cos θ
sinθ = sin θ
-- rot = R(-θ), rot' = R'(-θ), rot'' = R''(-θ)
-- NB: rot' is not the derivative of f(θ) = R(-θ)
rot = rotate cosθ -sinθ
rot' = rotate sinθ cosθ
rot90 = rotate 0 1
u, v, u_t, u_s, u_tt, u_ts, u_ss, v_t, v_s, v_tt, v_ts, v_ss :: T ( I i 2 )
u = rot p_t ^+^ θ_t *^ rot90 b ^+^ b_t
v = b_s
u_t = ( -θ_t *^ rot' p_t
^+^ rot p_tt
)
^+^
( θ_tt *^ rot90 b
^+^ θ_t *^ rot90 b_t
)
^+^ b_tt
u_s = θ_t *^ rot90 b_s ^+^ b_ts
u_tt = ( -( θ_t ^ 2 ) *^ rot p_t
^-^ θ_tt *^ rot' p_t
^-^ ( 2 * θ_t ) *^ rot' p_tt
^+^ rot p_ttt
)
^+^ ( θ_ttt *^ rot90 b
^+^ ( 2 * θ_tt ) *^ rot90 b_t
^+^ θ_t *^ rot90 b_tt
)
^+^ b_ttt
u_ts = θ_tt *^ rot90 b_s
^+^ θ_t *^ rot90 b_ts
^+^ b_tts
u_ss = θ_t *^ rot90 b_ss
^+^ b_tss
v_t = b_ts
v_s = b_ss
v_tt = b_tts
v_ts = b_tss
v_ss = b_sss
in ( unT $ p ^+^ rotate cosθ sinθ b
, D22 ( unT u ) u_t u_s u_tt u_ts u_ss
, D22 ( unT v ) v_t v_s v_tt v_ts v_ss
)
{-# INLINEABLE envelopeEquation #-}
cornerEnvelopeEquation
dpath@( D31 ( T -> p ) p' p'' p''' )
dcorner@(
Unrotated
( D31
( Corner ( T -> b ) t_s t_e )
( T ( Corner ( T -> b' ) t'_s t'_e ) )
( T ( Corner ( T -> b'' ) t''_s t''_e ) )
( T ( Corner ( T -> b''' ) t'''_s t'''_e ) )
) )
mbRotation =
let
dstroke = case mbRotation of
Nothing ->
D31 ( Corner ( unT $ p ^+^ b ) t_s t_e )
( T $ Corner ( unT $ p' ^+^ b' ) t'_s t'_e )
( T $ Corner ( unT $ p'' ^+^ b'' ) t''_s t''_e )
( T $ Corner ( unT $ p''' ^+^ b''' ) t'''_s t'''_e )
Just ( D31 θ ( T θ' ) ( T θ'' ) ( T θ''' ) ) ->
let cosθ = Ring.cos θ
sinθ = Ring.sin θ
rot90 = rotate 0 1
rot = rotate cosθ sinθ
rot' = rot . rot90
-- TODO: I shouldn't have to do this by hand.
--
-- c(t) = p(t) + rot(θ(t)) b(t)
--
-- abbreviate; let R = rot(θ(t)), R'=rot'(θ(t)), ...
--
-- c = p + R b
-- c' = p' + R b' + θ' R' b
-- c'' = p'' + R b'' + 2 θ' R' b' + θ'² R'' b + θ'' R' b
-- c''' = p''' + R b''' + 3 θ' R'' ( θ' b' + θ'' b )
-- + R' ( θ''' b + 3 θ'' b' + 3 θ' b'' )
-- + θ'³ R''' b
d0 v = rot v
d1 v v' = rot v' ^+^ θ' *^ rot' v
d2 v v' v'' = rot v''
^+^ ( 2 * θ' ) *^ rot' v'
^-^ ( θ' Ring.^ 2 ) *^ rot v
^+^ θ'' *^ rot' v
d3 v v' v'' v'''
= rot v'''
^-^ ( 3 * θ' ) *^ ( rot ( θ' *^ v' ^+^ θ'' *^ v ) )
^+^ rot' ( θ''' *^ v ^+^ ( 3 * θ'' ) *^ v' ^+^ ( 3 * θ' ) *^ v'' )
^-^ ( θ' Ring.^ 3 ) *^ rot' v
in
D31
( Corner
( unT $ p ^+^ d0 b )
( d0 t_s )
( d0 t_e )
)
( T $ Corner
( unT $ p' ^+^ d1 b b' )
( d1 t_s t'_s )
( d1 t_e t'_e )
)
( T $ Corner
( unT $ p'' ^+^ d2 b b' b'' )
( d2 t_s t'_s t''_s )
( d2 t_e t'_e t''_e )
)
( T $ Corner
( unT $ p''' ^+^ d3 b b' b'' b''' )
( d3 t_s t'_s t''_s t'''_s )
( d3 t_e t'_e t''_e t'''_e )
)
in
CornerStrokeDatum
{ dpath, dcorner, dstroke, mbRotation }
{-# INLINEABLE cornerEnvelopeEquation #-}
instance HasBézier 3 IsIV where
line co ( Segment a b :: Segment b ) =
D \ ( co -> t ) ->
D31 ( lerp @( T b ) t a b )
( a --> b )
origin
origin
{-# INLINEABLE line #-}
bezier2 co ( bez :: Quadratic.Bezier b ) =
D \ ( co -> t ) ->
D31 ( aabb bez ( `evaluateQuadratic` t ) )
( Quadratic.bezier' bez t )
( Quadratic.bezier'' bez )
origin
{-# INLINEABLE bezier2 #-}
bezier3 co ( bez :: Cubic.Bezier b ) =
D \ ( co -> t ) ->
D31 ( aabb bez ( `evaluateCubic` t ) )
( T $ aabb ( fmap unT $ Cubic.derivative ( fmap T bez ) ) ( `evaluateQuadratic` t ) )
( Cubic.bezier'' bez t )
( Cubic.bezier''' bez )
{-# INLINEABLE bezier3 #-}
{- Note [Computing Béziers over intervals]
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Here's how we evaluate a Bézier curve when the coefficients are intervals.
As we are using axis-aligned interval arithmetic, we reduce to the 1D situation.
Now, the formulas for Bézier curves (with value of the t parameter in [0,1])
are convex combinations
b_0(t) [p_0,q_0] + b_1(t) [p_1,q_1] + ... + b_n(t) [p_n,q_n]
-- Here b_0, ..., b_n are Bernstein polynomials.
This means that the minimum value attained by the Bézier curve as we vary
both the time parameter and the values of the points within their respective
intervals will occur when we take
b_0(t) p_0 + b_1(t) p_1 + ... + b_n(t) p_n
and the maximum when we take
b_0(t) q_0 + b_1(t) q_1 + ... + b_n(t) q_n
For each of these, we can compute a minimum/maximum by computing an axis-aligned
bounding box for the Bézier curve. This is done by computing the derivative
with respect to t and root-finding.
-}
-- | Evaluate a cubic Bézier curve, when both the coefficients and the
-- parameter are intervals.
evaluateCubic :: Cubic.Bezier 𝕀 -> 𝕀 -> 𝕀
evaluateCubic bez t =
-- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t ⊆ [0,1]
let inf_bez = fmap inf bez
sup_bez = fmap sup bez
mins = fmap (Cubic.bezier @( T Double ) inf_bez)
$ inf t :| ( sup t : filter ( ∈ t ) ( Cubic.extrema inf_bez ) )
maxs = fmap (Cubic.bezier @( T Double ) sup_bez)
$ inf t :| ( sup t : filter ( ∈ t ) ( Cubic.extrema sup_bez ) )
in 𝕀 ( minimum mins ) ( maximum maxs )
{-# INLINEABLE evaluateCubic #-}
-- | Evaluate a quadratic Bézier curve, when both the coefficients and the
-- parameter are intervals.
evaluateQuadratic :: Quadratic.Bezier 𝕀 -> 𝕀 -> 𝕀
evaluateQuadratic bez t =
-- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t ⊆ [0,1]
let inf_bez = fmap inf bez
sup_bez = fmap sup bez
mins = fmap (Quadratic.bezier @( T Double ) inf_bez)
$ inf t :| ( sup t : filter ( ∈ t ) ( Quadratic.extrema inf_bez ) )
maxs = fmap (Quadratic.bezier @( T Double ) sup_bez)
$ inf t :| ( sup t : filter ( ∈ t ) ( Quadratic.extrema sup_bez ) )
in 𝕀 ( minimum mins ) ( maximum maxs )
{-# INLINEABLE evaluateQuadratic #-}