# HG changeset patch # User Jaroslav Hajek # Date 1205495696 -3600 # Node ID 4586400dfda6ed0b289329a3b23bd9031b542477 # Parent 67403cfad8d7bfd7db1e65a66b1f53ba6cd21990 add pchip method to interp2 diff -r 67403cfad8d7 -r 4586400dfda6 scripts/ChangeLog --- a/scripts/ChangeLog Thu Mar 13 13:35:27 2008 -0400 +++ b/scripts/ChangeLog Fri Mar 14 12:54:56 2008 +0100 @@ -1,3 +1,9 @@ 2008-03-12 David Bateman + + * general/interp2.m: Added support for pchip bicubic interpolation. + Also simplified code and added support for natural extrapolation via + "extrap". + 2008-03-12 David Bateman * geometry/griddata3.m: Use griddatan and not griddata diff -r 67403cfad8d7 -r 4586400dfda6 scripts/general/interp2.m --- a/scripts/general/interp2.m Thu Mar 13 13:35:27 2008 -0400 +++ b/scripts/general/interp2.m Fri Mar 14 12:54:56 2008 +0100 @@ -87,6 +87,15 @@ ## * Rather big modification (XI,YI no longer need to be ## "meshgridded") to be consistent with the help message ## above and for compatibility. +## 2008-02-27 Jaroslav Hajek +## * Added support for "pchip" monotonicity-preserving +## interpolation via __pchip_deriv__ +## * Support "extrap" to do extrapolation. +## * some cleanup. Remove unneccessary meshgridding X and Y in +## nearest and linear methods, because only leading row and columns +## was used anyway. +## * duplicated demos based on peaks() function. +## function ZI = interp2 (varargin) Z = X = Y = XI = YI = n = []; @@ -138,7 +147,9 @@ function ZI = interp2 (varargin) if (!ischar (method)) error ("interp2 expected string 'method'"); endif - if (!isscalar (extrapval)) + if (ischar (extrapval) || strcmp (extrapval, "extrap")) + extrapval = []; + elseif (!isscalar (extrapval)) error ("interp2 expected n extrapval"); endif @@ -161,15 +172,18 @@ function ZI = interp2 (varargin) endif - if (strcmp (method, "linear") || strcmp (method, "nearest")) + if (strcmp (method, "linear") || strcmp (method, "nearest") ... + || strcmp (method, "pchip")) ## If X and Y vectors produce a grid from them if (isvector (X) && isvector (Y)) - [X, Y] = meshgrid (X, Y); - elseif (! size_equal (X, Y)) + X = X(:); Y = Y(:); + elseif (size_equal (X, Y)) + X = X(1,:)'; Y = Y(:,1); + else error ("X and Y must be matrices of same size"); endif - if (! size_equal (X, Z)) + if (columns (Z) != length (X) || rows (Z) != length (Y)) error ("X and Y size must match Z dimensions"); endif @@ -181,12 +195,8 @@ function ZI = interp2 (varargin) error ("XI and YI must be matrices of same size"); endif - shape = size (XI); - XI = reshape (XI, 1, prod (shape)); - YI = reshape (YI, 1, prod (shape)); - - xidx = lookup (X(1, 2:end-1), XI) + 1; - yidx = lookup (Y(2:end-1, 1), YI) + 1; + xidx = lookup (X(2:end-1), XI) + 1; + yidx = lookup (Y(2:end-1), YI) + 1; if (strcmp (method, "linear")) ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy @@ -202,41 +212,98 @@ function ZI = interp2 (varargin) idx = sub2ind (size (a), yidx, xidx); ## scale XI, YI values to a 1-spaced grid - Xsc = (XI - X(1, xidx)) ./ (X(1, xidx + 1) - X(1, xidx)); - Ysc = (YI - Y(yidx, 1)') ./ (Y(yidx + 1, 1) - Y(yidx, 1))'; + Xsc = (XI - X(xidx)) ./ (X(xidx + 1) - X(xidx)); + Ysc = (YI - Y(yidx)) ./ (Y(yidx + 1) - Y(yidx)); ## apply plane equation ZI = a(idx) + b(idx).*Xsc + c(idx).*Ysc + d(idx).*Xsc.*Ysc; elseif (strcmp (method, "nearest")) - xtable = X(1, :); - ytable = Y(:, 1)'; - ii = (XI - xtable(xidx) > xtable(xidx + 1) - XI); - jj = (YI - ytable(yidx) > ytable(yidx + 1) - YI); + ii = (XI - X(xidx) > X(xidx + 1) - XI); + jj = (YI - X(yidx) > X(yidx + 1) - YI); idx = sub2ind (size (Z), yidx+jj, xidx+ii); ZI = Z(idx); - endif - - ## set points outside the table to 'extrapval' - if (X (1, 1) < X (1, end)) - if (Y (1, 1) < Y (end, 1)) - ZI (XI < X(1,1) | XI > X(1,end) | YI < Y(1,1) | YI > Y(end,1)) = ... - extrapval; - else - ZI (XI < X(1,1) | XI > X(1,end) | YI < Y(end,1) | YI > Y(1,1)) = ... - extrapval; - endif - else - if (Y (1, 1) < Y (end, 1)) - ZI (XI < X(1,end) | XI > X(1,1) | YI < Y(1,1) | YI > Y(end,1)) = ... - extrapval; - else - ZI (XI < X(1,end) | XI > X(1,1) | YI < Y(end,1) | YI > Y(1,1)) = ... - extrapval; - endif - endif - - ZI = reshape (ZI, shape); + + elseif (strcmp (method, "pchip")) + + if (length (X) < 2 || length (Y) < 2) + error ("interp2: pchip2 requires at least 2 points in each dimension") + 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 + + if (! isempty (extrapval)) + ## set points outside the table to 'extrapval' + 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)) = ... + extrapval; + else + ZI (XI < X(1) | XI > X(end) | YI < Y(end) | YI > Y(1)) = ... + extrapval; + endif + else + if (Y (1) < Y (end)) + ZI (XI < X(end) | XI > X(1) | YI < Y(1) | YI > Y(end)) = ... + extrapval; + else + ZI (XI < X(1,end) | XI > X(1) | YI < Y(end) | YI > Y(1)) = ... + extrapval; + endif + endif + endif + else ## If X and Y vectors produce a grid from them @@ -285,6 +352,15 @@ endfunction %! 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); @@ -294,19 +370,64 @@ endfunction %! 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; diff -r 67403cfad8d7 -r 4586400dfda6 src/ChangeLog --- a/src/ChangeLog Thu Mar 13 13:35:27 2008 -0400 +++ b/src/ChangeLog Fri Mar 14 12:54:56 2008 +0100 @@ -1,3 +1,9 @@ 2008-03-13 John W. Eaton + + * DLD-FUNCTIONS/__pchip_deriv__.cc (F__pchip_deriv__): + added support for computing pchip derivatives along rows of matrix, + eliminated redundant copying. + 2008-03-13 John W. Eaton * ov-usr-fcn.cc (octave_user_function::octave_user_function): diff -r 67403cfad8d7 -r 4586400dfda6 src/DLD-FUNCTIONS/__pchip_deriv__.cc --- a/src/DLD-FUNCTIONS/__pchip_deriv__.cc Thu Mar 13 13:35:27 2008 -0400 +++ b/src/DLD-FUNCTIONS/__pchip_deriv__.cc Fri Mar 14 12:54:56 2008 +0100 @@ -1,6 +1,7 @@ /* Copyright (C) 2002, 2006, 2007 Kai Habel +Copyright (C) 2008 Jaroslav Hajek This file is part of Octave. @@ -20,6 +21,9 @@ along with Octave; see the file COPYING. */ +// Jaroslav Hajek, Feb 2008: handle row-wise derivatives, +// use const pointers to avoid unnecessary copying. + #ifdef HAVE_CONFIG_H #include #endif @@ -34,7 +38,7 @@ extern "C" extern "C" { F77_RET_T - F77_FUNC (dpchim, DPCHIM) (const octave_idx_type& n, double *x, double *f, + F77_FUNC (dpchim, DPCHIM) (const octave_idx_type& n, const double *x, const double *f, double *d, const octave_idx_type &incfd, octave_idx_type *ierr); } @@ -44,14 +48,16 @@ extern "C" DEFUN_DLD (__pchip_deriv__, args, , "-*- texinfo -*-\n\ address@hidden {Loadable Function} {} __pchip_deriv__ (@var{x}, @var{y})\n\ address@hidden {Loadable Function} {} __pchip_deriv__ (@var{x}, @var{y}, @var{dim})\n\ Undocumented internal function.\n\ @end deftypefn") { octave_value retval; const int nargin = args.length (); - if (nargin == 2) + bool rows = (nargin == 3 && args (2).uint_value() == 2); + + if (nargin >= 2) { ColumnVector xvec (args(0).vector_value ()); Matrix ymat (args(1).matrix_value ()); @@ -60,25 +66,29 @@ Undocumented internal function.\n\ octave_idx_type nyr = ymat.rows (); octave_idx_type nyc = ymat.columns (); - if (nx != nyr) + if (nx != (rows ? nyc : nyr)) { - error ("number of rows for x and y must match"); + error ("__pchip_deriv__: dimension mismatch"); return retval; } - ColumnVector dvec (nx), yvec (nx); + const double *yvec = ymat.data (); Matrix dmat (nyr, nyc); + double *dvec = dmat.fortran_vec (); octave_idx_type ierr; - const octave_idx_type incfd = 1; - for (int c = 0; c < nyc; c++) + const octave_idx_type incfd = rows ? nyr : 1; + const octave_idx_type inc = rows ? 1 : nyc; + + for (int i = (rows ? nyr : nyc); i > 0; i--) { - for (int r = 0; r < nx; r++) - yvec(r) = ymat(r,c); - F77_FUNC (dpchim, DPCHIM) (nx, xvec.fortran_vec (), - yvec.fortran_vec (), - dvec.fortran_vec (), incfd, &ierr); + F77_FUNC (dpchim, DPCHIM) (nx, xvec.data (), + yvec, + dvec, incfd, &ierr); + + yvec += inc; + dvec += inc; if (ierr < 0) { @@ -86,8 +96,6 @@ Undocumented internal function.\n\ return retval; } - for (int r=0; r