python - How to get the spline basis used by scipy.interpolate.splev -

i need evaluate b-splines in python. wrote code below works well.

import numpy np import scipy.interpolate si  def scipy_bspline(cv,n,degree):     """ bspline basis function         c        = list of control points.         n        = number of points on curve.         degree   = curve degree     """     # create range of u values     c = len(cv)     kv = np.array([0]*degree + range(c-degree+1) + [c-degree]*degree,dtype='int')     u  = np.linspace(0,c-degree,n)      # calculate result     arange = np.arange(n)     points = np.zeros((n,cv.shape[1]))     in xrange(cv.shape[1]):         points[arange,i] = si.splev(u, (kv,cv[:,i],degree))      return points 

giving 6 control points , asking evaluate 100k points on curve breeze:

# control points cv = np.array([[ 50.,  25., 0.],        [ 59.,  12., 0.],        [ 50.,  10., 0.],        [ 57.,   2., 0.],        [ 40.,   4., 0.],        [ 40.,   14., 0.]])  n = 100000  # 100k points degree = 3 # curve degree points_scipy = scipy_bspline(cv,n,degree) #cprofile clocks @ 0.012 seconds 

here's plot of "points_scipy":

now, make faster compute basis 100k points on curve, store in memory, , when need draw curve need multiply new control point positions stored basis new curve. prove point wrote function uses deboor's algorithm compute basis:

def basis(c, n, degree):     """ bspline basis function         c        = number of control points.         n        = number of points on curve.         degree   = curve degree     """     # create knot vector , range of samples on curve     kv = np.array([0]*degree + range(c-degree+1) + [c-degree]*degree,dtype='int') # knot vector     u  = np.linspace(0,c-degree,n) # samples range      # cox - deboor recursive function calculate basis     def coxdeboor(u, k, d):         # test end conditions         if (d == 0):             if (kv[k] <= u , u < kv[k+1]):                 return 1             return 0          den1 = kv[k+d] - kv[k]         den2 = 0         eq1  = 0         eq2  = 0          if den1 > 0:             eq1 = ((u-kv[k]) / den1) * coxdeboor(u,k,(d-1))          try:             den2 = kv[k+d+1] - kv[k+1]             if den2 > 0:                 eq2 = ((kv[k+d+1]-u) / den2) * coxdeboor(u,(k+1),(d-1))         except:             pass          return eq1 + eq2       # compute basis each point     b = np.zeros((n,c))     in xrange(n):         k in xrange(c):             b[i][k%c] += coxdeboor(u[i],k,degree)      b[n-1][-1] = 1      return b 

now lets use compute new basis, multiply control points , confirm we're getting same results splev:

b = basis(len(cv),n,degree) #5600011 function calls (600011 primitive calls) in 10.975 seconds points_basis = np.dot(b,cv) #3 function calls in 0.002 seconds print np.allclose(points_basis,points_scipy) # returns true 

my extremely slow function returned 100k basis values in 11 seconds, since these values need computed once, calculating points on curve ended being 6 times faster way doing via splev.

the fact able exact same results both method , splev leads me believe internally splev calculates basis do, except faster.

so goal find out how calculate basis fast, store in memory , use np.dot() compute new points on curve, , question is: possible use spicy.interpolate basis values (i presume) splev uses compute it's result? , if so, how?

[addendum]

following unutbu's , ev-br's useful insight on how scipy calculates spline's basis, looked fortran code , wrote equivalent best of abilities:

def fitpack_basis(c, n=100, d=3, rminoffset=0, rmaxoffset=0):     """ fitpack's spline basis function         c = number of control points.         n = number of points on curve.         d = curve degree     """     # create knot vector     kv = np.array([0]*d + range(c-d+1) + [c-d]*d, dtype='int')      # create sample range     u = np.linspace(rminoffset, rmaxoffset + c - d, n)  # samples range      # create buffers     b  = np.zeros((n,c)) # basis     bb = np.zeros((n,c)) # basis buffer     left  = np.clip(np.floor(u),0,c-d-1).astype(int)   # left  knot vector indices     right = left+d+1 # right knot vector indices      # go!     nrange = np.arange(n)     b[nrange,left] = 1.0     j in xrange(1, d+1):         crange = np.arange(j)[:,none]         bb[nrange,left+crange] = b[nrange,left+crange]                 b[nrange,left] = 0.0         in xrange(j):             f = bb[nrange,left+i] / (kv[right+i] - kv[right+i-j])             b[nrange,left+i] = b[nrange,left+i] + f * (kv[right+i] - u)             b[nrange,left+i+1] = f * (u - kv[right+i-j])      return b 

testing against unutbu's version of original basis function:

fb = fitpack_basis(c,n,d) #22 function calls in 0.044 seconds b = basis(c,n,d) #81 function calls (45 primitive calls) in 0.013 seconds  ~5 times faster print np.allclose(b,fb) # returns true 

my function 5 times slower, still relatively fast. lets me use sample ranges go beyond boundaries, of use in application. example:

print fitpack_basis(c,5,d,rminoffset=-0.1,rmaxoffset=.2) [[ 1.331  -0.3468  0.0159 -0.0002  0.      0.    ] [ 0.0208  0.4766  0.4391  0.0635  0.      0.    ] [ 0.      0.0228  0.4398  0.4959  0.0416  0.    ] [ 0.      0.      0.0407  0.3621  0.5444  0.0527] [ 0.      0.     -0.0013  0.0673 -0.794   1.728 ]] 

so reason use fitpack_basis since it's relatively fast. love suggestions improving performance , closer unutbu's version of original basis function wrote.

fitpack_basis uses double loop iteratively modifies elements in bb , b. don't see way use numpy vectorize these loops since value of bb , b @ each stage of iteration depends on values previous iterations. in situations cython can used improve performance of loops.

here cython-ized version of fitpack_basis, runs fast bspline_basis. main ideas used boost speed using cython declare type of every variable, , rewrite uses of numpy fancy indexing loops using plain integer indexing.

see this page instructions on how build code , run python.

import numpy np cimport numpy np cimport cython  ctypedef np.float64_t dtype_f ctypedef np.int64_t dtype_i  @cython.boundscheck(false) @cython.wraparound(false) @cython.nonecheck(false) def cython_fitpack_basis(int c, int n=100, int d=3,                           double rminoffset=0, double rmaxoffset=0):     """ fitpack's spline basis function         c = number of control points.         n = number of points on curve.         d = curve degree     """     cdef py_ssize_t i, j, k, l     cdef double f      # create knot vector     cdef np.ndarray[dtype_i, ndim=1] kv = np.array(         [0]*d + range(c-d+1) + [c-d]*d, dtype=np.int64)      # create sample range     cdef np.ndarray[dtype_f, ndim=1] u = np.linspace(         rminoffset, rmaxoffset + c - d, n)      # basis     cdef np.ndarray[dtype_f, ndim=2] b  = np.zeros((n,c))      # basis buffer     cdef np.ndarray[dtype_f, ndim=2] bb = np.zeros((n,c))      # left  knot vector indices     cdef np.ndarray[dtype_i, ndim=1] left = np.clip(np.floor(u), 0, c-d-1).astype(np.int64)        # right knot vector indices     cdef np.ndarray[dtype_i, ndim=1] right = left+d+1       k in range(n):         b[k, left[k]] = 1.0      j in range(1, d+1):         l in range(j):             k in range(n):                 bb[k, left[k] + l] = b[k, left[k] + l]                  b[k, left[k]] = 0.0         in range(j):             k in range(n):                 f = bb[k, left[k]+i] / (kv[right[k]+i] - kv[right[k]+i-j])                 b[k, left[k]+i] = b[k, left[k]+i] + f * (kv[right[k]+i] - u[k])                 b[k, left[k]+i+1] = f * (u[k] - kv[right[k]+i-j])     return b 

using timeit code benchmark it's performance,

import timeit import numpy np import cython_bspline cb import numpy_bspline nb  c = 6 n = 10**5 d = 3  fb = nb.fitpack_basis(c, n, d) bb = nb.bspline_basis(c, n, d)  cfb = cb.cython_fitpack_basis(c,n,d)   assert np.allclose(bb, fb)  assert np.allclose(cfb, fb)  # print(nb.fitpack_basis(c,5,d,rminoffset=-0.1,rmaxoffset=.2))  timing = dict() timing['nb.fitpack_basis'] = timeit.timeit(     stmt='nb.fitpack_basis(c, n, d)',      setup='from __main__ import nb, c, n, d',      number=10) timing['nb.bspline_basis'] = timeit.timeit(     stmt='nb.bspline_basis(c, n, d)',      setup='from __main__ import nb, c, n, d',      number=10) timing['cb.cython_fitpack_basis'] = timeit.timeit(     stmt='cb.cython_fitpack_basis(c, n, d)',      setup='from __main__ import cb, c, n, d',      number=10)  func_name, t in timing.items():     print "{:>25}: {:.4f}".format(func_name, t) 

it appears cython can make fitpack_basis code run fast (and perhaps bit faster) bspline_basis:

         nb.bspline_basis: 0.3322   cb.cython_fitpack_basis: 0.2939          nb.fitpack_basis: 0.9182