Mercurial > octave-nkf

comparison scripts/general/interp2.m @ 5837:55404f3b0da1

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