content: {-# LANGUAGE FlexibleInstances #-}
  {-# LANGUAGE GADTs #-}
  {-# LANGUAGE BangPatterns #-}
  {-# OPTIONS_GHC -fno-warn-orphans #-}
  module Graphics.Rasterific.CubicBezier
      ( cubicBezierCircle
      , cubicBezierFromPath
      , cubicBezierBreakAt
      , divideCubicBezier
      , clipCubicBezier
      , decomposeCubicBeziers
      , sanitizeCubicBezier
      , sanitizeCubicBezierFilling
      , offsetCubicBezier
      , flattenCubicBezier
      , cubicBezierLengthApproximation
      , cubicBezierBounds
      , cubicFromQuadraticBezier
      , isCubicBezierPoint
      ) where
  
  import Prelude hiding( or )
  
  import Control.Applicative( liftA2 )
  import Graphics.Rasterific.Linear
               ( V2( .. )
               , (^-^)
               , (^+^)
               , (^*)
               , norm
               , lerp
               )
  import Data.List( nub )
  import Graphics.Rasterific.Operators
  import Graphics.Rasterific.Types
  import Graphics.Rasterific.QuadraticFormula
  
  -- | Create a list of cubic bezier patch from a list of points.
  --
  -- > cubicBezierFromPath [a, b, c, d, e] = [CubicBezier a b c d]
  -- > cubicBezierFromPath [a, b, c, d, e, f, g] =
  -- >    [CubicBezier a b c d, CubicBezier d e f g]
  --
  cubicBezierFromPath :: [Point] -> [CubicBezier]
  cubicBezierFromPath (a:b:c:rest@(d:_)) =
      CubicBezier a b c d : cubicBezierFromPath rest
  cubicBezierFromPath _ = []
  
  cubicBezierLengthApproximation :: CubicBezier -> Float
  cubicBezierLengthApproximation (CubicBezier a _ _ d) =
      norm $ d ^-^ a
  
  -- | Represent a circle of radius 1 centered on 0 of
  -- a cubic bezier curve.
  cubicBezierCircle :: [CubicBezier]
  cubicBezierCircle =
      [ CubicBezier (V2 0 1) (V2 c 1) (V2 1 c) (V2 1 0)
      , CubicBezier (V2 1 0) (V2 1 (-c)) (V2 c (-1)) (V2 0 (-1))
      , CubicBezier (V2 0 (-1)) (V2 (-c) (-1)) (V2 (-1) (-c)) (V2 (-1) 0)
      , CubicBezier (V2 (-1) 0) (V2 (-1) c) (V2 (-c) 1) (V2 0 1)
      ]
    where c = 0.551915024494 -- magic constant? magic constant.
  
  straightLine :: Point -> Point -> CubicBezier
  straightLine a b = CubicBezier a p p b
    where p = a `midPoint` b
  
  isSufficientlyFlat :: Float -- ^ Tolerance
                     -> CubicBezier
                     -> Bool
  isSufficientlyFlat tol (CubicBezier a b c d) =
      x + y <= tolerance
    where u = (b ^* 3) ^-^ (a ^* 2) ^-^ d
          v = (c ^* 3) ^-^ (d ^* 2) ^-^ a
          (^*^) = liftA2 (*)
          V2 x y = vmax (u ^*^ u) (v ^*^ v)
          tolerance = 16 * tol * tol
  
  splitCubicBezier :: CubicBezier -> (Point, Point, Point, Point, Point, Point)
  {-# INLINE splitCubicBezier #-}
  splitCubicBezier (CubicBezier a b c d) = (ab, bc, cd, abbc, bccd, abbcbccd)
    where
      --                     BC
      --         B X----------X---------X C
      --    ^     /      ___/   \___     \     ^
      --   u \   /   __X------X------X_   \   / v
      --      \ /___/ ABBC       BCCD  \___\ /
      --    AB X/                          \X CD
      --      /                              \
      --     /                                \
      --    /                                  \
      -- A X                                    X D
      ab = a `midPoint` b
      bc = b `midPoint` c
      cd = c `midPoint` d
  
      abbc = ab `midPoint` bc
      bccd = bc `midPoint` cd
      abbcbccd = abbc `midPoint` bccd
  
  flattenCubicBezier :: CubicBezier -> Container Primitive
  flattenCubicBezier bezier@(CubicBezier a _ _ d)
      | isSufficientlyFlat 1 bezier = pure $ CubicBezierPrim bezier
      | otherwise =
          flattenCubicBezier (CubicBezier a ab abbc abbcbccd) <>
              flattenCubicBezier (CubicBezier abbcbccd bccd cd d)
    where
      (ab, _bc, cd, abbc, bccd, abbcbccd) = splitCubicBezier bezier
  
  --               3                    2            2                  3
  -- x(t) = (1 - t) ∙x     + 3∙t∙(1 - t) ∙x     + 3∙t ∙(1 - t)∙x     + t ∙x
  --                   0                    1                    2          3
  --
  --               3                    2            2                  3
  -- y(t) = (1 - t) ∙y     + 3∙t∙(1 - t) ∙y     + 3∙t ∙(1 - t)∙y     + t ∙y
  --                   0                    1                    2          3
  
  -- Other representation:
  --                3                2        2              3
  -- B(t) = x(1 - t)  + 3∙y∙t∙(1 - t)  + 3∙z∙t ∙(1 - t) + w∙t
  
  
  -- | Represent the cubic bezier curve as a vector ready
  -- for matrix multiplication
  data CachedBezier = CachedBezier
      { _cachedA :: {-# UNPACK #-} !Float
      , _cachedB :: {-# UNPACK #-} !Float
      , _cachedC :: {-# UNPACK #-} !Float
      , _cachedD :: {-# UNPACK #-} !Float
      }
  
  cacheBezier :: CubicBezier -> (CachedBezier, CachedBezier)
  cacheBezier (CubicBezier p0@(V2 x0 y0) p1 p2 p3) =
      (CachedBezier x0 bX cX dX, CachedBezier y0 bY cY dY)
    where
     V2 bX bY = p1 ^* 3 ^-^ p0 ^* 3
     V2 cX cY = p2 ^* 3 ^-^ p1 ^* 6 + p0 ^* 3
     V2 dX dY = p3 ^-^ p2 ^* 3 ^+^ p1 ^* 3 ^-^ p0
  
  cachedBezierAt :: CachedBezier -> Float -> Float
  cachedBezierAt (CachedBezier a b c d) t =
      a + b * t + c * tSquare + tCube * d
    where
      tSquare = t * t
      tCube = tSquare * t
  
  cachedBezierDerivative :: CachedBezier -> QuadraticFormula Float
  cachedBezierDerivative (CachedBezier _ b c d) =
      QuadraticFormula (3 * d) (2 * c) b
  
  -- | Find the coefficient of the extremum points
  extremums :: CachedBezier -> [Float]
  extremums cached =
    [ root | root <- formulaRoots $ cachedBezierDerivative cached
           , 0 <= root && root <= 1.0 ]
  
  extremumPoints :: (CachedBezier, CachedBezier) -> [Point]
  extremumPoints (onX, onY) = toPoints <$> nub (extremums onX <> extremums onY)
    where toPoints at = V2 (cachedBezierAt onX at) (cachedBezierAt onY at)
  
  cubicBezierBounds :: CubicBezier -> [Point]
  cubicBezierBounds bez@(CubicBezier p0 _ _ p3) =
      p0 : p3 : extremumPoints (cacheBezier bez)
  
  offsetCubicBezier :: Float -> CubicBezier -> Container Primitive
  offsetCubicBezier offset bezier@(CubicBezier a b c d)
      | isSufficientlyFlat 1 bezier =
          pure . CubicBezierPrim $ CubicBezier shiftedA shiftedB shiftedC shiftedD
      | otherwise =
          recurse (CubicBezier a ab abbc abbcbccd) <>
              recurse (CubicBezier abbcbccd bccd cd d)
    where
      recurse = offsetCubicBezier offset
  
      u = a `normal` b
      v = c `normal` d
  
      --                     BC
      --         B X----------X---------X C
      --    ^     /      ___/   \___     \     ^
      --   u \   /   __X------X------X_   \   / v
      --      \ /___/ ABBC       BCCD  \___\ /
      --    AB X/                          \X CD
      --      /                              \
      --     /                                \
      --    /                                  \
      -- A X                                    X D
      (ab, bc, cd, abbc, bccd, abbcbccd) = splitCubicBezier bezier
  
      w = ab `normal` bc
      x = bc `normal` cd
  
      shiftedA = a ^+^ (u ^* offset)
      shiftedD = d ^+^ (v ^* offset)
  
      {-shiftedABBCBCCD = abbcbccd ^+^ (w ^* offset)-}
      shiftedB = b ^+^ (w ^* offset)
      shiftedC = c ^+^ (x ^* offset)
  
  -- | Clamp the cubic bezier curve inside a rectangle
  -- given in parameter.
  clipCubicBezier
      :: Point   -- ^ Point representing the "minimal" point for cliping
      -> Point  -- ^ Point representing the "maximal" point for cliping
      -> CubicBezier -- ^ The cubic bezier curve to be clamped
      -> Container Primitive
  clipCubicBezier mini maxi bezier@(CubicBezier a b c d)
      -- If we are in the range bound, return the curve
      -- unaltered
      | insideX && insideY = pure $ CubicBezierPrim bezier
      -- If one of the component is outside, clamp
      -- the components on the boundaries and output a
      -- straight line on this boundary. Useful for the
      -- filing case, to clamp the polygon drawing on
      -- the edge
      | outsideX || outsideY =
          pure . CubicBezierPrim $ clampedA `straightLine` clampedD
      -- Not completly inside nor outside, just divide
      -- and conquer.
      | otherwise =
          recurse (CubicBezier a ab abbc m) <>
              recurse (CubicBezier m bccd cd d)
    where -- Minimal & maximal dimension of the bezier curve
          bmin = vmin a . vmin b $ vmin c d
          bmax = vmax a . vmax b $ vmax c d
  
          recurse = clipCubicBezier mini maxi
  
          clamper = clampPoint mini maxi
          clampedA = clamper a
          clampedD = clamper d
  
          V2 insideX insideY = mini ^<=^ bmin ^&&^ bmax ^<=^ maxi
          V2 outsideX outsideY = bmax ^<=^ mini ^||^ maxi ^<=^ bmin
  
          --                     BC
          --         B X----------X---------X C
          --          /      ___/   \___     \
          --         /   __X------X------X_   \
          --        /___/ ABBC       BCCD  \___\
          --    AB X/                          \X CD
          --      /                              \
          --     /                                \
          --    /                                  \
          -- A X                                    X D
          (ab, _bc, cd, abbc, bccd, abbcbccd) = splitCubicBezier bezier
  
          edgeSeparator = vabs (abbcbccd ^-^ mini) ^<^ vabs (abbcbccd ^-^ maxi)
          edge = vpartition edgeSeparator mini maxi
          m = vpartition (vabs (abbcbccd ^-^ edge) ^< 0.1) edge abbcbccd
  
  divideCubicBezier :: CubicBezier -> (CubicBezier, CubicBezier)
  divideCubicBezier bezier@(CubicBezier a _ _ d) = (left, right) where
    left = CubicBezier a ab abbc abbcbccd
    right = CubicBezier abbcbccd bccd cd d
    (ab, _bc, cd, abbc, bccd, abbcbccd) = splitCubicBezier bezier
  
  -- | Will subdivide the bezier from 0 to coeff and coeff to 1
  cubicBezierBreakAt :: CubicBezier -> Float
                     -> (CubicBezier, CubicBezier)
  cubicBezierBreakAt (CubicBezier a b c d) val =
      (CubicBezier a ab abbc abbcbccd, CubicBezier abbcbccd bccd cd d)
    where
      ab = lerp val b a
      bc = lerp val c b
      cd = lerp val d c
  
      abbc = lerp val bc ab
      bccd = lerp val cd bc
      abbcbccd = lerp val bccd abbc
  
  decomposeCubicBeziers :: CubicBezier -> Producer EdgeSample
  decomposeCubicBeziers cb@(CubicBezier a b c d)
     -- handle case of self closed bezier curve
    | not (a `isDistingableFrom` d) && ((a `isDistingableFrom` b) || (a `isDistingableFrom` c)) =
      let (l, r) = cubicBezierBreakAt cb 0.5 in
      decomposeCubicBeziers l . decomposeCubicBeziers r 
  decomposeCubicBeziers (CubicBezier (V2 aRx aRy) (V2 bRx bRy) (V2 cRx cRy) (V2 dRx dRy)) =
      go aRx aRy bRx bRy cRx cRy dRx dRy where
    go ax ay _bx _by _cx _cy dx dy cont | insideX && insideY =
      let !px = fromIntegral $ min floorAx floorDx
          !py = fromIntegral $ min floorAy floorDy
          !w = px + 1 - (dx `middle` ax)
          !h = dy - ay
      in
      EdgeSample (px + 0.5) (py + 0.5) (w * h) h : cont
      where
        floorAx, floorAy :: Int
        !floorAx = floor ax
        !floorAy = floor ay
  
        !floorDx = floor dx
        !floorDy = floor dy
  
        !insideX =
            floorAx == floorDx || ceiling ax == (ceiling dx :: Int)
        !insideY =
            floorAy == floorDy || ceiling ay == (ceiling dy :: Int)
  
  
    go !ax !ay !bx !by !cx !cy !dx !dy cont =
       go ax ay abx aby abbcx abbcy mx my $
          go mx my bccdx bccdy cdx cdy dx dy cont
      where
        --                     BC
        --         B X----------X---------X C
        --          /      ___/   \___     \
        --         /   __X------X------X_   \
        --        /___/ ABBC       BCCD  \___\
        --    AB X/                          \X CD
        --      /                              \
        --     /                                \
        --    /                                  \
        -- A X                                    X D
        !abx = ax `middle` bx
        !aby = ay `middle` by
        !bcx = bx `middle` cx
        !bcy = by `middle` cy
        !cdx = cx `middle` dx
        !cdy = cy `middle` dy
        !abbcx = abx `middle` bcx
        !abbcy = aby `middle` bcy
        !bccdx = bcx `middle` cdx
        !bccdy = bcy `middle` cdy
  
        !abbcbccdx = abbcx `middle` bccdx
        !abbcbccdy = abbcy `middle` bccdy
  
        !mx | abs (abbcbccdx - mini) < 0.1 = mini
            | abs (abbcbccdx - maxi) < 0.1 = maxi
            | otherwise = abbcbccdx
              where !mini = fromIntegral (floor abbcbccdx :: Int)
                    !maxi = fromIntegral (ceiling abbcbccdx :: Int)
  
        !my | abs (abbcbccdy - mini) < 0.1 = mini
            | abs (abbcbccdy - maxi) < 0.1 = maxi
            | otherwise = abbcbccdy
              where !mini = fromIntegral (floor abbcbccdy :: Int)
                    !maxi = fromIntegral (ceiling abbcbccdy :: Int)
  
  isCubicBezierPoint :: CubicBezier -> Bool
  isCubicBezierPoint (CubicBezier a b c d) =
    not $ a `isDistingableFrom` b || 
          b `isDistingableFrom` c ||
          c `isDistingableFrom` d
  
  sanitizeCubicBezier :: CubicBezier -> Container Primitive
  sanitizeCubicBezier bezier@(CubicBezier a b c d)
    | a `isDistingableFrom` b &&
      c `isDistingableFrom` d =
         pure . CubicBezierPrim $ bezier
    | ac `isDistingableFrom` b &&
       bd `isDistingableFrom` c =
        pure . CubicBezierPrim $ bezier
    | ac `isDistingableFrom` b =
        pure . CubicBezierPrim $ CubicBezier a ac c d
    | bd `isDistingableFrom` c =
        pure . CubicBezierPrim $ CubicBezier a b bd d
    | otherwise = mempty
      where ac = a `midPoint` c
            bd = a `midPoint` d
  
  sanitizeCubicBezierFilling :: CubicBezier -> Container Primitive
  sanitizeCubicBezierFilling bezier@(CubicBezier a b c d)
    | isDegenerate a || isDegenerate b ||
      isDegenerate c || isDegenerate d = mempty
    | otherwise = pure $ CubicBezierPrim bezier
  
  cubicFromQuadraticBezier :: Bezier -> CubicBezier
  cubicFromQuadraticBezier (Bezier p0 p1 p2) = CubicBezier p0 pa pb p2 where
    pa = p0 ^+^ (p1 ^-^ p0) ^* (2 / 3)
    pb = p2 ^+^ (p1 ^-^ p2) ^* (2 / 3)
  
  