3D-интерполяция массивов NumPy без SciPy

Question

Я пишу плагин для приложения, которое включает NumPy в бинарный дистрибутив, но не SciPy. Мой плагин должен интерполировать данные из одной обычной трехмерной сетки в другую обычную трехмерную сетку. Работая из исходного кода, это можно сделать очень эффективно, используя scipy.ndimage или, если у пользователя не установлен SciPy, сгенерировано переплетение .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 над каждой точкой в ​​трехмерном пространстве. Есть ли способ выполнить какую-то магию нарезки или векторизации, чтобы ускорить ее? Спасибо.

Answer 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

Answer 2

Большое спасибо за этот пост и за то, что следил за ним.Я основательно основываюсь на вашей векторизации, чтобы повысить скорость (по крайней мере, с данными, с которыми я работаю)!

Я работаю с корреляцией изображений, и поэтому я интерполирую множество наборовразличных координат в одном и том же input_array .

К сожалению, я сделал это немного сложнее, но если я смогу объяснить, что я сделал, дополнительное усложнение должно а) оправдать себя и b) стало ясно.Ваша последняя строка ( output = ) по-прежнему требует значительного количества поиска в непоследовательных местах в 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 /5259470/3d-interpolyatsiya-massivov-numpy-bez-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

Я не разделял координаты xy и z, как указано выше, потому что не представляется полезным их впоследствии объединять.В приведенном выше коде может быть что-то, что предполагает кубические данные (N = M = P), но я так не думаю ...

Дайте мне знать, что вы думаете!