Трехмерная интерполяция массивов NumPy без SciPy

Я пишу плагин для приложения, которое включает NumPy в двоичный дистрибутив, но не SciPy. Мой плагин должен интерполировать данные из одной обычной 3D-сетки в другую обычную 3D-сетку. При запуске из исходников это можно сделать очень эффективно, используя scipy.ndimage или, если у пользователя не установлен SciPy, сгенерированный мной файл weave .pyd . К сожалению, ни один из этих вариантов недоступен, если пользователь запускает двоичный файл.

Я написал простую процедуру трилинейной интерполяции на Python, которая дает правильный результат, но для размеров массива я использование занимает много времени (~ 5 минут). Мне интересно, есть ли способ ускорить его, используя только функциональность NumPy. Как и scipy.ndimage.map_coordinates , для интерполяции требуется трехмерный входной массив и массив с координатами 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

Очевидно, причина такой медленной работы функции в цикл для по каждой точке в трехмерном пространстве. Есть ли способ выполнить какую-то магию нарезки или векторизации, чтобы ускорить это? Спасибо.