Compute the B-spline basis of a Bivariate spline

Question

I need to compute uv queries on a bivariate spline in the B-spline basis. With this answer I have a good function (copied below) that leverages scipy.dfitpack.bispeu to get the results I need.

import numpy as np
import scipy.interpolate as si

def fitpack_bispeu(cv, u, v, count_u, count_v, degree_u, degree_v):
    # cv = grid of control vertices
    # u,v = list of u,v component queries
    # count_u, count_v = grid counts along the u and v directions
    # degree_u, degree_v = curve degree along the u and v directions
    
    # Calculate knot vectors for both u and v
    tck_u = np.clip(np.arange(count_u+degree_u+1)-degree_u,0,count_u-degree_u) # knot vector in the u direction
    tck_v = np.clip(np.arange(count_v+degree_v+1)-degree_v,0,count_v-degree_v) # knot vector in the v direction

    # Compute queries
    positions = np.empty((u.shape[0], cv.shape[1]))
    for i in range(cv.shape[1]):
        positions[:, i] = si.dfitpack.bispeu(tck_u, tck_v, cv[:,i], degree_u, degree_v, u, v)[0]

    return positions

This function works, but it occurred to me I could get better performance by calculating the bivariate basis ahead of time and then just get my result via a dot product. Here's what i wrote to compute the basis.

def basis_bispeu(cv, u, v, count_u, count_v, degree_u, degree_v):

    # Calculate knot vectors for both u and v
    tck_u = np.clip(np.arange(count_u+degree_u+1)-degree_u,0,count_u-degree_u) # knot vector in the u direction
    tck_v = np.clip(np.arange(count_v+degree_v+1)-degree_v,0,count_v-degree_v) # knot vector in the v direction

    # Compute basis for each control vertex
    basis = np.empty((u.shape[0], cv.shape[0]))
    cv_   = np.identity(len(cv))
    
    for i in range(cv.shape[0]):
        basis[:,i] = si.dfitpack.bispeu(tck_u, tck_v, cv_[i], degree_u, degree_v, u, v)[0]


    return basis

lets compare and profile with cProfile:

# A test grid of control vertices
cv = np.array([[-0.5       , -0.  ,        0.5       ],
               [-0.5       , -0.  ,        0.33333333],
               [-0.5       , -0.  ,        0.        ],
               [-0.5       ,  0.  ,       -0.33333333],
               [-0.5       ,  0.  ,       -0.5       ],
               [-0.16666667,  1.  ,        0.5       ],
               [-0.16666667, -0.  ,        0.33333333],
               [-0.16666667,  0.5 ,        0.        ],
               [-0.16666667,  0.5 ,       -0.33333333],
               [-0.16666667,  0.  ,       -0.5       ],
               [ 0.16666667, -0.  ,        0.5       ],
               [ 0.16666667, -0.  ,        0.33333333],
               [ 0.16666667, -0.  ,        0.        ],
               [ 0.16666667,  0.  ,       -0.33333333],
               [ 0.16666667,  0.  ,       -0.5       ],
               [ 0.5       , -0.  ,        0.5       ],
               [ 0.5       , -0.  ,        0.33333333],
               [ 0.5       , -0.5 ,        0.        ],
               [ 0.5       ,  0.  ,       -0.33333333],
               [ 0.5       ,  0.  ,       -0.5       ]])

count_u = 4
count_v = 5
degree_u = 3
degree_v = 3

n = 10**6 # make 1 million random queries
u = np.random.random(n) * (count_u-degree_u) 
v = np.random.random(n) * (count_v-degree_v) 

# get the result from fitpack_bispeu
result_bispeu = fitpack_bispeu(cv,u,v,count_u,count_v,degree_u,degree_v) # 0.482 seconds

# precompute the basis for the same grid
basis = basis_bispeu(cv,u,v,count_u,count_v,degree_u,degree_v) # 2.124 seconds

# get results via dot product
result_basis = np.dot(basis,cv) # 0.028 seconds (17x faster than fitpack_bispeu)

# all close?
print np.allclose(result_basis, result_bispeu) # True

With a 17x speed increase, pre-calculating the basis seems like the way to go, but basis_bispeu is rather slow.

QUESTION

Is there a faster way to compute the basis of a bivariate spline? I know of deBoor's algorithm to compute a similar basis on a curve. Are there similar algorithms for bivariates that once written with numba or cython could yield better performance?

Otherwise can the basis_bispeu function above be improved to compute the basis faster? Maybe there are built in numpy functions I'm not aware of that could help.