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

【问题标题】：How to obtain the gradient of a scipy interpolant directly?如何直接获得scipy插值的梯度？
【发布时间】：2011-10-17 19:08:36
【问题描述】：

我有一个较大的 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 的内部来获得它们。

【问题讨论】：

好问题！我认为从filtered_volume 中的样条系数获得梯度应该是相当直接的，但恐怕我没有比你更好的方法了......你可能想问在 scipy 邮件列表中也是如此。

标签： scipy gradient interpolation volume spline

【解决方案1】：

啊哈：根据numpy代码中引用的classic paper on splines，n阶样条线及其导数由

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

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

【讨论】：