## Copyright (C) 2000-2012 Kai Habel
## Copyright (C) 2009 Jaroslav Hajek
##
## 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} =} interp2 (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi})
## @deftypefnx {Function File} {@var{zi} =} interp2 (@var{Z}, @var{xi}, @var{yi})
## @deftypefnx {Function File} {@var{zi} =} interp2 (@var{Z}, @var{n})
## @deftypefnx {Function File} {@var{zi} =} interp2 (@dots{}, @var{method})
## @deftypefnx {Function File} {@var{zi} =} interp2 (@dots{}, @var{method}, @var{extrapval})
##
## Two-dimensional interpolation.  @var{x}, @var{y} and @var{z} describe a
## surface function.  If @var{x} and @var{y} are vectors their length
## must correspondent to the size of @var{z}.  @var{x} and @var{y} must be
## monotonic.  If they are matrices they must have the @code{meshgrid}
## format.
##
## @table @code
## @item interp2 (@var{x}, @var{y}, @var{Z}, @var{xi}, @var{yi}, @dots{})
## Returns a matrix corresponding to the points described by the
## matrices @var{xi}, @var{yi}.
##
## If the last argument is a string, the interpolation method can
## be specified.  The method can be 'linear', 'nearest' or 'cubic'.
## If it is omitted 'linear' interpolation is assumed.
##
## @item interp2 (@var{z}, @var{xi}, @var{yi})
## Assumes @code{@var{x} = 1:rows (@var{z})} and @code{@var{y} =
## 1:columns (@var{z})}
##
## @item interp2 (@var{z}, @var{n})
## Interleaves the matrix @var{z} n-times.  If @var{n} is omitted a value
## of @code{@var{n} = 1} is assumed.
## @end table
##
## The variable @var{method} defines the method to use for the
## interpolation.  It can take one of the following values
##
## @table @asis
## @item 'nearest'
## Return the nearest neighbor.
##
## @item 'linear'
## Linear interpolation from nearest neighbors.
##
## @item 'pchip'
## Piecewise cubic Hermite interpolating polynomial.
##
## @item 'cubic'
## Cubic interpolation from four nearest neighbors.
##
## @item 'spline'
## Cubic spline interpolation---smooth first and second derivatives
## throughout the curve.
## @end table
##
## If a scalar value @var{extrapval} is defined as the final value, then
## values outside the mesh as set to this value.  Note that in this case
## @var{method} must be defined as well.  If @var{extrapval} is not
## defined then NA is assumed.
##
## @seealso{interp1}
## @end deftypefn

## Author:      Kai Habel <kai.habel@gmx.de>
## 2005-03-02 Thomas Weber <weber@num.uni-sb.de>
##     * Add test cases
## 2005-03-02 Paul Kienzle <pkienzle@users.sf.net>
##     * Simplify
## 2005-04-23 Dmitri A. Sergatskov <dasergatskov@gmail.com>
##     * Modified demo and test for new gnuplot interface
## 2005-09-07 Hoxide <hoxide_dirac@yahoo.com.cn>
##     * Add bicubic interpolation method
##     * Fix the eat line bug when the last element of XI or YI is
##       negative or zero.
## 2005-11-26 Pierre Baldensperger <balden@libertysurf.fr>
##     * Rather big modification (XI,YI no longer need to be
##       "meshgridded") to be consistent with the help message
##       above and for compatibility.

function ZI = interp2 (varargin)
  Z = X = Y = XI = YI = n = [];
  method = "linear";
  extrapval = NA;

  switch (nargin)
    case 1
      Z = varargin{1};
      n = 1;
    case 2
      if (ischar (varargin{2}))
        [Z, method] = deal (varargin{:});
        n = 1;
      else
        [Z, n] = deal (varargin{:});
      endif
    case 3
      if (ischar (varargin{3}))
        [Z, n, method] = deal (varargin{:});
      else
        [Z, XI, YI] = deal (varargin{:});
      endif
    case 4
      if (ischar (varargin{4}))
        [Z, XI, YI, method] = deal (varargin{:});
      else
        [Z, n, method, extrapval] = deal (varargin{:});
      endif
    case 5
      if (ischar (varargin{4}))
        [Z, XI, YI, method, extrapval] = deal (varargin{:});
      else
        [X, Y, Z, XI, YI] = deal (varargin{:});
      endif
    case 6
        [X, Y, Z, XI, YI, method] = deal (varargin{:});
    case 7
        [X, Y, Z, XI, YI, method, extrapval] = deal (varargin{:});
    otherwise
      print_usage ();
  endswitch

  ## Type checking.
  if (!ismatrix (Z))
    error ("interp2: Z must be a matrix");
  endif
  if (!isempty (n) && !isscalar (n))
    error ("interp2: N must be a scalar");
  endif
  if (!ischar (method))
    error ("interp2: METHOD must be a string");
  endif
  if (ischar (extrapval) || strcmp (extrapval, "extrap"))
    extrapval = [];
  elseif (!isscalar (extrapval))
    error ("interp2: EXTRAPVAL must be a scalar");
  endif

  ## Define X, Y, XI, YI if needed
  [zr, zc] = size (Z);
  if (isempty (X))
    X = 1:zc;
    Y = 1:zr;
  endif
  if (! isnumeric (X) || ! isnumeric (Y))
    error ("interp2: X, Y must be numeric matrices");
  endif
  if (! isempty (n))
    ## Calculate the interleaved input vectors.
    p = 2^n;
    XI = (p:p*zc)/p;
    YI = (p:p*zr)'/p;
  endif
  if (! isnumeric (XI) || ! isnumeric (YI))
    error ("interp2: XI, YI must be numeric");
  endif


  if (strcmp (method, "linear") || strcmp (method, "nearest") ...
      || strcmp (method, "pchip"))

    ## If X and Y vectors produce a grid from them
    if (isvector (X) && isvector (Y))
      X = X(:); Y = Y(:);
    elseif (size_equal (X, Y))
      X = X(1,:)'; Y = Y(:,1);
    else
      error ("interp2: X and Y must be matrices of same size");
    endif
    if (columns (Z) != length (X) || rows (Z) != length (Y))
      error ("interp2: X and Y size must match the dimensions of Z");
    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))
      [XI, YI] = meshgrid (XI, YI);
    elseif (! size_equal (XI, YI))
      error ("interp2: XI and YI must be matrices of equal size");
    endif

    ## if XI, YI are vectors, X and Y should share their orientation.
    if (rows (XI) == 1)
      if (rows (X) != 1)
        X = X.';
      endif
      if (rows (Y) != 1)
        Y = Y.';
      endif
    elseif (columns (XI) == 1)
      if (columns (X) != 1)
        X = X.';
      endif
      if (columns (Y) != 1)
        Y = Y.';
      endif
    endif

    xidx = lookup (X, XI, "lr");
    yidx = lookup (Y, YI, "lr");

    if (strcmp (method, "linear"))
      ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy
      ##
      ## a-b
      ## | |
      ## c-d
      a = Z(1:(zr - 1), 1:(zc - 1));
      b = Z(1:(zr - 1), 2:zc) - a;
      c = Z(2:zr, 1:(zc - 1)) - a;
      d = Z(2:zr, 2:zc) - a - b - c;

      ## scale XI, YI values to a 1-spaced grid
      Xsc = (XI - X(xidx)) ./ (diff (X)(xidx));
      Ysc = (YI - Y(yidx)) ./ (diff (Y)(yidx));

      ## Get 2D index.
      idx = sub2ind (size (a), yidx, xidx);
      ## We can dispose of the 1D indices at this point to save memory.
      clear xidx yidx;

      ## apply plane equation
      ZI = a(idx) + b(idx).*Xsc + c(idx).*Ysc + d(idx).*Xsc.*Ysc;

    elseif (strcmp (method, "nearest"))
      ii = (XI - X(xidx) >= X(xidx + 1) - XI);
      jj = (YI - Y(yidx) >= Y(yidx + 1) - YI);
      idx = sub2ind (size (Z), yidx+jj, xidx+ii);
      ZI = Z(idx);

    elseif (strcmp (method, "pchip"))

      if (length (X) < 2 || length (Y) < 2)
        error ("interp2: pchip2 requires at least 2 points in each dimension");
      endif

      ## first order derivatives
      DX = __pchip_deriv__ (X, Z, 2);
      DY = __pchip_deriv__ (Y, Z, 1);
      ## Compute mixed derivatives row-wise and column-wise, use the average.
      DXY = (__pchip_deriv__ (X, DY, 2) + __pchip_deriv__ (Y, DX, 1))/2;

      ## do the bicubic interpolation
      hx = diff (X); hx = hx(xidx);
      hy = diff (Y); hy = hy(yidx);

      tx = (XI - X(xidx)) ./ hx;
      ty = (YI - Y(yidx)) ./ hy;

      ## construct the cubic hermite base functions in x, y

      ## formulas:
      ## b{1,1} =    ( 2*t.^3 - 3*t.^2     + 1);
      ## b{2,1} = h.*(   t.^3 - 2*t.^2 + t    );
      ## b{1,2} =    (-2*t.^3 + 3*t.^2        );
      ## b{2,2} = h.*(   t.^3 -   t.^2        );

      ## optimized equivalents of the above:
      t1 = tx.^2;
      t2 = tx.*t1 - t1;
      xb{2,2} = hx.*t2;
      t1 = t2 - t1;
      xb{2,1} = hx.*(t1 + tx);
      t2 += t1;
      xb{1,2} = -t2;
      xb{1,1} = t2 + 1;

      t1 = ty.^2;
      t2 = ty.*t1 - t1;
      yb{2,2} = hy.*t2;
      t1 = t2 - t1;
      yb{2,1} = hy.*(t1 + ty);
      t2 += t1;
      yb{1,2} = -t2;
      yb{1,1} = t2 + 1;

      ZI = zeros (size (XI));
      for i = 1:2
        for j = 1:2
          zidx = sub2ind (size (Z), yidx+(j-1), xidx+(i-1));
          ZI += xb{1,i} .* yb{1,j} .*   Z(zidx);
          ZI += xb{2,i} .* yb{1,j} .*  DX(zidx);
          ZI += xb{1,i} .* yb{2,j} .*  DY(zidx);
          ZI += xb{2,i} .* yb{2,j} .* DXY(zidx);
        endfor
      endfor

    endif

    if (! isempty (extrapval))
      ## set points outside the table to 'extrapval'
      if (X (1) < X (end))
        if (Y (1) < Y (end))
          ZI (XI < X(1,1) | XI > X(end) | YI < Y(1,1) | YI > Y(end)) = ...
                  extrapval;
        else
          ZI (XI < X(1) | XI > X(end) | YI < Y(end) | YI > Y(1)) = ...
                  extrapval;
        endif
      else
        if (Y (1) < Y (end))
          ZI (XI < X(end) | XI > X(1) | YI < Y(1) | YI > Y(end)) = ...
                  extrapval;
        else
          ZI (XI < X(1,end) | XI > X(1) | YI < Y(end) | YI > Y(1)) = ...
                  extrapval;
        endif
      endif
    endif

  else

    ## Check dimensions of X and Y
    if (isvector (X) && isvector (Y))
      X = X(:).';
      Y = Y(:);
      if (!isequal ([length(Y), length(X)], size(Z)))
        error ("interp2: X and Y size must match the dimensions of Z");
      endif
    elseif (!size_equal (X, Y))
      error ("interp2: X and Y must be matrices of equal size");
      if (! size_equal (X, Z))
        error ("interp2: X and Y size must match the dimensions of Z");
      endif
    endif

    ## Check dimensions of XI and YI
    if (isvector (XI) && isvector (YI) && ! size_equal (XI, YI))
      XI = XI(:).';
      YI = YI(:);
      [XI, YI] = meshgrid (XI, YI);
    elseif (! size_equal (XI, YI))
      error ("interp2: XI and YI must be matrices of equal size");
    endif

    if (strcmp (method, "cubic"))
      if (isgriddata (XI) && isgriddata (YI'))
        ZI = bicubic (X, Y, Z, XI (1, :), YI (:, 1), extrapval);
      elseif (isgriddata (X) && isgriddata (Y'))
        ## Allocate output
        ZI = zeros (size (X));

        ## Find inliers
        inside = !(XI < X (1) | XI > X (end) | YI < Y (1) | YI > Y (end));

        ## Scale XI and YI to match indices of Z
        XI = (columns (Z) - 1) * (XI - X (1)) / (X (end) - X (1)) + 1;
        YI = (rows (Z) - 1) * (YI - Y (1)) / (Y (end) - Y (1)) + 1;

        ## Start the real work
        K = floor (XI);
        L = floor (YI);

        ## Coefficients
        AY1  = bc ((YI - L + 1));
        AX1  = bc ((XI - K + 1));
        AY0  = bc ((YI - L + 0));
        AX0  = bc ((XI - K + 0));
        AY_1 = bc ((YI - L - 1));
        AX_1 = bc ((XI - K - 1));
        AY_2 = bc ((YI - L - 2));
        AX_2 = bc ((XI - K - 2));

        ## Perform interpolation
        sz = size(Z);
        ZI = AY_2 .* AX_2 .* Z (sym_sub2ind (sz, L+2, K+2)) ...
           + AY_2 .* AX_1 .* Z (sym_sub2ind (sz, L+2, K+1)) ...
           + AY_2 .* AX0  .* Z (sym_sub2ind (sz, L+2, K))   ...
           + AY_2 .* AX1  .* Z (sym_sub2ind (sz, L+2, K-1)) ...
           + AY_1 .* AX_2 .* Z (sym_sub2ind (sz, L+1, K+2)) ...
           + AY_1 .* AX_1 .* Z (sym_sub2ind (sz, L+1, K+1)) ...
           + AY_1 .* AX0  .* Z (sym_sub2ind (sz, L+1, K))   ...
           + AY_1 .* AX1  .* Z (sym_sub2ind (sz, L+1, K-1)) ...
           + AY0  .* AX_2 .* Z (sym_sub2ind (sz, L,   K+2)) ...
           + AY0  .* AX_1 .* Z (sym_sub2ind (sz, L,   K+1)) ...
           + AY0  .* AX0  .* Z (sym_sub2ind (sz, L,   K))   ...
           + AY0  .* AX1  .* Z (sym_sub2ind (sz, L,   K-1)) ...
           + AY1  .* AX_2 .* Z (sym_sub2ind (sz, L-1, K+2)) ...
           + AY1  .* AX_1 .* Z (sym_sub2ind (sz, L-1, K+1)) ...
           + AY1  .* AX0  .* Z (sym_sub2ind (sz, L-1, K))   ...
           + AY1  .* AX1  .* Z (sym_sub2ind (sz, L-1, K-1));
        ZI (!inside) = extrapval;

      else
        error ("interp2: input data must have `meshgrid' format");
      endif

    elseif (strcmp (method, "spline"))
      if (isgriddata (XI) && isgriddata (YI'))
        ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrapval,
                        "spline");
      else
        error ("interp2: input data must have `meshgrid' format");
      endif
    else
      error ("interp2: interpolation METHOD not recognized");
    endif

  endif
endfunction

function b = isgriddata (X)
  d1 = diff (X, 1, 1);
  b = all (d1 (:) == 0);
endfunction

## Compute the bicubic interpolation coefficients
function o = bc(x)
  x = abs(x);
  o = zeros(size(x));
  idx1 = (x < 1);
  idx2 = !idx1 & (x < 2);
  o(idx1) = 1 - 2.*x(idx1).^2 + x(idx1).^3;
  o(idx2) = 4 - 8.*x(idx2) + 5.*x(idx2).^2 - x(idx2).^3;
endfunction

## This version of sub2ind behaves as if the data was symmetrically padded
function ind = sym_sub2ind(sz, Y, X)
  Y (Y < 1) = 1 - Y (Y < 1);
  while (any (Y (:) > 2 * sz (1)))
    Y (Y > 2 * sz (1)) = round (Y (Y > 2 * sz (1)) / 2);
  endwhile
  Y (Y > sz (1)) = 1 + 2 * sz (1) - Y (Y > sz (1));
  X (X < 1) = 1 - X (X < 1);
  while (any (X (:) > 2 * sz (2)))
    X (X > 2 * sz (2)) = round (X (X > 2 * sz (2)) / 2);
  endwhile
  X (X > sz (2)) = 1 + 2 * sz (2) - X (X > sz (2));
  ind = sub2ind(sz, Y, X);
endfunction


%!demo
%! A=[13,-1,12;5,4,3;1,6,2];
%! x=[0,1,4]; 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,'linear'));
%! [x,y] = meshgrid(x,y);
%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;

%!demo
%! [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,'linear'));
%! [x,y] = meshgrid(x,y);
%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;

%!demo
%! A=[13,-1,12;5,4,3;1,6,2];
%! x=[0,1,4]; 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,'nearest'));
%! [x,y] = meshgrid(x,y);
%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;

%!demo
%! [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,'nearest'));
%! [x,y] = meshgrid(x,y);
%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;

%!demo
%! 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,'pchip'));
%! [x,y] = meshgrid(x,y);
%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;

%!demo
%! [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,'pchip'));
%! [x,y] = meshgrid(x,y);
%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;

%!demo
%! 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;

%!demo
%! [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;

%!demo
%! 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;

%!demo
%! [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,'spline'));
%! [x,y] = meshgrid(x,y);
%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;

%!test % simple test
%!  x = [1,2,3];
%!  y = [4,5,6,7];
%!  [X, Y] = meshgrid(x,y);
%!  Orig = X.^2 + Y.^3;
%!  xi = [1.2,2, 1.5];
%!  yi = [6.2, 4.0, 5.0]';
%!
%!  Expected = ...
%!    [243,   245.4,  243.9;
%!      65.6,  68,     66.5;
%!     126.6, 129,    127.5];
%!  Result = interp2(x,y,Orig, xi, yi);
%!
%!  assert(Result, Expected, 1000*eps);

%!test % 2^n form
%!  x = [1,2,3];
%!  y = [4,5,6,7];
%!  [X, Y] = meshgrid(x,y);
%!  Orig = X.^2 + Y.^3;
%!  xi = [1:0.25:3]; yi = [4:0.25:7]';
%!  Expected = interp2(x,y,Orig, xi, yi);
%!  Result = interp2(Orig,2);
%!
%!  assert(Result, Expected, 10*eps);

%!test % matrix slice
%!  A = eye(4);
%!  assert(interp2(A,[1:4],[1:4]),[1,1,1,1]);

%!test % non-gridded XI,YI
%!  A = eye(4);
%!  assert(interp2(A,[1,2;3,4],[1,3;2,4]),[1,0;0,1]);

%!test % for values outside of boundaries
%!  x = [1,2,3];
%!  y = [4,5,6,7];
%!  [X, Y] = meshgrid(x,y);
%!  Orig = X.^2 + Y.^3;
%!  xi = [0,4];
%!  yi = [3,8]';
%!  assert(interp2(x,y,Orig, xi, yi),[NA,NA;NA,NA]);
%!  assert(interp2(x,y,Orig, xi, yi,'linear', 0),[0,0;0,0]);

%!test % for values at boundaries
%!  A=[1,2;3,4];
%!  x=[0,1];
%!  y=[2,3]';
%!  assert(interp2(x,y,A,x,y,'linear'), A);
%!  assert(interp2(x,y,A,x,y,'nearest'), A);

%!test % for Matlab-compatible rounding for 'nearest'
%! X = meshgrid (1:4);
%! assert (interp2 (X, 2.5, 2.5, 'nearest'), 3);

%!shared z, zout, tol
%!  z = [1 3 5; 3 5 7; 5 7 9];
%!  zout = [1 2 3 4 5; 2 3 4 5 6; 3 4 5 6 7; 4 5 6 7 8; 5 6 7 8 9];
%!  tol = 2 * eps;
%!assert (interp2 (z), zout, tol);
%!assert (interp2 (z, "linear"), zout, tol);
%!assert (interp2 (z, "pchip"), zout, tol);
%!assert (interp2 (z, "cubic"), zout, 10 * tol);
%!assert (interp2 (z, "spline"), zout, tol);
%!assert (interp2 (z, [2 3 1], [2 2 2]', "linear"), repmat ([5, 7, 3], [3, 1]), tol) 
%!assert (interp2 (z, [2 3 1], [2 2 2]', "pchip"), repmat ([5, 7, 3], [3, 1]), tol) 
%!assert (interp2 (z, [2 3 1], [2 2 2]', "cubic"), repmat ([5, 7, 3], [3, 1]), 10 * tol) 
%!assert (interp2 (z, [2 3 1], [2 2 2]', "spline"), repmat ([5, 7, 3], [3, 1]), tol) 
%!assert (interp2 (z, [2 3 1], [2 2 2], "linear"), [5 7 3], tol);
%!assert (interp2 (z, [2 3 1], [2 2 2], "pchip"), [5 7 3], tol);
%!assert (interp2 (z, [2 3 1], [2 2 2], "cubic"), [5 7 3], 10 * tol);
%!assert (interp2 (z, [2 3 1], [2 2 2], "spline"), [5 7 3], tol);