使用numpy和/或scipy插值三维体积

Question

我非常沮丧，因为几个小时后，我似乎无法在python中完成一个看似简单的3D插值。在Matlab中，我所要做的就是

Vi = interp3(x,y,z,V,xi,yi,zi)

使用scipy的ndimage.map_坐标或其他numpy方法与此完全等效的是什么？

谢谢

Answer 1

基本上，ndimage.map_coordinates在“索引”坐标（也就是“体素”或“像素”坐标）中工作。它的界面一开始看起来有点笨拙，但它确实给了你很多灵活性。

如果要指定类似于matlab的interp3的插值坐标，则需要将intput坐标转换为“索引”坐标。

还有一个额外的问题map_coordinates总是保留输出中输入数组的数据类型。如果您插入一个整数数组，您将得到整数输出，这可能是您想要的，也可能不是您想要的。对于下面的代码片段，我假设您总是希望浮点输出。（如果你不这样做，实际上更简单。）

今晚晚些时候我将尝试添加更多的解释（这是相当密集的代码）。

总而言之，我所拥有的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()

Answer 2

在scipy 0.14或更高版本中，有一个新函数^{}与interp3非常相似。

MATLAB命令Vi = interp3(x,y,z,V,xi,yi,zi)将转换为如下内容：

from numpy import array
from scipy.interpolate import RegularGridInterpolator as rgi
my_interpolating_function = rgi((x,y,z), V)
Vi = my_interpolating_function(array([xi,yi,zi]).T)

这里有一个完整的例子来演示两者；它将帮助您理解它们之间的确切区别。。。

MATLAB代码：

x = linspace(1,4,11);
y = linspace(4,7,22);
z = linspace(7,9,33);
V = zeros(22,11,33);
for i=1:11
    for j=1:22
        for k=1:33
            V(j,i,k) = 100*x(i) + 10*y(j) + z(k);
        end
    end
end
xq = [2,3];
yq = [6,5];
zq = [8,7];
Vi = interp3(x,y,z,V,xq,yq,zq);

结果是Vi=[268 357]，这确实是这两点的值(2,6,8)和(3,5,7)。

SCIPY代码：

from scipy.interpolate import RegularGridInterpolator
from numpy import linspace, zeros, array
x = linspace(1,4,11)
y = linspace(4,7,22)
z = linspace(7,9,33)
V = zeros((11,22,33))
for i in range(11):
    for j in range(22):
        for k in range(33):
            V[i,j,k] = 100*x[i] + 10*y[j] + z[k]
fn = RegularGridInterpolator((x,y,z), V)
pts = array([[2,6,8],[3,5,7]])
print(fn(pts))

又是[268,357]。所以你可以看到一些细微的差别：Scipy使用x，y，z索引顺序，而MATLAB使用y，x，z索引顺序（奇怪地）；在Scipy中，你在单独的步骤中定义一个函数，当你调用它时，坐标被分组为（x1，y1，z1），（x2，y2，z2），。。。而matlab使用（x1，x2，…）（y1，y2，…）（z1，z2，…）。

除此之外，两者相似，同样易于使用。

Answer 3

与MATLAB的interp3等价的精确将使用scipy的^{}进行一次性插值：

import numpy as np
from scipy.interpolate import interpn

Vi = interpn((x,y,z), V, np.array([xi,yi,zi]).T)

MATLAB和scipy的默认方法都是线性插值，这可以通过method参数进行更改。注意，对于三维及三维以上，只有线性和最近邻插值才支持interpn，而MATLAB也支持三次和样条插值。

当在同一个网格上进行多个插值调用时，最好使用插值对象^{}，就像在接受的答案above中一样。interpn在内部使用RegularGridInterpolator。