bugGNU Octave - Bugs: bug #30587, bicubic can't handle non-gridded...

 
 

bug #30587: bicubic can't handle non-gridded targets / bug in the bicubic interpolation algorithm

Submitter:  Thorsten Meyer <tmeyier>
Submitted:  Thu 29 Jul 2010 08:53:54 AM UTC
   
 
Category:  Libraries Severity:  2 - Minor
Priority:  5 - Normal Item Group:  Feature Request
Status:  Duplicate Assigned to:  None
Originator Name:  Thorsten Meyer Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 29 Oct 2010 07:56:57 AM UTC, comment #7: 

This is being treated in <a href=https://savannah.gnu.org/bugs/?31352>#31352</a> so I'm closing this bug as a duplicate (even though this bug report was first)

D.



David Bateman <dbateman>
Group Member
Mon 30 Aug 2010 07:01:33 AM UTC, comment #6: 

Yes, both implementation (bicubic and interp2) are based on the paper I already quoted, but the choice of a free parameter and boundary conditions are different. Bicubic is more similar to matlab, but, as Thorsten noted, it works only for equispaced x and y data. As I said, I think I can fix their behavior, but I'm quite busy for a couple of weeks.

Marco

Marco Caliari <caliari>
Group Member
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'.

Thorsten Meyer <tmeyier>
Mon 23 Aug 2010 09:10:46 AM UTC, comment #4: 

I've disabled "bicubic" call from interp2 entirely, and I'm softening this bug to a feature request. If anyone feels up to fixing bicubic to handle non-gridded targets, please do.

Jaroslav Hajek <highegg>
Mon 16 Aug 2010 08:19:59 AM UTC, comment #3: 

I think the behavior where column x row vectors are automatically meshgridded is taken from Matlab. Note that an implementation doesn't actually need to form the meshgrid at all. So, as far as I understand, the bug is that bicubic can't handle non-meshgrid inputs. Given that bicubic is available as stand-alone function, unless someone extends bicubic I suggest that interp2 won't call bicubic at all, to keep the results consistent. Picking faster/more accurate code for special cases is OK, but two fundamentally different methods shouldn't be used under the same name.

Jaroslav Hajek <highegg>
Fri 13 Aug 2010 04:06:52 PM UTC, comment #2: 

Dear Marco,

thanks for your feedback.

1. You mention the case where xi and yi are vectors of the same length and orientation. At the moment, interp2 passes these to bicubic.m and there they will be used to define a meshgrid on which to interpolate. You say that the two vectors should be interpreted as scattered vectors instead (giving only a vector of interpolated values). If I understand you correctly, interp2 would:
 - interpret xi, yi as vectors of scattered data if xi and yi are of same length and orientation
 - use xi, yi to define a meshgrid when they have different orientation and/or different length

I personally would prefer an explicit option to force interp2 to treat xi, yi vectors this or the other way. What does matlab do in these cases?

2. You ask, why I should give xi, yi as matrices, if they represent scattered data. Well, actually, I was trying to understand the existing code and found only this way to get the bicubic interpolation algorithm within interp2 to be used. I also don't think, that such a use would be particularly useful (once issue 1 above has been resolved in some way).

Regards

Thorsten

Thorsten Meyer <tmeyier>
Wed 11 Aug 2010 02:50:17 PM UTC, comment #1: 

Dear Thorsten,

the bicubic.m interpolation is "better". It corresponds to

R. Keys, (1981). "Cubic convolution interpolation for digital image processing". IEEE Transactions on Signal Processing, Acoustics, Speech, and Signal Processing 29

In principle, I can fix the implementation inside interp2.m in order to match the implementation in bicubic.m (I can do for sure if x and y are both uniformly spaced vectors). First, there is something else to fix: if xi and yi are vectors of the same length and orientation (scattered data), then the algorithm inside interp2.m should be used. I don't understand your input: why should you give the points xi and yi as matrices not coming from a meshgrid?

Marco

Marco Caliari <caliari>
Group Member
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?

Thorsten Meyer <tmeyier>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Digest:
   bug dependencies.

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by dbateman (Posted a comment)
  • -email is unavailable- added by highegg (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by tmeyier (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2010-10-29 dbateman StatusPostponed Duplicate
        Open/ClosedOpen Closed
        Dependencies- Depends on bugs #31352
    2010-08-28 tmeyier Summarybicubic can\'t handle non-gridded targets bicubic can't handle non-gridded targets / bug in the bicubic interpolation algorithm
    2010-08-23 highegg Severity4 - Important 2 - Minor
        Item GroupIncorrect Result Feature Request
        StatusNone Postponed
        Summarybuggy implementation of bicubic interpolation in interp2/bicubic bicubic can't handle non-gridded targets

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code