|
7017
|
1 ## Copyright (C) 2000, 2006, 2007 Kai Habel |
|
5837
|
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. |
|
5837
|
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/>. |
|
5837
|
18 |
|
|
19 ## -*- texinfo -*- |
|
|
20 ## @deftypefn {Function File} {@var{zi}=} interp2 (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}) |
|
|
21 ## @deftypefnx {Function File} {@var{zi}=} interp2 (@var{Z}, @var{xi}, @var{yi}) |
|
|
22 ## @deftypefnx {Function File} {@var{zi}=} interp2 (@var{Z}, @var{n}) |
|
|
23 ## @deftypefnx {Function File} {@var{zi}=} interp2 (@dots{}, @var{method}) |
|
|
24 ## @deftypefnx {Function File} {@var{zi}=} interp2 (@dots{}, @var{method}, @var{extrapval}) |
|
|
25 ## |
|
|
26 ## Two-dimensional interpolation. @var{x}, @var{y} and @var{z} describe a |
|
|
27 ## surface function. If @var{x} and @var{y} are vectors their length |
|
6248
|
28 ## must correspondent to the size of @var{z}. @var{x} and @var{Yy} must be |
|
5837
|
29 ## monotonic. If they are matrices they must have the @code{meshgrid} |
|
|
30 ## format. |
|
|
31 ## |
|
|
32 ## @table @code |
|
|
33 ## @item interp2 (@var{x}, @var{y}, @var{Z}, @var{xi}, @var{yi}, @dots{}) |
|
|
34 ## Returns a matrix corresponding to the points described by the |
|
|
35 ## matrices @var{XI}, @var{YI}. |
|
|
36 ## |
|
|
37 ## If the last argument is a string, the interpolation method can |
|
|
38 ## be specified. The method can be 'linear', 'nearest' or 'cubic'. |
|
|
39 ## If it is omitted 'linear' interpolation is assumed. |
|
|
40 ## |
|
|
41 ## @item interp2 (@var{z}, @var{xi}, @var{yi}) |
|
|
42 ## Assumes @code{@var{x} = 1:rows (@var{z})} and @code{@var{y} = |
|
|
43 ## 1:columns (@var{z})} |
|
|
44 ## |
|
|
45 ## @item interp2 (@var{z}, @var{n}) |
|
7001
|
46 ## Interleaves the Matrix @var{z} n-times. If @var{n} is omitted a value |
|
5837
|
47 ## of @code{@var{n} = 1} is assumed. |
|
|
48 ## @end table |
|
|
49 ## |
|
|
50 ## The variable @var{method} defines the method to use for the |
|
|
51 ## interpolation. It can take one of the values |
|
|
52 ## |
|
|
53 ## @table @asis |
|
|
54 ## @item 'nearest' |
|
6218
|
55 ## Return the nearest neighbor. |
|
5837
|
56 ## @item 'linear' |
|
6218
|
57 ## Linear interpolation from nearest neighbors. |
|
5837
|
58 ## @item 'pchip' |
|
6218
|
59 ## Piece-wise cubic hermite interpolating polynomial (not implemented yet). |
|
5837
|
60 ## @item 'cubic' |
|
6218
|
61 ## Cubic interpolation from four nearest neighbors. |
|
5837
|
62 ## @item 'spline' |
|
|
63 ## Cubic spline interpolation--smooth first and second derivatives |
|
6702
|
64 ## throughout the curve. |
|
5837
|
65 ## @end table |
|
|
66 ## |
|
|
67 ## If a scalar value @var{extrapval} is defined as the final value, then |
|
|
68 ## values outside the mesh as set to this value. Note that in this case |
|
|
69 ## @var{method} must be defined as well. If @var{extrapval} is not |
|
6742
|
70 ## defined then NA is assumed. |
|
5837
|
71 ## |
|
|
72 ## @seealso{interp1} |
|
|
73 ## @end deftypefn |
|
|
74 |
|
|
75 ## Author: Kai Habel <kai.habel@gmx.de> |
|
|
76 ## 2005-03-02 Thomas Weber <weber@num.uni-sb.de> |
|
|
77 ## * Add test cases |
|
|
78 ## 2005-03-02 Paul Kienzle <pkienzle@users.sf.net> |
|
|
79 ## * Simplify |
|
|
80 ## 2005-04-23 Dmitri A. Sergatskov <dasergatskov@gmail.com> |
|
|
81 ## * Modified demo and test for new gnuplot interface |
|
|
82 ## 2005-09-07 Hoxide <hoxide_dirac@yahoo.com.cn> |
|
|
83 ## * Add bicubic interpolation method |
|
5838
|
84 ## * Fix the eat line bug when the last element of XI or YI is |
|
|
85 ## negative or zero. |
|
5837
|
86 ## 2005-11-26 Pierre Baldensperger <balden@libertysurf.fr> |
|
|
87 ## * Rather big modification (XI,YI no longer need to be |
|
|
88 ## "meshgridded") to be consistent with the help message |
|
|
89 ## above and for compatibility. |
|
|
90 |
|
|
91 function ZI = interp2 (varargin) |
|
5838
|
92 Z = X = Y = XI = YI = n = []; |
|
5837
|
93 method = "linear"; |
|
6742
|
94 extrapval = NA; |
|
5837
|
95 |
|
5838
|
96 switch (nargin) |
|
|
97 case 1 |
|
|
98 Z = varargin{1}; |
|
|
99 case 2 |
|
|
100 if (ischar (varargin{2})) |
|
|
101 [Z, method] = deal (varargin{:}); |
|
|
102 else |
|
|
103 [Z, n] = deal (varargin{:}); |
|
|
104 endif |
|
|
105 case 3 |
|
|
106 if (ischar (varargin{3})) |
|
|
107 [Z, n, method] = deal (varargin{:}); |
|
|
108 else |
|
|
109 [Z, XI, YI] = deal (varargin{:}); |
|
|
110 endif |
|
|
111 case 4 |
|
|
112 if (ischar (varargin{4})) |
|
|
113 [Z, XI, YI, method] = deal (varargin{:}); |
|
|
114 else |
|
|
115 [Z, n, method, extrapval] = deal (varargin{:}); |
|
|
116 endif |
|
|
117 case 5 |
|
|
118 if (ischar (varargin{4})) |
|
|
119 [Z, XI, YI, method, extrapval] = deal (varargin{:}); |
|
|
120 else |
|
|
121 [X, Y, Z, XI, YI] = deal (varargin{:}); |
|
|
122 endif |
|
|
123 case 6 |
|
|
124 [X, Y, Z, XI, YI, method] = deal (varargin{:}); |
|
|
125 case 7 |
|
|
126 [X, Y, Z, XI, YI, method, extrapval] = deal (varargin{:}); |
|
|
127 otherwise |
|
|
128 print_usage (); |
|
5837
|
129 endswitch |
|
|
130 |
|
|
131 ## Type checking. |
|
5838
|
132 if (!ismatrix (Z)) |
|
|
133 error ("interp2 expected matrix Z"); |
|
5837
|
134 endif |
|
5838
|
135 if (!isempty (n) && !isscalar (n)) |
|
|
136 error ("interp2 expected scalar n"); |
|
5837
|
137 endif |
|
5838
|
138 if (!ischar (method)) |
|
|
139 error ("interp2 expected string 'method'"); |
|
5837
|
140 endif |
|
5838
|
141 if (!isscalar (extrapval)) |
|
|
142 error ("interp2 expected n extrapval"); |
|
5837
|
143 endif |
|
|
144 |
|
5838
|
145 ## Define X, Y, XI, YI if needed |
|
5837
|
146 [zr, zc] = size (Z); |
|
5838
|
147 if (isempty (X)) |
|
|
148 X = 1:zc; |
|
|
149 Y = 1:zr; |
|
5837
|
150 endif |
|
5838
|
151 if (! isnumeric (X) || ! isnumeric (Y)) |
|
|
152 error ("interp2 expected numeric X, Y"); |
|
5837
|
153 endif |
|
5838
|
154 if (! isempty (n)) |
|
|
155 p = 2^n; |
|
|
156 XI = (p:p*zc)/p; |
|
|
157 YI = (p:p*zr)'/p; |
|
5837
|
158 endif |
|
5838
|
159 if (! isnumeric (XI) || ! isnumeric (YI)) |
|
|
160 error ("interp2 expected numeric XI, YI"); |
|
5837
|
161 endif |
|
|
162 |
|
6702
|
163 |
|
|
164 if (strcmp (method, "linear") || strcmp (method, "nearest")) |
|
|
165 |
|
|
166 ## If X and Y vectors produce a grid from them |
|
|
167 if (isvector (X) && isvector (Y)) |
|
|
168 [X, Y] = meshgrid (X, Y); |
|
|
169 elseif (! size_equal (X, Y)) |
|
|
170 error ("X and Y must be matrices of same size"); |
|
|
171 endif |
|
|
172 if (! size_equal (X, Z)) |
|
|
173 error ("X and Y size must match Z dimensions"); |
|
|
174 endif |
|
|
175 |
|
|
176 ## If Xi and Yi are vectors of different orientation build a grid |
|
|
177 if ((rows (XI) == 1 && columns (YI) == 1) |
|
|
178 || (columns (XI) == 1 && rows (YI) == 1)) |
|
|
179 [XI, YI] = meshgrid (XI, YI); |
|
|
180 elseif (! size_equal (XI, YI)) |
|
|
181 error ("XI and YI must be matrices of same size"); |
|
|
182 endif |
|
5837
|
183 |
|
6702
|
184 shape = size (XI); |
|
|
185 XI = reshape (XI, 1, prod (shape)); |
|
|
186 YI = reshape (YI, 1, prod (shape)); |
|
|
187 |
|
|
188 xidx = lookup (X(1, 2:end-1), XI) + 1; |
|
|
189 yidx = lookup (Y(2:end-1, 1), YI) + 1; |
|
5837
|
190 |
|
6702
|
191 if (strcmp (method, "linear")) |
|
|
192 ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy |
|
|
193 ## |
|
|
194 ## a-b |
|
|
195 ## | | |
|
|
196 ## c-d |
|
|
197 a = Z(1:(zr - 1), 1:(zc - 1)); |
|
|
198 b = Z(1:(zr - 1), 2:zc) - a; |
|
|
199 c = Z(2:zr, 1:(zc - 1)) - a; |
|
|
200 d = Z(2:zr, 2:zc) - a - b - c; |
|
5837
|
201 |
|
6702
|
202 idx = sub2ind (size (a), yidx, xidx); |
|
|
203 |
|
|
204 ## scale XI, YI values to a 1-spaced grid |
|
|
205 Xsc = (XI - X(1, xidx)) ./ (X(1, xidx + 1) - X(1, xidx)); |
|
|
206 Ysc = (YI - Y(yidx, 1)') ./ (Y(yidx + 1, 1) - Y(yidx, 1))'; |
|
|
207 |
|
|
208 ## apply plane equation |
|
|
209 ZI = a(idx) + b(idx).*Xsc + c(idx).*Ysc + d(idx).*Xsc.*Ysc; |
|
5837
|
210 |
|
6702
|
211 elseif (strcmp (method, "nearest")) |
|
|
212 xtable = X(1, :); |
|
|
213 ytable = Y(:, 1)'; |
|
|
214 ii = (XI - xtable(xidx) > xtable(xidx + 1) - XI); |
|
|
215 jj = (YI - ytable(yidx) > ytable(yidx + 1) - YI); |
|
|
216 idx = sub2ind (size (Z), yidx+jj, xidx+ii); |
|
|
217 ZI = Z(idx); |
|
|
218 endif |
|
5837
|
219 |
|
6742
|
220 ## set points outside the table to 'extrapval' |
|
6979
|
221 if (X (1, 1) < X (1, end)) |
|
|
222 if (Y (1, 1) < Y (end, 1)) |
|
|
223 ZI (XI < X(1,1) | XI > X(1,end) | YI < Y(1,1) | YI > Y(end,1)) = ... |
|
|
224 extrapval; |
|
|
225 else |
|
|
226 ZI (XI < X(1,1) | XI > X(1,end) | YI < Y(end,1) | YI > Y(1,1)) = ... |
|
|
227 extrapval; |
|
|
228 endif |
|
|
229 else |
|
|
230 if (Y (1, 1) < Y (end, 1)) |
|
|
231 ZI (XI < X(1,end) | XI > X(1,1) | YI < Y(1,1) | YI > Y(end,1)) = ... |
|
|
232 extrapval; |
|
|
233 else |
|
|
234 ZI (XI < X(1,end) | XI > X(1,1) | YI < Y(end,1) | YI > Y(1,1)) = ... |
|
|
235 extrapval; |
|
|
236 endif |
|
|
237 endif |
|
|
238 |
|
6702
|
239 ZI = reshape (ZI, shape); |
|
|
240 else |
|
5837
|
241 |
|
6702
|
242 ## If X and Y vectors produce a grid from them |
|
|
243 if (isvector (X) && isvector (Y)) |
|
|
244 X = X(:).'; |
|
|
245 Y = Y(:); |
|
|
246 if (!isequal ([length(X), length(Y)], size(Z))) |
|
|
247 error ("X and Y size must match Z dimensions"); |
|
|
248 endif |
|
|
249 elseif (!size_equal (X, Y)) |
|
|
250 error ("X and Y must be matrices of same size"); |
|
|
251 if (! size_equal (X, Z)) |
|
|
252 error ("X and Y size must match Z dimensions"); |
|
|
253 endif |
|
|
254 endif |
|
5837
|
255 |
|
6702
|
256 ## If Xi and Yi are vectors of different orientation build a grid |
|
|
257 if ((rows (XI) == 1 && columns (YI) == 1) |
|
|
258 || (columns (XI) == 1 && rows (YI) == 1)) |
|
|
259 ## Do nothing |
|
|
260 elseif (! size_equal (XI, YI)) |
|
|
261 error ("XI and YI must be matrices of same size"); |
|
|
262 endif |
|
5837
|
263 |
|
6702
|
264 ## FIXME bicubic/__splinen__ don't handle arbitrary XI, YI |
|
|
265 if (strcmp (method, "cubic")) |
|
|
266 ZI = bicubic (X, Y, Z, XI(1,:), YI(:,1), extrapval); |
|
5837
|
267 |
|
6702
|
268 elseif (strcmp (method, "spline")) |
|
|
269 ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrapval, |
|
|
270 "spline"); |
|
|
271 else |
|
|
272 error ("interpolation method not recognized"); |
|
|
273 endif |
|
5837
|
274 |
|
6702
|
275 endif |
|
5837
|
276 endfunction |
|
|
277 |
|
|
278 %!demo |
|
|
279 %! A=[13,-1,12;5,4,3;1,6,2]; |
|
|
280 %! x=[0,1,4]; y=[10,11,12]; |
|
|
281 %! xi=linspace(min(x),max(x),17); |
|
|
282 %! yi=linspace(min(y),max(y),26)'; |
|
|
283 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'linear')); |
|
|
284 %! [x,y] = meshgrid(x,y); |
|
|
285 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
|
|
286 |
|
|
287 %!demo |
|
|
288 %! A=[13,-1,12;5,4,3;1,6,2]; |
|
|
289 %! x=[0,1,4]; y=[10,11,12]; |
|
|
290 %! xi=linspace(min(x),max(x),17); |
|
|
291 %! yi=linspace(min(y),max(y),26)'; |
|
|
292 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'nearest')); |
|
|
293 %! [x,y] = meshgrid(x,y); |
|
|
294 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
|
|
295 |
|
6702
|
296 %!demo |
|
5837
|
297 %! A=[13,-1,12;5,4,3;1,6,2]; |
|
|
298 %! x=[0,1,2]; y=[10,11,12]; |
|
|
299 %! xi=linspace(min(x),max(x),17); |
|
6702
|
300 %! yi=linspace(min(y),max(y),26)'; |
|
5837
|
301 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic')); |
|
|
302 %! [x,y] = meshgrid(x,y); |
|
|
303 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
|
|
304 |
|
6702
|
305 %!demo |
|
|
306 %! A=[13,-1,12;5,4,3;1,6,2]; |
|
|
307 %! x=[0,1,2]; y=[10,11,12]; |
|
|
308 %! xi=linspace(min(x),max(x),17); |
|
|
309 %! yi=linspace(min(y),max(y),26)'; |
|
|
310 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'spline')); |
|
|
311 %! [x,y] = meshgrid(x,y); |
|
|
312 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
|
|
313 |
|
5837
|
314 %!test % simple test |
|
|
315 %! x = [1,2,3]; |
|
|
316 %! y = [4,5,6,7]; |
|
|
317 %! [X, Y] = meshgrid(x,y); |
|
|
318 %! Orig = X.^2 + Y.^3; |
|
|
319 %! xi = [1.2,2, 1.5]; |
|
|
320 %! yi = [6.2, 4.0, 5.0]'; |
|
|
321 %! |
|
|
322 %! Expected = ... |
|
|
323 %! [243, 245.4, 243.9; |
|
|
324 %! 65.6, 68, 66.5; |
|
|
325 %! 126.6, 129, 127.5]; |
|
|
326 %! Result = interp2(x,y,Orig, xi, yi); |
|
|
327 %! |
|
|
328 %! assert(Result, Expected, 1000*eps); |
|
|
329 |
|
|
330 %!test % 2^n form |
|
|
331 %! x = [1,2,3]; |
|
|
332 %! y = [4,5,6,7]; |
|
|
333 %! [X, Y] = meshgrid(x,y); |
|
|
334 %! Orig = X.^2 + Y.^3; |
|
|
335 %! xi = [1:0.25:3]; yi = [4:0.25:7]'; |
|
|
336 %! Expected = interp2(x,y,Orig, xi, yi); |
|
|
337 %! Result = interp2(Orig,2); |
|
|
338 %! |
|
|
339 %! assert(Result, Expected, 10*eps); |
|
|
340 |
|
|
341 %!test % matrix slice |
|
|
342 %! A = eye(4); |
|
|
343 %! assert(interp2(A,[1:4],[1:4]),[1,1,1,1]); |
|
|
344 |
|
|
345 %!test % non-gridded XI,YI |
|
|
346 %! A = eye(4); |
|
|
347 %! assert(interp2(A,[1,2;3,4],[1,3;2,4]),[1,0;0,1]); |
|
|
348 |
|
|
349 %!test % for values outside of boundaries |
|
|
350 %! x = [1,2,3]; |
|
|
351 %! y = [4,5,6,7]; |
|
|
352 %! [X, Y] = meshgrid(x,y); |
|
|
353 %! Orig = X.^2 + Y.^3; |
|
|
354 %! xi = [0,4]; |
|
|
355 %! yi = [3,8]'; |
|
6742
|
356 %! assert(interp2(x,y,Orig, xi, yi),[NA,NA;NA,NA]); |
|
5837
|
357 %! assert(interp2(x,y,Orig, xi, yi,'linear', 0),[0,0;0,0]); |
|
|
358 |
|
|
359 %!test % for values at boundaries |
|
|
360 %! A=[1,2;3,4]; |
|
|
361 %! x=[0,1]; |
|
|
362 %! y=[2,3]'; |
|
|
363 %! assert(interp2(x,y,A,x,y,'linear'), A); |
|
|
364 %! assert(interp2(x,y,A,x,y,'nearest'), A); |
|
|
365 |
