|
基本上，ndimage.map_coordinates是在“索引”坐标（也就是“体素”或“像素”坐标）下工作的。它的界面一开始看起来有点笨重，但它确实给予了你“很多"灵活性。如果你想指定类似于matlab的interp3的插值坐标，那么你需要把你的输入坐标转换成“索引”坐标。map_coordinates还有一个额外的问题，就是它总是在输出中保留输入数组的dtype。如果你插入一个整型数组，你会得到整型输出，这可能是你想要的，也可能不是。对于下面的代码片段，我假设你总是想要浮点输出。（如果你不想要，实际上它更简单。）我将在今晚晚些时候尝试添加更多的解释（这是相当密集的代码）。总而言之，我所拥有的interp3函数比您所需要的要复杂得多。但是，它应该或多或少地复制了我记忆中的interp3的行为（忽略interp3(data, zoom_factor)的“缩放”功能，scipy.ndimage.zoom处理该功能）。 import numpy as np
from scipy.ndimage import map_coordinates
def main():
data = np.arange(5*4*3).reshape(5,4,3)
x = np.linspace(5, 10, data.shape[0])
y = np.linspace(10, 20, data.shape[1])
z = np.linspace(-100, 0, data.shape[2])
# Interpolate at a single point
print interp3(x, y, z, data, 7.5, 13.2, -27)
# Interpolate a region of the x-y plane at z=-25
xi, yi = np.mgrid[6:8:10j, 13:18:10j]
print interp3(x, y, z, data, xi, yi, -25 * np.ones_like(xi))
def interp3(x, y, z, v, xi, yi, zi,**kwargs):
"""Sample a 3D array "v" with pixel corner locations at "x","y","z" at the
points in "xi", "yi", "zi" using linear interpolation. Additional kwargs
are passed on to ``scipy.ndimage.map_coordinates``."""
def index_coords(corner_locs, interp_locs):
index = np.arange(len(corner_locs))
if np.all(np.diff(corner_locs) < 0):
corner_locs, index = corner_locs[::-1], index[::-1]
return np.interp(interp_locs, corner_locs, index)
orig_shape = np.asarray(xi).shape
xi, yi, zi = np.atleast_1d(xi, yi, zi)
for arr in [xi, yi, zi]:
arr.shape = -1
output = np.empty(xi.shape, dtype=float)
coords = [index_coords(*item) for item in zip([x, y, z], [xi, yi, zi])]
map_coordinates(v, coords, order=1, output=output,**kwargs)
return output.reshape(orig_shape)
main()
|