如何直接获取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(xh,y,z))/(2*h)。 （是的，我知道我可以将额外的抽头数量减少到 3 次并进行“单侧”差异，但不对称性会让我烦恼。）

我的直觉是应该有一种更直接的方法来获取梯度 但我还不知道足够的样条数学来弄清楚，或者理解什么 Scipy 实现的内部正在进行： scipy/scipy/ndimage /src/ni_interpolation.c。

有没有比中心差分更好的“更直接”获得梯度的方法？最好是允许使用现有功能而不是侵入 Scipy 的内部结构来获取它们。

Answer 1

啊哈：根据经典论文在 numpy 代码中引用的样条上，n 阶样条及其导数通过

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

使用 SciPy 的样条插值 I 相关联可以通过维护低阶预过滤量并针对每个导数查询几次来获得我的导数。这意味着添加相当数量的内存（可能与缓存的“主”卷竞争），但大概对低阶样条的评估速度更快，所以对我来说，它是否会比中央差分总体上更快还是不明显使用我目前正在做的小偏移量。还没试过。