用numpy和scipy插入三维体积

Question

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

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

这是什么使用SciPy的的ndimage.map_coordinate或其他numpy的方法完全等效？

由于

Answer 1

在SciPy的0.14或更高版本，有一个新的功能scipy.interpolate.RegularGridInterpolator这与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]这的确是在这两点​​3210和(3,5,7)值。

SciPy的CODE：

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 2

基本上，ndimage.map_coordinates作品在 “索引” 的坐标（也称为 “体素” 或 “像素” 坐标）。它的界面起初似乎有点笨拙，但它确实给你提供了很多灵活性。

如果要指定类似于matlab的interp3的插补坐标，那么您需要将输入坐标转换为“索引”坐标。

还有一个额外的皱纹map_coordinates始终保留在输出中的输入数组的dtype。如果你插入一个整数数组，你会得到整数输出，这可能是也可能不是你想要的。对于下面的代码片段，我假设你总是需要浮点输出。 （如果你不这样做，这实际上更简单。）

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

总之，我所拥有的interp3函数比它可能需要用于确切目的更复杂。 Howver，我记得它应该或多或少复制的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 3

的确切相当于MATLAB的interp3将使用SciPy的的interpn用于一次性插补：

import numpy as np 
from scipy.interpolate import interpn 

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

两个MATLAB和SciPy的默认方法是线性内插，并且这可以与method来改变论据。请注意，对于3维和更高维度，只有线性和最近邻插值支持interpn，与支持三次和样条插值的MATLAB不同。

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