interpolate over a multi-dimensional function

11 views (last 30 days)
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.
So instead of doing:
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?

Answers (3)

the cyclist
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.)
  1 Comment
Lior Shapira
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.

Sign in to comment.


Alec
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.
  1 Comment
the cyclist
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.

Sign in to comment.


the cyclist
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);
  3 Comments
the cyclist
the cyclist on 7 Apr 2011
Alec, interesting stuff. I have to admit that I have not thought it thoroughly through, but my gut does not agree with you that the planes need to be "far apart" to "get away with it". By construction, there is no variation in what I defined as the z-axis, so I think you will theoretically get strictly the same solution regardless of the (arbitrary) distance between the z-planes. Although there may be some computer-y round-off issues, I guess.
Regarding performance, I did a quick bit of testing, and on the largest arrays that did not bring my machine to a whimpering halt, the vectorized method was about twice as fast. This is definitely a "your mileage may vary" situation. Lior, I suggest you do some test runs with your method and mine, and see what is faster.
Lior Shapira
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.

Sign in to comment.

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!