无需SciPy的NumPy数组3D插值

Question

我正在为一个应用程序编写插件，这个应用程序的二进制版本里包含了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

显然，函数运行这么慢的原因是因为在3D空间中对每个点进行的for循环。有没有什么办法可以通过切片或向量化的技巧来加速这个过程呢？谢谢。

Answer 1

非常感谢你发的这个帖子，也感谢你后续的跟进。我在你的向量化基础上进行了自由发挥，给它加了一点速度（至少在我处理的数据上是这样）！

我正在做图像相关性分析，因此需要在同一个input_array中插值许多不同坐标的集合。

不幸的是，我把事情弄得有点复杂，不过如果我能解释我做的事情，这些额外的复杂性应该会自我证明并且变得清晰。你最后一行的output =仍然需要在input_array中查找不连续的地方，这样会比较慢。

假设我的三维数据是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），但我觉得不是……

告诉我你的想法！

Answer 2

结果发现，把它变成向量形式其实非常简单。

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