|
6823
|
1 ## Copyright (C) 2007 David Bateman <dbateman@free.fr> |
|
|
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{yi} =} griddatan (@var{x}, @var{y}, @var{xi}, @var{method}, @var{options}) |
|
|
22 ## |
|
|
23 ## Generate a regular mesh from irregular data using interpolation. |
|
|
24 ## The function is defined by @code{@var{y} = f (@var{x})}. |
|
|
25 ## The interpolation points are all @var{xi}. |
|
|
26 ## |
|
6826
|
27 ## The interpolation method can be @code{"nearest"} or @code{"linear"}. |
|
|
28 ## If method is omitted it defaults to @code{"linear"}. |
|
6823
|
29 ## @seealso{griddata, delaunayn} |
|
|
30 ## @end deftypefn |
|
|
31 |
|
6826
|
32 function yi = griddatan (x, y, xi, method, varargin) |
|
|
33 |
|
6823
|
34 if (nargin == 3) |
|
6826
|
35 method = "linear"; |
|
6823
|
36 endif |
|
|
37 if (nargin < 3) |
|
6826
|
38 print_usage (); |
|
6823
|
39 endif |
|
|
40 |
|
6826
|
41 if (ischar (method)) |
|
|
42 method = tolower (method); |
|
6823
|
43 endif |
|
|
44 |
|
6826
|
45 [m, n] = size (x); |
|
|
46 [mi, ni] = size (xi); |
|
6823
|
47 |
|
6826
|
48 if (n != ni || size (y, 1) != m || size (y, 2) != 1) |
|
6823
|
49 error ("griddatan: dimensional mismatch"); |
|
|
50 endif |
|
|
51 |
|
|
52 ## triangulate data |
|
|
53 ## tri = delaunayn(x, varargin{:}); |
|
6826
|
54 tri = delaunayn (x); |
|
6823
|
55 |
|
6826
|
56 yi = nan (mi, 1); |
|
6823
|
57 |
|
6826
|
58 if (strcmp (method, "nearest")) |
|
6823
|
59 ## search index of nearest point |
|
6826
|
60 idx = dsearchn (x, tri, xi); |
|
|
61 valid = !isnan (idx); |
|
6823
|
62 yi(valid) = y(idx(valid)); |
|
|
63 |
|
6826
|
64 elseif (strcmp (method, "linear")) |
|
6823
|
65 ## search for every point the enclosing triangle |
|
6826
|
66 [tri_list, bary_list] = tsearchn (x, tri, xi); |
|
6823
|
67 |
|
|
68 ## only keep the points within triangles. |
|
6826
|
69 valid = !isnan (tri_list); |
|
|
70 tri_list = tri_list(!isnan (tri_list)); |
|
|
71 bary_list = bary_list(!isnan (tri_list), :); |
|
|
72 nr_t = rows (tri_list); |
|
6823
|
73 |
|
|
74 ## assign x,y for each point of simplex |
|
6826
|
75 xt = reshape (x(tri(tri_list,:),:), [nr_t, n+1, n]); |
|
6823
|
76 yt = y(tri(tri_list,:)); |
|
|
77 |
|
|
78 ## Use barycentric coordinate of point to calculate yi |
|
|
79 yi(valid) = sum (y(tri(tri_list,:)) .* bary_list, 2); |
|
|
80 |
|
|
81 else |
|
6826
|
82 error ("griddatan: unknown interpolation method"); |
|
6823
|
83 endif |
|
|
84 |
|
|
85 endfunction |
|
|
86 |
|
|
87 %!test |
|
|
88 %! [xx,yy] = meshgrid(linspace(-1,1,32)); |
|
|
89 %! xi = [xx(:), yy(:)]; |
|
|
90 %! x = (2 * rand(100,2) - 1); |
|
|
91 %! x = [x;1,1;1,-1;-1,-1;-1,1]; |
|
|
92 %! y = sin(2*(sum(x.^2,2))); |
|
|
93 %! zz = griddatan(x,y,xi,'linear'); |
|
|
94 %! zz2 = griddata(x(:,1),x(:,2),y,xi(:,1),xi(:,2),'linear'); |
|
|
95 %! assert (zz, zz2, 1e-10) |
|
|
96 |
|
|
97 %!test |
|
|
98 %! [xx,yy] = meshgrid(linspace(-1,1,32)); |
|
|
99 %! xi = [xx(:), yy(:)]; |
|
|
100 %! x = (2 * rand(100,2) - 1); |
|
|
101 %! x = [x;1,1;1,-1;-1,-1;-1,1]; |
|
|
102 %! y = sin(2*(sum(x.^2,2))); |
|
|
103 %! zz = griddatan(x,y,xi,'nearest'); |
|
|
104 %! zz2 = griddata(x(:,1),x(:,2),y,xi(:,1),xi(:,2),'nearest'); |
|
|
105 %! assert (zz, zz2, 1e-10) |
