2011-08-01 25 views
6

我有一個很大的標量值的3D numpy數組(如果必須的話,可以稱它爲「音量」)。我想插入一個平滑的標量場在這個連續的不規則的,而不是所有的已知前期非整數xyz座標。如何直接獲取scipy插值的漸變?

現在SciPy的對此的支持只是優秀:我篩選與

filtered_volume = scipy.ndimage.interpolation.spline_filter(volume) 

體積和調用

scipy.ndimage.interpolation.map_coordinates(
    filtered_volume, 
    [[z],[y],[x]], 
    prefilter=False) 

爲(X,Y,Z)的利益得到明顯討人喜歡(平滑等)插值。

到目前爲止這麼好。但是,我的應用程序還需要插值字段的局部導數。目前,我通過中心差分法得到這些數據:我還在6個額外點上採樣該音量(這至少可以用map_coordinates的一個呼叫來完成),並計算例如來自(i(x+h,y,z)-i(x-h,y,z))/(2*h)的x導數。 (是的,我知道我可以將額外的水龍頭數量減少到3,並做「單邊」差異,但不對稱會令我煩惱。)

我的直覺是,應該有一種更直接的方式獲得梯度 但我不知道足夠的樣條數學(還沒有)弄清楚,或者瞭解什麼是 在Scipy實現的內部進行:scipy/scipy/ndimage/src/ni_interpolation.c

有沒有比中心差分「更直接」獲得我的梯度的更好方法?最好能夠使用現有功能獲得它們,而不是在Scipy的內部進行攻擊。

+1

奇妙的問題!我認爲從「filtered_volume」中的樣條係數獲得梯度應該是相當直接的,但恐怕我沒有比你更準確的想法......你可能會問也在scipy郵件列表中。 –

回答

1

啊哈:根據在numpy的代碼引用了classic paper on splines,正的一級樣條及其衍生物被

n   n-1   n-1 
dB (x)/dx = B (x+1/2) - B (x-1/2) 

因此,使用SciPy的的樣條插值我可以也通過保持較低的讓我的衍生品相關預先過濾的體積,並且每個派生物查詢幾次。這意味着需要添加相當數量的內存(可能與高速緩存的「主要」卷相競爭),但推測較低階樣條的評估速度更快,所以對於它而言整體速度是否會比中心差異更快並不明顯使用我目前正在做的小偏移量。還沒有嘗試過。