splines 0.1 → 0.3
raw patch · 11 files changed
+419/−85 lines, 11 filesdep +polynomialdep +vector
Dependencies added: polynomial, vector
Files
- splines.cabal +9/−4
- src/Math/NURBS.hs +5/−4
- src/Math/Spline.hs +1/−1
- src/Math/Spline/BSpline.hs +17/−13
- src/Math/Spline/BSpline/Internal.hs +20/−18
- src/Math/Spline/BSpline/Reference.hs +93/−0
- src/Math/Spline/BezierCurve.hs +12/−10
- src/Math/Spline/Class.hs +3/−2
- src/Math/Spline/ISpline.hs +8/−8
- src/Math/Spline/Knots.hs +242/−12
- src/Math/Spline/MSpline.hs +9/−13
splines.cabal view
@@ -1,5 +1,5 @@ name: splines-version: 0.1+version: 0.3 stability: provisional cabal-version: >= 1.6@@ -19,18 +19,23 @@ and/or feature requests are welcome. source-repository head- type: darcs- location: http://code.haskell.org/~mokus/splines/+ type: git+ location: git://github.com/mokus0/splines.git Library hs-source-dirs: src exposed-modules: Math.Spline Math.Spline.BezierCurve Math.Spline.BSpline+ Math.Spline.BSpline.Reference Math.Spline.Class Math.Spline.ISpline Math.Spline.Knots Math.Spline.MSpline Math.NURBS other-modules: Math.Spline.BSpline.Internal- build-depends: base >= 3 && < 5, containers, vector-space+ build-depends: base >= 3 && < 5,+ containers,+ polynomial,+ vector,+ vector-space
src/Math/NURBS.hs view
@@ -7,6 +7,7 @@ , splitNURBS ) where +import qualified Data.Vector as V import Data.VectorSpace import Math.Spline.Class (Spline, toBSpline) import Math.Spline.BSpline.Internal@@ -28,19 +29,19 @@ nurbs :: (VectorSpace v, Scalar v ~ w, VectorSpace w, Scalar w ~ w)- => Knots (Scalar v) -> [(w, v)] -> NURBS v+ => Knots (Scalar v) -> V.Vector (w, v) -> NURBS v nurbs kts cps = NURBS (bSpline kts cps) -- |Constructs the homogeneous-coordinates B-spline that corresponds to this -- NURBS curve nurbsAsSpline (NURBS spline) = spline - { controlPoints = map homogenize (controlPoints spline) }+ { controlPoints = V.map homogenize (controlPoints spline) } where homogenize (w,v) = (w, w *^ v) -- |Constructs the NURBS curve corresponding to a homogeneous-coordinates B-spline splineAsNURBS spline = NURBS spline - { controlPoints = map unHomogenize (controlPoints spline) }+ { controlPoints = V.map unHomogenize (controlPoints spline) } where unHomogenize (w,v) = (w, recip w *^ v) @@ -67,7 +68,7 @@ nurbsKnotVector :: Scalar v ~ Scalar (Scalar v) => NURBS v -> Knots (Scalar v) nurbsKnotVector (NURBS spline) = knotVector spline -nurbsControlPoints :: NURBS v -> [(Scalar v, v)]+nurbsControlPoints :: NURBS v -> V.Vector (Scalar v, v) nurbsControlPoints (NURBS spline) = controlPoints spline splitNURBS :: (VectorSpace v, Scalar v ~ w,
src/Math/Spline.hs view
@@ -1,5 +1,5 @@ module Math.Spline- ( Spline(..)+ ( Spline(..), ControlPoints(..) , Knots, mkKnots, knots
src/Math/Spline/BSpline.hs view
@@ -13,21 +13,25 @@ import Data.Maybe (fromMaybe) import Data.VectorSpace+import qualified Data.Vector as V -- |@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)+bSpline :: Knots (Scalar a) -> V.Vector a -> BSpline a+bSpline kts cps+ | V.null cps = error "bSpline: no control points"+ | otherwise = fromMaybe+ (error "bSpline: too few knots")+ (maybeSpline kts cps) -maybeSpline :: Knots (Scalar a) -> [a] -> Maybe (BSpline a)+maybeSpline :: Knots (Scalar a) -> V.Vector a -> Maybe (BSpline a) maybeSpline kts cps | n > m = Nothing | otherwise = Just (Spline (m - n) kts cps) where- n = length cps+ n = V.length cps m = numKnots kts - 1 deriving instance (Eq (Scalar v), Eq v) => Eq (BSpline v)@@ -44,27 +48,27 @@ :: (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => BSpline v -> BSpline v differentiateBSpline spline | numKnots ks < 2 = spline- | numKnots ks == 2 = bSpline ks [zeroV]+ | numKnots ks == 2 = bSpline ks (V.singleton zeroV) | otherwise = bSpline ks' ds' where ks' = mkKnots . init . tail $ ts- ds' = zipWith (*^) (tail cs) (zipWith (^-^) (tail ds) ds)+ ds' = V.zipWith (*^) (V.tail cs) (V.zipWith (^-^) (V.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]+ cs = V.fromList [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')+integrateBSpline spline = bSpline (mkKnots ts') (V.scanl (^+^) zeroV ds') where- ds' = zipWith (*^) cs (controlPoints spline)+ ds' = V.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]+ cs = V.fromList [(t1 - t0) / fromIntegral p | (t0,t1) <- spans p ts] spans n xs = zip xs (drop n xs) @@ -83,9 +87,9 @@ dss = deBoor spline t us0 = takeWhile (<t) us ++ replicate (p+1) t- ds0 = trimTo (drop (p+1) us0) (map head dss)+ ds0 = V.fromList (trimTo (drop (p+1) us0) (map V.head dss)) us1 = replicate (p+1) t ++ dropWhile (<=t) us- ds1 = reverse (trimTo (drop (p+1) us1) (map last dss))+ ds1 = V.reverse (V.fromList (trimTo (drop (p+1) us1) (map V.last dss))) trimTo list xs = zipWith const xs list
src/Math/Spline/BSpline/Internal.hs view
@@ -4,22 +4,23 @@ import Math.Spline.Knots -import Data.List (zipWith4) import Data.Monoid+import Data.Vector as V import Data.VectorSpace+import Prelude as P -data BSpline v = Spline+data BSpline t = Spline { degree :: !Int- , knotVector :: Knots (Scalar v)- , controlPoints :: [v]+ , knotVector :: Knots (Scalar t)+ , controlPoints :: Vector t } mapControlPoints f spline = spline- { controlPoints = map f (controlPoints spline)+ { controlPoints = V.map f (controlPoints spline) , knotVector = knotVector spline } -evalBSpline spline = head . last . deBoor spline+evalBSpline spline = V.head . P.last . deBoor spline -- |Insert one knot into a 'BSpline' insertKnot@@ -27,39 +28,40 @@ 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)+ , controlPoints = V.zipWith4 (interp x) us (V.drop p us) ds (V.tail ds) } where- us = knots (knotVector spline)+ us = knotsVector (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+extend vec+ | V.null vec = V.empty+ | otherwise = V.cons (V.head vec) (V.snoc vec (V.last vec)) deBoor spline x = go us (controlPoints spline) where- us = knots (knotVector spline)+ us = knotsVector (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+ uHi = V.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'+ go us ds + | V.null ds = []+ | otherwise = ds : go uLo ds' where- ds' = zipWith4 (interp x) uLo uHi- ds (tail ds)+ uLo = V.tail us+ ds' = V.zipWith4 (interp x) uLo uHi+ ds (V.tail ds) interp x x0 x1 y0 y1 | x < x0 = y0
+ src/Math/Spline/BSpline/Reference.hs view
@@ -0,0 +1,93 @@+{-# LANGUAGE ParallelListComp #-}+-- |Reference implementation of B-Splines; very inefficient but \"obviously\"+-- correct.+module Math.Spline.BSpline.Reference+ ( bases+ , basisFunctions+ , basisPolynomials+ , basisPolynomialsAt+ ) where++import Math.Spline.Knots+import Math.Polynomial (Poly)+import qualified Math.Polynomial as Poly++ind True = 1+ind False = 0++bases :: (Fractional a, Ord a) => Knots a -> a -> [[a]]+bases kts x = coxDeBoor interp initial kts+ where+ initial = + [ ind (t_j <= x && x < t_jp1)+ | (t_j, t_jp1) <- knotSpans kts 1+ ]+ interp t_j d0 b_nm1_j t_jpnp1 d1 b_nm1_jp1+ = (if d0 == 0 then 0 else (x - t_j) / d0) * b_nm1_j+ + (if d1 == 0 then 0 else (t_jpnp1 - x) / d1) * b_nm1_jp1++-- Alternate version constructing table of functions rather than computing+-- table of values+basisFunctions :: (Fractional a, Ord a) => Knots a -> [[a -> a]]+basisFunctions kts = coxDeBoor interp initial kts+ where+ initial = + [ \x -> ind (t_j <= x && x < t_jp1)+ | (t_j, t_jp1) <- knotSpans kts 1+ ]+ interp t_j d0 b_nm1_j t_jpnp1 d1 b_nm1_jp1 x+ = (if d0 == 0 then 0 else (x - t_j) / d0) * b_nm1_j x+ + (if d1 == 0 then 0 else (t_jpnp1 - x) / d1) * b_nm1_jp1 x++-- compute all the basis polynomials for a knot vector, ordered by knot span.+basisPolynomials :: (Fractional a, Ord a) => Knots a -> [[[Poly a]]]+basisPolynomials kts+ | isEmpty kts = []+ | otherwise = [basisPolynomialsAt kts kt | kt <- init (distinctKnots kts)]++-- compute all the basis polynomials for the knot span containing a given location.+basisPolynomialsAt :: (Fractional a, Ord a) => Knots a -> a -> [[Poly a]]+basisPolynomialsAt kts x = coxDeBoor interp initial kts+ where+ indPoly True = Poly.one+ indPoly False = Poly.zero+ + initial = + [ indPoly (t_j <= x && x < t_jp1)+ | (t_j, t_jp1) <- knotSpans kts 1+ ]+ interp t_j d0 b_nm1_j t_jpnp1 d1 b_nm1_jp1+ = (if d0 == 0 then Poly.zero else (Poly.x - Poly.constPoly t_j) / d0) * b_nm1_j+ + (if d1 == 0 then Poly.zero else (Poly.constPoly t_jpnp1 - Poly.x) / d1) * b_nm1_jp1+ where+ infixl 6 +, -+ p + q = Poly.addPoly p q+ p - q = p + (Poly.negatePoly q)+ + infixl 7 *, /+ p * q = Poly.multPoly p q+ p / s = Poly.scalePoly (recip s) p++-- This is a straightforward implementation of the Cox-De Boor recursion scheme+-- generalized in a slightly strange way; the initial vector is a parameter +-- and the actual computation of the recursion step is a function parameter.+-- The purpose is to allow the same recursion to be applied when computing basis+-- function values and basis polynomials.+coxDeBoor interp initial kts = table+ where+ ts = knots kts+ table = initial :+ [ [ interp t_j d0 b_nm1_j t_jpnp1 d1 b_nm1_jp1+ | (b_nm1_j, b_nm1_jp1) <- spans 1 prevBasis+ | (d0, d1) <- spans 1 (spanDiffs n ts)+ | (t_j, t_jpnp1) <- spans (n+1) ts+ ]+ | prevBasis <- takeWhile (not . null) table+ | n <- [1..]+ ]++spans :: Int -> [a] -> [(a,a)]+spans = spansWith (,)+spanDiffs :: Num a => Int -> [a] -> [a]+spanDiffs = spansWith subtract+spansWith f n ts = zipWith f ts (drop n ts)
src/Math/Spline/BezierCurve.hs view
@@ -9,17 +9,18 @@ import Math.Spline.Knots import Control.Applicative+import qualified Data.Vector as V import Data.VectorSpace -- |A BezierCurve curve on @0 <= x <= 1@.-data BezierCurve v = BezierCurve !Int [v] deriving (Eq, Ord)+data BezierCurve t = BezierCurve !Int !(V.Vector t) deriving (Eq, Ord) -- |Construct a Bezier curve from a list of control points. The degree -- of the curve is one less than the number of control points.-bezierCurve :: [v] -> BezierCurve v+bezierCurve :: V.Vector t -> BezierCurve t bezierCurve cs- | null cs = error "bezierCurve: no control points given"- | otherwise = BezierCurve (length cs - 1) cs+ | V.null cs = error "bezierCurve: no control points given"+ | otherwise = BezierCurve (V.length cs - 1) cs instance Show v => Show (BezierCurve v) where showsPrec p (BezierCurve _ cs) = showParen (p>10)@@ -29,7 +30,7 @@ instance (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline BezierCurve v where splineDomain (BezierCurve _ _) = Just (0,1)- evalSpline (BezierCurve _ cs) = head . last . deCasteljau cs+ evalSpline (BezierCurve _ cs) = V.head . last . deCasteljau cs splineDegree (BezierCurve p _) = p knotVector (BezierCurve p _) = fromList [(0, p+1), (1, p+1)] toBSpline = bSpline <$> knotVector <*> controlPoints@@ -37,9 +38,10 @@ instance Spline BezierCurve v => ControlPoints BezierCurve v where controlPoints (BezierCurve _ cs) = cs -deCasteljau :: VectorSpace v => [v] -> Scalar v -> [[v]]-deCasteljau [] t = []-deCasteljau cs t = cs : deCasteljau (zipWith interp cs (tail cs)) t+deCasteljau :: VectorSpace v => V.Vector v -> Scalar v -> [V.Vector v]+deCasteljau cs t+ | V.null cs = []+ | otherwise = cs : deCasteljau (V.zipWith interp cs (V.tail cs)) t where interp x0 x1 = lerp x0 x1 t @@ -52,6 +54,6 @@ -- splitBezierCurve :: VectorSpace v => BezierCurve v -> Scalar v -> (BezierCurve v, BezierCurve v) splitBezierCurve (BezierCurve n cs) t = - ( BezierCurve n (map head css)- , BezierCurve n (reverse (map last css))+ ( BezierCurve n (V.fromList (map V.head css))+ , BezierCurve n (V.reverse (V.fromList (map V.last css))) ) where css = deCasteljau cs t
src/Math/Spline/Class.hs view
@@ -5,6 +5,7 @@ import Math.Spline.Knots import qualified Math.Spline.BSpline.Internal as BSpline +import qualified Data.Vector as V import Data.VectorSpace -- |A spline is a piecewise polynomial vector-valued function. The necessary@@ -32,10 +33,10 @@ toBSpline :: s v -> BSpline.BSpline v class Spline s v => ControlPoints s v where- controlPoints :: s v -> [v]+ controlPoints :: s v -> V.Vector v instance (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline BSpline.BSpline v where- evalSpline spline = head . last . BSpline.deBoor spline+ evalSpline spline = V.head . last . BSpline.deBoor spline splineDegree = BSpline.degree knotVector = BSpline.knotVector toBSpline = id
src/Math/Spline/ISpline.hs view
@@ -12,7 +12,7 @@ import Math.Spline.BSpline import Math.Spline.Class import Math.Spline.Knots-+import qualified Data.Vector as V import Data.VectorSpace -- |The I-Spline basis functions are the integrals of the M-splines, or@@ -23,7 +23,7 @@ data ISpline v = ISpline { iSplineDegree :: !Int , iSplineKnotVector :: Knots (Scalar v)- , iSplineControlPoints :: [v]+ , iSplineControlPoints :: !(V.Vector v) } deriving instance (Eq (Scalar v), Eq v) => Eq (ISpline v)@@ -41,19 +41,19 @@ -- 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@).-iSpline :: Knots (Scalar a) -> [a] -> ISpline a+iSpline :: Knots (Scalar a) -> V.Vector a -> ISpline a iSpline kts cps | n > m = error "iSpline: too few knots" | otherwise = ISpline (m - n) kts cps where- n = length cps+ n = V.length cps m = numKnots kts - 1 instance (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline ISpline v where splineDegree = (1 +) . iSplineDegree knotVector spline = mkKnots (head ts : ts ++ [last ts]) where ts = knots (iSplineKnotVector spline)- toBSpline spline = bSpline (knotVector spline) (scanl (^+^) zeroV cs)+ toBSpline spline = bSpline (knotVector spline) (V.scanl (^+^) zeroV cs) where cs = iSplineControlPoints spline instance Spline ISpline v => ControlPoints ISpline v where@@ -63,12 +63,12 @@ toISpline = fromBSpline . toBSpline fromBSpline spline- | head ds == zeroV - && numKnots ks >= 2 = iSpline (mkKnots (init (tail ts))) (tail ds')+ | V.head ds == zeroV + && numKnots ks >= 2 = iSpline (mkKnots (init (tail ts))) (V.tail ds') | otherwise = iSpline (mkKnots (init ts )) ds' where ks = knotVector spline ts = knots ks ds = controlPoints spline - ds' = zipWith (^-^) ds (zeroV:ds)+ ds' = V.zipWith (^-^) ds (V.cons zeroV ds)
src/Math/Spline/Knots.hs view
@@ -1,26 +1,56 @@+{-# LANGUAGE TypeFamilies #-} module Math.Spline.Knots ( Knots+ , empty, isEmpty+ , knot, multipleKnot , mkKnots, fromList - , knots, numKnots- , toList, distinctKnots, numDistinctKnots+ , numKnots, lookupKnot+ , toList, numDistinctKnots, lookupDistinctKnot + , knots, knotsVector+ , distinctKnots, distinctKnotsVector+ + , toMap+ , fromMap+ + , toVector+ , fromVector+ + , splitLookup+ , takeKnots, dropKnots, splitKnotsAt+ , takeDistinctKnots, dropDistinctKnots, splitDistinctKnotsAt+ + , maxMultiplicity , knotMultiplicity, setKnotMultiplicity + , fromAscList, fromDistinctAscList+ , valid+ + , knotSpan+ , knotsInSpan+ , knotSpans , knotDomain+ + , uniform ) where -import Prelude hiding (sum)-import Data.Foldable (Foldable(foldMap), sum)+import Prelude hiding (sum, maximum)+import Control.Monad (guard)+import Data.Foldable (Foldable(foldMap), sum, maximum) import qualified Data.Map as M import Data.Monoid (Monoid(..)) import Data.Maybe (fromMaybe)+import qualified Data.Set as S (Set)+import qualified Data.Vector as V+import Data.VectorSpace -- |Knot vectors - multisets of points in a 1-dimensional space. data Knots a = Knots !Int (M.Map a Int) deriving (Eq, Ord) instance Show a => Show (Knots a) where+ showsPrec p ks@(Knots 0 _) = showString "empty" showsPrec p ks@(Knots 1 _) = showParen (p > 10) ( showString "knot " . showsPrec 11 (head $ knots ks)@@ -31,7 +61,7 @@ ) instance (Ord a) => Monoid (Knots a) where- mempty = Knots 0 M.empty+ mempty = empty mappend (Knots n1 v1) (Knots n2 v2) = Knots (n1 + n2) (M.filter (/=0) (M.unionWith (+) v1 v2)) @@ -39,6 +69,14 @@ foldMap f = foldMap f . knots +-- |An empty knot vector+empty :: Knots a+empty = Knots 0 M.empty++isEmpty :: Knots a -> Bool+isEmpty (Knots 0 _) = True+isEmpty _ = False+ -- |Create a knot vector consisting of one knot. knot :: Ord a => a -> Knots a knot x = multipleKnot x 1@@ -59,11 +97,41 @@ fromList ks = Knots (sum kMap) kMap where kMap = M.fromListWith (+) (filter ((>0).snd) ks) +-- |Create a knot vector consisting of all the knots and corresponding +-- multiplicities in a list ordered by the knots' 'Ord' instance. The+-- ordering precondition is not checked.+fromAscList :: Eq k => [(k, Int)] -> Knots k+fromAscList ks = Knots (sum kMap) kMap+ where kMap = M.fromAscListWith (+) (filter ((>0).snd) ks)++-- |Create a knot vector consisting of all the knots and corresponding +-- multiplicities in a list ordered by the knots' 'Ord' instance with no+-- duplicates. The preconditions are not checked.+fromDistinctAscList :: [(k, Int)] -> Knots k+fromDistinctAscList ks = Knots (sum kMap) kMap+ where kMap = M.fromDistinctAscList (filter ((>0).snd) ks)++fromMap :: M.Map k Int -> Knots k+fromMap ks = Knots (sum kMap) kMap+ where+ kMap = mFilter (>0) ks+ -- filter is monotonic, I have no idea why M.filter requires Ord on the key+ mFilter p = M.fromDistinctAscList . filter (p.snd) . M.toAscList++fromVector :: Ord k => V.Vector (k,Int) -> Knots k+fromVector = fromList . V.toList+ -- |Returns a list of all distinct knots in ascending order along with -- their multiplicities. toList :: Knots k -> [(k, Int)]-toList (Knots _ ks) = M.toList ks+toList = M.toList . toMap +toVector :: Knots k -> V.Vector (k, Int)+toVector = V.fromList . toList++toMap :: Knots k -> M.Map k Int+toMap (Knots _ ks) = ks+ -- |Returns the number of knots (not necessarily distinct) in a knot vector. numKnots :: Knots t -> Int numKnots (Knots n _) = n@@ -72,14 +140,138 @@ numDistinctKnots :: Knots t -> Int numDistinctKnots (Knots _ ks) = M.size ks +maxMultiplicity :: Knots t -> Int+maxMultiplicity (Knots 0 _) = 0+maxMultiplicity (Knots _ ks) = maximum ks++lookupKnot :: Int -> Knots a -> Maybe a+lookupKnot k kts+ | k < 0 = Nothing+ | k < numKnots kts = fmap fst mbKt+ | otherwise = Nothing+ where (_, mbKt, _) = splitLookup k kts++lookupDistinctKnot :: Int -> Knots a -> Maybe a+lookupDistinctKnot k (Knots _ ks)+ | k < 0 = Nothing+ | k < M.size ks = Just (fst (M.elemAt k ks))+ | otherwise = Nothing++-- |@splitLookup n kts@: Split a knot vector @kts@ into 3 parts @(pre, mbKt, post)@+-- such that:+-- +-- * All the keys in @pre@, @mbKt@ (viewed as a knot vector of either 0+-- or 1 knot), and @post@ are disjoint and ordered+-- * Putting the 3 parts back together yields exactly the original knot vector+-- * The @n@'th knot, if one exists, will be in @mbKt@ along with its multiplicity+--+splitLookup :: Int -> Knots a -> (Knots a, Maybe (a, Int), Knots a)+splitLookup k (Knots n ks) = scan 0 M.empty n ks+ where+ -- The general plan: iteratively pull the smallest knot out of "post",+ -- either moving it to "pre" or terminating by returning it along with+ -- current values of "pre" and "post"+ + -- invariants:+ -- nPre = sum pre+ -- nPost = sum post+ -- M.union pre post = ks+ -- every key in pre < every key in post+ scan nPre pre nPost post+ | nPost <= 0 = (Knots nPre pre, Nothing, Knots nPost post)+ | nPre + m > k = (Knots nPre pre, Just kt, Knots nNewPost newPost)+ | otherwise = scan (nPre + m) (pre `ascSnoc` kt) nNewPost newPost+ where+ Just (kt@(x,m), newPost) = M.minViewWithKey post+ nNewPost = nPost - m+ done x = (Knots nPre pre, x, Knots nPost post)++-- Prepend or append an element to a map, without checking the precondition+-- that the new pair's key is less than (greater than, resp.) all keys in +-- the map.+ascCons x m = M.fromDistinctAscList (x : M.toAscList m)+ascSnoc m x = M.fromDistinctAscList (M.toAscList m ++ [x])++-- Prepend or append an knot to a knot vector, without checking the+-- precondition that the new knot's location is less than (greater than,+-- resp.) all knots in the vector.+ascConsKnot (_,0) kts = kts+ascConsKnot kt@(k,m) (Knots n ks) = Knots (n+m) (kt `ascCons` ks)++ascSnocKnot kts (_,0) = kts+ascSnocKnot (Knots n ks) kt@(k,m) = Knots (n+m) (ks `ascSnoc` kt)++clamp lo hi = max lo . min hi++dropKnots :: Int -> Knots a -> Knots a+dropKnots k kts = fromMaybe post $ do+ (x,xAvail) <- mbKt+ let xWanted = numKnots kts - (numKnots post + k)+ + return ((x, clamp 0 xAvail xWanted) `ascConsKnot` post)+ where+ (pre, mbKt, post) = splitLookup k kts++takeKnots :: Int -> Knots a -> Knots a+takeKnots k kts = fromMaybe pre $ do+ (x,xAvail) <- mbKt+ let xWanted = k - numKnots pre+ + return (pre `ascSnocKnot` (x, clamp 0 xAvail xWanted))+ where+ (pre, mbKt, post) = splitLookup k kts++splitKnotsAt :: Int -> Knots a -> (Knots a, Knots a)+splitKnotsAt k kts = fromMaybe (pre, post) $ do+ (x,xAvail) <- mbKt+ let xWanted = k - numKnots pre+ xTaken = clamp 0 xAvail xWanted+ + return ( pre `ascSnocKnot` (x,xTaken)+ , (x, xAvail - xTaken) `ascConsKnot` post+ )+ where+ (pre, mbKt, post) = splitLookup k kts+++takeDistinctKnots :: Int -> Knots a -> Knots a+takeDistinctKnots k (Knots n ks) = Knots (sum kMap) kMap+ where+ kMap = M.fromDistinctAscList (take k (M.toAscList ks))++dropDistinctKnots :: Int -> Knots a -> Knots a+dropDistinctKnots k (Knots n ks) = Knots (sum kMap) kMap+ where+ kMap = M.fromDistinctAscList (drop k (M.toAscList ks))++splitDistinctKnotsAt :: Int -> Knots a -> (Knots a, Knots a)+splitDistinctKnotsAt k (Knots n ks) = (Knots sz1 kMap1, Knots (n - sz1) kMap2)+ where+ (ks1, ks2) = splitAt k (M.toAscList ks)+ kMap1 = M.fromDistinctAscList ks1+ kMap2 = M.fromDistinctAscList ks2+ sz1 = sum kMap1+ -- |Returns a list of all knots (not necessarily distinct) of a knot vector in ascending order knots :: Knots t -> [t] knots (Knots _ ks) = concat [replicate n k | (k,n) <- M.toAscList ks] +-- |Returns a vector of all knots (not necessarily distinct) of a knot vector in ascending order+knotsVector :: Knots t -> V.Vector t+knotsVector (Knots _ ks) = V.concat [V.replicate n k | (k,n) <- M.toAscList ks]+ -- |Returns a list of all distinct knots of a knot vector in ascending order distinctKnots :: Knots t -> [t] distinctKnots (Knots _ ks) = M.keys ks +-- |Returns a vector of all distinct knots of a knot vector in ascending order+distinctKnotsVector :: Knots t -> V.Vector t+distinctKnotsVector = V.fromList . distinctKnots++-- |Returns a 'S.Set' of all distinct knots of a knot vector+distinctKnotsSet :: Knots k -> S.Set k+distinctKnotsSet (Knots _ ks) = M.keysSet ks+ -- |Looks up the multiplicity of a knot (which is 0 if the point is not a knot) knotMultiplicity :: (Ord k) => k -> Knots k -> Int knotMultiplicity k (Knots _ ks) = fromMaybe 0 (M.lookup k ks)@@ -93,16 +285,54 @@ where n' = knotMultiplicity k (Knots m ks) --- |@knotDomain kts p@ return the domain of a B-spline or NURBS with knot+-- |Check the internal consistency of a knot vector+valid :: Ord k => Knots k -> Bool+valid (Knots n ks) = and+ [ M.valid ks+ , n == sum ks+ , all (>0) (M.elems ks)+ ]++-- |@knotSpan kts i j@ returns the knot span extending from the @i@'th knot+-- to the @j@'th knot, if @i <= j@ and both knots exist.+knotSpan :: Knots a -> Int -> Int -> Maybe (a, a)+knotSpan ks i j = do+ guard (i <= j)+ lo <- lookupKnot i ks+ hi <- lookupKnot j ks+ return (lo,hi)++-- |@knotsInSpan kts i j@ returns the knots in the knot span extending from+-- the @i@'th knot to the @j@'th knot+knotsInSpan :: Knots a -> Int -> Int -> Knots a+knotsInSpan kts i j = takeKnots (j - i) (dropKnots i kts)++-- |@knotSpans kts width@ returns all knot spans of a given width in+-- ascending order.+--+-- For example, @knotSpans (mkKnots [1..5]) 2@ yields @[(1,3), (2,4), (3,5)]@.+knotSpans :: Knots a -> Int -> [(a,a)]+knotSpans ks w+ | w <= 0 = error "knotSpans: width must be positive"+ | otherwise = zip kts (drop w kts)+ where kts = knots ks++-- |@knotDomain kts p@ returns the domain of a B-spline or NURBS with knot -- vector @kts@ and degree @p@. This is the subrange spanned by all -- except the first and last @p@ knots. Outside this domain, the spline -- does not have a complete basis set. De Boor's algorithm assumes that -- the basis functions sum to 1, which is only true on this range, and so -- this is also precisely the domain on which de Boor's algorithm is valid. knotDomain :: Knots a -> Int -> Maybe (a,a)-knotDomain ks@(Knots n _) p - | n > 2*p = Just (head (drop p kts), head (drop p (reverse kts)))- | otherwise = Nothing- where- kts = knots ks+knotDomain ks@(Knots n _) p = knotSpan ks p (n-p-1) +-- |@uniform deg nPts (lo,hi)@ constructs a uniformly-spaced knot vector over+-- the interval from @lo@ to @hi@ which, when used to construct a B-spline +-- with @nPts@ control points will yield a clamped spline with degree @deg@.+uniform :: (Ord s, Fractional s) => Int -> Int -> (s,s) -> Knots s+uniform deg nPts (lo,hi) = ends `mappend` internal+ where+ ends = fromList [(lo,deg), (hi,deg)]+ n = nPts + deg - numKnots ends+ f i = (fromIntegral i * lo + fromIntegral (n - i) * hi) / fromIntegral n+ internal = mkKnots [f i | i <- [0..n]]
src/Math/Spline/MSpline.hs view
@@ -12,7 +12,7 @@ import Math.Spline.BSpline import Math.Spline.Class import Math.Spline.Knots-+import qualified Data.Vector as V import Data.VectorSpace -- |M-Splines are B-splines normalized so that the integral of each basis @@ -20,7 +20,7 @@ data MSpline v = MSpline { mSplineDegree :: !Int , mSplineKnotVector :: Knots (Scalar v)- , mSplineControlPoints :: [v]+ , mSplineControlPoints :: !(V.Vector v) } deriving instance (Eq (Scalar v), Eq v) => Eq (MSpline v)@@ -38,15 +38,15 @@ -- 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@).-mSpline :: Knots (Scalar a) -> [a] -> MSpline a+mSpline :: Knots (Scalar a) -> V.Vector a -> MSpline a mSpline kts cps | n > m = error "mSpline: too few knots" | otherwise = MSpline (m - n) kts cps where- n = length cps+ n = V.length cps m = numKnots kts - 1 -spans n xs = zip xs (drop n xs)+spans n xs = V.zip xs (V.drop n xs) instance (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline MSpline v where splineDegree = mSplineDegree@@ -54,10 +54,8 @@ toBSpline (MSpline p ks cs) = bSpline ks cs' where n = p + 1; n' = fromIntegral n- cs' = [ (n' / (t1 - t0)) *^ c - | c <- cs- | (t0, t1) <- spans n (knots ks)- ]+ cs' = V.zipWith f cs (spans n (V.fromList (knots ks)))+ f c (t0, t1) = ((n' / (t1 - t0)) *^ c) instance Spline MSpline v => ControlPoints MSpline v where controlPoints = mSplineControlPoints@@ -69,7 +67,5 @@ where n = splineDegree spline + 1; n' = fromIntegral n ks = knotVector spline- cs = [ ((t1 - t0) / n') *^ c- | c <- controlPoints spline- | (t0, t1) <- spans n (knots ks)- ]+ cs = V.zipWith f (controlPoints spline) (spans n (V.fromList (knots ks)))+ f c (t0, t1) = ((t1 - t0) / n') *^ c