Quadratic.hs

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

module Math.Bezier.Quadratic
  ( Bezier(..)
  , bezier, bezier', bezier''
  , curvature, squaredCurvature, signedCurvature
  , subdivide, restrict
  , ddist, closestPoint
  , interpolate
  , extrema
  , speedCriticalPoints
  )
  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 )

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

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

-- | Points defining a quadratic Bézier curve (Bernstein form).
--
-- @ p0 @ and @ p2 @ are endpoints, whereas @ p1 @ is a control point.
data Bezier p
  = Bezier
    { p0, p1, p2 :: !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) =
    show p1 ++ "--" ++ show p2 ++ "->" ++ show p3

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

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

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

-- | Curvature of a quadratic 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 quadratic 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
    sq_nm_g' :: r
    !sq_nm_g' = squaredNorm @v g'
{-# INLINEABLE squaredCurvature #-}

-- | Signed curvature of a planar quadratic 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

-- | Subdivide a quadratic 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 pt, Bezier pt r1 p2 )
  where
    pt, q1, r1 :: p
    !q1 = lerp @v t p0 p1
    !r1 = lerp @v t p1 p2
    !pt = lerp @v t q1 r1
{-# INLINEABLE subdivide #-}

-- | Restrict a quadratic 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 quadratic Bézier curve.
ddist :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> p -> [ r ]
ddist ( Bezier {..} ) c = [ a3, a2, a1, a0 ]
  where
    v, v', v'' :: v
    !v   = c  --> p0
    !v'  = p0 --> p1
    !v'' = p1 --> p0 ^+^ p1 --> p2

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

-- | Finds the closest point to a given point on a quadratic 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 )
{-# INLINEABLE closestPoint #-}

-- | Interpolation of a quadratic Bézier control point, given path points and Bézier parameter.
--
-- That is, given points `p0`, `p2`, and `q`, and parameter \( 0 < t < 1 \),
-- this function finds the unique control point `p1` such that
-- the quadratic Bézier curve with parameters `p0`, `p1`, `p2` passes through
-- the point `q` at time `t`.
interpolate :: forall v r p. ( Torsor v p, Module r v, Fractional r ) => p -> p -> r -> p -> Bezier p
interpolate p0 p2 t q = Bezier {..}
  where
    p1 :: p
    !p1 = (   ( 0.5 * ( t - 1 ) / t ) *^ ( q --> p0 :: v )
          ^+^ ( 0.5 * t / ( t - 1 ) ) *^ ( q --> p2 :: v )
          ) • q
{-# INLINEABLE interpolate #-}

-- | Critical points of the speed @|B′(t)|@ in @(0,1)@ for a quadratic Bézier.
--
-- These are the roots of @d\/dt |B′(t)|² = 0@, i.e. where the hodograph
-- @H(t) = d₀ + t·f@ is perpendicular to its (constant) derivative @f = d₁ − d₀@.
speedCriticalPoints
  :: forall diff ptData
  .  ( Torsor diff ptData, Module Double diff, Inner Double diff )
  => Bezier ptData -> [ Double ]
speedCriticalPoints ( Bezier {..} ) =
    let
      d0, d1 :: diff
      !d0 = p0 --> p1
      !d1 = p1 --> p2
      !f  = d1 ^-^ d0
      !ff = f ^.^ f
      !tc = -(f ^.^ d0) / ff
    in [ tc | ff > 0, tc > 0, tc < 1 ]
{-# INLINEABLE speedCriticalPoints #-}

-- | Extremal values of the Bézier parameter for a quadratic Bézier curve.
extrema :: Fractional r => Bezier r -> [ r ]
extrema ( Bezier {..} ) = [ t ]
  where
    !t = ( p0 - p1 ) / ( p0 - 2 * p1 + p2 )
{-# INLINEABLE extrema #-}