function ZI = newbicubic (X, Y, Z, XI, YI)
##
## bicubic interpolation
## based on 
## R. Keys, (1981). "Cubic convolution interpolation for digital image 
## processing". IEEE Transactions on Signal Processing, Acoustics,
## Speech, and Signal Processing 29
##
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{zi}] =} newbicubic (@var{x},
## @var{y}, @var{z}, @var{xi}, @var{yi})
  [Ny, Nx] = size (Z);
  hx = X(2) - X(1);
  hy = Y(2) - Y(1);
  X = [X(1) - hx, X, X(Nx) + hx];
  Y = [Y(1) - hy, Y, Y(Ny) + hy];
## edge extrapolation   
  C = [0, 3 * Z(1, :) - 3 * Z(2, :) + Z(3, :), 0; ...
       3 * Z(:, 1) - 3 * Z(:, 2) + Z(:, 3), Z, ...
       3 * Z(:, Nx) - 3 * Z(:, Nx - 1) + Z(:, Nx - 2); ...
       0, 3 * Z(Ny, :) - 3 * Z(Ny - 1, :) + Z(Ny - 2, :), 0];
## corner extrapolation
  C(1, 1) = 3 * C(2, 1) - 3 * C(3, 1) + C(4, 1);
  C(1, Nx + 2) = 3 * C(2, Nx + 2) - 3 * C(3, Nx + 2)+ C(4, Nx + 2);
  C(Ny + 2, 1) = 3 * C(Ny + 1, 1)- 3 * C(Ny, 1) + C(Ny - 1, 1);
  C(Ny + 2, Nx + 2) = 3  * C(Ny + 1, Nx + 2) - ...
      3 * C(Ny, Nx + 2) + C(Ny - 1, Nx + 2);
  XIc = XI(:);
  YIc = YI(:);
  Cx = cubkernel ((X(ones (length (XIc), 1),:) - ...
		   XIc(:, ones(1, length (X)))) / hx);
  Cy = cubkernel ((Y(ones (length (YIc), 1),:) - ...
		  YIc(:, ones(1, length (Y)))) / hy);
  if (size (XI, 1) == 1 && size (YI, 2) == 1)
% ZI corresponds to meshgrid (XI,YI)    
    ZI = Cy * C * Cx.';
  else
    ZI =  (reshape (sum ((Cy * C) .* Cx, 2), size (XI)));
  endif
endfunction

function u = cubkernel (s)
  a = - 0.5;
  s = abs (s);
  idx1 = 0 <= s & s < 1;
  idx2 = 1 < s & s < 2;
  u =  zeros (size (s));
  s1 = s (idx1);
  u (idx1) = ((a + 2) * s1 - (a + 3)) .* s1 .^ 2 + 1;
  s2 = s(idx2);
  u (idx2) = a * (((s2 - 5) .* s2 + 8) .* s2 - 4);
endfunction
