```{-# LANGUAGE FlexibleContexts #-}
module Math.Spline.BSpline.Internal
(BSpline(..), mapControlPoints, evalBSpline, insertKnot, deBoor) where

import Math.Spline.Knots

import Data.List (zipWith4)
import Data.Monoid
import Data.VectorSpace

data BSpline v = Spline
{ degree        :: !Int
, knotVector    :: Knots (Scalar v)
, controlPoints :: [v]
}

mapControlPoints f spline = spline
{ controlPoints = map f (controlPoints spline)
, knotVector = knotVector spline
}

evalBSpline spline = head . last . deBoor spline

-- |Insert one knot into a 'BSpline'
insertKnot
:: (VectorSpace a, Ord (Scalar a), Fractional (Scalar a)) =>
BSpline a -> Scalar a -> BSpline a
insertKnot spline x = spline
{ knotVector    = knotVector spline `mappend` knot x
, controlPoints = zipWith4 (interp x) us (drop p us) ds (tail ds)
}
where
us = knots (knotVector spline)
p  = degree spline
ds = extend (controlPoints spline)

-- duplicate the endpoints of a list; for example,
-- extend [1..5] -> [1,1,2,3,4,5,5]
extend []       = []
extend (x:xs)   = x : extend' x xs
where   extend' x []      = [x,x]
extend' x (x':xs) = x:   extend' x' xs

deBoor spline x = go us (controlPoints spline)
where
us = knots (knotVector spline)

-- Upper endpoints of the intervals are the same for
-- each row in the table (they just line up differently
-- with the lower endpoints):
uHi = drop (degree spline + 1) us

-- On each pass, the lower endpoints of the
-- interpolation intervals advance and the new
-- coefficients are given by linear interpolation
-- on the current intervals:
go       _ [] = []
go (_:uLo) ds = ds : go uLo ds'
where
ds' = zipWith4 (interp x) uLo uHi
ds (tail ds)

interp x x0 x1 y0 y1
|  x <  x0  = y0
|  x >= x1  = y1
| otherwise = lerp y0 y1 a
where
a = (x - x0) / (x1 - x0)

```