interpolate over a multi-dimensional function
6 次查看(过去 30 天)
显示 更早的评论
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?
0 个评论
回答(3 个)
the cyclist
2011-4-6
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.)
Alec
2011-4-7
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 个评论
the cyclist
2011-4-7
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
2011-4-7
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 个评论
the cyclist
2011-4-7
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!