Mercurial > octave-libtiff

comparison scripts/geometry/griddata.m @ 28210:bb929d5a34cb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
griddata.m: Added support for "v4" biharmonic spline interpolation method (bug #33539). * griddata.m: Add v4 method algorithm, documentation notes, and BIST tests. * griddatan.m: Updated method error messages and added BIST tests. * NEWS: Announce support for "v4" method in Matlab compatibility section.
author Nicholas R. Jankowski <jankowskin@asme.org>
date Thu, 09 Apr 2020 16:21:05 -0400
parents 37ee7c5cc6b2
children 0e0e0de09f1e
comparison
equal deleted inserted replaced
28209:63f5c721cfde 28210:bb929d5a34cb
44 ## For 3-D interpolation, the function is defined by 44 ## For 3-D interpolation, the function is defined by
45 ## @code{@var{v} = f (@var{x}, @var{y}, @var{z})}, and the interpolation points 45 ## @code{@var{v} = f (@var{x}, @var{y}, @var{z})}, and the interpolation points
46 ## are specified by @var{xi}, @var{yi}, @var{zi}. 46 ## are specified by @var{xi}, @var{yi}, @var{zi}.
47 ## 47 ##
48 ## For both 2-D and 3-D cases, the interpolation method can be 48 ## For both 2-D and 3-D cases, the interpolation method can be
49 ## @qcode{"nearest"} or @qcode{"linear"}. If method is omitted it defaults 49 ## @qcode{"linear"} or @qcode{"nearest"}. For 2-D cases only, the @qcode{"v4"}
50 ## to @qcode{"linear"}. 50 ## method is also available which implements a biharmonic spline interpolation.
51 ## If method is omitted it defaults to @qcode{"linear"}.
51 ## 52 ##
52 ## For 3-D interpolation, the optional argument @var{options} is passed 53 ## For 3-D interpolation, the optional argument @var{options} is passed
53 ## directly to Qhull when computing the Delaunay triangulation used for 54 ## directly to Qhull when computing the Delaunay triangulation used for
54 ## interpolation. See @code{delaunayn} for information on the defaults and 55 ## interpolation. See @code{delaunayn} for information on the defaults and
55 ## how to pass different values. 56 ## how to pass different values.
120 121
121 ## Triangulate data. 122 ## Triangulate data.
122 tri = delaunay (x, y); 123 tri = delaunay (x, y);
123 zi = NaN (size (xi)); 124 zi = NaN (size (xi));
124 125
125 if (any (strcmp (method, {"cubic", "natural", "v4"}))) 126 if (any (strcmp (method, {"cubic", "natural"})))
126 ## FIXME: implement missing interpolation methods. 127 ## FIXME: implement missing interpolation methods.
127 error ('griddata: "%s" interpolation not yet implemented', method); 128 error ('griddata: "%s" interpolation not yet implemented', method);
129
130 elseif (strcmp (method, "linear"))
131 ## Search for every point the enclosing triangle.
132 tri_list = tsearch (x, y, tri, xi(:), yi(:));
133
134 ## Only keep the points within triangles.
135 valid = ! isnan (tri_list);
136 tri_list = tri_list(valid);
137 nr_t = rows (tri_list);
138
139 tri = tri(tri_list,:);
140
141 ## Assign x,y,z for each point of triangle.
142 x1 = x(tri(:,1));
143 x2 = x(tri(:,2));
144 x3 = x(tri(:,3));
145
146 y1 = y(tri(:,1));
147 y2 = y(tri(:,2));
148 y3 = y(tri(:,3));
149
150 z1 = z(tri(:,1));
151 z2 = z(tri(:,2));
152 z3 = z(tri(:,3));
153
154 ## Calculate norm vector.
155 N = cross ([x2-x1, y2-y1, z2-z1], [x3-x1, y3-y1, z3-z1]);
156 ## Normalize.
157 N = diag (norm (N, "rows")) \ N;
158
159 ## Calculate D of plane equation: Ax+By+Cz+D = 0
160 D = -(N(:,1) .* x1 + N(:,2) .* y1 + N(:,3) .* z1);
161
162 ## Calculate zi by solving plane equation for xi, yi.
163 zi(valid) = -(N(:,1).*xi(:)(valid) + N(:,2).*yi(:)(valid) + D) ./ N(:,3);
128 164
129 elseif (strcmp (method, "nearest")) 165 elseif (strcmp (method, "nearest"))
130 ## Search index of nearest point. 166 ## Search index of nearest point.
131 idx = dsearch (x, y, tri, xi, yi); 167 idx = dsearch (x, y, tri, xi, yi);
132 valid = ! isnan (idx); 168 valid = ! isnan (idx);
133 zi(valid) = z(idx(valid)); 169 zi(valid) = z(idx(valid));
134 170
135 elseif (strcmp (method, "linear")) 171 elseif (strcmp (method, "v4"))
136 ## Search for every point the enclosing triangle. 172 ## Use Biharmonic Spline Interpolation Green's Function method.
137 tri_list = tsearch (x, y, tri, xi(:), yi(:)); 173 ## Compatible with Matlab v4 interpolation method, based on
138 174 ## D. Sandwell 1987 and Deng & Tang 2011.
139 ## Only keep the points within triangles. 175
140 valid = ! isnan (tri_list); 176 ## The free space Green Function which solves the two dimensional
141 tri_list = tri_list(valid); 177 ## Biharmonic PDE
142 nr_t = rows (tri_list); 178 ##
143 179 ## Delta(Delta(G(X))) = delta(X)
144 tri = tri(tri_list,:); 180 ##
145 181 ## for a point source yields
146 ## Assign x,y,z for each point of triangle. 182 ##
147 x1 = x(tri(:,1)); 183 ## G(X) = |X|^2 * (ln|X|-1) / (8 * pi)
148 x2 = x(tri(:,2)); 184 ##
149 x3 = x(tri(:,3)); 185 ## A N-point Biharmonic Interpolation at the point X is given by
150 186 ##
151 y1 = y(tri(:,1)); 187 ## z(X) = sum_j_N (alpha_j * G(X-Xj))
152 y2 = y(tri(:,2)); 188 ## = sum_j_N (alpha_j * G(rj))
153 y3 = y(tri(:,3)); 189 ##
154 190 ## in which the coefficients alpha_j are the unknowns. rj is the
155 z1 = z(tri(:,1)); 191 ## Euclidian distance between X and Xj.
156 z2 = z(tri(:,2)); 192 ## From N datapoints {zi, Xi} an equation system can be formed:
157 z3 = z(tri(:,3)); 193 ##
158 194 ## zi(Xi) = sum_j_N (alpha_j * G(Xi-Xj))
159 ## Calculate norm vector. 195 ## = sum_j_N (alpha_j * G(rij))
160 N = cross ([x2-x1, y2-y1, z2-z1], [x3-x1, y3-y1, z3-z1]); 196 ##
161 ## Normalize. 197 ## Its inverse yields the unknowns alpha_j.
162 N = diag (norm (N, "rows")) \ N; 198
163 199 ## Step1: Solve for weight coefficients alpha_j depending on the
164 ## Calculate D of plane equation: Ax+By+Cz+D = 0 200 ## Euclidian distances and the training data set {x,y,z}
165 D = -(N(:,1) .* x1 + N(:,2) .* y1 + N(:,3) .* z1); 201 r = sqrt ((x - x.').^2 + (y - y.').^2); # size N^2
166 202 D = (r.^2) .* (log (r) - 1);
167 ## Calculate zi by solving plane equation for xi, yi. 203 D(isnan (D)) = 0; # Fix Green Function for r=0
168 zi(valid) = -(N(:,1).*xi(:)(valid) + N(:,2).*yi(:)(valid) + D) ./ N(:,3); 204 alpha_j = D \ z;
205
206 ## Step2 - Use alphas and Green's functions to get interpolated points.
207 ## Use dim3 projection for vectorized calculation to avoid loops. Memory
208 ## usage is proportional to Ni x N.
209 ## FIXME: if this approach is too memory intensive, revert portion to loop
210 x = permute (x, [3, 2, 1]);
211 y = permute (y, [3, 2, 1]);
212 alpha_j = permute (alpha_j, [3, 2, 1]);
213 r_i = sqrt ((xi - x).^2 + (yi - y).^2); # size Ni x N
214 Di = (r_i.^2) .* (log (r_i) - 1);
215 Di(isnan (Di)) = 0; # Fix Green's Function for r==0
216 zi = sum (Di .* alpha_j, 3);
169 217
170 else 218 else
171 error ('griddata: unknown interpolation METHOD: "%s"', method); 219 error ('griddata: unknown interpolation METHOD: "%s"', method);
172 endif 220 endif
173 221
218 %! mesh (xx, yy, zz); 266 %! mesh (xx, yy, zz);
219 %! title ({"non-uniform grid sampled at 1,000 points", 267 %! title ({"non-uniform grid sampled at 1,000 points",
220 %! 'method = "nearest neighbor"'}); 268 %! 'method = "nearest neighbor"'});
221 269
222 %!testif HAVE_QHULL 270 %!testif HAVE_QHULL
223 %! [xx,yy] = meshgrid (linspace (-1, 1, 32)); 271 %! [xx, yy] = meshgrid (linspace (-1, 1, 32));
224 %! x = xx(:); 272 %! x = xx(:);
225 %! x = x + 10*(2*round (rand (size (x))) - 1) * eps; 273 %! x = x + 10*(2*round (rand (size (x))) - 1) * eps;
226 %! y = yy(:); 274 %! y = yy(:);
227 %! y = y + 10*(2*round (rand (size (y))) - 1) * eps; 275 %! y = y + 10*(2*round (rand (size (y))) - 1) * eps;
228 %! z = sin (2*(x.^2 + y.^2)); 276 %! z = sin (2*(x.^2 + y.^2));
229 %! zz = griddata (x,y,z,xx,yy, "linear"); 277 %! zz = griddata (x,y,z,xx,yy, "linear");
230 %! zz2 = sin (2*(xx.^2 + yy.^2)); 278 %! zz2 = sin (2*(xx.^2 + yy.^2));
279 %! zz2(isnan (zz)) = NaN;
280 %! assert (zz, zz2, 100*eps);
281
282 %!testif HAVE_QHULL
283 %! [xx, yy] = meshgrid (linspace (-1, 1, 5));
284 %! x = xx(:);
285 %! x = x + 10*(2*round (rand (size (x))) - 1) * eps;
286 %! y = yy(:);
287 %! y = y + 10*(2*round (rand (size (y))) - 1) * eps;
288 %! z = 2*(x.^2 + y.^2);
289 %! zz = griddata (x,y,z,xx,yy, "v4");
290 %! zz2 = 2*(xx.^2 + yy.^2);
291 %! zz2(isnan (zz)) = NaN;
292 %! assert (zz, zz2, 100*eps);
293
294 %!testif HAVE_QHULL
295 %! [xx, yy] = meshgrid (linspace (-1, 1, 5));
296 %! x = xx(:);
297 %! x = x + 10*(2*round (rand (size (x))) - 1) * eps;
298 %! y = yy(:);
299 %! y = y + 10*(2*round (rand (size (y))) - 1) * eps;
300 %! z = 2*(x.^2 + y.^2);
301 %! zz = griddata (x,y,z,xx,yy, "nearest");
302 %! zz2 = 2*(xx.^2 + yy.^2);
231 %! zz2(isnan (zz)) = NaN; 303 %! zz2(isnan (zz)) = NaN;
232 %! assert (zz, zz2, 100*eps); 304 %! assert (zz, zz2, 100*eps);
233 305
234 ## Test input validation 306 ## Test input validation
235 %!error griddata () 307 %!error griddata ()
246 %!error <the columns and rows of Z> griddata (1:4, 1:3, ones (3,5), 1:3, 1:3) 318 %!error <the columns and rows of Z> griddata (1:4, 1:3, ones (3,5), 1:3, 1:3)
247 %!error <XI and YI .* matrices of same size> griddata (1:3, 1:3, 1:3, 1:4, 1:3) 319 %!error <XI and YI .* matrices of same size> griddata (1:3, 1:3, 1:3, 1:4, 1:3)
248 %!error <XI and YI .* matrices of same size> griddata (1:3, 1:3, 1:3, 1:3, 1:4) 320 %!error <XI and YI .* matrices of same size> griddata (1:3, 1:3, 1:3, 1:3, 1:4)
249 %!error <"cubic" .* not yet implemented> griddata (1,2,3,4,5, "cubic") 321 %!error <"cubic" .* not yet implemented> griddata (1,2,3,4,5, "cubic")
250 %!error <"natural" .* not yet implemented> griddata (1,2,3,4,5, "natural") 322 %!error <"natural" .* not yet implemented> griddata (1,2,3,4,5, "natural")
251 %!error <"v4" .* not yet implemented> griddata (1,2,3,4,5, "v4")
252 %!error <unknown interpolation METHOD: "foobar"> griddata (1,2,3,4,5, "foobar") 323 %!error <unknown interpolation METHOD: "foobar"> griddata (1,2,3,4,5, "foobar")