/octave-3.6.2/scripts/general/interp2.m

# · Objective C · 610 lines · 547 code · 63 blank · 0 comment · 81 complexity · 7b378071460a9ccb2d26057d4b279319 MD5 · raw file 

    

  1. ## Copyright (C) 2000-2012 Kai Habel
    2. 

  2. ## Copyright (C) 2009 Jaroslav Hajek
    3. 

  3. ##
    4. 

  4. ## This file is part of Octave.
    5. 

  5. ##
    6. 

  6. ## Octave is free software; you can redistribute it and/or modify it
    7. 

  7. ## under the terms of the GNU General Public License as published by
    8. 

  8. ## the Free Software Foundation; either version 3 of the License, or (at
    9. 

  9. ## your option) any later version.
    10. 

  10. ##
    11. 

  11. ## Octave is distributed in the hope that it will be useful, but
    12. 

  12. ## WITHOUT ANY WARRANTY; without even the implied warranty of
    13. 

  13. ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    14. 

  14. ## General Public License for more details.
    15. 

  15. ##
    16. 

  16. ## You should have received a copy of the GNU General Public License
    17. 

  17. ## along with Octave; see the file COPYING.  If not, see
    18. 

  18. ## <http://www.gnu.org/licenses/>.
    19. 

  19. 

  20. ## -*- texinfo -*-
    21. 

  21. ## @deftypefn  {Function File} {@var{zi} =} interp2 (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi})
    22. 

  22. ## @deftypefnx {Function File} {@var{zi} =} interp2 (@var{Z}, @var{xi}, @var{yi})
    23. 

  23. ## @deftypefnx {Function File} {@var{zi} =} interp2 (@var{Z}, @var{n})
    24. 

  24. ## @deftypefnx {Function File} {@var{zi} =} interp2 (@dots{}, @var{method})
    25. 

  25. ## @deftypefnx {Function File} {@var{zi} =} interp2 (@dots{}, @var{method}, @var{extrapval})
    26. 

  26. ##
    27. 

  27. ## Two-dimensional interpolation.  @var{x}, @var{y} and @var{z} describe a
    28. 

  28. ## surface function.  If @var{x} and @var{y} are vectors their length
    29. 

  29. ## must correspondent to the size of @var{z}.  @var{x} and @var{y} must be
    30. 

  30. ## monotonic.  If they are matrices they must have the @code{meshgrid}
    31. 

  31. ## format.
    32. 

  32. ##
    33. 

  33. ## @table @code
    34. 

  34. ## @item interp2 (@var{x}, @var{y}, @var{Z}, @var{xi}, @var{yi}, @dots{})
    35. 

  35. ## Returns a matrix corresponding to the points described by the
    36. 

  36. ## matrices @var{xi}, @var{yi}.
    37. 

  37. ##
    38. 

  38. ## If the last argument is a string, the interpolation method can
    39. 

  39. ## be specified.  The method can be 'linear', 'nearest' or 'cubic'.
    40. 

  40. ## If it is omitted 'linear' interpolation is assumed.
    41. 

  41. ##
    42. 

  42. ## @item interp2 (@var{z}, @var{xi}, @var{yi})
    43. 

  43. ## Assumes @code{@var{x} = 1:rows (@var{z})} and @code{@var{y} =
    44. 

  44. ## 1:columns (@var{z})}
    45. 

  45. ##
    46. 

  46. ## @item interp2 (@var{z}, @var{n})
    47. 

  47. ## Interleaves the matrix @var{z} n-times.  If @var{n} is omitted a value
    48. 

  48. ## of @code{@var{n} = 1} is assumed.
    49. 

  49. ## @end table
    50. 

  50. ##
    51. 

  51. ## The variable @var{method} defines the method to use for the
    52. 

  52. ## interpolation.  It can take one of the following values
    53. 

  53. ##
    54. 

  54. ## @table @asis
    55. 

  55. ## @item 'nearest'
    56. 

  56. ## Return the nearest neighbor.
    57. 

  57. ##
    58. 

  58. ## @item 'linear'
    59. 

  59. ## Linear interpolation from nearest neighbors.
    60. 

  60. ##
    61. 

  61. ## @item 'pchip'
    62. 

  62. ## Piecewise cubic Hermite interpolating polynomial.
    63. 

  63. ##
    64. 

  64. ## @item 'cubic'
    65. 

  65. ## Cubic interpolation from four nearest neighbors.
    66. 

  66. ##
    67. 

  67. ## @item 'spline'
    68. 

  68. ## Cubic spline interpolation---smooth first and second derivatives
    69. 

  69. ## throughout the curve.
    70. 

  70. ## @end table
    71. 

  71. ##
    72. 

  72. ## If a scalar value @var{extrapval} is defined as the final value, then
    73. 

  73. ## values outside the mesh as set to this value.  Note that in this case
    74. 

  74. ## @var{method} must be defined as well.  If @var{extrapval} is not
    75. 

  75. ## defined then NA is assumed.
    76. 

  76. ##
    77. 

  77. ## @seealso{interp1}
    78. 

  78. ## @end deftypefn
    79. 

  79. 

  80. ## Author:      Kai Habel <kai.habel@gmx.de>
    81. 

  81. ## 2005-03-02 Thomas Weber <weber@num.uni-sb.de>
    82. 

  82. ##     * Add test cases
    83. 

  83. ## 2005-03-02 Paul Kienzle <pkienzle@users.sf.net>
    84. 

  84. ##     * Simplify
    85. 

  85. ## 2005-04-23 Dmitri A. Sergatskov <dasergatskov@gmail.com>
    86. 

  86. ##     * Modified demo and test for new gnuplot interface
    87. 

  87. ## 2005-09-07 Hoxide <hoxide_dirac@yahoo.com.cn>
    88. 

  88. ##     * Add bicubic interpolation method
    89. 

  89. ##     * Fix the eat line bug when the last element of XI or YI is
    90. 

  90. ##       negative or zero.
    91. 

  91. ## 2005-11-26 Pierre Baldensperger <balden@libertysurf.fr>
    92. 

  92. ##     * Rather big modification (XI,YI no longer need to be
    93. 

  93. ##       "meshgridded") to be consistent with the help message
    94. 

  94. ##       above and for compatibility.
    95. 

  95. 

  96. function ZI = interp2 (varargin)
    97. 

  97.   Z = X = Y = XI = YI = n = [];
    98. 

  98.   method = "linear";
    99. 

  99.   extrapval = NA;
    100. 

  100. 

  101.   switch (nargin)
    102. 

  102.     case 1
    103. 

  103.       Z = varargin{1};
    104. 

  104.       n = 1;
    105. 

  105.     case 2
    106. 

  106.       if (ischar (varargin{2}))
    107. 

  107.         [Z, method] = deal (varargin{:});
    108. 

  108.         n = 1;
    109. 

  109.       else
    110. 

  110.         [Z, n] = deal (varargin{:});
    111. 

  111.       endif
    112. 

  112.     case 3
    113. 

  113.       if (ischar (varargin{3}))
    114. 

  114.         [Z, n, method] = deal (varargin{:});
    115. 

  115.       else
    116. 

  116.         [Z, XI, YI] = deal (varargin{:});
    117. 

  117.       endif
    118. 

  118.     case 4
    119. 

  119.       if (ischar (varargin{4}))
    120. 

  120.         [Z, XI, YI, method] = deal (varargin{:});
    121. 

  121.       else
    122. 

  122.         [Z, n, method, extrapval] = deal (varargin{:});
    123. 

  123.       endif
    124. 

  124.     case 5
    125. 

  125.       if (ischar (varargin{4}))
    126. 

  126.         [Z, XI, YI, method, extrapval] = deal (varargin{:});
    127. 

  127.       else
    128. 

  128.         [X, Y, Z, XI, YI] = deal (varargin{:});
    129. 

  129.       endif
    130. 

  130.     case 6
    131. 

  131.         [X, Y, Z, XI, YI, method] = deal (varargin{:});
    132. 

  132.     case 7
    133. 

  133.         [X, Y, Z, XI, YI, method, extrapval] = deal (varargin{:});
    134. 

  134.     otherwise
    135. 

  135.       print_usage ();
    136. 

  136.   endswitch
    137. 

  137. 

  138.   ## Type checking.
    139. 

  139.   if (!ismatrix (Z))
    140. 

  140.     error ("interp2: Z must be a matrix");
    141. 

  141.   endif
    142. 

  142.   if (!isempty (n) && !isscalar (n))
    143. 

  143.     error ("interp2: N must be a scalar");
    144. 

  144.   endif
    145. 

  145.   if (!ischar (method))
    146. 

  146.     error ("interp2: METHOD must be a string");
    147. 

  147.   endif
    148. 

  148.   if (ischar (extrapval) || strcmp (extrapval, "extrap"))
    149. 

  149.     extrapval = [];
    150. 

  150.   elseif (!isscalar (extrapval))
    151. 

  151.     error ("interp2: EXTRAPVAL must be a scalar");
    152. 

  152.   endif
    153. 

  153. 

  154.   ## Define X, Y, XI, YI if needed
    155. 

  155.   [zr, zc] = size (Z);
    156. 

  156.   if (isempty (X))
    157. 

  157.     X = 1:zc;
    158. 

  158.     Y = 1:zr;
    159. 

  159.   endif
    160. 

  160.   if (! isnumeric (X) || ! isnumeric (Y))
    161. 

  161.     error ("interp2: X, Y must be numeric matrices");
    162. 

  162.   endif
    163. 

  163.   if (! isempty (n))
    164. 

  164.     ## Calculate the interleaved input vectors.
    165. 

  165.     p = 2^n;
    166. 

  166.     XI = (p:p*zc)/p;
    167. 

  167.     YI = (p:p*zr)'/p;
    168. 

  168.   endif
    169. 

  169.   if (! isnumeric (XI) || ! isnumeric (YI))
    170. 

  170.     error ("interp2: XI, YI must be numeric");
    171. 

  171.   endif
    172. 

  172. 

  173. 

  174.   if (strcmp (method, "linear") || strcmp (method, "nearest") ...
    175. 

  175.       || strcmp (method, "pchip"))
    176. 

  176. 

  177.     ## If X and Y vectors produce a grid from them
    178. 

  178.     if (isvector (X) && isvector (Y))
    179. 

  179.       X = X(:); Y = Y(:);
    180. 

  180.     elseif (size_equal (X, Y))
    181. 

  181.       X = X(1,:)'; Y = Y(:,1);
    182. 

  182.     else
    183. 

  183.       error ("interp2: X and Y must be matrices of same size");
    184. 

  184.     endif
    185. 

  185.     if (columns (Z) != length (X) || rows (Z) != length (Y))
    186. 

  186.       error ("interp2: X and Y size must match the dimensions of Z");
    187. 

  187.     endif
    188. 

  188. 

  189.     ## If Xi and Yi are vectors of different orientation build a grid
    190. 

  190.     if ((rows (XI) == 1 && columns (YI) == 1)
    191. 

  191.         || (columns (XI) == 1 && rows (YI) == 1))
    192. 

  192.       [XI, YI] = meshgrid (XI, YI);
    193. 

  193.     elseif (! size_equal (XI, YI))
    194. 

  194.       error ("interp2: XI and YI must be matrices of equal size");
    195. 

  195.     endif
    196. 

  196. 

  197.     ## if XI, YI are vectors, X and Y should share their orientation.
    198. 

  198.     if (rows (XI) == 1)
    199. 

  199.       if (rows (X) != 1)
    200. 

  200.         X = X.';
    201. 

  201.       endif
    202. 

  202.       if (rows (Y) != 1)
    203. 

  203.         Y = Y.';
    204. 

  204.       endif
    205. 

  205.     elseif (columns (XI) == 1)
    206. 

  206.       if (columns (X) != 1)
    207. 

  207.         X = X.';
    208. 

  208.       endif
    209. 

  209.       if (columns (Y) != 1)
    210. 

  210.         Y = Y.';
    211. 

  211.       endif
    212. 

  212.     endif
    213. 

  213. 

  214.     xidx = lookup (X, XI, "lr");
    215. 

  215.     yidx = lookup (Y, YI, "lr");
    216. 

  216. 

  217.     if (strcmp (method, "linear"))
    218. 

  218.       ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy
    219. 

  219.       ##
    220. 

  220.       ## a-b
    221. 

  221.       ## | |
    222. 

  222.       ## c-d
    223. 

  223.       a = Z(1:(zr - 1), 1:(zc - 1));
    224. 

  224.       b = Z(1:(zr - 1), 2:zc) - a;
    225. 

  225.       c = Z(2:zr, 1:(zc - 1)) - a;
    226. 

  226.       d = Z(2:zr, 2:zc) - a - b - c;
    227. 

  227. 

  228.       ## scale XI, YI values to a 1-spaced grid
    229. 

  229.       Xsc = (XI - X(xidx)) ./ (diff (X)(xidx));
    230. 

  230.       Ysc = (YI - Y(yidx)) ./ (diff (Y)(yidx));
    231. 

  231. 

  232.       ## Get 2D index.
    233. 

  233.       idx = sub2ind (size (a), yidx, xidx);
    234. 

  234.       ## We can dispose of the 1D indices at this point to save memory.
    235. 

  235.       clear xidx yidx;
    236. 

  236. 

  237.       ## apply plane equation
    238. 

  238.       ZI = a(idx) + b(idx).*Xsc + c(idx).*Ysc + d(idx).*Xsc.*Ysc;
    239. 

  239. 

  240.     elseif (strcmp (method, "nearest"))
    241. 

  241.       ii = (XI - X(xidx) >= X(xidx + 1) - XI);
    242. 

  242.       jj = (YI - Y(yidx) >= Y(yidx + 1) - YI);
    243. 

  243.       idx = sub2ind (size (Z), yidx+jj, xidx+ii);
    244. 

  244.       ZI = Z(idx);
    245. 

  245. 

  246.     elseif (strcmp (method, "pchip"))
    247. 

  247. 

  248.       if (length (X) < 2 || length (Y) < 2)
    249. 

  249.         error ("interp2: pchip2 requires at least 2 points in each dimension");
    250. 

  250.       endif
    251. 

  251. 

  252.       ## first order derivatives
    253. 

  253.       DX = __pchip_deriv__ (X, Z, 2);
    254. 

  254.       DY = __pchip_deriv__ (Y, Z, 1);
    255. 

  255.       ## Compute mixed derivatives row-wise and column-wise, use the average.
    256. 

  256.       DXY = (__pchip_deriv__ (X, DY, 2) + __pchip_deriv__ (Y, DX, 1))/2;
    257. 

  257. 

  258.       ## do the bicubic interpolation
    259. 

  259.       hx = diff (X); hx = hx(xidx);
    260. 

  260.       hy = diff (Y); hy = hy(yidx);
    261. 

  261. 

  262.       tx = (XI - X(xidx)) ./ hx;
    263. 

  263.       ty = (YI - Y(yidx)) ./ hy;
    264. 

  264. 

  265.       ## construct the cubic hermite base functions in x, y
    266. 

  266. 

  267.       ## formulas:
    268. 

  268.       ## b{1,1} =    ( 2*t.^3 - 3*t.^2     + 1);
    269. 

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

  270.       ## b{1,2} =    (-2*t.^3 + 3*t.^2        );
    271. 

  271.       ## b{2,2} = h.*(   t.^3 -   t.^2        );
    272. 

  272. 

  273.       ## optimized equivalents of the above:
    274. 

  274.       t1 = tx.^2;
    275. 

  275.       t2 = tx.*t1 - t1;
    276. 

  276.       xb{2,2} = hx.*t2;
    277. 

  277.       t1 = t2 - t1;
    278. 

  278.       xb{2,1} = hx.*(t1 + tx);
    279. 

  279.       t2 += t1;
    280. 

  280.       xb{1,2} = -t2;
    281. 

  281.       xb{1,1} = t2 + 1;
    282. 

  282. 

  283.       t1 = ty.^2;
    284. 

  284.       t2 = ty.*t1 - t1;
    285. 

  285.       yb{2,2} = hy.*t2;
    286. 

  286.       t1 = t2 - t1;
    287. 

  287.       yb{2,1} = hy.*(t1 + ty);
    288. 

  288.       t2 += t1;
    289. 

  289.       yb{1,2} = -t2;
    290. 

  290.       yb{1,1} = t2 + 1;
    291. 

  291. 

  292.       ZI = zeros (size (XI));
    293. 

  293.       for i = 1:2
    294. 

  294.         for j = 1:2
    295. 

  295.           zidx = sub2ind (size (Z), yidx+(j-1), xidx+(i-1));
    296. 

  296.           ZI += xb{1,i} .* yb{1,j} .*   Z(zidx);
    297. 

  297.           ZI += xb{2,i} .* yb{1,j} .*  DX(zidx);
    298. 

  298.           ZI += xb{1,i} .* yb{2,j} .*  DY(zidx);
    299. 

  299.           ZI += xb{2,i} .* yb{2,j} .* DXY(zidx);
    300. 

  300.         endfor
    301. 

  301.       endfor
    302. 

  302. 

  303.     endif
    304. 

  304. 

  305.     if (! isempty (extrapval))
    306. 

  306.       ## set points outside the table to 'extrapval'
    307. 

  307.       if (X (1) < X (end))
    308. 

  308.         if (Y (1) < Y (end))
    309. 

  309.           ZI (XI < X(1,1) | XI > X(end) | YI < Y(1,1) | YI > Y(end)) = ...
    310. 

  310.                   extrapval;
    311. 

  311.         else
    312. 

  312.           ZI (XI < X(1) | XI > X(end) | YI < Y(end) | YI > Y(1)) = ...
    313. 

  313.                   extrapval;
    314. 

  314.         endif
    315. 

  315.       else
    316. 

  316.         if (Y (1) < Y (end))
    317. 

  317.           ZI (XI < X(end) | XI > X(1) | YI < Y(1) | YI > Y(end)) = ...
    318. 

  318.                   extrapval;
    319. 

  319.         else
    320. 

  320.           ZI (XI < X(1,end) | XI > X(1) | YI < Y(end) | YI > Y(1)) = ...
    321. 

  321.                   extrapval;
    322. 

  322.         endif
    323. 

  323.       endif
    324. 

  324.     endif
    325. 

  325. 

  326.   else
    327. 

  327. 

  328.     ## Check dimensions of X and Y
    329. 

  329.     if (isvector (X) && isvector (Y))
    330. 

  330.       X = X(:).';
    331. 

  331.       Y = Y(:);
    332. 

  332.       if (!isequal ([length(Y), length(X)], size(Z)))
    333. 

  333.         error ("interp2: X and Y size must match the dimensions of Z");
    334. 

  334.       endif
    335. 

  335.     elseif (!size_equal (X, Y))
    336. 

  336.       error ("interp2: X and Y must be matrices of equal size");
    337. 

  337.       if (! size_equal (X, Z))
    338. 

  338.         error ("interp2: X and Y size must match the dimensions of Z");
    339. 

  339.       endif
    340. 

  340.     endif
    341. 

  341. 

  342.     ## Check dimensions of XI and YI
    343. 

  343.     if (isvector (XI) && isvector (YI) && ! size_equal (XI, YI))
    344. 

  344.       XI = XI(:).';
    345. 

  345.       YI = YI(:);
    346. 

  346.       [XI, YI] = meshgrid (XI, YI);
    347. 

  347.     elseif (! size_equal (XI, YI))
    348. 

  348.       error ("interp2: XI and YI must be matrices of equal size");
    349. 

  349.     endif
    350. 

  350. 

  351.     if (strcmp (method, "cubic"))
    352. 

  352.       if (isgriddata (XI) && isgriddata (YI'))
    353. 

  353.         ZI = bicubic (X, Y, Z, XI (1, :), YI (:, 1), extrapval);
    354. 

  354.       elseif (isgriddata (X) && isgriddata (Y'))
    355. 

  355.         ## Allocate output
    356. 

  356.         ZI = zeros (size (X));
    357. 

  357. 

  358.         ## Find inliers
    359. 

  359.         inside = !(XI < X (1) | XI > X (end) | YI < Y (1) | YI > Y (end));
    360. 

  360. 

  361.         ## Scale XI and YI to match indices of Z
    362. 

  362.         XI = (columns (Z) - 1) * (XI - X (1)) / (X (end) - X (1)) + 1;
    363. 

  363.         YI = (rows (Z) - 1) * (YI - Y (1)) / (Y (end) - Y (1)) + 1;
    364. 

  364. 

  365.         ## Start the real work
    366. 

  366.         K = floor (XI);
    367. 

  367.         L = floor (YI);
    368. 

  368. 

  369.         ## Coefficients
    370. 

  370.         AY1  = bc ((YI - L + 1));
    371. 

  371.         AX1  = bc ((XI - K + 1));
    372. 

  372.         AY0  = bc ((YI - L + 0));
    373. 

  373.         AX0  = bc ((XI - K + 0));
    374. 

  374.         AY_1 = bc ((YI - L - 1));
    375. 

  375.         AX_1 = bc ((XI - K - 1));
    376. 

  376.         AY_2 = bc ((YI - L - 2));
    377. 

  377.         AX_2 = bc ((XI - K - 2));
    378. 

  378. 

  379.         ## Perform interpolation
    380. 

  380.         sz = size(Z);
    381. 

  381.         ZI = AY_2 .* AX_2 .* Z (sym_sub2ind (sz, L+2, K+2)) ...
    382. 

  382.            + AY_2 .* AX_1 .* Z (sym_sub2ind (sz, L+2, K+1)) ...
    383. 

  383.            + AY_2 .* AX0  .* Z (sym_sub2ind (sz, L+2, K))   ...
    384. 

  384.            + AY_2 .* AX1  .* Z (sym_sub2ind (sz, L+2, K-1)) ...
    385. 

  385.            + AY_1 .* AX_2 .* Z (sym_sub2ind (sz, L+1, K+2)) ...
    386. 

  386.            + AY_1 .* AX_1 .* Z (sym_sub2ind (sz, L+1, K+1)) ...
    387. 

  387.            + AY_1 .* AX0  .* Z (sym_sub2ind (sz, L+1, K))   ...
    388. 

  388.            + AY_1 .* AX1  .* Z (sym_sub2ind (sz, L+1, K-1)) ...
    389. 

  389.            + AY0  .* AX_2 .* Z (sym_sub2ind (sz, L,   K+2)) ...
    390. 

  390.            + AY0  .* AX_1 .* Z (sym_sub2ind (sz, L,   K+1)) ...
    391. 

  391.            + AY0  .* AX0  .* Z (sym_sub2ind (sz, L,   K))   ...
    392. 

  392.            + AY0  .* AX1  .* Z (sym_sub2ind (sz, L,   K-1)) ...
    393. 

  393.            + AY1  .* AX_2 .* Z (sym_sub2ind (sz, L-1, K+2)) ...
    394. 

  394.            + AY1  .* AX_1 .* Z (sym_sub2ind (sz, L-1, K+1)) ...
    395. 

  395.            + AY1  .* AX0  .* Z (sym_sub2ind (sz, L-1, K))   ...
    396. 

  396.            + AY1  .* AX1  .* Z (sym_sub2ind (sz, L-1, K-1));
    397. 

  397.         ZI (!inside) = extrapval;
    398. 

  398. 

  399.       else
    400. 

  400.         error ("interp2: input data must have `meshgrid' format");
    401. 

  401.       endif
    402. 

  402. 

  403.     elseif (strcmp (method, "spline"))
    404. 

  404.       if (isgriddata (XI) && isgriddata (YI'))
    405. 

  405.         ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrapval,
    406. 

  406.                         "spline");
    407. 

  407.       else
    408. 

  408.         error ("interp2: input data must have `meshgrid' format");
    409. 

  409.       endif
    410. 

  410.     else
    411. 

  411.       error ("interp2: interpolation METHOD not recognized");
    412. 

  412.     endif
    413. 

  413. 

  414.   endif
    415. 

  415. endfunction
    416. 

  416. 

  417. function b = isgriddata (X)
    418. 

  418.   d1 = diff (X, 1, 1);
    419. 

  419.   b = all (d1 (:) == 0);
    420. 

  420. endfunction
    421. 

  421. 

  422. ## Compute the bicubic interpolation coefficients
    423. 

  423. function o = bc(x)
    424. 

  424.   x = abs(x);
    425. 

  425.   o = zeros(size(x));
    426. 

  426.   idx1 = (x < 1);
    427. 

  427.   idx2 = !idx1 & (x < 2);
    428. 

  428.   o(idx1) = 1 - 2.*x(idx1).^2 + x(idx1).^3;
    429. 

  429.   o(idx2) = 4 - 8.*x(idx2) + 5.*x(idx2).^2 - x(idx2).^3;
    430. 

  430. endfunction
    431. 

  431. 

  432. ## This version of sub2ind behaves as if the data was symmetrically padded
    433. 

  433. function ind = sym_sub2ind(sz, Y, X)
    434. 

  434.   Y (Y < 1) = 1 - Y (Y < 1);
    435. 

  435.   while (any (Y (:) > 2 * sz (1)))
    436. 

  436.     Y (Y > 2 * sz (1)) = round (Y (Y > 2 * sz (1)) / 2);
    437. 

  437.   endwhile
    438. 

  438.   Y (Y > sz (1)) = 1 + 2 * sz (1) - Y (Y > sz (1));
    439. 

  439.   X (X < 1) = 1 - X (X < 1);
    440. 

  440.   while (any (X (:) > 2 * sz (2)))
    441. 

  441.     X (X > 2 * sz (2)) = round (X (X > 2 * sz (2)) / 2);
    442. 

  442.   endwhile
    443. 

  443.   X (X > sz (2)) = 1 + 2 * sz (2) - X (X > sz (2));
    444. 

  444.   ind = sub2ind(sz, Y, X);
    445. 

  445. endfunction
    446. 

  446. 

  447. 

  448. %!demo
    449. 

  449. %! A=[13,-1,12;5,4,3;1,6,2];
    450. 

  450. %! x=[0,1,4]; y=[10,11,12];
    451. 

  451. %! xi=linspace(min(x),max(x),17);
    452. 

  452. %! yi=linspace(min(y),max(y),26)';
    453. 

  453. %! mesh(xi,yi,interp2(x,y,A,xi,yi,'linear'));
    454. 

  454. %! [x,y] = meshgrid(x,y);
    455. 

  455. %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
    456. 

  456. 

  457. %!demo
    458. 

  458. %! [x,y,A] = peaks(10);
    459. 

  459. %! x = x(1,:)'; y = y(:,1);
    460. 

  460. %! xi=linspace(min(x),max(x),41);
    461. 

  461. %! yi=linspace(min(y),max(y),41)';
    462. 

  462. %! mesh(xi,yi,interp2(x,y,A,xi,yi,'linear'));
    463. 

  463. %! [x,y] = meshgrid(x,y);
    464. 

  464. %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
    465. 

  465. 

  466. %!demo
    467. 

  467. %! A=[13,-1,12;5,4,3;1,6,2];
    468. 

  468. %! x=[0,1,4]; y=[10,11,12];
    469. 

  469. %! xi=linspace(min(x),max(x),17);
    470. 

  470. %! yi=linspace(min(y),max(y),26)';
    471. 

  471. %! mesh(xi,yi,interp2(x,y,A,xi,yi,'nearest'));
    472. 

  472. %! [x,y] = meshgrid(x,y);
    473. 

  473. %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
    474. 

  474. 

  475. %!demo
    476. 

  476. %! [x,y,A] = peaks(10);
    477. 

  477. %! x = x(1,:)'; y = y(:,1);
    478. 

  478. %! xi=linspace(min(x),max(x),41);
    479. 

  479. %! yi=linspace(min(y),max(y),41)';
    480. 

  480. %! mesh(xi,yi,interp2(x,y,A,xi,yi,'nearest'));
    481. 

  481. %! [x,y] = meshgrid(x,y);
    482. 

  482. %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
    483. 

  483. 

  484. %!demo
    485. 

  485. %! A=[13,-1,12;5,4,3;1,6,2];
    486. 

  486. %! x=[0,1,2]; y=[10,11,12];
    487. 

  487. %! xi=linspace(min(x),max(x),17);
    488. 

  488. %! yi=linspace(min(y),max(y),26)';
    489. 

  489. %! mesh(xi,yi,interp2(x,y,A,xi,yi,'pchip'));
    490. 

  490. %! [x,y] = meshgrid(x,y);
    491. 

  491. %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
    492. 

  492. 

  493. %!demo
    494. 

  494. %! [x,y,A] = peaks(10);
    495. 

  495. %! x = x(1,:)'; y = y(:,1);
    496. 

  496. %! xi=linspace(min(x),max(x),41);
    497. 

  497. %! yi=linspace(min(y),max(y),41)';
    498. 

  498. %! mesh(xi,yi,interp2(x,y,A,xi,yi,'pchip'));
    499. 

  499. %! [x,y] = meshgrid(x,y);
    500. 

  500. %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
    501. 

  501. 

  502. %!demo
    503. 

  503. %! A=[13,-1,12;5,4,3;1,6,2];
    504. 

  504. %! x=[0,1,2]; y=[10,11,12];
    505. 

  505. %! xi=linspace(min(x),max(x),17);
    506. 

  506. %! yi=linspace(min(y),max(y),26)';
    507. 

  507. %! mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic'));
    508. 

  508. %! [x,y] = meshgrid(x,y);
    509. 

  509. %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
    510. 

  510. 

  511. %!demo
    512. 

  512. %! [x,y,A] = peaks(10);
    513. 

  513. %! x = x(1,:)'; y = y(:,1);
    514. 

  514. %! xi=linspace(min(x),max(x),41);
    515. 

  515. %! yi=linspace(min(y),max(y),41)';
    516. 

  516. %! mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic'));
    517. 

  517. %! [x,y] = meshgrid(x,y);
    518. 

  518. %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
    519. 

  519. 

  520. %!demo
    521. 

  521. %! A=[13,-1,12;5,4,3;1,6,2];
    522. 

  522. %! x=[0,1,2]; y=[10,11,12];
    523. 

  523. %! xi=linspace(min(x),max(x),17);
    524. 

  524. %! yi=linspace(min(y),max(y),26)';
    525. 

  525. %! mesh(xi,yi,interp2(x,y,A,xi,yi,'spline'));
    526. 

  526. %! [x,y] = meshgrid(x,y);
    527. 

  527. %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
    528. 

  528. 

  529. %!demo
    530. 

  530. %! [x,y,A] = peaks(10);
    531. 

  531. %! x = x(1,:)'; y = y(:,1);
    532. 

  532. %! xi=linspace(min(x),max(x),41);
    533. 

  533. %! yi=linspace(min(y),max(y),41)';
    534. 

  534. %! mesh(xi,yi,interp2(x,y,A,xi,yi,'spline'));
    535. 

  535. %! [x,y] = meshgrid(x,y);
    536. 

  536. %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
    537. 

  537. 

  538. %!test % simple test
    539. 

  539. %!  x = [1,2,3];
    540. 

  540. %!  y = [4,5,6,7];
    541. 

  541. %!  [X, Y] = meshgrid(x,y);
    542. 

  542. %!  Orig = X.^2 + Y.^3;
    543. 

  543. %!  xi = [1.2,2, 1.5];
    544. 

  544. %!  yi = [6.2, 4.0, 5.0]';
    545. 

  545. %!
    546. 

  546. %!  Expected = ...
    547. 

  547. %!    [243,   245.4,  243.9;
    548. 

  548. %!      65.6,  68,     66.5;
    549. 

  549. %!     126.6, 129,    127.5];
    550. 

  550. %!  Result = interp2(x,y,Orig, xi, yi);
    551. 

  551. %!
    552. 

  552. %!  assert(Result, Expected, 1000*eps);
    553. 

  553. 

  554. %!test % 2^n form
    555. 

  555. %!  x = [1,2,3];
    556. 

  556. %!  y = [4,5,6,7];
    557. 

  557. %!  [X, Y] = meshgrid(x,y);
    558. 

  558. %!  Orig = X.^2 + Y.^3;
    559. 

  559. %!  xi = [1:0.25:3]; yi = [4:0.25:7]';
    560. 

  560. %!  Expected = interp2(x,y,Orig, xi, yi);
    561. 

  561. %!  Result = interp2(Orig,2);
    562. 

  562. %!
    563. 

  563. %!  assert(Result, Expected, 10*eps);
    564. 

  564. 

  565. %!test % matrix slice
    566. 

  566. %!  A = eye(4);
    567. 

  567. %!  assert(interp2(A,[1:4],[1:4]),[1,1,1,1]);
    568. 

  568. 

  569. %!test % non-gridded XI,YI
    570. 

  570. %!  A = eye(4);
    571. 

  571. %!  assert(interp2(A,[1,2;3,4],[1,3;2,4]),[1,0;0,1]);
    572. 

  572. 

  573. %!test % for values outside of boundaries
    574. 

  574. %!  x = [1,2,3];
    575. 

  575. %!  y = [4,5,6,7];
    576. 

  576. %!  [X, Y] = meshgrid(x,y);
    577. 

  577. %!  Orig = X.^2 + Y.^3;
    578. 

  578. %!  xi = [0,4];
    579. 

  579. %!  yi = [3,8]';
    580. 

  580. %!  assert(interp2(x,y,Orig, xi, yi),[NA,NA;NA,NA]);
    581. 

  581. %!  assert(interp2(x,y,Orig, xi, yi,'linear', 0),[0,0;0,0]);
    582. 

  582. 

  583. %!test % for values at boundaries
    584. 

  584. %!  A=[1,2;3,4];
    585. 

  585. %!  x=[0,1];
    586. 

  586. %!  y=[2,3]';
    587. 

  587. %!  assert(interp2(x,y,A,x,y,'linear'), A);
    588. 

  588. %!  assert(interp2(x,y,A,x,y,'nearest'), A);
    589. 

  589. 

  590. %!test % for Matlab-compatible rounding for 'nearest'
    591. 

  591. %! X = meshgrid (1:4);
    592. 

  592. %! assert (interp2 (X, 2.5, 2.5, 'nearest'), 3);
    593. 

  593. 

  594. %!shared z, zout, tol
    595. 

  595. %!  z = [1 3 5; 3 5 7; 5 7 9];
    596. 

  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. 

  597. %!  tol = 2 * eps;
    598. 

  598. %!assert (interp2 (z), zout, tol);
    599. 

  599. %!assert (interp2 (z, "linear"), zout, tol);
    600. 

  600. %!assert (interp2 (z, "pchip"), zout, tol);
    601. 

  601. %!assert (interp2 (z, "cubic"), zout, 10 * tol);
    602. 

  602. %!assert (interp2 (z, "spline"), zout, tol);
    603. 

  603. %!assert (interp2 (z, [2 3 1], [2 2 2]', "linear"), repmat ([5, 7, 3], [3, 1]), tol) 
    604. 

  604. %!assert (interp2 (z, [2 3 1], [2 2 2]', "pchip"), repmat ([5, 7, 3], [3, 1]), tol) 
    605. 

  605. %!assert (interp2 (z, [2 3 1], [2 2 2]', "cubic"), repmat ([5, 7, 3], [3, 1]), 10 * tol) 
    606. 

  606. %!assert (interp2 (z, [2 3 1], [2 2 2]', "spline"), repmat ([5, 7, 3], [3, 1]), tol) 
    607. 

  607. %!assert (interp2 (z, [2 3 1], [2 2 2], "linear"), [5 7 3], tol);
    608. 

  608. %!assert (interp2 (z, [2 3 1], [2 2 2], "pchip"), [5 7 3], tol);
    609. 

  609. %!assert (interp2 (z, [2 3 1], [2 2 2], "cubic"), [5 7 3], 10 * tol);
    610. 

  610. %!assert (interp2 (z, [2 3 1], [2 2 2], "spline"), [5 7 3], tol);
    611.