Mercurial > octave-nkf
view scripts/general/interp3.m @ 18840:4a4edf0f2077 nkf-ready
fix LLVM 3.4 build (bug #41061)
* configure.ac: Call new functions OCTAVE_LLVM_RAW_FD_OSTREAM_API and
OCTAVE_LLVM_LEGACY_PASSMANAGER_API, check for Verifier.h header file
* m4/acinclude.m4 (OCTAVE_LLVM_RAW_FD_OSTREAM_API): New function to
detect correct raw_fd_ostream API
* m4/acinclude.m4 (OCTAVE_LLVM_LEGACY_PASSMANAGER_API): New function
to detect legacy passmanager API
* libinterp/corefcn/jit-util.h: Use legacy passmanager namespace if
necessary
* libinterp/corefcn/pt-jit.h (class tree_jit): Use legacy passmanager
class if necessary
* libinterp/corefcn/pt-jit.cc: Include appropriate header files
* libinterp/corefcn/pt-jit.cc (tree_jit::initialize): Use legacy
passmanager if necessary
* libinterp/corefcn/pt-jit.cc (tree_jit::optimize): Use correct API
* libinterp/corefcn/jit-typeinfo.cc: Include appropriate header file
author | Stefan Mahr <dac922@gmx.de> |
---|---|
date | Sun, 11 May 2014 02:28:33 +0200 |
parents | 0ede4dbb37f1 |
children | 0e1f5a750d00 |
line wrap: on
line source
## Copyright (C) 2007-2013 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} =} 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{n}) ## @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}) ## ## Three-dimensional interpolation. ## ## Interpolate reference data @var{x}, @var{y}, @var{z}, @var{v} to determine ## @var{vi} at the coordinates @var{xi}, @var{yi}, @var{zi}. The reference ## data @var{x}, @var{y}, @var{z} can be matrices, as returned by ## @code{meshgrid}, in which case the sizes of ## @var{x}, @var{y}, @var{z}, and @var{v} must be equal. If @var{x}, @var{y}, ## @var{z} are vectors describing a cubic grid then ## @code{length (@var{x}) == columns (@var{v})}, ## @code{length (@var{y}) == rows (@var{v})}, ## and @code{length (@var{z}) == size (@var{v}, 3)}. In either case the input ## data must be strictly monotonic. ## ## If called without @var{x}, @var{y}, @var{z}, and just a single reference ## data matrix @var{v}, the 3-D region ## @code{@var{x} = 1:columns (@var{v}), @var{y} = 1:rows (@var{v}), ## @var{z} = 1:size (@var{v}, 3)} is assumed. ## This saves memory if the grid is regular and the distance between points is ## not important. ## ## If called with a single reference data matrix @var{v} and a refinement ## value @var{n}, then perform interpolation over a 3-D grid where each original ## interval has been recursively subdivided @var{n} times. This results in ## @code{2^@var{n}-1} additional points for every interval in the original ## grid. If @var{n} is omitted a value of 1 is used. As an example, the ## interval [0,1] with @code{@var{n}==2} results in a refined interval with ## points at [0, 1/4, 1/2, 3/4, 1]. ## ## The interpolation @var{method} is one of: ## ## @table @asis ## @item @qcode{"nearest"} ## Return the nearest neighbor. ## ## @item @qcode{"linear"} (default) ## Linear interpolation from nearest neighbors. ## ## @item @qcode{"pchip"} ## Piecewise cubic Hermite interpolating polynomial---shape-preserving ## interpolation with smooth first derivative (not implemented yet). ## ## @item @qcode{"cubic"} ## Cubic interpolation (same as @qcode{"pchip"} [not implemented yet]). ## ## @item @qcode{"spline"} ## Cubic spline interpolation---smooth first and second derivatives ## throughout the curve. ## @end table ## ## If @var{extrapval} is a number, then replace values beyond the endpoints ## with that number. When unspecified, @var{extrapval} defaults to @code{NA}. ## Note that if @var{extrapval} is used, @var{method} must be specified as well. ## @seealso{interp1, interp2, interpn, meshgrid} ## @end deftypefn ## FIXME: Need to validate N argument (maybe change interpn). ## FIXME: Need to add support for 'pchip' method (maybe change interpn). ## FIXME: Need to add support for "extrap" string value (maybe change interpn). function vi = interp3 (varargin) method = "linear"; extrapval = NA; nargs = nargin; if (nargin < 1 || ! isnumeric (varargin{1})) print_usage (); endif if (ischar (varargin{end})) method = varargin{end}; nargs--; elseif (nargs > 1 && ischar (varargin{end-1})) ## FIXME: No support for "extrap" string if (! isnumeric (varargin{end}) || ! isscalar (varargin{end})) error ("interp3: EXTRAPVAL must be a numeric scalar"); endif extrapval = varargin{end}; method = varargin{end-1}; nargs -= 2; endif if (method(1) == "*") warning ("interp3: ignoring unsupported '*' flag to METHOD"); method(1) = []; endif if (nargs < 3) ## Calling form interp3 (v) OR interp3 (v, n) v = varargin{1}; if (ndims (v) != 3) error ("interp3: V must be a 3-D array of values"); endif n = varargin(2:nargs); v = permute (v, [2, 1, 3]); vi = ipermute (interpn (v, n{:}, method, extrapval), [2, 1, 3]); elseif (nargs == 4 && ! isvector (varargin{1})) ## Calling form interp3 (v, xi, yi, zi) v = varargin{1}; if (ndims (v) != 3) error ("interp3: V must be a 3-D array of values"); endif xi = varargin(2:4); if (any (! cellfun (@isvector, xi))) ## Meshgridded values rather than vectors if (! size_equal (xi{:})) error ("interp3: XI, YI, and ZI dimensions must be equal"); endif for i = 1 : 3 xi{i} = permute (xi{i}, [2, 1, 3]); endfor endif v = permute (v, [2, 1, 3]); vi = ipermute (interpn (v, xi{:}, method, extrapval), [2, 1, 3]); elseif (nargs == 7) ## Calling form interp3 (x, y, z, v, xi, yi, zi) v = varargin{4}; if (ndims (v) != 3) error ("interp3: V must be a 3-D array of values"); endif x = varargin(1:3); if (any (! cellfun (@isvector, x))) ## Meshgridded values rather than vectors if (! size_equal (x{:}, v)) error ("interp3: X, Y, Z, and V dimensions must be equal"); endif for i = 1 : 3 x{i} = permute (x{i}, [2, 1, 3]); endfor endif xi = varargin(5:7); if (any (! cellfun (@isvector, xi))) ## Meshgridded values rather than vectors if (! size_equal (xi{:})) error ("interp3: XI, YI, and ZI dimensions must be equal"); endif for i = 1 : 3 xi{i} = permute (xi{i}, [2, 1, 3]); endfor endif v = permute (v, [2, 1, 3]); vi = ipermute (interpn (x{:}, v, xi{:}, method, extrapval), [2, 1, 3]); else error ("interp3: wrong number or incorrectly formatted input arguments"); endif endfunction %% FIXME: Need some demo blocks here to show off the function like interp2.m. %!test # basic test %! x = y = z = -1:1; y = y + 2; %! 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; yi = yi + 2.1; %! [xxi, yyi, zzi] = meshgrid (xi, yi, zi); %! vi = interp3 (x, y, z, v, xxi, yyi, zzi); %! [xxi, yyi, zzi] = ndgrid (yi, xi, zi); %! vi2 = interpn (y, x, z, v, xxi, yyi, zzi); %! assert (vi, vi2, 10*eps); %!test # meshgridded xi, yi, zi %! x = z = 1:2; y = 1:3; %! v = ones ([3,2,2]); v(:,2,1) = [7;5;4]; v(:,1,2) = [2;3;5]; %! xi = zi = .6:1.6; yi = 1; %! [xxi3, yyi3, zzi3] = meshgrid (xi, yi, zi); %! [xxi, yyi, zzi] = ndgrid (yi, xi, zi); %! vi = interp3 (x, y, z, v, xxi3, yyi3, zzi3, "nearest"); %! vi2 = interpn (y, x, z, v, xxi, yyi, zzi, "nearest"); %! assert (vi, vi2); %!test # vector xi, yi, zi %! x = z = 1:2; y = 1:3; %! v = ones ([3,2,2]); v(:,2,1) = [7;5;4]; v(:,1,2) = [2;3;5]; %! xi = zi = .6:1.6; yi = 1; %! vi = interp3 (x, y, z, v, xi, yi, zi, "nearest"); %! vi2 = interpn (y, x, z, v, yi, xi, zi,"nearest"); %! assert (vi, vi2); %!test # vector xi+1 with extrap value %! x = z = 1:2; y = 1:3; %! v = ones ([3,2,2]); v(:,2,1) = [7;5;4]; v(:,1,2) = [2;3;5]; %! xi = zi = .6:1.6; yi = 1; %! vi = interp3 (x, y, z, v, xi+1, yi, zi, "nearest", 3); %! vi2 = interpn (y, x, z, v, yi, xi+1, zi, "nearest", 3); %! assert (vi, vi2); %!test # input value matrix--no x,y,z %! x = z = 1:2; y = 1:3; %! v = ones ([3,2,2]); v(:,2,1) = [7;5;4]; v(:,1,2) = [2;3;5]; %! xi = zi = .6:1.6; yi = 1; %! vi = interp3 (v, xi, yi, zi, "nearest"); %! vi2 = interpn (v, yi, xi, zi,"nearest"); %! assert (vi, vi2); %!test # input value matrix--no x,y,z, with extrap value %! x = z = 1:2; y = 1:3; %! v = ones ([3,2,2]); v(:,2,1) = [7;5;4]; v(:,1,2) = [2;3;5]; %! xi = zi = .6:1.6; yi = 1; %! vi = interp3 (v, xi, yi, zi, "nearest", 3); %! vi2 = interpn (v, yi, xi, zi, "nearest", 3); %! assert (vi, vi2); %!shared z, zout, tol %! 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); %! end %! tol = 10 * eps; %! %!assert (interp3 (z), zout, tol) %!assert (interp3 (z, "linear"), zout, tol) %!assert (interp3 (z, "spline"), zout, tol) %% Test input validation %!error interp3 () %!error interp3 ({1}) %!error <EXTRAPVAL must be a numeric scalar> interp3 (1,2,3,4,1,2,3,"linear", {1}) %!error <EXTRAPVAL must be a numeric scalar> interp3 (1,2,3,4,1,2,3,"linear", ones (2,2)) %!warning <ignoring unsupported '\*' flag> interp3 (rand (3,3,3), 1, "*linear"); %!error <V must be a 3-D array> interp3 (rand (2,2)) %!error <V must be a 3-D array> interp3 (rand (2,2), 1,1,1) %!error <XI, YI, and ZI dimensions must be equal> interp3 (rand (2,2,2), 1,1, ones (2,2)) %!error <V must be a 3-D array> interp3 (1:2, 1:2, 1:2, rand (2,2), 1,1,1) %!error <X, Y, Z, and V dimensions must be equal> interp3 (ones(1,2,2), ones(2,2,2), ones(2,2,2), rand (2,2,2), 1,1,1) %!error <XI, YI, and ZI dimensions must be equal> interp3 (1:2, 1:2, 1:2, rand (2,2,2), 1,1, ones (2,2))