module Math.Spline.BSpline
( BSpline
, bSpline
, evalBSpline
, insertKnot
, splitBSpline
, differentiateBSpline, integrateBSpline
) where
import Math.Spline.Knots
import Math.Spline.BSpline.Internal
import Data.Maybe (fromMaybe)
import Data.VectorSpace
bSpline :: Knots (Scalar a) -> [a] -> BSpline a
bSpline _ [] = error "bSpline: no control points"
bSpline kts cps = fromMaybe (error "bSpline: too few knots") (maybeSpline kts cps)
maybeSpline :: Knots (Scalar a) -> [a] -> Maybe (BSpline a)
maybeSpline kts cps
| n > m = Nothing
| otherwise = Just (Spline (m n) kts cps)
where
n = length cps
m = numKnots kts 1
deriving instance (Eq (Scalar v), Eq v) => Eq (BSpline v)
deriving instance (Ord (Scalar v), Ord v) => Ord (BSpline v)
instance (Show (Scalar v), Show v) => Show (BSpline v) where
showsPrec p (Spline _ kts cps) = showParen (p>10)
( showString "bSpline "
. showsPrec 11 kts
. showChar ' '
. showsPrec 11 cps
)
differentiateBSpline
:: (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => BSpline v -> BSpline v
differentiateBSpline spline
| numKnots ks < 2 = spline
| numKnots ks == 2 = bSpline ks [zeroV]
| otherwise = bSpline ks' ds'
where
ks' = mkKnots . init . tail $ ts
ds' = zipWith (*^) (tail cs) (zipWith (^-^) (tail ds) ds)
ks = knotVector spline; ts = knots ks
ds = controlPoints spline
p = degree spline
cs = [fromIntegral p / (t1 t0) | (t0,t1) <- spans p ts]
integrateBSpline
:: (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => BSpline v -> BSpline v
integrateBSpline spline = bSpline (mkKnots ts') (scanl (^+^) zeroV ds')
where
ds' = zipWith (*^) cs (controlPoints spline)
ts = knots (knotVector spline)
ts' = head ts : ts ++ [last ts]
p = degree spline + 1
cs = [(t1 t0) / fromIntegral p | (t0,t1) <- spans p ts]
spans n xs = zip xs (drop n xs)
splitBSpline
:: (VectorSpace v, Ord (Scalar v), Fractional (Scalar v)) =>
BSpline v -> Scalar v -> Maybe (BSpline v, BSpline v)
splitBSpline spline@(Spline p kv ds) t
| inDomain = Just (Spline p (mkKnots us0) ds0, Spline p (mkKnots us1) ds1)
| otherwise = Nothing
where
inDomain = case knotDomain kv p of
Nothing -> False
Just (t0, t1) -> t >= t0 || t <= t1
us = knots kv
dss = deBoor spline t
us0 = takeWhile (<t) us ++ replicate (p+1) t
ds0 = trimTo (drop (p+1) us0) (map head dss)
us1 = replicate (p+1) t ++ dropWhile (<=t) us
ds1 = reverse (trimTo (drop (p+1) us1) (map last dss))
trimTo list xs = zipWith const xs list