Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project!

## [Octave-cvsupdate] octave-forge/main/general interp2.m,1.3,1.4

 [Octave-cvsupdate] octave-forge/main/general interp2.m,1.3,1.4 From: Paul Kienzle - 2005-03-03 04:47:04 ```Update of /cvsroot/octave/octave-forge/main/general In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv17520 Modified Files: interp2.m Log Message: [Thomas Weber and Paul Kienzle] simplify and add test cases Index: interp2.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/general/interp2.m,v retrieving revision 1.3 retrieving revision 1.4 diff -u -d -r1.3 -r1.4 --- interp2.m 24 May 2004 11:34:03 -0000 1.3 +++ interp2.m 3 Mar 2005 04:46:55 -0000 1.4 @@ -21,8 +21,8 @@ ## @deftypefnx {Function File} {@var{zi}=} interp2 (... , '@var{method}') ## 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}. If they are matrices they -## must have the 'meshgrid' format. +## must correspondent to the size of @var{Z}. @var{X} and @var{Y} must be +## monotonic. If they are matrices they must have the 'meshgrid' format. ## ## ZI = interp2 (X, Y, Z, XI, YI, ...) returns a matrix corresponding ## to the points described by the matrices @var{XI}, @var{YI}. @@ -41,103 +41,72 @@ ## @end deftypefn ## Author: Kai Habel +## 2005-03-02 Thomas Weber +## * Add test cases +## 2005-03-02 Paul Kienzle +## * Simplify -function ZI = interp2 (X, Y, Z, XI, YI, method) - - if (nargin > 6 || nargin == 0) - usage ("interp2 (X, Y, Z, XI, YI, method)"); - endif +function ZI = interp2 (varargin) + Z = X = Y = xi = yi = []; + n = 1; + method = "linear"; - if (nargin > 4) - if (is_vector (X) && is_vector (Y)) - [rz, cz] = size (Z); - if (rz != length (Y) || cz != length (X)) - error ("length of X and Y must match the size of Z"); - endif - [X, Y] = meshgrid (X, Y); - elseif !( (size (X) == size (Y)) && (size (X) == size (Z)) ) - error ("X,Y and Z must be matrices of same size"); - endif - - elseif (((nargin == 4) || (nargin == 3)) && !isstr (Z)) - - if (nargin == 4) - if (isstr (XI)) - method = XI; - else - usage("interp2 (z,xi,yi,'format'"); - endif + switch nargin + case 1 + Z = varargin{1}; + case 2 + if isstr(varargin{2}) + [Z,method] = deal(varargin{:}); + else + [Z,n] = deal(varargin{:}); endif - XI = Y; - YI = Z; - Z = X; - [X, Y] = meshgrid(1:columns(Z), 1:rows(Z)); - else - if (nargin == 1) - n = 1; - elseif (nargin == 2) - if (isstr (Y)) - method = Y; - n = 1; - elseif (is_scalar (Y)) - n = Y; - endif + case 3 + if isstr(varargin{3}) + [Z,n,method] = deal(varargin{:}); else - n = Y; - if (isstr (Z)) - method = Z; - endif + [Z,XI,YI] = deal(varargin{:}); endif - Z = X; - [zr, zc] = size (Z); - [X, Y] = meshgrid (1:zc, 1:zr); - xi = linspace (1, zc, pow2 (n) * (zc - 1) + 1); - yi = linspace (1, zr, pow2 (n) * (zr - 1) + 1); - [XI, YI] = meshgrid (xi, yi); - endif - - if (! exist ("method")) - method = "linear"; - endif + case 4 + [Z,XI,YI,method] = deal(varargin{:}); + case 5 + [X,Y,Z,XI,YI] = deal(varargin{:}); + case 6 + [X,Y,Z,XI,YI,method] = deal(varargin{:}); + otherwise + usage ("interp2 (X, Y, Z, XI, YI, method)"); + endswitch - xtable = X(1, :); - ytable = Y(:, 1); + % Type checking. + if !ismatrix(Z), error("interp2 expected matrix Z"); endif + if !isscalar(n), error("interp2 expected scalar n"); endif + if !isnumeric(X) || !isnumeric(Y), error("interp2 expected numeric X,Y"); endif + if !isnumeric(XI) || !isnumeric(YI), error("interp2 expected numeric XI,YI"); endif + if !isstr(method), error("interp2 expected string 'method'"); endif - if (is_vector (XI) && is_vector (YI)) + % Define X,Y,XI,YI if needed + [zr, zc] = size (Z); + if isempty(X), X=[1:zc]; Y=[1:zr]; endif + if isempty(XI), p=2^n; XI=[p:p*zc]/p; YI=[p:p*zr]/p; endif + + % Are X and Y vectors? => produce a grid from them + if isvector (X) && isvector (Y) + [X, Y] = meshgrid (X, Y); + elseif ! all(size (X) == size (Y)) + error ("X and Y must be matrices of same size"); + endif + if any(size (X) != size (Z)) + error ("X and Y size must match Z dimensions"); + endif + + % Are XI and YI vectors? => produce a grid from them + if isvector (XI) && isvector (YI) [XI, YI] = meshgrid (XI, YI); - elseif (! (size (XI) == size (YI))) + elseif any(size (XI) != size (YI)) error ("XI and YI must be matrices of same size"); endif - - ytlen = length (ytable); - xtlen = length (xtable); - - ## get x index of values in XI - xtable(xtlen) *= (1 + eps); - xtable(xtlen) > XI(1, :); - [m, n] = sort ([xtable(:); XI(1, :)']); - o = cumsum (n <= xtlen); - xidx = o([find(n > xtlen)])'; - - ## get y index of values in YI - ytable(ytlen) *= (1 + eps); - [m, n]=sort ([ytable(:); YI(:, 1)]); - o = cumsum (n <= ytlen); - yidx = o([find(n > ytlen)]); - - [zr, zc] = size (Z); - - ## mark values outside the lookup table - xfirst_val = find (XI(1,:) < xtable(1)); - xlast_val = find (XI(1,:) > xtable(xtlen)); - yfirst_val = find (YI(:,1) < ytable(1)); - ylast_val = find (YI(:,1) > ytable(ytlen)); - - ## set value outside the table preliminary to min max index - yidx(yfirst_val) = 1; - xidx(xfirst_val) = 1; - yidx(ylast_val) = zr - 1; - xidx(xlast_val) = zc - 1; + + xidx = lookup(X(1, 2:end-1), XI(1,:))' + 1; + yidx = lookup(Y(2:end-1, 1), YI(:,1)) + 1; if strcmp (method, "linear") ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy @@ -145,32 +114,31 @@ ## 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; + 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(yidx, xidx)) ./ (X(yidx, xidx + 1) - X(yidx, xidx)); - Ysc = (YI .- Y(yidx, xidx)) ./ (Y(yidx + 1, xidx) - Y(yidx, xidx)); + Xsc = (XI - X(yidx, xidx)) ./ (X(yidx, xidx + 1) - X(yidx, xidx)); + Ysc = (YI - Y(yidx, xidx)) ./ (Y(yidx + 1, xidx) - Y(yidx, xidx)); ## apply plane equation - ZI = a(yidx, xidx) .+ b(yidx, xidx) .* Xsc \ - .+ c(yidx, xidx) .* Ysc .+ d(yidx, xidx) .* Xsc .* Ysc; + ZI = a(yidx, xidx) + b(yidx, xidx) .* Xsc \ + + c(yidx, xidx) .* Ysc + d(yidx, xidx) .* Xsc .* Ysc; elseif strcmp (method, "nearest") - i = XI(1, :) - xtable(xidx) > xtable(xidx + 1) - XI(1, :); - j = YI(:, 1) - ytable(yidx) > ytable(yidx + 1) - YI(:, 1); - ZI = Z(yidx + j, xidx + i); + xtable = X(1, :); + ytable = Y(:, 1); + i = (XI(1, :) - xtable(xidx) > xtable(xidx + 1) - XI(1, :)); + j = (YI(:, 1) - ytable(yidx) > ytable(yidx + 1) - YI(:, 1)); + ZI = Z(yidx + j(:), xidx + i(:)); else error ("interpolation method not (yet) supported"); endif ## set points outside the table to NaN - if (! (isempty (xfirst_val) && isempty (xlast_val))) - ZI(:, [xfirst_val, xlast_val]) = NaN; - endif - if (! (isempty (yfirst_val) && isempty (ylast_val))) - ZI([yfirst_val; ylast_val], :) = NaN; - endif + ZI(:, XI(1,:) < X(1,1) | XI(1,:) > X(1,end)) = NaN; + ZI(YI(:,1) < Y(1,1) | YI(:,1) > Y(end,1), :) = NaN; + endfunction %!demo @@ -199,3 +167,49 @@ %! mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic')); %! [x,y] = meshgrid(x,y); gset nohidden3d; %! 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.2,2, 1.5]; +%! yi = [6.2, 4.0, 5.0]'; +%! +%! Expected = \ +%! [243, 245.4, 243.9; +%! 65.6, 68, 66.5; +%! 126.6, 129, 127.5]; +%! Result = interp2(x,y,Orig, xi, yi); +%! +%! assert(Result, Expected, 1000*eps); + +%!test +%! % 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]'; +%! Expected = [nan,nan; nan,nan]; +%! Result = interp2(x,y,Orig, xi, yi); +%! +%! assert(Result, Expected); + +%!test +%! % test for values at boundaries +%! A=[1,2;3,4]; +%! x=[0,1]; +%! y=[2,3]; +%! xi=x; +%! yi=y; +%! Expected = A; +%! +%! Result = interp2(x,y,A,xi,yi,'linear'); +%! assert(Result, Expected); +%! +%! Result = interp2(x,y,A,xi,yi,'nearest'); +%! assert(Result, Expected); + ```

 [Octave-cvsupdate] octave-forge/main/general interp2.m,1.3,1.4 From: Paul Kienzle - 2005-03-03 04:47:04 ```Update of /cvsroot/octave/octave-forge/main/general In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv17520 Modified Files: interp2.m Log Message: [Thomas Weber and Paul Kienzle] simplify and add test cases Index: interp2.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/general/interp2.m,v retrieving revision 1.3 retrieving revision 1.4 diff -u -d -r1.3 -r1.4 --- interp2.m 24 May 2004 11:34:03 -0000 1.3 +++ interp2.m 3 Mar 2005 04:46:55 -0000 1.4 @@ -21,8 +21,8 @@ ## @deftypefnx {Function File} {@var{zi}=} interp2 (... , '@var{method}') ## 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}. If they are matrices they -## must have the 'meshgrid' format. +## must correspondent to the size of @var{Z}. @var{X} and @var{Y} must be +## monotonic. If they are matrices they must have the 'meshgrid' format. ## ## ZI = interp2 (X, Y, Z, XI, YI, ...) returns a matrix corresponding ## to the points described by the matrices @var{XI}, @var{YI}. @@ -41,103 +41,72 @@ ## @end deftypefn ## Author: Kai Habel +## 2005-03-02 Thomas Weber +## * Add test cases +## 2005-03-02 Paul Kienzle +## * Simplify -function ZI = interp2 (X, Y, Z, XI, YI, method) - - if (nargin > 6 || nargin == 0) - usage ("interp2 (X, Y, Z, XI, YI, method)"); - endif +function ZI = interp2 (varargin) + Z = X = Y = xi = yi = []; + n = 1; + method = "linear"; - if (nargin > 4) - if (is_vector (X) && is_vector (Y)) - [rz, cz] = size (Z); - if (rz != length (Y) || cz != length (X)) - error ("length of X and Y must match the size of Z"); - endif - [X, Y] = meshgrid (X, Y); - elseif !( (size (X) == size (Y)) && (size (X) == size (Z)) ) - error ("X,Y and Z must be matrices of same size"); - endif - - elseif (((nargin == 4) || (nargin == 3)) && !isstr (Z)) - - if (nargin == 4) - if (isstr (XI)) - method = XI; - else - usage("interp2 (z,xi,yi,'format'"); - endif + switch nargin + case 1 + Z = varargin{1}; + case 2 + if isstr(varargin{2}) + [Z,method] = deal(varargin{:}); + else + [Z,n] = deal(varargin{:}); endif - XI = Y; - YI = Z; - Z = X; - [X, Y] = meshgrid(1:columns(Z), 1:rows(Z)); - else - if (nargin == 1) - n = 1; - elseif (nargin == 2) - if (isstr (Y)) - method = Y; - n = 1; - elseif (is_scalar (Y)) - n = Y; - endif + case 3 + if isstr(varargin{3}) + [Z,n,method] = deal(varargin{:}); else - n = Y; - if (isstr (Z)) - method = Z; - endif + [Z,XI,YI] = deal(varargin{:}); endif - Z = X; - [zr, zc] = size (Z); - [X, Y] = meshgrid (1:zc, 1:zr); - xi = linspace (1, zc, pow2 (n) * (zc - 1) + 1); - yi = linspace (1, zr, pow2 (n) * (zr - 1) + 1); - [XI, YI] = meshgrid (xi, yi); - endif - - if (! exist ("method")) - method = "linear"; - endif + case 4 + [Z,XI,YI,method] = deal(varargin{:}); + case 5 + [X,Y,Z,XI,YI] = deal(varargin{:}); + case 6 + [X,Y,Z,XI,YI,method] = deal(varargin{:}); + otherwise + usage ("interp2 (X, Y, Z, XI, YI, method)"); + endswitch - xtable = X(1, :); - ytable = Y(:, 1); + % Type checking. + if !ismatrix(Z), error("interp2 expected matrix Z"); endif + if !isscalar(n), error("interp2 expected scalar n"); endif + if !isnumeric(X) || !isnumeric(Y), error("interp2 expected numeric X,Y"); endif + if !isnumeric(XI) || !isnumeric(YI), error("interp2 expected numeric XI,YI"); endif + if !isstr(method), error("interp2 expected string 'method'"); endif - if (is_vector (XI) && is_vector (YI)) + % Define X,Y,XI,YI if needed + [zr, zc] = size (Z); + if isempty(X), X=[1:zc]; Y=[1:zr]; endif + if isempty(XI), p=2^n; XI=[p:p*zc]/p; YI=[p:p*zr]/p; endif + + % Are X and Y vectors? => produce a grid from them + if isvector (X) && isvector (Y) + [X, Y] = meshgrid (X, Y); + elseif ! all(size (X) == size (Y)) + error ("X and Y must be matrices of same size"); + endif + if any(size (X) != size (Z)) + error ("X and Y size must match Z dimensions"); + endif + + % Are XI and YI vectors? => produce a grid from them + if isvector (XI) && isvector (YI) [XI, YI] = meshgrid (XI, YI); - elseif (! (size (XI) == size (YI))) + elseif any(size (XI) != size (YI)) error ("XI and YI must be matrices of same size"); endif - - ytlen = length (ytable); - xtlen = length (xtable); - - ## get x index of values in XI - xtable(xtlen) *= (1 + eps); - xtable(xtlen) > XI(1, :); - [m, n] = sort ([xtable(:); XI(1, :)']); - o = cumsum (n <= xtlen); - xidx = o([find(n > xtlen)])'; - - ## get y index of values in YI - ytable(ytlen) *= (1 + eps); - [m, n]=sort ([ytable(:); YI(:, 1)]); - o = cumsum (n <= ytlen); - yidx = o([find(n > ytlen)]); - - [zr, zc] = size (Z); - - ## mark values outside the lookup table - xfirst_val = find (XI(1,:) < xtable(1)); - xlast_val = find (XI(1,:) > xtable(xtlen)); - yfirst_val = find (YI(:,1) < ytable(1)); - ylast_val = find (YI(:,1) > ytable(ytlen)); - - ## set value outside the table preliminary to min max index - yidx(yfirst_val) = 1; - xidx(xfirst_val) = 1; - yidx(ylast_val) = zr - 1; - xidx(xlast_val) = zc - 1; + + xidx = lookup(X(1, 2:end-1), XI(1,:))' + 1; + yidx = lookup(Y(2:end-1, 1), YI(:,1)) + 1; if strcmp (method, "linear") ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy @@ -145,32 +114,31 @@ ## 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; + 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(yidx, xidx)) ./ (X(yidx, xidx + 1) - X(yidx, xidx)); - Ysc = (YI .- Y(yidx, xidx)) ./ (Y(yidx + 1, xidx) - Y(yidx, xidx)); + Xsc = (XI - X(yidx, xidx)) ./ (X(yidx, xidx + 1) - X(yidx, xidx)); + Ysc = (YI - Y(yidx, xidx)) ./ (Y(yidx + 1, xidx) - Y(yidx, xidx)); ## apply plane equation - ZI = a(yidx, xidx) .+ b(yidx, xidx) .* Xsc \ - .+ c(yidx, xidx) .* Ysc .+ d(yidx, xidx) .* Xsc .* Ysc; + ZI = a(yidx, xidx) + b(yidx, xidx) .* Xsc \ + + c(yidx, xidx) .* Ysc + d(yidx, xidx) .* Xsc .* Ysc; elseif strcmp (method, "nearest") - i = XI(1, :) - xtable(xidx) > xtable(xidx + 1) - XI(1, :); - j = YI(:, 1) - ytable(yidx) > ytable(yidx + 1) - YI(:, 1); - ZI = Z(yidx + j, xidx + i); + xtable = X(1, :); + ytable = Y(:, 1); + i = (XI(1, :) - xtable(xidx) > xtable(xidx + 1) - XI(1, :)); + j = (YI(:, 1) - ytable(yidx) > ytable(yidx + 1) - YI(:, 1)); + ZI = Z(yidx + j(:), xidx + i(:)); else error ("interpolation method not (yet) supported"); endif ## set points outside the table to NaN - if (! (isempty (xfirst_val) && isempty (xlast_val))) - ZI(:, [xfirst_val, xlast_val]) = NaN; - endif - if (! (isempty (yfirst_val) && isempty (ylast_val))) - ZI([yfirst_val; ylast_val], :) = NaN; - endif + ZI(:, XI(1,:) < X(1,1) | XI(1,:) > X(1,end)) = NaN; + ZI(YI(:,1) < Y(1,1) | YI(:,1) > Y(end,1), :) = NaN; + endfunction %!demo @@ -199,3 +167,49 @@ %! mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic')); %! [x,y] = meshgrid(x,y); gset nohidden3d; %! 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.2,2, 1.5]; +%! yi = [6.2, 4.0, 5.0]'; +%! +%! Expected = \ +%! [243, 245.4, 243.9; +%! 65.6, 68, 66.5; +%! 126.6, 129, 127.5]; +%! Result = interp2(x,y,Orig, xi, yi); +%! +%! assert(Result, Expected, 1000*eps); + +%!test +%! % 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]'; +%! Expected = [nan,nan; nan,nan]; +%! Result = interp2(x,y,Orig, xi, yi); +%! +%! assert(Result, Expected); + +%!test +%! % test for values at boundaries +%! A=[1,2;3,4]; +%! x=[0,1]; +%! y=[2,3]; +%! xi=x; +%! yi=y; +%! Expected = A; +%! +%! Result = interp2(x,y,A,xi,yi,'linear'); +%! assert(Result, Expected); +%! +%! Result = interp2(x,y,A,xi,yi,'nearest'); +%! assert(Result, Expected); + ```