非常感谢您的帖子,并为其跟进。我已经自由地使用了您的向量化,以便为其提供另一个速度提升(至少对于我正在处理的数据而言)!
我正在处理图像相关性,因此我需要在相同的输入数组中插值多组不同的坐标。
不幸的是,我让它变得有点复杂,但如果我能解释我所做的事情,额外的复杂性应该会证明其合理性并变得清晰。您的最后一行(output =)仍然需要在输入数组的非连续位置进行大量查找,这使得它相对较慢。
假设我的三维数据长为NxMxP。我决定做以下事情:如果我可以获得一个(8 x(NxMxP))的预计算灰度值点及其最近邻居的矩阵,并且我也可以计算一个((NxMxP)X8)的系数矩阵(您上面示例中的第一个系数为(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 https://dev59.com/7Gw15IYBdhLWcg3wv-NM
# 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),但我不这么认为...
你觉得呢?