scipy - 如何直接获得scipy插值的梯度？

Question

我有一个较大的 3D numpy 标量值数组(如果必须，可以将其称为“卷”)。我想在一系列不规则的序列上插入一个平滑的标量场，而不是全部
已知的前期非整数 xyz 坐标。

现在 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 的内部功能来获取它们的方法。

Answer 1

啊哈:根据classic paper on splines在 numpy 代码中引用，n 阶样条及其导数与

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

因此，使用 SciPy 的样条插值，我还可以通过维护低阶预过滤体积并查询每个导数几次来获得我的导数。这意味着添加相当数量的内存(可能与缓存的“主要”卷竞争)，但大概对低阶样条的评估更快，所以我不清楚它是否会比中心差分更快或不整体使用我目前正在做的小偏移量。还没试过。