Mercurial > octave

comparison scripts/geometry/griddata.m @ 28214:0e0e0de09f1e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
griddata.m: Overhaul function. * griddata.m: Rewrite documentation for clarity. Place all input validation before calculations. Validate METHOD input more precisely. Don't calculate Delaunay triangulation for "v4" method as it is unnecessary. Update BIST tests.
author Rik <rik@octave.org>
date Mon, 13 Apr 2020 18:07:28 -0700
parents bb929d5a34cb
children f5644ccd1df5
comparison
equal deleted inserted replaced
28213:bb50407f9c60 28214:0e0e0de09f1e
29 ## @deftypefnx {} {[@var{xi}, @var{yi}, @var{zi}] =} griddata (@dots{}) 29 ## @deftypefnx {} {[@var{xi}, @var{yi}, @var{zi}] =} griddata (@dots{})
30 ## @deftypefnx {} {@var{vi} =} griddata (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi}) 30 ## @deftypefnx {} {@var{vi} =} griddata (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi})
31 ## @deftypefnx {} {@var{vi} =} griddata (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi}, @var{method}) 31 ## @deftypefnx {} {@var{vi} =} griddata (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi}, @var{method})
32 ## @deftypefnx {} {@var{vi} =} griddata (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi}, @var{method}, @var{options}) 32 ## @deftypefnx {} {@var{vi} =} griddata (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi}, @var{method}, @var{options})
33 ## 33 ##
34 ## Generate a regular mesh from irregular data using 2-D or 3-D interpolation. 34 ## Interpolate irregular 2-D and 3-D source data at specified points.
35 ## 35 ##
36 ## For 2-D interpolation, the function is defined by 36 ## For 2-D interpolation, the inputs @var{x} and @var{y} define the points
37 ## @code{@var{z} = f (@var{x}, @var{y})}. Inputs 37 ## where the function @code{@var{z} = f (@var{x}, @var{y})} is evaluated.
38 ## @code{@var{x}, @var{y}, @var{z}} are vectors of the same length or 38 ## The inputs @var{x}, @var{y}, @var{z} are either vectors of the same length,
39 ## @code{@var{x}, @var{y}} are vectors and @code{@var{z}} is a matrix. 39 ## or the unequal vectors @var{x}, @var{y} are expanded to a 2-D grid with
40 ## 40 ## @code{meshgrid} and @var{z} is a 2-D matrix matching the resulting size of
41 ## The interpolation points are all @code{(@var{xi}, @var{yi})}. If @var{xi}, 41 ## the X-Y grid.
42 ## @var{yi} are vectors then they are made into a 2-D mesh. 42 ##
43 ## 43 ## The interpolation points are (@var{xi}, @var{yi}). If, and only if,
44 ## For 3-D interpolation, the function is defined by 44 ## @var{xi} is a row vector and @var{yi} is a column vector, then
45 ## @code{@var{v} = f (@var{x}, @var{y}, @var{z})}, and the interpolation points 45 ## @code{meshgrid} will be used to create a mesh of interpolation points.
46 ## are specified by @var{xi}, @var{yi}, @var{zi}. 46 ##
47 ## 47 ## For 3-D interpolation, the inputs @var{x}, @var{y}, and @var{z} define the
48 ## For both 2-D and 3-D cases, the interpolation method can be 48 ## points where the function @code{@var{v} = f (@var{x}, @var{y}, @var{z})}
49 ## @qcode{"linear"} or @qcode{"nearest"}. For 2-D cases only, the @qcode{"v4"} 49 ## is evaluated. The inputs @var{x}, @var{y}, @var{z} are either vectors of
50 ## method is also available which implements a biharmonic spline interpolation. 50 ## the same length, or if they are of unequal length, then they are expanded to
51 ## If method is omitted it defaults to @qcode{"linear"}. 51 ## a 3-D grid with @code{meshgrid}. The size of the input @var{v} must match
52 ## the size of the original data, either as a vector or a matrix.
53 ##
54 ## The optional input interpolation @var{method} can be @qcode{"nearest"},
55 ## @qcode{"linear"}, or for 2-D data @qcode{"v4"}. When the method is
56 ## @qcode{"nearest"}, the output @var{vi} will be the closest point in the
57 ## original data (@var{x}, @var{y}, @var{z}) to the query point (@var{xi},
58 ## @var{yi}, @var{zi}). When the method is @qcode{"linear"}, the output
59 ## @var{vi} will be a linear interpolation between the two closest points in
60 ## the original source data in each dimension. For 2-D cases only, the
61 ## @qcode{"v4"} method is also available which implements a biharmonic spline
62 ## interpolation. If @var{method} is omitted or empty, it defaults to
63 ## @qcode{"linear"}.
52 ## 64 ##
53 ## For 3-D interpolation, the optional argument @var{options} is passed 65 ## For 3-D interpolation, the optional argument @var{options} is passed
54 ## directly to Qhull when computing the Delaunay triangulation used for 66 ## directly to Qhull when computing the Delaunay triangulation used for
55 ## interpolation. See @code{delaunayn} for information on the defaults and 67 ## interpolation. See @code{delaunayn} for information on the defaults and
56 ## how to pass different values. 68 ## how to pass different values.
69 ##
70 ## Programming Notes: If the input is complex the real and imaginary parts
71 ## are interpolated separately. Interpolation is normally based on a
72 ## Delaunay triangulation. Any query values outside the convex hull of the
73 ## input points will return @code{NaN}. However, the @qcode{"v4"} method does
74 ## not use the triangulation and will return values outside the original data
75 ## (extrapolation).
57 ## @seealso{griddata3, griddatan, delaunay} 76 ## @seealso{griddata3, griddatan, delaunay}
58 ## @end deftypefn 77 ## @end deftypefn
59
60 ## Algorithm: xi and yi are not "meshgridded" if both are vectors
61 ## of the same size (for compatibility)
62 78
63 function [rx, ry, rz] = griddata (x, y, z, varargin) 79 function [rx, ry, rz] = griddata (x, y, z, varargin)
64 80
65 if (nargin < 5) 81 if (nargin < 5)
66 print_usage (); 82 print_usage ();
69 if (nargin > 6) 85 if (nargin > 6)
70 ## Current 2D implementation has nargin max = 6, since no triangulation 86 ## Current 2D implementation has nargin max = 6, since no triangulation
71 ## options are passed to the 2D algorithm. 3D algorithm requires nargin >=7 87 ## options are passed to the 2D algorithm. 3D algorithm requires nargin >=7
72 88
73 if (nargout > 1) 89 if (nargout > 1)
74 error ("griddata: only one output argument valid for 3D interpolation"); 90 error ("griddata: only one output argument valid for 3-D interpolation");
75 endif 91 endif
76 rx = griddata3 (x, y, z, varargin{:}); 92 rx = griddata3 (x, y, z, varargin{:});
77 93
78 else 94 else
79 ## for nargin 5 or 6, assign varargin terms to variables for 2D algorithm 95 ## for nargin 5 or 6, assign varargin terms to variables for 2D algorithm
80 xi = varargin{1}; 96 xi = varargin{1};
81 yi = varargin{2}; 97 yi = varargin{2};
82
83 if (nargin == 6)
84 if (! ischar (varargin{3}))
85 error ("griddata: unknown interpolation METHOD");
86 endif
87 method = tolower (varargin{3});
88 else
89 method = "linear";
90 endif
91 98
92 ## Meshgrid if x and y are vectors but z is matrix 99 ## Meshgrid if x and y are vectors but z is matrix
93 if (isvector (x) && isvector (y) && all ([numel(y), numel(x)] == size (z))) 100 if (isvector (x) && isvector (y) && all ([numel(y), numel(x)] == size (z)))
94 [x, y] = meshgrid (x, y); 101 [x, y] = meshgrid (x, y);
95 endif 102 endif
100 endif 107 endif
101 elseif (! size_equal (x, y, z)) 108 elseif (! size_equal (x, y, z))
102 error ("griddata: lengths of X, Y must match the columns and rows of Z"); 109 error ("griddata: lengths of X, Y must match the columns and rows of Z");
103 endif 110 endif
104 111
105 ## Meshgrid xi and yi if they are a row and column vector. 112 ## Meshgrid xi and yi if they are a row and column vector, but not
106 if (rows (xi) == 1 && columns (yi) == 1) 113 ## if they are simply vectors of the same size (for compatibility).
114 if (isrow (xi) && iscolumn (yi))
107 [xi, yi] = meshgrid (xi, yi); 115 [xi, yi] = meshgrid (xi, yi);
108 elseif (isvector (xi) && isvector (yi)) 116 elseif (isvector (xi) && isvector (yi))
109 ## Otherwise, convert to column vectors 117 ## Otherwise, convert to column vectors
110 xi = xi(:); 118 xi = xi(:);
111 yi = yi(:); 119 yi = yi(:);
113 121
114 if (! size_equal (xi, yi)) 122 if (! size_equal (xi, yi))
115 error ("griddata: XI and YI must be vectors or matrices of same size"); 123 error ("griddata: XI and YI must be vectors or matrices of same size");
116 endif 124 endif
117 125
126 if (nargin == 6)
127 method = varargin{3};
128 if (isempty (method))
129 method = "linear";
130 elseif (! ischar (method))
131 error ("griddata: METHOD must be a string");
132 endif
133 method = tolower (method);
134
135 if (any (strcmp (method, {"linear", "nearest", "v4"})))
136 ## Do nothing, these are implemented methods
137 elseif (any (strcmp (method, {"cubic", "natural"})))
138 ## FIXME: implement missing interpolation methods.
139 error ('griddata: "%s" interpolation not yet implemented', method);
140 else
141 error ('griddata: unknown interpolation METHOD: "%s"', method);
142 endif
143 else
144 method = "linear";
145 endif
146
118 x = x(:); 147 x = x(:);
119 y = y(:); 148 y = y(:);
120 z = z(:); 149 z = z(:);
121 150
122 ## Triangulate data. 151 ## Triangulate data.
123 tri = delaunay (x, y); 152 if (! strcmp (method, "v4"))
153 tri = delaunay (x, y);
154 endif
124 zi = NaN (size (xi)); 155 zi = NaN (size (xi));
125 156
126 if (any (strcmp (method, {"cubic", "natural"}))) 157 if (strcmp (method, "linear"))
127 ## FIXME: implement missing interpolation methods.
128 error ('griddata: "%s" interpolation not yet implemented', method);
129
130 elseif (strcmp (method, "linear"))
131 ## Search for every point the enclosing triangle. 158 ## Search for every point the enclosing triangle.
132 tri_list = tsearch (x, y, tri, xi(:), yi(:)); 159 tri_list = tsearch (x, y, tri, xi(:), yi(:));
133 160
134 ## Only keep the points within triangles. 161 ## Only keep the points within triangles.
135 valid = ! isnan (tri_list); 162 valid = ! isnan (tri_list);
180 ## 207 ##
181 ## for a point source yields 208 ## for a point source yields
182 ## 209 ##
183 ## G(X) = |X|^2 * (ln|X|-1) / (8 * pi) 210 ## G(X) = |X|^2 * (ln|X|-1) / (8 * pi)
184 ## 211 ##
185 ## A N-point Biharmonic Interpolation at the point X is given by 212 ## An N-point Biharmonic Interpolation at the point X is given by
186 ## 213 ##
187 ## z(X) = sum_j_N (alpha_j * G(X-Xj)) 214 ## z(X) = sum_j_N (alpha_j * G(X-Xj))
188 ## = sum_j_N (alpha_j * G(rj)) 215 ## = sum_j_N (alpha_j * G(rj))
189 ## 216 ##
190 ## in which the coefficients alpha_j are the unknowns. rj is the 217 ## in which the coefficients alpha_j are the unknowns. rj is the
191 ## Euclidian distance between X and Xj. 218 ## Euclidean distance between X and Xj.
192 ## From N datapoints {zi, Xi} an equation system can be formed: 219 ## From N datapoints {zi, Xi} an equation system can be formed:
193 ## 220 ##
194 ## zi(Xi) = sum_j_N (alpha_j * G(Xi-Xj)) 221 ## zi(Xi) = sum_j_N (alpha_j * G(Xi-Xj))
195 ## = sum_j_N (alpha_j * G(rij)) 222 ## = sum_j_N (alpha_j * G(rij))
196 ## 223 ##
197 ## Its inverse yields the unknowns alpha_j. 224 ## Its inverse yields the unknowns alpha_j.
198 225
199 ## Step1: Solve for weight coefficients alpha_j depending on the 226 ## Step1: Solve for weight coefficients alpha_j depending on the
200 ## Euclidian distances and the training data set {x,y,z} 227 ## Euclidean distances and the training data set {x,y,z}
201 r = sqrt ((x - x.').^2 + (y - y.').^2); # size N^2 228 r = sqrt ((x - x.').^2 + (y - y.').^2); # size N^2
202 D = (r.^2) .* (log (r) - 1); 229 D = (r.^2) .* (log (r) - 1);
203 D(isnan (D)) = 0; # Fix Green Function for r=0 230 D(isnan (D)) = 0; # Fix Green Function for r=0
204 alpha_j = D \ z; 231 alpha_j = D \ z;
205 232
206 ## Step2 - Use alphas and Green's functions to get interpolated points. 233 ## Step2 - Use alphas and Green's functions to get interpolated points.
207 ## Use dim3 projection for vectorized calculation to avoid loops. Memory 234 ## Use dim3 projection for vectorized calculation to avoid loops.
208 ## usage is proportional to Ni x N. 235 ## Memory usage is proportional to Ni x N.
209 ## FIXME: if this approach is too memory intensive, revert portion to loop 236 ## FIXME: if this approach is too memory intensive, revert portion to loop
210 x = permute (x, [3, 2, 1]); 237 x = permute (x, [3, 2, 1]);
211 y = permute (y, [3, 2, 1]); 238 y = permute (y, [3, 2, 1]);
212 alpha_j = permute (alpha_j, [3, 2, 1]); 239 alpha_j = permute (alpha_j, [3, 2, 1]);
213 r_i = sqrt ((xi - x).^2 + (yi - y).^2); # size Ni x N 240 r_i = sqrt ((xi - x).^2 + (yi - y).^2); # size Ni x N
214 Di = (r_i.^2) .* (log (r_i) - 1); 241 Di = (r_i.^2) .* (log (r_i) - 1);
215 Di(isnan (Di)) = 0; # Fix Green's Function for r==0 242 Di(isnan (Di)) = 0; # Fix Green's Function for r==0
216 zi = sum (Di .* alpha_j, 3); 243 zi = sum (Di .* alpha_j, 3);
217 244
218 else
219 error ('griddata: unknown interpolation METHOD: "%s"', method);
220 endif 245 endif
221 246
222 if (nargout > 1) 247 if (nargout > 1)
223 rx = xi; 248 rx = xi;
224 ry = yi; 249 ry = yi;
308 %!error griddata (1) 333 %!error griddata (1)
309 %!error griddata (1,2) 334 %!error griddata (1,2)
310 %!error griddata (1,2,3) 335 %!error griddata (1,2,3)
311 %!error griddata (1,2,3,4) 336 %!error griddata (1,2,3,4)
312 %!error <only one output argument> [xi,yi] = griddata (1,2,3,4,5,6,7) 337 %!error <only one output argument> [xi,yi] = griddata (1,2,3,4,5,6,7)
313 %!error <unknown interpolation METHOD> griddata (1,2,3,4,5, {"linear"})
314 %!error <vectors of the same length> griddata (1:4, 1:3, 1:3, 1:3, 1:3) 338 %!error <vectors of the same length> griddata (1:4, 1:3, 1:3, 1:3, 1:3)
315 %!error <vectors of the same length> griddata (1:3, 1:4, 1:3, 1:3, 1:3) 339 %!error <vectors of the same length> griddata (1:3, 1:4, 1:3, 1:3, 1:3)
316 %!error <vectors of the same length> griddata (1:3, 1:3, 1:4, 1:3, 1:3) 340 %!error <vectors of the same length> griddata (1:3, 1:3, 1:4, 1:3, 1:3)
317 %!error <the columns and rows of Z> griddata (1:4, 1:3, ones (4,4), 1:3, 1:3) 341 %!error <the columns and rows of Z> griddata (1:4, 1:3, ones (4,4), 1:3, 1:3)
318 %!error <the columns and rows of Z> griddata (1:4, 1:3, ones (3,5), 1:3, 1:3) 342 %!error <the columns and rows of Z> griddata (1:4, 1:3, ones (3,5), 1:3, 1:3)
319 %!error <XI and YI .* matrices of same size> griddata (1:3, 1:3, 1:3, 1:4, 1:3) 343 %!error <XI and YI .* matrices of same size> griddata (1:3, 1:3, 1:3, 1:4, 1:3)
320 %!error <XI and YI .* matrices of same size> griddata (1:3, 1:3, 1:3, 1:3, 1:4) 344 %!error <XI and YI .* matrices of same size> griddata (1:3, 1:3, 1:3, 1:3, 1:4)
345 %!error <METHOD must be a string> griddata (1,2,3,4,5, {"linear"})
321 %!error <"cubic" .* not yet implemented> griddata (1,2,3,4,5, "cubic") 346 %!error <"cubic" .* not yet implemented> griddata (1,2,3,4,5, "cubic")
322 %!error <"natural" .* not yet implemented> griddata (1,2,3,4,5, "natural") 347 %!error <"natural" .* not yet implemented> griddata (1,2,3,4,5, "natural")
323 %!error <unknown interpolation METHOD: "foobar"> griddata (1,2,3,4,5, "foobar") 348 %!error <unknown interpolation METHOD: "foobar"> griddata (1,2,3,4,5, "foobar")