Mercurial > hg > octave-avbm
changeset 6702:b2391d403ed2
[project @ 2007-06-12 21:39:26 by dbateman]
author | dbateman |
---|---|
date | Tue, 12 Jun 2007 21:39:27 +0000 |
parents | 3933e0693fe0 |
children | 31c8d115f25d |
files | doc/ChangeLog doc/interpreter/hashing.txi doc/interpreter/interp.txi doc/interpreter/octave.texi doc/interpreter/system.txi scripts/ChangeLog scripts/general/__splinen__.m scripts/general/bicubic.m scripts/general/interp1.m scripts/general/interp2.m scripts/general/interp3.m scripts/general/interpft.m scripts/general/interpn.m scripts/polynomial/spline.m src/ChangeLog src/DLD-FUNCTIONS/__lin_interpn__.cc src/DLD-FUNCTIONS/interpn.cc src/Makefile.in |
diffstat | 18 files changed, 925 insertions(+), 364 deletions(-) [+] |
line wrap: on
line diff
--- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,3 +1,16 @@ +2007-06-12 David Bateman <dbateman@free.fr> + + * interpreter/interp.txi: Split into two section and document + interp3 and the differences in the treatement of the dimensions + between interpn and interp3. + * hashing.txi: Remove. + * system.txi: Move it here as a subsection. Include explanation + and example. + * interpreter/octave.texi: Add sections for the Interpolation + chapter. Remove references to Hashing chapter and hashing.texi, + and subsections for hashing to system utilities chapter. + + 2007-06-12 2007-06-10 Søren Hauberg <hauberg@gmail.com> * interpreter/diffeq.txi: Note that x-dot is the derivative of x.
--- a/doc/interpreter/hashing.txi +++ b/doc/interpreter/hashing.txi @@ -1,8 +0,0 @@ -@c Copyright (C) 2007 John W. Eaton -@c This is part of the Octave manual. -@c For copying conditions, see the file gpl.texi. - -@node Hashing Functions -@chapter Hashing Functions - -@DOCSTRING(md5sum)
--- a/doc/interpreter/interp.txi +++ b/doc/interpreter/interp.txi @@ -5,18 +5,96 @@ @node Interpolation @chapter Interpolation -Octave provides the following functions for interpolation. +@menu +* One-dimensional Interpolation:: +* Multi-dimensional Interpolation:: +@end menu + +@node One-dimensional Interpolation +@section One-dimensional Interpolation @DOCSTRING(interp1) -@DOCSTRING(interp2) +Fourier interpolation, is a resampling technique where a signal is +converted to the frequency domain, padded with zeros and then +reconverted to the time domain. @DOCSTRING(interpft) -@DOCSTRING(interpn) +There are two significant limitations on Fourier interpolation. Firstly, +the function signal is assumed to be periodic, and so no periodic +signals will be poorly represented at the edges. Secondly, both the +signal and its interpolation are required to be sampled at equispaced +points. An example of the use of @code{interpft} is -@DOCSTRING(bicubic) +@example +@group +t = 0 : 0.3 : pi; dt = t(2)-t(1); +n = length (t); k = 100; +ti = t(1) + [0 : k-1]*dt*n/k; +y = sin (4*t + 0.3) .* cos (3*t - 0.1); +yp = sin (4*ti + 0.3) .* cos (3*ti - 0.1); +plot (ti, yp, 'g', ti, interp1(t, y, ti, 'spline'), 'b', ... + ti, interpft (y, k), 'c', t, y, 'r+'); +legend ('sin(4t+0.3)cos(3t-0.1','spline','interpft','data'); +@end group +@end example + +which demonstrates the poor behavior of Fourier interpolation for non +periodic functions. + +In additional the support function @code{spline} and @code{lookup} that +underlie the @code{interp1} function can be called directly. @DOCSTRING(spline) +The @code{lookup} is used by other interpolation function to identify +the points of the original data that are closest to the current point +of interest. + @DOCSTRING(lookup) + +@node Multi-dimensional Interpolation +@section Multi-dimensional Interpolation + +There are three multi-dimensional interpolation function in Octave, with +similar capabilities. + +@DOCSTRING(interp2) + +@DOCSTRING(interp3) + +@DOCSTRING(interpn) + +A significant difference between @code{interpn} and the other two +multidimensional interpolation function is the fashion in which the +dimensions are treated. For @code{interp2} and @code{interp3}, the 'y' +axis is considered to be the columns of the matrix, whereas the 'x' +axis corresponds to the rows the the array. As Octave indexes arrays in +column major order, the first dimension of any array is the columns, and +so @code{interpn} effectively reverses the 'x' and 'y' dimensions. +Consider the example + +@example +@group +x = y = z = -1:1; +f = @@(x,y,z) x.^2 - y - z.^2; +[xx, yy, zz] = meshgrid (x, y, z); +v = f (xx,yy,zz); +xi = yi = zi = -1:0.5:1; +[xxi, yyi, zzi] = meshgrid (xi, yi, zi); +vi = interp3(x, y, z, v, xxi, yyi, zzi); +[xxi, yyi, zzi] = ndgrid (xi, yi, zi); +vi2 = interpn(x, y, z, v, xxi, yyi, zzi); +@end group +@end example + +@noindent +where @code{vi} and @code{vi2} are identical. The reversal of the +dimensions is treated in the @code{meshgrid} and @code{ndgrid} functions +respectively. + +In additional the support function @code{bicubic} that underlies the +cubic interpolation of @code{interp2} function can be called directly. + +@DOCSTRING(bicubic)
--- a/doc/interpreter/octave.texi +++ b/doc/interpreter/octave.texi @@ -158,7 +158,6 @@ * Polynomial Manipulations:: * Interpolation:: * Geometry:: -* Hashing Functions:: * Control Theory:: * Signal Processing:: * Image Processing:: @@ -437,7 +436,9 @@ * Models:: * Distributions:: -Hashing Functions +Interpolation +* One-dimensional Interpolation:: +* Multi-dimensional Interpolation:: Control Theory @@ -488,6 +489,7 @@ * Password Database Functions:: * Group Database Functions:: * System Information:: +* Hashing Functions:: Packages @@ -587,7 +589,6 @@ @include poly.texi @include interp.texi @include geometry.texi -@include hashing.texi @include control.texi @include signal.texi @include image.texi
--- a/doc/interpreter/system.txi +++ b/doc/interpreter/system.txi @@ -23,6 +23,7 @@ * Password Database Functions:: * Group Database Functions:: * System Information:: +* Hashing Functions:: @end menu @node Timing Utilities @@ -398,3 +399,34 @@ @DOCSTRING(octave_config_info) @DOCSTRING(getrusage) + +@node Hashing Functions +@section Hashing Functions + +It is often necessary to find if two strings or files are +identical. This might be done by comparing them character by character +and looking for differences. However, this can be slow, and so comparing +a hash of the string or file can be a rapid way of finding if the files +differ. + +Another use of the hashing function is to check for file integrity. The +user can check the hash of the file against a known value and find if +the file they have is the same as the one that the original hash was +produced with. + +Octave supplies the @code{md5sum} function to perfrom MD5 hashes on +strings and files. An example of the use of @code{md5sum} function might +be + +@example +@group +if exist (file, "file") + hash = md5sum (file); +else + # Treat the variable "file" as a string + hash = md5sum (file, true); +endif +@end group +@end example + +@DOCSTRING(md5sum)
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,21 @@ +2007-06-12 David Bateman <dbateman@free.fr> + + * general/interp1.m: Change examples to use new graphics + interface. + * general/__splinen__.m: New support function for N-dimensional + spline interpolation. + * general/bicubic.m: Allow definition of extrapolation + value. Adapt tests to use new graphics interface + * general/interp2.m: Call __splinen__ for 2-D spline + interpolation. Make the lookup table code only be called for + linear and nearest methods. + * general/interpn.m: New function for N-dimensional, linear, nearest + and spline interpolation. + * general/interp3.m: New function for 3-dimensional, linear, nearest + and spline interpolation. + * polynomial/spline.m: Change examples to use new graphics + interface. + 2007-06-12 Steve M. Robbins <steve@sumost.ca> * statistics/tests/wilcoxon_test.m: Error if N <= 25.
new file mode 100644 --- /dev/null +++ b/scripts/general/__splinen__.m @@ -0,0 +1,47 @@ +## Copyright (C) 2007 David Bateman +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{yi} = } __splinen__ (@var{x}, @var{y}, @var{xi}) +## Internal support function for multi-dimensional splines. +## @end deftypefn + +## FIXME: Allow arbitrary grids.. + +function yi = __splinen__ (x, y, xi, extrapval, f) + if (nargin != 5) + error ("Incorrect number of arguments"); + endif + + if (!iscell (x) || length(x) < ndims(y) || any (! cellfun (@isvector, x)) || + !iscell (xi) || length(xi) < ndims(y) || any (! cellfun (@isvector, xi))) + error ("%s: non gridded data or dimensions inconsistent", f); + endif + yi = y; + for i = length(x):-1:1 + yi = spline (x{i}, yi, xi{i}).'; + endfor + + [xi{:}] = ndgrid (xi{:}); + 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; +endfunction
--- a/scripts/general/bicubic.m +++ b/scripts/general/bicubic.m @@ -18,32 +18,37 @@ ## 02110-1301, USA. ## -*- texinfo -*- -## @deftypefn {Function File} {@var{zi}=} bicubic (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}) +## @deftypefn {Function File} {@var{zi}=} bicubic (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{extrapval}) ## ## Return a matrix @var{zi} corresponding to the bicubic ## interpolations at @var{xi} and @var{yi} of the data supplied -## as @var{x}, @var{y} and @var{z}. +## as @var{x}, @var{y} and @var{z}. Points outside the grid are set +## to @var{extrapval} ## -## For further information please see bicubic.pdf available at -## @url{http://wiki.woodpecker.org.cn/moin/Octave/Bicubic} +## See @url{http://wiki.woodpecker.org.cn/moin/Octave/Bicubic} +## for further information. ## @seealso{interp2} ## @end deftypefn ## Bicubic interpolation method. ## Author: Hoxide Ma <hoxide_dirac@yahoo.com.cn> -function F = bicubic (X, Y, Z, XI, YI, spline_alpha) +function F = bicubic (X, Y, Z, XI, YI, extrapval, spline_alpha) - if (nargin < 1 || nargin > 6) + if (nargin < 1 || nargin > 7) print_usage (); endif - if (nargin == 6 && prod (size (spline_alpha)) == 1) + if (nargin == 7 && isscalar(spline_alpha)) a = spline_alpha else a = 0.5; endif + if (nargin < 6) + extrapval = NaN; + endif + if (nargin <= 2) ## bicubic (Z) or bicubic (Z, 2) if (nargin == 1) @@ -177,12 +182,12 @@ p(int,inds+2) .* cs2 + p(int,inds+3) .* cs3); endfor - ## set points outside the table to NaN + ## set points outside the table to extrapval if (! (isempty (xfirst_ind) && isempty (xlast_ind))) - F(:, [xfirst_ind, xlast_ind]) = NaN; + F(:, [xfirst_ind, xlast_ind]) = extrapval; endif if (! (isempty (yfirst_ind) && isempty (ylast_ind))) - F([yfirst_ind; ylast_ind], :) = NaN; + F([yfirst_ind; ylast_ind], :) = extrapval; endif endfunction @@ -191,7 +196,7 @@ %! A=[13,-1,12;5,4,3;1,6,2]; %! x=[0,1,4]+10; y=[-10,-9,-8]; %! xi=linspace(min(x),max(x),17); -%! yi=linspace(min(y),max(y),26); +%! yi=linspace(min(y),max(y),26)'; %! mesh(xi,yi,bicubic(x,y,A,xi,yi)); %! [x,y] = meshgrid(x,y); %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
--- a/scripts/general/interp1.m +++ b/scripts/general/interp1.m @@ -70,8 +70,9 @@ ## spl=interp1(xp,yp,xf,'spline'); ## cub=interp1(xp,yp,xf,'cubic'); ## near=interp1(xp,yp,xf,'nearest'); -## plot(xf,yf,';original;',xf,lin,';linear;',xf,spl,';spline;',... -## xf,cub,';cubic;',xf,near,';nearest;',xp,yp,'*;;'); +## plot(xf,yf,"r",xf,lin,"g",xf,spl,"b", ... +## xf,cub,"c",xf,near,"m",xp,yp,"r*"); +## legend ("original","linear","spline","cubic","nearest") ## @end group ## @end example ## @@ -327,8 +328,8 @@ %! spl=interp1(xp,yp,xf,"spline"); %! cub=interp1(xp,yp,xf,"pchip"); %! near=interp1(xp,yp,xf,"nearest"); -%! plot(xf,yf,";original;",xf,near,";nearest;",xf,lin,";linear;",... -%! xf,cub,";pchip;",xf,spl,";spline;",xp,yp,"*;;"); +%! plot(xf,yf,"r",xf,near,"g",xf,lin,"b",xf,cub,"c",xf,spl,"m",xp,yp,"r*"); +%! legend ("original","nearest","linear","pchip","spline") %! %-------------------------------------------------------- %! % confirm that interpolated function matches the original @@ -339,8 +340,8 @@ %! spl=interp1(xp,yp,xf,"*spline"); %! cub=interp1(xp,yp,xf,"*cubic"); %! near=interp1(xp,yp,xf,"*nearest"); -%! plot(xf,yf,";*original;",xf,near,";*nearest;",xf,lin,";*linear;",... -%! xf,cub,";*cubic;",xf,spl,";*spline;",xp,yp,"*;;"); +%! plot(xf,yf,"r",xf,near,"g",xf,lin,"b",xf,cub,"c",xf,spl,"m",xp,yp,"r*"); +%! legend ("*original","*nearest","*linear","*cubic","*spline") %! %-------------------------------------------------------- %! % confirm that interpolated function matches the original
--- a/scripts/general/interp2.m +++ b/scripts/general/interp2.m @@ -62,7 +62,7 @@ ## Cubic interpolation from four nearest neighbors. ## @item 'spline' ## Cubic spline interpolation--smooth first and second derivatives -## throughout the curve (not implemented yet). +## throughout the curve. ## @end table ## ## If a scalar value @var{extrapval} is defined as the final value, then @@ -161,75 +161,102 @@ error ("interp2 expected numeric XI, YI"); endif - ## 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)) - error ("X and Y must be matrices of same size"); - endif - if (! size_equal (X, Z)) - error ("X and Y size must match Z dimensions"); - endif + + if (strcmp (method, "linear") || strcmp (method, "nearest")) + + ## 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)) + error ("X and Y must be matrices of same size"); + endif + if (! size_equal (X, Z)) + error ("X and Y size must match Z dimensions"); + endif + + ## If Xi and Yi are vectors of different orientation build a grid + if ((rows (XI) == 1 && columns (YI) == 1) + || (columns (XI) == 1 && rows (YI) == 1)) + [XI, YI] = meshgrid (XI, YI); + elseif (! size_equal (XI, YI)) + error ("XI and YI must be matrices of same size"); + endif - ## If Xi and Yi are vectors of different orientation build a grid - if ((rows (XI) == 1 && columns (YI) == 1) - || (columns (XI) == 1 && rows (YI) == 1)) - [XI, YI] = meshgrid (XI, YI); - elseif (! size_equal (XI, YI)) - 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; - shape = size (XI); - XI = reshape (XI, 1, prod (shape)); - YI = reshape (YI, 1, prod (shape)); + if (strcmp (method, "linear")) + ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy + ## + ## 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; + d = Z(2:zr, 2:zc) - a - b - c; - xidx = lookup (X(1, 2:end-1), XI) + 1; - yidx = lookup (Y(2:end-1, 1), YI) + 1; + 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))'; + + ## apply plane equation + ZI = a(idx) + b(idx).*Xsc + c(idx).*Ysc + d(idx).*Xsc.*Ysc; - if (strcmp (method, "linear")) - ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy - ## - ## 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; - d = Z(2:zr, 2:zc) - a - b - c; + 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); + idx = sub2ind (size (Z), yidx+jj, xidx+ii); + ZI = Z(idx); + endif - 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))'; - - ## apply plane equation - ZI = a(idx) + b(idx).*Xsc + c(idx).*Ysc + d(idx).*Xsc.*Ysc; + ## set points outside the table to NaN + ZI (XI < X(1,1) | XI > X(1,end) | YI < Y(1,1) | YI > Y(end,1)) = extrapval; + ZI = reshape (ZI, shape); + else - 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); - idx = sub2ind (size (Z), yidx+jj, xidx+ii); - ZI = Z(idx); - - elseif (strcmp (method, "cubic")) - ## FIXME bicubic doesn't handle arbitrary XI, YI - ZI = bicubic (X, Y, Z, XI(1,:), YI(:,1)); + ## If X and Y vectors produce a grid from them + if (isvector (X) && isvector (Y)) + X = X(:).'; + Y = Y(:); + if (!isequal ([length(X), length(Y)], size(Z))) + error ("X and Y size must match Z dimensions"); + endif + elseif (!size_equal (X, Y)) + error ("X and Y must be matrices of same size"); + if (! size_equal (X, Z)) + error ("X and Y size must match Z dimensions"); + endif + endif - elseif (strcmp (method, "spline")) - ## FIXME Implement 2-D (or in fact ND) spline interpolation - error ("interp2, spline interpolation not yet implemented"); + ## If Xi and Yi are vectors of different orientation build a grid + if ((rows (XI) == 1 && columns (YI) == 1) + || (columns (XI) == 1 && rows (YI) == 1)) + ## Do nothing + elseif (! size_equal (XI, YI)) + error ("XI and YI must be matrices of same size"); + endif - else - error ("interpolation method not recognized"); - endif + ## FIXME bicubic/__splinen__ don't handle arbitrary XI, YI + if (strcmp (method, "cubic")) + ZI = bicubic (X, Y, Z, XI(1,:), YI(:,1), extrapval); - ## set points outside the table to NaN - ZI (XI < X(1,1) | XI > X(1,end) | YI < Y(1,1) | YI > Y(end,1)) = extrapval; - ZI = reshape (ZI, shape); + elseif (strcmp (method, "spline")) + ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrapval, + "spline"); + else + error ("interpolation method not recognized"); + endif + endif endfunction %!demo @@ -250,15 +277,24 @@ %! [x,y] = meshgrid(x,y); %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; -%!#demo +%!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); +%! 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 +%! 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; + %!test % simple test %! x = [1,2,3]; %! y = [4,5,6,7];
new file mode 100644 --- /dev/null +++ b/scripts/general/interp3.m @@ -0,0 +1,114 @@ +## Copyright (C) 2007 David Bateman +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{vi} =} interp3 (@var{x}, @var{y},@var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi}) +## @deftypefnx {Function File} {@var{vi} =} interp3 (@var{v}, @var{xi}, @var{yi}, @var{zi}) +## @deftypefnx {Function File} {@var{vi} =} interp3 (@var{v}, @var{m}) +## @deftypefnx {Function File} {@var{vi} =} interp3 (@var{v}) +## @deftypefnx {Function File} {@var{vi} =} interp3 (@dots{}, @var{method}) +## @deftypefnx {Function File} {@var{vi} =} interp3 (@dots{}, @var{method}, @var{extrapval}) +## +## Perform 3-dimensional interpolation. Each element of then 3-dimensional +## array @var{v} represents a value at a location given by the parameters +## @var{x}, @var{y}, and @var{z}. The parameters @var{x}, @var{x}, and +## @var{z} are either 3-dimensional arrays of the same size as the array +## @var{v} in the 'meshgrid' format or vectors. The parameters @var{xi}, etc +## respect a similar format to @var{x}, etc, and they represent the points +## at which the array @var{vi} is interpolated. +## +## If @var{x}, @var{y}, @var{z} are ommitted, they are assumed to be +## @code{x = 1 : size (@var{v}, 2)}, @code{y = 1 : size (@var{v}, 1)} and +## @code{z = 1 : size (@var{v}, 3)}. If @var{m} is specified, then +## the interpolation adds a point half way between each of the interplation +## 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: +## +## @table @asis +## @item 'nearest' +## Return the nearest neighbour. +## @item 'linear' +## Linear interpolation from nearest neighbours. +## @item 'cubic' +## Cubic interpolation from four nearest neighbours (not implemented yet). +## @item 'spline' +## Cubic spline interpolation--smooth first and second derivatives +## throughout the curve. +## @end table +## +## The default method is 'linear'. +## +## If @var{extrap} is the string 'extrap', then extrapolate values beyond +## the endpoints. If @var{extrap} is a number, replace values beyond the +## endpoints with that number. If @var{extrap} is missing, assume NaN. +## @seealso{interp1, interp2, spline, meshgrid} +## @end deftypefn + +function vi = interp3 (varargin) + method = "linear"; + extrapval = NaN; + nargs = nargin; + + if (nargin < 1) + print_usage (); + endif + + if (ischar (varargin {end})) + method = varargin {end}; + nargs = nargs - 1; + elseif (ischar (varargin {end - 1})) + if (! isnumeric (vargin {end}) || ! isscalar (vargin {end})) + error ("extrapal is expected to be a numeric scalar"); + endif + method = varargin {end - 1}; + nargs = nargs - 2; + endif + + if (nargs < 3 || (nargs == 4 && ! isvector (varargin {1}) && + nargs == (ndims (varargin {1}) + 1))) + v = varargin {1}; + if (ndims (v) != 3) + error ("expect 3-dimensional array of values"); + endif + vi = ipermute (interpn (permute(varargin, [1, 3, 2, 4]){:}), [2, 1, 3]); + elseif (nargs == 7 && nargs == (2 * ndims (varargin {ceil (nargs / 2)})) + 1) + v = varargin {4}; + if (ndims (v) != 3) + error ("expect 3-dimensional array of values"); + endif + vi = ipermute (interpn (permute(varargin, [2, 1, 3, 4, 6, 5, 7]){:}), + [2, 1, 3]); + else + error ("wrong number or incorrectly formatted input arguments"); + endif +endfunction + +%!test +%! x = y = z = -1:1; +%! f = @(x,y,z) x.^2 - y - z.^2; +%! [xx, yy, zz] = meshgrid (x, y, z); +%! v = f (xx,yy,zz); +%! xi = yi = zi = -1:0.5:1; +%! [xxi, yyi, zzi] = meshgrid (xi, yi, zi); +%! vi = interp3(x, y, z, v, xxi, yyi, zzi); +%! [xxi, yyi, zzi] = ndgrid (xi, yi, zi); +%! vi2 = interpn(x, y, z, v, xxi, yyi, zzi); +%! assert (vi, vi2);
--- a/scripts/general/interpft.m +++ b/scripts/general/interpft.m @@ -101,10 +101,9 @@ %! ti = t(1) + [0 : k-1]*dt*n/k; %! y = sin (4*t + 0.3) .* cos (3*t - 0.1); %! yp = sin (4*ti + 0.3) .* cos (3*ti - 0.1); -%! plot (ti, yp, 'g;sin(4t+0.3)cos(3t-0.1);', ... -%! ti, interp1 (t, y, ti, 'spline'), 'b;spline;', ... -%! ti, interpft (y ,k), 'c;interpft;', ... -%! t, y, 'r+;data;'); +%! plot (ti, yp, 'g', ti, interp1(t, y, ti, 'spline'), 'b', ... +%! ti, interpft (y, k), 'c', t, y, 'r+'); +%! legend ('sin(4t+0.3)cos(3t-0.1','spline','interpft','data'); %!shared n,y %! x = [0:10]'; y = sin(x); n = length (x);
new file mode 100644 --- /dev/null +++ b/scripts/general/interpn.m @@ -0,0 +1,218 @@ +## Copyright (C) 2007 David Bateman +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{vi} =} interpn (@var{x1}, @var{x2}, @dots{}, @var{v}, @var{y1}, @var{y2}, @dots{}) +## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}, @var{y1}, @var{y2}, @dots{}) +## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}, @var{m}) +## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}) +## @deftypefnx {Function File} {@var{vi} =} interpn (@dots{}, @var{method}) +## @deftypefnx {Function File} {@var{vi} =} interpn (@dots{}, @var{method}, @var{extrapval}) +## +## Perform @var{n}-dimensional interpolation, where @var{n} is at least two. +## Each element of then @var{n}-dimensional array @var{v} represents a value +## at a location given by the parameters @var{x1}, @var{x2}, @dots{}, @var{xn}. +## The parameters @var{x1}, @var{x2}, @dots{}, @var{xn} are either +## @var{n}-dimensional arrays of the same size as the array @var{v} in +## the 'ndgrid' format or vectors. The parameters @var{y1}, etc respect a +## similar format to @var{x1}, etc, and they represent the points at which +## the array @var{vi} is interpolated. +## +## If @var{x1}, @dots{}, @var{xn} are ommitted, they are assumed to be +## @code{x1 = 1 : size (@var{v}, 1)}, etc. If @var{m} is specified, then +## the interpolation adds a point half way between each of the interplation +## 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: +## +## @table @asis +## @item 'nearest' +## Return the nearest neighbour. +## @item 'linear' +## Linear interpolation from nearest neighbours. +## @item 'cubic' +## Cubic interpolation from four nearest neighbours (not implemented yet). +## @item 'spline' +## Cubic spline interpolation--smooth first and second derivatives +## throughout the curve. +## @end table +## +## The default method is 'linear'. +## +## If @var{extrap} is the string 'extrap', then extrapolate values beyond +## the endpoints. If @var{extrap} is a number, replace values beyond the +## endpoints with that number. If @var{extrap} is missing, assume NaN. +## @seealso{interp1, interp2, spline, ndgrid} +## @end deftypefn + +function vi = interpn (varargin) + + method = "linear"; + extrapval = NaN; + nargs = nargin; + + if (nargin < 1) + print_usage (); + endif + + if (ischar (varargin {end})) + method = varargin {end}; + nargs = nargs - 1; + elseif (ischar (varargin {end - 1})) + if (! isnumeric (vargin {end}) || ! isscalar (vargin {end})) + error ("extrapal is expected to be a numeric scalar"); + endif + method = varargin {end - 1}; + nargs = nargs - 2; + endif + + if (nargs < 3) + v = varargin {1}; + m = 1; + if (nargs == 2) + m = varargin {2}; + if (! isnumeric (m) || ! isscalar (m) || floor (m) != m) + error ("m is expected to be a integer scalar"); + endif + endif + sz = size (v); + nd = ndims (v); + x = cell (1, nd); + y = cell (1, nd); + for i = 1 : nd; + x{i} = 1 : sz(i); + y{i} = 1 : (1 / (2 ^ m)) : sz(i); + endfor + elseif (! isvector (varargin {1}) && nargs == (ndims (varargin {1}) + 1)) + v = varargin {1}; + sz = size (v); + nd = ndims (v); + x = cell (1, nd); + y = varargin (2 : nargs); + for i = 1 : nd; + x{i} = 1 : sz(i); + endfor + elseif (rem (nargs, 2) == 1 && nargs == + (2 * ndims (varargin {ceil (nargs / 2)})) + 1) + nv = ceil (nargs / 2); + v = varargin {nv}; + sz = size (v); + nd = ndims (v); + x = varargin (1 : (nv - 1)); + y = varargin ((nv + 1) : nargs); + else + error ("wrong number or incorrectly formatted input arguments"); + endif + + if (any (! cellfun (@isvector, x))) + for i = 2 : nd + if (! size_equal (x{1}, x{i}) || ! size_equal (x{i}, v)) + error ("dimensional mismatch"); + endif + idx (1 : nd) = {1}; + idx (i) = ":"; + x{i} = x{i}(idx{:}); + endfor + idx (1 : nd) = {1}; + idx (1) = ":"; + x{1} = x{1}(idx{:}); + endif + + if (strcmp (method, "linear") || strcmp (method, "nearest")) + if (all (cellfun (@isvector, y))) + [y{:}] = ndgrid (y{:}); + endif + elseif (any (! cellfun (@isvector, x))) + for i = 1 : nd + idx (1 : nd) = {1}; + idx (i) = ":"; + y{i} = y{i}(idx{:}); + endfor + endif + + method = tolower (method); + if (strcmp (method, "linear")) + vi = __lin_interpn__ (x{:}, v, y{:}); + vi (vi == NaN) = extrapval; + elseif (strcmp (method, "nearest")) + yshape = size (y{1}); + yidx = cell (1, nd); + for i = 1 : nd + y{i} = y{i}(:); + yidx{i} = lookup (x{i}(2:end-1), y{i}) + 1; + endfor + idx = cell (1,nd); + for i = 1 : nd + idx {i} = yidx{i} + (y{i} - x{i}(yidx{i}).' > ... + x{i}(yidx{i} + 1).' - y{i}); + endfor + vi = v (sub2ind (sz, idx{:})); + idx = zeros (prod(yshape),1); + for i = 1 : nd + idx |= y{i} < min (x{i}(:)) | y{i} > max (x{i}(:)); + endfor + vi(idx) = extrapval; + vi = reshape (vi, yshape); + elseif (strcmp (method, "spline")) + vi = __splinen__ (x, v, y, extrapval, "interpn"); + elseif (strcmp (method, "cubic")) + error ("cubic interpolation not yet implemented"); + else + error ("unrecognized interpolation method"); + endif + +endfunction + +%!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,interpn(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,interpn(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,interpn(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,interpn(x,y,A.',xi,yi,"spline").'); +%! [x,y] = meshgrid(x,y); +%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; +
--- a/scripts/polynomial/spline.m +++ b/scripts/polynomial/spline.m @@ -212,9 +212,8 @@ %! x = 0:10; y = sin(x); %! xspline = 0:0.1:10; yspline = spline(x,y,xspline); %! title("spline fit to points from sin(x)"); -%! plot(xspline,sin(xspline),";original;",... -%! xspline,yspline,"-;interpolation;",... -%! x,y,"+;interpolation points;"); +%! plot(xspline,sin(xspline),"r",xspline,yspline,"g-",x,y,"b+"); +%! legend("original","interpolation","interpolation points"); %! %-------------------------------------------------------- %! % confirm that interpolated function matches the original
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,10 @@ +2007-06-12 David Bateman <dbateman@free.fr> + + * DLD-FUNCTIONS/interpn.cc: Remove it. + * DLD-FUNCTIONS/__lin_interpn__.cc: Move it. This is now a support + function of interpn.m. + * Makefile.in (DLD_XSRC): Remove interpn.cc and add __lin_interpn__.cc. + 2007-06-07 David Bateman <dbateman@free.fr> * ov-fcn-handles.cc (octave_fcn_handle::save_hdf5): More care that
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/__lin_interpn__.cc @@ -0,0 +1,261 @@ +/* + +Copyright (C) 2007 Alexander Barth + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +02110-1301, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "dNDArray.h" + +#include "defun-dld.h" +#include "error.h" +#include "oct-obj.h" + +// equivalent to isvector.m + +bool +isvector (const NDArray& array) +{ + const dim_vector dv = array.dims (); + return dv.length () == 2 && (dv(0) == 1 || dv(1) == 1); +} + +// lookup a value in a sorted table (lookup.m) +octave_idx_type +lookup (const double *x, octave_idx_type n, double y) +{ + octave_idx_type j, j0, j1; + + if (y > x[n-1] || y < x[0]) + return -1; + +#ifdef EXHAUSTIF + for (j = 0; j < n - 1; j++) + { + if (x[j] <= y && y <= x[j+1]) + return j; + } +#else + j0 = 0; + j1 = n - 1; + + while (true) + { + j = (j0+j1)/2; + + if (y <= x[j+1]) + { + if (x[j] <= y) + return j; + + j1 = j; + } + + if (x[j] <= y) + j0 = j; + } + +#endif +} + +// n-dimensional linear interpolation + +void +lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale, + octave_idx_type Ni, double extrapval, const double **x, + const double *v, const double **y, double *vi) +{ + bool out = false; + int bit; + + OCTAVE_LOCAL_BUFFER (double, coef, 2*n); + OCTAVE_LOCAL_BUFFER (octave_idx_type, index, n); + + // loop over all points + for (octave_idx_type m = 0; m < Ni; m++) + { + // loop over all dimensions + for (int i = 0; i < n; i++) + { + index[i] = lookup (x[i], size[i], y[i][m]); + out = index[i] == -1; + + if (out) + break; + else + { + octave_idx_type j = index[i]; + coef[2*i+1] = (y[i][m] - x[i][j])/(x[i][j+1] - x[i][j]); + coef[2*i] = 1 - coef[2*i+1]; + } + } + + + if (out) + vi[m] = extrapval; + else + { + vi[m] = 0; + + // loop over all corners of hypercube (1<<n = 2^n) + for (int i = 0; i < (1 << n); i++) + { + double c = 1; + octave_idx_type l = 0; + + // loop over all dimensions + for (int j = 0; j < n; j++) + { + // test if the jth bit in i is set + bit = i >> j & 1; + l += scale[j] * (index[j] + bit); + c *= coef[2*j+bit]; + } + + vi[m] += c * v[l]; + } + } + } +} + +DEFUN_DLD (__lin_interpn__, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{vi} =} __lin_interpn__ (@var{x1}, @var{x2}, @dots{}, @var{xn}, @var{v}, @var{y1}, @var{y2}, @dots{}, @var{yn})\n\ +Perform @var{n}-dimensional interpolation. Each element of then\n\ +@var{n}-dimensional array @var{v} represents a value at a location\n\ +given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters\n\ +@var{x1}, @var{x2}, @dots{}, @var{xn} are either @var{n}-dimensional\n\ +arrays of the same size as the array @var{v} in the \"ndgrid\" format\n\ +or vectors. The parameters @var{y1}, @var{y2}, @dots{}, @var{yn} are\n\ +all @var{n}-dimensional arrays of the same size and represent the\n\ +points at which the array @var{vi} is interpolated.\n\ +\n\ +This function only performs linear interpolation.\n\ +@seealso{interp1, interp2, ndgrid}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin < 2 || nargin % 2 == 0) + { + print_usage (); + return retval; + } + + // dimension of the problem + int n = (nargin-1)/2; + + OCTAVE_LOCAL_BUFFER (NDArray, X, n); + OCTAVE_LOCAL_BUFFER (NDArray, Y, n); + + OCTAVE_LOCAL_BUFFER (const double *, x, n); + OCTAVE_LOCAL_BUFFER (const double *, y, n); + OCTAVE_LOCAL_BUFFER (octave_idx_type, scale, n); + OCTAVE_LOCAL_BUFFER (octave_idx_type, size, n); + + const NDArray V = args(n).array_value (); + NDArray Vi = NDArray (args(n+1).dims ()); + + if (error_state) + { + print_usage (); + return retval; + } + + const double *v = V.data (); + double *vi = Vi.fortran_vec (); + octave_idx_type Ni = Vi.numel (); + + double extrapval = octave_NaN; + + for (int i = 0; i < n; i++) + { + X[i] = args(i).array_value (); + Y[i] = args(n+i+1).array_value (); + + if (error_state) + { + print_usage (); + return retval; + } + + y[i] = Y[i].data (); + size[i] = V.dims()(i); + + if (Y[0].dims () != Y[i].dims ()) + { + error ("interpn: incompatible size of argument number %d", n+i+2); + return retval; + } + } + + // offset in memory of each dimension + + scale[0] = 1; + + for (int i = 1; i < n; i++) + scale[i] = scale[i-1] * size[i-1]; + + // tests if X[0] is a vector, if yes, assume that all elements of X are + // in the ndgrid format. + + if (! isvector (X[0])) + { + for (int i = 0; i < n; i++) + { + if (X[i].dims () != V.dims ()) + { + error ("interpn: incompatible size of argument number %d", i+1); + return retval; + } + else + { + NDArray tmp = NDArray (dim_vector (size[i], 1)); + + for (octave_idx_type j = 0; j < size[i]; j++) + tmp(j) = X[i](scale[i]*j); + + X[i] = tmp; + } + } + } + + for (int i = 0; i < n; i++) + { + if (! isvector (X[i]) && X[i].numel () != size[i]) + { + error ("interpn: incompatible size of argument number %d", i+1); + return retval; + } + else + x[i] = X[i].data (); + } + + lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi); + + retval = Vi; + + return retval; +}
deleted file mode 100644 --- a/src/DLD-FUNCTIONS/interpn.cc +++ /dev/null @@ -1,261 +0,0 @@ -/* - -Copyright (C) 2007 Alexander Barth - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 2, or (at your option) any -later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, write to the Free -Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA -02110-1301, USA. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include "dNDArray.h" - -#include "defun-dld.h" -#include "error.h" -#include "oct-obj.h" - -// equivalent to isvector.m - -bool -isvector (const NDArray& array) -{ - const dim_vector dv = array.dims (); - return dv.length () == 2 && (dv(0) == 1 || dv(1) == 1); -} - -// lookup a value in a sorted table (lookup.m) -octave_idx_type -lookup (const double *x, octave_idx_type n, double y) -{ - octave_idx_type j, j0, j1; - - if (y > x[n-1] || y < x[0]) - return -1; - -#ifdef EXHAUSTIF - for (j = 0; j < n - 1; j++) - { - if (x[j] <= y && y <= x[j+1]) - return j; - } -#else - j0 = 0; - j1 = n - 1; - - while (true) - { - j = (j0+j1)/2; - - if (y <= x[j+1]) - { - if (x[j] <= y) - return j; - - j1 = j; - } - - if (x[j] <= y) - j0 = j; - } - -#endif -} - -// n-dimensional linear interpolation - -void -lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale, - octave_idx_type Ni, double extrapval, const double **x, - const double *v, const double **y, double *vi) -{ - bool out = false; - int bit; - - OCTAVE_LOCAL_BUFFER (double, coef, 2*n); - OCTAVE_LOCAL_BUFFER (octave_idx_type, index, n); - - // loop over all points - for (octave_idx_type m = 0; m < Ni; m++) - { - // loop over all dimensions - for (int i = 0; i < n; i++) - { - index[i] = lookup (x[i], size[i], y[i][m]); - out = index[i] == -1; - - if (out) - break; - else - { - octave_idx_type j = index[i]; - coef[2*i+1] = (y[i][m] - x[i][j])/(x[i][j+1] - x[i][j]); - coef[2*i] = 1 - coef[2*i+1]; - } - } - - - if (out) - vi[m] = extrapval; - else - { - vi[m] = 0; - - // loop over all corners of hypercube (1<<n = 2^n) - for (int i = 0; i < (1 << n); i++) - { - double c = 1; - octave_idx_type l = 0; - - // loop over all dimensions - for (int j = 0; j < n; j++) - { - // test if the jth bit in i is set - bit = i >> j & 1; - l += scale[j] * (index[j] + bit); - c *= coef[2*j+bit]; - } - - vi[m] += c * v[l]; - } - } - } -} - -DEFUN_DLD (interpn, args, , - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{vi} =} interpn (@var{x1}, @var{x2}, @dots{}, @var{xn}, @var{v}, @var{y1}, @var{y2}, @dots{}, @var{yn})\n\ -Perform @var{n}-dimensional interpolation. Each element of then\n\ -@var{n}-dimensional array @var{v} represents a value at a location\n\ -given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters\n\ -@var{x1}, @var{x2}, @dots{}, @var{xn} are either @var{n}-dimensional\n\ -arrays of the same size as the array @var{v} in the \"ndgrid\" format\n\ -or vectors. The parameters @var{y1}, @var{y2}, @dots{}, @var{yn} are\n\ -all @var{n}-dimensional arrays of the same size and represent the\n\ -points at which the array @var{vi} is interpolated.\n\ -\n\ -This function only performs linear interpolation.\n\ -@seealso{interp1, interp2, ndgrid}\n\ -@end deftypefn") -{ - octave_value retval; - - int nargin = args.length (); - - if (nargin < 2 || nargin % 2 == 0) - { - print_usage (); - return retval; - } - - // dimension of the problem - int n = (nargin-1)/2; - - OCTAVE_LOCAL_BUFFER (NDArray, X, n); - OCTAVE_LOCAL_BUFFER (NDArray, Y, n); - - OCTAVE_LOCAL_BUFFER (const double *, x, n); - OCTAVE_LOCAL_BUFFER (const double *, y, n); - OCTAVE_LOCAL_BUFFER (octave_idx_type, scale, n); - OCTAVE_LOCAL_BUFFER (octave_idx_type, size, n); - - const NDArray V = args(n).array_value (); - NDArray Vi = NDArray (args(n+1).dims ()); - - if (error_state) - { - print_usage (); - return retval; - } - - const double *v = V.data (); - double *vi = Vi.fortran_vec (); - octave_idx_type Ni = Vi.numel (); - - double extrapval = octave_NaN; - - for (int i = 0; i < n; i++) - { - X[i] = args(i).array_value (); - Y[i] = args(n+i+1).array_value (); - - if (error_state) - { - print_usage (); - return retval; - } - - y[i] = Y[i].data (); - size[i] = V.dims()(i); - - if (Y[0].dims () != Y[i].dims ()) - { - error ("interpn: incompatible size of argument number %d", n+i+2); - return retval; - } - } - - // offset in memory of each dimension - - scale[0] = 1; - - for (int i = 1; i < n; i++) - scale[i] = scale[i-1] * size[i-1]; - - // tests if X[0] is a vector, if yes, assume that all elements of X are - // in the ndgrid format. - - if (! isvector (X[0])) - { - for (int i = 0; i < n; i++) - { - if (X[i].dims () != V.dims ()) - { - error ("interpn: incompatible size of argument number %d", i+1); - return retval; - } - else - { - NDArray tmp = NDArray (dim_vector (size[i], 1)); - - for (octave_idx_type j = 0; j < size[i]; j++) - tmp(j) = X[i](scale[i]*j); - - X[i] = tmp; - } - } - } - - for (int i = 0; i < n; i++) - { - if (! isvector (X[i]) && X[i].numel () != size[i]) - { - error ("interpn: incompatible size of argument number %d", i+1); - return retval; - } - else - x[i] = X[i].data (); - } - - lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi); - - retval = Vi; - - return retval; -}
--- a/src/Makefile.in +++ b/src/Makefile.in @@ -52,12 +52,13 @@ dassl.cc det.cc dispatch.cc eig.cc expm.cc fft.cc fft2.cc \ fftn.cc fftw.cc filter.cc find.cc fsolve.cc \ gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \ - givens.cc hess.cc interpn.cc inv.cc kron.cc lpsolve.cc lsode.cc \ + givens.cc hess.cc inv.cc kron.cc lpsolve.cc lsode.cc \ lu.cc luinc.cc matrix_type.cc md5sum.cc minmax.cc pinv.cc qr.cc \ quad.cc qz.cc rand.cc regexp.cc schur.cc sort.cc sparse.cc \ spchol.cc spdet.cc spfind.cc spkron.cc splu.cc spparms.cc spqr.cc \ sqrtm.cc svd.cc syl.cc symrcm.cc time.cc urlwrite.cc __contourc__.cc \ - __gnuplot_raw__.l __glpk__.cc __pchip_deriv__.cc __qp__.cc + __gnuplot_raw__.l __glpk__.cc __lin_interpn__.cc __pchip_deriv__.cc \ + __qp__.cc DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))