Cubic.hs

{-# LANGUAGE AllowAmbiguousTypes  #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE UndecidableInstances #-}

module Math.Bezier.Cubic
  ( Bezier(..)
  , fromQuadratic
  , bezier, bezier', bezier'', bezier'''
  , derivative
  , curvature, squaredCurvature, signedCurvature
  , subdivide, restrict
  , ddist, closestPoint
  , drag, selfIntersectionParameters
  , extrema
  , speedCriticalPoints
  )
  where

-- base
import Data.Complex
  ( realPart, Complex )
import Data.Function
  ( on )
import Data.List
  ( sortOn, nubBy )
import Data.List.NonEmpty
  ( NonEmpty(..) )
import Data.Monoid
  ( Ap(..) )
import Data.Semigroup
  ( ArgMin, Min(..), Arg(..) )
import GHC.Generics
  ( Generic, Generic1
  , Generically(..), Generically1(..)
  )

-- acts
import Data.Act
  ( Act(..)
  , Torsor
    ( (-->) )
  )

-- deepseq
import Control.DeepSeq
  ( NFData, NFData1 )

-- fp-ieee
import Numeric.Floating.IEEE.NaN
  ( RealFloatNaN )

-- groups
import Data.Group
  ( Group )

-- groups-generic
import Data.Group.Generics
  ()

-- primitive
import Data.Primitive.Types
  ( Prim )

-- brush-strokes
import qualified Math.Bezier.Quadratic as Quadratic
  ( Bezier(..), bezier )
import Math.Epsilon
  ( epsilon )
import Math.Module
  ( Module (..)
  , lerp
  , Inner((^.^)), norm, squaredNorm
  , Cross((×))
  )
import Math.Roots
  ( realRoots, solveQuadratic, roots )
import Math.Linear
  ( ℝ(..), T(..) )
import qualified Math.Ring as Ring

--------------------------------------------------------------------------------

-- | Points defining a cubic Bézier curve (Bernstein form).
--
-- @ p0 @ and @ p3 @ are endpoints, whereas @ p1 @ and @ p2 @ are control points.
data Bezier p
  = Bezier
    { p0, p1, p2, p3 :: !p }
  deriving stock ( Generic, Generic1, Functor, Foldable, Traversable )
  deriving ( Semigroup, Monoid, Group )
    via Generically ( Bezier p )
  deriving Applicative
    via Generically1 Bezier
  deriving anyclass ( NFData, NFData1 )

deriving via Ap Bezier p
  instance {-# OVERLAPPING #-} Act v p => Act v ( Bezier p )
deriving via Ap Bezier ( T b )
  instance Module r ( T b ) => Module r ( T ( Bezier b ) )

instance Show p => Show (Bezier p) where
  show (Bezier p1 p2 p3 p4) =
    show p1 ++ "--" ++ show p2 ++ "--" ++ show p3 ++ "->" ++ show p4

-- | Degree raising: convert a quadratic Bézier curve to a cubic Bézier curve.
fromQuadratic :: forall v r p. ( Torsor v p, Module r v, Fractional r ) => Quadratic.Bezier p -> Bezier p
fromQuadratic ( Quadratic.Bezier { p0 = q0, p1 = q1, p2 = q2 } ) = Bezier {..}
  where
    !p0 = q0
    !p1 = lerp @v (2/3) q0 q1
    !p2 = lerp @v (1/3) q1 q2
    !p3 = q2
{-# INLINEABLE fromQuadratic #-}

-- | Cubic Bézier curve.
bezier :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> p
bezier ( Bezier {..} ) t =
  lerp @v t
    ( Quadratic.bezier @v ( Quadratic.Bezier p0 p1 p2 ) t )
    ( Quadratic.bezier @v ( Quadratic.Bezier p1 p2 p3 ) t )
{-# INLINEABLE bezier #-}

-- | The derivative of a Cubic Bézier curve, as a quadratic Bézier curve.
derivative :: ( Group v, Module r v ) => Bezier v -> Quadratic.Bezier v
derivative ( Bezier {..} ) = ( Ring.fromInteger 3 *^ )
                          <$> Quadratic.Bezier ( p0 --> p1 ) ( p1 --> p2 ) ( p2 --> p3 )
{-# INLINEABLE derivative #-}

-- | Derivative of a cubic Bézier curve.
bezier' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
bezier' ( Bezier {..} )
  = ( Ring.fromInteger 3 *^ )
  . Quadratic.bezier @v ( Quadratic.Bezier ( p0 --> p1 ) ( p1 --> p2 ) ( p2 --> p3 ) )
{-# INLINEABLE bezier' #-}

-- | Second derivative of a cubic Bézier curve.
bezier'' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
bezier'' ( Bezier {..} ) t
  = ( Ring.fromInteger 6 *^ )
  $ lerp @v t
      ( p1 --> p0 ^+^ p1 --> p2 )
      ( p2 --> p1 ^+^ p2 --> p3 )
{-# INLINEABLE bezier'' #-}

-- | 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 ) )
{-# INLINEABLE bezier''' #-}

-- | 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
{-# INLINEABLE curvature #-}

-- | Square of curvature of a cubic Bézier curve.
squaredCurvature :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> r -> r
squaredCurvature bez t
  | sq_nm_g' < epsilon
  = 1 / 0
  | otherwise
  = ( sq_nm_g' * squaredNorm @v g'' - ( g' ^.^ g'' ) ^ ( 2 :: Int ) )
  / ( sq_nm_g' ^ ( 3 :: Int ) )
  where
    g', g'' :: v
    !g'  = bezier'  @v bez t
    !g'' = bezier'' @v bez t
    sq_nm_g' :: r
    !sq_nm_g' = squaredNorm @v g'
{-# INLINEABLE squaredCurvature #-}

-- | Signed curvature of a planar cubic Bézier curve.
signedCurvature :: Bezier ( ℝ 2 ) -> Double -> Double
signedCurvature bez t = ( g' × g'' ) / norm g' ^ ( 3 :: Int )
  where
    g', g'' :: T ( ℝ 2 )
    !g'  = bezier'  @( T ( ℝ 2 ) ) bez t
    !g'' = bezier'' @( T ( ℝ 2 ) ) bez t

-- | Subdivide a cubic Bézier curve into two parts.
subdivide :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> ( Bezier p, Bezier p )
subdivide ( Bezier {..} ) t = ( Bezier p0 q1 q2 pt, Bezier pt r1 r2 p3 )
  where
    pt, s, q1, q2, r1, r2 :: p
    !q1 = lerp @v t p0 p1
    !s  = lerp @v t p1 p2
    !r2 = lerp @v t p2 p3
    !q2 = lerp @v t q1 s
    !r1 = lerp @v t s  r2
    !pt = lerp @v t q2 r1
{-# INLINEABLE subdivide #-}

-- | Restrict a cubic Bézier curve to a sub-interval, re-parametrising
-- to \( [0,1] \).
restrict :: forall v r p. ( Torsor v p, Ring.Field r, Module r v ) => Bezier p -> ( r , r ) -> Bezier p
restrict bez ( a, b ) = fst $ ( flip ( subdivide @v ) b' ) $ snd $ subdivide @v bez a
  where
    b' = ( b Ring.- a ) Ring./ ( Ring.fromInteger 1 Ring.- a )
  -- TODO: this could be made more efficient.
  -- See e.g. "https://math.stackexchange.com/questions/4172835/cubic-b%C3%A9zier-spline-multiple-split"
  -- or the paper "On the numerical condition of Bernstein-Bézier subdivision process".
{-# INLINEABLE restrict #-}

-- | Polynomial coefficients of the derivative of the distance to a cubic Bézier curve.
ddist :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> p -> [ r ]
ddist ( Bezier {..} ) c = [ a5, a4, a3, a2, a1, a0 ]
  where
    v, v', v'', v''' :: v
    !v    = c  --> p0
    !v'   = p0 --> p1
    !v''  = p1 --> p0 ^+^ p1 --> p2
    !v''' = p0 --> p3 ^+^ 3 *^ ( p2 --> p1 )

    a0, a1, a2, a3, a4, a5 :: r
    !a0 = v ^.^ v'
    !a1 = 3 * squaredNorm v' + 2 * v ^.^ v''
    !a2 = 9 * v' ^.^ v'' + v ^.^ v'''
    !a3 = 6 * squaredNorm v'' + 4 * v' ^.^ v'''
    !a4 = 5 * v'' ^.^ v'''
    !a5 = squaredNorm v'''
{-# INLINEABLE ddist #-}

-- | Finds the closest point to a given point on a cubic Bézier curve.
closestPoint
  :: forall v r p. ( Torsor v p, Inner r v, RealFloat r, Prim r, NFData r )
  => Bezier p -> p -> ArgMin r ( r, p )
closestPoint pts c = pickClosest ( 0 :| 1 : distRoots )
  where
    distRoots :: [ r ]
    distRoots = filter ( \ r -> r > 0 && r < 1 ) ( realRoots 50 $ ddist @v pts c )

    pickClosest :: NonEmpty r -> ArgMin r ( r, p )
    pickClosest ( s :| ss ) = go s q nm0 ss
      where
        q :: p
        !q  = bezier @v pts s
        nm0 :: r
        !nm0 = squaredNorm ( c --> q :: v )
        go t p nm []  = Min ( Arg nm ( t, p ) )
        go t p nm ( t' : ts )
          | nm' < nm  = go t' p' nm' ts
          | otherwise = go t  p  nm  ts
          where
            p' :: p
            !p' = bezier @v pts t'
            nm' :: r
            !nm' = squaredNorm ( c --> p' :: v )
{-# INLINEABLE closestPoint #-}

-- | Drag a cubic Bézier curve to pass through a given point.
--
-- Given a cubic Bézier curve, a time \( 0 < t < 1 \) and a point `q`,
-- modifies the control points to make the curve pass through `q` at time `t`.
--
-- Affects the two control points depending on how far along the dragged point is.
-- For instance, dragging near the middle moves both control points equally,
-- while dragging near an endpoint will mostly affect the control point associated with that endpoint.
drag :: forall v r p. ( Torsor v p, Module r v, Fractional r ) => Bezier p -> r -> p -> Bezier p
drag ( Bezier {..} ) t q = Bezier { p0, p1 = p1', p2 = p2', p3 }
  where
    v0, v1, v2, v3, delta :: v
    !v0 = q --> p0
    !v1 = q --> p1
    !v2 = q --> p2
    !v3 = q --> p3
    !delta = ( recip $ t * ( -3 + t * ( 9 + t * ( -12 + 6 * t ) ) ) )
          *^ bezier @v ( Bezier v0 v1 v2 v3 ) t
    p1', p2' :: p
    !p1' = ( ( 1 - t ) *^ delta ) • p1
    !p2' = (     t     *^ delta ) • p2
{-# INLINEABLE drag #-}

-- | Compute parameter values for the self-intersection of a planar cubic Bézier curve, if such exist.
--
-- The parameter values might lie outside the interval [0,1],
-- indicating a self-intersection of the extended curve.
--
-- Formula taken from:
--  "A Basis for the Implicit Representation of Planar Rational Cubic Bézier Curves"
--      – Oliver J. D. Barrowclough, 2016
selfIntersectionParameters :: Bezier ( ℝ 2 ) -> [ Double ]
selfIntersectionParameters ( Bezier {..} ) = solveQuadratic c0 c1 c2
  where
    areaConstant :: ℝ 2 -> ℝ 2 -> ℝ 2 -> Double
    areaConstant ( ℝ2 x1 y1 ) ( ℝ2 x2 y2 ) ( ℝ2 x3 y3 ) =
      x1 * ( y2 - y3 ) + x2 * ( y3 - y1 ) + x3 * ( y1 - y2 )
    l0, l1, l2, l3, f1, f2, f3, c0, c1, c2 :: Double
    !l0 = areaConstant p3 p2 p1
    !l1 = areaConstant p2 p3 p0
    !l2 = areaConstant p1 p0 p3
    !l3 = areaConstant p0 p1 p2
    !f1 = 3 * ( l1 * l1 - 3 * l0 * l2 )
    !f2 = 3 * ( l2 * l2 - 3 * l1 * l3 )
    !f3 = 3 * ( 9 * l0 * l3 - l1 * l2 )
    !c0 = f2
    !c1 = f3 - 2 * f2
    !c2 = f1 + f2 - f3

-- | Critical points of the speed @|B′(t)|@ in @(0,1)@ for a cubic Bézier.
--
-- These are the roots of @d\/dt |B′(t)|²\/18 = 0@, a cubic polynomial whose
-- coefficients come from expanding @H(t)·H′(t) = 0@ where
-- @H(t) = d₀ + 2t·f + t²·e@ is the (scaled) hodograph,
-- @f = d₁ − d₀@ and @e = d₀ − 2d₁ + d₂@ its first and second differences.
speedCriticalPoints
  :: forall diff ptData
  .  ( Torsor diff ptData, Module Double diff, Inner Double diff )
  => Bezier ptData -> [ Complex Double ]
speedCriticalPoints ( Bezier {..} ) =
  let
    d0, d1, d2 :: diff
    !d0 = p0 --> p1
    !d1 = p1 --> p2
    !d2 = p2 --> p3
    !f  = d1 ^-^ d0
    !e  = ( d0 ^-^ ( 2 *^ d1 ) ) ^+^ d2
    !c3 = e ^.^ e
    !c2 = 3 * ( f ^.^ e )
    !c1 = 2 * ( f ^.^ f ) + ( d0 ^.^ e )
    !c0 = d0 ^.^ f
  in  nubBy ( (==) `on` realPart )
    $ sortOn realPart
    $ filter ( \ t -> realPart t > 0 && realPart t < 1 )
    $ roots epsilon 64 [c3, c2, c1, c0]
{-# INLINEABLE speedCriticalPoints #-}

-- | Extremal values of the Bézier parameter for a cubic Bézier curve.
extrema :: RealFloatNaN r => Bezier r -> [ r ]
extrema ( Bezier {..} ) = solveQuadratic c b a
  where
    !a = p3 - 3 * p2 + 3 * p1 - p0
    !b = 2 * ( p0 - 2 * p1 + p2 )
    !c = p1 - p0
{-# INLINEABLE extrema #-}
{-# SPECIALISE extrema @Float  #-}
{-# SPECIALISE extrema @Double #-}