Mercurial > hg > octave-lyh
annotate scripts/general/interp2.m @ 8347:fa78cb8d8a5c
corrections for typos
Here is a patch with some corrections for typos and missing/extra
words in the manual.
changeset: 8347:34fd1d1c2294
user: Brian Gough <bjg@gnu.org>
date: Wed Nov 26 11:00:15 2008 -0500
summary: [docs] can not => cannot
| author | Brian Gough<bjg@network-theory.co.uk> |
|---|---|
| date | Thu, 27 Nov 2008 10:28:24 +0100 |
| parents | 4fbaba9abec1 |
| children | 5aafd61e4c3c |
| rev | line source |
|---|---|
| 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 | |
|
8347
fa78cb8d8a5c
corrections for typos
Brian Gough<bjg@network-theory.co.uk>
parents:
7671
diff
changeset
|
28 ## must correspondent to the size of @var{z}. @var{x} and @var{y} 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 | |
|
7671
4fbaba9abec1
implement compiled binary lookup
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
188 xidx = lookup (X(1, :), XI, "lr"); |
|
4fbaba9abec1
implement compiled binary lookup
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
189 yidx = lookup (Y(:, 1), YI, "lr"); |
| 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 |