Mercurial > octave-nkf
changeset 19692:5a59c0e1203d
Modified the "extrap" option for interp2, interp3, and interpn (bug #44002).
* interp2.m: allow only "spline" method to set extrapolation values, improved
input validation, updated documentation, testcases added and modified
* interp3.m: allow only "spline" method to set extrapolation values, improved
input validation, updated documentation (no need to implement "pchip"),
testcases added and modified
* interpn.m: allow only "spline" method to set extrapolation values, improved
input validation, updated documentation, testcases added and modified
* __splinen__.m: only apply external supplied extrapolation value
author | Kai T. Ohlhus <k.ohlhus@gmail.com> |
---|---|
date | Thu, 05 Feb 2015 17:58:20 +0100 |
parents | 0cdda69dc2b4 |
children | a9516bc4c55c |
files | scripts/general/interp2.m scripts/general/interp3.m scripts/general/interpn.m scripts/general/private/__splinen__.m |
diffstat | 4 files changed, 167 insertions(+), 211 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/general/interp2.m Tue Feb 03 21:42:32 2015 -0500 +++ b/scripts/general/interp2.m Thu Feb 05 17:58:20 2015 +0100 @@ -70,96 +70,65 @@ ## throughout the curve. ## @end table ## -## If @var{extrap} is the string @qcode{"extrap"}, then extrapolate values -## beyond the endpoints using the current @var{method}. If @var{extrap} is a -## number, then replace values beyond the endpoints with that number. When -## unspecified, @var{extrap} defaults to @code{NA}. Note that if @var{extrap} -## is used, @var{method} must be specified as well. +## @var{extrap} is a scalar number. It replaces values beyond the endpoints +## with @var{extrap}. Note that if @var{extrapval} 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 -## 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) -## FIXME: Need better input validation. -## E.g, interp2 (1,1,1) => A(I): index out of bounds - -function ZI = interp2 (varargin) + narginchk (1, 7); + nargs = nargin; Z = X = Y = XI = YI = n = []; method = "linear"; - extrap = NA; + extrap = []; - switch (nargin) + ## 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--; + 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 - if (ischar (varargin{2})) - [Z, method] = deal (varargin{:}); - n = 1; - else - [Z, n] = deal (varargin{:}); - endif + [Z, n] = deal (varargin{1:nargs}); 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, extrap] = deal (varargin{:}); - endif + [Z, XI, YI] = deal (varargin{1:nargs}); case 5 - if (ischar (varargin{4})) - [Z, XI, YI, method, extrap] = 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, extrap] = deal (varargin{:}); + [X, Y, Z, XI, YI] = deal (varargin{1:nargs}); otherwise print_usage (); endswitch ## Type checking - if (! (ismatrix (Z) && ndims (Z) == 2)) + if (! isnumeric (Z) || isscalar (Z) || ! ismatrix (Z) || ndims (Z) != 2) 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 - if (! ischar (method)) - error ("interp2: METHOD must be a string"); - elseif (method(1) == "*") - warning ("interp2: ignoring unsupported '*' flag to METHOD"); - method(1) = []; - endif - if (isnumeric (extrap) && isscalar (extrap)) - ## Typical case - elseif (strcmp (extrap, "extrap")) - extrap = []; - else - error ('interp2: EXTRAP must be a numeric scalar or "extrap"'); - endif ## Define X, Y, XI, YI if needed [zr, zc] = size (Z); @@ -324,27 +293,6 @@ endif - if (! isempty (extrap)) - ## set points outside the table to 'extrap' - 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 - endif - else ## Check dimensions of XI and YI @@ -354,77 +302,35 @@ error ("interp2: XI and YI must be matrices of equal size"); endif - ## FIXME: Previously used algorithm for cubic. - ## This produced results within a few eps of "spline". - ## Matlab compatibility requires "cubic" to be a C1 algorithm - ## equivalent to "pchip" so this was commented out 2014/03/30. - ## This can be removed completely in the future if no problems are - ## encountered. - #{ - if (strcmp (method, "cubic")) - if (isgriddata (XI) && isgriddata (YI')) - ZI = bicubic (X, Y, Z, XI (1, :), YI (:, 1), extrap); - 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) = extrap; - - else - error ("interp2: input data must have 'meshgrid' format"); - endif - #} - if (strcmp (method, "spline")) if (isgriddata (XI) && isgriddata (YI')) - ZI = __splinen__ ({Y, X}, Z, {YI(:,1), XI(1,:)}, extrap, - "spline"); + ZI = __splinen__ ({Y, X}, Z, {YI(:,1), XI(1,:)}, extrap, "spline"); else error ("interp2: XI, YI must have uniform spacing ('meshgrid' format)"); endif - else - error ("interp2: unrecognized interpolation method '%s'", method); endif + return; # spline doesn't need NA extrapolation value (MATLAB compatibility) + + endif + + ## extrapolation 'extrap' + if (isempty (extrap)) + extrap = NA; + 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 @@ -616,8 +522,11 @@ %! 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", "extrap"), [1,17;468,484]); -%! assert (interp2 (x,y,orig, xi, yi,"nearest", "extrap"), orig([1,end],[1,end])); +%! 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]); + %!test # for values at boundaries %! A = [1,2;3,4]; @@ -654,31 +563,32 @@ %!assert (interp2 (z, [2 3 1], [2 2 2], "spline"), [5 7 3], 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 <Z must be a 2-D matrix> interp2 ({1}) +%!error <Z must be a 2-D matrix> interp2 (1,1,1) %!error <Z must be a 2-D matrix> interp2 (ones (2,2,2)) -%!error <N must be an integer .= 0> interp2 (1, ones (2)) -%!error <N must be an integer .= 0> interp2 (1, -1) -%!error <N must be an integer .= 0> interp2 (1, 1.5) -%!error <METHOD must be a string> interp2 (1, 1, 1, 1, 1, 2) +%!error <N must be an integer .= 0> interp2 (ones (2), ones (2)) +%!error <N must be an integer .= 0> interp2 (ones (2), -1) +%!error <N must be an integer .= 0> interp2 (ones (2), 1.5) %!warning <ignoring unsupported '\*' flag> interp2 (rand (3,3), 1, "*linear"); -%!error <EXTRAP must be a numeric scalar or "extrap"> interp2 (1, 1, 1, 1, 1, 'linear', {1}) -%!error <EXTRAP must be a numeric scalar or "extrap"> interp2 (1, 1, 1, 1, 1, 'linear', ones (2,2)) -%!error <EXTRAP must be a numeric scalar or "extrap"> interp2 (1, 1, 1, 1, 1, 'linear', "abc") -%!error <X, Y must be numeric matrices> interp2 ({1}, 1, 1, 1, 1) -%!error <X, Y must be numeric matrices> interp2 (1, {1}, 1, 1, 1) -%!error <XI, YI must be numeric> interp2 (1, 1, 1, {1}, 1) -%!error <XI, YI must be numeric> interp2 (1, 1, 1, 1, {1}) -%!error <X and Y must be matrices of equal size> interp2 (ones(2,2), 1, 1, 1, 1) -%!error <X and Y must be matrices of equal size> interp2 (ones(2,2), ones(2,3), 1, 1, 1) +%!error <EXTRAP must be a numeric scalar> interp2 (1, 1, 1, 1, 1, 'linear', {1}) +%!error <EXTRAP must be a numeric scalar> interp2 (1, 1, 1, 1, 1, 'linear', ones (2,2)) +%!error <EXTRAP must be a numeric scalar> interp2 (1, 1, 1, 1, 1, 'linear', "abc") +%!error <EXTRAP must be a numeric scalar> interp2 (1, 1, 1, 1, 1, 'linear', "extrap") +%!error <X, Y must be numeric matrices> interp2 ({1}, 1, ones (2), 1, 1) +%!error <X, Y must be numeric matrices> interp2 (1, {1}, ones (2), 1, 1) +%!error <XI, YI must be numeric> interp2 (1, 1, ones (2), {1}, 1) +%!error <XI, YI must be numeric> interp2 (1, 1, ones (2), 1, {1}) +%!error <X and Y must be matrices of equal size> interp2 (ones(2,2), 1, ones (2), 1, 1) +%!error <X and Y must be matrices of equal size> interp2 (ones(2,2), ones(2,3), ones (2), 1, 1) %!error <X and Y size must match the dimensions of Z> interp2 (1:3, 1:3, ones (3,2), 1, 1) %!error <X and Y size must match the dimensions of Z> interp2 (1:2, 1:2, ones (3,2), 1, 1) %!error <X must be strictly monotonic> interp2 ([1 0 2], 1:3, ones (3,3), 1, 1) %!error <Y must be strictly monotonic> interp2 (1:3, [1 0 2], ones (3,3), 1, 1) -%!error <XI and YI must be matrices of equal size> interp2 (1, 1, 1, ones(2,2), 1) -%!error <XI and YI must be matrices of equal size> interp2 (1, 1, 1, 1, ones(2,2)) -%!error <pchip requires at least 2 points> interp2 (1, 1, 1, 1, 1, "pchip") -%!error <cubic requires at least 2 points> interp2 (1, 1, 1, 1, 1, "cubic") -%!error <XI, YI must have uniform spacing> interp2 (1, 1, 1, [1 2 4], [1 2 3], "spline") -%!error <XI, YI must have uniform spacing> interp2 (1, 1, 1, [1 2 3], [1 2 4], "spline") -%!error <unrecognized interpolation method 'foobar'> interp2 (1, 1, 1, 1, 1, "foobar") +%!error <XI and YI must be matrices of equal size> interp2 (1:2, 1:2, ones (2), ones(2,2), 1) +%!error <XI and YI must be matrices of equal size> interp2 (1:2, 1:2, ones (2), 1, ones(2,2)) +%!error <XI, YI must have uniform spacing> interp2 (1:2, 1:2, ones (2), [1 2 4], [1 2 3], "spline") +%!error <XI, YI must have uniform spacing> interp2 (1:2, 1:2, ones (2), [1 2 3], [1 2 4], "spline") +%!error interp2 (1, 1, 1, 1, 1, "foobar")
--- a/scripts/general/interp3.m Tue Feb 03 21:42:32 2015 -0500 +++ b/scripts/general/interp3.m Thu Feb 05 17:58:20 2015 +0100 @@ -61,55 +61,56 @@ ## @item @qcode{"linear"} (default) ## Linear interpolation from nearest neighbors. ## -## @item @qcode{"pchip"} +## @item @qcode{"cubic"} ## Piecewise cubic Hermite interpolating polynomial---shape-preserving ## interpolation with smooth first derivative (not implemented yet). ## -## @item @qcode{"cubic"} -## Cubic interpolation (same as @qcode{"pchip"} [not implemented yet]). -## ## @item @qcode{"spline"} ## Cubic spline interpolation---smooth first and second derivatives ## throughout the curve. ## @end table ## -## If @var{extrapval} is a number, then replace values beyond the endpoints -## with that number. When unspecified, @var{extrapval} defaults to @code{NA}. -## Note that if @var{extrapval} is used, @var{method} must be specified as well. +## @var{extrapval} is a scalar number. It replaces values beyond the endpoints +## with @var{extrapval}. Note that if @var{extrapval} is used, @var{method} +## must be specified as well. If @var{extrapval} is omitted and the +## @var{method} is @qcode{"spline"}, then the extrapolated values of the +## @qcode{"spline"} are used. Otherwise the default @var{extrapval} value for +## any other @var{method} is @qcode{"NA"}. ## @seealso{interp1, interp2, interpn, meshgrid} ## @end deftypefn -## FIXME: Need to validate N argument (maybe change interpn). -## FIXME: Need to add support for 'pchip' method (maybe change interpn). -## FIXME: Need to add support for "extrap" string value (maybe change interpn). +## FIXME: Need to add support for 'cubic' method (maybe change interpn). function vi = interp3 (varargin) + narginchk (1,9); + method = "linear"; - extrapval = NA; + extrapval = []; nargs = nargin; if (nargin < 1 || ! isnumeric (varargin{1})) print_usage (); endif - if (ischar (varargin{end})) - method = varargin{end}; - nargs--; - elseif (nargs > 1 && ischar (varargin{end-1})) - ## FIXME: No support for "extrap" string + if (nargs > 1 && ischar (varargin{end-1})) if (! isnumeric (varargin{end}) || ! isscalar (varargin{end})) error ("interp3: EXTRAPVAL must be a numeric scalar"); endif extrapval = varargin{end}; method = varargin{end-1}; nargs -= 2; + elseif (ischar (varargin{end})) + method = varargin{end}; + nargs--; endif if (method(1) == "*") warning ("interp3: ignoring unsupported '*' flag to METHOD"); method(1) = []; endif + method = validatestring (method, ... + {"nearest", "linear", "cubic", "spline"}); if (nargs < 3) ## Calling form interp3 (v) OR interp3 (v, n) @@ -119,7 +120,12 @@ endif n = varargin(2:nargs); v = permute (v, [2, 1, 3]); - vi = ipermute (interpn (v, n{:}, method, extrapval), [2, 1, 3]); + if (isempty (extrapval)) + vi = interpn (v, n{:}, method); + else + vi = interpn (v, n{:}, method, extrapval); + endif + vi = ipermute (vi, [2, 1, 3]); elseif (nargs == 4 && ! isvector (varargin{1})) ## Calling form interp3 (v, xi, yi, zi) @@ -138,7 +144,12 @@ endfor endif v = permute (v, [2, 1, 3]); - vi = ipermute (interpn (v, xi{:}, method, extrapval), [2, 1, 3]); + if (isempty (extrapval)) + vi = interpn (v, xi{:}, method); + else + vi = interpn (v, xi{:}, method, extrapval); + endif + vi = ipermute (vi, [2, 1, 3]); elseif (nargs == 7) ## Calling form interp3 (x, y, z, v, xi, yi, zi) @@ -167,7 +178,12 @@ endfor endif v = permute (v, [2, 1, 3]); - vi = ipermute (interpn (x{:}, v, xi{:}, method, extrapval), [2, 1, 3]); + if (isempty (extrapval)) + vi = interpn (x{:}, v, xi{:}, method); + else + vi = interpn (x{:}, v, xi{:}, method, extrapval); + endif + vi = ipermute (vi, [2, 1, 3]); else error ("interp3: wrong number or incorrectly formatted input arguments"); @@ -232,6 +248,19 @@ %! vi2 = interpn (v, yi, xi, zi, "nearest", 3); %! assert (vi, vi2); +%!test # extrapolation +%! X=[0,0.5,1]; Y=X; Z=X; +%! V = zeros (3,3,3); +%! V(:,:,1) = [1 3 5; 3 5 7; 5 7 9]; +%! V(:,:,2) = V(:,:,1) + 2; +%! V(:,:,3) = V(:,:,2) + 2; +%! tol = 10 * eps; +%! x=[-0.1,0,0.1]; y=x; z=x; +%! assert(interp3(X,Y,Z,V,x,y,z,"spline"), [-0.2, 1.0, 2.2]',tol); +%! assert(interp3(X,Y,Z,V,x,y,z,"linear"), [NA, 1.0, 2.2]',tol); +%! assert(interp3(X,Y,Z,V,x,y,z,"spline", 0), [0, 1.0, 2.2]',tol); +%! assert(interp3(X,Y,Z,V,x,y,z,"linear", 0), [0, 1.0, 2.2]',tol); + %!shared z, zout, tol %! z = zeros (3, 3, 3); %! zout = zeros (5, 5, 5);
--- a/scripts/general/interpn.m Tue Feb 03 21:42:32 2015 -0500 +++ b/scripts/general/interpn.m Thu Feb 05 17:58:20 2015 +0100 @@ -39,17 +39,21 @@ ## points. This process is performed @var{m} times. If only @var{v} is ## specified, then @var{m} is assumed to be @code{1}. ## -## Method is one of: +## The interpolation @var{method} is one of: ## ## @table @asis ## @item @qcode{"nearest"} ## Return the nearest neighbor. ## -## @item @qcode{"linear"} +## @item @qcode{"linear"} (default) ## Linear interpolation from nearest neighbors. ## +## @item @qcode{"pchip"} +## Piecewise cubic Hermite interpolating polynomial---shape-preserving +## interpolation with smooth first derivative (not implemented yet). +## ## @item @qcode{"cubic"} -## Cubic interpolation from four nearest neighbors (not implemented yet). +## Cubic interpolation (same as @qcode{"pchip"} [not implemented yet]). ## ## @item @qcode{"spline"} ## Cubic spline interpolation---smooth first and second derivatives @@ -58,38 +62,43 @@ ## ## The default method is @qcode{"linear"}. ## -## If @var{extrapval} is the scalar value, use it to replace the values -## beyond the endpoints with that number. If @var{extrapval} is missing, -## assume NA. -## @seealso{interp1, interp2, spline, ndgrid} +## @var{extrapval} is a scalar number. It replaces values beyond the endpoints +## with @var{extrapval}. Note that if @var{extrapval} is used, @var{method} +## must be specified as well. If @var{extrapval} is omitted and the +## @var{method} is @qcode{"spline"}, then the extrapolated values of the +## @qcode{"spline"} are used. Otherwise the default @var{extrapval} value for +## any other @var{method} is @qcode{"NA"}. +## @seealso{interp1, interp2, interp3, spline, ndgrid} ## @end deftypefn function vi = interpn (varargin) method = "linear"; - extrapval = NA; + extrapval = []; nargs = nargin; if (nargin < 1 || ! isnumeric (varargin{1})) print_usage (); endif - if (ischar (varargin{end})) - method = varargin{end}; - nargs -= 1; - elseif (nargs > 1 && ischar (varargin{end - 1})) + if (nargs > 1 && ischar (varargin{end-1})) if (! isnumeric (varargin{end}) || ! isscalar (varargin{end})) - error ("interpn: extrapal is expected to be a numeric scalar"); + error ("interpn: EXTRAPVAL must be a numeric scalar"); endif - method = varargin{end - 1}; extrapval = varargin{end}; + method = varargin{end-1}; nargs -= 2; + elseif (ischar (varargin{end})) + method = varargin{end}; + nargs--; endif if (method(1) == "*") warning ("interpn: ignoring unsupported '*' flag to METHOD"); method(1) = []; endif + method = validatestring (method, ... + {"nearest", "linear", "pchip", "cubic", "spline"}); if (nargs < 3) v = varargin{1}; @@ -159,6 +168,9 @@ if (strcmp (method, "linear")) vi = __lin_interpn__ (x{:}, v, y{:}); + if (isempty (extrapval)) + extrapval = NA; + endif vi(isna (vi)) = extrapval; elseif (strcmp (method, "nearest")) yshape = size (y{1}); @@ -176,6 +188,9 @@ for i = 1 : nd idx |= y{i} < min (x{i}(:)) | y{i} > max (x{i}(:)); endfor + if (isempty (extrapval)) + extrapval = NA; + endif vi(idx) = extrapval; vi = reshape (vi, yshape); elseif (strcmp (method, "spline"))
--- a/scripts/general/private/__splinen__.m Tue Feb 03 21:42:32 2015 -0500 +++ b/scripts/general/private/__splinen__.m Thu Feb 05 17:58:20 2015 +0100 @@ -38,10 +38,12 @@ endfor [xi{:}] = ndgrid (cellfun (@(x) x(:), xi, "uniformoutput", false){:}); - idx = zeros (size (xi{1})); - for i = 1 : length (x) - idx |= xi{i} < min (x{i}(:)) | xi{i} > max (x{i}(:)); - endfor - yi(idx) = extrapval; + if (!isempty (extrapval)) + idx = zeros (size (xi{1})); + for i = 1 : length (x) + idx |= xi{i} < min (x{i}(:)) | xi{i} > max (x{i}(:)); + endfor + yi(idx) = extrapval; + endif endfunction