Thread Subject:

3d interpolation

HI all-

I have a three-D array 120*100*100, I wish to interpolate the data along the third dimension so that the end result is a matrix 120*100*199. The third dimension is equally spaced. The ouput array would have interpolated values at points between each of the rows and columns and the original data should not be touched (a simple average of nearest neighbors would be fine).

interp3 works across all elements as aa=interp3(X,3), and doesnt work.

It would be better than what I am thinking (although i am brain dead at this moment)
something like...
for i= 1: length(U-1)

Uaaa=(U(:,:,i)+U(:,:,i+1))/2; % not actually interp.. just element average

U_new=cat(3, U(:,:,i), Uaaa, U(:,:,i+1));...
end

any help would be appreciated

thanks
chris

"Christian Reiss" <christian.reiss@noaa.gov> wrote in message <hvgnbc$ohq$1@fred.mathworks.com>...
> HI all-
>
> I have a three-D array 120*100*100, I wish to interpolate the data along the third dimension so that the end result is a matrix 120*100*199. The third dimension is equally spaced. The ouput array would have interpolated values at points between each of the rows and columns and the original data should not be touched (a simple average of nearest neighbors would be fine).
>
> interp3 works across all elements as aa=interp3(X,3), and doesnt work.
>
> It would be better than what I am thinking (although i am brain dead at this moment)
> something like...
> for i= 1: length(U-1)
>
> Uaaa=(U(:,:,i)+U(:,:,i+1))/2; % not actually interp.. just element average
>
> U_new=cat(3, U(:,:,i), Uaaa, U(:,:,i+1));...
> end
>
> any help would be appreciated
>
> thanks
> chris

%A hint; if you're still looking.
A = cat(3,ones(3),2*ones(3),3*ones(3));
[xold yold zold] = ndgrid(1:3,1:3,1:2:5);
[xnew ynew znew] = ndgrid(1:3,1:3,1:5);

griddata3(xold,yold,zold,A,xnew,ynew,znew)

"Sean " <sean.dewolski@nospamplease.umit.maine.edu> wrote in message <i02ar8$7on$1@fred.mathworks.com>...

>
> %A hint; if you're still looking.
> A = cat(3,ones(3),2*ones(3),3*ones(3));
> [xold yold zold] = ndgrid(1:3,1:3,1:2:5);
> [xnew ynew znew] = ndgrid(1:3,1:3,1:5);
>
> griddata3(xold,yold,zold,A,xnew,ynew,znew)
==============

This will be quite slow in higher dimensions.

Because MATLAB's interpolation functions are rather slow, it is probably best to do this upsampling by sparse matrix multiplication (or get a faster interpolation routine from the FEX).

Below is a speed comparison of the method Bruno described using interp1 (in a duplicate thread) to essentially the same method replacing interp1 by a sparse matrix multiplication.


U=rand( 120,100,100);

%Method interp1
tic;

    p = size(U,3);
    Unew = permute(U,[3 1 2]);
    Unew = interp1(1:p,Unew,1:0.5:p);
    Unew = ipermute(Unew,[3 1 2]);

toc;
%Elapsed time is 0.340892 seconds.

 

%Method: sparse matrix
tic;

    [nx,ny,N]=size(U);
    A=[speye(N),spdiags(ones(N,2),-1:0,N,N-1)/2];
    A(:,[1:2:end,2:2:end])=A;
    A=A.';
    
    p = size(U,3);
    Unew = reshape(U,[],p).';
    Unew = A*Unew;
    Unew2=reshape(Unew.',nx,ny,[]);

toc;
%Elapsed time is 0.067520 seconds.





Also, as i;llustated below you can abbreviate the syntax of the sparse matrix approach using my KronProd tool

http://www.mathworks.com/matlabcentral/fileexchange/25969-efficient-object-oriented-kronecker-product-manipulation

You can also save the KronProd object for later reuse on arrays of the same size.



Z=KronProd({1,1,A},[1,2,3],[nx,ny,N]); %make KronProd object

tic; %Assume KronProd object has been saved previously for reuse

    Unew=Z*U;

toc;
%Elapsed time is 0.042638 seconds.