1 ## Copyright (C) 2000-2012 Kai Habel
2 ## Copyright (C) 2009 Jaroslav Hajek
4 ## This file is part of Octave.
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.
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.
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/>.
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})
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}
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}.
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.
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})}
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.
51 ## The variable @var{method} defines the method to use for the
52 ## interpolation. It can take one of the following values
56 ## Return the nearest neighbor.
59 ## Linear interpolation from nearest neighbors.
62 ## Piecewise cubic Hermite interpolating polynomial.
65 ## Cubic interpolation from four nearest neighbors.
68 ## Cubic spline interpolation---smooth first and second derivatives
69 ## throughout the curve.
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.
80 ## Author: Kai Habel <kai.habel@gmx.de>
81 ## 2005-03-02 Thomas Weber <weber@num.uni-sb.de>
83 ## 2005-03-02 Paul Kienzle <pkienzle@users.sf.net>
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
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.
96 function ZI = interp2 (varargin)
97 Z = X = Y = XI = YI = n = [];
106 if (ischar (varargin{2}))
107 [Z, method] = deal (varargin{:});
110 [Z, n] = deal (varargin{:});
113 if (ischar (varargin{3}))
114 [Z, n, method] = deal (varargin{:});
116 [Z, XI, YI] = deal (varargin{:});
119 if (ischar (varargin{4}))
120 [Z, XI, YI, method] = deal (varargin{:});
122 [Z, n, method, extrapval] = deal (varargin{:});
125 if (ischar (varargin{4}))
126 [Z, XI, YI, method, extrapval] = deal (varargin{:});
128 [X, Y, Z, XI, YI] = deal (varargin{:});
131 [X, Y, Z, XI, YI, method] = deal (varargin{:});
133 [X, Y, Z, XI, YI, method, extrapval] = deal (varargin{:});
140 error ("interp2: Z must be a matrix");
142 if (!isempty (n) && !isscalar (n))
143 error ("interp2: N must be a scalar");
145 if (!ischar (method))
146 error ("interp2: METHOD must be a string");
148 if (ischar (extrapval) || strcmp (extrapval, "extrap"))
150 elseif (!isscalar (extrapval))
151 error ("interp2: EXTRAPVAL must be a scalar");
154 ## Define X, Y, XI, YI if needed
160 if (! isnumeric (X) || ! isnumeric (Y))
161 error ("interp2: X, Y must be numeric matrices");
164 ## Calculate the interleaved input vectors.
169 if (! isnumeric (XI) || ! isnumeric (YI))
170 error ("interp2: XI, YI must be numeric");
174 if (strcmp (method, "linear") || strcmp (method, "nearest") ...
175 || strcmp (method, "pchip"))
177 ## If X and Y vectors produce a grid from them
178 if (isvector (X) && isvector (Y))
180 elseif (size_equal (X, Y))
181 X = X(1,:)'; Y = Y(:,1);
183 error ("interp2: X and Y must be matrices of same size");
185 if (columns (Z) != length (X) || rows (Z) != length (Y))
186 error ("interp2: X and Y size must match the dimensions of Z");
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");
197 ## if XI, YI are vectors, X and Y should share their orientation.
205 elseif (columns (XI) == 1)
206 if (columns (X) != 1)
209 if (columns (Y) != 1)
214 xidx = lookup (X, XI, "lr");
215 yidx = lookup (Y, YI, "lr");
217 if (strcmp (method, "linear"))
218 ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy
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;
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));
233 idx = sub2ind (size (a), yidx, xidx);
234 ## We can dispose of the 1D indices at this point to save memory.
237 ## apply plane equation
238 ZI = a(idx) + b(idx).*Xsc + c(idx).*Ysc + d(idx).*Xsc.*Ysc;
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);
246 elseif (strcmp (method, "pchip"))
248 if (length (X) < 2 || length (Y) < 2)
249 error ("interp2: pchip2 requires at least 2 points in each dimension");
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;
258 ## do the bicubic interpolation
259 hx = diff (X); hx = hx(xidx);
260 hy = diff (Y); hy = hy(yidx);
262 tx = (XI - X(xidx)) ./ hx;
263 ty = (YI - Y(yidx)) ./ hy;
265 ## construct the cubic hermite base functions in x, y
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 );
273 ## optimized equivalents of the above:
278 xb{2,1} = hx.*(t1 + tx);
287 yb{2,1} = hy.*(t1 + ty);
292 ZI = zeros (size (XI));
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);
305 if (! isempty (extrapval))
306 ## set points outside the table to 'extrapval'
309 ZI (XI < X(1,1) | XI > X(end) | YI < Y(1,1) | YI > Y(end)) = ...
312 ZI (XI < X(1) | XI > X(end) | YI < Y(end) | YI > Y(1)) = ...
317 ZI (XI < X(end) | XI > X(1) | YI < Y(1) | YI > Y(end)) = ...
320 ZI (XI < X(1,end) | XI > X(1) | YI < Y(end) | YI > Y(1)) = ...
328 ## Check dimensions of X and Y
329 if (isvector (X) && isvector (Y))
332 if (!isequal ([length(Y), length(X)], size(Z)))
333 error ("interp2: X and Y size must match the dimensions of Z");
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");
342 ## Check dimensions of XI and YI
343 if (isvector (XI) && isvector (YI) && ! size_equal (XI, 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");
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'))
356 ZI = zeros (size (X));
359 inside = !(XI < X (1) | XI > X (end) | YI < Y (1) | YI > Y (end));
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;
365 ## Start the real work
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));
379 ## Perform interpolation
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;
400 error ("interp2: input data must have `meshgrid' format");
403 elseif (strcmp (method, "spline"))
404 if (isgriddata (XI) && isgriddata (YI'))
405 ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrapval,
408 error ("interp2: input data must have `meshgrid' format");
411 error ("interp2: interpolation METHOD not recognized");
417 function b = isgriddata (X)
419 b = all (d1 (:) == 0);
422 ## Compute the bicubic interpolation coefficients
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;
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);
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);
443 X (X > sz (2)) = 1 + 2 * sz (2) - X (X > sz (2));
444 ind = sub2ind(sz, Y, X);
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;
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;
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;
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;
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;
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;
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;
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;
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;
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;
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]';
547 %! [243, 245.4, 243.9;
549 %! 126.6, 129, 127.5];
550 %! Result = interp2(x,y,Orig, xi, yi);
552 %! assert(Result, Expected, 1000*eps);
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);
563 %! assert(Result, Expected, 10*eps);
565 %!test % matrix slice
567 %! assert(interp2(A,[1:4],[1:4]),[1,1,1,1]);
569 %!test % non-gridded XI,YI
571 %! assert(interp2(A,[1,2;3,4],[1,3;2,4]),[1,0;0,1]);
573 %!test % for values outside of boundaries
576 %! [X, Y] = meshgrid(x,y);
577 %! Orig = X.^2 + Y.^3;
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]);
583 %!test % for values at boundaries
587 %! assert(interp2(x,y,A,x,y,'linear'), A);
588 %! assert(interp2(x,y,A,x,y,'nearest'), A);
590 %!test % for Matlab-compatible rounding for 'nearest'
591 %! X = meshgrid (1:4);
592 %! assert (interp2 (X, 2.5, 2.5, 'nearest'), 3);
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];
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);