1 ## Copyright (C) 2007-2012 David Bateman
3 ## This file is part of Octave.
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.
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.
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/>.
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})
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.
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}.
46 ## Return the nearest neighbor.
49 ## Linear interpolation from nearest neighbors.
52 ## Cubic interpolation from four nearest neighbors (not implemented yet).
55 ## Cubic spline interpolation---smooth first and second derivatives
56 ## throughout the curve.
59 ## The default method is 'linear'.
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,
64 ## @seealso{interp1, interp2, spline, ndgrid}
67 function vi = interpn (varargin)
73 if (nargin < 1 || ! isnumeric (varargin{1}))
77 if (ischar (varargin{end}))
78 method = varargin{end};
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");
84 method = varargin{end - 1};
85 extrapval = varargin{end};
93 if (ischar (varargin{2}))
95 elseif (isnumeric (m) && isscalar (m) && fix (m) == m)
107 y{i} = 1 : (1 / (2 ^ m)) : sz(i);
110 [y{:}] = ndgrid (y{:});
111 elseif (! isvector (varargin{1}) && nargs == (ndims (varargin{1}) + 1))
116 y = varargin (2 : nargs);
120 elseif (rem (nargs, 2) == 1 && nargs ==
121 (2 * ndims (varargin{ceil (nargs / 2)})) + 1)
122 nv = ceil (nargs / 2);
126 x = varargin (1 : (nv - 1));
127 y = varargin ((nv + 1) : nargs);
129 error ("interpn: wrong number or incorrectly formatted input arguments");
132 if (any (! cellfun ("isvector", x)))
134 if (! size_equal (x{1}, x{i}) || ! size_equal (x{i}, v))
135 error ("interpn: dimensional mismatch");
139 x{i} = x{i}(idx{:})(:);
143 x{1} = x{1}(idx{:})(:);
146 method = tolower (method);
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{:});
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});
163 yidx{i} = lookup (x{i}, y{i}, "lr");
167 idx{i} = yidx{i} + (y{i} - x{i}(yidx{i})(:) >= x{i}(yidx{i} + 1)(:) - y{i});
169 vi = v (sub2ind (sz, idx{:}));
170 idx = zeros (prod (yshape), 1);
172 idx |= y{i} < min (x{i}(:)) | y{i} > max (x{i}(:));
175 vi = reshape (vi, yshape);
176 elseif (strcmp (method, "spline"))
177 if (any (! cellfun ("isvector", y)))
179 if (! size_equal (y{1}, y{i}))
180 error ("interpn: dimensional mismatch");
191 vi = __splinen__ (x, v, y, extrapval, "interpn");
193 if (size_equal (y{:}))
201 vi = vi (cellfun (@(x) sub2ind (size(vi), x{:}), idx));
202 vi = reshape (vi, size(y{1}));
204 elseif (strcmp (method, "cubic"))
205 error ("interpn: cubic interpolation not yet implemented");
207 error ("interpn: unrecognized interpolation METHOD");
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;
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;
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;
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;
251 %! f = @(x,y,z) x.^2 - y - z.^2;
252 %! [xx, yy, zz] = meshgrid (x, y, z);
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,:,:)));
261 %! [x,y,z] = ndgrid(0:2);
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)
271 %! [x, y, z] = ndgrid (0:2, 1:4, 2:6);
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)
284 %! [x, y, z] = ndgrid (xi, yi, zi);
286 %! fi = interpn (x, y, z, f, xi, yi, zi, "nearest");
287 %! assert (fi, x + y + z)
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)
294 %!test % for Matlab-compatible rounding for 'nearest'
295 %! X = meshgrid (1:4);
296 %! assert (interpn (X, 2.5, 2.5, 'nearest'), 3);
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;
305 %! zout(:,:,n) = [1 2 3 4 5;
309 %! 5 6 7 8 9] + (n-1);
312 %!assert (interpn (z), zout, tol)
313 %!assert (interpn (z, "linear"), zout, tol)
314 %!assert (interpn (z, "spline"), zout, tol)