Mercurial > octave-nkf

changeset 19692:5a59c0e1203d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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