我需要构建一个 3DB-spline表面并在各种参数坐标处对其进行多次采样。我发现最接近的解决方案是使用bisplev,它需要一个tck由 计算的输入bsplprep。不幸的是,我无法使用该tck组件,因为它会生成一个通过控制点的表面,而我想要的是一个在 中计算的表面B-spline basis。因此,我手动构建了可用于生成所需表面的tck输入。bsplev
B-spline
bisplev
tck
bsplprep
basis
bsplev
不幸的是,我不知道如何在不使用 2 个嵌套循环的情况下做到这一点:每个uv查询一个,每个空间组件一个。后者是可以接受的,但前者在处理非常大的查询数组时非常慢。
uv
代码如下:
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 次。
我确实尝试for j in xrange(u.shape[0]):用一个bisplev调用来替换循环,一次性给它 u 和 v 数组,但这会引发ValueError: Invalid input dataat scipy\interpolate\_fitpack_impl.py", line 1048, in bisplev。
for j in xrange(u.shape[0]):
ValueError: Invalid input data
scipy\interpolate\_fitpack_impl.py", line 1048, in bisplev
有没有办法摆脱这两者,或者至少摆脱uv查询循环并uv在单个矢量化操作中执行所有查询?
简短回答:替换
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)
和
for i in xrange(cv.shape[1]): position[:, i] = si.dfitpack.bispeu(u_kv, v_kv, cv[:,i], uDegree, vDegree, u, v)[0]
bisplev确实接受数组作为si.bisplev(u, v, tck),但它将它们解释为定义 xy 网格。因此,u和v必须按升序排序,并且将对所有对执行评估(u[j], v[k]),输出为值的 2D 数组。这不是您想要的;对评估次数进行平方可能很糟糕,并且从返回的 2D 数组中提取您真正想要的值并不容易(它们不一定在对角线上,因为您的 u、v 一开始就没有排序)。
si.bisplev(u, v, tck)
u
v
(u[j], v[k])
但是SmoothBivariateSpline 的调用方法包含一个布尔参数grid,将其设置为 False 使其只对(u[j], v[j])样条线进行求值。向量 u, v 不再需要排序,但现在它们必须具有相同的大小。
grid
(u[j], v[j])
但是你正在准备自己的tck。这提供了两种方法:尝试SmoothBivariateSpline使用手工制作的 tck 进行实例化;或者阅读其调用方法的源代码并执行当参数grid设置为 False 时它所做的事情。我采用了后一种方法。
SmoothBivariateSpline