######################################################################## ## ## Copyright (C) 2000-2022 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or . ## ## 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 ## . ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {@var{zi} =} interp2 (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}) ## @deftypefnx {} {@var{zi} =} interp2 (@var{z}, @var{xi}, @var{yi}) ## @deftypefnx {} {@var{zi} =} interp2 (@var{z}, @var{n}) ## @deftypefnx {} {@var{zi} =} interp2 (@var{z}) ## @deftypefnx {} {@var{zi} =} interp2 (@dots{}, @var{method}) ## @deftypefnx {} {@var{zi} =} interp2 (@dots{}, @var{method}, @var{extrap}) ## ## Two-dimensional interpolation. ## ## Interpolate reference data @var{x}, @var{y}, @var{z} to determine @var{zi} ## at the coordinates @var{xi}, @var{yi}. The reference data @var{x}, @var{y} ## can be matrices, as returned by @code{meshgrid}, in which case the sizes of ## @var{x}, @var{y}, and @var{z} must be equal. If @var{x}, @var{y} are ## vectors describing a grid then @code{length (@var{x}) == columns (@var{z})} ## and @code{length (@var{y}) == rows (@var{z})}. In either case the input ## data must be strictly monotonic. ## ## If called without @var{x}, @var{y}, and just a single reference data matrix ## @var{z}, the 2-D region ## @code{@var{x} = 1:columns (@var{z}), @var{y} = 1:rows (@var{z})} is assumed. ## This saves memory if the grid is regular and the distance between points is ## not important. ## ## If called with a single reference data matrix @var{z} and a refinement ## value @var{n}, then perform interpolation over a grid where each original ## interval has been recursively subdivided @var{n} times. This results in ## @code{2^@var{n}-1} additional points for every interval in the original ## grid. If @var{n} is omitted a value of 1 is used. As an example, the ## interval [0,1] with @code{@var{n}==2} results in a refined interval with ## points at [0, 1/4, 1/2, 3/4, 1]. ## ## The interpolation @var{method} is one of: ## ## @table @asis ## @item @qcode{"nearest"} ## Return the nearest neighbor. ## ## @item @qcode{"linear"} (default) ## Linear interpolation from nearest neighbors. ## ## @item @qcode{"pchip"} ## Piecewise cubic Hermite interpolating polynomial---shape-preserving ## interpolation with smooth first derivative. ## ## @item @qcode{"cubic"} ## Cubic interpolation (same as @qcode{"pchip"}). ## ## @item @qcode{"spline"} ## Cubic spline interpolation---smooth first and second derivatives ## throughout the curve. ## @end table ## ## @var{extrap} is a scalar number. It replaces values beyond the endpoints ## with @var{extrap}. Note that if @var{extrap} is used, @var{method} must ## be specified as well. If @var{extrap} is omitted and the @var{method} is ## @qcode{"spline"}, then the extrapolated values of the @qcode{"spline"} are ## used. Otherwise the default @var{extrap} value for any other @var{method} ## is @qcode{"NA"}. ## @seealso{interp1, interp3, interpn, meshgrid} ## @end deftypefn function ZI = interp2 (varargin) narginchk (1, 7); nargs = nargin; Z = X = Y = XI = YI = n = []; method = "linear"; extrap = []; ## Check for method and extrap if (nargs > 1 && ischar (varargin{end-1})) if (! isnumeric (varargin{end}) || ! isscalar (varargin{end})) error ("interp2: EXTRAP must be a numeric scalar"); endif extrap = varargin{end}; method = varargin{end-1}; nargs -= 2; elseif (ischar (varargin{end})) method = varargin{end}; nargs -= 1; endif if (method(1) == "*") warning ("interp2: ignoring unsupported '*' flag to METHOD"); method(1) = []; endif method = validatestring (method, ... {"nearest", "linear", "pchip", "cubic", "spline"}); ## Read numeric input switch (nargs) case 1 Z = varargin{1}; n = 1; case 2 [Z, n] = deal (varargin{1:nargs}); case 3 [Z, XI, YI] = deal (varargin{1:nargs}); case 5 [X, Y, Z, XI, YI] = deal (varargin{1:nargs}); otherwise print_usage (); endswitch ## Type checking if (! isnumeric (Z) || isscalar (Z) || ! ismatrix (Z)) error ("interp2: Z must be a 2-D matrix"); endif if (! isempty (n) && ! (isscalar (n) && n >= 0 && n == fix (n))) error ("interp2: N must be an integer >= 0"); 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 (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 equal size"); endif if (columns (Z) != length (X) || rows (Z) != length (Y)) error ("interp2: X and Y size must match the dimensions of Z"); endif dx = diff (X); if (all (dx < 0)) X = flipud (X); Z = fliplr (Z); elseif (any (dx <= 0)) error ("interp2: X must be strictly monotonic"); endif dy = diff (Y); if (all (dy < 0)) Y = flipud (Y); Z = flipud (Z); elseif (any (dy <= 0)) error ("interp2: Y must be strictly monotonic"); endif if (any (strcmp (method, {"nearest", "linear", "pchip"}))) ## If Xi and Yi are vectors of different orientation build a grid if ((isrow (XI) && iscolumn (YI)) || (iscolumn (XI) && isrow (YI))) [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 (isrow (XI)) if (rows (X) != 1) X = X.'; endif if (rows (Y) != 1) Y = Y.'; endif elseif (iscolumn (XI)) 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); ## Dispose of the 1D indices at this point to save memory. clear xidx yidx; ## Apply plane equation ## Handle case where idx and coefficients are both vectors and resulting ## coeff(idx) follows orientation of coeff, rather than that of idx. forient = @(x) reshape (x, size (idx)); ZI = forient (a(idx)) ... + forient (b(idx)) .* Xsc ... + forient (c(idx)) .* Ysc ... + forient (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: %s requires at least 2 points in each dimension", method); 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 else ## Check dimensions of XI and YI if (isvector (XI) && isvector (YI) && ! size_equal (XI, YI)) XI = XI(:).'; YI = YI(:); elseif (! size_equal (XI, YI)) error ("interp2: XI and YI must be matrices of equal size"); endif if (strcmp (method, "spline")) if (isgriddata (XI) && isgriddata (YI')) ZI = __splinen__ ({Y, X}, Z, {YI(:,1), XI(1,:)}, extrap, "spline"); else error ("interp2: XI, YI must have uniform spacing ('meshgrid' format)"); endif return; # spline doesn't need NA extrapolation value (MATLAB compatibility) elseif (strcmp (method, "cubic")) ## check if interpolation points are a meshgrid and reduce to vectors if possible if (size_equal (XI, YI) && all (all (XI(1, :) == XI & YI(:, 1) == YI))) XI = XI(1, :); YI = YI(:, 1); endif ## make X a row vector X = X.'; ## quadratic padding + additional zeros for the special case of copying ## the last line (like x=1:5, xi=5, requires to have indexes 6 and 7) row_1 = 3*Z(1, :, :) - 3*Z(2, :, :) + Z(3, :, :); row_end = 3*Z(end, :, :) - 3*Z(end-1, :, :) + Z(end-2, :, :); ZI = [3*row_1(:, 1, :) - 3*row_1(:, 2, :) + row_1(:, 3, :), ... row_1, ... 3*row_1(:, end, :) - 3*row_1(:, end-1, :) + row_1(:, end-2, :), ... 0; # 3*Z(:, 1, :) - 3*Z(:, 2, :) + Z(:, 3, :), ... Z, ... 3*Z(:, end, :) - 3*Z(:, end-1, :) + Z(:, end-2, :), ... zeros(rows (Z), 1, size (Z, 3)); # 3*row_end(:, 1, :) - 3*row_end(:, 2, :) + row_end(:, 3, :), ... row_end, ... 3*row_end(:, end, :) - 3*row_end(:, end-1, :) + row_end(:, end-2, :), ... 0; zeros(1, columns (Z) + 3, size (Z, 3))]; ## interpolate if (isrow (XI) && iscolumn (YI)) ZI = conv_interp_vec (ZI, Y, YI, @cubic, [-2, 2], 1); ZI = conv_interp_vec (ZI, X, XI, @cubic, [-2, 2], 2); else ZI = conv_interp_pairs (ZI, X, Y, XI, YI, @cubic, [-2, 2]); endif endif endif ## extrapolation 'extrap' if (isempty (extrap)) if (iscomplex (Z)) extrap = NA + i*NA; else extrap = NA; endif endif 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)) = extrap; else ZI(XI < X(1) | XI > X(end) | YI < Y(end) | YI > Y(1)) = extrap; endif else if (Y(1) < Y(end)) ZI(XI < X(end) | XI > X(1) | YI < Y(1) | YI > Y(end)) = extrap; else ZI(XI < X(1,end) | XI > X(1) | YI < Y(end) | YI > Y(1)) = extrap; endif endif endfunction function b = isgriddata (X) d1 = diff (X, 1, 1); b = ! any (d1(:) != 0); endfunction ## cubic convolution kernel with a = -0.5 for MATLAB compatibility. function w = cubic (h) absh = abs (h); absh01 = absh <= 1; absh12 = absh <= 2 & ~absh01; w = (1.5 * absh.^3 - 2.5 * absh.^2 + 1) .* absh01 + ... ## for |h| <= 1 (-0.5 * absh.^3 + 2.5 * absh.^2 - 4 * absh + 2) .* absh12; ## for 1 < |h| <= 2 end ## bicubic interpolation of full matrix in one direction with vector function out = conv_interp_vec (Z, XY, XIYI, kernel, kernel_bounds, axis) ## allocate output out_shape = size (Z); out_shape(axis) = length (XIYI); out = zeros (out_shape); ## get indexes and distances h spread = abs (XY(1) - XY(2)); idx = lookup (XY, XIYI, "l"); h = (XIYI - XY(idx)) / spread; idx += -kernel_bounds(1) - 1; ## apply padding for indexes ## interpolate for shift = kernel_bounds(1)+1:kernel_bounds(2) if axis == 1 out += Z(idx + shift, :, :) .* kernel (shift - h); else out += Z(:, idx + shift, :) .* kernel (shift - h); endif endfor endfunction ## bicubic interpolation of arbitrary XI-YI-pairs function out = conv_interp_pairs (Z, X, Y, XI, YI, kernel, kernel_bounds) spread_x = abs (X(1, 1) - X(1, 2)); spread_y = abs (Y(1, 1) - Y(2, 1)); idx_x = lookup (X, XI, "l"); idx_y = lookup (Y, YI, "l"); h_x = (XI - reshape (X(idx_x), size (idx_x))) / spread_x; h_y = (YI - reshape (Y(idx_y), size (idx_y))) / spread_y; idx_x += -kernel_bounds(1) - 1; ## apply padding for indexes idx_y += -kernel_bounds(1) - 1; ## apply padding for indexes out = zeros (rows (YI), columns (XI), size (Z, 3)); for shift_x = kernel_bounds(1)+1:kernel_bounds(2) ## interpolate in y-direction out_y = zeros (rows (YI), columns (XI), size (Z, 3)); for shift_y = kernel_bounds(1)+1:kernel_bounds(2) pixels = Z(sub2ind (size (Z), idx_y + shift_y, idx_x + shift_x)); out_y += pixels .* kernel(shift_y - h_y); endfor ## interpolate in x-direction out += out_y .* kernel (shift_x - h_x); endfor endfunction %!demo %! clf; %! colormap ("default"); %! 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 %! clf; %! colormap ("default"); %! [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 %! clf; %! colormap ("default"); %! 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 %! clf; %! colormap ("default"); %! [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 %! clf; %! colormap ("default"); %! 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 %! clf; %! colormap ("default"); %! [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 %! clf; %! colormap ("default"); %! 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 %! clf; %! colormap ("default"); %! [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 %! clf; %! colormap ("default"); %! 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 %! clf; %! colormap ("default"); %! [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; %!shared x, y, orig, xi, yi, expected %!test # simple test %! x = [1,2,3]; %! y = [4,5,6,7]; %! [X, Y] = meshgrid (x, y); %! orig = X.^2 + Y.^3; %! xi = [0.8, 1.2, 2.0, 1.5]; %! yi = [6.2, 4.0, 5.0, 7.1]'; %! %! # check nearest neighbor %! expected = ... %! [NA, 217, 220, 220; %! NA, 65, 68, 68; %! NA, 126, 129, 129; %! NA, NA, NA, NA]; %! result = interp2 (x, y, orig, xi, yi, "nearest"); %! assert (result, expected); %! %! # check invariance of translation %! result = interp2 (x+3, y-7, orig, xi+3, yi-7, "nearest"); %! assert (result, expected); %! %! # check invariance of scaling %! result = interp2 (x*3, y*(-7), orig, xi*3, yi*(-7), "nearest"); %! assert (result, expected); %! %! # check interpolation with index pairs %! result = interp2 (x, y, orig, xi(2:4), yi(1:3)', "nearest"); %! assert (result, expected(sub2ind(size(expected), 1:3, 2:4))); %! %! # check bilinear interpolation %! expected = ... %! [NA, 243, 245.4, 243.9; %! NA, 65.6, 68, 66.5; %! NA, 126.6, 129, 127.5; %! NA, NA, NA, NA]; %! result = interp2 (x, y, orig, xi, yi); %! assert (result, expected, 1000*eps); %! %! # check invariance of translation %! result = interp2 (x+3, y-7, orig, xi+3, yi-7); %! assert (result, expected, 1000*eps); %! %! # check invariance of scaling %! result = interp2 (x*3, y*(-7), orig, xi*3, yi*(-7)); %! assert (result, expected, 1000*eps); %! %! # check interpolation with index pairs %! result = interp2 (x, y, orig, xi(2:4), yi(1:3)'); %! assert (result, expected(sub2ind(size(expected), 1:3, 2:4)), 1000*eps); %! %! # check spline interpolation %! expected = ... %! [238.968, 239.768, 242.328, 240.578; %! 64.64, 65.44, 68, 66.25; %! 125.64, 126.44, 129, 127.25; %! 358.551, 359.351, 361.911, 360.161]; %! result = interp2 (x, y, orig, xi, yi, "spline"); %! assert (result, expected, 1000*eps); %! %! # check invariance of translation %! result = interp2 (x+3, y-7, orig, xi+3, yi-7, "spline"); %! assert (result, expected, 1000*eps); %! %! # check invariance of scaling %! result = interp2 (x*3, y*(-7), orig, xi*3, yi*(-7), "spline"); %! assert (result, expected, 1000*eps); %! %!xtest %! # FIXME: spline interpolation does not support index pairs, Matlab does. %! result = interp2 (x, y, orig, xi(2:4), yi(1:3)', "spline"); %! assert (result, expected(sub2ind(size(expected), 1:3, 2:4)), 1000*eps); %! %!test %! # check bicubic interpolation %! expected = ... %! [NA, 239.96, 242.52, 240.77; %! NA, 65.44, 68, 66.25; %! NA, 126.44, 129, 127.25; %! NA, NA, NA, NA]; %! result = interp2 (x, y, orig, xi, yi, "cubic"); %! assert (result, expected, 10000*eps); %! %! # check invariance of translation %! result = interp2 (x+3, y-7, orig, xi+3, yi-7, "cubic"); %! assert (result, expected, 10000*eps); %! %! # check invariance of scaling %! result = interp2 (x*3, y*(-7), orig, xi*3, yi*(-7), "cubic"); %! assert (result, expected, 10000*eps); %! %! # check interpolation with index pairs %! result = interp2 (x, y, orig, xi(2:4), yi(1:3)', "cubic"); %! assert (result, expected(sub2ind(size(expected), 1:3, 2:4)), 10000*eps); ## Test that interpolating a complex matrix is equivalent to interpolating its ## real and imaginary parts separately. %!test %! %! xi = 3:0.25:5; %! yi = (1:0.25:2)'; %! orig = rand (5, 4) + i * rand (5, 4); %! for method = {"nearest", "linear", "cubic", "spline"} %! interp_complex = interp2 (orig, xi, yi, method{1}); %! interp_real = interp2 (real (orig), xi, yi, method{1}); %! interp_imag = interp2 (imag (orig), xi, yi, method{1}); %! assert (real (interp_complex), interp_real) %! assert (imag (interp_complex), interp_imag) %! endfor %!test # 2^n refinement 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]); %! assert (interp2 (x,y,orig, xi, yi,"linear", 2), [2,2;2,2]); %! assert (interp2 (x,y,orig, xi, yi,"spline", 2), [2,2;2,2]); %! assert (interp2 (x,y,orig, xi, yi,"linear", 0+1i), [0+1i,0+1i;0+1i,0+1i]); %! assert (interp2 (x,y,orig, xi, yi,"spline"), [27,43;512,528]); %! assert (interp2 (x,y,orig, xi, yi,"cubic"), [NA,NA;NA,NA]); %! assert (interp2 (x,y,orig, xi, yi,"cubic", 2), [2,2;2,2]); %!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); ## re-order monotonically decreasing %!assert <*41838> (interp2 ([1 2 3], [3 2 1], magic (3), 2.5, 3), 3.5) %!assert <*41838> (interp2 ([3 2 1], [1 2 3], magic (3), 1.5, 1), 3.5) ## Linear interpretation with vector XI doesn't lead to matrix output %!assert <*49506> (interp2 ([2 3], [2 3 4], [1 2; 3 4; 5 6], [2 3], 3, "linear"), [3 4]) %!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, 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]), 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], tol) %!assert (interp2 (z, [2 3 1], [2 2 2], "spline"), [5 7 3], tol) %!assert (interp2 (z, [3; 3; 3], [2; 3; 1], "linear"), [7; 9; 5], tol) %!assert (interp2 (z, [3; 3; 3], [2; 3; 1], "pchip"), [7; 9; 5], tol) %!assert (interp2 (z, [3; 3; 3], [2; 3; 1], "cubic"), [7; 9; 5], tol) %!xtest %! # FIXME: single column yields single row with spline interpolation (numbers are correct) %! assert (interp2 (z, [3; 3; 3], [2; 3; 1], "spline"), [7; 9; 5], tol) ## Test input validation %!error interp2 (1, 1, 1, 1, 1, 2) # only 5 numeric inputs %!error interp2 (1, 1, 1, 1, 1, 2, 2) # only 5 numeric inputs %!error interp2 ({1}) %!error interp2 (1,1,1) %!error interp2 (ones (2,2,2)) %!error interp2 (ones (2), ones (2)) %!error interp2 (ones (2), -1) %!error interp2 (ones (2), 1.5) %!warning interp2 (rand (3,3), 1, "*linear"); %!error interp2 (1, 1, 1, 1, 1, "linear", {1}) %!error interp2 (1, 1, 1, 1, 1, "linear", ones (2,2)) %!error interp2 (1, 1, 1, 1, 1, "linear", "abc") %!error interp2 (1, 1, 1, 1, 1, "linear", "extrap") %!error interp2 ({1}, 1, ones (2), 1, 1) %!error interp2 (1, {1}, ones (2), 1, 1) %!error interp2 (1, 1, ones (2), {1}, 1) %!error interp2 (1, 1, ones (2), 1, {1}) %!error interp2 (ones (2,2), 1, ones (2), 1, 1) %!error interp2 (ones (2,2), ones (2,3), ones (2), 1, 1) %!error interp2 (1:3, 1:3, ones (3,2), 1, 1) %!error interp2 (1:2, 1:2, ones (3,2), 1, 1) %!error interp2 ([1 0 2], 1:3, ones (3,3), 1, 1) %!error interp2 (1:3, [1 0 2], ones (3,3), 1, 1) %!error interp2 (1:2, 1:2, ones (2), ones (2,2), 1) %!error interp2 (1:2, 1:2, ones (2), 1, ones (2,2)) %!error interp2 (1:2, 1:2, ones (2), [1 2 4], [1 2 3], "spline") %!error interp2 (1:2, 1:2, ones (2), [1 2 3], [1 2 4], "spline") %!error interp2 (1, 1, 1, 1, 1, "foobar")