Mercurial > octave-nkf

comparison scripts/general/interp2.m @ 10549:95c3e38098bf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Untabify .m scripts
author Rik <code@nomad.inbox5.com>
date Fri, 23 Apr 2010 11:28:50 -0700
parents 62bb59f927b1
children d1978e7364ad
comparison
equal deleted inserted replaced
10548:479536c5bb10 10549:95c3e38098bf
71 ## defined then NA is assumed. 71 ## defined then NA is assumed.
72 ## 72 ##
73 ## @seealso{interp1} 73 ## @seealso{interp1}
74 ## @end deftypefn 74 ## @end deftypefn
75 75
76 ## Author: Kai Habel <kai.habel@gmx.de> 76 ## Author: Kai Habel <kai.habel@gmx.de>
77 ## 2005-03-02 Thomas Weber <weber@num.uni-sb.de> 77 ## 2005-03-02 Thomas Weber <weber@num.uni-sb.de>
78 ## * Add test cases 78 ## * Add test cases
79 ## 2005-03-02 Paul Kienzle <pkienzle@users.sf.net> 79 ## 2005-03-02 Paul Kienzle <pkienzle@users.sf.net>
80 ## * Simplify 80 ## * Simplify
81 ## 2005-04-23 Dmitri A. Sergatskov <dasergatskov@gmail.com> 81 ## 2005-04-23 Dmitri A. Sergatskov <dasergatskov@gmail.com>
97 switch (nargin) 97 switch (nargin)
98 case 1 98 case 1
99 Z = varargin{1}; 99 Z = varargin{1};
100 case 2 100 case 2
101 if (ischar (varargin{2})) 101 if (ischar (varargin{2}))
102 [Z, method] = deal (varargin{:}); 102 [Z, method] = deal (varargin{:});
103 else 103 else
104 [Z, n] = deal (varargin{:}); 104 [Z, n] = deal (varargin{:});
105 endif 105 endif
106 case 3 106 case 3
107 if (ischar (varargin{3})) 107 if (ischar (varargin{3}))
108 [Z, n, method] = deal (varargin{:}); 108 [Z, n, method] = deal (varargin{:});
109 else 109 else
110 [Z, XI, YI] = deal (varargin{:}); 110 [Z, XI, YI] = deal (varargin{:});
111 endif 111 endif
112 case 4 112 case 4
113 if (ischar (varargin{4})) 113 if (ischar (varargin{4}))
114 [Z, XI, YI, method] = deal (varargin{:}); 114 [Z, XI, YI, method] = deal (varargin{:});
115 else 115 else
116 [Z, n, method, extrapval] = deal (varargin{:}); 116 [Z, n, method, extrapval] = deal (varargin{:});
117 endif 117 endif
118 case 5 118 case 5
119 if (ischar (varargin{4})) 119 if (ischar (varargin{4}))
120 [Z, XI, YI, method, extrapval] = deal (varargin{:}); 120 [Z, XI, YI, method, extrapval] = deal (varargin{:});
121 else 121 else
122 [X, Y, Z, XI, YI] = deal (varargin{:}); 122 [X, Y, Z, XI, YI] = deal (varargin{:});
123 endif 123 endif
124 case 6 124 case 6
125 [X, Y, Z, XI, YI, method] = deal (varargin{:}); 125 [X, Y, Z, XI, YI, method] = deal (varargin{:});
126 case 7 126 case 7
127 [X, Y, Z, XI, YI, method, extrapval] = deal (varargin{:}); 127 [X, Y, Z, XI, YI, method, extrapval] = deal (varargin{:});
128 otherwise 128 otherwise
129 print_usage (); 129 print_usage ();
130 endswitch 130 endswitch
131 131
132 ## Type checking. 132 ## Type checking.
179 error ("X and Y size must match Z dimensions"); 179 error ("X and Y size must match Z dimensions");
180 endif 180 endif
181 181
182 ## If Xi and Yi are vectors of different orientation build a grid 182 ## If Xi and Yi are vectors of different orientation build a grid
183 if ((rows (XI) == 1 && columns (YI) == 1) 183 if ((rows (XI) == 1 && columns (YI) == 1)
184 || (columns (XI) == 1 && rows (YI) == 1)) 184 || (columns (XI) == 1 && rows (YI) == 1))
185 [XI, YI] = meshgrid (XI, YI); 185 [XI, YI] = meshgrid (XI, YI);
186 elseif (! size_equal (XI, YI)) 186 elseif (! size_equal (XI, YI))
187 error ("XI and YI must be matrices of same size"); 187 error ("XI and YI must be matrices of same size");
188 endif 188 endif
189 189
237 ZI = Z(idx); 237 ZI = Z(idx);
238 238
239 elseif (strcmp (method, "pchip")) 239 elseif (strcmp (method, "pchip"))
240 240
241 if (length (X) < 2 || length (Y) < 2) 241 if (length (X) < 2 || length (Y) < 2)
242 error ("interp2: pchip2 requires at least 2 points in each dimension") 242 error ("interp2: pchip2 requires at least 2 points in each dimension")
243 endif 243 endif
244 244
245 ## first order derivatives 245 ## first order derivatives
246 DX = __pchip_deriv__ (X, Z, 2); 246 DX = __pchip_deriv__ (X, Z, 2);
247 DY = __pchip_deriv__ (Y, Z, 1); 247 DY = __pchip_deriv__ (Y, Z, 1);
282 yb{1,2} = -t2; 282 yb{1,2} = -t2;
283 yb{1,1} = t2 + 1; 283 yb{1,1} = t2 + 1;
284 284
285 ZI = zeros (size (XI)); 285 ZI = zeros (size (XI));
286 for i = 1:2 286 for i = 1:2
287 for j = 1:2 287 for j = 1:2
288 zidx = sub2ind (size (Z), yidx+(j-1), xidx+(i-1)); 288 zidx = sub2ind (size (Z), yidx+(j-1), xidx+(i-1));
289 ZI += xb{1,i} .* yb{1,j} .* Z(zidx); 289 ZI += xb{1,i} .* yb{1,j} .* Z(zidx);
290 ZI += xb{2,i} .* yb{1,j} .* DX(zidx); 290 ZI += xb{2,i} .* yb{1,j} .* DX(zidx);
291 ZI += xb{1,i} .* yb{2,j} .* DY(zidx); 291 ZI += xb{1,i} .* yb{2,j} .* DY(zidx);
292 ZI += xb{2,i} .* yb{2,j} .* DXY(zidx); 292 ZI += xb{2,i} .* yb{2,j} .* DXY(zidx);
293 endfor 293 endfor
294 endfor 294 endfor
295 295
296 endif 296 endif
297 297
298 if (! isempty (extrapval)) 298 if (! isempty (extrapval))
299 ## set points outside the table to 'extrapval' 299 ## set points outside the table to 'extrapval'
300 if (X (1) < X (end)) 300 if (X (1) < X (end))
301 if (Y (1) < Y (end)) 301 if (Y (1) < Y (end))
302 ZI (XI < X(1,1) | XI > X(end) | YI < Y(1,1) | YI > Y(end)) = ... 302 ZI (XI < X(1,1) | XI > X(end) | YI < Y(1,1) | YI > Y(end)) = ...
303 extrapval; 303 extrapval;
304 else 304 else
305 ZI (XI < X(1) | XI > X(end) | YI < Y(end) | YI > Y(1)) = ... 305 ZI (XI < X(1) | XI > X(end) | YI < Y(end) | YI > Y(1)) = ...
306 extrapval; 306 extrapval;
307 endif 307 endif
308 else 308 else
309 if (Y (1) < Y (end)) 309 if (Y (1) < Y (end))
310 ZI (XI < X(end) | XI > X(1) | YI < Y(1) | YI > Y(end)) = ... 310 ZI (XI < X(end) | XI > X(1) | YI < Y(1) | YI > Y(end)) = ...
311 extrapval; 311 extrapval;
312 else 312 else
313 ZI (XI < X(1,end) | XI > X(1) | YI < Y(end) | YI > Y(1)) = ... 313 ZI (XI < X(1,end) | XI > X(1) | YI < Y(end) | YI > Y(1)) = ...
314 extrapval; 314 extrapval;
315 endif 315 endif
316 endif 316 endif
317 endif 317 endif
318 318
319 else 319 else
320 320
321 ## If X and Y vectors produce a grid from them 321 ## If X and Y vectors produce a grid from them
322 if (isvector (X) && isvector (Y)) 322 if (isvector (X) && isvector (Y))
323 X = X(:).'; 323 X = X(:).';
324 Y = Y(:); 324 Y = Y(:);
325 if (!isequal ([length(X), length(Y)], size(Z))) 325 if (!isequal ([length(X), length(Y)], size(Z)))
326 error ("X and Y size must match Z dimensions"); 326 error ("X and Y size must match Z dimensions");
327 endif 327 endif
328 elseif (!size_equal (X, Y)) 328 elseif (!size_equal (X, Y))
329 error ("X and Y must be matrices of same size"); 329 error ("X and Y must be matrices of same size");
330 if (! size_equal (X, Z)) 330 if (! size_equal (X, Z))
331 error ("X and Y size must match Z dimensions"); 331 error ("X and Y size must match Z dimensions");
332 endif 332 endif
333 endif 333 endif
334 334
335 ## If Xi and Yi are vectors of different orientation build a grid 335 ## If Xi and Yi are vectors of different orientation build a grid
336 if ((rows (XI) == 1 && columns (YI) == 1) 336 if ((rows (XI) == 1 && columns (YI) == 1)
337 || (columns (XI) == 1 && rows (YI) == 1)) 337 || (columns (XI) == 1 && rows (YI) == 1))
338 ## Do nothing 338 ## Do nothing
339 elseif (! size_equal (XI, YI)) 339 elseif (! size_equal (XI, YI))
340 error ("XI and YI must be matrices of same size"); 340 error ("XI and YI must be matrices of same size");
341 endif 341 endif
342 342
394 endif 394 endif
395 395
396 elseif (strcmp (method, "spline")) 396 elseif (strcmp (method, "spline"))
397 if (isgriddata (XI) && isgriddata (YI')) 397 if (isgriddata (XI) && isgriddata (YI'))
398 ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrapval, 398 ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrapval,
399 "spline"); 399 "spline");
400 else 400 else
401 error ("interp2: input data must have `meshgrid' format"); 401 error ("interp2: input data must have `meshgrid' format");
402 endif 402 endif
403 else 403 else
404 error ("interpolation method not recognized"); 404 error ("interpolation method not recognized");