# interpolate over a multi-dimensional function

26 views (last 30 days)
Lior Shapira on 6 Apr 2011
I want to use interpn to interpolate the values of function V:R^4->R^3. I.e. size(V) = [C C C C 3], but interpn requires that the last dimension be of size 1.
anew = interpn(C, M, Y, K, V, a(:,1), a(:,2), a(:,3), a(:,4));
I have to do this:
anew = zeros(size(a_cmyk));
for i=1:4
anew(:,i) = interpn(C, M, Y, K, V(:,:,:,:,i),a(:,1), a(:,2), a(:,3), a(:,4));
end
And it really slows me down. Am I doing something wrong? Is there a way around this?

the cyclist on 6 Apr 2011
I think you may be mistaken, about the last dimension of V needing to be 1. In the example in the help file, their V is length 6 in the final dimension, and it works fine. (They do use a loop over the final dimension in that example, but only for the visualization.)
Lior Shapira on 7 Apr 2011
Here's a simple test:
c = 1:100;
m = 1:100;
y = 1:100;
V = rand([100 100 100 4]);
[C M Y] = ndgrid(c,m,y);
vals = interpn(C,M,Y,V, 1:50,1:50,1:50);
It doesn't work, returns an error:
??? Error using ==> interpn at 155
Wrong number of input arguments or some dimension of V is less than 2.

Alec on 7 Apr 2011
It seems interpn only works with scalar functions. If you issue
open interpn
you can see where it parses the input it expects that you have as many X1,X2,... as dimensions in V. So V has to be a scalar function: R^n --> R. You'll have to interpolate each coordinate function independently.
In your original post, you write that V is a function from R^4 --> R^3, but then in the code V seems to be a function from R^4 --> R^4 (you're looping from 1 to 4 rather than 1 to 3). Was that just a typo?
In the example you posted in the comment. V is a function from R^3 to R^4. So you'll have to interpolate each of the 4 coordinate functions separately, perhaps like this:
c = 1:100;
m = 1:100;
y = 1:100;
V = rand([100 100 100 4]);
[C M Y] = ndgrid(c,m,y);
vals = zeros([size(1:50) 4]);
for ii = 1:4
vals(:,:,ii) = interpn(C,M,Y,V(:,:,:,ii), 1:50,1:50,1:50);
end
(Note: vals here has size(vals) = [ 1 50 4] because the 50 interpolated points are defined on a "grid" that only has one row )
As far as performance, the for loop isn't really a "bad" for loop in terms of matlab optimization. Presuming the number of coordinates is much much smaller than the number of data points or interpolated points.
the cyclist on 7 Apr 2011
If I am not mistaken, this is pretty much the same coding strategy that Lior included with her original question.
See my answer, posted (time-wise) just after yours, in which I think I have the correct vectorization of the last dimension. Performance wise, I am not sure if it will be faster or not, but it is worth a shot.

the cyclist on 7 Apr 2011
Lior,
Putting this as a new answer, because I did not have a perfect understanding of what you were trying to do before.
I am going to put this code in without explanation, because I think you are going to understand what it is doing. If you do not, add a comment, and I will try to explain a bit what is happening. The code that does what you want is the one commented "Four planes".
c = 1:20;
m = 1:20;
y = 1:20;
z = 1:4;
V = rand([20 20 20 4]);
[C M Y Z] = ndgrid(c,m,y,z);
% One plane:
ci = 1:10;
mi = 1:10;
yi = 1:10;
zi = ones(1,10);
vals = interpn(C,M,Y,Z,V,ci,mi,yi,zi);
% Four planes:
ci = repmat(1:10,[1 4]);
mi = repmat(1:10,[1 4]);
yi = repmat(1:10,[1 4]);
zi = [ones(1,10), 2*ones(1,10), 3*ones(1,10), 4*ones(1,10)];
vals = interpn(C,M,Y,Z,V,ci,mi,yi,zi);
Lior Shapira on 11 Apr 2011
Thanks for the discussion guys. I'm running tests to see if this can help me. Also I'm writing the interpolation part as a MEX file, which will hopefully speed things up.

### Categories

Find more on Creating and Concatenating Matrices in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by