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

Submitted by:  Thorsten Meyer <tmeyier>
Submitted on:  Thu 29 Jul 2010 08:53:54 AM UTC  
 
Category: LibrariesSeverity: 2 - Minor
Priority: 5 - NormalItem Group: Feature Request
Status: DuplicateAssigned to: None
Originator Name: Thorsten MeyerOpen/Closed: Closed
Release: devOperating System: GNU/Linux

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

(Jump to the original submission 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>
Project 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>
Project 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>
Project 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 File(s):
   
   
Comment:
   

No files currently attached

 

Digest:
   bug dependencies.

Items that depend on this one: None found

 

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

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 8 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Fri 29 Oct 2010 07:56:57 AM UTCdbatemanStatusPostponed=>Duplicate
      Open/ClosedOpen=>Closed
      Dependencies-=>Depends on bugs #31352
    Sat 28 Aug 2010 02:58:45 PM UTCtmeyierSummarybicubic can\'t handle non-gridded targets=>bicubic can't handle non-gridded targets / bug in the bicubic interpolation algorithm
    Mon 23 Aug 2010 09:10:46 AM UTChigheggSeverity4 - 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.1-cleanup1