octave-antonio: scripts/general/interp2.m comparison

comparison scripts/general/interp2.m @ 11587:c792872f8942

all script files: untabify and strip trailing whitespace

author	John W. Eaton <jwe@octave.org>
date	Thu, 20 Jan 2011 17:35:29 -0500
parents	702dbd0c53f5
children	b0084095098e

comparison

equal deleted inserted replaced

-:12df7854fa7c
+:c792872f8942
 ## @deftypefnx {Function File} {@var{zi} =} interp2 (@dots{}, @var{method}, @var{extrapval})
 ##
 ## Two-dimensional interpolation.  @var{x}, @var{y} and @var{z} describe a
 ## surface function.  If @var{x} and @var{y} are vectors their length
 ## must correspondent to the size of @var{z}.  @var{x} and @var{y} must be
 ## monotonic.  If they are matrices they must have the @code{meshgrid}
 ## format.
 ##
 ## @table @code
 ## @item interp2 (@var{x}, @var{y}, @var{Z}, @var{xi}, @var{yi}, @dots{})
 ## Returns a matrix corresponding to the points described by the
 ## matrices @var{xi}, @var{yi}.
 ##
 ## If the last argument is a string, the interpolation method can
 ## be specified.  The method can be 'linear', 'nearest' or 'cubic'.
 ## If it is omitted 'linear' interpolation is assumed.
 ##
 ## @item interp2 (@var{z}, @var{xi}, @var{yi})
 ## Assumes @code{@var{x} = 1:rows (@var{z})} and @code{@var{y} =
 ## 1:columns (@var{z})}
 ##
 ## @item interp2 (@var{z}, @var{n})
 ## Interleaves the matrix @var{z} n-times.  If @var{n} is omitted a value
 ## of @code{@var{n} = 1} is assumed.
 ## @end table
 ##
 ## The variable @var{method} defines the method to use for the
 ## interpolation.  It can take one of the following values
 ##
 ## @table @asis
 ## @item 'nearest'
 ## Return the nearest neighbor.
 ##
 ## Cubic spline interpolation--smooth first and second derivatives
 ## throughout the curve.
 ## @end table
 ##
 ## If a scalar value @var{extrapval} is defined as the final value, then
 ## values outside the mesh as set to this value.  Note that in this case
 ## @var{method} must be defined as well.  If @var{extrapval} is not
 ## defined then NA is assumed.
 ##
 ## @seealso{interp1}
 ## @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
 if (ischar (varargin{4}))
 [Z, XI, YI, method, extrapval] = 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, extrapval] = deal (varargin{:});
 otherwise
 print_usage ();
 endswitch
 ## Type checking.
 if (!ismatrix (Z))
 error ("interp2: Z must be a matrix");
 endif
 if (!isempty (n) && !isscalar (n))
 error ("interp2: N must be a scalar");
 endif
 if (!ischar (method))
 error ("interp2: METHOD must be a string");
 endif
 if (ischar (extrapval) || strcmp (extrapval, "extrap"))
 extrapval = [];
 elseif (!isscalar (extrapval))
 error ("interp2: EXTRAPVAL must be a scalar");
 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))
 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 (strcmp (method, "linear") || strcmp (method, "nearest") ...
 || strcmp (method, "pchip"))
 ## 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;
 if (0 && isgriddata (XI) && isgriddata (YI'))
 ZI = bicubic (X, Y, Z, XI (1, :), YI (:, 1), extrapval);
 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  .* 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) = extrapval;
 else
 error ("interp2: input data must have `meshgrid' format");
 endif
 elseif (strcmp (method, "spline"))
 if (isgriddata (XI) && isgriddata (YI'))
 ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrapval,
 "spline");
 else
 error ("interp2: input data must have `meshgrid' format");
 endif
 else
 %! 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
 %! [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
 %! 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
 %! [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
 %! 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
 %! [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
 %! 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
 %! [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
 %! 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
 %! [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;
 %!test % simple test
 %!  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]);
 %!  assert(interp2(x,y,Orig, xi, yi),[NA,NA;NA,NA]);
 %!  assert(interp2(x,y,Orig, xi, yi,'linear', 0),[0,0;0,0]);
 %!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'

Mercurial > octave-antonio

comparison scripts/general/interp2.m @ 11587:c792872f8942