diff --git a/benchmark/DeBoor.hs b/benchmark/DeBoor.hs new file mode 100644 --- /dev/null +++ b/benchmark/DeBoor.hs @@ -0,0 +1,63 @@ +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE FlexibleContexts #-} + +import Criterion.Main + +import Data.List (find) +import qualified Data.Vector.Generic.Safe as V +import qualified Data.Vector.Safe as BV +import qualified Data.Vector.Unboxed.Safe as UV +import Math.Spline.BSpline +import Math.Spline.BSpline.Reference +import Math.Spline.Knots +import Math.Polynomial + +import Debug.Trace +import Control.Monad + +kts = mkKnots $ [0,0,0] ++ [0..10] ++ [11,11,11] +ctPts = map sin [0..12] + +unboxedSpline :: BSpline UV.Vector Double +unboxedSpline = bSpline kts (V.fromList ctPts) + +boxedSpline :: BSpline BV.Vector Double +boxedSpline = bSpline kts (V.fromList ctPts) + + +intervalPoly :: [(Double, Double, Poly Double)] +intervalPoly = map f $ zip3 dkts (tail dkts) (basisPolynomials kts) + where + dkts = distinctKnots kts + f (begin, end, k) = (begin, end, ) . sumPolys $ zipWith scalePoly ctPts (k !! 3) + +applyDeBoor :: V.Vector v Double => BSpline v Double -> Double -> Double +applyDeBoor s = evalBSpline s + +{-# SPECIALIZE applyNaturalDeBoor :: BSpline BV.Vector Double -> Double -> Double #-} +{-# SPECIALIZE applyNaturalDeBoor :: BSpline UV.Vector Double -> Double -> Double #-} +applyNaturalDeBoor :: V.Vector v Double => BSpline v Double -> Double -> Double +applyNaturalDeBoor s = evalNaturalBSpline s + +applyPoly :: Double -> Double +applyPoly x = maybe 0 (\(_,_,p) -> evalPoly p x) $ find (\(b,e,_) -> x >= b && x < e) intervalPoly + +applyAndSum :: (Double -> Double) -> [Double] -> Double +applyAndSum f = sum . map f + +main = defaultMain + [ bgroup "Boxed" + [ bench "deBoor 1000" $ whnf (applyAndSum (applyDeBoor boxedSpline)) [0,0.01..10] + , bench "deBoor 10000" $ whnf (applyAndSum (applyDeBoor boxedSpline)) [0,0.001..10] + , bench "natural 1000" $ whnf (applyAndSum (applyNaturalDeBoor boxedSpline)) [0,0.01..10] + , bench "natural 10000" $ whnf (applyAndSum (applyNaturalDeBoor boxedSpline)) [0,0.001..10] + ] + , bgroup "Unboxed" + [ bench "deBoor 1000" $ whnf (applyAndSum (applyDeBoor unboxedSpline)) [0,0.01..10] + , bench "deBoor 10000" $ whnf (applyAndSum (applyDeBoor unboxedSpline)) [0,0.001..10] + , bench "natural 1000" $ whnf (applyAndSum (applyNaturalDeBoor unboxedSpline)) [0,0.01..10] + , bench "natural 10000" $ whnf (applyAndSum (applyNaturalDeBoor unboxedSpline)) [0,0.001..10] + ] + , bench "poly 1000" $ whnf (applyAndSum applyPoly) [0,0.01..10] + , bench "poly 10000" $ whnf (applyAndSum applyPoly) [0,0.001..10] + ] diff --git a/splines.cabal b/splines.cabal --- a/splines.cabal +++ b/splines.cabal @@ -1,8 +1,8 @@ name: splines -version: 0.3 +version: 0.5.0.1 stability: provisional -cabal-version: >= 1.6 +cabal-version: >= 1.9.2 build-type: Simple author: James Cook @@ -24,11 +24,13 @@ Library hs-source-dirs: src + ghc-options: -Wall exposed-modules: Math.Spline Math.Spline.BezierCurve Math.Spline.BSpline Math.Spline.BSpline.Reference Math.Spline.Class + Math.Spline.Hermite Math.Spline.ISpline Math.Spline.Knots Math.Spline.MSpline @@ -37,5 +39,31 @@ build-depends: base >= 3 && < 5, containers, polynomial, + vector >= 0.8, + vector-space + +Test-Suite splines-test + type: exitcode-stdio-1.0 + hs-source-dirs: test + main-is: Main.hs + + build-depends: base >= 3 && <5, + containers, + polynomial, + splines, + test-framework, + test-framework-quickcheck2, + QuickCheck >= 2, vector, vector-space + +Benchmark splines-bench + type: exitcode-stdio-1.0 + hs-source-dirs: benchmark + main-is: DeBoor.hs + + build-depends: base >= 3 && < 5, + criterion, + polynomial, + splines, + vector diff --git a/src/Math/NURBS.hs b/src/Math/NURBS.hs --- a/src/Math/NURBS.hs +++ b/src/Math/NURBS.hs @@ -8,13 +8,13 @@ ) where import qualified Data.Vector as V -import Data.VectorSpace +import Data.VectorSpace hiding (project) import Math.Spline.Class (Spline, toBSpline) import Math.Spline.BSpline.Internal import Math.Spline.BSpline import Math.Spline.Knots -newtype NURBS v = NURBS (BSpline (Scalar v, v)) +newtype NURBS v = NURBS (BSpline V.Vector (Scalar v, v)) deriving instance (Eq v, Eq (Scalar v), Eq (Scalar (Scalar v))) => Eq (NURBS v) deriving instance (Ord v, Ord (Scalar v), Ord (Scalar (Scalar v))) => Ord (NURBS v) @@ -34,16 +34,18 @@ -- |Constructs the homogeneous-coordinates B-spline that corresponds to this -- NURBS curve +nurbsAsSpline :: VectorSpace v => NURBS v -> BSpline V.Vector (Scalar v, v) nurbsAsSpline (NURBS spline) = spline { controlPoints = V.map homogenize (controlPoints spline) } where - homogenize (w,v) = (w, w *^ v) + homogenize (w,v) = (w, v ^* w) -- |Constructs the NURBS curve corresponding to a homogeneous-coordinates B-spline +splineAsNURBS :: (VectorSpace v, Fractional (Scalar v)) => BSpline V.Vector (Scalar v, v) -> NURBS v splineAsNURBS spline = NURBS spline { controlPoints = V.map unHomogenize (controlPoints spline) } where - unHomogenize (w,v) = (w, recip w *^ v) + unHomogenize (w,v) = (w, v ^/ w) evalNURBS @@ -51,7 +53,7 @@ VectorSpace w, Scalar w ~ w, Fractional w, Ord w) => NURBS v -> w -> v -evalNURBS nurbs = project . evalBSpline (nurbsAsSpline nurbs) +evalNURBS f = project . evalBSpline (nurbsAsSpline f) where project (w,v) = recip w *^ v @@ -75,6 +77,6 @@ VectorSpace w, Scalar w ~ w, Ord w, Fractional w) => NURBS v -> Scalar v -> Maybe (NURBS v, NURBS v) -splitNURBS nurbs t = do - (s0, s1) <- splitBSpline (nurbsAsSpline nurbs) t +splitNURBS f t = do + (s0, s1) <- splitBSpline (nurbsAsSpline f) t return (splineAsNURBS s0, splineAsNURBS s1) diff --git a/src/Math/Spline.hs b/src/Math/Spline.hs --- a/src/Math/Spline.hs +++ b/src/Math/Spline.hs @@ -7,6 +7,7 @@ , BSpline, bSpline , MSpline, mSpline, toMSpline , ISpline, iSpline, toISpline + , CSpline, cSpline ) where import Math.Spline.Class @@ -15,3 +16,4 @@ import Math.Spline.BSpline import Math.Spline.MSpline import Math.Spline.ISpline +import Math.Spline.Hermite 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 @@ -1,11 +1,13 @@ -{-# LANGUAGE MultiParamTypeClasses, StandaloneDeriving, FlexibleContexts, UndecidableInstances, TypeFamilies, ParallelListComp #-} +{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, TypeFamilies, ParallelListComp #-} module Math.Spline.BSpline ( BSpline , bSpline , evalBSpline + , evalNaturalBSpline , insertKnot , splitBSpline , differentiateBSpline, integrateBSpline + , deBoor ) where import Math.Spline.Knots @@ -13,20 +15,19 @@ import Data.Maybe (fromMaybe) import Data.VectorSpace -import qualified Data.Vector as V +import qualified Data.Vector.Generic 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) -> 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) +bSpline :: V.Vector v a => Knots (Scalar a) -> v a -> BSpline v a +bSpline kts cps = fromMaybe + (error "bSpline: too few knots") + (maybeSpline kts cps) -maybeSpline :: Knots (Scalar a) -> V.Vector a -> Maybe (BSpline a) +-- not exported; precondition: n > 0 +maybeSpline :: V.Vector v a => Knots (Scalar a) -> v a -> Maybe (BSpline v a) maybeSpline kts cps | n > m = Nothing | otherwise = Just (Spline (m - n) kts cps) @@ -34,23 +35,18 @@ n = V.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 + :: (VectorSpace a, Fractional (Scalar a), Ord (Scalar a), V.Vector v a + , V.Vector v (Scalar a)) => BSpline v a -> BSpline v a differentiateBSpline spline - | numKnots ks < 2 = spline - | numKnots ks == 2 = bSpline ks (V.singleton zeroV) - | otherwise = bSpline ks' ds' + | V.null ds = error "differentiateBSpline: no control points" + | m < 1 = spline + | p == 0 = bSpline ks (V.replicate n zeroV) + | otherwise = bSpline ks' ds' where + n = V.length ds + m = numKnots ks - 1 + ks' = mkKnots . init . tail $ ts ds' = V.zipWith (*^) (V.tail cs) (V.zipWith (^-^) (V.tail ds) ds) @@ -58,10 +54,12 @@ ds = controlPoints spline p = degree spline - cs = V.fromList [fromIntegral p / (t1 - t0) | (t0,t1) <- spans p ts] + cs = V.fromList [ if t1 /= t0 then fromIntegral p / (t1 - t0) else 0 | (t0,t1) <- spans p ts] integrateBSpline - :: (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => BSpline v -> BSpline v + :: (VectorSpace a, Fractional (Scalar a), Ord (Scalar a), V.Vector v a + , V.Vector v (Scalar a)) => + BSpline v a -> BSpline v a integrateBSpline spline = bSpline (mkKnots ts') (V.scanl (^+^) zeroV ds') where ds' = V.zipWith (*^) cs (controlPoints spline) @@ -70,26 +68,31 @@ p = degree spline + 1 cs = V.fromList [(t1 - t0) / fromIntegral p | (t0,t1) <- spans p ts] +spans :: Int -> [a] -> [(a,a)] spans n xs = zip xs (drop n xs) +-- |Split a B-spline at the specified point (which must be inside the spline's domain), +-- returning two disjoint splines, the sum of which is equal to the original. The domain +-- of the first will be below the split point and the domain of the second will be above. 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) + :: ( VectorSpace a, Ord (Scalar a), Fractional (Scalar a), V.Vector v a + , V.Vector v (Scalar a)) => + BSpline v a -> Scalar a -> Maybe (BSpline v a, BSpline v a) +splitBSpline spline@(Spline p kv _) t + | inDomain = Just (Spline p us0 ds0, Spline p us1 ds1) | otherwise = Nothing where inDomain = case knotDomain kv p of Nothing -> False - Just (t0, t1) -> t >= t0 || t <= t1 + Just (t0, t1) -> t >= t0 && t <= t1 - us = knots kv + (lt, _, gt) = splitFind t kv dss = deBoor spline t - us0 = takeWhile ( Eq (BSpline v a) +deriving instance (Ord (Scalar a), Ord (v a)) => Ord (BSpline v a) +instance (Show (Scalar a), Show a, Show (v a)) => Show (BSpline v a) where + showsPrec p (Spline _ kts cps) = showParen (p>10) + ( showString "bSpline " + . showsPrec 11 kts + . showChar ' ' + . showsPrec 11 cps + ) + +mapControlPoints :: (Scalar a ~ Scalar b, V.Vector v a, V.Vector v b) => (a -> b) -> BSpline v a -> BSpline v b +{-# SPECIALIZE mapControlPoints :: (Scalar a ~ Scalar b) => (a -> b) + -> BSpline BV.Vector a -> BSpline BV.Vector b #-} mapControlPoints f spline = spline { controlPoints = V.map f (controlPoints spline) , knotVector = knotVector spline } -evalBSpline spline = V.head . P.last . deBoor spline +-- |Evaluate a B-spline at the given point. This uses a slightly modified version of +-- de Boor's algorithm which is only strictly correct inside the domain of the spline. +-- Unlike the standard algorithm, the basis functions always sum to 1, even outside the +-- domain of the spline. This is mainly useful for \"clamped\" splines - the values at +-- or outside the endpoints will always be the value of the nearest control point. +-- +-- For a standard implementation of de Boor's algorithm, see 'evalNaturalBSpline'. +-- For a (much slower) strictly mathematically correct evaluation, see 'evalReferenceBSpline'. +evalBSpline :: ( VectorSpace a, Fractional (Scalar a), Ord (Scalar a) + , V.Vector v a, V.Vector v (Scalar a)) + => BSpline v a -> Scalar a -> a +{-# SPECIALIZE evalBSpline :: ( VectorSpace a, Fractional (Scalar a) + , Ord (Scalar a)) => BSpline BV.Vector a -> Scalar a -> a#-} +evalBSpline spline + | V.null (controlPoints spline) = zeroV + | otherwise = V.head . P.last . deBoor spline --- |Insert one knot into a 'BSpline' +-- | Evaluate a B-spline at the given point. This uses de Boor's algorithm, which is +-- only strictly correct inside the domain of the spline. +-- +-- For a (much slower) strictly mathematically correct evaluation, see 'evalReferenceBSpline'. +{-# INLINE evalNaturalBSpline #-} +{-# SPECIALIZE evalNaturalBSpline :: (Fractional (Scalar a), Ord (Scalar a), VectorSpace a) + => BSpline BV.Vector a -> Scalar a -> a #-} +{-# SPECIALIZE evalNaturalBSpline :: (Fractional (Scalar a), Ord (Scalar a), VectorSpace a, UV.Unbox a) + => BSpline UV.Vector a -> Scalar a -> a #-} +evalNaturalBSpline :: (Fractional (Scalar a), Ord (Scalar a), VectorSpace a, V.Vector v a) + => BSpline v a -> Scalar a -> a +evalNaturalBSpline spline x = runST $ do + let Spline p kvec cps = slice spline x + kts = V.tail (knotsVector kvec) + s = p - 1 - V.length (V.takeWhile (x==) kts) + + ds <- V.thaw cps + + sequence_ + [ do + -- No need to check whether u0 < x < u1: + -- x > u0 is guaranteed by choice of 's' + -- x < u1 is guaranteed by 'slice' + let !u0 = kts V.! (j + i) + !u1 = kts V.! (j + p) + !du = u1 - u0 + !a = if du <= 0 then 1 else (x - u0) / du + + d0 <- MV.read ds j + d1 <- MV.read ds (j + 1) + MV.write ds j (lerp d0 d1 a) + | i <- [0 .. s] + , j <- [0 .. s - i] + ] + + MV.read ds 0 + +-- |Insert one knot into a 'BSpline' without changing the spline's shape. insertKnot - :: (VectorSpace a, Ord (Scalar a), Fractional (Scalar a)) => - BSpline a -> Scalar a -> BSpline a + :: (VectorSpace a, Ord (Scalar a), Fractional (Scalar a), V.Vector v a, V.Vector v (Scalar a)) => + BSpline v a -> Scalar a -> BSpline v a +{-# SPECIALIZE insertKnot :: (VectorSpace a, Ord (Scalar a), Fractional (Scalar a)) => + BSpline BV.Vector a -> Scalar a -> BSpline BV.Vector a #-} insertKnot spline x = spline { knotVector = knotVector spline `mappend` knot x , controlPoints = V.zipWith4 (interp x) us (V.drop p us) ds (V.tail ds) } where - us = knotsVector (knotVector spline) + us = V.convert $ 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] +{-# INLINE extend #-} +extend :: V.Vector v t => v t -> v t extend vec | V.null vec = V.empty - | otherwise = V.cons (V.head vec) (V.snoc vec (V.last vec)) + | otherwise = V.cons (V.head vec) (V.snoc vec (V.last vec)) -deBoor spline x = go us (controlPoints spline) +{-# SPECIALIZE deBoor :: (Fractional (Scalar a), Ord (Scalar a), VectorSpace a) + => BSpline BV.Vector a -> Scalar a -> [BV.Vector a] #-} +{-# SPECIALIZE deBoor :: (Fractional (Scalar a), Ord (Scalar a), VectorSpace a, + UV.Unbox a, UV.Unbox (Scalar a)) + => BSpline UV.Vector a -> Scalar a -> [UV.Vector a] #-} + +-- | The table from de Boor's algorithm, calculated for the entire spline. If that is not necessary +-- (for example, if you are only evaluating the spline), then use 'slice' on the spline first. +-- 'splitBSpline' currently uses the whole table. It is probably not necessary there, but it +-- greatly simplifies the definition and makes the similarity to splitting Bezier curves very obvious. +deBoor :: (Fractional (Scalar a), Ord (Scalar a), VectorSpace a, V.Vector v a, V.Vector v (Scalar a)) + => BSpline v a -> Scalar a -> [v a] +deBoor spline x = go us0 (controlPoints spline) where - us = knotsVector (knotVector spline) - + us0 = V.convert $ 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 = V.drop (degree spline + 1) us - - -- On each pass, the lower endpoints of the - -- interpolation intervals advance and the new + uHi = V.drop (degree spline + 1) us0 + + -- 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 us ds + go us ds | V.null ds = [] | otherwise = ds : go uLo ds' where @@ -63,10 +158,64 @@ ds' = V.zipWith4 (interp x) uLo uHi ds (V.tail ds) +interp :: (Fractional (Scalar v), Ord (Scalar v), VectorSpace v) + => Scalar v -> Scalar v -> Scalar v -> v -> v -> v interp x x0 x1 y0 y1 | x < x0 = y0 | x >= x1 = y1 | otherwise = lerp y0 y1 a where a = (x - x0) / (x1 - x0) + +-- "slice" a spline to contain only those knots and control points that +-- actually influence the value at 'x'. +-- +-- It should be true for any valid BSpline that: +-- degree (slice f x) == degree f +-- slice (slice f x) x == slice f x +-- {x in domain of f} => {x in domain of slice f x} +-- {x in domain of f} => evalBSpline (slice f x) x == evalBSpline f x +{-# INLINE slice #-} +slice :: (Num (Scalar a), Ord (Scalar a), AdditiveGroup a, V.Vector v a) + => BSpline v a -> Scalar a -> BSpline v a +slice spline x = spline + { knotVector = stakeKnots (n + n) . sdropKnots (l - n) $ knotVector spline + , controlPoints = vtake n . vdrop (l - n) $ controlPoints spline + } + where + l = maybe 0 id $ V.findIndex (> x) us + n = degree spline + 1 + + us = knotsVector (knotVector spline) + +-- Try to take n, but if there's not enough, pad the rest with 0s +vtake :: (V.Vector v t, AdditiveGroup t) => Int -> v t -> v t +{-# SPECIALIZE vtake :: AdditiveGroup t => Int -> BV.Vector t -> BV.Vector t #-} +vtake n xs + | n <= V.length xs = V.take n xs + | otherwise = xs V.++ V.replicate (n - V.length xs) zeroV + +-- Try to drop n, but if n is negative, pad the beginning with 0s +vdrop :: (V.Vector v t, AdditiveGroup t) => Int -> v t -> v t +{-# SPECIALIZE vdrop :: AdditiveGroup t => Int -> BV.Vector t -> BV.Vector t #-} +vdrop n xs + | n >= 0 = V.drop n xs + | otherwise = V.replicate (-n) zeroV V.++ xs + +-- Try to take n knots, but if there aren't enough, increase the multiplicity of the last knot +stakeKnots :: (Num k, Ord k) => Int -> Knots k -> Knots k +stakeKnots n kts + | n <= nKts = takeKnots n kts + | otherwise = case maxKnot kts of + Nothing -> multipleKnot 0 (n - nKts) + Just (k, m) -> setKnotMultiplicity k (m + n - nKts) kts + where nKts = numKnots kts + +-- Try to drop n knots, but if n is negative, increase the multiplicity of the first knot by @abs n@ +sdropKnots :: (Num k, Ord k) => Int -> Knots k -> Knots k +sdropKnots n kts + | n >= 0 = dropKnots n kts + | otherwise = case minKnot kts of + Nothing -> multipleKnot 0 (-n) + Just (k, m) -> setKnotMultiplicity k (m - n) kts diff --git a/src/Math/Spline/BSpline/Reference.hs b/src/Math/Spline/BSpline/Reference.hs --- a/src/Math/Spline/BSpline/Reference.hs +++ b/src/Math/Spline/BSpline/Reference.hs @@ -1,4 +1,6 @@ +{-# LANGUAGE TypeFamilies #-} {-# LANGUAGE ParallelListComp #-} +{-# LANGUAGE FlexibleContexts #-} -- |Reference implementation of B-Splines; very inefficient but \"obviously\" -- correct. module Math.Spline.BSpline.Reference @@ -6,15 +8,40 @@ , basisFunctions , basisPolynomials , basisPolynomialsAt + , evalReferenceBSpline + , fitPolyToBSplineAt ) where +import qualified Data.Vector.Generic as V +import Data.VectorSpace (VectorSpace, Scalar, (^*), sumV) import Math.Spline.Knots +import Math.Spline.BSpline.Internal import Math.Polynomial (Poly) import qualified Math.Polynomial as Poly +-- | This is a fairly slow function which computes the value of a B-spline at a given point, +-- using the mathematical definition of B-splines. It is mainly for testing purposes, as a +-- reference against which the other evaluation functions are checked. +evalReferenceBSpline :: (VectorSpace a, Fractional (Scalar a), Ord (Scalar a), V.Vector v a) + => BSpline v a -> Scalar a -> a +evalReferenceBSpline (Spline deg kts cps) x = + sumV (zipWith (^*) (V.toList cps) (bases kts x !! deg)) + +-- | This is a fairly slow function which computes one polynomial segment of a B-spline (the +-- one containing the given point), using the mathematical definition of B-splines. It is +-- mainly for testing purposes, as a reference against which the other evaluation functions +-- are checked. +fitPolyToBSplineAt :: (Fractional a, Ord a, Scalar a ~ a, V.Vector v a) + => BSpline v a -> a -> Poly a +fitPolyToBSplineAt (Spline deg kts cps) x = + Poly.sumPolys (zipWith Poly.scalePoly (V.toList cps) (basisPolynomialsAt kts x !! deg)) + +ind :: Num a => Bool -> a ind True = 1 ind False = 0 +-- | The values of all the B-spline basis functions for the given knot vector at the given +-- point, ordered by degree; \"b_{i,j}(x)\" is @bases kts x !! i !! j@. bases :: (Fractional a, Ord a) => Knots a -> a -> [[a]] bases kts x = coxDeBoor interp initial kts where @@ -26,8 +53,8 @@ = (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 +-- | All the B-spline basis functions for the given knot vector at the given +-- point, ordered by degree; \"b_{i,j}\" is @basisFunctions kts x !! i !! j@. basisFunctions :: (Fractional a, Ord a) => Knots a -> [[a -> a]] basisFunctions kts = coxDeBoor interp initial kts where @@ -39,13 +66,15 @@ = (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. +-- | All the B-spline basis polynomials for the given knot vector, ordered first +-- by knot span and then by degree. 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. +-- | All the B-spline basis polynomials for the given knot vector at the given +-- point, ordered by degree; \"b_{i,j}\" is @basisPolynomialsAt kts x !! i !! j@. basisPolynomialsAt :: (Fractional a, Ord a) => Knots a -> a -> [[Poly a]] basisPolynomialsAt kts x = coxDeBoor interp initial kts where @@ -57,22 +86,23 @@ | (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 + = (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 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 + 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 +-- | 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 :: Num a => (a -> a -> b -> a -> a -> b -> b) -> [b] -> Knots a -> [[b]] coxDeBoor interp initial kts = table where ts = knots kts @@ -90,4 +120,6 @@ spans = spansWith (,) spanDiffs :: Num a => Int -> [a] -> [a] spanDiffs = spansWith subtract + +spansWith :: (a -> a -> b) -> Int -> [a] -> [b] 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 @@ -12,7 +12,7 @@ import qualified Data.Vector as V import Data.VectorSpace --- |A BezierCurve curve on @0 <= x <= 1@. +-- |A Bezier curve on @0 <= x <= 1@. data BezierCurve t = BezierCurve !Int !(V.Vector t) deriving (Eq, Ord) -- |Construct a Bezier curve from a list of control points. The degree 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 @@ -1,4 +1,5 @@ {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, FlexibleContexts, UndecidableInstances #-} +{-# LANGUAGE OverlappingInstances, IncoherentInstances #-} module Math.Spline.Class where import Control.Applicative @@ -6,6 +7,7 @@ import qualified Math.Spline.BSpline.Internal as BSpline import qualified Data.Vector as V +import qualified Data.Vector.Generic as G import Data.VectorSpace -- |A spline is a piecewise polynomial vector-valued function. The necessary @@ -30,16 +32,27 @@ knotVector :: s v -> Knots (Scalar v) knotVector = knotVector . toBSpline - toBSpline :: s v -> BSpline.BSpline v + toBSpline :: s v -> BSpline.BSpline V.Vector v +-- TODO: this class should probably go away. all it really does is overload something that doesn't really have any implementation-independent semantics (or does it?). class Spline s v => ControlPoints s v where controlPoints :: s v -> V.Vector v -instance (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline BSpline.BSpline v where - evalSpline spline = V.head . last . BSpline.deBoor spline +instance (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline (BSpline.BSpline V.Vector) v where + evalSpline = BSpline.evalBSpline splineDegree = BSpline.degree knotVector = BSpline.knotVector toBSpline = id -instance Spline BSpline.BSpline v => ControlPoints BSpline.BSpline v where +instance ( VectorSpace a, Fractional (Scalar a), Ord (Scalar a), G.Vector v a + , G.Vector v (Scalar a)) => Spline (BSpline.BSpline v) a where + evalSpline = BSpline.evalBSpline + splineDegree = BSpline.degree + knotVector = BSpline.knotVector + toBSpline (BSpline.Spline deg ks ctp) = BSpline.Spline deg ks (G.convert $ ctp) + +instance Spline (BSpline.BSpline V.Vector) a => ControlPoints (BSpline.BSpline V.Vector) a where controlPoints = BSpline.controlPoints + +instance (Spline (BSpline.BSpline v) a, G.Vector v a) => ControlPoints (BSpline.BSpline v) a where + controlPoints = V.convert . BSpline.controlPoints diff --git a/src/Math/Spline/Hermite.hs b/src/Math/Spline/Hermite.hs new file mode 100644 --- /dev/null +++ b/src/Math/Spline/Hermite.hs @@ -0,0 +1,89 @@ +{-# LANGUAGE ParallelListComp #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +module Math.Spline.Hermite + ( CSpline, cSpline + , evalSpline + ) where + +import Data.List +import Data.Ord +import Math.Polynomial hiding (x) +import Math.Spline.BSpline +import Math.Spline.Class +import Math.Spline.Knots +import qualified Data.Vector as V +import Data.VectorSpace + +-- | Cubic Hermite splines. These are cubic splines defined by a +-- sequence of control points and derivatives at those points. +newtype CSpline a = CSpline [(Scalar a,a,a)] + +-- | Cubic splines specified by a list of control points, +-- where each control point is given by a triple of parameter value, +-- position of the spline at that parameter value, +-- and derivative of the spline at that parameter value. +cSpline :: Ord (Scalar a) => [(Scalar a,a,a)] -> CSpline a +cSpline = CSpline . sortBy (comparing fst3) + where fst3 (a,_,_) = a + +h00, h10, h01, h11 :: (Num a, Eq a) => Poly a +h00 = poly BE [ 2,-3, 0, 1] +h10 = poly BE [ 1,-2, 1, 0] +h01 = poly BE [-2, 3, 0, 0] +h11 = poly BE [ 1,-1, 0, 0] + +evalHermite + :: (Eq (Scalar v), Num (Scalar v), VectorSpace v) => + v -> v -> v -> v -> Scalar v -> v +evalHermite y0 m0 y1 m1 x = sumV + [ evalPoly h x *^ p + | p <- [ y0, m0, y1, m1] + | h <- [h00, h10, h01, h11] + ] + +evalCSpline + :: (VectorSpace v, Ord (Scalar v), Fractional (Scalar v)) => + CSpline v -> Scalar v -> v +evalCSpline (CSpline cps) = loop cps + where + loop [] _ = zeroV + loop [(_,y,_)] _ = y + loop ((x0,y0,m0):rest@((x1,y1,m1):_)) x + | x <= x0 = y0 + | x < x1 = let dx = x1 - x0 + in evalHermite y0 (m0 ^* dx) y1 (m1 ^* dx) ((x - x0) / dx) + | otherwise = loop rest x + +instance (VectorSpace a, Fractional (Scalar a), Ord (Scalar a)) => Spline CSpline a where + splineDegree _ = 3 + + splineDomain (CSpline []) = Nothing + splineDomain (CSpline cps) = Just (head xs, last xs) where (xs, _, _) = unzip3 cps + + evalSpline = evalCSpline + + -- TODO: check. Also work out a more compact translation, taking advantage of the + -- known continuity on the interior knots. It should be possible to work out an + -- equivalent b-spline with only 'n+4' knots. If that translation isn't ill- + -- conditioned, it might be a good thing to implement. + toBSpline (CSpline cSpl) = bSpline kts (V.fromList cps) + where + kts = fromList [(x,4) | x <- xs] + + (xs, _, _) = unzip3 cSpl + cps = concat + [ [ y0 + , y0 ^+^ dy0 ^* dx3 + , y1 ^-^ dy1 ^* dx3 + , y1 + ] + | ((x0, y0, dy0), (x1, y1, dy1)) <- spans 1 cSpl + , let dx3 = (x1 - x0) / 3 + ] + +spans :: Int -> [a] -> [(a,a)] +spans n xs = zip xs (drop n xs) + 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 @@ -62,6 +62,8 @@ toISpline :: (Spline s v, Eq v) => s v -> ISpline v toISpline = fromBSpline . toBSpline +fromBSpline :: (Eq v, VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) + => BSpline V.Vector v -> ISpline v fromBSpline spline | V.head ds == zeroV && numKnots ks >= 2 = iSpline (mkKnots (init (tail ts))) (V.tail 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 @@ -10,7 +10,9 @@ , toList, numDistinctKnots, lookupDistinctKnot , knots, knotsVector - , distinctKnots, distinctKnotsVector + , distinctKnots, multiplicities + , distinctKnotsVector, multiplicitiesVector + , distinctKnotsSet , toMap , fromMap @@ -25,6 +27,8 @@ , maxMultiplicity , knotMultiplicity, setKnotMultiplicity + , splitFind + , fromAscList, fromDistinctAscList , valid @@ -34,27 +38,26 @@ , knotDomain , uniform + + , minKnot + , maxKnot ) where import Prelude hiding (sum, maximum) +import Control.Arrow ((***)) import Control.Monad (guard) -import Data.Foldable (Foldable(foldMap), sum, maximum) +import Data.Foldable (Foldable(foldMap), maximum) +import Data.List (sortBy, sort, unfoldr) import qualified Data.Map as M import Data.Monoid (Monoid(..)) -import Data.Maybe (fromMaybe) -import qualified Data.Set as S (Set) +import Data.Ord +import qualified Data.Set as S (Set, fromAscList) 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) +newtype Knots a = Knots (V.Vector a) 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) - ) showsPrec p ks = showParen (p > 10) ( showString "mkKnots " . showsPrec 11 (knots ks) @@ -62,20 +65,19 @@ instance (Ord a) => Monoid (Knots a) where mempty = empty - mappend (Knots n1 v1) (Knots n2 v2) = - Knots (n1 + n2) (M.filter (/=0) (M.unionWith (+) v1 v2)) + mappend (Knots v1) (Knots v2) = + Knots . V.fromList . sort . V.toList $ v1 V.++ v2 instance Foldable Knots where - foldMap f = foldMap f . knots + foldMap f = foldMap f . knotsVector -- |An empty knot vector empty :: Knots a -empty = Knots 0 M.empty +empty = Knots V.empty isEmpty :: Knots a -> Bool -isEmpty (Knots 0 _) = True -isEmpty _ = False +isEmpty (Knots v) = V.null v -- |Create a knot vector consisting of one knot. knot :: Ord a => a -> Knots a @@ -83,79 +85,74 @@ -- |Create a knot vector consisting of one knot with the specified multiplicity. multipleKnot :: Ord a => a -> Int -> Knots a -multipleKnot k n - | n <= 0 = Knots 0 (M.empty) - | otherwise = Knots n (M.singleton k n) +multipleKnot k n = Knots $ V.replicate n k -- |Create a knot vector consisting of all the knots in a list. mkKnots :: (Ord a) => [a] -> Knots a -mkKnots ks = fromList (map (\k -> (k,1)) ks) +mkKnots = Knots . V.fromList . sort -- |Create a knot vector consisting of all the knots and corresponding -- multiplicities in a list. fromList :: (Ord k) => [(k, Int)] -> Knots k -fromList ks = Knots (sum kMap) kMap - where kMap = M.fromListWith (+) (filter ((>0).snd) ks) +fromList ks = Knots v + where v = V.concat . map (\(k, mult) -> V.replicate mult k) . + sortBy (comparing fst) $ 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) +fromAscList ks = Knots v + where v = V.concat . map (\(k, mult) -> V.replicate mult k) + $ 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) +fromDistinctAscList :: Eq k => [(k, Int)] -> Knots k +fromDistinctAscList = fromAscList -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 +fromMap :: Eq k => M.Map k Int -> Knots k +fromMap = fromAscList . 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 = M.toList . toMap +toList :: Eq k => Knots k -> [(k, Int)] +toList = unfoldr $ \kts -> do + kt <- minKnot kts + return (kt, dropKnots (snd kt) kts) -toVector :: Knots k -> V.Vector (k, Int) -toVector = V.fromList . toList +toVector :: Eq k => Knots k -> V.Vector (k, Int) +toVector = V.unfoldr $ \kts -> do + kt <- minKnot kts + return (kt, dropKnots (snd kt) kts) -toMap :: Knots k -> M.Map k Int -toMap (Knots _ ks) = ks +toMap :: Ord k => Knots k -> M.Map k Int +toMap = M.fromAscListWith (+) . toList -- |Returns the number of knots (not necessarily distinct) in a knot vector. numKnots :: Knots t -> Int -numKnots (Knots n _) = n +numKnots (Knots v) = V.length v -- |Returns the number of distinct knots in a knot vector. -numDistinctKnots :: Knots t -> Int -numDistinctKnots (Knots _ ks) = M.size ks +numDistinctKnots :: Eq t => Knots t -> Int +numDistinctKnots = V.length . distinctKnotsVector -maxMultiplicity :: Knots t -> Int -maxMultiplicity (Knots 0 _) = 0 -maxMultiplicity (Knots _ ks) = maximum ks +maxMultiplicity :: Ord t => Knots t -> Int +maxMultiplicity kts@(Knots v) + | V.length v == 0 = 0 + | otherwise = maximum $ toMap kts 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 +lookupKnot k (Knots kts) = do + guard (0 <= k && k < V.length kts) + V.indexM kts k -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 +lookupDistinctKnot :: Eq a => Int -> Knots a -> Maybe a +lookupDistinctKnot k kts = lookupKnot k . Knots $ distinctKnotsVector kts -- |@splitLookup n kts@: Split a knot vector @kts@ into 3 parts @(pre, mbKt, post)@ -- such that: @@ -165,133 +162,84 @@ -- * 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 +splitLookup :: Int -> Knots a -> (Knots a, Maybe a, Knots a) +splitLookup k (Knots v) + | V.null gt = (Knots lt, Nothing, Knots V.empty) + | otherwise = (Knots lt, Just $ V.head gt, Knots $ V.tail gt) + where + (lt, gt) = V.splitAt k v 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 +dropKnots k (Knots v) = Knots $ V.drop k v 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 +takeKnots k (Knots v) = Knots $ V.take k v 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 +splitKnotsAt k (Knots v) = Knots *** Knots $ V.splitAt k v +-- |Count the number of knots less than the n'th distinct knot. +findDistinctKnot :: Eq a => Int -> Knots a -> Int +findDistinctKnot n = V.last . V.take (1 + max 0 n) . V.scanl (+) 0 . multiplicitiesVector -takeDistinctKnots :: Int -> Knots a -> Knots a -takeDistinctKnots k (Knots n ks) = Knots (sum kMap) kMap - where - kMap = M.fromDistinctAscList (take k (M.toAscList ks)) +takeDistinctKnots :: (Ord a) => Int -> Knots a -> Knots a +takeDistinctKnots k kts = takeKnots (findDistinctKnot k kts) kts -dropDistinctKnots :: Int -> Knots a -> Knots a -dropDistinctKnots k (Knots n ks) = Knots (sum kMap) kMap - where - kMap = M.fromDistinctAscList (drop k (M.toAscList ks)) +dropDistinctKnots :: (Ord a) => Int -> Knots a -> Knots a +dropDistinctKnots k kts = dropKnots (findDistinctKnot k kts) kts -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 +splitDistinctKnotsAt :: (Ord a, Eq a) => Int -> Knots a -> (Knots a, Knots a) +splitDistinctKnotsAt k kts = splitKnotsAt (findDistinctKnot k kts) kts -- |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] +knots = V.toList . knotsVector -- |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] +knotsVector (Knots v) = v -- |Returns a list of all distinct knots of a knot vector in ascending order -distinctKnots :: Knots t -> [t] -distinctKnots (Knots _ ks) = M.keys ks +distinctKnots :: Eq t => Knots t -> [t] +distinctKnots = map fst . toList +multiplicities :: Eq t => Knots t -> [Int] +multiplicities = map snd . toList + -- |Returns a vector of all distinct knots of a knot vector in ascending order -distinctKnotsVector :: Knots t -> V.Vector t -distinctKnotsVector = V.fromList . distinctKnots +distinctKnotsVector :: Eq t => Knots t -> V.Vector t +distinctKnotsVector = V.map fst . toVector +multiplicitiesVector :: Eq a => Knots a -> V.Vector Int +multiplicitiesVector = V.map snd . toVector + -- |Returns a 'S.Set' of all distinct knots of a knot vector -distinctKnotsSet :: Knots k -> S.Set k -distinctKnotsSet (Knots _ ks) = M.keysSet ks +distinctKnotsSet :: Eq k => Knots k -> S.Set k +distinctKnotsSet (Knots k) = S.fromAscList $ V.toList k -- |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) +knotMultiplicity k (Knots ks) = V.length $ V.elemIndices k ks -- |Returns a new knot vector with the given knot set to the specified -- multiplicity and all other knots unchanged. setKnotMultiplicity :: Ord k => k -> Int -> Knots k -> Knots k -setKnotMultiplicity k n (Knots m ks) - | n <= 0 = Knots (m - n') (M.delete k ks) - | otherwise = Knots (m + n - n') (M.insert k n ks) - where - n' = knotMultiplicity k (Knots m ks) +setKnotMultiplicity k n kts@(Knots v) + | n <= 0 = Knots (V.filter (/= k) v) + | otherwise = Knots $ V.concat [lt, V.replicate n k, gt] + where (Knots lt, _, Knots gt) = splitFind k kts +splitFind :: Ord k => k -> Knots k -> (Knots k, Knots k, Knots k) +splitFind k (Knots v) = (Knots lt, Knots eq, Knots gt) + where + (lt, xs) = V.span ( Knots k -> Bool -valid (Knots n ks) = and - [ M.valid ks - , n == sum ks - , all (>0) (M.elems ks) - ] +valid :: (Ord k, Num k) => Knots k -> Bool +valid (Knots v) + | V.null v = True + | otherwise = V.and $ V.zipWith (>=) (V.tail v) v -- |@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. @@ -324,7 +272,7 @@ -- 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 = knotSpan ks p (n-p-1) +knotDomain ks@(Knots v) p = knotSpan ks p (V.length v-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 @@ -336,3 +284,14 @@ n = nPts + deg - numKnots ends f i = (fromIntegral i * lo + fromIntegral (n - i) * hi) / fromIntegral n internal = mkKnots [f i | i <- [0..n]] + +{-# INLINE minKnot #-} +minKnot :: (Eq a) => Knots a -> Maybe (a, Int) +minKnot (Knots v) + | V.null v = Nothing + | otherwise = Just (kt, V.length (V.takeWhile (kt ==) v)) + where kt = V.head v + +{-# INLINE maxKnot #-} +maxKnot :: Eq a => Knots a -> Maybe (a, Int) +maxKnot (Knots v) = minKnot (Knots (V.reverse v)) 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 @@ -46,6 +46,7 @@ n = V.length cps m = numKnots kts - 1 +spans :: Int -> V.Vector a -> V.Vector (a,a) spans n xs = V.zip xs (V.drop n xs) instance (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline MSpline v where @@ -63,6 +64,8 @@ toMSpline :: Spline s v => s v -> MSpline v toMSpline = fromBSpline . toBSpline +fromBSpline :: (VectorSpace a, Fractional (Scalar a), Ord (Scalar a)) + => BSpline (V.Vector) a -> MSpline a fromBSpline spline = mSpline ks cs where n = splineDegree spline + 1; n' = fromIntegral n diff --git a/test/Main.hs b/test/Main.hs new file mode 100644 --- /dev/null +++ b/test/Main.hs @@ -0,0 +1,14 @@ +#!/usr/bin/env runhaskell +module Main where + +import Test.Framework (defaultMain, testGroup) + +import Tests.BSpline.Reference (referenceBSplineTests) +import Tests.BSpline (bSplineTests) +import Tests.Knots (knotsTests) + +main = defaultMain + [ testGroup "Math.Spline.BSpline.Reference" referenceBSplineTests + , testGroup "Math.Spline.BSpline" bSplineTests + , testGroup "Math.Spline.Knots" knotsTests + ]