6

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

4

1 回答 1

1

啊哈：根据numpy代码中引用的关于样条的经典论文，n阶样条及其导数与

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

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

于 2011-08-16T13:50:55.013 回答