我正在为二进制发行版中包含NumPy的应用程序编写插件，但不是SciPy。我的插件需要将数据从一个常规3D网格插入另一个常规3D网格。从源代码运行，这可以非常有效地使用 scipy.ndimage 或者，如果用户未安装SciPy，则生成编织 .pyd 我写过。不幸的是，如果用户正在运行二进制文件，那么这两个选项都不可用。
我写了一个简单的 三线插值 在python中提供正确结果的例程，但对于我正在使用的数组大小，需要很长时间（约5分钟）。我想知道是否有办法只使用NumPy中的功能加快速度。就像 scipy.ndimage.map_coordinates，它需要一个3D输入数组和一个数组，每个点的x，y和z坐标都要进行插值。
def trilinear_interp(input_array, indices):
"""Evaluate the input_array data at the indices given"""
output = np.empty(indices[0].shape)
x_indices = indices[0]
y_indices = indices[1]
z_indices = indices[2]
for i in np.ndindex(x_indices.shape):
x0 = np.floor(x_indices[i])
y0 = np.floor(y_indices[i])
z0 = np.floor(z_indices[i])
x1 = x0 + 1
y1 = y0 + 1
z1 = z0 + 1
#Check if xyz1 is beyond array boundary:
if x1 == input_array.shape[0]:
x1 = x0
if y1 == input_array.shape[1]:
y1 = y0
if z1 == input_array.shape[2]:
z1 = z0
x = x_indices[i] - x0
y = y_indices[i] - y0
z = z_indices[i] - z0
output[i] = (input_array[x0,y0,z0]*(1-x)*(1-y)*(1-z) +
input_array[x1,y0,z0]*x*(1-y)*(1-z) +
input_array[x0,y1,z0]*(1-x)*y*(1-z) +
input_array[x0,y0,z1]*(1-x)*(1-y)*z +
input_array[x1,y0,z1]*x*(1-y)*z +
input_array[x0,y1,z1]*(1-x)*y*z +
input_array[x1,y1,z0]*x*y*(1-z) +
input_array[x1,y1,z1]*x*y*z)
return output
显然，功能如此缓慢的原因是 for 循环遍历3D空间中的每个点。有没有办法执行某种切片或矢量化魔术来加速它？谢谢。
事实证明，它很容易被矢量化。
output = np.empty(indices[0].shape)
x_indices = indices[0]
y_indices = indices[1]
z_indices = indices[2]
x0 = x_indices.astype(np.integer)
y0 = y_indices.astype(np.integer)
z0 = z_indices.astype(np.integer)
x1 = x0 + 1
y1 = y0 + 1
z1 = z0 + 1
#Check if xyz1 is beyond array boundary:
x1[np.where(x1==input_array.shape[0])] = x0.max()
y1[np.where(y1==input_array.shape[1])] = y0.max()
z1[np.where(z1==input_array.shape[2])] = z0.max()
x = x_indices - x0
y = y_indices - y0
z = z_indices - z0
output = (input_array[x0,y0,z0]*(1-x)*(1-y)*(1-z) +
input_array[x1,y0,z0]*x*(1-y)*(1-z) +
input_array[x0,y1,z0]*(1-x)*y*(1-z) +
input_array[x0,y0,z1]*(1-x)*(1-y)*z +
input_array[x1,y0,z1]*x*(1-y)*z +
input_array[x0,y1,z1]*(1-x)*y*z +
input_array[x1,y1,z0]*x*y*(1-z) +
input_array[x1,y1,z1]*x*y*z)
return output
非常感谢这篇文章，并对此进行了跟进。我已经将自己基于你的矢量化，以便给它提供另一个速度提升（至少我正在使用的数据）！
我正在处理图像相关，因此我在同一个内插多组不同的坐标 input_array。
不幸的是我让它变得有点复杂，但是如果我能解释我做了什么，额外的复杂性应该a）证明自己是正确的，并且b）变得清晰。你的最后一行（ 输出= ）仍然需要在非连续的地方查找相当数量的 input_array，使它相对缓慢。
假设我的3D数据长度为NxMxP。我决定做以下事情：如果我能得到一个（8 x（NxMxP））矩阵的预先计算的灰度值为一个点及其最近的邻居，我也可以计算（（NxMxP）X 8）矩阵系数（你上面例子中的第一个系数是（x-1）（y-1）（z-1））然后我可以一起乘以一起，然后回家免费！
对我来说一个很好的奖励是我可以预先计算灰色矩阵并回收它！
这是一个代码示例（从两个不同的函数粘贴，所以可能不是开箱即用，但应该作为一个很好的灵感来源）：
def trilinear_interpolator_speedup( input_array, coords ):
input_array_precut_2x2x2 = numpy.zeros( (input_array.shape[0]-1, input_array.shape[1]-1, input_array.shape[2]-1, 8 ), dtype=DATA_DTYPE )
input_array_precut_2x2x2[ :, :, :, 0 ] = input_array[ 0:new_dimension-1, 0:new_dimension-1, 0:new_dimension-1 ]
input_array_precut_2x2x2[ :, :, :, 1 ] = input_array[ 1:new_dimension , 0:new_dimension-1, 0:new_dimension-1 ]
input_array_precut_2x2x2[ :, :, :, 2 ] = input_array[ 0:new_dimension-1, 1:new_dimension , 0:new_dimension-1 ]
input_array_precut_2x2x2[ :, :, :, 3 ] = input_array[ 0:new_dimension-1, 0:new_dimension-1, 1:new_dimension ]
input_array_precut_2x2x2[ :, :, :, 4 ] = input_array[ 1:new_dimension , 0:new_dimension-1, 1:new_dimension ]
input_array_precut_2x2x2[ :, :, :, 5 ] = input_array[ 0:new_dimension-1, 1:new_dimension , 1:new_dimension ]
input_array_precut_2x2x2[ :, :, :, 6 ] = input_array[ 1:new_dimension , 1:new_dimension , 0:new_dimension-1 ]
input_array_precut_2x2x2[ :, :, :, 7 ] = input_array[ 1:new_dimension , 1:new_dimension , 1:new_dimension ]
# adapted from from http://stackoverflow.com/questions/6427276/3d-interpolation-of-numpy-arrays-without-scipy
# 2012.03.02 - heavy modifications, to vectorise the final calculation... it is now superfast.
# - the checks are now removed in order to go faster...
# IMPORTANT: Input array is a pre-split, 8xNxMxO array.
# input coords could contain indexes at non-integer values (it's kind of the idea), whereas the coords_0 and coords_1 are integer values.
if coords.max() > min(input_array.shape[0:3])-1 or coords.min() < 0:
# do some checks to bring back the extremeties
# Could check each parameter in x y and z separately, but I know I get cubic data...
coords[numpy.where(coords>min(input_array.shape[0:3])-1)] = min(input_array.shape[0:3])-1
coords[numpy.where(coords<0 )] = 0
# for NxNxN data, coords[0].shape = N^3
output_array = numpy.zeros( coords[0].shape, dtype=DATA_DTYPE )
# a big array to hold all the coefficients for the trilinear interpolation
all_coeffs = numpy.zeros( (8,coords.shape[1]), dtype=DATA_DTYPE )
# the "floored" coordinates x, y, z
coords_0 = coords.astype(numpy.integer)
# all the above + 1 - these define the top left and bottom right (highest and lowest coordinates)
coords_1 = coords_0 + 1
# make the input coordinates "local"
coords = coords - coords_0
# Calculate one minus these values, in order to be able to do a one-shot calculation
# of the coefficients.
one_minus_coords = 1 - coords
# calculate those coefficients.
all_coeffs[0] = (one_minus_coords[0])*(one_minus_coords[1])*(one_minus_coords[2])
all_coeffs[1] = (coords[0]) *(one_minus_coords[1])*(one_minus_coords[2])
all_coeffs[2] = (one_minus_coords[0])* (coords[1]) *(one_minus_coords[2])
all_coeffs[3] = (one_minus_coords[0])*(one_minus_coords[1])* (coords[2])
all_coeffs[4] = (coords[0]) *(one_minus_coords[1])* (coords[2])
all_coeffs[5] = (one_minus_coords[0])* (coords[1]) * (coords[2])
all_coeffs[6] = (coords[0]) * (coords[1]) *(one_minus_coords[2])
all_coeffs[7] = (coords[0]) * (coords[1]) * (coords[2])
# multiply 8 greyscale values * 8 coefficients, and sum them across the "8 coefficients" direction
output_array = ( input_array[ coords_0[0], coords_0[1], coords_0[2] ].T * all_coeffs ).sum( axis=0 )
# and return it...
return output_array
我没有像上面那样分割x y和z坐标，因为之后重新合并它们似乎没有用。上面的代码中可能有一些假设立方数据（N = M = P），但我不这么认为......
让我知道你的想法！