Python 使用 numpy 和或 scipy 插入 3D 体积
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/21836067/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
interpolate 3D volume with numpy and or scipy
提问by user1301295
I am extremely frustrated because after several hours I can't seem to be able to do a seemingly easy 3D interpolation in python. In Matlab all I had to do was
我非常沮丧，因为几个小时后我似乎无法在 python 中进行看似简单的 3D 插值。在 Matlab 中，我所要做的就是
Vi = interp3(x,y,z,V,xi,yi,zi)
What is the exact equivalent of this using scipy's ndimage.map_coordinate or other numpy methods?
使用 scipy 的 ndimage.map_coordinate 或其他 numpy 方法的确切等价物是什么？
Thanks
谢谢
采纳答案by Steve Byrnes
In scipy 0.14 or later, there is a new function scipy.interpolate.RegularGridInterpolatorwhich closely resembles interp3.
在 scipy 0.14 或更高版本中，有一个scipy.interpolate.RegularGridInterpolator与interp3.
The MATLAB command Vi = interp3(x,y,z,V,xi,yi,zi)would translate to something like:
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)
Here is a full example demonstrating both; it will help you understand the exact differences...
这是一个完整的示例，演示了两者；它将帮助您了解确切的差异......
MATLAB CODE:
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);
The result is Vi=[268 357]which is indeed the value at those two points (2,6,8)and (3,5,7).
结果是Vi=[268 357]这确实是这两个点的值(2,6,8)和(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))
Again it's [268,357]. So you see some slight differences: Scipy uses x,y,z index order while MATLAB uses y,x,z (strangely); In Scipy you define a function in a separate step and when you call it, the coordinates are grouped like (x1,y1,z1),(x2,y2,z2),... while matlab uses (x1,x2,...),(y1,y2,...),(z1,z2,...).
又是[268,357]。所以你会看到一些细微的差异：Scipy 使用 x,y,z 索引顺序，而 MATLAB 使用 y,x,z（奇怪）；在 Scipy 中，您在单独的步骤中定义一个函数，当您调用它时，坐标分组为 (x1,y1,z1),(x2,y2,z2),... 而 matlab 使用 (x1,x2,.. .),(y1,y2,...),(z1,z2,...)。
Other than that, the two are similar and equally easy to use.
除此之外，两者相似且同样易于使用。
回答by Joe Kington
Basically, ndimage.map_coordinatesworks in "index" coordinates (a.k.a. "voxel" or "pixel" coordinates). The interface to it seems a bit clunky at first, but it does give you a lotof flexibility.
基本上，ndimage.map_coordinates适用于“索引”坐标（又名“体素”或“像素”坐标）。它的界面起初似乎有点笨拙，但它确实为您提供了很大的灵活性。
If you want to specify the interpolated coordinates similar to matlab's interp3, then you'll need to convert your intput coordinates into "index" coordinates.
如果要指定类似于 matlab 的插值坐标interp3，则需要将输入坐标转换为“索引”坐标。
There's also the additional wrinkle that map_coordinatesalways preserves the dtype of the input array in the output. If you interpolate an integer array, you'll get integer output, which may or may not be what you want. For the code snippet below, I'll assume that you always want floating point output. (If you don't, it's actually simpler.)
还有一个额外的皱纹，它map_coordinates总是在输出中保留输入数组的 dtype。如果你插入一个整数数组，你会得到整数输出，这可能是也可能不是你想要的。对于下面的代码片段，我假设您总是想要浮点输出。（如果你不这样做，它实际上更简单。）
I'll try to add more explanation later tonight (this is rather dense code).
今晚晚些时候我会尝试添加更多解释（这是相当密集的代码）。
All in all, the interp3function I have is more complex than it may need to be for your exact purposes. Howver, it should more or less replicate the behavior of interp3as I remember it (ignoring the "zooming" functionality of interp3(data, zoom_factor), which scipy.ndimage.zoomhandles.)
总而言之，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()
回答by buzjwa
The exactequivalent to MATLAB's interp3would be using scipy's interpnfor one-off interpolation:
该确切相当于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)
The default method for both MATLAB and scipy is linear interpolation, and this can be changed with the methodargument. Note that only linear and nearest-neighbor interpolation is supported by interpnfor 3 dimensions and above, unlike MATLAB which supports cubic and spline interpolation as well.
MATLAB 和 scipy 的默认方法是线性插值，这可以通过method参数进行更改。请注意，interpn对于 3 维及以上，仅支持线性和最近邻插值，这与 MATLAB 不同，它也支持三次和样条插值。
When making multiple interpolation calls on the same grid it is preferable to use the interpolation object RegularGridInterpolator, as in the accepted answer above. interpnuses RegularGridInterpolatorinternally.
在同一个网格上进行多个插值调用时，最好使用插值对象RegularGridInterpolator，如上面已接受的答案。内部interpn使用RegularGridInterpolator。
