public
Created

Derivative of B-spline in Python

  • Download Gist
splder.py
Python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
r"""
Show how to compute a derivative spline.
 
Scipy's splines are represented in terms of the standard B-spline basis
functions. In short, a spline of degree ``k`` is represented in terms of the
knots ``t`` and coefficients ``c`` by:
 
.. math::
 
s(x) = \sum_{j=-\infty}^\infty c_{j} B^k_{j}(x)
\\
B^0_i(x) = 1, \text{if $t_i \le x < t_{i+1}$, otherwise $0$,}
\\
B^k_i(x) = \frac{x - t_i}{t_{i+k} - t_i} B^{k-1}_i(x)
+ \frac{t_{i+k+1} - x}{t_{i+k+1} - t_{i+1}} B^{k-1}_{i+1}(x)
 
where :math:`c_j` is nonzero only for ``0 <= j <= N`` and the first
and last `k` knots are at the same position and only used to set up
boundary conditions; terms with vanishing denominators are omitted.
 
One can follow standard spline textbooks here, and work out that by
differentiating this, one obtains another B-spline, of one order
lower:
 
.. math::
 
s'(x) = \sum_{j=-\infty}^\infty d_j B^{k-1}_j(x)
\\
d_j = k \frac{c_{j+1} - c_{j}}{t_{j+k+1} - t_{j+1}}
 
Care needs to be paid at the endpoints: the first knot
should be thrown away since the order is reduced by one.
 
"""
import numpy as np
from scipy.interpolate import UnivariateSpline
 
import matplotlib.pyplot as plt
 
 
class MyUnivariateSpline(UnivariateSpline):
@classmethod
def _from_tck(cls, t, c, k):
self = cls.__new__(cls)
self._eval_args = t, c, k
#_data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
self._data = [None,None,None,None,None,k,None,len(t),t,
c,None,None,None,None]
return self
 
def derivative_spline(self):
"""
Compute the derivative spline of a spline in FITPACK tck
representation.
"""
t, c, k = self._eval_args
if k <= 0:
raise ValueError("Cannot differentiate order 0 spline")
# Compute the denominator in the differentiation formula.
dt = t[k+1:-1] - t[1:-k-1]
# Compute the new coefficients
d = (c[1:-1-k] - c[:-2-k]) * k / dt
# Adjust knots
t2 = t[1:-1]
# Pad coefficient array to same size as knots (FITPACK convention)
d = np.r_[d, [0]*k]
# Done, return a new spline
return self._from_tck(t2, d, k-1)
 
def main():
x = np.linspace(0, 20, 200)
y = np.sin(x)
 
# Fit a spline (must use order 4, since roots() works only for order 3)
s1 = MyUnivariateSpline(x, y, s=0, k=3+1)
 
# Compute the derivative spline of the fitted spline
s2 = s1.derivative_spline()
 
# Roots (works only with order 3 splines)
r = s2.roots()
print "minima/maxima: x =", r/np.pi
print "min/max values: y =", np.sin(r)
 
# Plot and compare the derivative
xx = np.linspace(0, 20, 5000)
plt.plot(xx, np.cos(xx), 'c-', # exact result
xx, s2(xx), 'k--', # derivative spline
xx, np.sin(xx), 'm:', # original curve
r, s1(r), 'k.', # minima/maxima
)
plt.ylim(-1.1, 1.1)
plt.savefig('out.png', dpi=96)
plt.show()
 
if __name__ == "__main__":
main()

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.