[Octave-bug-tracker] [bug #30587] bicubic can't handle non-gridded targets / bug in the bicubic interpolation algorithm

[Thorsten Meyer]

Update of bug #30587 (project octave):

                 Summary: bicubic can't handle non-gridded targets => bicubic
can't handle non-gridded targets / bug in the bicubic interpolation algorithm

    _______________________________________________________

Follow-up 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'. 


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?30587>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/