Path: news.mathworks.com!not-for-mail
From: "us " <us@neurol.unizh.ch>
Newsgroups: comp.soft-sys.matlab
Subject: Re: vectorizing interp2 for use in 3D arrays?
Date: Wed, 22 Aug 2007 18:39:07 +0000 (UTC)
Organization: Universit&#228;tsSpital Z&#252;rich
Lines: 85
Message-ID: <fahvsb$ku5$1@fred.mathworks.com>
References: <f93kfr$ppb$1@fred.mathworks.com> <f9sdak$oci$1@fred.mathworks.com> <f9sqfn$rpb$1@fred.mathworks.com> <fahepe$bil$1@fred.mathworks.com>
Reply-To: "us " <us@neurol.unizh.ch>
NNTP-Posting-Host: webapp-02-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1187807947 21445 172.30.248.37 (22 Aug 2007 18:39:07 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Wed, 22 Aug 2007 18:39:07 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 11
Xref: news.mathworks.com comp.soft-sys.matlab:425027


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