2014-03-18 27 views
5

我正在使用UnivariateSpline爲我有的一些數據構造分段多項式。然後我想在其他程序中使用這些樣條曲線(在C或FORTRAN中),所以我想了解生成的樣條曲線背後的公式。從UnivariateSpline對象獲取樣條方程

這裏是我的代碼:

import numpy as np 
import scipy as sp 
from scipy.interpolate import UnivariateSpline 
import matplotlib.pyplot as plt 
import bisect 

data = np.loadtxt('test_C12H26.dat') 
Tmid = 800.0 
print "Tmid", Tmid 
nmid = bisect.bisect(data[:,0],Tmid) 
fig = plt.figure() 
plt.plot(data[:,0], data[:,7],ls='',marker='o',markevery=20) 
npts = len(data[:,0]) 
#print "npts", npts 
w = np.ones(npts) 
w[0] = 100 
w[nmid] = 100 
w[npts-1] = 100 
spline1 = UnivariateSpline(data[:nmid,0],data[:nmid,7],s=1,w=w[:nmid]) 
coeffs = spline1.get_coeffs() 
print coeffs 
print spline1.get_knots() 
print spline1.get_residual() 
print coeffs[0] + coeffs[1] * (data[0,0] - data[0,0]) \ 
       + coeffs[2] * (data[0,0] - data[0,0])**2 \ 
       + coeffs[3] * (data[0,0] - data[0,0])**3, \ 
     data[0,7] 
print coeffs[0] + coeffs[1] * (data[nmid,0] - data[0,0]) \ 
       + coeffs[2] * (data[nmid,0] - data[0,0])**2 \ 
       + coeffs[3] * (data[nmid,0] - data[0,0])**3, \ 
     data[nmid,7] 

print Tmid,data[-1,0] 
spline2 = UnivariateSpline(data[nmid-1:,0],data[nmid-1:,7],s=1,w=w[nmid-1:]) 
print spline2.get_coeffs() 
print spline2.get_knots() 
print spline2.get_residual() 
plt.plot(data[:,0],spline1(data[:,0])) 
plt.plot(data[:,0],spline2(data[:,0])) 
plt.savefig('test.png') 

splines

這裏是導致情節。我相信每個區間都有有效的樣條線,但看起來像我的樣條方程不正確......我找不到它在scipy文檔中應該是什麼的任何參考。有人知道嗎?謝謝 !

回答

5

scipy documentation對於如何獲取係數和手動生成樣條曲線沒有任何說法。但是,可以從現有的B樣條文獻中找出如何做到這一點。下面的函數bspleval顯示如何通過乘以具有最高階基函數的係數和相加來構造B樣條基函數(在代碼中的矩陣B),從中可以很容易地生成樣條曲線:

def bspleval(x, knots, coeffs, order, debug=False): 
    ''' 
    Evaluate a B-spline at a set of points. 

    Parameters 
    ---------- 
    x : list or ndarray 
     The set of points at which to evaluate the spline. 
    knots : list or ndarray 
     The set of knots used to define the spline. 
    coeffs : list of ndarray 
     The set of spline coefficients. 
    order : int 
     The order of the spline. 

    Returns 
    ------- 
    y : ndarray 
     The value of the spline at each point in x. 
    ''' 

    k = order 
    t = knots 
    m = alen(t) 
    npts = alen(x) 
    B = zeros((m-1,k+1,npts)) 

    if debug: 
     print('k=%i, m=%i, npts=%i' % (k, m, npts)) 
     print('t=', t) 
     print('coeffs=', coeffs) 

    ## Create the zero-order B-spline basis functions. 
    for i in range(m-1): 
     B[i,0,:] = float64(logical_and(x >= t[i], x < t[i+1])) 

    if (k == 0): 
     B[m-2,0,-1] = 1.0 

    ## Next iteratively define the higher-order basis functions, working from lower order to higher. 
    for j in range(1,k+1): 
     for i in range(m-j-1): 
      if (t[i+j] - t[i] == 0.0): 
       first_term = 0.0 
      else: 
       first_term = ((x - t[i])/(t[i+j] - t[i])) * B[i,j-1,:] 

      if (t[i+j+1] - t[i+1] == 0.0): 
       second_term = 0.0 
      else: 
       second_term = ((t[i+j+1] - x)/(t[i+j+1] - t[i+1])) * B[i+1,j-1,:] 

      B[i,j,:] = first_term + second_term 
     B[m-j-2,j,-1] = 1.0 

    if debug: 
     plt.figure() 
     for i in range(m-1): 
      plt.plot(x, B[i,k,:]) 
     plt.title('B-spline basis functions') 

    ## Evaluate the spline by multiplying the coefficients with the highest-order basis functions. 
    y = zeros(npts) 
    for i in range(m-k-1): 
     y += coeffs[i] * B[i,k,:] 

    if debug: 
     plt.figure() 
     plt.plot(x, y) 
     plt.title('spline curve') 
     plt.show() 

    return(y) 

爲了舉例說明Scipy現有的單變量樣條函數如何使用它,以下是一個示例腳本。這需要輸入數據,並使用Scipy的功能和麪向對象的方法進行樣條擬合。從兩者中的任何一個取係數和結點,並將它們用作我們手動計算的例程bspleval的輸入,我們再現與它們相同的曲線。請注意,手動評估曲線和Scipy評估方法之間的差異非常小,幾乎可以肯定是浮點噪聲。

x = array([-273.0, -176.4, -79.8, 16.9, 113.5, 210.1, 306.8, 403.4, 500.0]) 
y = array([2.25927498e-53, 2.56028619e-03, 8.64512988e-01, 6.27456769e+00, 1.73894734e+01, 
     3.29052124e+01, 5.14612316e+01, 7.20531200e+01, 9.40718450e+01]) 

x_nodes = array([-273.0, -263.5, -234.8, -187.1, -120.3, -34.4, 70.6, 194.6, 337.8, 500.0]) 
y_nodes = array([2.25927498e-53, 3.83520726e-46, 8.46685318e-11, 6.10568083e-04, 1.82380809e-01, 
       2.66344008e+00, 1.18164677e+01, 3.01811501e+01, 5.78812583e+01, 9.40718450e+01]) 

## Now get scipy's spline fit. 
k = 3 
tck = splrep(x_nodes, y_nodes, k=k, s=0) 
knots = tck[0] 
coeffs = tck[1] 
print('knot points=', knots) 
print('coefficients=', coeffs) 

## Now try scipy's object-oriented version. The result is exactly the same as "tck": the knots are the 
## same and the coeffs are the same, they are just queried in a different way. 
uspline = UnivariateSpline(x_nodes, y_nodes, s=0) 
uspline_knots = uspline.get_knots() 
uspline_coeffs = uspline.get_coeffs() 

## Here are scipy's native spline evaluation methods. Again, "ytck" and "y_uspline" are exactly equal. 
ytck = splev(x, tck) 
y_uspline = uspline(x) 
y_knots = uspline(knots) 

## Now let's try our manually-calculated evaluation function. 
y_eval = bspleval(x, knots, coeffs, k, debug=False) 

plt.plot(x, ytck, label='tck') 
plt.plot(x, y_uspline, label='uspline') 
plt.plot(x, y_eval, label='manual') 

## Next plot the knots and nodes. 
plt.plot(x_nodes, y_nodes, 'ko', markersize=7, label='input nodes')   ## nodes 
plt.plot(knots, y_knots, 'mo', markersize=5, label='tck knots') ## knots 
plt.xlim((-300.0,530.0)) 
plt.legend(loc='best', prop={'size':14}) 

plt.figure() 
plt.title('difference') 
plt.plot(x, ytck-y_uspline, label='tck-uspl') 
plt.plot(x, ytck-y_eval, label='tck-manual') 
plt.legend(loc='best', prop={'size':14}) 

plt.show() 

enter image description here enter image description here

2

通過get_coeffs給出的係數B樣條(基本樣條)係數,說明如下:B-spline (Wikipedia)

也許無論你將使用其他程序/語言都有的實現。提供結點位置和係數,你應該全部設置。