使用NumPy进行3D插值，无需使用SciPy

Question

我正在编写一个插件，该应用程序在其二进制分发中包含NumPy，但不包括SciPy。我的插件需要将数据从一个常规的3D网格插值到另一个常规的3D网格。如果从源代码运行，则可以使用scipy.ndimage或者我编写的weave生成的.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

显然，这个函数变得如此缓慢的原因是在3D空间中循环遍历每个点的for循环。是否有任何方法可以执行某种切片或向量化魔法来加速它？谢谢。

- Stephen Terry

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

- Stephen Terry

Answer 2

非常感谢您的帖子，并为其跟进。我已经自由地使用了您的向量化，以便为其提供另一个速度提升（至少对于我正在处理的数据而言）！ 我正在处理图像相关性，因此我需要在相同的输入数组中插值多组不同的坐标。 不幸的是，我让它变得有点复杂，但如果我能解释我所做的事情，额外的复杂性应该会证明其合理性并变得清晰。您的最后一行（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），但我不这么认为...

你觉得呢？

- EdwardAndo