【发布时间】：2016-01-30 16:24:04
【问题描述】：
给定一个定义样条曲线的控制顶点列表和查询值列表（0 = 行开始，1= 行结束，0.5= 中途，0.25= 中途等...）我想找到（尽可能快速有效地）样条上这些查询的 3D 坐标。
我试图找到一些内置的 scipy 但失败了。所以我写了一个函数来用蛮力方法解决这个问题：
- 将样条细分为n个细分
- 添加细分以比较查询
- 为每个查询找到适当的细分步骤并推断位置
下面的代码运行良好，但我很想知道是否有更快/更有效的方法来计算我需要的东西，或者更好的一些已经内置的 scipy 我可能错过了。
这是我的功能：
import numpy as np
import scipy.interpolate as interpolate
def uQuery(cv,u,steps=100,projection=True):
''' Brute force point query on spline
cv = list of spline control vertices
u = list of queries (0-1)
steps = number of curve subdivisions (higher value = more precise result)
projection = method by wich we get the final result
- True : project a query onto closest spline segments.
this gives good results but requires a high step count
- False: modulates the parametric samples and recomputes new curve with splev.
this can give better results with fewer samples.
definitely works better (and cheaper) when dealing with b-splines (not in this examples)
'''
u = np.clip(u,0,1) # Clip u queries between 0 and 1
# Create spline points
samples = np.linspace(0,1,steps)
tck,u_=interpolate.splprep(cv.T,s=0.0)
p = np.array(interpolate.splev(samples,tck)).T
# at first i thought that passing my query list to splev instead
# of np.linspace would do the trick, but apparently not.
# Approximate spline length by adding all the segments
p_= np.diff(p,axis=0) # get distances between segments
m = np.sqrt((p_*p_).sum(axis=1)) # segment magnitudes
s = np.cumsum(m) # cumulative summation of magnitudes
s/=s[-1] # normalize distances using its total length
# Find closest index boundaries
s = np.insert(s,0,0) # prepend with 0 for proper index matching
i0 = (s.searchsorted(u,side='left')-1).clip(min=0) # Find closest lowest boundary position
i1 = i0+1 # upper boundary will be the next up
# Return projection on segments for each query
if projection:
return ((p[i1]-p[i0])*((u-s[i0])/(s[i1]-s[i0]))[:,None])+p[i0]
# Else, modulate parametric samples and and pass back to splev
mod = (((u-s[i0])/(s[i1]-s[i0]))/steps)+samples[i0]
return np.array(interpolate.splev(mod,tck)).T
这是一个用法示例：
import matplotlib.pyplot as plt
cv = np.array([[ 50., 25., 0.],
[ 59., 12., 0.],
[ 50., 10., 0.],
[ 57., 2., 0.],
[ 40., 4., 0.],
[ 40., 14., 0.]])
# Lets plot a few queries
u = [0.,0.2,0.3,0.5,1.0]
steps = 10000 # The more subdivisions the better
x,y,z = uQuery(cv,u,steps).T
fig, ax = plt.subplots()
ax.plot(x, y, 'bo')
for i, txt in enumerate(u):
ax.annotate(' u=%s'%txt, (x[i],y[i]))
# Plot the curve we're sampling
tck,u_=interpolate.splprep(cv.T,s=0.0)
x,y,z = np.array(interpolate.splev(np.linspace(0,1,1000),tck))
plt.plot(x,y,'k-',label='Curve')
# Plot control points
p = cv.T
plt.scatter(p[0],p[1],s=80, facecolors='none', edgecolors='r',label='Control Points')
plt.minorticks_on()
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(35, 70)
plt.ylim(0, 30)
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
【问题讨论】：