如何直接获得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)并计算例如来自的x导数(i(x+h,y,z)-i(x-h,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 的样条插值，我还可以通过维护低阶预过滤体积并针对每个导数查询几次来获得导数。这意味着添加相当数量的内存（可能与缓存的“主”卷竞争），但大概对低阶样条的评估速度更快，所以对我来说，它是否会比中央差分总体上更快还是不明显使用我目前正在做的小偏移量。还没试过。