Interpolacja 3D macierzy NumPy bez SciPy

Question

Piszę wtyczkę dla aplikacji, która zawiera NumPy w dystrybucji binarnej, ale nie SciPy. Moja wtyczka musi interpolować dane z jednej regularnej siatki 3D do innej regularnej siatki 3D. Uruchamiając go ze źródła, można to zrobić bardzo skutecznie, używając scipy.ndimage lub, jeśli użytkownik nie ma zainstalowanej SciPy, splot wygenerowany przez .pyd, który napisałem. Niestety, żadna z tych opcji nie jest dostępna, jeśli użytkownik uruchamia plik binarny.

Napisałem prostą rutynę trilinear interpolation w pythonie, która daje poprawny wynik, ale dla rozmiarów tablic, których używam, zajmuje dużo czasu (~ 5 minut). Zastanawiam się, czy istnieje sposób, aby przyspieszyć za pomocą tylko funkcji w NumPy. Podobnie jak scipy.ndimage.map_coordinates, pobiera ona tablicę wejściową 3D oraz tablicę ze współrzędnymi x, yiz dla każdego interpolowanego punktu.

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 

Oczywiście powodu funkcja jest tak powolny, to pętla for przez każdy punkt w przestrzeni 3D. Czy istnieje sposób na wykonanie magii krojenia lub wektoryzacji, aby przyspieszyć grę? Dzięki.

Answer 1

Okazuje się, że jest to w zawstydzeniu łatwe do wektoryzacji.

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

Dziękuję bardzo za ten wpis i za jego wykonanie. Opierałem się swobodnie o twoją wektoryzację, aby nadać jej kolejny wzrost prędkości (przynajmniej z danymi, z którymi pracuję)!

Pracuję z korelacją obrazu, a zatem I do interpolacji wielu zestawów o różnych współrzędnych w tym samym input_array.

Niestety, zrobiłem to trochę bardziej skomplikowanym, ale jeśli mogę wyjaśnić, co zrobiłem, dodatkowa komplikacja powinna a) usprawiedliwić się i b) stać się jasnym. Twoja ostatnia linia (output =) nadal wymaga sporego patrzenia w niesekwencyjne miejsca w input_array, przez co jest względnie powolna.

Załóżmy, że moje dane 3D mają długość NxMxP. Postanowiłem wykonać następujące czynności: Jeśli mogę uzyskać macierz (8 x (NxMxP)) obliczonych wstępnie wartości szarości dla punktu i jego najbliższych sąsiadów, a także mogę obliczyć macierz ((NxMxP) X 8) współczynniki (twój pierwszy współczynnik w powyższym przykładzie to (x-1) (y-1) (z-1)), wtedy mogę po prostu pomnożyć to razem i być w domu za darmo!

Dobrą premią dla mnie jest to, że mogę wstępnie obliczyć szare matryce i poddać je recyklingowi!

Oto nieco próbki kodu (wklejony z dwóch różnych funkcji, więc może nie działać po wyjęciu z pudełka, ale powinna służyć jako dobre źródło inspiracji):

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 

nie dzielonego Współrzędne xy i z, jak wyżej, ponieważ nie wydaje się użyteczne ich późniejsze odtworzenie. W powyższym kodzie może być coś, co zakłada dane sześcienne (N = M = P), ale nie sądzę ...

Daj mi znać, co myślisz!