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
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
3D interpolation of NumPy arrays without SciPy