interpolate 3D volume with numpy and or scipy

Question

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?

Thanks

Answer 1

In scipy 0.14 or later, there is a new function scipy.interpolate.RegularGridInterpolator which closely resembles interp3.

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

MATLAB CODE:

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 (2,6,8) and (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))

Again it's [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.