Mercurial > hg > octave-lyh
annotate scripts/general/interp2.m @ 10780:6e7590d003dc
Fix bugs in interp2 (bug #29601)
author | Thorsten Meyer <thorsten.meyier@gmx.de> |
---|---|
date | Fri, 09 Jul 2010 19:56:21 +0200 |
parents | d1978e7364ad |
children | 3140cb7a05a1 |
rev | line source |
---|---|
8920 | 1 ## Copyright (C) 2000, 2006, 2007, 2008 Kai Habel |
2 ## Copyright (C) 2009 Jaroslav Hajek | |
5837 | 3 ## |
4 ## This file is part of Octave. | |
5 ## | |
6 ## Octave is free software; you can redistribute it and/or modify it | |
7 ## under the terms of the GNU General Public License as published by | |
7016 | 8 ## the Free Software Foundation; either version 3 of the License, or (at |
9 ## your option) any later version. | |
5837 | 10 ## |
11 ## Octave is distributed in the hope that it will be useful, but | |
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
14 ## General Public License for more details. | |
15 ## | |
16 ## You should have received a copy of the GNU General Public License | |
7016 | 17 ## along with Octave; see the file COPYING. If not, see |
18 ## <http://www.gnu.org/licenses/>. | |
5837 | 19 |
20 ## -*- texinfo -*- | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
21 ## @deftypefn {Function File} {@var{zi} =} interp2 (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}) |
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
22 ## @deftypefnx {Function File} {@var{zi} =} interp2 (@var{Z}, @var{xi}, @var{yi}) |
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
23 ## @deftypefnx {Function File} {@var{zi} =} interp2 (@var{Z}, @var{n}) |
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
24 ## @deftypefnx {Function File} {@var{zi} =} interp2 (@dots{}, @var{method}) |
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
25 ## @deftypefnx {Function File} {@var{zi} =} interp2 (@dots{}, @var{method}, @var{extrapval}) |
5837 | 26 ## |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
27 ## Two-dimensional interpolation. @var{x}, @var{y} and @var{z} describe a |
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
28 ## surface function. If @var{x} and @var{y} are vectors their length |
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
29 ## must correspondent to the size of @var{z}. @var{x} and @var{y} must be |
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
30 ## monotonic. If they are matrices they must have the @code{meshgrid} |
5837 | 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 | |
8478
5aafd61e4c3c
[docs] @var{XI}, @var{YI}. => @var{xi}, @var{yi}.
Brian Gough <bjg@gnu.org>
parents:
8347
diff
changeset
|
36 ## matrices @var{xi}, @var{yi}. |
5837 | 37 ## |
38 ## If the last argument is a string, the interpolation method can | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
39 ## be specified. The method can be 'linear', 'nearest' or 'cubic'. |
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
40 ## If it is omitted 'linear' interpolation is assumed. |
5837 | 41 ## |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
42 ## @item interp2 (@var{z}, @var{xi}, @var{yi}) |
5837 | 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}) | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
47 ## Interleaves the matrix @var{z} n-times. If @var{n} is omitted a value |
5837 | 48 ## of @code{@var{n} = 1} is assumed. |
49 ## @end table | |
50 ## | |
51 ## The variable @var{method} defines the method to use for the | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
52 ## interpolation. It can take one of the following values |
5837 | 53 ## |
54 ## @table @asis | |
55 ## @item 'nearest' | |
6218 | 56 ## Return the nearest neighbor. |
5837 | 57 ## @item 'linear' |
6218 | 58 ## Linear interpolation from nearest neighbors. |
5837 | 59 ## @item 'pchip' |
9755
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
60 ## Piece-wise cubic hermite interpolating polynomial. |
5837 | 61 ## @item 'cubic' |
6218 | 62 ## Cubic interpolation from four nearest neighbors. |
5837 | 63 ## @item 'spline' |
64 ## Cubic spline interpolation--smooth first and second derivatives | |
6702 | 65 ## throughout the curve. |
5837 | 66 ## @end table |
67 ## | |
68 ## If a scalar value @var{extrapval} is defined as the final value, then | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
69 ## values outside the mesh as set to this value. Note that in this case |
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
70 ## @var{method} must be defined as well. If @var{extrapval} is not |
6742 | 71 ## defined then NA is assumed. |
5837 | 72 ## |
73 ## @seealso{interp1} | |
74 ## @end deftypefn | |
75 | |
10549 | 76 ## Author: Kai Habel <kai.habel@gmx.de> |
5837 | 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 | |
5838 | 85 ## * Fix the eat line bug when the last element of XI or YI is |
86 ## negative or zero. | |
5837 | 87 ## 2005-11-26 Pierre Baldensperger <balden@libertysurf.fr> |
88 ## * Rather big modification (XI,YI no longer need to be | |
89 ## "meshgridded") to be consistent with the help message | |
90 ## above and for compatibility. | |
91 | |
92 function ZI = interp2 (varargin) | |
5838 | 93 Z = X = Y = XI = YI = n = []; |
5837 | 94 method = "linear"; |
6742 | 95 extrapval = NA; |
5837 | 96 |
5838 | 97 switch (nargin) |
98 case 1 | |
99 Z = varargin{1}; | |
100 case 2 | |
101 if (ischar (varargin{2})) | |
10549 | 102 [Z, method] = deal (varargin{:}); |
5838 | 103 else |
10549 | 104 [Z, n] = deal (varargin{:}); |
5838 | 105 endif |
106 case 3 | |
107 if (ischar (varargin{3})) | |
10549 | 108 [Z, n, method] = deal (varargin{:}); |
5838 | 109 else |
10549 | 110 [Z, XI, YI] = deal (varargin{:}); |
5838 | 111 endif |
112 case 4 | |
113 if (ischar (varargin{4})) | |
10549 | 114 [Z, XI, YI, method] = deal (varargin{:}); |
5838 | 115 else |
10549 | 116 [Z, n, method, extrapval] = deal (varargin{:}); |
5838 | 117 endif |
118 case 5 | |
119 if (ischar (varargin{4})) | |
10549 | 120 [Z, XI, YI, method, extrapval] = deal (varargin{:}); |
5838 | 121 else |
10549 | 122 [X, Y, Z, XI, YI] = deal (varargin{:}); |
5838 | 123 endif |
124 case 6 | |
10549 | 125 [X, Y, Z, XI, YI, method] = deal (varargin{:}); |
5838 | 126 case 7 |
10549 | 127 [X, Y, Z, XI, YI, method, extrapval] = deal (varargin{:}); |
5838 | 128 otherwise |
129 print_usage (); | |
5837 | 130 endswitch |
131 | |
132 ## Type checking. | |
5838 | 133 if (!ismatrix (Z)) |
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
134 error ("interp2: expected matrix Z"); |
5837 | 135 endif |
5838 | 136 if (!isempty (n) && !isscalar (n)) |
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
137 error ("interp2: expected scalar n"); |
5837 | 138 endif |
5838 | 139 if (!ischar (method)) |
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
140 error ("interp2: expected string 'method'"); |
5837 | 141 endif |
8712
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
142 if (ischar (extrapval) || strcmp (extrapval, "extrap")) |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
143 extrapval = []; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
144 elseif (!isscalar (extrapval)) |
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
145 error ("interp2: expected n extrapval"); |
5837 | 146 endif |
147 | |
5838 | 148 ## Define X, Y, XI, YI if needed |
5837 | 149 [zr, zc] = size (Z); |
5838 | 150 if (isempty (X)) |
151 X = 1:zc; | |
152 Y = 1:zr; | |
5837 | 153 endif |
5838 | 154 if (! isnumeric (X) || ! isnumeric (Y)) |
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
155 error ("interp2: expected numeric X, Y"); |
5837 | 156 endif |
5838 | 157 if (! isempty (n)) |
158 p = 2^n; | |
159 XI = (p:p*zc)/p; | |
160 YI = (p:p*zr)'/p; | |
5837 | 161 endif |
5838 | 162 if (! isnumeric (XI) || ! isnumeric (YI)) |
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
163 error ("interp2: expected numeric XI, YI"); |
5837 | 164 endif |
165 | |
6702 | 166 |
8712
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
167 if (strcmp (method, "linear") || strcmp (method, "nearest") ... |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
168 || strcmp (method, "pchip")) |
6702 | 169 |
170 ## If X and Y vectors produce a grid from them | |
171 if (isvector (X) && isvector (Y)) | |
8712
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
172 X = X(:); Y = Y(:); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
173 elseif (size_equal (X, Y)) |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
174 X = X(1,:)'; Y = Y(:,1); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
175 else |
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
176 error ("interp2: X and Y must be matrices of same size"); |
6702 | 177 endif |
8712
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
178 if (columns (Z) != length (X) || rows (Z) != length (Y)) |
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
179 error ("interp2: X and Y size must match Z dimensions"); |
6702 | 180 endif |
181 | |
182 ## If Xi and Yi are vectors of different orientation build a grid | |
183 if ((rows (XI) == 1 && columns (YI) == 1) | |
10549 | 184 || (columns (XI) == 1 && rows (YI) == 1)) |
6702 | 185 [XI, YI] = meshgrid (XI, YI); |
186 elseif (! size_equal (XI, YI)) | |
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
187 error ("interp2: XI and YI must be matrices of same size"); |
6702 | 188 endif |
5837 | 189 |
8712
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
190 ## if XI, YI are vectors, X and Y should share their orientation. |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
191 if (rows (XI) == 1) |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
192 if (rows (X) != 1) |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
193 X = X.'; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
194 endif |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
195 if (rows (Y) != 1) |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
196 Y = Y.'; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
197 endif |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
198 elseif (columns (XI) == 1) |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
199 if (columns (X) != 1) |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
200 X = X.'; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
201 endif |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
202 if (columns (Y) != 1) |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
203 Y = Y.'; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
204 endif |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
205 endif |
6702 | 206 |
8712
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
207 xidx = lookup (X, XI, "lr"); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
208 yidx = lookup (Y, YI, "lr"); |
5837 | 209 |
6702 | 210 if (strcmp (method, "linear")) |
211 ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy | |
212 ## | |
213 ## a-b | |
214 ## | | | |
215 ## c-d | |
216 a = Z(1:(zr - 1), 1:(zc - 1)); | |
217 b = Z(1:(zr - 1), 2:zc) - a; | |
218 c = Z(2:zr, 1:(zc - 1)) - a; | |
219 d = Z(2:zr, 2:zc) - a - b - c; | |
5837 | 220 |
6702 | 221 ## scale XI, YI values to a 1-spaced grid |
10326
ade59ae10e80
slight optimization in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
10320
diff
changeset
|
222 Xsc = (XI - X(xidx)) ./ (diff (X)(xidx)); |
ade59ae10e80
slight optimization in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
10320
diff
changeset
|
223 Ysc = (YI - Y(yidx)) ./ (diff (Y)(yidx)); |
6702 | 224 |
10328
2210d3070543
further memory optimization in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
10326
diff
changeset
|
225 ## Get 2D index. |
2210d3070543
further memory optimization in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
10326
diff
changeset
|
226 idx = sub2ind (size (a), yidx, xidx); |
2210d3070543
further memory optimization in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
10326
diff
changeset
|
227 ## We can dispose of the 1D indices at this point to save memory. |
2210d3070543
further memory optimization in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
10326
diff
changeset
|
228 clear xidx yidx |
2210d3070543
further memory optimization in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
10326
diff
changeset
|
229 |
6702 | 230 ## apply plane equation |
231 ZI = a(idx) + b(idx).*Xsc + c(idx).*Ysc + d(idx).*Xsc.*Ysc; | |
5837 | 232 |
6702 | 233 elseif (strcmp (method, "nearest")) |
10427
62bb59f927b1
scripts/general/interp2.m, scripts/general/interpn.m: For nearest neighbour interpolation ceil (instead of floor) at the center of the data intervals to be compatible with Matlab. Add test.
Soren Hauberg <hauberg@gmail.com>
parents:
10328
diff
changeset
|
234 ii = (XI - X(xidx) >= X(xidx + 1) - XI); |
62bb59f927b1
scripts/general/interp2.m, scripts/general/interpn.m: For nearest neighbour interpolation ceil (instead of floor) at the center of the data intervals to be compatible with Matlab. Add test.
Soren Hauberg <hauberg@gmail.com>
parents:
10328
diff
changeset
|
235 jj = (YI - Y(yidx) >= Y(yidx + 1) - YI); |
6702 | 236 idx = sub2ind (size (Z), yidx+jj, xidx+ii); |
237 ZI = Z(idx); | |
8712
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
238 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
239 elseif (strcmp (method, "pchip")) |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
240 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
241 if (length (X) < 2 || length (Y) < 2) |
10549 | 242 error ("interp2: pchip2 requires at least 2 points in each dimension") |
8712
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
243 endif |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
244 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
245 ## first order derivatives |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
246 DX = __pchip_deriv__ (X, Z, 2); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
247 DY = __pchip_deriv__ (Y, Z, 1); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
248 ## Compute mixed derivatives row-wise and column-wise, use the average. |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
249 DXY = (__pchip_deriv__ (X, DY, 2) + __pchip_deriv__ (Y, DX, 1))/2; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
250 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
251 ## do the bicubic interpolation |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
252 hx = diff (X); hx = hx(xidx); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
253 hy = diff (Y); hy = hy(yidx); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
254 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
255 tx = (XI - X(xidx)) ./ hx; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
256 ty = (YI - Y(yidx)) ./ hy; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
257 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
258 ## construct the cubic hermite base functions in x, y |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
259 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
260 ## formulas: |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
261 ## b{1,1} = ( 2*t.^3 - 3*t.^2 + 1); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
262 ## b{2,1} = h.*( t.^3 - 2*t.^2 + t ); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
263 ## b{1,2} = (-2*t.^3 + 3*t.^2 ); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
264 ## b{2,2} = h.*( t.^3 - t.^2 ); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
265 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
266 ## optimized equivalents of the above: |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
267 t1 = tx.^2; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
268 t2 = tx.*t1 - t1; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
269 xb{2,2} = hx.*t2; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
270 t1 = t2 - t1; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
271 xb{2,1} = hx.*(t1 + tx); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
272 t2 += t1; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
273 xb{1,2} = -t2; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
274 xb{1,1} = t2 + 1; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
275 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
276 t1 = ty.^2; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
277 t2 = ty.*t1 - t1; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
278 yb{2,2} = hy.*t2; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
279 t1 = t2 - t1; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
280 yb{2,1} = hy.*(t1 + ty); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
281 t2 += t1; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
282 yb{1,2} = -t2; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
283 yb{1,1} = t2 + 1; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
284 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
285 ZI = zeros (size (XI)); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
286 for i = 1:2 |
10549 | 287 for j = 1:2 |
288 zidx = sub2ind (size (Z), yidx+(j-1), xidx+(i-1)); | |
289 ZI += xb{1,i} .* yb{1,j} .* Z(zidx); | |
290 ZI += xb{2,i} .* yb{1,j} .* DX(zidx); | |
291 ZI += xb{1,i} .* yb{2,j} .* DY(zidx); | |
292 ZI += xb{2,i} .* yb{2,j} .* DXY(zidx); | |
293 endfor | |
8712
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
294 endfor |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
295 |
6702 | 296 endif |
5837 | 297 |
8712
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
298 if (! isempty (extrapval)) |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
299 ## set points outside the table to 'extrapval' |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
300 if (X (1) < X (end)) |
10549 | 301 if (Y (1) < Y (end)) |
302 ZI (XI < X(1,1) | XI > X(end) | YI < Y(1,1) | YI > Y(end)) = ... | |
303 extrapval; | |
304 else | |
305 ZI (XI < X(1) | XI > X(end) | YI < Y(end) | YI > Y(1)) = ... | |
306 extrapval; | |
307 endif | |
6979 | 308 else |
10549 | 309 if (Y (1) < Y (end)) |
310 ZI (XI < X(end) | XI > X(1) | YI < Y(1) | YI > Y(end)) = ... | |
311 extrapval; | |
312 else | |
313 ZI (XI < X(1,end) | XI > X(1) | YI < Y(end) | YI > Y(1)) = ... | |
314 extrapval; | |
315 endif | |
6979 | 316 endif |
317 endif | |
318 | |
6702 | 319 else |
5837 | 320 |
10780
6e7590d003dc
Fix bugs in interp2 (bug #29601)
Thorsten Meyer <thorsten.meyier@gmx.de>
parents:
10635
diff
changeset
|
321 ## Check dimensions of X and Y |
6702 | 322 if (isvector (X) && isvector (Y)) |
323 X = X(:).'; | |
324 Y = Y(:); | |
325 if (!isequal ([length(X), length(Y)], size(Z))) | |
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
326 error ("interp2: X and Y size must match Z dimensions"); |
6702 | 327 endif |
328 elseif (!size_equal (X, Y)) | |
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
329 error ("interp2: X and Y must be matrices of same size"); |
6702 | 330 if (! size_equal (X, Z)) |
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
331 error ("interp2: X and Y size must match Z dimensions"); |
6702 | 332 endif |
333 endif | |
5837 | 334 |
10780
6e7590d003dc
Fix bugs in interp2 (bug #29601)
Thorsten Meyer <thorsten.meyier@gmx.de>
parents:
10635
diff
changeset
|
335 ## Check dimensions of XI and YI |
6e7590d003dc
Fix bugs in interp2 (bug #29601)
Thorsten Meyer <thorsten.meyier@gmx.de>
parents:
10635
diff
changeset
|
336 if (isvector (XI) && isvector (YI)) |
6e7590d003dc
Fix bugs in interp2 (bug #29601)
Thorsten Meyer <thorsten.meyier@gmx.de>
parents:
10635
diff
changeset
|
337 XI = XI(:).'; |
6e7590d003dc
Fix bugs in interp2 (bug #29601)
Thorsten Meyer <thorsten.meyier@gmx.de>
parents:
10635
diff
changeset
|
338 YI = YI(:); |
6702 | 339 elseif (! size_equal (XI, YI)) |
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
340 error ("interp2: XI and YI must be matrices of same size"); |
6702 | 341 endif |
5837 | 342 |
6702 | 343 ## FIXME bicubic/__splinen__ don't handle arbitrary XI, YI |
344 if (strcmp (method, "cubic")) | |
9755
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
345 if (isgriddata (XI) && isgriddata (YI')) |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
346 ZI = bicubic (X, Y, Z, XI (1, :), YI (:, 1), extrapval); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
347 elseif (isgriddata (X) && isgriddata (Y')) |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
348 ## Allocate output |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
349 ZI = zeros (size (X)); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
350 |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
351 ## Find inliers |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
352 inside = !(XI < X (1) | XI > X (end) | YI < Y (1) | YI > Y (end)); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
353 |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
354 ## Scale XI and YI to match indices of Z |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
355 XI = (columns (Z) - 1) * (XI - X (1)) / (X (end) - X (1)) + 1; |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
356 YI = (rows (Z) - 1) * (YI - Y (1)) / (Y (end) - Y (1)) + 1; |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
357 |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
358 ## Start the real work |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
359 K = floor (XI); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
360 L = floor (YI); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
361 |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
362 ## Coefficients |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
363 AY1 = bc ((YI - L + 1)); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
364 AX1 = bc ((XI - K + 1)); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
365 AY0 = bc ((YI - L + 0)); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
366 AX0 = bc ((XI - K + 0)); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
367 AY_1 = bc ((YI - L - 1)); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
368 AX_1 = bc ((XI - K - 1)); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
369 AY_2 = bc ((YI - L - 2)); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
370 AX_2 = bc ((XI - K - 2)); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
371 |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
372 ## Perform interpolation |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
373 sz = size(Z); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
374 ZI = AY_2 .* AX_2 .* Z (sym_sub2ind (sz, L+2, K+2)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
375 + AY_2 .* AX_1 .* Z (sym_sub2ind (sz, L+2, K+1)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
376 + AY_2 .* AX0 .* Z (sym_sub2ind (sz, L+2, K)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
377 + AY_2 .* AX1 .* Z (sym_sub2ind (sz, L+2, K-1)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
378 + AY_1 .* AX_2 .* Z (sym_sub2ind (sz, L+1, K+2)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
379 + AY_1 .* AX_1 .* Z (sym_sub2ind (sz, L+1, K+1)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
380 + AY_1 .* AX0 .* Z (sym_sub2ind (sz, L+1, K)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
381 + AY_1 .* AX1 .* Z (sym_sub2ind (sz, L+1, K-1)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
382 + AY0 .* AX_2 .* Z (sym_sub2ind (sz, L, K+2)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
383 + AY0 .* AX_1 .* Z (sym_sub2ind (sz, L, K+1)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
384 + AY0 .* AX0 .* Z (sym_sub2ind (sz, L, K)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
385 + AY0 .* AX1 .* Z (sym_sub2ind (sz, L, K-1)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
386 + AY1 .* AX_2 .* Z (sym_sub2ind (sz, L-1, K+2)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
387 + AY1 .* AX_1 .* Z (sym_sub2ind (sz, L-1, K+1)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
388 + AY1 .* AX0 .* Z (sym_sub2ind (sz, L-1, K)) ... |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
389 + AY1 .* AX1 .* Z (sym_sub2ind (sz, L-1, K-1)); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
390 ZI (!inside) = extrapval; |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
391 |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
392 else |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
393 error ("interp2: input data must have `meshgrid' format"); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
394 endif |
5837 | 395 |
6702 | 396 elseif (strcmp (method, "spline")) |
9755
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
397 if (isgriddata (XI) && isgriddata (YI')) |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
398 ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrapval, |
10549 | 399 "spline"); |
9755
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
400 else |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
401 error ("interp2: input data must have `meshgrid' format"); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
402 endif |
6702 | 403 else |
10635
d1978e7364ad
Print name of function in error() string messages.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
404 error ("interp2: interpolation method not recognized"); |
6702 | 405 endif |
5837 | 406 |
6702 | 407 endif |
5837 | 408 endfunction |
409 | |
9755
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
410 function b = isgriddata (X) |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
411 d1 = diff (X, 1, 1); |
10780
6e7590d003dc
Fix bugs in interp2 (bug #29601)
Thorsten Meyer <thorsten.meyier@gmx.de>
parents:
10635
diff
changeset
|
412 b = all (d1 (:) == 0); |
9755
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
413 endfunction |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
414 |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
415 ## Compute the bicubic interpolation coefficients |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
416 function o = bc(x) |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
417 x = abs(x); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
418 o = zeros(size(x)); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
419 idx1 = (x < 1); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
420 idx2 = !idx1 & (x < 2); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
421 o(idx1) = 1 - 2.*x(idx1).^2 + x(idx1).^3; |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
422 o(idx2) = 4 - 8.*x(idx2) + 5.*x(idx2).^2 - x(idx2).^3; |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
423 endfunction |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
424 |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
425 ## This version of sub2ind behaves as if the data was symmetrically padded |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
426 function ind = sym_sub2ind(sz, Y, X) |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
427 Y (Y < 1) = 1 - Y (Y < 1); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
428 while (any (Y (:) > 2 * sz (1))) |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
429 Y (Y > 2 * sz (1)) = round (Y (Y > 2 * sz (1)) / 2); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
430 endwhile |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
431 Y (Y > sz (1)) = 1 + 2 * sz (1) - Y (Y > sz (1)); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
432 X (X < 1) = 1 - X (X < 1); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
433 while (any (X (:) > 2 * sz (2))) |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
434 X (X > 2 * sz (2)) = round (X (X > 2 * sz (2)) / 2); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
435 endwhile |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
436 X (X > sz (2)) = 1 + 2 * sz (2) - X (X > sz (2)); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
437 ind = sub2ind(sz, Y, X); |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
438 endfunction |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
439 |
4f4873f6f873
general/interp2.m: improved error checking and support for bicubic interpolation when X and Y are meshgrid format
Soren Hauberg <hauberg@gmail.com>
parents:
9399
diff
changeset
|
440 |
5837 | 441 %!demo |
442 %! A=[13,-1,12;5,4,3;1,6,2]; | |
443 %! x=[0,1,4]; y=[10,11,12]; | |
444 %! xi=linspace(min(x),max(x),17); | |
445 %! yi=linspace(min(y),max(y),26)'; | |
446 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'linear')); | |
447 %! [x,y] = meshgrid(x,y); | |
448 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; | |
449 | |
450 %!demo | |
8712
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
451 %! [x,y,A] = peaks(10); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
452 %! x = x(1,:)'; y = y(:,1); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
453 %! xi=linspace(min(x),max(x),41); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
454 %! yi=linspace(min(y),max(y),41)'; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
455 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'linear')); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
456 %! [x,y] = meshgrid(x,y); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
457 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
458 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
459 %!demo |
5837 | 460 %! A=[13,-1,12;5,4,3;1,6,2]; |
461 %! x=[0,1,4]; y=[10,11,12]; | |
462 %! xi=linspace(min(x),max(x),17); | |
463 %! yi=linspace(min(y),max(y),26)'; | |
464 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'nearest')); | |
465 %! [x,y] = meshgrid(x,y); | |
466 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; | |
467 | |
6702 | 468 %!demo |
8712
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
469 %! [x,y,A] = peaks(10); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
470 %! x = x(1,:)'; y = y(:,1); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
471 %! xi=linspace(min(x),max(x),41); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
472 %! yi=linspace(min(y),max(y),41)'; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
473 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'nearest')); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
474 %! [x,y] = meshgrid(x,y); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
475 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
476 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
477 %!demo |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
478 %! A=[13,-1,12;5,4,3;1,6,2]; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
479 %! x=[0,1,2]; y=[10,11,12]; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
480 %! xi=linspace(min(x),max(x),17); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
481 %! yi=linspace(min(y),max(y),26)'; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
482 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'pchip')); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
483 %! [x,y] = meshgrid(x,y); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
484 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
485 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
486 %!demo |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
487 %! [x,y,A] = peaks(10); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
488 %! x = x(1,:)'; y = y(:,1); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
489 %! xi=linspace(min(x),max(x),41); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
490 %! yi=linspace(min(y),max(y),41)'; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
491 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'pchip')); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
492 %! [x,y] = meshgrid(x,y); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
493 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
494 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
495 %!demo |
5837 | 496 %! A=[13,-1,12;5,4,3;1,6,2]; |
497 %! x=[0,1,2]; y=[10,11,12]; | |
498 %! xi=linspace(min(x),max(x),17); | |
6702 | 499 %! yi=linspace(min(y),max(y),26)'; |
5837 | 500 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic')); |
501 %! [x,y] = meshgrid(x,y); | |
502 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; | |
503 | |
6702 | 504 %!demo |
8712
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
505 %! [x,y,A] = peaks(10); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
506 %! x = x(1,:)'; y = y(:,1); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
507 %! xi=linspace(min(x),max(x),41); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
508 %! yi=linspace(min(y),max(y),41)'; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
509 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic')); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
510 %! [x,y] = meshgrid(x,y); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
511 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
512 |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
513 %!demo |
6702 | 514 %! A=[13,-1,12;5,4,3;1,6,2]; |
515 %! x=[0,1,2]; y=[10,11,12]; | |
516 %! xi=linspace(min(x),max(x),17); | |
517 %! yi=linspace(min(y),max(y),26)'; | |
518 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'spline')); | |
519 %! [x,y] = meshgrid(x,y); | |
520 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; | |
521 | |
8712
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
522 %!demo |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
523 %! [x,y,A] = peaks(10); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
524 %! x = x(1,:)'; y = y(:,1); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
525 %! xi=linspace(min(x),max(x),41); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
526 %! yi=linspace(min(y),max(y),41)'; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
527 %! mesh(xi,yi,interp2(x,y,A,xi,yi,'spline')); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
528 %! [x,y] = meshgrid(x,y); |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
529 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; |
010e15c7be9a
support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents:
8479
diff
changeset
|
530 |
5837 | 531 %!test % simple test |
532 %! x = [1,2,3]; | |
533 %! y = [4,5,6,7]; | |
534 %! [X, Y] = meshgrid(x,y); | |
535 %! Orig = X.^2 + Y.^3; | |
536 %! xi = [1.2,2, 1.5]; | |
537 %! yi = [6.2, 4.0, 5.0]'; | |
538 %! | |
539 %! Expected = ... | |
540 %! [243, 245.4, 243.9; | |
541 %! 65.6, 68, 66.5; | |
542 %! 126.6, 129, 127.5]; | |
543 %! Result = interp2(x,y,Orig, xi, yi); | |
544 %! | |
545 %! assert(Result, Expected, 1000*eps); | |
546 | |
547 %!test % 2^n form | |
548 %! x = [1,2,3]; | |
549 %! y = [4,5,6,7]; | |
550 %! [X, Y] = meshgrid(x,y); | |
551 %! Orig = X.^2 + Y.^3; | |
552 %! xi = [1:0.25:3]; yi = [4:0.25:7]'; | |
553 %! Expected = interp2(x,y,Orig, xi, yi); | |
554 %! Result = interp2(Orig,2); | |
555 %! | |
556 %! assert(Result, Expected, 10*eps); | |
557 | |
558 %!test % matrix slice | |
559 %! A = eye(4); | |
560 %! assert(interp2(A,[1:4],[1:4]),[1,1,1,1]); | |
561 | |
562 %!test % non-gridded XI,YI | |
563 %! A = eye(4); | |
564 %! assert(interp2(A,[1,2;3,4],[1,3;2,4]),[1,0;0,1]); | |
565 | |
566 %!test % for values outside of boundaries | |
567 %! x = [1,2,3]; | |
568 %! y = [4,5,6,7]; | |
569 %! [X, Y] = meshgrid(x,y); | |
570 %! Orig = X.^2 + Y.^3; | |
571 %! xi = [0,4]; | |
572 %! yi = [3,8]'; | |
6742 | 573 %! assert(interp2(x,y,Orig, xi, yi),[NA,NA;NA,NA]); |
5837 | 574 %! assert(interp2(x,y,Orig, xi, yi,'linear', 0),[0,0;0,0]); |
575 | |
576 %!test % for values at boundaries | |
577 %! A=[1,2;3,4]; | |
578 %! x=[0,1]; | |
579 %! y=[2,3]'; | |
580 %! assert(interp2(x,y,A,x,y,'linear'), A); | |
581 %! assert(interp2(x,y,A,x,y,'nearest'), A); | |
582 | |
10427
62bb59f927b1
scripts/general/interp2.m, scripts/general/interpn.m: For nearest neighbour interpolation ceil (instead of floor) at the center of the data intervals to be compatible with Matlab. Add test.
Soren Hauberg <hauberg@gmail.com>
parents:
10328
diff
changeset
|
583 %!test % for Matlab-compatible rounding for 'nearest' |
62bb59f927b1
scripts/general/interp2.m, scripts/general/interpn.m: For nearest neighbour interpolation ceil (instead of floor) at the center of the data intervals to be compatible with Matlab. Add test.
Soren Hauberg <hauberg@gmail.com>
parents:
10328
diff
changeset
|
584 %! X = meshgrid (1:4); |
62bb59f927b1
scripts/general/interp2.m, scripts/general/interpn.m: For nearest neighbour interpolation ceil (instead of floor) at the center of the data intervals to be compatible with Matlab. Add test.
Soren Hauberg <hauberg@gmail.com>
parents:
10328
diff
changeset
|
585 %! assert (interp2 (X, 2.5, 2.5, 'nearest'), 3); |
62bb59f927b1
scripts/general/interp2.m, scripts/general/interpn.m: For nearest neighbour interpolation ceil (instead of floor) at the center of the data intervals to be compatible with Matlab. Add test.
Soren Hauberg <hauberg@gmail.com>
parents:
10328
diff
changeset
|
586 |