Mercurial > forge

comparison main/general/interp2.m @ 0:6b33357c7561 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial revision
author pkienzle
date Wed, 10 Oct 2001 19:54:49 +0000
parents
children 3165e038dba4
comparison
equal deleted inserted replaced
-1:000000000000 0:6b33357c7561
1 ## Copyright (C) 2000 Kai Habel
2 ##
3 ## This program is free software; you can redistribute it and/or modify
4 ## it under the terms of the GNU General Public License as published by
5 ## the Free Software Foundation; either version 2 of the License, or
6 ## (at your option) any later version.
7 ##
8 ## This program is distributed in the hope that it will be useful,
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 ## GNU General Public License for more details.
12 ##
13 ## You should have received a copy of the GNU General Public License
14 ## along with this program; if not, write to the Free Software
15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16
17 ## -*- texinfo -*-
18 ## @deftypefn {Function File} {@var{zi}=} interp2 (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi})
19 ## @deftypefnx {Function File} {@var{zi}=} interp2 (@var{Z}, @var{xi}, @var{yi})
20 ## @deftypefnx {Function File} {@var{zi}=} interp2 (@var{Z}, @var{n})
21 ## @deftypefnx {Function File} {@var{zi}=} interp2 (... , '@var{method}')
22 ## Two-dimensional interpolation. @var{X},@var{Y} and @var{Z} describe a
23 ## surface function. If @var{X} and @var{Y} are vectors their length
24 ## must correspondent to the size of @var{Z}. If they are matrices they
25 ## must have the 'meshgrid' format.
26 ##
27 ## ZI = interp2 (X, Y, Z, XI, YI, ...) returns a matrix corresponding
28 ## to the points described by the matrices @var{XI}, @var{YI}.
29 ##
30 ## If the last argument is a string, the interpolation method can
31 ## be specified. At the moment only 'linear' and 'nearest' methods
32 ## are provided. If it is omitted 'linear' interpolation is
33 ## assumed.
34 ##
35 ## ZI = interp2 (Z, XI, YI) assumes X = 1:rows(Z) and Y = 1:columns(Z)
36 ##
37 ## ZI = interp2 (Z, n) interleaves the Matrix Z n-times. If n is ommited
38 ## n = 1 is assumed
39 ##
40 ## @seealso{interp1}
41 ## @end deftypefn
42
43 ## Author: Kai Habel <kai.habel@gmx.de>
44
45 function ZI = interp2 (X, Y, Z, XI, YI, method)
46
47 if (nargin > 6 || nargin == 0)
48 usage ("interp2 (X, Y, Z, XI, YI, method)");
49 endif
50
51 if (nargin > 4)
52 if (is_vector (X) && is_vector (Y))
53 [rz, cz] = size (Z);
54 if (rz != length (Y) || cz != length (X))
55 error ("length of X and Y must match the size of Z");
56 endif
57 [X, Y] = meshgrid (X, Y);
58 elseif !( (size (X) == size (Y)) && (size (X) == size (Z)) )
59 error ("X,Y and Z must be matrices of same size");
60 endif
61 endif
62
63 if (((nargin == 4) || (nargin == 3)) && !isstr (Z))
64
65 if (nargin == 4)
66 if (isstr (XI))
67 method = XI;
68 else
69 usage("interp2 (z,xi,yi,'format'");
70 endif
71 endif
72 XI = Y;
73 YI = Z;
74 Z = X;
75 [X, Y] = meshgrid(1:columns(Z), 1:rows(Z));
76 else
77 if (nargin == 1)
78 n = 1;
79 elseif (nargin == 2)
80 if (isstr (Y))
81 method = Y;
82 n = 1;
83 elseif (is_scalar (Y))
84 n = Y;
85 endif
86 else
87 n = Y;
88 if (isstr (Z))
89 method = Z;
90 endif
91 endif
92 Z = X;
93 [zr, zc] = size (Z);
94 [X, Y] = meshgrid (1:zc, 1:zr);
95 xi = linspace (1, zc, pow2 (n) * (zc - 1) + 1);
96 yi = linspace (1, zr, pow2 (n) * (zr - 1) + 1);
97 [XI, YI] = meshgrid (xi, yi);
98 endif
99
100 if (! exist ("method"))
101 method = "linear";
102 endif
103
104 xtable = X(1, :);
105 ytable = Y(:, 1);
106
107 if (is_vector (XI) && is_vector (YI))
108 [XI, YI] = meshgrid (XI, YI);
109 elseif (! (size (XI) == size (YI)))
110 error ("XI and YI must be matrices of same size");
111 endif
112
113 ytlen = length (ytable);
114 xtlen = length (xtable);
115
116 ## get x index of values in XI
117 xtable(xtlen) *= (1 + eps);
118 xtable(xtlen) > XI(1, :);
119 [m, n] = sort ([xtable(:); XI(1, :)']);
120 o = cumsum (n <= xtlen);
121 xidx = o([find (n > xtlen)])';
122
123 ## get y index of values in YI
124 ytable(ytlen) *= (1 + eps);
125 [m, n]=sort ([ytable(:); YI(:, 1)]);
126 o = cumsum (n <= ytlen);
127 yidx = o([find (n > ytlen)]);
128
129 [zr, zc] = size (Z);
130
131 ## mark values outside the lookup table
132 xfirst_val = find (XI(1,:) < xtable(1));
133 xlast_val = find (XI(1,:) > xtable(xtlen));
134 yfirst_val = find (YI(:,1) < ytable(1));
135 ylast_val = find (YI(:,1) > ytable(ytlen));
136
137 ## set value outside the table preliminary to min max index
138 yidx(yfirst_val) = 1;
139 xidx(xfirst_val) = 1;
140 yidx(ylast_val) = zr - 1;
141 xidx(xlast_val) = zc - 1;
142
143 if strcmp (method, "linear")
144 ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy
145 ##
146 ## a-b
147 ## | |
148 ## c-d
149 a = Z(1:zr - 1, 1:zc - 1);
150 b = Z(1:zr - 1, 2:zc) - a;
151 c = Z(2:zr, 1:zc - 1) - a;
152 d = Z(2:zr, 2:zc) - a - b - c;
153
154 ## scale XI,YI values to a 1-spaced grid
155 Xsc = (XI .- X(yidx, xidx)) ./ (X(yidx, xidx + 1) - X(yidx, xidx));
156 Ysc = (YI .- Y(yidx, xidx)) ./ (Y(yidx + 1, xidx) - Y(yidx, xidx));
157 ## apply plane equation
158 ZI = a(yidx, xidx) .+ b(yidx, xidx) .* Xsc \
159 .+ c(yidx, xidx) .* Ysc .+ d(yidx, xidx) .* Xsc .* Ysc;
160 elseif strcmp (method, "nearest")
161 i = XI(1, :) - xtable(xidx) > xtable(xidx + 1) - XI(1, :);
162 j = YI(:, 1) - ytable(yidx) > ytable(yidx + 1) - YI(:, 1);
163 ZI = Z(yidx + j, xidx + i);
164 else
165 error ("interpolation method not (yet) supported");
166 endif
167
168 ## set points outside the table to NaN
169 if (! (isempty (xfirst_val) && isempty (xlast_val)))
170 ZI(:, [xfirst_val, xlast_val]) = NaN;
171 endif
172 if (! (isempty (yfirst_val) && isempty (ylast_val)))
173 ZI([yfirst_val; ylast_val], :) = NaN;
174 endif
175 endfunction
176
177 %!demo
178 %! A=[13,-1,12;5,4,3;1,6,2];
179 %! x=[0,1,4]; y=[10,11,12];
180 %! xi=linspace(min(x),max(x),17);
181 %! yi=linspace(min(y),max(y),26);
182 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'linear'));
183 %! [x,y] = meshgrid(x,y); gset nohidden3d;
184 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
185
186 %!demo
187 %! A=[13,-1,12;5,4,3;1,6,2];
188 %! x=[0,1,4]; y=[10,11,12];
189 %! xi=linspace(min(x),max(x),17);
190 %! yi=linspace(min(y),max(y),26);
191 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'nearest'));
192 %! [x,y] = meshgrid(x,y); gset nohidden3d;
193 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
194
195 %!#demo
196 %! A=[13,-1,12;5,4,3;1,6,2];
197 %! x=[0,1,2]; y=[10,11,12];
198 %! xi=linspace(min(x),max(x),17);
199 %! yi=linspace(min(y),max(y),26);
200 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic'));
201 %! [x,y] = meshgrid(x,y); gset nohidden3d;
202 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;