3D-интерполяция массивов NumPy без SciPy
Я пишу плагин для приложения, которое включает NumPy в двоичном дистрибутиве, но не SciPy. Мой плагин должен интерполировать данные из одной регулярной трехмерной сетки в другую обычную трехмерную сетку. Исходя из источника, это можно сделать очень эффективно с помощью scipy.ndimage или, если у пользователя нет установленного SciPy, сгенерировано weave .pyd. К сожалению, ни одна из этих опций недоступна, если пользователь запускает двоичный файл.
Я написал простую процедуру трилинейную интерполяцию в python, которая дает правильный результат, но для размеров массивов, которые я использую, принимает долгое время (~ 5 минут). Мне интересно, есть ли способ ускорить его, используя только функциональность в NumPy. Точно так же, как scipy.ndimage.map_coordinates, он принимает трехмерный входной массив и массив с координатами 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 для каждой точки в трехмерном пространстве. Есть ли способ выполнить какую-то магию для нарезки или векторизации, чтобы ускорить ее? Спасибо.
Ответы
Ответ 1
Получается, что смущающе легко его векторизовать.
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
Ответ 2
Большое спасибо за этот пост и за то, что он это сделал. Я основывал себя на вашей векторности, чтобы дать ей еще одно ускорение скорости (по крайней мере, с данными, с которыми я работаю)!
Я работаю с корреляцией изображения, и поэтому я должен интерполировать множество наборов разных координат в одном и том же input_array.
К сожалению, я сделал это немного сложнее, но если я могу объяснить, что я сделал, дополнительное усложнение должно: a) оправдать себя и b) стать ясным. Ваша последняя строка (output =) по-прежнему требует значительного поиска в непоследовательных местах в input_array, что делает ее относительно медленной.
Предположим, что мои 3D-данные длинны NxMxP. Я решил сделать следующее: если я могу получить (8 x (NxMxP)) матрицу предварительно вычисленных greyvalues для точки и ее ближайших соседей, и я также могу вычислить ((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/info/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 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), но я так не думаю...
Дайте мне знать, что вы думаете!
