|
6823
|
1 ## Copyright (C) 1999,2000 Kai Habel |
|
|
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{zi} =} griddata (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{method}) |
|
|
22 ## @deftypefnx {Function File} {[@var{xi}, @var{yi}, @var{zi}] =} griddata (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{method}) |
|
|
23 ## |
|
|
24 ## Generate a regular mesh from irregular data using interpolation. |
|
|
25 ## The function is defined by @code{@var{z} = f (@var{x}, @var{y})}. |
|
|
26 ## The interpolation points are all @code{(@var{xi}, @var{yi})}. If |
|
|
27 ## @var{xi}, @var{yi} are vectors then they are made into a 2D mesh. |
|
|
28 ## |
|
6826
|
29 ## The interpolation method can be @code{"nearest"}, @code{"cubic"} or |
|
|
30 ## @code{"linear"}. If method is omitted it defaults to @code{"linear"}. |
|
6823
|
31 ## @seealso{delaunay} |
|
|
32 ## @end deftypefn |
|
|
33 |
|
|
34 ## Author: Kai Habel <kai.habel@gmx.de> |
|
|
35 ## Adapted-by: Alexander Barth <barth.alexander@gmail.com> |
|
|
36 ## xi and yi are not "meshgridded" if both are vectors |
|
|
37 ## of the same size (for compatibility) |
|
|
38 |
|
6826
|
39 function [rx, ry, rz] = griddata (x, y, z, xi, yi, method) |
|
6823
|
40 |
|
|
41 if (nargin == 5) |
|
6826
|
42 method = "linear"; |
|
6823
|
43 endif |
|
|
44 if (nargin < 5 || nargin > 7) |
|
6826
|
45 print_usage (); |
|
6823
|
46 endif |
|
|
47 |
|
6826
|
48 if (ischar (method)) |
|
|
49 method = tolower (method); |
|
6823
|
50 endif |
|
6826
|
51 if (! all (size (x) == size (y) & size (x) == size (z))) |
|
|
52 error ("griddata: x, y, and z must be vectors of same length"); |
|
6823
|
53 endif |
|
|
54 |
|
|
55 ## meshgrid xi and yi if they are vectors unless they |
|
|
56 ## are vectors of the same length |
|
6826
|
57 if (isvector (xi) && isvector (yi) && numel (xi) != numel (yi)) |
|
|
58 [xi, yi] = meshgrid (xi, yi); |
|
6823
|
59 endif |
|
|
60 |
|
6826
|
61 if (any (size (xi) != size (yi))) |
|
|
62 error ("griddata: xi and yi must be vectors or matrices of same size"); |
|
6823
|
63 endif |
|
|
64 |
|
6826
|
65 [nr, nc] = size (xi); |
|
6823
|
66 |
|
|
67 ## triangulate data |
|
6826
|
68 tri = delaunay (x, y); |
|
|
69 zi = nan (size (xi)); |
|
6823
|
70 |
|
6826
|
71 if (strcmp (method, "cubic")) |
|
|
72 error ("griddata: cubic interpolation not yet implemented") |
|
6823
|
73 |
|
6826
|
74 elseif (strcmp (method, "nearest")) |
|
6823
|
75 ## search index of nearest point |
|
6826
|
76 idx = dsearch (x, y, tri, xi, yi); |
|
|
77 valid = !isnan (idx); |
|
6823
|
78 zi(valid) = z(idx(valid)); |
|
|
79 |
|
6826
|
80 elseif (strcmp (method, "linear")) |
|
6823
|
81 ## search for every point the enclosing triangle |
|
6826
|
82 tri_list= tsearch (x, y, tri, xi(:), yi(:)); |
|
6823
|
83 |
|
|
84 ## only keep the points within triangles. |
|
6826
|
85 valid = !isnan (reshape (tri_list, size (xi))); |
|
|
86 tri_list = tri_list(!isnan (tri_list)); |
|
|
87 nr_t = rows (tri_list); |
|
6823
|
88 |
|
|
89 ## assign x,y,z for each point of triangle |
|
6826
|
90 x1 = x(tri(tri_list,1)); |
|
|
91 x2 = x(tri(tri_list,2)); |
|
|
92 x3 = x(tri(tri_list,3)); |
|
|
93 |
|
|
94 y1 = y(tri(tri_list,1)); |
|
|
95 y2 = y(tri(tri_list,2)); |
|
|
96 y3 = y(tri(tri_list,3)); |
|
|
97 |
|
|
98 z1 = z(tri(tri_list,1)); |
|
|
99 z2 = z(tri(tri_list,2)); |
|
|
100 z3 = z(tri(tri_list,3)); |
|
6823
|
101 |
|
|
102 ## calculate norm vector |
|
6826
|
103 N = cross ([x2-x1, y2-y1, z2-z1], [x3-x1, y3-y1, z3-z1]); |
|
|
104 N_norm = sqrt (sumsq (N, 2)); |
|
6823
|
105 N = N ./ N_norm(:,[1,1,1]); |
|
|
106 |
|
|
107 ## calculate D of plane equation |
|
|
108 ## Ax+By+Cz+D=0; |
|
|
109 D = -(N(:,1) .* x1 + N(:,2) .* y1 + N(:,3) .* z1); |
|
|
110 |
|
|
111 ## calculate zi by solving plane equation for xi,yi |
|
|
112 zi(valid) = -(N(:,1).*xi(valid) + N(:,2).*yi(valid) + D ) ./ N(:,3); |
|
|
113 |
|
|
114 else |
|
6826
|
115 error ("griddata: unknown interpolation method"); |
|
6823
|
116 endif |
|
|
117 |
|
6826
|
118 if (nargout == 3) |
|
6823
|
119 rx = xi; |
|
|
120 ry = yi; |
|
|
121 rz = zi; |
|
6826
|
122 elseif (nargout == 1) |
|
6823
|
123 rx = zi; |
|
6826
|
124 elseif (nargout == 0) |
|
|
125 mesh (xi, yi, zi); |
|
6823
|
126 endif |
|
|
127 endfunction |
|
|
128 |
|
|
129 %!test |
|
|
130 %! [xx,yy]=meshgrid(linspace(-1,1,32)); |
|
|
131 %! x = xx(:); |
|
|
132 %! x = x + 10 * (2 * round(rand(size(x))) - 1) * eps; |
|
|
133 %! y = yy(:); |
|
|
134 %! y = y + 10 * (2 * round(rand(size(y))) - 1) * eps; |
|
|
135 %! z = sin(2*(x.^2+y.^2)); |
|
|
136 %! zz = griddata(x,y,z,xx,yy,'linear'); |
|
|
137 %! zz2 = sin(2*(xx.^2+yy.^2)); |
|
|
138 %! zz2(isnan(zz)) = NaN; |
|
|
139 %! assert (zz, zz2, 100 * eps) |
|
|
140 |
|
|
141 %!demo |
|
|
142 %! x=2*rand(100,1)-1; |
|
|
143 %! y=2*rand(size(x))-1; |
|
|
144 %! z=sin(2*(x.^2+y.^2)); |
|
|
145 %! [xx,yy]=meshgrid(linspace(-1,1,32)); |
|
|
146 %! griddata(x,y,z,xx,yy); |
|
|
147 %! title('nonuniform grid sampled at 100 points'); |
|
|
148 |
|
|
149 %!demo |
|
|
150 %! x=2*rand(1000,1)-1; |
|
|
151 %! y=2*rand(size(x))-1; |
|
|
152 %! z=sin(2*(x.^2+y.^2)); |
|
|
153 %! [xx,yy]=meshgrid(linspace(-1,1,32)); |
|
|
154 %! griddata(x,y,z,xx,yy); |
|
|
155 %! title('nonuniform grid sampled at 1000 points'); |
|
|
156 |
|
|
157 %!demo |
|
|
158 %! x=2*rand(1000,1)-1; |
|
|
159 %! y=2*rand(size(x))-1; |
|
|
160 %! z=sin(2*(x.^2+y.^2)); |
|
|
161 %! [xx,yy]=meshgrid(linspace(-1,1,32)); |
|
|
162 %! griddata(x,y,z,xx,yy,'nearest'); |
|
|
163 %! title('nonuniform grid sampled at 1000 points with nearest neighbor'); |
