1
## Copyright (C) 2007-2013 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} =} 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})
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 @qcode{"meshgrid"} format or vectors. The parameters
32
## @var{xi}, etc. respect a similar format to @var{x}, etc., and they
33
## represent the points at which the array @var{vi} is interpolated.
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}.
45
## @item @qcode{"nearest"}
46
## Return the nearest neighbor.
48
## @item @qcode{"linear"}
49
## Linear interpolation from nearest neighbors.
51
## @item @qcode{"cubic"}
52
## Cubic interpolation from four nearest neighbors (not implemented yet).
54
## @item @qcode{"spline"}
55
## Cubic spline interpolation---smooth first and second derivatives
56
## throughout the curve.
59
## The default method is @qcode{"linear"}.
61
## If @var{extrap} is the string @qcode{"extrap"}, then extrapolate values
62
## beyond the endpoints. If @var{extrap} is a number, replace values beyond
63
## the endpoints with that number. If @var{extrap} is missing, assume NA.
64
## @seealso{interp1, interp2, spline, meshgrid}
67
function vi = interp3 (varargin)
72
if (nargin < 1 || ! isnumeric (varargin{1}))
76
if (ischar (varargin{end}))
77
method = varargin{end};
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");
83
extrapval = varargin{end};
84
method = varargin{end-1};
88
if (nargs < 3 || (nargs == 4 && ! isvector (varargin{1})
89
&& nargs == (ndims (varargin{1}) + 1)))
92
error ("interp3: expect 3-dimensional array of values");
94
x = varargin (2:nargs);
95
if (any (! cellfun (@isvector, x)))
97
if (! size_equal (x{1}, x{i}))
98
error ("interp3: dimensional mismatch");
100
x{i} = permute (x{i}, [2, 1, 3]);
102
x{1} = permute (x{1}, [2, 1, 3]);
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)
109
error ("interp3: expect 3-dimensional array of values");
112
if (any (! cellfun (@isvector, x)))
114
if (! size_equal (x{1}, x{i}) || ! size_equal (x{i}, v))
115
error ("interp3: dimensional mismatch");
117
x{i} = permute (x{i}, [2, 1, 3]);
119
x{1} = permute (x{1}, [2, 1, 3]);
122
if (any (! cellfun (@isvector, y)))
124
if (! size_equal (y{1}, y{i}))
125
error ("interp3: dimensional mismatch");
127
y{i} = permute (y{i}, [2, 1, 3]);
129
y{1} = permute (y{1}, [2, 1, 3]);
131
v = permute (v, [2, 1, 3]);
132
vi = ipermute (interpn (x{:}, v, y{:}, method, extrapval), [2, 1, 3]);
134
error ("interp3: wrong number or incorrectly formatted input arguments");
140
%! x = y = z = -1:1; y = y + 2;
141
%! f = @(x,y,z) x.^2 - y - z.^2;
142
%! [xx, yy, zz] = meshgrid (x, y, z);
144
%! xi = yi = zi = -1:0.5:1; yi = yi + 2.1;
145
%! [xxi, yyi, zzi] = meshgrid (xi, yi, zi);
146
%! vi = interp3 (x, y, z, v, xxi, yyi, zzi);
147
%! [xxi, yyi, zzi] = ndgrid (yi, xi, zi);
148
%! vi2 = interpn (y, x, z, v, xxi, yyi, zzi);
150
%! assert (vi, vi2, tol);
153
%! x=z=1:2; y=1:3;xi=zi=.6:1.6; yi=1; v=ones([3,2,2]); v(:,2,1)=[7 ;5;4]; v(:,1,2)=[2 ;3;5];
154
%! [xxi3, yyi3, zzi3] = meshgrid (xi, yi, zi);
155
%! [xxi, yyi, zzi] = ndgrid (yi, xi, zi);
156
%! vi = interp3 (x, y, z, v, xxi3, yyi3, zzi3, "nearest");
157
%! vi2 = interpn (y, x, z, v, xxi, yyi, zzi,"nearest");
161
%! x=z=1:2; y=1:3;xi=zi=.6:1.6; yi=1; v=ones([3,2,2]); v(:,2,1)=[7 ;5;4]; v(:,1,2)=[2 ;3;5];
162
%! vi = interp3 (x, y, z, v, xi+1, yi, zi, "nearest",3);
163
%! vi2 = interpn (y, x, z, v, yi, xi+1, zi,"nearest", 3);
167
%! x=z=1:2; y=1:3;xi=zi=.6:1.6; yi=1; v=ones([3,2,2]); v(:,2,1)=[7 ;5;4]; v(:,1,2)=[2 ;3;5];
168
%! vi = interp3 (x, y, z, v, xi, yi, zi, "nearest");
169
%! vi2 = interpn (y, x, z, v, yi, xi, zi,"nearest");
173
%! x=z=1:2; y=1:3;xi=zi=.6:1.6; yi=1; v=ones([3,2,2]); v(:,2,1)=[7 ;5;4]; v(:,1,2)=[2 ;3;5];
174
%! vi = interp3 (v, xi, yi, zi, "nearest",3);
175
%! vi2 = interpn (v, yi, xi, zi,"nearest", 3);
179
%! xi=zi=.6:1.6; yi=1; v=ones([3,2,2]); v(:,2,1)=[7 ;5;4]; v(:,1,2)=[2 ;3;5];
180
%! vi = interp3 (v, xi, yi, zi, "nearest");
181
%! vi2 = interpn (v, yi, xi, zi,"nearest");
184
%!shared z, zout, tol
185
%! z = zeros (3, 3, 3);
186
%! zout = zeros (5, 5, 5);
187
%! z(:,:,1) = [1 3 5; 3 5 7; 5 7 9];
188
%! z(:,:,2) = z(:,:,1) + 2;
189
%! z(:,:,3) = z(:,:,2) + 2;
191
%! zout(:,:,n) = [1 2 3 4 5;
195
%! 5 6 7 8 9] + (n-1);
198
%!assert (interp3 (z), zout, tol)
199
%!assert (interp3 (z, "linear"), zout, tol)
200
%!assert (interp3 (z, "spline"), zout, tol)