interpoler le volume 3D avec numpy et ou scipy

Je suis extrêmement frustré parce que, après plusieurs heures, je n'arrive pas à être en mesure de faire un semblant de 3D facile d'interpolation en python. Dans Matlab tout ce que j'avais à faire était de

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

Ce qui est l'équivalent exact de cette aide scipy est ndimage.map_coordinate ou d'autres numpy méthodes?

Grâce

source d'informationauteur user1301295 | 2014-02-17

  1. 18

    Dans scipy 0.14 ou tard, il y est une nouvelle fonction scipy.interpoler.RegularGridInterpolator qui ressemble étroitement à interp3.

    L'MATLAB commande Vi = interp3(x,y,z,V,xi,yi,zi) serait traduit en quelque chose comme:

    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)

    Voici un exemple complet démontrant à la fois; il va vous aider à comprendre les différences exactes...

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

    Le résultat est Vi=[268 357] qui est en effet la valeur de ces deux points (2,6,8) et (3,5,7).

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

    De nouveau, c'est [268,357]. Donc, vous voyez quelques différences: Scipy utilise x,y,z ordre d'index alors que MATLAB utilise y,x,z (étrangement); Dans Scipy vous de définir une fonction dans une étape distincte et quand vous l'appelez, les coordonnées sont regroupés comme (x1,y1,z1),(x2,y2,z2),... alors que matlab utilise (x1,x2,...),(y1,y2,...),(z1,z2,...).

    Autre que cela, les deux sont similaires et tout aussi facile à utiliser.

  2. 2

    Fondamentalement, ndimage.map_coordinates fonctionne en "index" de coordonnées (un.k.un. "voxel" ou "pixel" coordonnées). L'interface semble un peu maladroit au début, mais il ne vous donner un beaucoup de flexibilité.

    Si vous souhaitez spécifier l'interpolation des coordonnées similaire à matlab interp3alors vous aurez besoin de convertir vos entrées coordonnées dans "index" coordonnées.

    Il y a aussi les rides que map_coordinates conserve toujours la dtype du tableau d'entrée à la sortie. Si vous interpoler un tableau d'entiers, vous obtiendrez un nombre entier de sortie, qui peut ou peut ne pas être ce que vous voulez. Pour l'extrait de code ci-dessous, je vais supposer que vous en voulez toujours flottant point de sortie. (Si vous n'en avez pas, c'est effectivement plus simple.)

    Je vais essayer d'ajouter plus d'explications plus tard ce soir (ce qui est assez dense code).

    Dans l'ensemble, le interp3 fonction de ce que j'ai est plus complexe qu'il faudra peut-être pour exactement à vos fins. Howver, il devrait plus ou moins reproduire le comportement de interp3 que je m'en souviens (en ignorant le "zoom" de la fonctionnalité de interp3(data, zoom_factor)qui scipy.ndimage.zoom poignées.) 

    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. 2

    La exacte équivalent à MATLAB interp3 serait à l'aide de scipy est interpn pour un hors d'interpolation:

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

    La méthode par défaut pour les deux MATLAB et scipy est l'interpolation linéaire, et cela peut être modifié avec la method argument. Notez que seulement linéaire et le plus proche voisin de l'interpolation est pris en charge par interpn pour les 3 dimensions et au-dessus, à la différence de MATLAB qui prend en charge cubes et interpolation spline.

    Lors de la prise de plusieurs interpolation des appels sur la même grille, il est préférable d'utiliser l'interpolation objet RegularGridInterpolatorcomme dans la accepté de répondre à ci-dessus. interpn utilise RegularGridInterpolator en interne.