SciPyを使用しないNumPy配列の3D補間

Question

SciPyではなく、バイナリディストリビューションにNumPyを含むアプリケーション用のプラグインを作成しています。私のプラグインは、1つの通常の3Dグリッドから別の通常の3Dグリッドにデータを補間する必要があります。ソースから実行すると、これはscipy.ndimageを使用して非常に効率的に行うことができます。ユーザーがSciPyをインストールしていない場合は、.pydと書かれています。残念ながら、ユーザーがバイナリを実行している場合、これらのオプションはどちらも使用できません。

私は正しい結果を与える単純なtrilinear interpolationルーチンを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空間内の各点でループすることです。スピードアップのためにスライシングやベクトル化の魔法を実行する方法はありますか？ありがとう。

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。

残念ながら私はもう少し複雑にしましたが、私が行ったことを説明できる場合、余分な合併症はa）それ自体を正当化し、b）明確になるべきです。あなたの最後の行（出力=）は、input_arrayの非連続的な場所でかなりの量のルックアップを必要とし、比較的遅くなります。

私の3DデータがNxMxP長いとします。私は次のことをすることにしました：ポイントとそれに最も近い隣人のための事前計算されたグレイレベルの（8x（NxMxP））マトリックスを得ることができ、また（（NxMxP）X 8） （上の例の最初の係数は（x-1）（y-1）（z-1）です）、次に私はちょうど一緒に乗算して家に帰ることができます！

私にとって良いボーナスは、グレーマトリックスをあらかじめ計算してリサイクルできることです。ここで

は、コードのサンプルビットである（2つの異なる機能から貼り付け、その箱から出して動作しない場合がありますが、インスピレーションの良い情報源として役立つはずです）：

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 

私は分割しませんでしたxyとz座標は上記のように後でそれらを修復するために有用ではないように思われたからです。上記のコードにキュービックデータ（N = M = P）を仮定したものがあるかもしれませんが、そうは思わないでしょう...

私はあなたの考えをお知らせください！