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
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
basisPolynomials :: (Fractional a, Ord a) => Knots a -> [[[Poly a]]]
basisPolynomials kts
| isEmpty kts = []
| otherwise = [basisPolynomialsAt kts kt | kt <- init (distinctKnots kts)]
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
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)