octave_packages/m/general/interp2.m

   1 ## Copyright (C) 2000-2012 Kai Habel
   2 ## Copyright (C) 2009 Jaroslav Hajek
   3 ##
   4 ## This file is part of Octave.
   5 ##
   6 ## Octave is free software; you can redistribute it and/or modify it
   7 ## under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation; either version 3 of the License, or (at
   9 ## your option) any later version.
  10 ##
  11 ## Octave is distributed in the hope that it will be useful, but
  12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ## General Public License for more details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License
  17 ## along with Octave; see the file COPYING.  If not, see
  18 ## <http://www.gnu.org/licenses/>.
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn  {Function File} {@var{zi} =} interp2 (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi})
  22 ## @deftypefnx {Function File} {@var{zi} =} interp2 (@var{Z}, @var{xi}, @var{yi})
  23 ## @deftypefnx {Function File} {@var{zi} =} interp2 (@var{Z}, @var{n})
  24 ## @deftypefnx {Function File} {@var{zi} =} interp2 (@dots{}, @var{method})
  25 ## @deftypefnx {Function File} {@var{zi} =} interp2 (@dots{}, @var{method}, @var{extrapval})
  26 ##
  27 ## Two-dimensional interpolation.  @var{x}, @var{y} and @var{z} describe a
  28 ## surface function.  If @var{x} and @var{y} are vectors their length
  29 ## must correspondent to the size of @var{z}.  @var{x} and @var{y} must be
  30 ## monotonic.  If they are matrices they must have the @code{meshgrid}
  31 ## format.
  32 ##
  33 ## @table @code
  34 ## @item interp2 (@var{x}, @var{y}, @var{Z}, @var{xi}, @var{yi}, @dots{})
  35 ## Returns a matrix corresponding to the points described by the
  36 ## matrices @var{xi}, @var{yi}.
  37 ##
  38 ## If the last argument is a string, the interpolation method can
  39 ## be specified.  The method can be 'linear', 'nearest' or 'cubic'.
  40 ## If it is omitted 'linear' interpolation is assumed.
  41 ##
  42 ## @item interp2 (@var{z}, @var{xi}, @var{yi})
  43 ## Assumes @code{@var{x} = 1:rows (@var{z})} and @code{@var{y} =
  44 ## 1:columns (@var{z})}
  45 ##
  46 ## @item interp2 (@var{z}, @var{n})
  47 ## Interleaves the matrix @var{z} n-times.  If @var{n} is omitted a value
  48 ## of @code{@var{n} = 1} is assumed.
  49 ## @end table
  50 ##
  51 ## The variable @var{method} defines the method to use for the
  52 ## interpolation.  It can take one of the following values
  53 ##
  54 ## @table @asis
  55 ## @item 'nearest'
  56 ## Return the nearest neighbor.
  57 ##
  58 ## @item 'linear'
  59 ## Linear interpolation from nearest neighbors.
  60 ##
  61 ## @item 'pchip'
  62 ## Piecewise cubic Hermite interpolating polynomial.
  63 ##
  64 ## @item 'cubic'
  65 ## Cubic interpolation from four nearest neighbors.
  66 ##
  67 ## @item 'spline'
  68 ## Cubic spline interpolation---smooth first and second derivatives
  69 ## throughout the curve.
  70 ## @end table
  71 ##
  72 ## If a scalar value @var{extrapval} is defined as the final value, then
  73 ## values outside the mesh as set to this value.  Note that in this case
  74 ## @var{method} must be defined as well.  If @var{extrapval} is not
  75 ## defined then NA is assumed.
  76 ##
  77 ## @seealso{interp1}
  78 ## @end deftypefn
  79
  80 ## Author:      Kai Habel <kai.habel@gmx.de>
  81 ## 2005-03-02 Thomas Weber <weber@num.uni-sb.de>
  82 ##     * Add test cases
  83 ## 2005-03-02 Paul Kienzle <pkienzle@users.sf.net>
  84 ##     * Simplify
  85 ## 2005-04-23 Dmitri A. Sergatskov <dasergatskov@gmail.com>
  86 ##     * Modified demo and test for new gnuplot interface
  87 ## 2005-09-07 Hoxide <hoxide_dirac@yahoo.com.cn>
  88 ##     * Add bicubic interpolation method
  89 ##     * Fix the eat line bug when the last element of XI or YI is
  90 ##       negative or zero.
  91 ## 2005-11-26 Pierre Baldensperger <balden@libertysurf.fr>
  92 ##     * Rather big modification (XI,YI no longer need to be
  93 ##       "meshgridded") to be consistent with the help message
  94 ##       above and for compatibility.
  95
  96 function ZI = interp2 (varargin)
  97   Z = X = Y = XI = YI = n = [];
  98   method = "linear";
  99   extrapval = NA;
 100
 101   switch (nargin)
 102     case 1
 103       Z = varargin{1};
 104       n = 1;
 105     case 2
 106       if (ischar (varargin{2}))
 107         [Z, method] = deal (varargin{:});
 108         n = 1;
 109       else
 110         [Z, n] = deal (varargin{:});
 111       endif
 112     case 3
 113       if (ischar (varargin{3}))
 114         [Z, n, method] = deal (varargin{:});
 115       else
 116         [Z, XI, YI] = deal (varargin{:});
 117       endif
 118     case 4
 119       if (ischar (varargin{4}))
 120         [Z, XI, YI, method] = deal (varargin{:});
 121       else
 122         [Z, n, method, extrapval] = deal (varargin{:});
 123       endif
 124     case 5
 125       if (ischar (varargin{4}))
 126         [Z, XI, YI, method, extrapval] = deal (varargin{:});
 127       else
 128         [X, Y, Z, XI, YI] = deal (varargin{:});
 129       endif
 130     case 6
 131         [X, Y, Z, XI, YI, method] = deal (varargin{:});
 132     case 7
 133         [X, Y, Z, XI, YI, method, extrapval] = deal (varargin{:});
 134     otherwise
 135       print_usage ();
 136   endswitch
 137
 138   ## Type checking.
 139   if (!ismatrix (Z))
 140     error ("interp2: Z must be a matrix");
 141   endif
 142   if (!isempty (n) && !isscalar (n))
 143     error ("interp2: N must be a scalar");
 144   endif
 145   if (!ischar (method))
 146     error ("interp2: METHOD must be a string");
 147   endif
 148   if (ischar (extrapval) || strcmp (extrapval, "extrap"))
 149     extrapval = [];
 150   elseif (!isscalar (extrapval))
 151     error ("interp2: EXTRAPVAL must be a scalar");
 152   endif
 153
 154   ## Define X, Y, XI, YI if needed
 155   [zr, zc] = size (Z);
 156   if (isempty (X))
 157     X = 1:zc;
 158     Y = 1:zr;
 159   endif
 160   if (! isnumeric (X) || ! isnumeric (Y))
 161     error ("interp2: X, Y must be numeric matrices");
 162   endif
 163   if (! isempty (n))
 164     ## Calculate the interleaved input vectors.
 165     p = 2^n;
 166     XI = (p:p*zc)/p;
 167     YI = (p:p*zr)'/p;
 168   endif
 169   if (! isnumeric (XI) || ! isnumeric (YI))
 170     error ("interp2: XI, YI must be numeric");
 171   endif
 172
 173
 174   if (strcmp (method, "linear") || strcmp (method, "nearest") ...
 175       || strcmp (method, "pchip"))
 176
 177     ## If X and Y vectors produce a grid from them
 178     if (isvector (X) && isvector (Y))
 179       X = X(:); Y = Y(:);
 180     elseif (size_equal (X, Y))
 181       X = X(1,:)'; Y = Y(:,1);
 182     else
 183       error ("interp2: X and Y must be matrices of same size");
 184     endif
 185     if (columns (Z) != length (X) || rows (Z) != length (Y))
 186       error ("interp2: X and Y size must match the dimensions of Z");
 187     endif
 188
 189     ## If Xi and Yi are vectors of different orientation build a grid
 190     if ((rows (XI) == 1 && columns (YI) == 1)
 191         || (columns (XI) == 1 && rows (YI) == 1))
 192       [XI, YI] = meshgrid (XI, YI);
 193     elseif (! size_equal (XI, YI))
 194       error ("interp2: XI and YI must be matrices of equal size");
 195     endif
 196
 197     ## if XI, YI are vectors, X and Y should share their orientation.
 198     if (rows (XI) == 1)
 199       if (rows (X) != 1)
 200         X = X.';
 201       endif
 202       if (rows (Y) != 1)
 203         Y = Y.';
 204       endif
 205     elseif (columns (XI) == 1)
 206       if (columns (X) != 1)
 207         X = X.';
 208       endif
 209       if (columns (Y) != 1)
 210         Y = Y.';
 211       endif
 212     endif
 213
 214     xidx = lookup (X, XI, "lr");
 215     yidx = lookup (Y, YI, "lr");
 216
 217     if (strcmp (method, "linear"))
 218       ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy
 219       ##
 220       ## a-b
 221       ## | |
 222       ## c-d
 223       a = Z(1:(zr - 1), 1:(zc - 1));
 224       b = Z(1:(zr - 1), 2:zc) - a;
 225       c = Z(2:zr, 1:(zc - 1)) - a;
 226       d = Z(2:zr, 2:zc) - a - b - c;
 227
 228       ## scale XI, YI values to a 1-spaced grid
 229       Xsc = (XI - X(xidx)) ./ (diff (X)(xidx));
 230       Ysc = (YI - Y(yidx)) ./ (diff (Y)(yidx));
 231
 232       ## Get 2D index.
 233       idx = sub2ind (size (a), yidx, xidx);
 234       ## We can dispose of the 1D indices at this point to save memory.
 235       clear xidx yidx;
 236
 237       ## apply plane equation
 238       ZI = a(idx) + b(idx).*Xsc + c(idx).*Ysc + d(idx).*Xsc.*Ysc;
 239
 240     elseif (strcmp (method, "nearest"))
 241       ii = (XI - X(xidx) >= X(xidx + 1) - XI);
 242       jj = (YI - Y(yidx) >= Y(yidx + 1) - YI);
 243       idx = sub2ind (size (Z), yidx+jj, xidx+ii);
 244       ZI = Z(idx);
 245
 246     elseif (strcmp (method, "pchip"))
 247
 248       if (length (X) < 2 || length (Y) < 2)
 249         error ("interp2: pchip2 requires at least 2 points in each dimension");
 250       endif
 251
 252       ## first order derivatives
 253       DX = __pchip_deriv__ (X, Z, 2);
 254       DY = __pchip_deriv__ (Y, Z, 1);
 255       ## Compute mixed derivatives row-wise and column-wise, use the average.
 256       DXY = (__pchip_deriv__ (X, DY, 2) + __pchip_deriv__ (Y, DX, 1))/2;
 257
 258       ## do the bicubic interpolation
 259       hx = diff (X); hx = hx(xidx);
 260       hy = diff (Y); hy = hy(yidx);
 261
 262       tx = (XI - X(xidx)) ./ hx;
 263       ty = (YI - Y(yidx)) ./ hy;
 264
 265       ## construct the cubic hermite base functions in x, y
 266
 267       ## formulas:
 268       ## b{1,1} =    ( 2*t.^3 - 3*t.^2     + 1);
 269       ## b{2,1} = h.*(   t.^3 - 2*t.^2 + t    );
 270       ## b{1,2} =    (-2*t.^3 + 3*t.^2        );
 271       ## b{2,2} = h.*(   t.^3 -   t.^2        );
 272
 273       ## optimized equivalents of the above:
 274       t1 = tx.^2;
 275       t2 = tx.*t1 - t1;
 276       xb{2,2} = hx.*t2;
 277       t1 = t2 - t1;
 278       xb{2,1} = hx.*(t1 + tx);
 279       t2 += t1;
 280       xb{1,2} = -t2;
 281       xb{1,1} = t2 + 1;
 282
 283       t1 = ty.^2;
 284       t2 = ty.*t1 - t1;
 285       yb{2,2} = hy.*t2;
 286       t1 = t2 - t1;
 287       yb{2,1} = hy.*(t1 + ty);
 288       t2 += t1;
 289       yb{1,2} = -t2;
 290       yb{1,1} = t2 + 1;
 291
 292       ZI = zeros (size (XI));
 293       for i = 1:2
 294         for j = 1:2
 295           zidx = sub2ind (size (Z), yidx+(j-1), xidx+(i-1));
 296           ZI += xb{1,i} .* yb{1,j} .*   Z(zidx);
 297           ZI += xb{2,i} .* yb{1,j} .*  DX(zidx);
 298           ZI += xb{1,i} .* yb{2,j} .*  DY(zidx);
 299           ZI += xb{2,i} .* yb{2,j} .* DXY(zidx);
 300         endfor
 301       endfor
 302
 303     endif
 304
 305     if (! isempty (extrapval))
 306       ## set points outside the table to 'extrapval'
 307       if (X (1) < X (end))
 308         if (Y (1) < Y (end))
 309           ZI (XI < X(1,1) | XI > X(end) | YI < Y(1,1) | YI > Y(end)) = ...
 310                   extrapval;
 311         else
 312           ZI (XI < X(1) | XI > X(end) | YI < Y(end) | YI > Y(1)) = ...
 313                   extrapval;
 314         endif
 315       else
 316         if (Y (1) < Y (end))
 317           ZI (XI < X(end) | XI > X(1) | YI < Y(1) | YI > Y(end)) = ...
 318                   extrapval;
 319         else
 320           ZI (XI < X(1,end) | XI > X(1) | YI < Y(end) | YI > Y(1)) = ...
 321                   extrapval;
 322         endif
 323       endif
 324     endif
 325
 326   else
 327
 328     ## Check dimensions of X and Y
 329     if (isvector (X) && isvector (Y))
 330       X = X(:).';
 331       Y = Y(:);
 332       if (!isequal ([length(Y), length(X)], size(Z)))
 333         error ("interp2: X and Y size must match the dimensions of Z");
 334       endif
 335     elseif (!size_equal (X, Y))
 336       error ("interp2: X and Y must be matrices of equal size");
 337       if (! size_equal (X, Z))
 338         error ("interp2: X and Y size must match the dimensions of Z");
 339       endif
 340     endif
 341
 342     ## Check dimensions of XI and YI
 343     if (isvector (XI) && isvector (YI) && ! size_equal (XI, YI))
 344       XI = XI(:).';
 345       YI = YI(:);
 346       [XI, YI] = meshgrid (XI, YI);
 347     elseif (! size_equal (XI, YI))
 348       error ("interp2: XI and YI must be matrices of equal size");
 349     endif
 350
 351     if (strcmp (method, "cubic"))
 352       if (isgriddata (XI) && isgriddata (YI'))
 353         ZI = bicubic (X, Y, Z, XI (1, :), YI (:, 1), extrapval);
 354       elseif (isgriddata (X) && isgriddata (Y'))
 355         ## Allocate output
 356         ZI = zeros (size (X));
 357
 358         ## Find inliers
 359         inside = !(XI < X (1) | XI > X (end) | YI < Y (1) | YI > Y (end));
 360
 361         ## Scale XI and YI to match indices of Z
 362         XI = (columns (Z) - 1) * (XI - X (1)) / (X (end) - X (1)) + 1;
 363         YI = (rows (Z) - 1) * (YI - Y (1)) / (Y (end) - Y (1)) + 1;
 364
 365         ## Start the real work
 366         K = floor (XI);
 367         L = floor (YI);
 368
 369         ## Coefficients
 370         AY1  = bc ((YI - L + 1));
 371         AX1  = bc ((XI - K + 1));
 372         AY0  = bc ((YI - L + 0));
 373         AX0  = bc ((XI - K + 0));
 374         AY_1 = bc ((YI - L - 1));
 375         AX_1 = bc ((XI - K - 1));
 376         AY_2 = bc ((YI - L - 2));
 377         AX_2 = bc ((XI - K - 2));
 378
 379         ## Perform interpolation
 380         sz = size(Z);
 381         ZI = AY_2 .* AX_2 .* Z (sym_sub2ind (sz, L+2, K+2)) ...
 382            + AY_2 .* AX_1 .* Z (sym_sub2ind (sz, L+2, K+1)) ...
 383            + AY_2 .* AX0  .* Z (sym_sub2ind (sz, L+2, K))   ...
 384            + AY_2 .* AX1  .* Z (sym_sub2ind (sz, L+2, K-1)) ...
 385            + AY_1 .* AX_2 .* Z (sym_sub2ind (sz, L+1, K+2)) ...
 386            + AY_1 .* AX_1 .* Z (sym_sub2ind (sz, L+1, K+1)) ...
 387            + AY_1 .* AX0  .* Z (sym_sub2ind (sz, L+1, K))   ...
 388            + AY_1 .* AX1  .* Z (sym_sub2ind (sz, L+1, K-1)) ...
 389            + AY0  .* AX_2 .* Z (sym_sub2ind (sz, L,   K+2)) ...
 390            + AY0  .* AX_1 .* Z (sym_sub2ind (sz, L,   K+1)) ...
 391            + AY0  .* AX0  .* Z (sym_sub2ind (sz, L,   K))   ...
 392            + AY0  .* AX1  .* Z (sym_sub2ind (sz, L,   K-1)) ...
 393            + AY1  .* AX_2 .* Z (sym_sub2ind (sz, L-1, K+2)) ...
 394            + AY1  .* AX_1 .* Z (sym_sub2ind (sz, L-1, K+1)) ...
 395            + AY1  .* AX0  .* Z (sym_sub2ind (sz, L-1, K))   ...
 396            + AY1  .* AX1  .* Z (sym_sub2ind (sz, L-1, K-1));
 397         ZI (!inside) = extrapval;
 398
 399       else
 400         error ("interp2: input data must have `meshgrid' format");
 401       endif
 402
 403     elseif (strcmp (method, "spline"))
 404       if (isgriddata (XI) && isgriddata (YI'))
 405         ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrapval,
 406                         "spline");
 407       else
 408         error ("interp2: input data must have `meshgrid' format");
 409       endif
 410     else
 411       error ("interp2: interpolation METHOD not recognized");
 412     endif
 413
 414   endif
 415 endfunction
 416
 417 function b = isgriddata (X)
 418   d1 = diff (X, 1, 1);
 419   b = all (d1 (:) == 0);
 420 endfunction
 421
 422 ## Compute the bicubic interpolation coefficients
 423 function o = bc(x)
 424   x = abs(x);
 425   o = zeros(size(x));
 426   idx1 = (x < 1);
 427   idx2 = !idx1 & (x < 2);
 428   o(idx1) = 1 - 2.*x(idx1).^2 + x(idx1).^3;
 429   o(idx2) = 4 - 8.*x(idx2) + 5.*x(idx2).^2 - x(idx2).^3;
 430 endfunction
 431
 432 ## This version of sub2ind behaves as if the data was symmetrically padded
 433 function ind = sym_sub2ind(sz, Y, X)
 434   Y (Y < 1) = 1 - Y (Y < 1);
 435   while (any (Y (:) > 2 * sz (1)))
 436     Y (Y > 2 * sz (1)) = round (Y (Y > 2 * sz (1)) / 2);
 437   endwhile
 438   Y (Y > sz (1)) = 1 + 2 * sz (1) - Y (Y > sz (1));
 439   X (X < 1) = 1 - X (X < 1);
 440   while (any (X (:) > 2 * sz (2)))
 441     X (X > 2 * sz (2)) = round (X (X > 2 * sz (2)) / 2);
 442   endwhile
 443   X (X > sz (2)) = 1 + 2 * sz (2) - X (X > sz (2));
 444   ind = sub2ind(sz, Y, X);
 445 endfunction
 446
 447
 448 %!demo
 449 %! A=[13,-1,12;5,4,3;1,6,2];
 450 %! x=[0,1,4]; y=[10,11,12];
 451 %! xi=linspace(min(x),max(x),17);
 452 %! yi=linspace(min(y),max(y),26)';
 453 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'linear'));
 454 %! [x,y] = meshgrid(x,y);
 455 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 456
 457 %!demo
 458 %! [x,y,A] = peaks(10);
 459 %! x = x(1,:)'; y = y(:,1);
 460 %! xi=linspace(min(x),max(x),41);
 461 %! yi=linspace(min(y),max(y),41)';
 462 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'linear'));
 463 %! [x,y] = meshgrid(x,y);
 464 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 465
 466 %!demo
 467 %! A=[13,-1,12;5,4,3;1,6,2];
 468 %! x=[0,1,4]; y=[10,11,12];
 469 %! xi=linspace(min(x),max(x),17);
 470 %! yi=linspace(min(y),max(y),26)';
 471 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'nearest'));
 472 %! [x,y] = meshgrid(x,y);
 473 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 474
 475 %!demo
 476 %! [x,y,A] = peaks(10);
 477 %! x = x(1,:)'; y = y(:,1);
 478 %! xi=linspace(min(x),max(x),41);
 479 %! yi=linspace(min(y),max(y),41)';
 480 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'nearest'));
 481 %! [x,y] = meshgrid(x,y);
 482 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 483
 484 %!demo
 485 %! A=[13,-1,12;5,4,3;1,6,2];
 486 %! x=[0,1,2]; y=[10,11,12];
 487 %! xi=linspace(min(x),max(x),17);
 488 %! yi=linspace(min(y),max(y),26)';
 489 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'pchip'));
 490 %! [x,y] = meshgrid(x,y);
 491 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 492
 493 %!demo
 494 %! [x,y,A] = peaks(10);
 495 %! x = x(1,:)'; y = y(:,1);
 496 %! xi=linspace(min(x),max(x),41);
 497 %! yi=linspace(min(y),max(y),41)';
 498 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'pchip'));
 499 %! [x,y] = meshgrid(x,y);
 500 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 501
 502 %!demo
 503 %! A=[13,-1,12;5,4,3;1,6,2];
 504 %! x=[0,1,2]; y=[10,11,12];
 505 %! xi=linspace(min(x),max(x),17);
 506 %! yi=linspace(min(y),max(y),26)';
 507 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic'));
 508 %! [x,y] = meshgrid(x,y);
 509 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 510
 511 %!demo
 512 %! [x,y,A] = peaks(10);
 513 %! x = x(1,:)'; y = y(:,1);
 514 %! xi=linspace(min(x),max(x),41);
 515 %! yi=linspace(min(y),max(y),41)';
 516 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic'));
 517 %! [x,y] = meshgrid(x,y);
 518 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 519
 520 %!demo
 521 %! A=[13,-1,12;5,4,3;1,6,2];
 522 %! x=[0,1,2]; y=[10,11,12];
 523 %! xi=linspace(min(x),max(x),17);
 524 %! yi=linspace(min(y),max(y),26)';
 525 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'spline'));
 526 %! [x,y] = meshgrid(x,y);
 527 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 528
 529 %!demo
 530 %! [x,y,A] = peaks(10);
 531 %! x = x(1,:)'; y = y(:,1);
 532 %! xi=linspace(min(x),max(x),41);
 533 %! yi=linspace(min(y),max(y),41)';
 534 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'spline'));
 535 %! [x,y] = meshgrid(x,y);
 536 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 537
 538 %!test % simple test
 539 %!  x = [1,2,3];
 540 %!  y = [4,5,6,7];
 541 %!  [X, Y] = meshgrid(x,y);
 542 %!  Orig = X.^2 + Y.^3;
 543 %!  xi = [1.2,2, 1.5];
 544 %!  yi = [6.2, 4.0, 5.0]';
 545 %!
 546 %!  Expected = ...
 547 %!    [243,   245.4,  243.9;
 548 %!      65.6,  68,     66.5;
 549 %!     126.6, 129,    127.5];
 550 %!  Result = interp2(x,y,Orig, xi, yi);
 551 %!
 552 %!  assert(Result, Expected, 1000*eps);
 553
 554 %!test % 2^n form
 555 %!  x = [1,2,3];
 556 %!  y = [4,5,6,7];
 557 %!  [X, Y] = meshgrid(x,y);
 558 %!  Orig = X.^2 + Y.^3;
 559 %!  xi = [1:0.25:3]; yi = [4:0.25:7]';
 560 %!  Expected = interp2(x,y,Orig, xi, yi);
 561 %!  Result = interp2(Orig,2);
 562 %!
 563 %!  assert(Result, Expected, 10*eps);
 564
 565 %!test % matrix slice
 566 %!  A = eye(4);
 567 %!  assert(interp2(A,[1:4],[1:4]),[1,1,1,1]);
 568
 569 %!test % non-gridded XI,YI
 570 %!  A = eye(4);
 571 %!  assert(interp2(A,[1,2;3,4],[1,3;2,4]),[1,0;0,1]);
 572
 573 %!test % for values outside of boundaries
 574 %!  x = [1,2,3];
 575 %!  y = [4,5,6,7];
 576 %!  [X, Y] = meshgrid(x,y);
 577 %!  Orig = X.^2 + Y.^3;
 578 %!  xi = [0,4];
 579 %!  yi = [3,8]';
 580 %!  assert(interp2(x,y,Orig, xi, yi),[NA,NA;NA,NA]);
 581 %!  assert(interp2(x,y,Orig, xi, yi,'linear', 0),[0,0;0,0]);
 582
 583 %!test % for values at boundaries
 584 %!  A=[1,2;3,4];
 585 %!  x=[0,1];
 586 %!  y=[2,3]';
 587 %!  assert(interp2(x,y,A,x,y,'linear'), A);
 588 %!  assert(interp2(x,y,A,x,y,'nearest'), A);
 589
 590 %!test % for Matlab-compatible rounding for 'nearest'
 591 %! X = meshgrid (1:4);
 592 %! assert (interp2 (X, 2.5, 2.5, 'nearest'), 3);
 593
 594 %!shared z, zout, tol
 595 %!  z = [1 3 5; 3 5 7; 5 7 9];
 596 %!  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];
 597 %!  tol = 2 * eps;
 598 %!assert (interp2 (z), zout, tol);
 599 %!assert (interp2 (z, "linear"), zout, tol);
 600 %!assert (interp2 (z, "pchip"), zout, tol);
 601 %!assert (interp2 (z, "cubic"), zout, 10 * tol);
 602 %!assert (interp2 (z, "spline"), zout, tol);
 603 %!assert (interp2 (z, [2 3 1], [2 2 2]', "linear"), repmat ([5, 7, 3], [3, 1]), tol)
 604 %!assert (interp2 (z, [2 3 1], [2 2 2]', "pchip"), repmat ([5, 7, 3], [3, 1]), tol)
 605 %!assert (interp2 (z, [2 3 1], [2 2 2]', "cubic"), repmat ([5, 7, 3], [3, 1]), 10 * tol)
 606 %!assert (interp2 (z, [2 3 1], [2 2 2]', "spline"), repmat ([5, 7, 3], [3, 1]), tol)
 607 %!assert (interp2 (z, [2 3 1], [2 2 2], "linear"), [5 7 3], tol);
 608 %!assert (interp2 (z, [2 3 1], [2 2 2], "pchip"), [5 7 3], tol);
 609 %!assert (interp2 (z, [2 3 1], [2 2 2], "cubic"), [5 7 3], 10 * tol);
 610 %!assert (interp2 (z, [2 3 1], [2 2 2], "spline"), [5 7 3], tol);