[FIXED] interpolate 3D volume with numpy and or scipy

Issue

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

Solution

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.

Answered By - Steve Byrnes