{-# LANGUAGE AllowAmbiguousTypes   #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE DeriveTraversable     #-}
{-# LANGUAGE DerivingStrategies    #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RecordWildCards       #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeApplications      #-}
{-# LANGUAGE UndecidableInstances  #-}

module Math.Bezier.Quadratic
  ( Bezier(..)
  , bezier, bezier'
  , subdivide
  , closestPoint
  )
  where

-- base
import Data.List.NonEmpty
  ( NonEmpty(..) )
import GHC.Generics
  ( Generic )

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

-- MetaBrush
import Math.Module
  ( Module (..)
  , lerp
  , Inner(..), squaredNorm
  )
import Math.RealRoots
  ( realRoots )

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

-- | Points defining a quadratic Bézier curve.
--
-- @ p0 @ and @ p2 @ are endpoints, whereas @ p1 @ is a control point.
data Bezier p
  = Bezier
    { p0 :: !p
    , p1 :: !p
    , p2 :: !p
    }
  deriving stock ( Show, Generic, Functor, Foldable, Traversable )

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

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

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

-- | 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 ) => Bezier p -> p -> ( r, p )
closestPoint pts@( Bezier { .. } ) c = pickClosest ( 0 :| 1 : roots )
  where
    roots :: [ r ]
    roots = filter ( \ r -> r > 0 && r < 1 ) ( realRoots [ a0, a1, a2, a3 ] )

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

    pickClosest :: NonEmpty 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 _  []  = ( 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 )