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