Sat 28 Aug 2010 02:58:45 PM UTC, comment #5:
Actually, the original bug report was about the implementation of bicubic interpolation. The handling of non-gridded data only came up as part of the discussion. Maybe it should be a separate feature request.
About bicubic interpolation:
I think there really is a bug in the implementation of the bicubic interpolation in bicubic.m.
For example:
x=[0,1,4]+10;
y=[-10,-9,-8];
[X,Y] = meshgrid(x,y);
B=X.^2 - 10 * (Y+9).^2 + X.*Y;
xi=linspace(min(x),max(x),17);
yi=linspace(min(y),max(y),26)';
mesh(xi,yi,bicubic(x,y,B,xi,yi));
hold on; plot3(X(:),Y(:),B(:),"b*"); hold off;
xlabel("x")
ylabel("y")
You can see, that the interpolation is not continuous in the 1st derivative in x direction (as it is supposed to be).
to make it even more clear:
xx=[0,1,20];
xxi=linspace(min(xx),max(xx),30);
[XX, Y] = meshgrid(xx, y);
mesh(xxi,yi,bicubic(xx,y,B,xxi,yi));
hold on; plot3(XX(:),Y(:),B(:),"b*"); hold off;
xlabel("xx")
ylabel("y")
pause
xxi2=linspace(min(xx),max(xx),1000);
plot(xxi2, bicubic(xx,y,B,xxi2,-9));
I can also produce the same problem in y-direction:
yy=[-10 -9 -1];
yyi=linspace(min(yy),max(yy),30);
[XX, YY] = meshgrid(xx, yy);
mesh(xxi,yyi,bicubic(xx,yy,B,xxi,yyi));
hold on; plot3(XX(:),YY(:),B(:),"b*"); hold off;
xlabel("xx")
ylabel("yy")
pause
yyi2=linspace(min(yy),max(yy),1000);
xxi2=0;
plot(yyi2, bicubic(xx,yy,B,xxi2,yyi2))
I think that bicubic does not deal correctly with unequally spaced x and y data. I've also checked with matlab and, indeed, only for equally spaced x or y vectors I get the same results.
Strangely, the implementation within interp2.m gives neither the result of bicubic nor the result of matlab:
octave:
bicubic(x,y,B,xi,yi)(1:2,1:2)
ans =
-10.0000 -10.4375
-7.6640 -8.0965
matlab:
bla=interp2(x,y,B,xi,yi,'cubic')(1:2,1:2);
bla(1:2,1:2)
ans =
-10.0000 -10.4844
-7.6640 -8.1471
octave (freshly built):
interp2(x,y',B,[xi(1:3); xi(1:3)],[yi(1) yi(1) yi(1); yi(2) yi(2) yi(3)], 'cubic')(1:2,1:2)
ans =
-10.0000 -9.2402
-8.4000 -7.6340
So, if I take matlab as a reference, both implementations are buggy.
Also note, that with your latest patch (e48a45b9a265), Jaroslav, the cubic demos of interp2 do no longer work and I had to use the above trick to actually get interp2 to calculate anything with 'cubic'.
|
Thu 29 Jul 2010 08:53:54 AM UTC, original submission:
while implementing regression tests for recent bug fixes in interp2, I came across the following:
A=[13,-1,12;13,-1,12;13,-1,12];
x=[0,1,2]; y=[10,11,12];
xi=linspace(min(x),max(x),17);
yi=linspace(min(y),max(y),26)';
interp2(x,y,A,xi,yi,'cubic')(1,:)
gives this:
ans =
Columns 1 through 6:
13.000000 9.773438 6.968750 4.585938 2.625000 1.085938
Columns 7 through 12:
-0.031250 -0.726562 -1.000000 -0.851562 -0.281250 0.710938
Columns 13 through 17:
2.125000 3.960938 6.218750 8.898437 12.000000
in this case, the bicubic interpolation in bicubic.m is used and it gives what I would naively expect, i.e., a parabolic fit to the three data points:
pp=polyfit([0 1 2], [13 -1 12], 2);
polyval(pp, linspace(0, 2, 17))
ans =
Columns 1 through 6:
13.000000 9.773438 6.968750 4.585937 2.625000 1.085937
Columns 7 through 12:
-0.031250 -0.726563 -1.000000 -0.851563 -0.281250 0.710937
Columns 13 through 17:
2.125000 3.960938 6.218750 8.898438 12.000000
Now, I wanted to try the other implementation of bicubic interpolation directly within interp2, which is used for non-meshgrid xi and yi:
interp2(x,y,A,[xi; xi], [10 * ones(size(xi)); [10 11*ones(1, 16)]],'cubic')(1,:)
And here, I get the following:
ans =
Columns 1 through 7:
13.00000 11.07227 8.89062 6.60742 4.37500 2.34570 0.67188
Columns 8 through 14:
-0.49414 -1.00000 -0.71484 0.28125 1.82422 3.75000 5.89453
Columns 15 through 17:
8.09375 10.18359 12.00000
The question is: which of the two implementations is right?
|