compute roots of envelope equation

This commit is contained in:
sheaf 2021-05-15 23:27:21 +02:00
parent df2469487d
commit 7431e8ba67
2 changed files with 456 additions and 0 deletions

View file

@ -83,6 +83,7 @@ library
exposed-modules:
Math.Bezier.Cubic
, Math.Bezier.Cubic.Fit
, Math.Bezier.Envelope
, Math.Bezier.Quadratic
, Math.Bezier.Spline
, Math.Bezier.Stroke

View file

@ -0,0 +1,455 @@
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
module Math.Bezier.Envelope where
-- acts
import Data.Act
( Torsor((-->)) )
-- deepseq
import Control.DeepSeq
( NFData )
-- primitive
import Data.Primitive.Types
( Prim )
-- MetaBrush
import Math.Roots
( realRoots )
import qualified Math.Bezier.Cubic as Cubic
( Bezier(..), bezier, bezier' )
import qualified Math.Bezier.Quadratic as Quadratic
( Bezier(..), bezier, bezier' )
import Math.Module
( Module((^+^),(*^)), lerp, cross )
import Math.Vector2D
( Point2D(..), Vector2D(..), Segment(..) )
--------------------------------------------------------------------------------
-- | Find the roots of the envelope equation for a family of cubic Bézier curves
-- varying along a cubic Bézier path.
--
-- \[ c(t,u) = p(t) + b(t,u), \]
--
-- where \( t \mapsto p(t) \) describes the underlying path,
-- and \( u \mapsto b(t_0,u) \) describes the brush shape at point \( t = t_0 \).
--
-- The envelope equation is then:
--
-- \[ \frac{\partial c}{\partial t} \cross \frac{\partial c}{\partial u} = 0. \]
--
-- Given \( t_0 \), this function returns a (possibly empty) list of values \( u_i \)
-- satisfying the envelope equation at \( (t_0, u_i) \).
--
-- The points \( c(t_0,u_i) \) are thus potential outline points on the contour stroked
-- by the brush as it moves along the path.
envelope33
:: forall r
. ( RealFloat r, Prim r, NFData r )
=> Cubic.Bezier ( Point2D r ) -> Cubic.Bezier ( Cubic.Bezier ( Point2D r ) ) -> r -> [ r ]
envelope33 path
( Cubic.Bezier
( Cubic.Bezier b00 b01 b02 b03 )
( Cubic.Bezier b10 b11 b12 b13 )
( Cubic.Bezier b20 b21 b22 b23 )
( Cubic.Bezier b30 b31 b32 b33 )
) t0 = realRoots 50 [ a5, a4, a3, a2, a1, a0 ]
where
-- Compute ∂p/∂t(t0).
dpdt :: Vector2D r
dpdt = Cubic.bezier' @( Vector2D r ) path t0
-- Compute ∂b/∂t(t0,u) using the Bernstein basis:
--
-- ∂b/∂t(t0,u) = Cubic.bezier ( Cubic.Bezier ct0 ct1 ct2 ct3 ) u.
ct0, ct1, ct2, ct3, dt0, dt1, dt2, dt3 :: Vector2D r
ct0 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b00 b10 b20 b30 ) t0
ct1 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b01 b11 b21 b31 ) t0
ct2 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b02 b12 b22 b32 ) t0
ct3 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b03 b13 b23 b33 ) t0
-- Add ∂p/∂t and convert the Bernstein representation to the monomial basis to obtain
--
-- ∂c/∂t(t0,u) = dt0 + u dt1 + u² dt2 + u³ dt3.
dt0 = ct0 ^+^ dpdt
dt1 = 3 *^ ( ct0 --> ct1 )
dt2 = 3 *^ ( ct1 --> ct0 ^+^ ct1 --> ct2 )
dt3 = ct0 --> ct3 ^+^ 3 *^ ( ct2 --> ct1 )
-- Compute ∂c/∂u(t0,u) using the Bernstein basis:
--
-- ∂c/∂u(t0,u) = Cubic.bezier' ( Cubic.Bezier cu0 cu1 cu2 cu3 ) u.
cu0, cu1, cu2, cu3 :: Point2D r
cu0 = Cubic.bezier @( Vector2D r ) ( Cubic.Bezier b00 b10 b20 b30 ) t0
cu1 = Cubic.bezier @( Vector2D r ) ( Cubic.Bezier b01 b11 b21 b31 ) t0
cu2 = Cubic.bezier @( Vector2D r ) ( Cubic.Bezier b02 b12 b22 b32 ) t0
cu3 = Cubic.bezier @( Vector2D r ) ( Cubic.Bezier b03 b13 b23 b33 ) t0
-- Convert the Bernstein representation to the monomial basis to obtain
--
-- ∂c/∂u(t0,u) = du0 + u du1 + u² du2.
du0, du1, du2 :: Vector2D r
du0 = cu0 --> cu1
du1 = 2 *^ ( cu1 --> cu0 ^+^ cu1 --> cu2 )
du2 = cu0 --> cu3 ^+^ 3 *^ ( cu2 --> cu1 )
-- Expand out the cross-product ∂c/∂t × ∂c/∂u to obtain the envelope equation:
--
-- a0 + a1 u + a2 u² + a3 u³ + a4 u⁴ + a5 u⁵.
a0, a1, a2, a3, a4, a5 :: r
a0 = dt0 `cross` du0
a1 = dt1 `cross` du0 + dt0 `cross` du1
a2 = dt2 `cross` du0 + dt1 `cross` du1 + dt0 `cross` du2
a3 = dt3 `cross` du0 + dt2 `cross` du1 + dt1 `cross` du2
a4 = dt3 `cross` du1 + dt2 `cross` du2
a5 = dt3 `cross` du2
-- | Find the roots of the envelope equation for a family of cubic Bézier curves
-- varying along a quadratic Bézier path.
--
-- See 'envelope33' for more information.
envelope23
:: forall r
. ( RealFloat r, Prim r, NFData r )
=> Quadratic.Bezier ( Point2D r ) -> Quadratic.Bezier ( Cubic.Bezier ( Point2D r ) ) -> r -> [ r ]
envelope23 path
( Quadratic.Bezier
( Cubic.Bezier b00 b01 b02 b03 )
( Cubic.Bezier b10 b11 b12 b13 )
( Cubic.Bezier b20 b21 b22 b23 )
) t0 = realRoots 50 [ a5, a4, a3, a2, a1, a0 ]
where
dpdt :: Vector2D r
dpdt = Quadratic.bezier' @( Vector2D r ) path t0
ct0, ct1, ct2, ct3, dt0, dt1, dt2, dt3 :: Vector2D r
ct0 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b00 b10 b20 ) t0
ct1 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b01 b11 b21 ) t0
ct2 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b02 b12 b22 ) t0
ct3 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b03 b13 b23 ) t0
dt0 = ct0 ^+^ dpdt
dt1 = 3 *^ ( ct0 --> ct1 )
dt2 = 3 *^ ( ct1 --> ct0 ^+^ ct1 --> ct2 )
dt3 = ct0 --> ct3 ^+^ 3 *^ ( ct2 --> ct1 )
cu0, cu1, cu2, cu3 :: Point2D r
cu0 = Quadratic.bezier @( Vector2D r ) ( Quadratic.Bezier b00 b10 b20 ) t0
cu1 = Quadratic.bezier @( Vector2D r ) ( Quadratic.Bezier b01 b11 b21 ) t0
cu2 = Quadratic.bezier @( Vector2D r ) ( Quadratic.Bezier b02 b12 b22 ) t0
cu3 = Quadratic.bezier @( Vector2D r ) ( Quadratic.Bezier b03 b13 b23 ) t0
du0, du1, du2 :: Vector2D r
du0 = cu0 --> cu1
du1 = 2 *^ ( cu1 --> cu0 ^+^ cu1 --> cu2 )
du2 = cu0 --> cu3 ^+^ 3 *^ ( cu2 --> cu1 )
a0, a1, a2, a3, a4, a5 :: r
a0 = dt0 `cross` du0
a1 = dt1 `cross` du0 + dt0 `cross` du1
a2 = dt2 `cross` du0 + dt1 `cross` du1 + dt0 `cross` du2
a3 = dt3 `cross` du0 + dt2 `cross` du1 + dt1 `cross` du2
a4 = dt3 `cross` du1 + dt2 `cross` du2
a5 = dt3 `cross` du2
-- | Find the roots of the envelope equation for a family of cubic Bézier curves
-- varying along a straight line path.
--
-- See 'envelope33' for more information.
envelope13
:: forall r
. ( RealFloat r, Prim r, NFData r )
=> Segment ( Point2D r ) -> Segment ( Cubic.Bezier ( Point2D r ) ) -> r -> [ r ]
envelope13 ( Segment p0 p1 )
( Segment
( Cubic.Bezier b00 b01 b02 b03 )
( Cubic.Bezier b10 b11 b12 b13 )
) t0 = realRoots 50 [ a5, a4, a3, a2, a1, a0 ]
where
dpdt :: Vector2D r
dpdt = p0 --> p1
ct0, ct1, ct2, ct3, dt0, dt1, dt2, dt3 :: Vector2D r
ct0 = b00 --> b10
ct1 = b01 --> b11
ct2 = b02 --> b12
ct3 = b03 --> b13
dt0 = ct0 ^+^ dpdt
dt1 = 3 *^ ( ct0 --> ct1 )
dt2 = 3 *^ ( ct1 --> ct0 ^+^ ct1 --> ct2 )
dt3 = ct0 --> ct3 ^+^ 3 *^ ( ct2 --> ct1 )
cu0, cu1, cu2, cu3 :: Point2D r
cu0 = lerp @( Vector2D r ) t0 b00 b10
cu1 = lerp @( Vector2D r ) t0 b01 b11
cu2 = lerp @( Vector2D r ) t0 b02 b12
cu3 = lerp @( Vector2D r ) t0 b03 b13
du0, du1, du2 :: Vector2D r
du0 = cu0 --> cu1
du1 = 2 *^ ( cu1 --> cu0 ^+^ cu1 --> cu2 )
du2 = cu0 --> cu3 ^+^ 3 *^ ( cu2 --> cu1 )
a0, a1, a2, a3, a4, a5 :: r
a0 = dt0 `cross` du0
a1 = dt1 `cross` du0 + dt0 `cross` du1
a2 = dt2 `cross` du0 + dt1 `cross` du1 + dt0 `cross` du2
a3 = dt3 `cross` du0 + dt2 `cross` du1 + dt1 `cross` du2
a4 = dt3 `cross` du1 + dt2 `cross` du2
a5 = dt3 `cross` du2
-- | Find the roots of the envelope equation for a family of quadratic Bézier curves
-- varying along a cubic Bézier path.
--
-- See 'envelope33' for more information.
envelope32
:: forall r
. ( RealFloat r, Prim r, NFData r )
=> Cubic.Bezier ( Point2D r ) -> Cubic.Bezier ( Quadratic.Bezier ( Point2D r ) ) -> r -> [ r ]
envelope32 path
( Cubic.Bezier
( Quadratic.Bezier b00 b01 b02 )
( Quadratic.Bezier b10 b11 b12 )
( Quadratic.Bezier b20 b21 b22 )
( Quadratic.Bezier b30 b31 b32 )
) t0 = realRoots 50 [ a3, a2, a1, a0 ]
where
dpdt :: Vector2D r
dpdt = Cubic.bezier' @( Vector2D r ) path t0
ct0, ct1, ct2, dt0, dt1, dt2 :: Vector2D r
ct0 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b00 b10 b20 b30 ) t0
ct1 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b01 b11 b21 b31 ) t0
ct2 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b02 b12 b22 b32 ) t0
dt0 = ct0 ^+^ dpdt
dt1 = 2 *^ ( ct0 --> ct1 )
dt2 = ct1 --> ct0 ^+^ ct1 --> ct2
cu0, cu1, cu2 :: Point2D r
cu0 = Cubic.bezier @( Vector2D r ) ( Cubic.Bezier b00 b10 b20 b30 ) t0
cu1 = Cubic.bezier @( Vector2D r ) ( Cubic.Bezier b01 b11 b21 b31 ) t0
cu2 = Cubic.bezier @( Vector2D r ) ( Cubic.Bezier b02 b12 b22 b32 ) t0
du0, du1 :: Vector2D r
du0 = cu0 --> cu1
du1 = cu1 --> cu0 ^+^ cu1 --> cu2
a0, a1, a2, a3 :: r
a0 = dt0 `cross` du0
a1 = dt1 `cross` du0 + dt0 `cross` du1
a2 = dt2 `cross` du0 + dt1 `cross` du1
a3 = dt2 `cross` du1
-- | Find the roots of the envelope equation for a family of quadratic Bézier curves
-- varying along a quadratic Bézier path.
--
-- See 'envelope33' for more information.
envelope22
:: forall r
. ( RealFloat r, Prim r, NFData r )
=> Quadratic.Bezier ( Point2D r ) -> Quadratic.Bezier ( Quadratic.Bezier ( Point2D r ) ) -> r -> [ r ]
envelope22 path
( Quadratic.Bezier
( Quadratic.Bezier b00 b01 b02 )
( Quadratic.Bezier b10 b11 b12 )
( Quadratic.Bezier b20 b21 b22 )
) t0 = realRoots 50 [ a3, a2, a1, a0 ]
where
dpdt :: Vector2D r
dpdt = Quadratic.bezier' @( Vector2D r ) path t0
ct0, ct1, ct2, dt0, dt1, dt2 :: Vector2D r
ct0 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b00 b10 b20 ) t0
ct1 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b01 b11 b21 ) t0
ct2 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b02 b12 b22 ) t0
dt0 = ct0 ^+^ dpdt
dt1 = 2 *^ ( ct0 --> ct1 )
dt2 = ct1 --> ct0 ^+^ ct1 --> ct2
cu0, cu1, cu2 :: Point2D r
cu0 = Quadratic.bezier @( Vector2D r ) ( Quadratic.Bezier b00 b10 b20 ) t0
cu1 = Quadratic.bezier @( Vector2D r ) ( Quadratic.Bezier b01 b11 b21 ) t0
cu2 = Quadratic.bezier @( Vector2D r ) ( Quadratic.Bezier b02 b12 b22 ) t0
du0, du1 :: Vector2D r
du0 = cu0 --> cu1
du1 = cu1 --> cu0 ^+^ cu1 --> cu2
a0, a1, a2, a3 :: r
a0 = dt0 `cross` du0
a1 = dt1 `cross` du0 + dt0 `cross` du1
a2 = dt2 `cross` du0 + dt1 `cross` du1
a3 = dt2 `cross` du1
-- | Find the roots of the envelope equation for a family of quadratic Bézier curves
-- varying along a straight line.
--
-- See 'envelope33' for more information.
envelope12
:: forall r
. ( RealFloat r, Prim r, NFData r )
=> Segment ( Point2D r ) -> Segment ( Quadratic.Bezier ( Point2D r ) ) -> r -> [ r ]
envelope12 ( Segment p0 p1 )
( Segment
( Quadratic.Bezier b00 b01 b02 )
( Quadratic.Bezier b10 b11 b12 )
) t0 = realRoots 50 [ a3, a2, a1, a0 ]
where
dpdt :: Vector2D r
dpdt = p0 --> p1
ct0, ct1, ct2, dt0, dt1, dt2 :: Vector2D r
ct0 = b00 --> b10
ct1 = b01 --> b11
ct2 = b02 --> b12
dt0 = ct0 ^+^ dpdt
dt1 = 2 *^ ( ct0 --> ct1 )
dt2 = ct1 --> ct0 ^+^ ct1 --> ct2
cu0, cu1, cu2 :: Point2D r
cu0 = lerp @( Vector2D r ) t0 b00 b10
cu1 = lerp @( Vector2D r ) t0 b01 b11
cu2 = lerp @( Vector2D r ) t0 b02 b12
du0, du1 :: Vector2D r
du0 = cu0 --> cu1
du1 = cu1 --> cu0 ^+^ cu1 --> cu2
a0, a1, a2, a3 :: r
a0 = dt0 `cross` du0
a1 = dt1 `cross` du0 + dt0 `cross` du1
a2 = dt2 `cross` du0 + dt1 `cross` du1
a3 = dt2 `cross` du1
-- | Find the roots of the envelope equation for a family of line segments
-- varying along a cubic Bézier curve.
--
-- See 'envelope33' for more information.
envelope31
:: forall r
. ( RealFloat r, Prim r, NFData r )
=> Cubic.Bezier ( Point2D r ) -> Cubic.Bezier ( Segment ( Point2D r ) ) -> r -> [ r ]
envelope31 path
( Cubic.Bezier
( Segment b00 b01 )
( Segment b10 b11 )
( Segment b20 b21 )
( Segment b30 b31 )
) t0 = [ - a1 / a0 ]
where
dpdt :: Vector2D r
dpdt = Cubic.bezier' @( Vector2D r ) path t0
ct0, ct1, dt0, dt1 :: Vector2D r
ct0 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b00 b10 b20 b30 ) t0
ct1 = Cubic.bezier' @( Vector2D r ) ( Cubic.Bezier b01 b11 b21 b31 ) t0
dt0 = ct0 ^+^ dpdt
dt1 = ct0 --> ct1
cu0, cu1 :: Point2D r
cu0 = Cubic.bezier @( Vector2D r ) ( Cubic.Bezier b00 b10 b20 b30 ) t0
cu1 = Cubic.bezier @( Vector2D r ) ( Cubic.Bezier b01 b11 b21 b31 ) t0
du0 :: Vector2D r
du0 = cu0 --> cu1
a0, a1 :: r
a0 = dt0 `cross` du0
a1 = dt1 `cross` du0
-- | Find the roots of the envelope equation for a family of line segments
-- varying along a quadratic Bézier curve.
--
-- See 'envelope33' for more information.
envelope21
:: forall r
. ( RealFloat r, Prim r, NFData r )
=> Quadratic.Bezier ( Point2D r ) -> Quadratic.Bezier ( Segment ( Point2D r ) ) -> r -> [ r ]
envelope21 path
( Quadratic.Bezier
( Segment b00 b01 )
( Segment b10 b11 )
( Segment b20 b21 )
) t0 = [ - a1 / a0 ]
where
dpdt :: Vector2D r
dpdt = Quadratic.bezier' @( Vector2D r ) path t0
ct0, ct1, dt0, dt1 :: Vector2D r
ct0 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b00 b10 b20 ) t0
ct1 = Quadratic.bezier' @( Vector2D r ) ( Quadratic.Bezier b01 b11 b21 ) t0
dt0 = ct0 ^+^ dpdt
dt1 = ct0 --> ct1
cu0, cu1 :: Point2D r
cu0 = Quadratic.bezier @( Vector2D r ) ( Quadratic.Bezier b00 b10 b20 ) t0
cu1 = Quadratic.bezier @( Vector2D r ) ( Quadratic.Bezier b01 b11 b21 ) t0
du0 :: Vector2D r
du0 = cu0 --> cu1
a0, a1 :: r
a0 = dt0 `cross` du0
a1 = dt1 `cross` du0
-- | Find the roots of the envelope equation for a family of line segments
-- varying along a straight line path.
--
-- See 'envelope33' for more information.
envelope11
:: forall r
. ( RealFloat r, Prim r, NFData r )
=> Segment ( Point2D r ) -> Segment ( Segment ( Point2D r ) ) -> r -> [ r ]
envelope11 ( Segment p0 p1 )
( Segment
( Segment b00 b01 )
( Segment b10 b11 )
) t0 = [ - a1 / a0 ]
where
dpdt :: Vector2D r
dpdt = p0 --> p1
ct0, ct1, dt0, dt1 :: Vector2D r
ct0 = b00 --> b10
ct1 = b01 --> b11
dt0 = ct0 ^+^ dpdt
dt1 = ct0 --> ct1
cu0, cu1 :: Point2D r
cu0 = lerp @( Vector2D r ) t0 b00 b10
cu1 = lerp @( Vector2D r ) t0 b01 b11
du0 :: Vector2D r
du0 = cu0 --> cu1
a0, a1 :: r
a0 = dt0 `cross` du0
a1 = dt1 `cross` du0