Mercurial > octave-nkf
changeset 8712:010e15c7be9a
support pchip method in interp2
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Fri, 14 Mar 2008 12:54:56 +0100 |
| parents | 822cab55ca85 |
| children | 02d4c764de61 |
| files | scripts/ChangeLog scripts/general/interp2.m src/ChangeLog src/DLD-FUNCTIONS/__pchip_deriv__.cc |
| diffstat | 4 files changed, 208 insertions(+), 62 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog Mon Feb 09 16:53:12 2009 -0500 +++ b/scripts/ChangeLog Fri Mar 14 12:54:56 2008 +0100 @@ -1,3 +1,9 @@ +2008-03-14 Jaroslav Hajek <highegg@gmail.com> + + * general/interp2.m: Added support for pchip bicubic interpolation. + Also simplified code and added support for natural extrapolation via + "extrap". + 2009-02-09 John W. Eaton <jwe@octave.org> * miscellaneous/Makefile.in (SOURCES): Include __xzip__.m in the list.
--- a/scripts/general/interp2.m Mon Feb 09 16:53:12 2009 -0500 +++ b/scripts/general/interp2.m Fri Mar 14 12:54:56 2008 +0100 @@ -1,4 +1,5 @@ ## Copyright (C) 2000, 2006, 2007 Kai Habel +## Copyright (C) 2009 Jaroslav Hajek <highegg@gmail.com> ## ## This file is part of Octave. ## @@ -138,7 +139,9 @@ 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 +164,18 @@ 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 +187,25 @@ 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)); + ## if XI, YI are vectors, X and Y should share their orientation. + if (rows (XI) == 1) + if (rows (X) != 1) + X = X.'; + endif + if (rows (Y) != 1) + Y = Y.'; + endif + elseif (columns (XI) == 1) + if (columns (X) != 1) + X = X.'; + endif + if (columns (Y) != 1) + Y = Y.'; + endif + endif - xidx = lookup (X(1, :), XI, "lr"); - yidx = lookup (Y(:, 1), YI, "lr"); + xidx = lookup (X, XI, "lr"); + yidx = lookup (Y, YI, "lr"); if (strcmp (method, "linear")) ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy @@ -202,41 +221,98 @@ 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 - Y(yidx) > Y(yidx + 1) - YI); idx = sub2ind (size (Z), yidx+jj, xidx+ii); ZI = Z(idx); + + 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 - ## 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; + 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 - 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; + 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 - ZI = reshape (ZI, shape); else ## If X and Y vectors produce a grid from them @@ -285,6 +361,15 @@ %! 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,6 +379,33 @@ %! 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); @@ -303,6 +415,15 @@ %! 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); @@ -311,6 +432,15 @@ %! [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];
--- a/src/ChangeLog Mon Feb 09 16:53:12 2009 -0500 +++ b/src/ChangeLog Fri Mar 14 12:54:56 2008 +0100 @@ -1,3 +1,9 @@ +2009-02-10 Jaroslav Hajek <highegg@gmail.com> + + * DLD-FUNCTIONS/__pchip_deriv__.cc (F__pchip_deriv__): + Add support for computing pchip derivatives along rows of matrix. + Eliminate redundant copying by using const args where appropriate. + 2009-02-09 John W. Eaton <jwe@octave.org> * version.h (OCTAVE_VERSION): Now 3.1.52.
--- a/src/DLD-FUNCTIONS/__pchip_deriv__.cc Mon Feb 09 16:53:12 2009 -0500 +++ 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, 2009 Jaroslav Hajek This file is part of Octave. @@ -20,6 +21,9 @@ */ +// Jaroslav Hajek, Feb 2008: handle row-wise derivatives, +// use const pointers to avoid unnecessary copying. + #ifdef HAVE_CONFIG_H #include <config.h> #endif @@ -34,12 +38,12 @@ 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); F77_RET_T - F77_FUNC (pchim, PCHIM) (const octave_idx_type& n, float *x, float *f, + F77_FUNC (pchim, PCHIM) (const octave_idx_type& n, const float *x, const float *f, float *d, const octave_idx_type &incfd, octave_idx_type *ierr); } @@ -49,14 +53,16 @@ DEFUN_DLD (__pchip_deriv__, args, , "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {} __pchip_deriv__ (@var{x}, @var{y})\n\ +@deftypefn {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) { if (args(0).is_single_type () || args(1).is_single_type ()) { @@ -67,34 +73,33 @@ 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; } - FloatColumnVector dvec (nx), yvec (nx); + const float *yvec = ymat.data (); FloatMatrix dmat (nyr, nyc); + float *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 (octave_idx_type i = (rows ? nyr : nyc); i > 0; i--) { - for (int r = 0; r < nx; r++) - yvec(r) = ymat(r,c); + F77_FUNC (pchim, PCHIM) (nx, xvec.data (), + yvec, dvec, incfd, &ierr); - F77_FUNC (pchim, PCHIM) (nx, xvec.fortran_vec (), - yvec.fortran_vec (), - dvec.fortran_vec (), incfd, &ierr); + yvec += inc; + dvec += inc; if (ierr < 0) { error ("PCHIM error: %i\n", ierr); return retval; } - - for (int r=0; r<nx; r++) - dmat(r,c) = dvec(r); } retval = dmat; @@ -108,34 +113,33 @@ 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 (octave_idx_type 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.data (), + yvec, dvec, incfd, &ierr); - F77_FUNC (dpchim, DPCHIM) (nx, xvec.fortran_vec (), - yvec.fortran_vec (), - dvec.fortran_vec (), incfd, &ierr); + yvec += inc; + dvec += inc; if (ierr < 0) { error ("DPCHIM error: %i\n", ierr); return retval; } - - for (int r=0; r<nx; r++) - dmat(r,c) = dvec(r); } retval = dmat;