在 B 样条基础上查询双变量样条上的多个点 [英] Querying multiple points on a bivariate spline in the B-spline basis

问题描述

我需要构建一个 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?

推荐答案

简答:替换

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])，输出是一个二维数组.这不是你想要的；求值次数的平方可能不好，并且从返回的二维数组中提取您实际想要的值并不容易(它们不一定在对角线上，因为您的 u、v 没有开始排序).

Explanation

bisplev does accept arrays as si.bisplev(u, v, tck), but it interprets them as defining an xy-grid. Hence, u and v must be sorted in ascending order, and evaluation will be performed on all pairs (u[j], v[k]), the output being a 2D array of values. This is not what you want; squaring the number of evaluations is probably bad, and it's not going to be easy to extract the values you actually want from the 2D array returned (they are not necessarily on the diagonal, since your u, v were not sorted to begin with).

但是 SmoothBivariateSpline 的调用方法 包含一个布尔参数grid，将which 设置为False 使其仅评估(u[j], v[j])<处的样条曲线/code> 对.向量 u、v 不再需要排序，但现在它们必须具有相同的大小.

But the call method of SmoothBivariateSpline includes a boolean parameter grid, setting which to False makes it just evaluate the spline at (u[j], v[j]) pairs. The vectors u, v no longer need to be sorted, but now they must be of the same size.

但是您正在准备自己的tck.这提供了两种方法:尝试使用手工 tck 实例化 SmoothBivariateSpline；或阅读 它的调用方法 并在参数grid 设置为False 时执行它的操作.我选择了后一种方法.

But you are preparing your own tck. This presents two approaches: seek to instantiate SmoothBivariateSpline with hand-made tck; or read the source of its call method and do what it does when the parameter grid is set to False. I went with the latter approach.

这篇关于在 B 样条基础上查询双变量样条上的多个点的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！