интерполировать трехмерный объем с помощью numpy и scipy

Ответ 1

В scipy 0.14 или новее есть новая функция scipy.interpolate.RegularGridInterpolator которая очень похожа на interp3.

Команда MATLAB Vi = interp3(x,y,z,V,xi,yi,zi) переведет нечто вроде:

from numpy import array
from scipy.interpolate import RegularGridInterpolator as rgi
my_interpolating_function = rgi((x,y,z), V)
Vi = my_interpolating_function(array([xi,yi,zi]).T)

Вот полный пример, демонстрирующий оба; это поможет вам понять точные различия...

КОД MATLAB:

x = linspace(1,4,11);
y = linspace(4,7,22);
z = linspace(7,9,33);
V = zeros(22,11,33);
for i=1:11
    for j=1:22
        for k=1:33
            V(j,i,k) = 100*x(i) + 10*y(j) + z(k);
        end
    end
end
xq = [2,3];
yq = [6,5];
zq = [8,7];
Vi = interp3(x,y,z,V,xq,yq,zq);

В результате получается Vi=[268 357] что действительно является значением в этих двух точках (2,6,8) и (3,5,7).

КОД SCIPY:

from scipy.interpolate import RegularGridInterpolator
from numpy import linspace, zeros, array
x = linspace(1,4,11)
y = linspace(4,7,22)
z = linspace(7,9,33)
V = zeros((11,22,33))
for i in range(11):
    for j in range(22):
        for k in range(33):
            V[i,j,k] = 100*x[i] + 10*y[j] + z[k]
fn = RegularGridInterpolator((x,y,z), V)
pts = array([[2,6,8],[3,5,7]])
print(fn(pts))

Снова это [268,357]. Итак, вы видите некоторые незначительные отличия: Scipy использует x, y, z индексный порядок, в то время как MATLAB использует y, x, z (странно); В Scipy вы определяете функцию на отдельном шаге, и когда вы ее вызываете, координаты группируются как (x1, y1, z1), (x2, y2, z2),... тогда как в matlab используются (x1, x2,....), (у1, у2,...), (z1, z2,...).

Помимо этого, они похожи и одинаково просты в использовании.

Ответ 2

В принципе, ndimage.map_coordinates работает в координатах "индекс" (так называемые "воксельные" или "пиксельные" координаты). Интерфейс к нему кажется немного неуклюжим, но это дает вам большую гибкость.

Если вы хотите указать интерполированные координаты, похожие на matlab interp3, тогда вам нужно будет преобразовать ваши intput-координаты в "индексные" координаты.

Там также дополнительная морщина, которую map_coordinates всегда сохраняет dtype входного массива в выходе. Если вы интерполируете целочисленный массив, вы получите целочисленный вывод, который может быть или не быть тем, что вы хотите. Для нижеприведенного фрагмента кода я предполагаю, что вам всегда нужен вывод с плавающей запятой. (Если вы этого не сделаете, это будет проще.)

Сегодня вечером я попытаюсь добавить больше объяснений (это довольно плотный код).

В целом, функция interp3 я имею, более сложна, чем может быть, для ваших конкретных целей. Как бы то ни было, он должен более или менее воспроизводить поведение interp3 как я его помню (игнорируя "масштабирование" функций interp3(data, zoom_factor), которые обрабатывает scipy.ndimage.zoom.)

import numpy as np
from scipy.ndimage import map_coordinates

def main():
    data = np.arange(5*4*3).reshape(5,4,3)

    x = np.linspace(5, 10, data.shape[0])
    y = np.linspace(10, 20, data.shape[1])
    z = np.linspace(-100, 0, data.shape[2])

    # Interpolate at a single point
    print interp3(x, y, z, data, 7.5, 13.2, -27)

    # Interpolate a region of the x-y plane at z=-25
    xi, yi = np.mgrid[6:8:10j, 13:18:10j]
    print interp3(x, y, z, data, xi, yi, -25 * np.ones_like(xi))

def interp3(x, y, z, v, xi, yi, zi, **kwargs):
    """Sample a 3D array "v" with pixel corner locations at "x","y","z" at the
    points in "xi", "yi", "zi" using linear interpolation. Additional kwargs
    are passed on to ''scipy.ndimage.map_coordinates''."""
    def index_coords(corner_locs, interp_locs):
        index = np.arange(len(corner_locs))
        if np.all(np.diff(corner_locs) < 0):
            corner_locs, index = corner_locs[::-1], index[::-1]
        return np.interp(interp_locs, corner_locs, index)

    orig_shape = np.asarray(xi).shape
    xi, yi, zi = np.atleast_1d(xi, yi, zi)
    for arr in [xi, yi, zi]:
        arr.shape = -1

    output = np.empty(xi.shape, dtype=float)
    coords = [index_coords(*item) for item in zip([x, y, z], [xi, yi, zi])]

    map_coordinates(v, coords, order=1, output=output, **kwargs)

    return output.reshape(orig_shape)

main()

Ответ 3

Точный эквивалент MATLAB interp3 будет использовать scipy interpn для однократной интерполяции:

import numpy as np
from scipy.interpolate import interpn

Vi = interpn((x,y,z), V, np.array([xi,yi,zi]).T)

Метод по умолчанию для MATLAB и scipy - это линейная интерполяция, и это можно изменить с помощью аргумента method. Обратите внимание, что только интерполяция с interpn линейных и ближайших соседей поддерживается методом interpn для 3-х измерений и выше, в отличие от MATLAB, который поддерживает также кубическую и сплайн-интерполяцию.

При выполнении нескольких вызовов интерполяции в одной и той же сетке предпочтительнее использовать объект интерполяции RegularGridInterpolator, как в принятом ответе выше. interpn использует RegularGridInterpolator внутренне.