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

-----------------------------------------------------------------------------
-- |
-- Module      :  Diagrams.CubicSpline.Boehm
--
-- 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)

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 b t x y = lerp ((t-a)/(b-a)) y 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 k = takeWhile ((==k) . length) . map (take k) . 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 k xs = replicate k (head xs) ++ xs ++ replicate k (last 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 { unPP :: Point v n, _knots :: [n] }

mkPolarPt :: Ord n => Point v n -> [n] -> PolarPt v n
mkPolarPt pt kts = PP pt (sort 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 k (PP pt1 kts1) (PP pt2 kts2)
= mkPolarPt
(affineCombo (head kts1) (last kts2) newKt pt1 pt2)
(newKt : drop 1 kts1)
where
newKt = kts2 !! 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 controls = beziers
where
n                            = length controls
numKnots                     = n + 2
knots                        = iterateN numKnots (+1/(fromIntegral numKnots - 1)) 0

-- The control points are P(a,b,c), P(b,c,d), P(c,d,e), and so on.
controls' = zipWith mkPolarPt (extend 2 controls) (windows 3 \$ extend 2 knots)

-- The bezier internal control points are affine combinations of
-- the spline control points.
bezierControls        = map combineC (windows 2 controls')
combineC [pabc, pbcd] = (combine 0 pabc pbcd, combine 1 pabc pbcd)
combineC _ = error "combineC must be called on a list of length 2"

-- The bezier end points are affine combinations of the bezier
-- control points.
bezierEnds                   = map combineE (windows 2 bezierControls)
combineE [(_,pabb),(pbbc,_)] = combine 0 pabb pbbc
combineE _ = error "combineE must be called on a list of length 2"

-- Finally, we actually put together the generated bezier segments.
beziers                      = zipWith mkBezier (drop 1 bezierControls) (windows 2 bezierEnds)
where
mkBezier (paab,pabb) [paaa,pbbb]
= FCubic (unPP paaa) (unPP paab) (unPP pabb) (unPP pbbb)
mkBezier _ _ = 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 = fromLocSegments . fixup . map fromFixedSeg . bsplineToBeziers
where
fixup []        = [] `at` origin
fixup (b1:rest) = (unLoc b1 : map unLoc rest) `at` loc b1
```