Interpolation rapide sur un tableau 3D

j'ai un tableau 3D que je dois interpoler sur un axe (la dernière dimension). Disons y.shape = (nx, ny, nz), je tiens à interpoler dans nz pour tout (nx, ny). Cependant, je veux interpoler pour une valeur différente dans chaque [i, j].

voici quelques exemples de code. Si je voulais interpoler à une seule valeur, dire new_z, je voudrais utiliser scipy.interpolate.interp1d comme ceci

# y is a 3D ndarray
# x is a 1D ndarray with the abcissa values
# new_z is a number
f = scipy.interpolate.interp1d(x, y, axis=-1, kind='linear')
result = f(new_z)

Cependant, pour ce problème ce que je veux réellement c'est interpoler à un différent new_z pour y[i, j]. Donc ce que je fais:

# y is a 3D ndarray
# x is a 1D ndarray with the abcissa values
# new_z is a 2D array
result = numpy.empty(y.shape[:-1])
for i in range(nx):
    for j in range(ny):
        f = scipy.interpolate.interp1d(x, y[i, j], axis=-1, kind='linear')
        result[i, j] = f(new_z[i, j])

malheureusement, avec plusieurs boucles, cela devient inefficace et lent. Est-il une meilleure façon de faire ce type d'interpolation? L'interpolation linéaire est suffisante. Une possibilité est de mettre en Cython, mais j'essayais d'éviter que parce que je veux avoir la possibilité de changer d'interpolation cubique et ne veulent pas le faire à la main en Cython.

10
demandé sur tiago 2012-11-21 11:53:21

6 réponses

Pour l'accélération élevée afin d'interpoler, vous pouvez appeler interp1d() une seule fois, puis utilisez le _spline attribut et la fonction de bas niveau _bspleval() dans le _fitpack module. Voici le code:

from scipy.interpolate import interp1d
import numpy as np

nx, ny, nz = 30, 40, 50
x = np.arange(0, nz, 1.0)
y = np.random.randn(nx, ny, nz)
new_x = np.random.random_integers(1, (nz-1)*10, size=(nx, ny))/10.0

def original_interpolation(x, y, new_x):
    result = np.empty(y.shape[:-1])
    for i in xrange(nx):
        for j in xrange(ny):
            f = interp1d(x, y[i, j], axis=-1, kind=3)
            result[i, j] = f(new_x[i, j])
    return result

def fast_interpolation(x, y, new_x):
    from scipy.interpolate._fitpack import _bspleval
    f = interp1d(x, y, axis=-1, kind=3)
    xj,cvals,k = f._spline
    result = np.empty_like(new_x)
    for (i, j), value in np.ndenumerate(new_x):
        result[i, j] = _bspleval(value, x, cvals[:, i, j], k, 0)
    return result

r1 = original_interpolation(x, y, new_x)
r2 = fast_interpolation(x, y, new_x)

>>> np.allclose(r1, r2)
True

%timeit original_interpolation(x, y, new_x)
%timeit fast_interpolation(x, y, new_x)
1 loops, best of 3: 3.78 s per loop
100 loops, best of 3: 15.4 ms per loop
6
répondu HYRY 2017-06-04 15:12:39

comme la suggestion de numpy ci-dessus était trop longue, je pourrais attendre donc voici la version de cython pour référence future. De quelques repères lâches il est environ 3000 fois plus rapide (certes, il est seulement interpolation linéaire et ne fait pas autant que interp1d mais c'est ok pour cette raison).

import numpy as N
cimport numpy as N
cimport cython

DTYPEf = N.float64
ctypedef N.float64_t DTYPEf_t

@cython.boundscheck(False) # turn of bounds-checking for entire function
@cython.wraparound(False)  # turn of bounds-checking for entire function
cpdef interp3d(N.ndarray[DTYPEf_t, ndim=1] x, N.ndarray[DTYPEf_t, ndim=3] y,
               N.ndarray[DTYPEf_t, ndim=2] new_x):
    """
    interp3d(x, y, new_x)

    Performs linear interpolation over the last dimension of a 3D array,
    according to new values from a 2D array new_x. Thus, interpolate
    y[i, j, :] for new_x[i, j].

    Parameters
    ----------
    x : 1-D ndarray (double type)
        Array containg the x (abcissa) values. Must be monotonically
        increasing.
    y : 3-D ndarray (double type)
        Array containing the y values to interpolate.
    x_new: 2-D ndarray (double type)
        Array with new abcissas to interpolate.

    Returns
    -------
    new_y : 3-D ndarray
        Interpolated values.
    """
    cdef int nx = y.shape[0]
    cdef int ny = y.shape[1]
    cdef int nz = y.shape[2]
    cdef int i, j, k
    cdef N.ndarray[DTYPEf_t, ndim=2] new_y = N.zeros((nx, ny), dtype=DTYPEf)

    for i in range(nx):
        for j in range(ny):
            for k in range(1, nz):
                 if x[k] > new_x[i, j]:
                     new_y[i, j] = (y[i, j, k] - y[i, j, k - 1]) * \
                  (new_x[i, j] - x[k-1]) / (x[k] - x[k - 1]) + y[i, j, k - 1]
                     break
    return new_y
3
répondu tiago 2012-11-21 13:51:28

je ne pense pas interp1d a une méthode pour faire cela rapidement, donc vous ne pouvez pas éviter la boucle ici.

Cython vous pouvez probablement encore éviter en codant l'interpolation linéaire en utilisant np.searchsorted, quelque chose de ce genre (pas testé):

def interp3d(x, y, new_x):
    assert x.ndim == 1 and y.ndim == 3 and new_x.ndim == 2
    assert y.shape[:2] == new_x.shape and x.shape == y.shape[2:]

    nx, ny = y.shape[:2]
    new_x = new_x.ravel()
    j = np.arange(len(new_x))
    k = np.searchsorted(x, new_x).clip(1, len(x) - 1)
    y = y.reshape(-1, x.shape[0])
    p = (new_x - x[k-1]) / (x[k] - x[k-1])
    result = (1 - p) * y[j,k-1] + p * y[j,k]
    return result.reshape(nx, ny)

N'aide pas avec interpolation cubique.

EDIT: en a fait une fonction et a corrigé les erreurs de type "off-by-one". Certains timings vs Cython (500x500x500 de la grille):

In [58]: %timeit interp3d(x, y, new_x)
10 loops, best of 3: 82.7 ms per loop

In [59]: %timeit cyfile.interp3d(x, y, new_x)
10 loops, best of 3: 86.3 ms per loop

In [60]: abs(interp3d(x, y, new_x) - cyfile.interp3d(x, y, new_x)).max()
Out[60]: 2.2204460492503131e-16

cependant, on peut argumenter que le code Cython est plus facile à lire.

3
répondu pv. 2012-11-21 22:21:37

@pv.'s réponse, et la vectorisation de la boucle interne, ce qui suit donne une accélération substantielle (EDIT: changé le numpy.tile à l'aide de numpy.lib.stride_tricks.as_strided):

import numpy
from scipy import interpolate

nx = 30
ny = 40
nz = 50

y = numpy.random.randn(nx, ny, nz)
x = numpy.float64(numpy.arange(0, nz))

# We select some locations in the range [0.1, nz-0.1]
new_z = numpy.random.random_integers(1, (nz-1)*10, size=(nx, ny))/10.0

# y is a 3D ndarray
# x is a 1D ndarray with the abcissa values
# new_z is a 2D array

def original_interpolation():
    result = numpy.empty(y.shape[:-1])
    for i in range(nx):
        for j in range(ny):
            f = interpolate.interp1d(x, y[i, j], axis=-1, kind='linear')
            result[i, j] = f(new_z[i, j])

    return result

grid_x, grid_y = numpy.mgrid[0:nx, 0:ny]
def faster_interpolation():
    flat_new_z = new_z.ravel()
    k = numpy.searchsorted(x, flat_new_z)
    k = k.reshape(nx, ny)

    lower_index = [grid_x, grid_y, k-1]
    upper_index = [grid_x, grid_y, k]

    tiled_x = numpy.lib.stride_tricks.as_strided(x, shape=(nx, ny, nz), 
        strides=(0, 0, x.itemsize))

    z_upper = tiled_x[upper_index]
    z_lower = tiled_x[lower_index]

    z_step = z_upper - z_lower
    z_delta = new_z - z_lower

    y_lower = y[lower_index]
    result = y_lower + z_delta * (y[upper_index] - y_lower)/z_step

    return result

# both should be the same (giving a small difference)
print numpy.max(
        numpy.abs(original_interpolation() - faster_interpolation()))

Que donne le temps sur mon ordinateur:

In [8]: timeit foo.original_interpolation()
10 loops, best of 3: 102 ms per loop

In [9]: timeit foo.faster_interpolation()
1000 loops, best of 3: 564 us per loop

Va nx = 300,ny = 300 et nz = 500, donne un 130 x speedup:

In [2]: timeit original_interpolation()
1 loops, best of 3: 8.27 s per loop

In [3]: timeit faster_interpolation()
10 loops, best of 3: 60.1 ms per loop

vous auriez besoin d'écrire votre propre algorithme pour l'interpolation cubique, mais cela ne devrait pas être si difficile.

2
répondu Henry Gomersall 2017-05-23 11:44:31

bien Qu'il y ait plusieurs bonnes réponses, ils font encore des interpolations de 250k dans un réseau fixe de 500 long:

j250k = np.searchsorted( X500, X250k )  # indices in [0, 500)

cela peut être accéléré avec un LUT, une table de recherche, avec par exemple des fentes de 5k:

lut = np.interp( np.arange(5000), X500, np.arange(500) ).round().astype(int)
xscale = (X - X.min()) * (5000 - 1) \
        / (X.max() - X.min()) 
j = lut.take( xscale.astype(int), mode="clip" )  # take(floats) in numpy 1.7 ?

#---------------------------------------------------------------------------
# X     |    |       | |             |
# j     0    1       2 3             4 ...
# LUT   |....|.......|.|.............|....  -> int j (+ offset in [0, 1) )
#---------------------------------------------------------------------------

searchsorted est assez rapide, temps ~ ln2 500, donc c'est probablement pas beaucoup plus vite.

Mais les LUTs sont très rapide en C, un simple vitesse de la mémoire et du compromis.

2
répondu denis 2013-02-08 10:16:35

Vous pouvez utiliser map_coordinates pour:

from numpy import random, meshgrid, arange
from scipy.ndimage import map_coordinates

(nx, ny, nz) = (4, 5, 6)
# some random array
A = random.rand(nx, ny, nz)

# random floating-point indices in [0, nz-1]
Z = random.rand(nx, ny)*(nz-1)

# regular integer indices of shape (nx,ny)
X, Y = meshgrid(arange(nx), arange(ny), indexing='ij')

coords = (X, Y, Z) # X, Y, and Z are of shape (nx, ny)

print map_coordinates(A, coords, order=1, cval=-999.)
1
répondu François 2012-11-22 10:10:56