Thread Subject: vectorizing interp2 for use in 3D arrays?

Hi,

I am trying to see if it is possible to avoid a loop in the
following problem regarding the use of interp2.

Suppose Z is an N-by-M-by-K 3D-array
I want to interpolate each of the K M-by-N 2D-arrays over a
finer grid.
The only solution I could think of is to use interp2 on a
loop over the K dimension as follows:

for k=1:K
     ZI(:,:,k) = interp2(X, Y, Z(:,:,k), XI, YI, '*linear');
end

Is there a better and faster alternative?

Thanks.

Lorenzo Garlappi

Lorenzo Garlappi:
<SNIP looking for <interpX>...

a hint:

     help interp3;
     help interpn;

us

"Lorenzo Garlappi" <lorenzo.garlappi@gmail.com> wrote in message
<f93kfr$ppb$1@fred.mathworks.com>...
> I am trying to see if it is possible to avoid a loop in the
> following problem regarding the use of interp2.
>
> Suppose Z is an N-by-M-by-K 3D-array
> I want to interpolate each of the K M-by-N 2D-arrays over a
> finer grid.
> The only solution I could think of is to use interp2 on a
> loop over the K dimension as follows:

Use of interp3 or interpn is for some reason
out of the question?

John

Thanks John,
I was worried about execution time and was under the
impression that interp3 would be slower than interp2 with
loop. I will experiment.

Lorenzo

"John D'Errico" <woodchips@rochester.rr.com> wrote in
message <f94503$565$1@fred.mathworks.com>...
> "Lorenzo Garlappi" <lorenzo.garlappi@gmail.com> wrote in
message
> <f93kfr$ppb$1@fred.mathworks.com>...
> > I am trying to see if it is possible to avoid a loop in the
> > following problem regarding the use of interp2.
> >
> > Suppose Z is an N-by-M-by-K 3D-array
> > I want to interpolate each of the K M-by-N 2D-arrays over a
> > finer grid.
> > The only solution I could think of is to use interp2 on a
> > loop over the K dimension as follows:
>
> Use of interp3 or interpn is for some reason
> out of the question?
>
> John
>
>

interp2 and interp3 are not doing the same job and will not
give you exactly the same results! You need to figure out if
you need 2D or 3D interpolation. (i.e. Do you need to
interpolate in the 3rd dimension or not?)

It seems that you don't need to interpolate in the 3rd
dimension. So the loop-interp2 combination is probably your
best choice!

However it seems that what you actually want to do is not
only interpolation but regridding. If your regridding ratio
can be a power of 2 you should probably use ZI =
interp2(Z,ntimes) instead of creating your XI, YI with
meshgrid. For my application it is 40% slower to use interp2
with meshgrid instead of "ntimes".


Aurelien


"Lorenzo Garlappi" <lorenzo.garlappi@gmail.com> wrote in
message <f93kfr$ppb$1@fred.mathworks.com>...
> Hi,
>
> I am trying to see if it is possible to avoid a loop in the
> following problem regarding the use of interp2.
>
> Suppose Z is an N-by-M-by-K 3D-array
> I want to interpolate each of the K M-by-N 2D-arrays over a
> finer grid.
> The only solution I could think of is to use interp2 on a
> loop over the K dimension as follows:
>
> for k=1:K
> ZI(:,:,k) = interp2(X, Y, Z(:,:,k), XI, YI, '*linear');
> end
>
> Is there a better and faster alternative?
>
> Thanks.
>
> Lorenzo Garlappi

Aurelien Stalder:
<SNIP down to good point...

> interp2 and interp3 are not doing the same job and will
not give you exactly the same results...

correct - albeit, down to FP minutes...
however, the problem is timing, as one can see running the
snippet below, which yields this on a
wintel system (p5/1.3g/512m/win2k.sp4/r2007a)

Elapsed time is 0.304378 seconds. <- loop
Elapsed time is 0.649021 seconds. <- 3d
     -0.46836 10.816 11.284 <- loop
     -0.46836 10.816 11.284 <- 3d
 -3.5527e-015 3.5527e-015 7.1054e-015 <- loop-3d

just a thought
us

% the data
     nr=10;
     nc=25;
     nz=50;
     m=10*rand(nr,nc,nz);
     x=1:nc;
     y=1:nr;
     z=1:nz;
     xi=1:.5:nc;
     yi=1:.5:nr;
     zi=ones(size(z));
% the loop
tic
     m1=zeros(numel(yi),numel(xi),numel(zi));
for i=z
     m1(:,:,i)=interp2(x,y.',m(:,:,1),xi,yi.','cubic');
end
toc
% the 3d
tic
     m2=interp3(x,y.',z,m,xi,yi.',zi,'cubic');
toc
% the result
% - equal within FP accuracy...
     r=m1-m2;
     disp([
          [min(m1(:)),max(m1(:)),range(m1(:))]
          [min(m2(:)),max(m2(:)),range(m2(:))]
          [min(r(:)),max(r(:)),range(r(:))]
     ]);

"us " <us@neurol.unizh.ch> wrote in message
> correct - albeit, down to FP minutes...

What are you trying to demonstrate?!?
I don't think the point is on showing that doing z-times 2D
interpolation on the same slice give the same result as
doing 3D interpolation on a volume constant along the 3rd
dimension.

The following modifications would give a little bit of sense
to your code:
> zi=ones(size(z));
This way, 3D interpolation is just happening on the 1st
slice. What you really want is:
=> zi=1:nz;

...
> for i=z
> m1(:,:,i)=interp2(x,y.',m(:,:,1),xi,yi.','cubic');

m(:,:,1) => m(:,:,i)
Otherwise you are always doing the same interpolation on the
first slice.

3D interpolation is not equivalent to successive 2D
interpolation!!! Depending on the data, both might give
similar results but they are nevertheless intrinsically
different.

"Aurelien Stalder" <aurelien.stalder.remove@uniklinik-freiburg.de> wrote in
message <fahepe$bil$1@fred.mathworks.com>...
> 3D interpolation is not equivalent to successive 2D
> interpolation!!! Depending on the data, both might give
> similar results but they are nevertheless intrinsically
> different.

Actually, in the case of some interpolation methods, they
ARE identical.

A 3-d tensor product linear interpolant is mathematically
equivalent to a linear interpolation applied to a pair of
2-d tensor product linear interpolants.

That equivalence should also apply to nearest neighbor
interpolants.

John

Aurelien Stalder:
<SNIP slighly premature rant...

> What are you trying to demonstrate?!?
> The following modifications would give a little bit
> of sense to your code...

well, you obviously didn't take the time to carefully read
what i tried to show:
- TIMING... (i admitt, i could have been more clear about
this point!)

however, you obviously didn't bother to run your own
suggestion...
let's do it then with a modified snippet - and - let's see
what i meant by <down to FP minutes>...

% the modified snippet
% the data
     format long; % <- may have to be reset!
     amode={ % check different interpolation modes
          'nearest'
'linear'
'cubic'
     };
     stp=.25;
     nr=5;
     nc=7;
     nz=15;
     m=10*rand(nr,nc,nz);
     x=1:nc;
     y=1:nr;
     z=1:nz;
     xi=1:stp:nc;
     yi=1:stp:nr;
     zi=1:nz;
% the loop
for mix=1:numel(amode)
tic
     m1=zeros(numel(yi),numel(xi),numel(zi));
for i=z
     m1(:,:,i)=interp2(x,y.',m(:,:,i),xi,yi.',amode{mix});
end
toc
% the 3d
tic
     m2=interp3(x,y.',z,m,xi,yi.',zi,amode{mix});
toc
% the result
% - equal within FP accuracy...
     r=m1-m2;
     rr=[
          [min(m1(:)),max(m1(:)),range(m1(:))]
          [min(m2(:)),max(m2(:)),range(m2(:))]
          [min(r(:)),max(r(:)),range(r(:))]
     ];
     disp(sprintf('mode <%s> output equal %-1d',...
                   amode{mix},isequal(m1,m2)));
     disp(rr);
end

% yields these exemplary results for one run
Elapsed time is 0.085091 seconds.
Elapsed time is 0.027607 seconds.
mode <nearest> output equal 1 <- EQUAL
   0.021618029388651 9.981778860424260 9.960160831035609
   0.021618029388651 9.981778860424260 9.960160831035609
                   0 0 0

Elapsed time is 0.104865 seconds.
Elapsed time is 0.036034 seconds.
mode <linear> output equal 1 <- EQUAL
   0.021618029388651 9.981778860424260 9.960160831035609
   0.021618029388651 9.981778860424260 9.960160831035609
                   0 0 0

Elapsed time is 0.120199 seconds.
Elapsed time is 0.154096 seconds.
mode <cubic> output equal 0 <- NOT EQUAL, but SMALL errors
  -0.884119795211513 10.768503311369718 11.652623106581231
  -0.884119795211513 10.768503311369717 11.652623106581229
  -0.000000000000005 0.000000000000009 0.000000000000014

just another thought
us

"John D'Errico" <woodchips@rochester.rr.com> wrote in
message <fahk4p$64l$1@fred.mathworks.com>...
> A 3-d tensor product linear interpolant is mathematically
> equivalent to a linear interpolation applied to a pair of
> 2-d tensor product linear interpolants.
>
> That equivalence should also apply to nearest neighbor
> interpolants.

Exact! The mathematical results are identical for those low
order interpolation methods and for any interpolation method
based on a separable filter. Nevertheless, the computation
algorithms (and hence computation time) are intrinsically
different between 1D, 2D, 3D or nD interpolation methods.

Generally, if you don't need interpolation in the 3rd
dimension, you should use successive 2D interpolation (as
Lorenzo initially suggested). Depending on the interpolation
implementation, the matrix size and other aspects, you may
find exceptions to this. As "us" demonstrated, when working
on a volume of size 5x7x15, interp3 can be faster than
interp2 in a loop.

Aurelien