octave_packages/m/general/interpn.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} =} interpn (@var{x1}, @var{x2}, @dots{}, @var{v}, @var{y1}, @var{y2}, @dots{})
  21 ## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}, @var{y1}, @var{y2}, @dots{})
  22 ## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}, @var{m})
  23 ## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v})
  24 ## @deftypefnx {Function File} {@var{vi} =} interpn (@dots{}, @var{method})
  25 ## @deftypefnx {Function File} {@var{vi} =} interpn (@dots{}, @var{method}, @var{extrapval})
  26 ##
  27 ## Perform @var{n}-dimensional interpolation, where @var{n} is at least two.
  28 ## Each element of the @var{n}-dimensional array @var{v} represents a value
  29 ## at a location given by the parameters @var{x1}, @var{x2}, @dots{}, @var{xn}.
  30 ## The parameters @var{x1}, @var{x2}, @dots{}, @var{xn} are either
  31 ## @var{n}-dimensional arrays of the same size as the array @var{v} in
  32 ## the 'ndgrid' format or vectors.  The parameters @var{y1}, etc. respect a
  33 ## similar format to @var{x1}, etc., and they represent the points at which
  34 ## the array @var{vi} is interpolated.
  35 ##
  36 ## If @var{x1}, @dots{}, @var{xn} are omitted, they are assumed to be
  37 ## @code{x1 = 1 : size (@var{v}, 1)}, etc.  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{extrapval} is the scalar value, use it to replace the values
  62 ## beyond the endpoints with that number.  If @var{extrapval} is missing,
  63 ## assume NA.
  64 ## @seealso{interp1, interp2, spline, ndgrid}
  65 ## @end deftypefn
  66
  67 function vi = interpn (varargin)
  68
  69   method = "linear";
  70   extrapval = NA;
  71   nargs = nargin;
  72
  73   if (nargin < 1 || ! isnumeric (varargin{1}))
  74     print_usage ();
  75   endif
  76
  77   if (ischar (varargin{end}))
  78     method = varargin{end};
  79     nargs = nargs - 1;
  80   elseif (nargs > 1 && ischar (varargin{end - 1}))
  81     if (! isnumeric (varargin{end}) || ! isscalar (varargin{end}))
  82       error ("interpn: extrapal is expected to be a numeric scalar");
  83     endif
  84     method = varargin{end - 1};
  85     extrapval = varargin{end};
  86     nargs = nargs - 2;
  87   endif
  88
  89   if (nargs < 3)
  90     v = varargin{1};
  91     m = 1;
  92     if (nargs == 2)
  93       if (ischar (varargin{2}))
  94         method = varargin{2};
  95       elseif (isnumeric (m) && isscalar (m) && fix (m) == m)
  96         m = varargin{2};
  97       else
  98         print_usage ();
  99       endif
 100     endif
 101     sz = size (v);
 102     nd = ndims (v);
 103     x = cell (1, nd);
 104     y = cell (1, nd);
 105     for i = 1 : nd;
 106       x{i} = 1 : sz(i);
 107       y{i} = 1 : (1 / (2 ^ m)) : sz(i);
 108     endfor
 109     y{1} = y{1}.';
 110     [y{:}] = ndgrid (y{:});
 111   elseif (! isvector (varargin{1}) && nargs == (ndims (varargin{1}) + 1))
 112     v = varargin{1};
 113     sz = size (v);
 114     nd = ndims (v);
 115     x = cell (1, nd);
 116     y = varargin (2 : nargs);
 117     for i = 1 : nd;
 118       x{i} = 1 : sz(i);
 119     endfor
 120   elseif (rem (nargs, 2) == 1 && nargs ==
 121           (2 * ndims (varargin{ceil (nargs / 2)})) + 1)
 122     nv = ceil (nargs / 2);
 123     v = varargin{nv};
 124     sz = size (v);
 125     nd = ndims (v);
 126     x = varargin (1 : (nv - 1));
 127     y = varargin ((nv + 1) : nargs);
 128   else
 129     error ("interpn: wrong number or incorrectly formatted input arguments");
 130   endif
 131
 132   if (any (! cellfun ("isvector", x)))
 133     for i = 2 : nd
 134       if (! size_equal (x{1}, x{i}) || ! size_equal (x{i}, v))
 135         error ("interpn: dimensional mismatch");
 136       endif
 137       idx (1 : nd) = {1};
 138       idx (i) = ":";
 139       x{i} = x{i}(idx{:})(:);
 140     endfor
 141     idx (1 : nd) = {1};
 142     idx (1) = ":";
 143     x{1} = x{1}(idx{:})(:);
 144   endif
 145
 146   method = tolower (method);
 147
 148   all_vectors = all (cellfun ("isvector", y));
 149   different_lengths = numel (unique (cellfun ("numel", y))) > 1;
 150   if (all_vectors && different_lengths)
 151     [foobar(1:numel(y)).y] = ndgrid (y{:});
 152     y = {foobar.y};
 153   endif
 154
 155   if (strcmp (method, "linear"))
 156     vi = __lin_interpn__ (x{:}, v, y{:});
 157     vi (isna (vi)) = extrapval;
 158   elseif (strcmp (method, "nearest"))
 159     yshape = size (y{1});
 160     yidx = cell (1, nd);
 161     for i = 1 : nd
 162       y{i} = y{i}(:);
 163       yidx{i} = lookup (x{i}, y{i}, "lr");
 164     endfor
 165     idx = cell (1,nd);
 166     for i = 1 : nd
 167       idx{i} = yidx{i} + (y{i} - x{i}(yidx{i})(:) >= x{i}(yidx{i} + 1)(:) - y{i});
 168     endfor
 169     vi = v (sub2ind (sz, idx{:}));
 170     idx = zeros (prod (yshape), 1);
 171     for i = 1 : nd
 172       idx |= y{i} < min (x{i}(:)) | y{i} > max (x{i}(:));
 173     endfor
 174     vi(idx) = extrapval;
 175     vi = reshape (vi, yshape);
 176   elseif (strcmp (method, "spline"))
 177     if (any (! cellfun ("isvector", y)))
 178       for i = 2 : nd
 179         if (! size_equal (y{1}, y{i}))
 180           error ("interpn: dimensional mismatch");
 181         endif
 182         idx (1 : nd) = {1};
 183         idx (i) = ":";
 184         y{i} = y{i}(idx{:});
 185       endfor
 186       idx (1 : nd) = {1};
 187       idx (1) = ":";
 188       y{1} = y{1}(idx{:});
 189     endif
 190
 191     vi = __splinen__ (x, v, y, extrapval, "interpn");
 192
 193     if (size_equal (y{:}))
 194       ly = length (y{1});
 195       idx = cell (1, ly);
 196       q = cell (1, nd);
 197       for i = 1 : ly
 198         q(:) = i;
 199         idx {i} = q;
 200       endfor
 201       vi = vi (cellfun (@(x) sub2ind (size(vi), x{:}), idx));
 202       vi = reshape (vi, size(y{1}));
 203     endif
 204   elseif (strcmp (method, "cubic"))
 205     error ("interpn: cubic interpolation not yet implemented");
 206   else
 207     error ("interpn: unrecognized interpolation METHOD");
 208   endif
 209
 210 endfunction
 211
 212 %!demo
 213 %! A=[13,-1,12;5,4,3;1,6,2];
 214 %! x=[0,1,4]; y=[10,11,12];
 215 %! xi=linspace(min(x),max(x),17);
 216 %! yi=linspace(min(y),max(y),26)';
 217 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"linear").');
 218 %! [x,y] = meshgrid(x,y);
 219 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 220
 221 %!demo
 222 %! A=[13,-1,12;5,4,3;1,6,2];
 223 %! x=[0,1,4]; y=[10,11,12];
 224 %! xi=linspace(min(x),max(x),17);
 225 %! yi=linspace(min(y),max(y),26)';
 226 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"nearest").');
 227 %! [x,y] = meshgrid(x,y);
 228 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 229
 230 %!#demo
 231 %! A=[13,-1,12;5,4,3;1,6,2];
 232 %! x=[0,1,2]; y=[10,11,12];
 233 %! xi=linspace(min(x),max(x),17);
 234 %! yi=linspace(min(y),max(y),26)';
 235 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"cubic").');
 236 %! [x,y] = meshgrid(x,y);
 237 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 238
 239 %!demo
 240 %! A=[13,-1,12;5,4,3;1,6,2];
 241 %! x=[0,1,2]; y=[10,11,12];
 242 %! xi=linspace(min(x),max(x),17);
 243 %! yi=linspace(min(y),max(y),26)';
 244 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"spline").');
 245 %! [x,y] = meshgrid(x,y);
 246 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 247
 248
 249 %!demo
 250 %! x = y = z = -1:1;
 251 %! f = @(x,y,z) x.^2 - y - z.^2;
 252 %! [xx, yy, zz] = meshgrid (x, y, z);
 253 %! v = f (xx,yy,zz);
 254 %! xi = yi = zi = -1:0.1:1;
 255 %! [xxi, yyi, zzi] = ndgrid (xi, yi, zi);
 256 %! vi = interpn(x, y, z, v, xxi, yyi, zzi, 'spline');
 257 %! mesh (yi, zi, squeeze (vi(1,:,:)));
 258
 259
 260 %!test
 261 %! [x,y,z] = ndgrid(0:2);
 262 %! f = x+y+z;
 263 %! assert (interpn(x,y,z,f,[.5 1.5],[.5 1.5],[.5 1.5]), [1.5, 4.5])
 264 %! assert (interpn(x,y,z,f,[.51 1.51],[.51 1.51],[.51 1.51],'nearest'), [3, 6])
 265 %! assert (interpn(x,y,z,f,[.5 1.5],[.5 1.5],[.5 1.5],'spline'), [1.5, 4.5])
 266 %! assert (interpn(x,y,z,f,x,y,z), f)
 267 %! assert (interpn(x,y,z,f,x,y,z,'nearest'), f)
 268 %! assert (interpn(x,y,z,f,x,y,z,'spline'), f)
 269
 270 %!test
 271 %! [x, y, z] = ndgrid (0:2, 1:4, 2:6);
 272 %! f = x + y + z;
 273 %! xi = [0.5 1.0 1.5];
 274 %! yi = [1.5 2.0 2.5 3.5];
 275 %! zi = [2.5 3.5 4.0 5.0 5.5];
 276 %! fi = interpn (x, y, z, f, xi, yi, zi);
 277 %! [xi, yi, zi] = ndgrid (xi, yi, zi);
 278 %! assert (fi, xi + yi + zi)
 279
 280 %!test
 281 %! xi = 0:2;
 282 %! yi = 1:4;
 283 %! zi = 2:6;
 284 %! [x, y, z] = ndgrid (xi, yi, zi);
 285 %! f = x + y + z;
 286 %! fi = interpn (x, y, z, f, xi, yi, zi, "nearest");
 287 %! assert (fi, x + y + z)
 288
 289 %!test
 290 %! [x,y,z] = ndgrid(0:2);
 291 %! f = x.^2+y.^2+z.^2;
 292 %! assert (interpn(x,y,-z,f,1.5,1.5,-1.5), 7.5)
 293
 294 %!test % for Matlab-compatible rounding for 'nearest'
 295 %! X = meshgrid (1:4);
 296 %! assert (interpn (X, 2.5, 2.5, 'nearest'), 3);
 297
 298 %!shared z, zout, tol
 299 %! z = zeros (3, 3, 3);
 300 %! zout = zeros (5, 5, 5);
 301 %! z(:,:,1) = [1 3 5; 3 5 7; 5 7 9];
 302 %! z(:,:,2) = z(:,:,1) + 2;
 303 %! z(:,:,3) = z(:,:,2) + 2;
 304 %! for n = 1:5
 305 %!   zout(:,:,n) = [1 2 3 4 5;
 306 %!                  2 3 4 5 6;
 307 %!                  3 4 5 6 7;
 308 %!                  4 5 6 7 8;
 309 %!                  5 6 7 8 9] + (n-1);
 310 %! end
 311 %! tol = 10 * eps;
 312 %!assert (interpn (z), zout, tol)
 313 %!assert (interpn (z, "linear"), zout, tol)
 314 %!assert (interpn (z, "spline"), zout, tol)