diff --git a/splines.cabal b/splines.cabal --- a/splines.cabal +++ b/splines.cabal @@ -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 diff --git a/src/Math/NURBS.hs b/src/Math/NURBS.hs --- a/src/Math/NURBS.hs +++ b/src/Math/NURBS.hs @@ -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, diff --git a/src/Math/Spline.hs b/src/Math/Spline.hs --- a/src/Math/Spline.hs +++ b/src/Math/Spline.hs @@ -1,5 +1,5 @@ module Math.Spline - ( Spline(..) + ( Spline(..), ControlPoints(..) , Knots, mkKnots, knots diff --git a/src/Math/Spline/BSpline.hs b/src/Math/Spline/BSpline.hs --- a/src/Math/Spline/BSpline.hs +++ b/src/Math/Spline/BSpline.hs @@ -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 ( 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 diff --git a/src/Math/Spline/BSpline/Reference.hs b/src/Math/Spline/BSpline/Reference.hs new file mode 100644 --- /dev/null +++ b/src/Math/Spline/BSpline/Reference.hs @@ -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) diff --git a/src/Math/Spline/BezierCurve.hs b/src/Math/Spline/BezierCurve.hs --- a/src/Math/Spline/BezierCurve.hs +++ b/src/Math/Spline/BezierCurve.hs @@ -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 diff --git a/src/Math/Spline/Class.hs b/src/Math/Spline/Class.hs --- a/src/Math/Spline/Class.hs +++ b/src/Math/Spline/Class.hs @@ -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 diff --git a/src/Math/Spline/ISpline.hs b/src/Math/Spline/ISpline.hs --- a/src/Math/Spline/ISpline.hs +++ b/src/Math/Spline/ISpline.hs @@ -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) diff --git a/src/Math/Spline/Knots.hs b/src/Math/Spline/Knots.hs --- a/src/Math/Spline/Knots.hs +++ b/src/Math/Spline/Knots.hs @@ -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]] diff --git a/src/Math/Spline/MSpline.hs b/src/Math/Spline/MSpline.hs --- a/src/Math/Spline/MSpline.hs +++ b/src/Math/Spline/MSpline.hs @@ -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