Mercurial > octave
view scripts/general/interpn.m @ 14237:11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Add clf() to all demos using plot features to get reproducibility.
Use 64 as input to all colormaps (jet (64)) to get reproducibility.
* bicubic.m, cell2mat.m, celldisp.m, cplxpair.m, interp1.m, interp2.m,
interpft.m, interpn.m, profile.m, profshow.m, convhull.m, delaunay.m,
griddata.m, inpolygon.m, voronoi.m, autumn.m, bone.m, contrast.m, cool.m,
copper.m, flag.m, gmap40.m, gray.m, hot.m, hsv.m, image.m, imshow.m, jet.m,
ocean.m, pink.m, prism.m, rainbow.m, spring.m, summer.m, white.m, winter.m,
condest.m, onenormest.m, axis.m, clabel.m, colorbar.m, comet.m, comet3.m,
compass.m, contour.m, contour3.m, contourf.m, cylinder.m, daspect.m,
ellipsoid.m, errorbar.m, ezcontour.m, ezcontourf.m, ezmesh.m, ezmeshc.m,
ezplot.m, ezplot3.m, ezpolar.m, ezsurf.m, ezsurfc.m, feather.m, fill.m,
fplot.m, grid.m, hold.m, isosurface.m, legend.m, loglog.m, loglogerr.m,
pareto.m, patch.m, pbaspect.m, pcolor.m, pie.m, pie3.m, plot3.m, plotmatrix.m,
plotyy.m, polar.m, quiver.m, quiver3.m, rectangle.m, refreshdata.m, ribbon.m,
rose.m, scatter.m, scatter3.m, semilogx.m, semilogxerr.m, semilogy.m,
semilogyerr.m, shading.m, slice.m, sombrero.m, stairs.m, stem.m, stem3.m,
subplot.m, surf.m, surfc.m, surfl.m, surfnorm.m, text.m, title.m, trimesh.m,
triplot.m, trisurf.m, uigetdir.m, uigetfile.m, uimenu.m, uiputfile.m,
waitbar.m, xlim.m, ylim.m, zlim.m, mkpp.m, pchip.m, polyaffine.m, spline.m,
bicgstab.m, cgs.m, gplot.m, pcg.m, pcr.m, treeplot.m, strtok.m, demo.m,
example.m, rundemos.m, speed.m, test.m, calendar.m, datestr.m, datetick.m,
weekday.m: Revamp %!demos to use Octave coding conventions on spacing, etc.
author | Rik <octave@nomad.inbox5.com> |
---|---|
date | Fri, 20 Jan 2012 12:59:53 -0800 |
parents | 72c96de7a403 |
children | c4fa5e0b6193 |
line wrap: on
line source
## Copyright (C) 2007-2012 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 3 of the License, 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, see ## <http://www.gnu.org/licenses/>. ## -*- 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 the @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 omitted, 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 interpolation ## 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 neighbor. ## ## @item 'linear' ## Linear interpolation from nearest neighbors. ## ## @item 'cubic' ## Cubic interpolation from four nearest neighbors (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{extrapval} is the scalar value, use it to replace the values ## beyond the endpoints with that number. If @var{extrapval} is missing, ## assume NA. ## @seealso{interp1, interp2, spline, ndgrid} ## @end deftypefn function vi = interpn (varargin) method = "linear"; extrapval = NA; nargs = nargin; if (nargin < 1 || ! isnumeric (varargin{1})) print_usage (); endif if (ischar (varargin{end})) method = varargin{end}; nargs = nargs - 1; elseif (nargs > 1 && ischar (varargin{end - 1})) if (! isnumeric (varargin{end}) || ! isscalar (varargin{end})) error ("interpn: extrapal is expected to be a numeric scalar"); endif method = varargin{end - 1}; extrapval = varargin{end}; nargs = nargs - 2; endif if (nargs < 3) v = varargin{1}; m = 1; if (nargs == 2) if (ischar (varargin{2})) method = varargin{2}; elseif (isnumeric (m) && isscalar (m) && fix (m) == m) m = varargin{2}; else print_usage (); 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 y{1} = y{1}.'; [y{:}] = ndgrid (y{:}); 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 ("interpn: 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 ("interpn: 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 method = tolower (method); all_vectors = all (cellfun ("isvector", y)); different_lengths = numel (unique (cellfun ("numel", y))) > 1; if (all_vectors && different_lengths) [foobar(1:numel(y)).y] = ndgrid (y{:}); y = {foobar.y}; endif if (strcmp (method, "linear")) vi = __lin_interpn__ (x{:}, v, y{:}); vi (isna (vi)) = 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}, y{i}, "lr"); 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")) if (any (! cellfun ("isvector", y))) for i = 2 : nd if (! size_equal (y{1}, y{i})) error ("interpn: dimensional mismatch"); endif idx (1 : nd) = {1}; idx (i) = ":"; y{i} = y{i}(idx{:}); endfor idx (1 : nd) = {1}; idx (1) = ":"; y{1} = y{1}(idx{:}); endif vi = __splinen__ (x, v, y, extrapval, "interpn"); if (size_equal (y{:})) ly = length (y{1}); idx = cell (1, ly); q = cell (1, nd); for i = 1 : ly q(:) = i; idx {i} = q; endfor vi = vi (cellfun (@(x) sub2ind (size(vi), x{:}), idx)); vi = reshape (vi, size(y{1})); endif elseif (strcmp (method, "cubic")) error ("interpn: cubic interpolation not yet implemented"); else error ("interpn: unrecognized interpolation METHOD"); endif endfunction %!demo %! clf; %! 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 %! clf; %! 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 # FIXME: Uncomment when support for "cubic" has been added %! clf; %! 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 %! clf; %! 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; %!demo %! clf; %! 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.1:1; %! [xxi, yyi, zzi] = ndgrid (xi, yi, zi); %! vi = interpn (x, y, z, v, xxi, yyi, zzi, 'spline'); %! mesh (yi, zi, squeeze (vi(1,:,:))); %!test %! [x,y,z] = ndgrid (0:2); %! f = x + y + z; %! assert (interpn (x,y,z,f,[.5 1.5],[.5 1.5],[.5 1.5]), [1.5, 4.5]); %! assert (interpn (x,y,z,f,[.51 1.51],[.51 1.51],[.51 1.51],"nearest"), [3, 6]); %! assert (interpn (x,y,z,f,[.5 1.5],[.5 1.5],[.5 1.5],"spline"), [1.5, 4.5]); %! assert (interpn (x,y,z,f,x,y,z), f); %! assert (interpn (x,y,z,f,x,y,z,"nearest"), f); %! assert (interpn (x,y,z,f,x,y,z,"spline"), f); %!test %! [x, y, z] = ndgrid (0:2, 1:4, 2:6); %! f = x + y + z; %! xi = [0.5 1.0 1.5]; yi = [1.5 2.0 2.5 3.5]; zi = [2.5 3.5 4.0 5.0 5.5]; %! fi = interpn (x, y, z, f, xi, yi, zi); %! [xi, yi, zi] = ndgrid (xi, yi, zi); %! assert (fi, xi + yi + zi); %!test %! xi = 0:2; yi = 1:4; zi = 2:6; %! [x, y, z] = ndgrid (xi, yi, zi); %! f = x + y + z; %! fi = interpn (x, y, z, f, xi, yi, zi, "nearest"); %! assert (fi, x + y + z); %!test %! [x,y,z] = ndgrid (0:2); %! f = x.^2 + y.^2 + z.^2; %! assert (interpn (x,y,-z,f,1.5,1.5,-1.5), 7.5); %!test # for Matlab-compatible rounding for "nearest" %! x = meshgrid (1:4); %! assert (interpn (x, 2.5, 2.5, "nearest"), 3); %!test %! z = zeros (3, 3, 3); %! zout = zeros (5, 5, 5); %! z(:,:,1) = [1 3 5; 3 5 7; 5 7 9]; %! z(:,:,2) = z(:,:,1) + 2; %! z(:,:,3) = z(:,:,2) + 2; %! for n = 1:5 %! zout(:,:,n) = [1 2 3 4 5; %! 2 3 4 5 6; %! 3 4 5 6 7; %! 4 5 6 7 8; %! 5 6 7 8 9] + (n-1); %! endfor %! tol = 10*eps; %! assert (interpn (z), zout, tol); %! assert (interpn (z, "linear"), zout, tol); %! assert (interpn (z, "spline"), zout, tol);