|
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 |
|
7016
|
7 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
8 ## your option) any later version. |
|
6702
|
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 |
|
7016
|
16 ## along with Octave; see the file COPYING. If not, see |
|
|
17 ## <http://www.gnu.org/licenses/>. |
|
6702
|
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 then @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 ## |
|
7001
|
36 ## If @var{x1}, @dots{}, @var{xn} are omitted, they are assumed to be |
|
6702
|
37 ## @code{x1 = 1 : size (@var{v}, 1)}, etc. If @var{m} is specified, then |
|
7001
|
38 ## the interpolation adds a point half way between each of the interpolation |
|
6702
|
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 neighbour. |
|
|
47 ## @item 'linear' |
|
|
48 ## Linear interpolation from nearest neighbours. |
|
|
49 ## @item 'cubic' |
|
|
50 ## Cubic interpolation from four nearest neighbours (not implemented yet). |
|
|
51 ## @item 'spline' |
|
|
52 ## Cubic spline interpolation--smooth first and second derivatives |
|
|
53 ## throughout the curve. |
|
|
54 ## @end table |
|
|
55 ## |
|
|
56 ## The default method is 'linear'. |
|
|
57 ## |
|
|
58 ## If @var{extrap} is the string 'extrap', then extrapolate values beyond |
|
|
59 ## the endpoints. If @var{extrap} is a number, replace values beyond the |
|
6742
|
60 ## endpoints with that number. If @var{extrap} is missing, assume NA. |
|
6702
|
61 ## @seealso{interp1, interp2, spline, ndgrid} |
|
|
62 ## @end deftypefn |
|
|
63 |
|
|
64 function vi = interpn (varargin) |
|
|
65 |
|
|
66 method = "linear"; |
|
6742
|
67 extrapval = NA; |
|
6702
|
68 nargs = nargin; |
|
|
69 |
|
|
70 if (nargin < 1) |
|
|
71 print_usage (); |
|
|
72 endif |
|
|
73 |
|
7208
|
74 if (ischar (varargin{end})) |
|
|
75 method = varargin{end}; |
|
6702
|
76 nargs = nargs - 1; |
|
7208
|
77 elseif (ischar (varargin{end - 1})) |
|
|
78 if (! isnumeric (varargin{end}) || ! isscalar (varargin{end})) |
|
6702
|
79 error ("extrapal is expected to be a numeric scalar"); |
|
|
80 endif |
|
7208
|
81 method = varargin{end - 1}; |
|
6702
|
82 nargs = nargs - 2; |
|
|
83 endif |
|
|
84 |
|
|
85 if (nargs < 3) |
|
7208
|
86 v = varargin{1}; |
|
6702
|
87 m = 1; |
|
|
88 if (nargs == 2) |
|
7208
|
89 m = varargin{2}; |
|
6702
|
90 if (! isnumeric (m) || ! isscalar (m) || floor (m) != m) |
|
|
91 error ("m is expected to be a integer scalar"); |
|
|
92 endif |
|
|
93 endif |
|
|
94 sz = size (v); |
|
|
95 nd = ndims (v); |
|
|
96 x = cell (1, nd); |
|
|
97 y = cell (1, nd); |
|
|
98 for i = 1 : nd; |
|
|
99 x{i} = 1 : sz(i); |
|
|
100 y{i} = 1 : (1 / (2 ^ m)) : sz(i); |
|
|
101 endfor |
|
7208
|
102 elseif (! isvector (varargin{1}) && nargs == (ndims (varargin{1}) + 1)) |
|
|
103 v = varargin{1}; |
|
6702
|
104 sz = size (v); |
|
|
105 nd = ndims (v); |
|
|
106 x = cell (1, nd); |
|
|
107 y = varargin (2 : nargs); |
|
|
108 for i = 1 : nd; |
|
|
109 x{i} = 1 : sz(i); |
|
|
110 endfor |
|
|
111 elseif (rem (nargs, 2) == 1 && nargs == |
|
7208
|
112 (2 * ndims (varargin{ceil (nargs / 2)})) + 1) |
|
6702
|
113 nv = ceil (nargs / 2); |
|
7208
|
114 v = varargin{nv}; |
|
6702
|
115 sz = size (v); |
|
|
116 nd = ndims (v); |
|
|
117 x = varargin (1 : (nv - 1)); |
|
|
118 y = varargin ((nv + 1) : nargs); |
|
|
119 else |
|
|
120 error ("wrong number or incorrectly formatted input arguments"); |
|
|
121 endif |
|
|
122 |
|
|
123 if (any (! cellfun (@isvector, x))) |
|
|
124 for i = 2 : nd |
|
|
125 if (! size_equal (x{1}, x{i}) || ! size_equal (x{i}, v)) |
|
|
126 error ("dimensional mismatch"); |
|
|
127 endif |
|
|
128 idx (1 : nd) = {1}; |
|
|
129 idx (i) = ":"; |
|
6721
|
130 x{i} = x{i}(idx{:})(:); |
|
6702
|
131 endfor |
|
|
132 idx (1 : nd) = {1}; |
|
|
133 idx (1) = ":"; |
|
6721
|
134 x{1} = x{1}(idx{:})(:); |
|
6702
|
135 endif |
|
|
136 |
|
7421
|
137 method = tolower (method); |
|
6702
|
138 |
|
|
139 if (strcmp (method, "linear")) |
|
|
140 vi = __lin_interpn__ (x{:}, v, y{:}); |
|
6742
|
141 vi (isna (vi)) = extrapval; |
|
6702
|
142 elseif (strcmp (method, "nearest")) |
|
|
143 yshape = size (y{1}); |
|
|
144 yidx = cell (1, nd); |
|
|
145 for i = 1 : nd |
|
|
146 y{i} = y{i}(:); |
|
|
147 yidx{i} = lookup (x{i}(2:end-1), y{i}) + 1; |
|
|
148 endfor |
|
|
149 idx = cell (1,nd); |
|
|
150 for i = 1 : nd |
|
7421
|
151 idx{i} = yidx{i} + (y{i} - x{i}(yidx{i}) > x{i}(yidx{i} + 1) - y{i}); |
|
6702
|
152 endfor |
|
|
153 vi = v (sub2ind (sz, idx{:})); |
|
|
154 idx = zeros (prod(yshape),1); |
|
|
155 for i = 1 : nd |
|
|
156 idx |= y{i} < min (x{i}(:)) | y{i} > max (x{i}(:)); |
|
|
157 endfor |
|
|
158 vi(idx) = extrapval; |
|
|
159 vi = reshape (vi, yshape); |
|
6721
|
160 elseif (strcmp (method, "spline")) |
|
7421
|
161 szi = size (y{1}); |
|
|
162 for i = 1 : nd |
|
|
163 y{i} = y{i}(:); |
|
|
164 endfor |
|
|
165 |
|
6702
|
166 vi = __splinen__ (x, v, y, extrapval, "interpn"); |
|
7421
|
167 |
|
|
168 ## get all diagonal elements of vi |
|
|
169 sc = sum ([1 cumprod(size (vi)(1:end-1))]); |
|
|
170 vi = vi(sc * [0:size(vi,1)-1] + 1); |
|
|
171 vi = reshape (vi,szi); |
|
6702
|
172 elseif (strcmp (method, "cubic")) |
|
|
173 error ("cubic interpolation not yet implemented"); |
|
|
174 else |
|
|
175 error ("unrecognized interpolation method"); |
|
|
176 endif |
|
|
177 |
|
|
178 endfunction |
|
|
179 |
|
|
180 %!demo |
|
|
181 %! A=[13,-1,12;5,4,3;1,6,2]; |
|
|
182 %! x=[0,1,4]; y=[10,11,12]; |
|
|
183 %! xi=linspace(min(x),max(x),17); |
|
7421
|
184 %! AI=linspace(min(y),max(y),26)'; |
|
6702
|
185 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"linear").'); |
|
|
186 %! [x,y] = meshgrid(x,y); |
|
|
187 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
|
|
188 |
|
|
189 %!demo |
|
|
190 %! A=[13,-1,12;5,4,3;1,6,2]; |
|
|
191 %! x=[0,1,4]; y=[10,11,12]; |
|
|
192 %! xi=linspace(min(x),max(x),17); |
|
|
193 %! yi=linspace(min(y),max(y),26)'; |
|
|
194 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"nearest").'); |
|
|
195 %! [x,y] = meshgrid(x,y); |
|
|
196 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
|
|
197 |
|
|
198 %!#demo |
|
|
199 %! A=[13,-1,12;5,4,3;1,6,2]; |
|
|
200 %! x=[0,1,2]; y=[10,11,12]; |
|
|
201 %! xi=linspace(min(x),max(x),17); |
|
|
202 %! yi=linspace(min(y),max(y),26)'; |
|
|
203 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"cubic").'); |
|
|
204 %! [x,y] = meshgrid(x,y); |
|
|
205 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
|
|
206 |
|
|
207 %!demo |
|
|
208 %! A=[13,-1,12;5,4,3;1,6,2]; |
|
|
209 %! x=[0,1,2]; y=[10,11,12]; |
|
|
210 %! xi=linspace(min(x),max(x),17); |
|
|
211 %! yi=linspace(min(y),max(y),26)'; |
|
|
212 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"spline").'); |
|
|
213 %! [x,y] = meshgrid(x,y); |
|
|
214 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
|
|
215 |
|
6721
|
216 |
|
|
217 %!demo |
|
|
218 %! x = y = z = -1:1; |
|
|
219 %! f = @(x,y,z) x.^2 - y - z.^2; |
|
|
220 %! [xx, yy, zz] = meshgrid (x, y, z); |
|
|
221 %! v = f (xx,yy,zz); |
|
|
222 %! xi = yi = zi = -1:0.1:1; |
|
|
223 %! [xxi, yyi, zzi] = ndgrid (xi, yi, zi); |
|
|
224 %! vi = interpn(x, y, z, v, xxi, yyi, zzi, 'spline'); |
|
|
225 %! mesh (yi, zi, squeeze (vi(1,:,:))); |
|
|
226 |
|
7421
|
227 |
|
|
228 %!test |
|
|
229 %! [x,y,z] = ndgrid(0:2); |
|
|
230 %! f = x+y+z; |
|
|
231 %! assert (interpn(x,y,z,f,[.5 1.5],[.5 1.5],[.5 1.5]), [1.5, 4.5]) |
|
|
232 %! assert (interpn(x,y,z,f,[.51 1.51],[.51 1.51],[.51 1.51],'nearest'), [3, 6]) |
|
|
233 %! assert (interpn(x,y,z,f,[.5 1.5],[.5 1.5],[.5 1.5],'spline'), [1.5, 4.5]) |
|
|
234 %! assert (interpn(x,y,z,f,x,y,z), f) |
|
|
235 %! assert (interpn(x,y,z,f,x,y,z,'nearest'), f) |
|
|
236 %! assert (interpn(x,y,z,f,x,y,z,'spline'), f) |
