    {-# LANGUAGE FlexibleContexts #-}
    {-# LANGUAGE TypeFamilies     #-}

    -----------------------------------------------------------------------------
    -- |
    -- Module      :  Diagrams.CubicSpline.Boehm
    -- Copyright   :  (c) 2015 diagrams-lib team (see LICENSE)
    -- License     :  BSD-style (see LICENSE)
    -- Maintainer  :  diagrams-discuss@googlegroups.com
    --
    -- Boehm's algorithm for converting a cubic B-spline into a sequence
    -- of cubic Bezier curves.
    --
    -- See
    --
    --   * Thomas W. Sederberg, /An Introduction to B-Spline Curves/,
    --     <http://web.archive.org/web/20120227050519/http://tom.cs.byu.edu/~455/bs.pdf>
    --
    --   * Lyle Ramshaw, /Blossoming: A Connect-the-Dots Approach to Splines/,
    --     <http://www.hpl.hp.com/techreports/Compaq-DEC/SRC-RR-19.pdf>
    --
    -----------------------------------------------------------------------------
    module Diagrams.CubicSpline.Boehm
           ( BSpline
           , bsplineToBeziers
           , bspline
           ) where

    import           Data.List          (sort, tails)
    import           Diagrams.Core      (N, Point, V, origin)
    import           Diagrams.Located   (at, loc, unLoc)
    import           Diagrams.Segment   (FixedSegment (..), fromFixedSeg)
    import           Diagrams.TrailLike (TrailLike, fromLocSegments)
    import           Diagrams.Util      (iterateN)
    import           Linear.Vector      (Additive, lerp)

    type BSpline v n = [Point v n]

    -- | @affineCombo a b t x y@ computes an affine combination of x and y
    --   which lies at parameter t, if x has parameter a and y has parameter b.
    --   The usual @lerp@ arises by giving x parameter 0 and y parameter 1.
    affineCombo :: (Additive f, Fractional a) => a -> a -> a -> f a -> f a -> f a
    affineCombo :: a -> a -> a -> f a -> f a -> f a
    affineCombo a :: a
    a b :: a
    b t :: a
    t x :: f a
    x y :: f a
    y = a -> f a -> f a -> f a
    forall (f :: * -> *) a.
    (Additive f, Num a) =>
    a -> f a -> f a -> f a
    lerp ((a
    ta -> a -> a
    forall a. Num a => a -> a -> a
    -a
    a)a -> a -> a
    forall a. Fractional a => a -> a -> a
    /(a
    ba -> a -> a
    forall a. Num a => a -> a -> a
    -a
    a)) f a
    y f a
    x

    -- | @windows k xs@ yields all the length-@k@ windows from @xs@, e.g.
    --   @windows 3 [a,b,c,d,e] == [[a,b,c], [b,c,d], [c,d,e]]@.
    windows :: Int -> [a] -> [[a]]
    windows :: Int -> [a] -> [[a]]
    windows k :: Int
    k = ([a] -> Bool) -> [[a]] -> [[a]]
    forall a. (a -> Bool) -> [a] -> [a]
    takeWhile ((Int -> Int -> Bool
    forall a. Eq a => a -> a -> Bool
    ==Int
    k) (Int -> Bool) -> ([a] -> Int) -> [a] -> Bool
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [a] -> Int
    forall (t :: * -> *) a. Foldable t => t a -> Int
    length) ([[a]] -> [[a]]) -> ([a] -> [[a]]) -> [a] -> [[a]]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . ([a] -> [a]) -> [[a]] -> [[a]]
    forall a b. (a -> b) -> [a] -> [b]
    map (Int -> [a] -> [a]
    forall a. Int -> [a] -> [a]
    take Int
    k) ([[a]] -> [[a]]) -> ([a] -> [[a]]) -> [a] -> [[a]]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [a] -> [[a]]
    forall a. [a] -> [[a]]
    tails

    -- | @extend k xs@ extends @xs@ on both ends by prepending @k@ copies
    -- of its head and appending @k@ copies of its last element.  For example,
    -- @extend 2 [1..5] == [1,1,1,2,3,4,5,5,5]@.
    extend :: Int -> [a] -> [a]
    extend :: Int -> [a] -> [a]
    extend k :: Int
    k xs :: [a]
    xs = Int -> a -> [a]
    forall a. Int -> a -> [a]
    replicate Int
    k ([a] -> a
    forall a. [a] -> a
    head [a]
    xs) [a] -> [a] -> [a]
    forall a. [a] -> [a] -> [a]
    ++ [a]
    xs [a] -> [a] -> [a]
    forall a. [a] -> [a] -> [a]
    ++ Int -> a -> [a]
    forall a. Int -> a -> [a]
    replicate Int
    k ([a] -> a
    forall a. [a] -> a
    last [a]
    xs)

    -- | A "polar point" is a point along with three knot values.
    --   We consider the "blossom" of a cubic spline, a 3-ary symmetric
    --   polynomial; a polar point consists of 3 values paired with the
    --   output of the blossom at those input values.  Blossoms have nice
    --   affine properties so this makes it easy to keep track of how
    --   points may be combined to yield other points of interest.
    --
    --   Invariant: knot values are in nondecreasing order.
    data PolarPt v n = PP { PolarPt v n -> Point v n
    unPP :: Point v n, PolarPt v n -> [n]
    _knots :: [n] }

    mkPolarPt :: Ord n => Point v n -> [n] -> PolarPt v n
    mkPolarPt :: Point v n -> [n] -> PolarPt v n
    mkPolarPt pt :: Point v n
    pt kts :: [n]
    kts = Point v n -> [n] -> PolarPt v n
    forall (v :: * -> *) n. Point v n -> [n] -> PolarPt v n
    PP Point v n
    pt ([n] -> [n]
    forall a. Ord a => [a] -> [a]
    sort [n]
    kts)

    -- | Precondition: the knots of the two polar points overlap, like abc
    --   and bcd.  The @Int@ should be 0 or 1, indicating which knot to
    --   replicate (0 means to replicate b, yielding bbc, 1 means to
    --   replicate c, yielding bcc).
    combine
      :: (Additive v, Fractional n, Ord n)
      => Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
    combine :: Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
    combine k :: Int
    k (PP pt1 :: Point v n
    pt1 kts1 :: [n]
    kts1) (PP pt2 :: Point v n
    pt2 kts2 :: [n]
    kts2)
      = Point v n -> [n] -> PolarPt v n
    forall n (v :: * -> *). Ord n => Point v n -> [n] -> PolarPt v n
    mkPolarPt
          (n -> n -> n -> Point v n -> Point v n -> Point v n
    forall (f :: * -> *) a.
    (Additive f, Fractional a) =>
    a -> a -> a -> f a -> f a -> f a
    affineCombo ([n] -> n
    forall a. [a] -> a
    head [n]
    kts1) ([n] -> n
    forall a. [a] -> a
    last [n]
    kts2) n
    newKt Point v n
    pt1 Point v n
    pt2)
          (n
    newKt n -> [n] -> [n]
    forall a. a -> [a] -> [a]
    : Int -> [n] -> [n]
    forall a. Int -> [a] -> [a]
    drop 1 [n]
    kts1)
      where
        newKt :: n
    newKt = [n]
    kts2 [n] -> Int -> n
    forall a. [a] -> Int -> a
    !! Int
    k

    -- | Convert a uniform cubic B-spline to a sequence of cubic beziers.
    --   (/Uniform/ refers to the fact that the knots are assumed to be
    --   evenly spaced, with no duplicates.)  The knots at the end are
    --   replicated so the cubic spline begins and ends at the first and
    --   last control points, tangent to the line from the end control
    --   point to the next.
    bsplineToBeziers
      :: (Additive v, Fractional n, Num n, Ord n)
      => BSpline v n
      -> [FixedSegment v n]
    bsplineToBeziers :: BSpline v n -> [FixedSegment v n]
    bsplineToBeziers controls :: BSpline v n
    controls = [FixedSegment v n]
    beziers
      where
        n :: Int
    n                            = BSpline v n -> Int
    forall (t :: * -> *) a. Foldable t => t a -> Int
    length BSpline v n
    controls
        numKnots :: Int
    numKnots                     = Int
    n Int -> Int -> Int
    forall a. Num a => a -> a -> a
    + 2
        knots :: [n]
    knots                        = Int -> (n -> n) -> n -> [n]
    forall a. Int -> (a -> a) -> a -> [a]
    iterateN Int
    numKnots (n -> n -> n
    forall a. Num a => a -> a -> a
    +1n -> n -> n
    forall a. Fractional a => a -> a -> a
    /(Int -> n
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral Int
    numKnots n -> n -> n
    forall a. Num a => a -> a -> a
    - 1)) 0

        -- The control points are P(a,b,c), P(b,c,d), P(c,d,e), and so on.
        controls' :: [PolarPt v n]
    controls' = (Point v n -> [n] -> PolarPt v n)
    -> BSpline v n -> [[n]] -> [PolarPt v n]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith Point v n -> [n] -> PolarPt v n
    forall n (v :: * -> *). Ord n => Point v n -> [n] -> PolarPt v n
    mkPolarPt (Int -> BSpline v n -> BSpline v n
    forall a. Int -> [a] -> [a]
    extend 2 BSpline v n
    controls) (Int -> [n] -> [[n]]
    forall a. Int -> [a] -> [[a]]
    windows 3 ([n] -> [[n]]) -> [n] -> [[n]]
    forall a b. (a -> b) -> a -> b
    $ Int -> [n] -> [n]
    forall a. Int -> [a] -> [a]
    extend 2 [n]
    knots)

        -- The bezier internal control points are affine combinations of
        -- the spline control points.
        bezierControls :: [(PolarPt v n, PolarPt v n)]
    bezierControls        = ([PolarPt v n] -> (PolarPt v n, PolarPt v n))
    -> [[PolarPt v n]] -> [(PolarPt v n, PolarPt v n)]
    forall a b. (a -> b) -> [a] -> [b]
    map [PolarPt v n] -> (PolarPt v n, PolarPt v n)
    forall (v :: * -> *) n.
    (Additive v, Fractional n, Ord n) =>
    [PolarPt v n] -> (PolarPt v n, PolarPt v n)
    combineC (Int -> [PolarPt v n] -> [[PolarPt v n]]
    forall a. Int -> [a] -> [[a]]
    windows 2 [PolarPt v n]
    controls')
        combineC :: [PolarPt v n] -> (PolarPt v n, PolarPt v n)
    combineC [pabc :: PolarPt v n
    pabc, pbcd :: PolarPt v n
    pbcd] = (Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
    forall (v :: * -> *) n.
    (Additive v, Fractional n, Ord n) =>
    Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
    combine 0 PolarPt v n
    pabc PolarPt v n
    pbcd, Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
    forall (v :: * -> *) n.
    (Additive v, Fractional n, Ord n) =>
    Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
    combine 1 PolarPt v n
    pabc PolarPt v n
    pbcd)
        combineC _ = [Char] -> (PolarPt v n, PolarPt v n)
    forall a. HasCallStack => [Char] -> a
    error "combineC must be called on a list of length 2"

        -- The bezier end points are affine combinations of the bezier
        -- control points.
        bezierEnds :: [PolarPt v n]
    bezierEnds                   = ([(PolarPt v n, PolarPt v n)] -> PolarPt v n)
    -> [[(PolarPt v n, PolarPt v n)]] -> [PolarPt v n]
    forall a b. (a -> b) -> [a] -> [b]
    map [(PolarPt v n, PolarPt v n)] -> PolarPt v n
    forall (v :: * -> *) n.
    (Additive v, Fractional n, Ord n) =>
    [(PolarPt v n, PolarPt v n)] -> PolarPt v n
    combineE (Int
    -> [(PolarPt v n, PolarPt v n)] -> [[(PolarPt v n, PolarPt v n)]]
    forall a. Int -> [a] -> [[a]]
    windows 2 [(PolarPt v n, PolarPt v n)]
    bezierControls)
        combineE :: [(PolarPt v n, PolarPt v n)] -> PolarPt v n
    combineE [(_,pabb :: PolarPt v n
    pabb),(pbbc :: PolarPt v n
    pbbc,_)] = Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
    forall (v :: * -> *) n.
    (Additive v, Fractional n, Ord n) =>
    Int -> PolarPt v n -> PolarPt v n -> PolarPt v n
    combine 0 PolarPt v n
    pabb PolarPt v n
    pbbc
        combineE _ = [Char] -> PolarPt v n
    forall a. HasCallStack => [Char] -> a
    error "combineE must be called on a list of length 2"

        -- Finally, we actually put together the generated bezier segments.
        beziers :: [FixedSegment v n]
    beziers                      = ((PolarPt v n, PolarPt v n) -> [PolarPt v n] -> FixedSegment v n)
    -> [(PolarPt v n, PolarPt v n)]
    -> [[PolarPt v n]]
    -> [FixedSegment v n]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith (PolarPt v n, PolarPt v n) -> [PolarPt v n] -> FixedSegment v n
    forall (v :: * -> *) n.
    (PolarPt v n, PolarPt v n) -> [PolarPt v n] -> FixedSegment v n
    mkBezier (Int -> [(PolarPt v n, PolarPt v n)] -> [(PolarPt v n, PolarPt v n)]
    forall a. Int -> [a] -> [a]
    drop 1 [(PolarPt v n, PolarPt v n)]
    bezierControls) (Int -> [PolarPt v n] -> [[PolarPt v n]]
    forall a. Int -> [a] -> [[a]]
    windows 2 [PolarPt v n]
    bezierEnds)
          where
            mkBezier :: (PolarPt v n, PolarPt v n) -> [PolarPt v n] -> FixedSegment v n
    mkBezier (paab :: PolarPt v n
    paab,pabb :: PolarPt v n
    pabb) [paaa :: PolarPt v n
    paaa,pbbb :: PolarPt v n
    pbbb]
              = Point v n
    -> Point v n -> Point v n -> Point v n -> FixedSegment v n
    forall (v :: * -> *) n.
    Point v n
    -> Point v n -> Point v n -> Point v n -> FixedSegment v n
    FCubic (PolarPt v n -> Point v n
    forall (v :: * -> *) n. PolarPt v n -> Point v n
    unPP PolarPt v n
    paaa) (PolarPt v n -> Point v n
    forall (v :: * -> *) n. PolarPt v n -> Point v n
    unPP PolarPt v n
    paab) (PolarPt v n -> Point v n
    forall (v :: * -> *) n. PolarPt v n -> Point v n
    unPP PolarPt v n
    pabb) (PolarPt v n -> Point v n
    forall (v :: * -> *) n. PolarPt v n -> Point v n
    unPP PolarPt v n
    pbbb)
            mkBezier _ _ = [Char] -> FixedSegment v n
    forall a. HasCallStack => [Char] -> a
    error "mkBezier must be called on a list of length 2"

        -- Note that the above algorithm works in any dimension but is
        -- very specific to *cubic* splines.  This can of course be
        -- generalized to higher degree splines but keeping track of
        -- everything gets a bit more complicated; to be honest I am not
        -- quite sure how to do it.

    -- | Generate a uniform cubic B-spline from the given control points.
    --   The spline starts and ends at the first and last control points,
    --   and is tangent to the line to the second(-to-last) control point.
    --   It does not necessarily pass through any of the other control
    --   points.
    --
    --   <<diagrams/src_Diagrams_CubicSpline_Boehm_bsplineEx.svg#diagram=bsplineEx&width=300>>
    --
    --   > pts = map p2 [(0,0), (2,3), (5,-2), (-4,1), (0,3)]
    --   > spot = circle 0.2 # fc blue # lw none
    --   > bsplineEx = mconcat
    --   >   [ position (zip pts (repeat spot))
    --   >   , bspline pts
    --   >   ]
    --   >   # frame 0.5

    bspline :: (TrailLike t, V t ~ v, N t ~ n) => BSpline v n -> t
    bspline :: BSpline v n -> t
    bspline = Located [Segment Closed v n] -> t
    forall t. TrailLike t => Located [Segment Closed (V t) (N t)] -> t
    fromLocSegments (Located [Segment Closed v n] -> t)
    -> (BSpline v n -> Located [Segment Closed v n])
    -> BSpline v n
    -> t
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [Located (Segment Closed v n)] -> Located [Segment Closed v n]
    forall a. (Additive (V a), Num (N a)) => [Located a] -> Located [a]
    fixup ([Located (Segment Closed v n)] -> Located [Segment Closed v n])
    -> (BSpline v n -> [Located (Segment Closed v n)])
    -> BSpline v n
    -> Located [Segment Closed v n]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (FixedSegment v n -> Located (Segment Closed v n))
    -> [FixedSegment v n] -> [Located (Segment Closed v n)]
    forall a b. (a -> b) -> [a] -> [b]
    map FixedSegment v n -> Located (Segment Closed v n)
    forall n (v :: * -> *).
    (Num n, Additive v) =>
    FixedSegment v n -> Located (Segment Closed v n)
    fromFixedSeg ([FixedSegment v n] -> [Located (Segment Closed v n)])
    -> (BSpline v n -> [FixedSegment v n])
    -> BSpline v n
    -> [Located (Segment Closed v n)]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . BSpline v n -> [FixedSegment v n]
    forall (v :: * -> *) n.
    (Additive v, Fractional n, Num n, Ord n) =>
    BSpline v n -> [FixedSegment v n]
    bsplineToBeziers
      where
        fixup :: [Located a] -> Located [a]
    fixup []        = [] [a] -> Point (V [a]) (N [a]) -> Located [a]
    forall a. a -> Point (V a) (N a) -> Located a
    `at` Point (V [a]) (N [a])
    forall (f :: * -> *) a. (Additive f, Num a) => Point f a
    origin
        fixup (b1 :: Located a
    b1:rest :: [Located a]
    rest) = (Located a -> a
    forall a. Located a -> a
    unLoc Located a
    b1 a -> [a] -> [a]
    forall a. a -> [a] -> [a]
    : (Located a -> a) -> [Located a] -> [a]
    forall a b. (a -> b) -> [a] -> [b]
    map Located a -> a
    forall a. Located a -> a
    unLoc [Located a]
    rest) [a] -> Point (V [a]) (N [a]) -> Located [a]
    forall a. a -> Point (V a) (N a) -> Located a
    `at` Located a -> Point (V a) (N a)
    forall a. Located a -> Point (V a) (N a)
    loc Located a
    b1