Nội suy 3D của các mảng NumPy mà không cần SciPy

Question

Tôi đang viết một plugin cho một ứng dụng bao gồm NumPy trong phân phối nhị phân, nhưng không phải là SciPy. Plugin của tôi cần phải nội suy dữ liệu từ một lưới 3D thông thường sang một lưới 3D thông thường khác. Chạy từ nguồn, điều này có thể được thực hiện rất hiệu quả bằng cách sử dụng scipy.ndimage hoặc, nếu người dùng không có SciPy được cài đặt, một dệt tạo ra .pyd Tôi đã viết. Thật không may, không phải của các tùy chọn có sẵn nếu người dùng đang chạy nhị phân.

Tôi đã viết một thói quen đơn giản trilinear interpolation trong python cho kết quả chính xác, nhưng đối với kích thước mảng tôi đang sử dụng, mất một thời gian dài (~ 5 phút). Tôi tự hỏi nếu có một cách để tăng tốc độ nó lên bằng cách sử dụng chỉ các chức năng trong NumPy. Cũng giống như scipy.ndimage.map_coordinates, Phải mất một mảng đầu vào 3D và một mảng có toạ độ x, y và z của mỗi điểm được nội suy.

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 

Rõ ràng lý do chức năng quá chậm trên mỗi điểm trong không gian 3D. Có cách nào để thực hiện một số loại cắt hoặc ma thuật vector hóa để tăng tốc nó? Cảm ơn.

Answer 1

Hóa ra thật dễ dàng để làm mờ nó.

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

Cảm ơn rất nhiều vì bài đăng này và đã theo dõi nó. Tôi đã dựa trên bản thân mình một cách tự do trên vectơ của bạn để tăng thêm tốc độ (ít nhất là với dữ liệu tôi đang làm việc)!

Tôi đang làm việc với tương quan hình ảnh và do đó tôi nội suy nhiều tập tọa độ khác nhau trong cùng một input_array.

Thật không may tôi đã làm cho nó phức tạp hơn một chút, nhưng nếu tôi có thể giải thích những gì tôi đã làm, các biến chứng thêm nên a) biện minh cho chính nó và b) trở nên rõ ràng. Dòng cuối cùng của bạn (đầu ra =) vẫn yêu cầu số tiền hợp lý để tra cứu ở những vị trí không tuần tự trong input_array, khiến cho nó tương đối chậm.

Giả sử dữ liệu 3D của tôi dài NxMxP. Tôi đã quyết định làm điều sau đây: Nếu tôi có thể nhận được ma trận (8 x (NxMxP)) của các giá trị được tính trước cho một điểm và các láng giềng gần nhất, và tôi cũng có thể tính toán ma trận ((NxMxP) X 8) hệ số (hệ số đầu tiên của bạn trong ví dụ trên là (x-1) (y-1) (z-1)) sau đó tôi chỉ có thể nhân với nhau và trở về nhà miễn phí!

Điểm thưởng tốt cho tôi là tôi có thể tính toán trước ma trận màu xám và tái chế nó!

Dưới đây là một chút mẫu mã (dán từ hai chức năng khác nhau, vì vậy có thể không làm việc ra khỏi hộp, nhưng phải phục vụ như là một nguồn tốt của cảm hứng):

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 

tôi không chia xy và z tọa độ như trên bởi vì nó không có vẻ hữu ích để remerge chúng sau đó. Có thể có điều gì đó trong đoạn mã trên giả định dữ liệu khối (N = M = P) nhưng tôi không nghĩ vậy ...

Hãy cho tôi biết suy nghĩ của bạn!