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

Question

Я очень расстроен, потому что через несколько часов я не могу сделать, казалось бы, легкую трехмерную интерполяцию в python. В Matlab все, что я должен был сделать

Vi = interp3(x,y,z,V,xi,yi,zi) 

Что является точным эквивалентом ndimage.map_coordinate это с помощью SciPy или других Numpy методов?

Благодаря

Answer 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, ...).

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

Answer 2

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

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

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

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

В целом, функция interp3 У меня есть более сложная задача, чем это возможно для ваших конкретных целей. Howver, он должен более или менее копирует поведение 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() 

Answer 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 для 3-х измерений и выше, в отличие от MATLAB, которая поддерживает также кубическую и сплайн-интерполяцию.

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