[Octave-bug-tracker] [bug #31352] interp2 broken for non-quadratic matrices when using cubic interpolation

[Thorsten Meyer]

Follow-up Comment #13, bug #31352 (project octave):

Dear Marco,

finally I have found the time to have a look at your latest version of
bicubic.m. Thanks for your good work.

To answer your question: yes, I could still construct a test case, where
bicubic fails but the old interp2 manages to deliver (on my machine):

  x = linspace (-3, 3, 8000);
  y = linspace (-5, 7, 9111)';
  A = (1/10000 * y.^4 + 3 * sin (y)) * (x.^5 - 10 * x.^3) - y * cos (x) *
10;

  xi = linspace (min(x), max(x), 19071);
  yi = linspace (min(y), max(y), 19071);

  Z1 = bicubic (x, y, A, xi, yi);

gives "memory exhausted" for me.

while
  Z2 = interp2_old (x, y, A, [xi; xi], [yi; yi], "cubic");
succeeds. Even this is possible:
  xi = linspace (min(x), max(x), 1907100);
  yi = linspace (min(y), max(y), 1907100);
  Z2 = interp2_old (x, y, A, [xi; xi], [yi; yi], "cubic");

I think the reason for this is that your bicubic function uses the same
amount of memory for the intermediate matrix
(ZI = ZI * Cxy) as in the full meshgrid case.

Then, I couldn't resist and wrote yet another version of bicubic: see the
attached file bicubic_tm.m. It is based on the bicubic code in the old
interp2.m. Consequently, it can cope with the above example:
  xi = linspace (min(x), max(x), 1907100);
  yi = linspace (min(y), max(y), 1907100);
  Z1 = bicubic_tm (x, y, A, xi, yi);

However, in case of full meshgrid interpolation, it is much slower than your
code:
  x = linspace (-3, 3, 8000);
  y = linspace (-5, 7, 1111)';
  A = (1/10000 * y.^4 + 3 * sin (y)) * (x.^5 - 10 * x.^3) - y * cos (x) *
10;

  xi = linspace (min(x), max(x), 1071);
  yi = linspace (min(y), max(y), 821)';

  tic; Z5 = bicubic (x, y, A, xi, yi); toc
Elapsed time is 0.1334 seconds.

  tic; Z6 = bicubic_tm (x, y, A, xi, yi); toc
Elapsed time is 3.118 seconds.

For interpolating vectors of points, bicubic_tm can be a lot faster:
  x = linspace (-3, 3, 8000);
  y = linspace (-5, 7, 9111)';
  A = (1/10000 * y.^4 + 3 * sin (y)) * (x.^5 - 10 * x.^3) - y * cos (x) *
10;
  xi = min(x) + rand(10000,1) * (max(x) - min(x));
  yi = min(y) + rand(10000,1) * (max(y) - min(y));
  tic; Z7 = bicubic (x, y, A, xi, yi); toc
Elapsed time is 3.193 seconds.
  tic; Z8 = bicubic_tm (x, y, A, xi, yi); toc
Elapsed time is 0.0745 seconds.

Maybe we should include two implementations in bicubic: one optimized for the
full meshgrid case the other for interpolating vectors of points. They could
be either chosen by an input parameter or selected automatically depending on
the xi, yi input (depending on how well the results of the two implementations
can be matched numerically).

A few additional details in the code:
1. bicubic_tm seems to be a bit more robust wrt rounding errors in the step
size:
  x = linspace (-3, 3, 8);
  y = linspace (-5, 7, 11)';
  A = (1/10000 * y.^4 + 3 * sin (y)) * (x.^5 - 10 * x.^3) - y * cos (x) *
10;
  xi = linspace (min(x), max(x), 101);
  yi = linspace (min(y), max(y), 81)';

  Z1 = bicubic (x, y, A, xi, yi);
  Z2 = bicubic (x+100000*pi, y, A, xi+100000*pi, yi);
  Z3 = bicubic_tm (x, y, A, xi, yi);
  Z4 = bicubic_tm (x+100000*pi, y, A, xi+100000*pi, yi);
  max(max(abs(Z1-Z2)))
ans = 8.2191e-09
  max(max(abs(Z1-Z3)))
ans = 4.7962e-13
  max(max(abs(Z3-Z4)))
ans =  2.6484e-09
  I suspect this is due to the different calculation of the step size (as
x(2)-x(1) in bicubic.m or effectively as (x(end)-x(1)/N in bicubic_tm.m).
 
2. on the other hand, bicubic_tm doesn't pass the accuracy tests in your
interp2.m. E.g:
  Z = (1:4)' * (1:5);
  x = 1:5;
  y = 1:4;
  Xi = [1.2 4.7];
  Yi = [1.8 2.3 3.1];
  Zi = [2.16 8.46
        2.76 10.81
        3.72 14.57];
  assert(bicubic_tm(x, y', Z, Xi, Yi'), Zi, -5*eps);
maximum relative error 1.64477e-15 exceeds tolerance 1.11022e-15
I am not sure why this is so (or if it is really a problem).

3. in bicubic you don't test if xi, yi are really meshgridded in case they
are equally sized matrices. I think this is good. So there will be the interp2
function with full input checking and the bicubic function that doesn't waste
time with extended input checks.

4. in your bicubic.m you reuse several variable names for different purposes
(ZI, Cxy, A, D). That made the code look rather confusing for me. Are you sure
that you can save memory usage by this?

regards

Thorsten

(file #22084)
    _______________________________________________________

Additional Item Attachment:

File name: bicubic_tm.m                   Size:6 KB


    _______________________________________________________

Reply to this item at:

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

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