## Copyright (C) 2010 Marco Caliari
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {@var{zi} =} bicubic (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi})
## @deftypefnx {Function File} {@var{zi} =} bicubic (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{extrapval})
## @deftypefnx {Function File} {@var{zi} =} bicubic (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{extrapval}, @var{a})
## 
## Return a matrix @var{zi} corresponding to the bicubic
## interpolation at @var{xi} and @var{yi} of the data supplied
## as @var{x}, @var{y} (equispaced meshgrid format or vectors of the same 
## length) and @var{z}.  
##
## Points outside the grid are set to @var{extrapval}. If it is not
## specified, NaN is used.
##
## The parameter @var{a} can be used to get different bicubic
## interpolations. If not specified, @var{a} is equal to -1/2.
## @end deftypef

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

function ZI = bicubic (X, Y, Z, XI, YI, extrapval, a)
  if (min (size (Z)) < 3)
    error ('bicubic: Z must be at least a 3 by 3 matrix');
  end
  if (nargin < 7)
    a = -0.5;
  endif
  if (nargin < 6)
    extrapval = NaN;
  endif
  [Ny, Nx] = size (Z);
  if (isvector (X) && isvector (Y))
    X = X(:);
    Y = Y(:);
  elseif (size_equal (X,Y))
    X = X(1,:).';
    Y = Y(:,1);
  endif
  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   
  ZI = [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
  ZI(1, 1) = 3 * ZI(2, 1) - 3 * ZI(3, 1) + ZI(4, 1);
  ZI(1, Nx + 2) = 3 * ZI(2, Nx + 2) - 3 * ZI(3, Nx + 2)+ ZI(4, Nx + 2);
  ZI(Ny + 2, 1) = 3 * ZI(Ny + 1, 1)- 3 * ZI(Ny, 1) + ZI(Ny - 1, 1);
  ZI(Ny + 2, Nx + 2) = 3  * ZI(Ny + 1, Nx + 2) - ...
      3 * ZI(Ny, Nx + 2) + ZI(Ny - 1, Nx + 2);
  XIc = XI(:).';
  YIc = YI(:).';
  if (length (YIc) * length (X) < length (XIc) * length (Y))
    indexes = max( min (floor ((rows (Y) - 1) * (YIc - Y(1)) / ...
			  (Y(end) - Y(1))), rows (Y) - 3), 1);
    A = [indexes; indexes + 1; indexes + 2; indexes + 3];
    D = cubkernel (Y(A) - YIc(ones(1, 4), :), hy, a);
    init = sub2ind([length(Y),length(YIc)], indexes, 1:length(YIc));
    Cxy = spalloc (length (Y),length (YIc),4 * length (Y));
    Cxy([init; init + 1; init + 2; init + 3](:)) = D(:);
    ZI = Cxy' * ZI;
    indexes = max (min (floor ((rows (X) - 1) * (XIc - X(1)) / ...
			  (X(end) - X(1))), rows (X) - 3), 1);
    A = [indexes; indexes + 1; indexes + 2; indexes + 3];
    D = cubkernel (X(A) - XIc(ones(1, 4), :), hx, a);
    init = sub2ind([length(X),length(XIc)], indexes, 1:length(XIc));
    Cxy = spalloc (length (X),length (XIc),4 * length (X));
    Cxy([init; init + 1; init + 2; init + 3](:)) = D(:);
    if ((size (XI, 1) == 1 && size (YI, 2) == 1) || ...
	(size (XI, 2) == 1 && size (YI, 1) == 1))
## ZI corresponds to meshgrid (XI,YI)    
      ZI = ZI * Cxy;
## set points outside X,Y to extrapval
      ZI (:, XI < X(2) | XI > X(Nx + 1)) = extrapval;
      ZI (YI < Y(2) | YI > Y(Ny + 1), :) = extrapval;
    else
      ZI = full (reshape (sum (ZI .* Cxy', 2), size (XI)));
      ZI (XI < X(2) | XI > X(Nx + 1) | YI < Y(2) | YI > Y(Ny + 1)) = extrapval;
    endif
  else
    indexes = max (min (floor ((rows (X) - 1) * (XIc - X(1)) / ...
			  (X(end) - X(1))), rows (X) - 3), 1);
    A = [indexes; indexes + 1; indexes + 2; indexes + 3];
    D = cubkernel (X(A) - XIc(ones(1, 4), :), hx, a);
    init = sub2ind([length(X),length(XIc)], indexes, 1:length(XIc));
    Cxy = spalloc (length (X),length (XIc),4 * length (X));
    Cxy([init; init + 1; init + 2; init + 3](:)) = D(:);
    ZI = ZI * Cxy;
    indexes = max (min (floor ((rows (Y) - 1) * (YIc - Y(1)) / ...
			  (Y(end) - Y(1))), rows (Y) - 3), 1);
    A = [indexes; indexes + 1; indexes + 2; indexes + 3];
    D = cubkernel (Y(A) - YIc(ones(1, 4), :), hy, a);
    init = sub2ind([length(Y),length(YIc)], indexes, 1:length(YIc));
    Cxy = spalloc (length (Y),length (YIc),4 * length (Y));
    Cxy([init; init + 1; init + 2; init + 3](:)) = D(:);
    if ((size (XI, 1) == 1 && size (YI, 2) == 1) || ...
	(size (XI, 2) == 1 && size (YI, 1) == 1))
## ZI corresponds to meshgrid (XI,YI)    
      ZI = Cxy' * ZI;
## set points outside X,Y to extrapval
      ZI (:, XI < X(2) | XI > X(Nx + 1)) = extrapval;
      ZI (YI < Y(2) | YI > Y(Ny + 1), :) = extrapval;
    else
      ZI = full (reshape (sum (ZI .* Cxy, 1), size (XI)));
      ZI (XI < X(2) | XI > X(Nx + 1) | YI < Y(2) | YI > Y(Ny + 1)) = extrapval;
    endif
  endif  
endfunction

function u = cubkernel (s, h, a)
  s = abs (s);
  idx1 = s < h;
  s1 = s (idx1) / h;
  idx2 = h < s & s < 2 * h;
  s2 = s (idx2) / h;
  m = size (s,1);
  n = size (s,2);
  u = spalloc (m, n, 4 * m);
  u (idx1) = ((a + 2) * s1 - (a + 3)) .* s1 .^ 2 + 1;
  u (idx2) = a * (((s2 - 5) .* s2 + 8) .* s2 - 4);
endfunction
