|
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 ## |
|
|
27 ## The interpolation method can be 'nearest' or 'linear'. If method is |
|
|
28 ## omitted it defaults to 'linear'. |
|
|
29 ## @seealso{griddata, delaunayn} |
|
|
30 ## @end deftypefn |
|
|
31 |
|
|
32 function [yi] = griddatan (x,y,xi,method,varargin) |
|
|
33 |
|
|
34 if (nargin == 3) |
|
|
35 method="linear"; |
|
|
36 endif |
|
|
37 if (nargin < 3) |
|
|
38 print_usage(); |
|
|
39 endif |
|
|
40 |
|
|
41 if (ischar(method)) |
|
|
42 method=tolower(method); |
|
|
43 endif |
|
|
44 |
|
|
45 [m, n] = size(x); |
|
|
46 [mi, ni] = size(xi); |
|
|
47 |
|
|
48 if (n != ni || size(y,1) != m || size(y,2) != 1) |
|
|
49 error ("griddatan: dimensional mismatch"); |
|
|
50 endif |
|
|
51 |
|
|
52 ## triangulate data |
|
|
53 ## tri = delaunayn(x, varargin{:}); |
|
|
54 tri = delaunayn(x); |
|
|
55 |
|
|
56 yi = nan(mi, 1); |
|
|
57 |
|
|
58 if strcmp(method,'nearest') |
|
|
59 ## search index of nearest point |
|
|
60 idx = dsearchn(x,tri,xi); |
|
|
61 valid = !isnan(idx); |
|
|
62 yi(valid) = y(idx(valid)); |
|
|
63 |
|
|
64 elseif strcmp(method,'linear') |
|
|
65 ## search for every point the enclosing triangle |
|
|
66 [tri_list, bary_list] = tsearchn(x,tri,xi); |
|
|
67 |
|
|
68 ## only keep the points within triangles. |
|
|
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); |
|
|
73 |
|
|
74 ## assign x,y for each point of simplex |
|
|
75 xt = reshape (x(tri(tri_list,:),:),[nr_t, n+1, n]); |
|
|
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 |
|
|
82 method |
|
|
83 error("griddatan: unknown interpolation method"); |
|
|
84 endif |
|
|
85 |
|
|
86 endfunction |
|
|
87 |
|
|
88 %!test |
|
|
89 %! [xx,yy] = meshgrid(linspace(-1,1,32)); |
|
|
90 %! xi = [xx(:), yy(:)]; |
|
|
91 %! x = (2 * rand(100,2) - 1); |
|
|
92 %! x = [x;1,1;1,-1;-1,-1;-1,1]; |
|
|
93 %! y = sin(2*(sum(x.^2,2))); |
|
|
94 %! zz = griddatan(x,y,xi,'linear'); |
|
|
95 %! zz2 = griddata(x(:,1),x(:,2),y,xi(:,1),xi(:,2),'linear'); |
|
|
96 %! assert (zz, zz2, 1e-10) |
|
|
97 |
|
|
98 %!test |
|
|
99 %! [xx,yy] = meshgrid(linspace(-1,1,32)); |
|
|
100 %! xi = [xx(:), yy(:)]; |
|
|
101 %! x = (2 * rand(100,2) - 1); |
|
|
102 %! x = [x;1,1;1,-1;-1,-1;-1,1]; |
|
|
103 %! y = sin(2*(sum(x.^2,2))); |
|
|
104 %! zz = griddatan(x,y,xi,'nearest'); |
|
|
105 %! zz2 = griddata(x(:,1),x(:,2),y,xi(:,1),xi(:,2),'nearest'); |
|
|
106 %! assert (zz, zz2, 1e-10) |
