在 B 样条基础上查询双变量样条上的多个点

我需要构建一个 3D B-spline 曲面并在各种参数坐标下对其进行多次采样.我找到的最接近的解决方案是使用 bisplev，它需要 tck 输入，由 bsplprep.不幸的是，我不能使用那个 tck 组件，因为它产生一个通过控制点的表面，而我想要的是在 B-spline 基础中计算的表面代码>.所以我手动构建了 tck 输入 bsplev 可以用来产生所需的表面.

I need to construct a 3D B-spline surface and sample it multiple times at various parametric coordinates. The nearest solution I found was to use bisplev, which expects a tck input calculated by bsplprep. Unfortunately i cannot use that tck component because it produces a surface that passes through the control points, while what I want is a surfaces computed in the B-spline basis. So I manually construct the tck input bsplev can use to produce the desired surface.

不幸的是，如果不使用 2 个嵌套循环，我无法弄清楚如何做到这一点:每个 uv 查询 1 个，每个空间组件一个.后者是可以接受的，但前者在处理非常大的查询数组时很慢.

Unfortunately, I can't figure out how to do this without using 2 nested loops: 1 for each uv query, and one for each spacial component. The latter is acceptable but the former is very slow when dealing with very large query arrays.

代码如下:

import numpy as np
import scipy.interpolate as si

def bivariate_bspline(cv,u,v,uCount,vCount,uDegree,vDegree):
    # cv = grid of control vertices
    # u,v = list of u,v component queries
    # uCount, vCount = number of control points along the u and v directions
    # uDegree, vDegree = curve degree along the u and v directions
    
    uMax = uCount-uDegree # Max u parameter
    vMax = vCount-vDegree # Max v parameter
    
    # Calculate knot vectors for both u and v
    u_kv = np.clip(np.arange(uCount+uDegree+1)-uDegree,0,uCount-uDegree) # knot vector in the u direction
    v_kv = np.clip(np.arange(vCount+vDegree+1)-vDegree,0,vCount-vDegree) # knot vector in the v direction
    
    # Compute queries
    position = np.empty((u.shape[0], cv.shape[1]))
    for i in xrange(cv.shape[1]):
        tck = (u_kv, v_kv, cv[:,i], uDegree,vDegree)
        
        for j in xrange(u.shape[0]):
            position[j,i] = si.bisplev(u[j],v[j], tck)

    return position

测试:

# 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       ]])

uCount = 4
vCount = 5
uDegree = 3
vDegree = 3

n = 10**4 # make 10k random queries
u = np.random.random(n) * (uCount-uDegree) 
v = np.random.random(n) * (vCount-vDegree) 
bivariate_bspline(cv,u,v,uCount,vCount,uDegree,vDegree) # will return n correct samples on a b-spline basis surface

速度测试:

import cProfile
cProfile.run('bivariate_bspline(cv,u,v,uCount,vCount,uDegree,vDegree)') # 0.929 seconds

对于 10k 个样本来说，将近 1 秒，其中 bisplev 调用占用了大部分计算时间，因为每个空间组件都调用了 10k 次.

So nearly 1 sec for 10k samples, where the bisplev call takes the lion's share of computation time because it's being called 10k times for each space component.

我确实尝试将 for j in xrange(u.shape[0]): 循环替换为单个 bisplev 调用，并为其提供 u 和 v 数组一次，但这会在 scipy\interpolate\_fitpack_impl.py"，第 1048 行，在 bisplev 中引发 ValueError: Invalid input data.

I did try to replace the for j in xrange(u.shape[0]): loop with a single bisplev call giving it the u and v arrays in one go, but that raises a ValueError: Invalid input data at scipy\interpolate\_fitpack_impl.py", line 1048, in bisplev.

有没有办法摆脱这两者，或者至少摆脱 uv 查询循环并在单个向量化操作中执行所有 uv 查询?

Is there a way to get rid both, or at minimum the uv query loop and do all the uv queries in a single vectorized operation?