I am extremely frustrated because after several hours I can't seem to be able to do a seemingly easy 3D interpolation in python. In Matlab all I had to do was
Vi = interp3(x,y,z,V,xi,yi,zi)
What is the exact equivalent of this using scipy's ndimage.map_coordinate or other numpy methods?
In scipy 0.14 or later, there is a new function
scipy.interpolate.RegularGridInterpolator which closely resembles
The MATLAB command
Vi = interp3(x,y,z,V,xi,yi,zi) would translate to something like:
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)
Here is a full example demonstrating both; it will help you understand the exact differences...
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);
The result is
Vi=[268 357] which is indeed the value at those two points
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]. So you see some slight differences: Scipy uses x,y,z index order while MATLAB uses y,x,z (strangely); In Scipy you define a function in a separate step and when you call it, the coordinates are grouped like (x1,y1,z1),(x2,y2,z2),... while matlab uses (x1,x2,...),(y1,y2,...),(z1,z2,...).
Other than that, the two are similar and equally easy to use.
ndimage.map_coordinates works in "index" coordinates (a.k.a. "voxel" or "pixel" coordinates). The interface to it seems a bit clunky at first, but it does give you a lot of flexibility.
If you want to specify the interpolated coordinates similar to matlab's
interp3, then you'll need to convert your intput coordinates into "index" coordinates.
There's also the additional wrinkle that
map_coordinates always preserves the dtype of the input array in the output. If you interpolate an integer array, you'll get integer output, which may or may not be what you want. For the code snippet below, I'll assume that you always want floating point output. (If you don't, it's actually simpler.)
I'll try to add more explanation later tonight (this is rather dense code).
All in all, the
interp3 function I have is more complex than it may need to be for your exact purposes. Howver, it should more or less replicate the behavior of
interp3 as I remember it (ignoring the "zooming" functionality of
interp3(data, zoom_factor), which
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) y = np.linspace(10, 20, data.shape) z = np.linspace(-100, 0, data.shape) # 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()