content: module Graphics.CubicBezier (
        CubicBezier (..)
      , firstDerivativeAt
      , secondDerivativeAt
      , splitAt
      , isClockwise
      , inflectionPoints
      , maxArcLength
  ) where
  
  import Prelude hiding (splitAt)
  
  import Linear
  import Control.Lens
  import Data.Complex
  
  import Graphics.Curve
  
  data CubicBezier = CubicBezier { _p1 :: V2 Double
                                 , _c1 :: V2 Double
                                 , _c2 :: V2 Double
                                 , _p2 :: V2 Double
                                 } deriving Show
  
  instance Curve CubicBezier where
      pointAt bezier t =  (1 - t) ** 3 *^ _p1 bezier +
                          (1 - t) ** 2 * 3 * t *^ _c1 bezier +
                          t ** 2 * (1 - t) * 3 *^ _c2 bezier +
                          t ** 3 *^ _p2 bezier
  
  firstDerivativeAt :: CubicBezier -> Double -> V2 Double
  firstDerivativeAt bezier t = (1 - t) ** 2 * 3 *^ (_c1 bezier - _p1 bezier) +
                               (1 - t) * t * 6 *^ (_c2 bezier - _c1 bezier) +
                               t * t * 3 *^ (_p2 bezier - _c2 bezier)
  
  secondDerivativeAt :: CubicBezier -> Double -> V2 Double
  secondDerivativeAt bezier t = (1 - t) * 6 *^ (_c2 bezier - 2 *^ _c1 bezier + _p1 bezier) +
                                t * 6 *^ (_p2 bezier - 2 *^ _c2 bezier + _c1 bezier)
  
  
  splitAt :: CubicBezier -> Double -> (CubicBezier, CubicBezier)
  splitAt bezier t = (CubicBezier (_p1 bezier) p0 p01 dp, CubicBezier dp p12 p2 (_p2 bezier))
      where
          p0 = _p1 bezier + t *^ (_c1 bezier - _p1 bezier)
          p1 = _c1 bezier + t *^ (_c2 bezier - _c1 bezier)
          p2 = _c2 bezier + t *^ (_p2 bezier - _c2 bezier)
  
          p01 = p0 + t *^ (p1 - p0)
          p12 = p1 + t *^ (p2 - p1)
  
          dp = p01 + t *^ (p12 - p01)
  
  isClockwise :: CubicBezier -> Bool
  isClockwise bezier = isClockwise4 (_p1 bezier) (_c1 bezier) (_c2 bezier) (_p2 bezier)
  
  inflectionPoints :: CubicBezier -> [Double]
  inflectionPoints bezier
      | a /= 0 = realInflectionPoints [t1, t2]
      | otherwise = realInflectionPoints [t]
      where
          pa = _c1 bezier - _p1 bezier
          pb = _c2 bezier - _c1 bezier - pa
          pc = _p2 bezier - _c2 bezier - pa - 2 *^ pb
  
          a = (pb ^. _x * pc ^. _y - pb ^. _y * pc ^. _x) :+ 0
          b = (pa ^. _x * pc ^. _y - pa ^. _y * pc ^. _x) :+ 0
          c = (pa ^. _x * pb ^. _y - pa ^. _y * pb ^. _x) :+ 0
  
          -- linear case
          t = -c / b
  
          -- quadratic case
          t1 = (-b + sqrt (b * b  - 4 * a * c)) / (2 * a)
          t2 = (-b - sqrt (b * b  - 4 * a * c)) / (2 * a)
  
  realInflectionPoints :: [Complex Double] -> [Double]
  realInflectionPoints = map realPart . filter isInflectionPoint
  
  isInflectionPoint :: Complex Double -> Bool
  isInflectionPoint c = imagPart c == 0 && realPart c > 0 && realPart c < 1
  
  maxArcLength :: CubicBezier -> Double
  maxArcLength bezier =
      distance (_p1 bezier) (_c1 bezier) +
      distance (_c1 bezier) (_c2 bezier) +
      distance (_c2 bezier) (_p2 bezier)
  