|
6702
|
1 ## Copyright (C) 2007 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 2, or (at your option) |
|
|
8 ## 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, write to the Free |
|
|
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
18 ## 02110-1301, USA. |
|
|
19 |
|
|
20 ## -*- texinfo -*- |
|
|
21 ## @deftypefn {Function File} {@var{vi} =} interpn (@var{x1}, @var{x2}, @dots{}, @var{v}, @var{y1}, @var{y2}, @dots{}) |
|
|
22 ## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}, @var{y1}, @var{y2}, @dots{}) |
|
|
23 ## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}, @var{m}) |
|
|
24 ## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}) |
|
|
25 ## @deftypefnx {Function File} {@var{vi} =} interpn (@dots{}, @var{method}) |
|
|
26 ## @deftypefnx {Function File} {@var{vi} =} interpn (@dots{}, @var{method}, @var{extrapval}) |
|
|
27 ## |
|
|
28 ## Perform @var{n}-dimensional interpolation, where @var{n} is at least two. |
|
|
29 ## Each element of then @var{n}-dimensional array @var{v} represents a value |
|
|
30 ## at a location given by the parameters @var{x1}, @var{x2}, @dots{}, @var{xn}. |
|
|
31 ## The parameters @var{x1}, @var{x2}, @dots{}, @var{xn} are either |
|
|
32 ## @var{n}-dimensional arrays of the same size as the array @var{v} in |
|
|
33 ## the 'ndgrid' format or vectors. The parameters @var{y1}, etc respect a |
|
|
34 ## similar format to @var{x1}, etc, and they represent the points at which |
|
|
35 ## the array @var{vi} is interpolated. |
|
|
36 ## |
|
|
37 ## If @var{x1}, @dots{}, @var{xn} are ommitted, they are assumed to be |
|
|
38 ## @code{x1 = 1 : size (@var{v}, 1)}, etc. If @var{m} is specified, then |
|
|
39 ## the interpolation adds a point half way between each of the interplation |
|
|
40 ## points. This process is performed @var{m} times. If only @var{v} is |
|
|
41 ## specified, then @var{m} is assumed to be @code{1}. |
|
|
42 ## |
|
|
43 ## Method is one of: |
|
|
44 ## |
|
|
45 ## @table @asis |
|
|
46 ## @item 'nearest' |
|
|
47 ## Return the nearest neighbour. |
|
|
48 ## @item 'linear' |
|
|
49 ## Linear interpolation from nearest neighbours. |
|
|
50 ## @item 'cubic' |
|
|
51 ## Cubic interpolation from four nearest neighbours (not implemented yet). |
|
|
52 ## @item 'spline' |
|
|
53 ## Cubic spline interpolation--smooth first and second derivatives |
|
|
54 ## throughout the curve. |
|
|
55 ## @end table |
|
|
56 ## |
|
|
57 ## The default method is 'linear'. |
|
|
58 ## |
|
|
59 ## If @var{extrap} is the string 'extrap', then extrapolate values beyond |
|
|
60 ## the endpoints. If @var{extrap} is a number, replace values beyond the |
|
6742
|
61 ## endpoints with that number. If @var{extrap} is missing, assume NA. |
|
6702
|
62 ## @seealso{interp1, interp2, spline, ndgrid} |
|
|
63 ## @end deftypefn |
|
|
64 |
|
|
65 function vi = interpn (varargin) |
|
|
66 |
|
|
67 method = "linear"; |
|
6742
|
68 extrapval = NA; |
|
6702
|
69 nargs = nargin; |
|
|
70 |
|
|
71 if (nargin < 1) |
|
|
72 print_usage (); |
|
|
73 endif |
|
|
74 |
|
|
75 if (ischar (varargin {end})) |
|
|
76 method = varargin {end}; |
|
|
77 nargs = nargs - 1; |
|
|
78 elseif (ischar (varargin {end - 1})) |
|
6722
|
79 if (! isnumeric (varargin {end}) || ! isscalar (varargin {end})) |
|
6702
|
80 error ("extrapal is expected to be a numeric scalar"); |
|
|
81 endif |
|
|
82 method = varargin {end - 1}; |
|
|
83 nargs = nargs - 2; |
|
|
84 endif |
|
|
85 |
|
|
86 if (nargs < 3) |
|
|
87 v = varargin {1}; |
|
|
88 m = 1; |
|
|
89 if (nargs == 2) |
|
|
90 m = varargin {2}; |
|
|
91 if (! isnumeric (m) || ! isscalar (m) || floor (m) != m) |
|
|
92 error ("m is expected to be a integer scalar"); |
|
|
93 endif |
|
|
94 endif |
|
|
95 sz = size (v); |
|
|
96 nd = ndims (v); |
|
|
97 x = cell (1, nd); |
|
|
98 y = cell (1, nd); |
|
|
99 for i = 1 : nd; |
|
|
100 x{i} = 1 : sz(i); |
|
|
101 y{i} = 1 : (1 / (2 ^ m)) : sz(i); |
|
|
102 endfor |
|
|
103 elseif (! isvector (varargin {1}) && nargs == (ndims (varargin {1}) + 1)) |
|
|
104 v = varargin {1}; |
|
|
105 sz = size (v); |
|
|
106 nd = ndims (v); |
|
|
107 x = cell (1, nd); |
|
|
108 y = varargin (2 : nargs); |
|
|
109 for i = 1 : nd; |
|
|
110 x{i} = 1 : sz(i); |
|
|
111 endfor |
|
|
112 elseif (rem (nargs, 2) == 1 && nargs == |
|
|
113 (2 * ndims (varargin {ceil (nargs / 2)})) + 1) |
|
|
114 nv = ceil (nargs / 2); |
|
|
115 v = varargin {nv}; |
|
|
116 sz = size (v); |
|
|
117 nd = ndims (v); |
|
|
118 x = varargin (1 : (nv - 1)); |
|
|
119 y = varargin ((nv + 1) : nargs); |
|
|
120 else |
|
|
121 error ("wrong number or incorrectly formatted input arguments"); |
|
|
122 endif |
|
|
123 |
|
|
124 if (any (! cellfun (@isvector, x))) |
|
|
125 for i = 2 : nd |
|
|
126 if (! size_equal (x{1}, x{i}) || ! size_equal (x{i}, v)) |
|
|
127 error ("dimensional mismatch"); |
|
|
128 endif |
|
|
129 idx (1 : nd) = {1}; |
|
|
130 idx (i) = ":"; |
|
6721
|
131 x{i} = x{i}(idx{:})(:); |
|
6702
|
132 endfor |
|
|
133 idx (1 : nd) = {1}; |
|
|
134 idx (1) = ":"; |
|
6721
|
135 x{1} = x{1}(idx{:})(:); |
|
6702
|
136 endif |
|
|
137 |
|
|
138 if (strcmp (method, "linear") || strcmp (method, "nearest")) |
|
|
139 if (all (cellfun (@isvector, y))) |
|
|
140 [y{:}] = ndgrid (y{:}); |
|
|
141 endif |
|
6721
|
142 elseif (any (! cellfun (@isvector, y))) |
|
6702
|
143 for i = 1 : nd |
|
|
144 idx (1 : nd) = {1}; |
|
|
145 idx (i) = ":"; |
|
6721
|
146 y{i} = y{i}(idx{:})(:).'; |
|
6702
|
147 endfor |
|
|
148 endif |
|
|
149 |
|
|
150 method = tolower (method); |
|
|
151 if (strcmp (method, "linear")) |
|
|
152 vi = __lin_interpn__ (x{:}, v, y{:}); |
|
6742
|
153 vi (isna (vi)) = extrapval; |
|
6702
|
154 elseif (strcmp (method, "nearest")) |
|
|
155 yshape = size (y{1}); |
|
|
156 yidx = cell (1, nd); |
|
|
157 for i = 1 : nd |
|
|
158 y{i} = y{i}(:); |
|
|
159 yidx{i} = lookup (x{i}(2:end-1), y{i}) + 1; |
|
|
160 endfor |
|
|
161 idx = cell (1,nd); |
|
|
162 for i = 1 : nd |
|
|
163 idx {i} = yidx{i} + (y{i} - x{i}(yidx{i}).' > ... |
|
|
164 x{i}(yidx{i} + 1).' - y{i}); |
|
|
165 endfor |
|
|
166 vi = v (sub2ind (sz, idx{:})); |
|
|
167 idx = zeros (prod(yshape),1); |
|
|
168 for i = 1 : nd |
|
|
169 idx |= y{i} < min (x{i}(:)) | y{i} > max (x{i}(:)); |
|
|
170 endfor |
|
|
171 vi(idx) = extrapval; |
|
|
172 vi = reshape (vi, yshape); |
|
6721
|
173 elseif (strcmp (method, "spline")) |
|
6702
|
174 vi = __splinen__ (x, v, y, extrapval, "interpn"); |
|
|
175 elseif (strcmp (method, "cubic")) |
|
|
176 error ("cubic interpolation not yet implemented"); |
|
|
177 else |
|
|
178 error ("unrecognized interpolation method"); |
|
|
179 endif |
|
|
180 |
|
|
181 endfunction |
|
|
182 |
|
|
183 %!demo |
|
|
184 %! A=[13,-1,12;5,4,3;1,6,2]; |
|
|
185 %! x=[0,1,4]; y=[10,11,12]; |
|
|
186 %! xi=linspace(min(x),max(x),17); |
|
|
187 %! yi=linspace(min(y),max(y),26)'; |
|
|
188 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"linear").'); |
|
|
189 %! [x,y] = meshgrid(x,y); |
|
|
190 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
|
|
191 |
|
|
192 %!demo |
|
|
193 %! A=[13,-1,12;5,4,3;1,6,2]; |
|
|
194 %! x=[0,1,4]; y=[10,11,12]; |
|
|
195 %! xi=linspace(min(x),max(x),17); |
|
|
196 %! yi=linspace(min(y),max(y),26)'; |
|
|
197 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"nearest").'); |
|
|
198 %! [x,y] = meshgrid(x,y); |
|
|
199 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
|
|
200 |
|
|
201 %!#demo |
|
|
202 %! A=[13,-1,12;5,4,3;1,6,2]; |
|
|
203 %! x=[0,1,2]; y=[10,11,12]; |
|
|
204 %! xi=linspace(min(x),max(x),17); |
|
|
205 %! yi=linspace(min(y),max(y),26)'; |
|
|
206 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"cubic").'); |
|
|
207 %! [x,y] = meshgrid(x,y); |
|
|
208 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
|
|
209 |
|
|
210 %!demo |
|
|
211 %! A=[13,-1,12;5,4,3;1,6,2]; |
|
|
212 %! x=[0,1,2]; y=[10,11,12]; |
|
|
213 %! xi=linspace(min(x),max(x),17); |
|
|
214 %! yi=linspace(min(y),max(y),26)'; |
|
|
215 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"spline").'); |
|
|
216 %! [x,y] = meshgrid(x,y); |
|
|
217 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
|
|
218 |
|
6721
|
219 |
|
|
220 %!demo |
|
|
221 %! x = y = z = -1:1; |
|
|
222 %! f = @(x,y,z) x.^2 - y - z.^2; |
|
|
223 %! [xx, yy, zz] = meshgrid (x, y, z); |
|
|
224 %! v = f (xx,yy,zz); |
|
|
225 %! xi = yi = zi = -1:0.1:1; |
|
|
226 %! [xxi, yyi, zzi] = ndgrid (xi, yi, zi); |
|
|
227 %! vi = interpn(x, y, z, v, xxi, yyi, zzi, 'spline'); |
|
|
228 %! mesh (yi, zi, squeeze (vi(1,:,:))); |
|
|
229 |
