python/scipy中的多元样条插值？

Question

是否有库模块或其他直接的方法在python中实现多元样条插值？

具体来说,我在一个规则间隔的三维网格上有一组标量数据,我需要在分散在整个域中的少量点进行插值.对于二维,我一直在使用scipy.interpolate.RectBivariateSpline,我实际上是在寻找三维数据的扩展.

我发现的N维插值例程还不够好:我更喜欢使用LinearNDInterpolator上的样条线来获得平滑度,而且我有太多的数据点(通常超过一百万),例如,径向基函数可以工作.

如果有人知道可以执行此操作的python库,或者我可以调用或移植的另一种语言,我真的很感激.

Answer 1

如果我正确理解您的问题,您的输入"观察"数据会定期网格化吗？

如果是这样,scipy.ndimage.map_coordinates那就完全是你想要的.

在第一次通过时有点难以理解,但实际上,您只需要为其提供一系列坐标,您希望在像素/体素/ n维索引坐标中插入网格的值.

作为2D示例:

import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt

# Note that the output interpolated coords will be the same dtype as your input
# data.  If we have an array of ints, and we want floating point precision in
# the output interpolated points, we need to cast the array as floats
data = np.arange(40).reshape((8,5)).astype(np.float)

# I'm writing these as row, column pairs for clarity...
coords = np.array([[1.2, 3.5], [6.7, 2.5], [7.9, 3.5], [3.5, 3.5]])
# However, map_coordinates expects the transpose of this
coords = coords.T

# The "mode" kwarg here just controls how the boundaries are treated
# mode='nearest' is _not_ nearest neighbor interpolation, it just uses the
# value of the nearest cell if the point lies outside the grid.  The default is
# to treat the values outside the grid as zero, which can cause some edge
# effects if you're interpolating points near the edge
# The "order" kwarg controls the order of the splines used. The default is 
# cubic splines, order=3
zi = ndimage.map_coordinates(data, coords, order=3, mode='nearest')

row, column = coords
nrows, ncols = data.shape
im = plt.imshow(data, interpolation='nearest', extent=[0, ncols, nrows, 0])
plt.colorbar(im)
plt.scatter(column, row, c=zi, vmin=data.min(), vmax=data.max())
for r, c, z in zip(row, column, zi):
    plt.annotate('%0.3f' % z, (c,r), xytext=(-10,10), textcoords='offset points',
            arrowprops=dict(arrowstyle='->'), ha='right')
plt.show()

要在n维中执行此操作,我们只需要传入适当大小的数组:

import numpy as np
from scipy import ndimage

data = np.arange(3*5*9).reshape((3,5,9)).astype(np.float)
coords = np.array([[1.2, 3.5, 7.8], [0.5, 0.5, 6.8]])
zi = ndimage.map_coordinates(data, coords.T)

就缩放和内存使用而言,map_coordinates如果您使用的顺序> 1(即非线性插值),将创建数组的过滤副本.如果您只想在非常少的点进行插值,这是一个相当大的开销.但是,它不会随着您要插入的数字点而增加.只要有足够的RAM用于输入数据阵列的单个临时副本,你就可以了.

如果您无法将数据副本存储在内存中,您可将)指定prefilter=False并order=1和使用线性内插,或b),所用过滤的版本替换您的原始数据ndimage.spline_filter,然后调用map_coordinates用prefilter=False.

即使你有足够的ram,如果你需要多次调用map_coordinates(例如交互式使用等),保持过滤的数据集可以是一个很大的加速.