{-# LANGUAGE MultiParamTypeClasses, StandaloneDeriving, FlexibleContexts, UndecidableInstances, TypeFamilies, ParallelListComp #-} 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 kts cps@ creates a B-spline with the given knot vector and control -- points. The degree is automatically inferred as the difference between the -- number of spans in the knot vector (@numKnots kts - 1@) and the number of -- control points (@length cps@). 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