2017-07-14 54 views
1

我想用B樣條來計算numpy中閉合曲線的曲率。我想評估樣條表示中的衍生物而不是數據以得到平滑的結果。但是下面的代碼返回一個錯誤:如何從splprep評估樣條曲線?

Traceback (most recent call last): 
     File "<stdin>", line 16, in <module> 
     File "/Users/jfl/anaconda/lib/python2.7/site-packages/scipy/interpolate/fitpack.py", line 1212, in splder 
     c = (c[1:-1-k] - c[:-2-k]) * k/dt 
    TypeError: unsupported operand type(s) for -: 'list' and 'list' 

shell returned 1 

所以我如果我使用正確的splder功能不知道...

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


t = np.linspace(0,2*np.pi,100) 
x = np.cos(t) 
y = np.sin(t) 


pts = np.vstack((x,y)) 
tck, u = intplt.splprep(pts, u=None,k=3, s=0.0, per=1) 
u_new = np.linspace(u.min(), u.max(), 1000) 
x_new, y_new = intplt.splev(u_new, tck, der=0) 

tck_der1 = intplt.splder(tck) 
tck_der2 = intplt.splder(tck_der1) 


xp, yp = intplt.splev(u_new, tck_der1, der=0) 
xpp, ypp = intplt.splev(u_new, tck_der2, der=0) 



plt.figure() 
plt.plot(x,y,".") 
plt.plot(x_new,y_new) 

plt.figure() 
curvature = np.abs(xp* ypp - yp* xpp)/np.power(xp** 2 + yp** 2, 3/2) 
plt.plot(u_new,curvature) 

plt.show() 

回答

1

splder只支持標樣條(插值功能y = f(x)splrep) ,而不是使用splprep獲得的矢量樣條。在任何情況下,你不需要它:只使用der參數splev

xp, yp = intplt.splev(u_new, tck, der=1) 
xpp, ypp = intplt.splev(u_new, tck, der=2) 

這不是衍生物的通用有限差分評價; splev實際上將使用樣條結構來計算它,這就是你想要做的。

我猜你試過以上,並不滿與輸出:

output

但這僅僅是Matplotlib是滑稽。繪圖窗口從+ 9.998e-1到+ 9.998e-1 + 0.0006。換句話說,曲率幾乎是恆定的(1),所以matplotlib放大了來自插值的不可避免的噪聲。只需設置合理的繪圖窗口,問題就會消失。

nice

(或使用類似x = 2*np.cos(t)爲你的榜樣,讓漂亮的非恆定的曲率。)

+0

哇好,我感到非常愚蠢的,我確實試過,但沒有注意到軸極限..無論如何,它很好地瞭解der方法的工作原理! – Jack