metabrush/src/splines/Math/Bezier/Cubic.hs

{-# LANGUAGE AllowAmbiguousTypes   #-}
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE DeriveTraversable     #-}
{-# LANGUAGE DerivingVia           #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE NamedFieldPuns        #-}
{-# LANGUAGE NegativeLiterals      #-}
{-# LANGUAGE RecordWildCards       #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE StandaloneDeriving    #-}
{-# LANGUAGE TypeApplications      #-}
{-# LANGUAGE UndecidableInstances  #-}

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

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

-- groups
import Data.Group
  ( Group )

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

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

-- MetaBrush
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 )
import Math.Vector2D
  ( Point2D(..), Vector2D(..) )

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

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

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

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

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

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

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

-- | 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
  = ( 6 *^ )
  $ lerp @v t
      ( p1 --> p0 ^+^ p1 --> p2 )
      ( p2 --> p1 ^+^ p2 --> p3 )

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

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

-- | Signed curvature of a planar cubic Bézier curve.
signedCurvature :: forall r. Floating r => Bezier ( Point2D r ) -> r -> r
signedCurvature bez t = ( g' `cross` g'' ) / norm g'
  where
    g', g'' :: Vector2D r
    g'  = bezier'  @( Vector2D r ) bez t
    g'' = bezier'' @( Vector2D r ) 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

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

-- | 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 : roots )
  where
    roots :: [ r ]
    roots = 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 )

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

-- | 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 :: forall r. RealFloat r => Bezier ( Point2D r ) -> [ r ]
selfIntersectionParameters ( Bezier {..} ) = solveQuadratic c0 c1 c2
  where
    areaConstant :: Point2D r -> Point2D r -> Point2D r -> r
    areaConstant ( Point2D x1 y1 ) ( Point2D x2 y2 ) ( Point2D x3 y3 ) =
      x1 * ( y2 - y3 ) + x2 * ( y3 - y1 ) + x3 * ( y1 - y2 )
    l0, l1, l2, l3, f1, f2, f3, c0, c1, c2 :: r
    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