由于 AttributeError，无法区分 3D 中的样条

Question

我正在尝试对一些数据进行平滑 B 样条拟合，我发现这在此处非常有用 post。但是，我不仅需要样条曲线，还需要它的导数，所以我尝试在示例中添加以下代码：

tck_der = interpolate.splder(tck, n=1)
x_der, y_der, z_der = interpolate.splev(u_fine, tck_der)

由于某些数据类型问题，这似乎无法正常工作。我得到以下回溯：

Traceback (most recent call last):
  File "interpolate_point_trace.py", line 31, in spline_example
    tck_der = interpolate.splder(tck, n=1)
  File "/home/user/anaconda3/lib/python3.7/site-packages/scipy/interpolate/fitpack.py", line 657, in splder
     return _impl.splder(tck, n)
   File "/home/user/anaconda3/lib/python3.7/site-packages/scipy/interpolate/_fitpack_impl.py", line 1206, in splder
     sh = (slice(None),) + ((None,)*len(c.shape[1:]))
 AttributeError: 'list' object has no attribute 'shape'

原因似乎是 tck 元组的第二个参数包含一个 numpy 数组列表。我认为将输入数据转换为 numpy 数组也会有所帮助，但它不会更改 tck.

的数据类型

此行为是否反映了 scipy 中的错误，或者输入是否格式错误？我尝试手动将列表转换为数组：

tck[1] = np.array(tck[1])

但是这个（这并不让我感到惊讶）也给出了一个错误：

ValueError: operands could not be broadcast together with shapes (0,8) (7,1)

知道问题出在哪里吗？我在 1D 样条曲线之前和上面使用过 scipy splder 函数工作得很好，所以我认为它与样条曲线在 3D 中是一条线有关。

--------编辑--------

这是一个最小的工作示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
from mpl_toolkits.mplot3d import Axes3D

total_rad = 10
z_factor = 3
noise = 0.1

num_true_pts = 200
s_true = np.linspace(0, total_rad, num_true_pts)
x_true = np.cos(s_true)
y_true = np.sin(s_true)
z_true = s_true / z_factor

num_sample_pts = 80
s_sample = np.linspace(0, total_rad, num_sample_pts)
x_sample = np.cos(s_sample) + noise * np.random.randn(num_sample_pts)
y_sample = np.sin(s_sample) + noise * np.random.randn(num_sample_pts)
z_sample = s_sample / z_factor + noise * np.random.randn(num_sample_pts)

tck, u = interpolate.splprep([x_sample, y_sample, z_sample], s=2)
x_knots, y_knots, z_knots = interpolate.splev(tck[0], tck)
u_fine = np.linspace(0, 1, num_true_pts)
x_fine, y_fine, z_fine = interpolate.splev(u_fine, tck)

# this is the part of the code I inserted: the line under this causes the crash
tck_der = interpolate.splder(tck, n=1)
x_der, y_der, z_der = interpolate.splev(u_fine, tck_der)

# end of the inserted code

fig2 = plt.figure(2)
ax3d = fig2.add_subplot(111, projection='3d')
ax3d.plot(x_true, y_true, z_true, 'b')
ax3d.plot(x_sample, y_sample, z_sample, 'r*')
ax3d.plot(x_knots, y_knots, z_knots, 'go')
ax3d.plot(x_fine, y_fine, z_fine, 'g')
fig2.show()
plt.show()

Answer 1

遇到同样的问题...

我通过使用 interpolate.splder(tck, n=1) 来规避错误，而是使用 interpolate.splev(spline_ev, tck, der=1)，其中 returns 点 spline_ev 的导数（参见 Scipy Doku）。

如果您需要样条曲线，我想您可以再次使用 interpolate.splprep()。

总的来说：

import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt

points = np.random.rand(10,2) * 10

(tck, u), fp, ier, msg = interpolate.splprep(points.T, s=0, k=3, full_output=True)

spline_ev = np.linspace(0.0, 1.0, 100, endpoint=True)

spline_points = interpolate.splev(spline_ev, tck)

# Calculate derivative
spline_der_points = interpolate.splev(spline_ev, tck, der=1)
spline_der = interpolate.splprep(spline_der_points.T, s=0, k=3, full_output=True)


# Plot the data and derivative
fig = plt.figure()

plt.plot(points[:,0], points[:,1], '.-', label="points")
plt.plot(spline_points[0], spline_points[1], '.-', label="tck")
plt.plot(spline_der_points[0], spline_der_points[1], '.-', label="tck_der")

#   Show tangent
plt.arrow(spline_points[0][23]-spline_der_points[0][23], spline_points[1][23]-spline_der_points[1][23], 2.0*spline_der_points[0][23], 2.0*spline_der_points[1][23])

plt.legend()

plt.show()

编辑：

我还在 Github 上开了一个 Issue，根据 ev-br interpolate.splprep 的用法已贬值，应该使用 make_interp_spline / BSpline 代替。

Answer 2

如其他答案所述，splprep 输出与 splder 不兼容，但与 splev 兼容。而后者可以评估导数。

然而，对于插值，有一种替代方法，它完全避免了splprep。我基本上是在 SciPy 问题跟踪器 (https://github.com/scipy/scipy/issues/10389):

上复制回复

这是复制 splprep 输出的示例。首先让我们理解 splprep 输出：

# start with the OP example
import numpy as np
from scipy import interpolate

points = np.random.rand(10,2) * 10

(tck, u), fp, ier, msg = interpolate.splprep(points.T, s=0, k=3, full_output=True)

# check the meaning of the `u` array: evaluation of the spline at `u`
# gives back the original points (up to a list/transpose)

xy = interpolate.splev(u, tck)
xy = np.asarray(xy)
np.allclose(xy.T, points)

接下来，让我们在没有splprep的情况下复制它。首先，构建 u 数组：曲线以参数方式表示，u 本质上是弧长的近似值。其他参数化也是可能的，但这里让我们坚持 splprep 所做的。从文档页面翻译伪代码，https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.splprep.html

vv = np.sum((points[1:, :] - points[:-1, :])**2, axis=1)
vv = np.sqrt(vv).cumsum()
vv/= vv[-1]
vv = np.r_[0, vv]

# check: 
np.allclose(u, vv)

现在，沿着参数曲线插值：points vs vv:

spl = interpolate.make_interp_spline(vv, points)

# check spl.t vs knots from splPrep
spl.t - tck[0]

结果 spl 是一个 BSpline 对象，您可以用通常的方式对其求值、求导等：

np.allclose(points, spl(vv))

# differentiate
spl_derivative = spl.derivative(vv)

由于 AttributeError，无法区分 3D 中的样条

Spline in 3D can not be differentiated due to an AttributeError

python

spline

scipy

derivative

编辑：