octave_packages/m/general/interp3.m

   1 ## Copyright (C) 2007-2012 David Bateman
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {@var{vi} =} interp3 (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi})
  21 ## @deftypefnx {Function File} {@var{vi} =} interp3 (@var{v}, @var{xi}, @var{yi}, @var{zi})
  22 ## @deftypefnx {Function File} {@var{vi} =} interp3 (@var{v}, @var{m})
  23 ## @deftypefnx {Function File} {@var{vi} =} interp3 (@var{v})
  24 ## @deftypefnx {Function File} {@var{vi} =} interp3 (@dots{}, @var{method})
  25 ## @deftypefnx {Function File} {@var{vi} =} interp3 (@dots{}, @var{method}, @var{extrapval})
  26 ##
  27 ## Perform 3-dimensional interpolation.  Each element of the 3-dimensional
  28 ## array @var{v} represents a value at a location given by the parameters
  29 ## @var{x}, @var{y}, and @var{z}.  The parameters @var{x}, @var{x}, and
  30 ## @var{z} are either 3-dimensional arrays of the same size as the array
  31 ## @var{v} in the 'meshgrid' format or vectors.  The parameters @var{xi}, etc.
  32 ## respect a similar format to @var{x}, etc., and they represent the points
  33 ## at which the array @var{vi} is interpolated.
  34 ##
  35 ## If @var{x}, @var{y}, @var{z} are omitted, they are assumed to be
  36 ## @code{x = 1 : size (@var{v}, 2)}, @code{y = 1 : size (@var{v}, 1)} and
  37 ## @code{z = 1 : size (@var{v}, 3)}.  If @var{m} is specified, then
  38 ## the interpolation adds a point half way between each of the interpolation
  39 ## points.  This process is performed @var{m} times.  If only @var{v} is
  40 ## specified, then @var{m} is assumed to be @code{1}.
  41 ##
  42 ## Method is one of:
  43 ##
  44 ## @table @asis
  45 ## @item 'nearest'
  46 ## Return the nearest neighbor.
  47 ##
  48 ## @item 'linear'
  49 ## Linear interpolation from nearest neighbors.
  50 ##
  51 ## @item 'cubic'
  52 ## Cubic interpolation from four nearest neighbors (not implemented yet).
  53 ##
  54 ## @item 'spline'
  55 ## Cubic spline interpolation---smooth first and second derivatives
  56 ## throughout the curve.
  57 ## @end table
  58 ##
  59 ## The default method is 'linear'.
  60 ##
  61 ## If @var{extrap} is the string 'extrap', then extrapolate values beyond
  62 ## the endpoints.  If @var{extrap} is a number, replace values beyond the
  63 ## endpoints with that number.  If @var{extrap} is missing, assume NA.
  64 ## @seealso{interp1, interp2, spline, meshgrid}
  65 ## @end deftypefn
  66
  67 function vi = interp3 (varargin)
  68   method = "linear";
  69   extrapval = NA;
  70   nargs = nargin;
  71
  72   if (nargin < 1 || ! isnumeric (varargin{1}))
  73     print_usage ();
  74   endif
  75
  76   if (ischar (varargin{end}))
  77     method = varargin{end};
  78     nargs = nargs - 1;
  79   elseif (nargs > 1 && ischar (varargin{end - 1}))
  80     if (! isnumeric (varargin{end}) || ! isscalar (varargin{end}))
  81       error ("interp3: extrapal is expected to be a numeric scalar");
  82     endif
  83     extrapval = varargin{end};
  84     method = varargin{end-1};
  85     nargs = nargs - 2;
  86   endif
  87
  88   if (nargs < 3 || (nargs == 4 && ! isvector (varargin{1})
  89                     && nargs == (ndims (varargin{1}) + 1)))
  90     v = varargin{1};
  91     if (ndims (v) != 3)
  92       error ("interp3: expect 3-dimensional array of values");
  93     endif
  94     x = varargin (2:end);
  95     if (any (! cellfun (@isvector, x)))
  96       for i = 2 : 3
  97         if (! size_equal (x{1}, x{i}))
  98           error ("interp3: dimensional mismatch");
  99         endif
 100         x{i} = permute (x{i}, [2, 1, 3]);
 101       endfor
 102       x{1} = permute (x{1}, [2, 1, 3]);
 103     endif
 104     v = permute (v, [2, 1, 3]);
 105     vi = ipermute (interpn (v, x{:}, method, extrapval), [2, 1, 3]);
 106   elseif (nargs == 7 && nargs == (2 * ndims (varargin{ceil (nargs / 2)})) + 1)
 107     v = varargin{4};
 108     if (ndims (v) != 3)
 109       error ("interp3: expect 3-dimensional array of values");
 110     endif
 111     x = varargin (1:3);
 112     if (any (! cellfun (@isvector, x)))
 113       for i = 2 : 3
 114         if (! size_equal (x{1}, x{i}) || ! size_equal (x{i}, v))
 115           error ("interp3: dimensional mismatch");
 116         endif
 117         x{i} = permute (x{i}, [2, 1, 3]);
 118       endfor
 119       x{1} = permute (x{1}, [2, 1, 3]);
 120     endif
 121     y = varargin (5:7);
 122     if (any (! cellfun (@isvector, y)))
 123       for i = 2 : 3
 124         if (! size_equal (y{1}, y{i}))
 125           error ("interp3: dimensional mismatch");
 126         endif
 127         y{i} = permute (y{i}, [2, 1, 3]);
 128       endfor
 129       y{1} = permute (y{1}, [2, 1, 3]);
 130     endif
 131     v = permute (v, [2, 1, 3]);
 132     vi = ipermute (interpn (x{:}, v, y{:}, method, extrapval), [2, 1, 3]);
 133   else
 134     error ("interp3: wrong number or incorrectly formatted input arguments");
 135   endif
 136 endfunction
 137
 138 %!test
 139 %! x = y = z = -1:1;
 140 %! f = @(x,y,z) x.^2 - y - z.^2;
 141 %! [xx, yy, zz] = meshgrid (x, y, z);
 142 %! v = f (xx,yy,zz);
 143 %! xi = yi = zi = -1:0.5:1;
 144 %! [xxi, yyi, zzi] = meshgrid (xi, yi, zi);
 145 %! vi = interp3(x, y, z, v, xxi, yyi, zzi);
 146 %! [xxi, yyi, zzi] = ndgrid (xi, yi, zi);
 147 %! vi2 = interpn(x, y, z, v, xxi, yyi, zzi);
 148 %! assert (vi, vi2);
 149
 150 %!shared z, zout, tol
 151 %! z = zeros (3, 3, 3);
 152 %! zout = zeros (5, 5, 5);
 153 %! z(:,:,1) = [1 3 5; 3 5 7; 5 7 9];
 154 %! z(:,:,2) = z(:,:,1) + 2;
 155 %! z(:,:,3) = z(:,:,2) + 2;
 156 %! for n = 1:5
 157 %!   zout(:,:,n) = [1 2 3 4 5;
 158 %!                  2 3 4 5 6;
 159 %!                  3 4 5 6 7;
 160 %!                  4 5 6 7 8;
 161 %!                  5 6 7 8 9] + (n-1);
 162 %! end
 163 %! tol = 10 * eps;
 164 %!assert (interp3 (z), zout, tol)
 165 %!assert (interp3 (z, "linear"), zout, tol)
 166 %!assert (interp3 (z, "spline"), zout, tol)