bugGNU Octave - Bugs: bug #29601, bugs in cubic and spline...

 
 

bug #29601: bugs in cubic and spline interpolation with interp2

Submitter:  Thorsten Meyer <tmeyier>
Submitted:  Sun 18 Apr 2010 09:03:51 PM UTC
   
 
Category:  None Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  tmeyier
Originator Name:  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

Tue 03 Aug 2010 07:33:13 PM UTC, comment #6: 

Patches applied.  All demos run now.  Closing bug.

Rik <rik5>
Group administrator
Mon 19 Jul 2010 07:26:32 AM UTC, comment #5: 

Sorry about the late reply.

If memory serves me correct that code comes from the 'image' package where interpolation is used when e.g. resizing images. When working with images you can assume that XI and YI correspond to pixel coordinates and as such will be constant in one direction and have a constant increment in the other direction. This was essentially what the code in question was ensuring. This check is, however, a bit too strong as the interpolation code does not require (I think) the increment to be constant; the input just have to be "meshgridded".

Søren Hauberg <hauberg>
Fri 09 Jul 2010 06:24:37 PM UTC, comment #4: 

Thanks for the code review. I have pushed the patch (changeset 6e7590d003dc).

I am not closing the bug yet, because I would like to understand what the part of the code is good for that caused the error in demo example 7 (i.e. interp2.m line 347 ff., see the original bug report below). Mercurial says, it originates from changeset 4f4873f6f873:

o  changeset:   9755:4f4873f6f873
|  user:        Soren Hauberg <hauberg@gmail.com>
|  date:        Tue Oct 20 20:22:10 2009 +0200
|  summary:     general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format

Does that perhaps ring a bell in your memory, Soren? Any hint is appreciated.

Thorsten Meyer <tmeyier>
Sun 04 Jul 2010 10:30:37 AM UTC, comment #3: 

Sorry about the late reply (the last half year has been crazy for me...). To be honest I am not sure if I am the author of the code in question (I might be; I just don't remember), but I'll do a review anyway :-)

  - ## If X and Y vectors produce a grid from them
  + ## Check dimensions of X and Y
  if (isvector (X) && isvector (Y))
  X = X(:).';
  Y = Y(:);

Seems like a reasonable change to me (the new comment is more descriptive).

  - ## If Xi and Yi are vectors of different orientation build a grid
  - if ((rows (XI) == 1 && columns (YI) == 1)
  - || (columns (XI) == 1 && rows (YI) == 1))
  - ## Do nothing
  + ## Check dimensions of XI and YI
  + if (isvector (XI) && isvector (YI))
  + XI = XI(:).';
  + YI = YI(:);

This seems simpler then the current code. I am not sure if the

  + XI = XI(:).';
  + YI = YI(:);

part is really needed, but I guess it is better to play it safe, so this seems fine to me.

  function b = isgriddata (X)
  d1 = diff (X, 1, 1);
  - d2 = diff (X, 1, 2);
  - b = all (d1 (:) == 0) & all (d2 (:) == d2 (1));
  + b = all (d1 (:) == 0);
  endfunction

This removes the check that the derivative of X is constant. I guess this is not a requirement for X to be gridded (consider e.g. the output of 'meshgrid (logspace (1, 10, 10))'). So, I think this change is fine as well.

In summary, I'd say check in the change.

Søren Hauberg <hauberg>
Wed 28 Apr 2010 07:15:06 AM UTC, comment #2: 

The failure of demo examples 7-10 is confirmed.  It would be best if Søren could evaluate the patch attached to this bug report and decide whether it is the correct strategy for resolving this bug.

Rik <rik5>
Group administrator
Mon 19 Apr 2010 08:52:02 PM UTC, comment #1: 

The attached patch fixes the fails of demo("interp2"). But now the cubic interpolation code within interp2 will never be used (as far as I can see). Could you have a look Søren (as you seem to be author of that code).

--- a/scripts/general/interp2.m Sun Apr 18 21:27:53 2010 +0200
+++ b/scripts/general/interp2.m Mon Apr 19 22:41:16 2010 +0200
@@ -318,7 +318,7 @@
 
   else
 
-    ## If X and Y vectors produce a grid from them
+    ## Check dimensions of X and Y
     if (isvector (X) && isvector (Y))
       X = X(:).';
       Y = Y(:);
@@ -332,10 +332,10 @@
       endif
     endif
 
-    ## If Xi and Yi are vectors of different orientation build a grid
-    if ((rows (XI) == 1 && columns (YI) == 1)
-       || (columns (XI) == 1 && rows (YI) == 1))
-      ## Do nothing
+    ## Check dimensions of XI and YI
+    if (isvector (XI) && isvector (YI))
+      XI = XI(:).';
+      YI = YI(:);
     elseif (! size_equal (XI, YI))
       error ("XI and YI must be matrices of same size");
     endif
@@ -409,8 +409,7 @@
 
 function b = isgriddata (X)
   d1 = diff (X, 1, 1);
-  d2 = diff (X, 1, 2);
-  b = all (d1 (:) == 0) & all (d2 (:) == d2 (1));
+  b = all (d1 (:) == 0);
 endfunction
 
 ## Compute the bicubic interpolation coefficients

Thorsten Meyer <tmeyier>
Sun 18 Apr 2010 09:03:51 PM UTC, original submission:  

There seem to be several bugs in the cubic and spline interpolation code of interp2.

Reproduce with:
octave> demo interp2
[...]
interp2 example 7:
 A=[13,-1,12;5,4,3;1,6,2];
 x=[0,1,2]; y=[10,11,12];
 xi=linspace(min(x),max(x),17);
 yi=linspace(min(y),max(y),26)';
 mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic'));
 [x,y] = meshgrid(x,y);
 hold on; plot3(x(:),y(:),A(:),"b*"); hold off;

interp2 example 7: failed
mx_el_or: nonconformant arguments (op1 is 1x17, op2 is 26x1)
error: called from:
error:   /home/thorsten/hg/octave/scripts/general/interp2.m at line 352, column 8

The offending line is
        inside = !(XI < X (1) | XI > X (end) | YI < Y (1) | YI > Y (end));

To me it seems that this works only if xi and yi are meshgrid data as well. But this case is handled in the preceding conditional. Strange... Also strange: the bicubic function can handle the input of the above example, i.e.
  mesh(xi,yi,bicubic(x,y,A,xi,yi));
works. So what is the purpose of the 2nd implementation of cubic interpolation in interp2.m?

next failure of the demos:
[...]
interp2 example 8:
 [x,y,A] = peaks(10);
 x = x(1,:)'; y = y(:,1);
 xi=linspace(min(x),max(x),41);
 yi=linspace(min(y),max(y),41)';
 mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic'));
 [x,y] = meshgrid(x,y);
 hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
interp2 example 8: failed
interp2: input data must have `meshgrid' format

Here, the data generated by peaks are not accepted as meshgrid data due to a floating point problem. Also, the function isgriddata (defined in interp2.m) seems to require not only equal lines in the grid matrix but also equally spaced values. Is that intended (or really necessary?).

and the last failure in the demos:
interp2 example 9:
 A=[13,-1,12;5,4,3;1,6,2];
 x=[0,1,2]; y=[10,11,12];
 xi=linspace(min(x),max(x),17);
 yi=linspace(min(y),max(y),26)';
 mesh(xi,yi,interp2(x,y,A,xi,yi,'spline'));
 [x,y] = meshgrid(x,y);
 hold on; plot3(x(:),y(:),A(:),"b*"); hold off;

interp2 example 9: failed
interp2: input data must have `meshgrid' format

Looking at the sources, the _splinen_ function that actually calculates the interpolation is called only with the first row/column of the x and y variables, so why require meshgrids here?

regards

Thorsten

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

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by hauberg (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by tmeyier
  • -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 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2010-08-03 rik5 StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2010-05-08 rik5 StatusConfirmed In Progress
    2010-04-28 rik5 StatusNone Confirmed
    2010-04-19 tmeyier Assigned toNone tmeyier
        Carbon-Copy- Added søren hauberg <soren@hauberg.org>

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code