interpolieren 3D-Volumen mit numpy und oder scipy

Question

Ich bin extrem frustriert, weil ich nach mehreren Stunden nicht in der Lage, eine scheinbar einfache 3D-Interpolation in Python zu tun scheinen kann. In Matlab alles, was ich tun musste, war

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

Was ist das genaue Äquivalent dieser ndimage.map_coordinate der mit scipy oder anderen numpy Methoden?

Dank

Answer 1

In scipy 0,14 oder höher, gibt es eine neue Funktion scipy.interpolate.RegularGridInterpolator das ähnelt interp3.

Der Befehl MATLAB Vi = interp3(x,y,z,V,xi,yi,zi) wie etwas übersetzen würde:

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) 

Hier ist ein vollständiges Beispiel demonstriert beide; es wird Ihnen die genauen Unterschiede ...

MATLAB-Code verstehen helfen:

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); 

Das Ergebnis ist Vi=[268 357], die in der Tat an diesen beiden Punkten (2,6,8) und (3,5,7) der Wert ist.

SciPy Code:

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)) 

Wieder ist es [268,357]. Sie sehen also einige kleine Unterschiede: Scipy verwendet x, y, z Indexreihenfolge, während MATLAB y, x, z (seltsamerweise) verwendet; In Scipy definieren Sie eine Funktion in einem separaten Schritt und wenn Sie sie aufrufen, werden die Koordinaten gruppiert wie (x1, y1, z1), (x2, y2, z2), ... während Matlab (x1, x2, ..) verwendet. .), (y1, y2, ...), (z1, z2, ...).

Ansonsten sind die beiden ähnlich und ebenso einfach zu bedienen.

Answer 2

Grundsätzlich ndimage.map_coordinates Arbeiten in "Index" Koordinaten (auch bekannt als "Voxel" oder "Pixel" Koordinaten). Die Schnittstelle zu ihr scheint zunächst ein wenig klobig, aber es gibt Ihnen eine Los von Flexibilität.

Wenn Sie die interpolierten Koordinaten ähnlich wie Matlab interp3 angeben möchten, dann müssen Sie Ihre eingegebenen Koordinaten in "Index" -Koordinaten konvertieren.

Es gibt auch die zusätzliche Falte, map_coordinates behält immer den dtype des Eingabearrays in der Ausgabe. Wenn Sie ein Integer-Array interpolieren, erhalten Sie eine Integer-Ausgabe, die möglicherweise nicht Ihren Vorstellungen entspricht. Für das folgende Codefragment nehme ich an, dass Sie immer Gleitkommaausgabe möchten. (Wenn nicht, ist es eigentlich einfacher.)

Ich werde versuchen, später noch mehr Erklärung hinzufügen (das ist ziemlich dichten Code).

Alles in allem ist die interp3 Funktion, die ich habe, komplexer, als es für Ihre genauen Zwecke sein muss. Howver, sollte das Verhalten von interp3 mehr oder weniger replizieren, wie ich es erinnere (ohne Berücksichtigung der „Zoom“ -Funktion von interp3(data, zoom_factor), die scipy.ndimage.zoom Griffe.)

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

Die genaue entspricht MATLAB interp3'S wäre mit scipy interpn für einmalige Interpolation:

import numpy as np 
from scipy.interpolate import interpn 

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

Die Standardmethode sowohl für MATLAB und scipy ist eine lineare Interpolation, und dies mit dem method geändert werden kann Streit. Beachten Sie, dass nur lineare und nächste Nachbar-Interpolation von interpn für 3 Dimensionen und darüber unterstützt wird, im Gegensatz zu MATLAB, das ebenfalls kubische und Spline-Interpolation unterstützt.

Wenn mehrere Interpolationsaufrufe auf dem gleichen Raster ausgeführt werden, ist es vorzuziehen, das Interpolationsobjekt RegularGridInterpolator zu verwenden, wie in der angenommenen Antwort above. interpn verwendet RegularGridInterpolator intern.