# alang9/splines forked from mokus0/splines

Faster version of deBoor. It assumes that the splines are `natural', …

`…i.e. they take the value 0 before and after the last spline.`
1 parent 15e96f5 commit bb7fe7f2f10fbd29ef08812c6f6381b0a67b4c0b committed Feb 6, 2012
Showing with 30 additions and 8 deletions.
1. +30 −8 src/Math/Spline/BSpline/Internal.hs
 @@ -5,6 +5,7 @@ module Math.Spline.BSpline.Internal import Math.Spline.Knots import Data.Monoid +import Data.Ratio import Data.Vector as V import Data.VectorSpace import Prelude as P @@ -40,29 +41,50 @@ insertKnot spline x = spline -- extend [1..5] -> [1,1,2,3,4,5,5] 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) +deBoor spline x = go uLo . vtake (deg + 1) . vdrop (l - deg) \$ + controlPoints spline where + l = maybe (-1) pred \$ V.findIndex (> x) us + deg = degree spline + zero = fromRational (0 % 1) + us = knotsVector (knotVector spline) - + uLo = stake (deg + 1) \$ sdrop (l - deg) \$ us -- 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 = sdrop (l + 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 us ds + go us ds | V.null ds = [] | otherwise = ds : go uLo ds' where uLo = V.tail us ds' = V.zipWith4 (interp x) uLo uHi ds (V.tail ds) + -- Try to take n, but if there's not enough, pan the rest with 0s + vtake n xs + | n <= V.length xs = V.take n xs + | otherwise = xs V.++ V.replicate (n - V.length xs) zeroV + stake n xs + | n <= V.length xs = V.take n xs + | otherwise = xs V.++ V.replicate (n - V.length xs) 0 + + -- Try to drop n, but if n is negative, pan the beginning with 0s + vdrop n xs + | n >= 0 = V.drop n xs + | otherwise = V.replicate (-n) zeroV V.++ xs + sdrop n xs + | n >= 0 = V.drop n xs + | otherwise = V.replicate (-n) 0 V.++ xs + interp x x0 x1 y0 y1 | x < x0 = y0 | x >= x1 = y1