# HG changeset patch # User Nicholas R. Jankowski # Date 1586463665 14400 # Node ID bb929d5a34cb2ebf83b6a3e2674f60e46e934a62 # Parent 63f5c721cfde6a2510f77659f9ebeb07a3893a33 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. diff -r 63f5c721cfde -r bb929d5a34cb NEWS --- a/NEWS Sat Apr 11 16:26:19 2020 -0700 +++ b/NEWS Thu Apr 09 16:21:05 2020 -0400 @@ -41,8 +41,9 @@ ### Matlab compatibility -- The function `griddata` now accepts 3-D inputs by passing data -directly to `griddata3`. +- The function `griddata` now implements the "v4" Biharmonic Spline +Interpolation method. In adddition, the function now accepts 3-D inputs +by passing the data to `griddata3`. - Coordinate transformation functions `cart2sph`, `sph2cart`, `cart2pol`, and `pol2cart` can now accept either row or column vectors diff -r 63f5c721cfde -r bb929d5a34cb scripts/geometry/griddata.m --- a/scripts/geometry/griddata.m Sat Apr 11 16:26:19 2020 -0700 +++ b/scripts/geometry/griddata.m Thu Apr 09 16:21:05 2020 -0400 @@ -46,8 +46,9 @@ ## are specified by @var{xi}, @var{yi}, @var{zi}. ## ## For both 2-D and 3-D cases, the interpolation method can be -## @qcode{"nearest"} or @qcode{"linear"}. If method is omitted it defaults -## to @qcode{"linear"}. +## @qcode{"linear"} or @qcode{"nearest"}. For 2-D cases only, the @qcode{"v4"} +## method is also available which implements a biharmonic spline interpolation. +## If method is omitted it defaults to @qcode{"linear"}. ## ## For 3-D interpolation, the optional argument @var{options} is passed ## directly to Qhull when computing the Delaunay triangulation used for @@ -122,16 +123,10 @@ tri = delaunay (x, y); zi = NaN (size (xi)); - if (any (strcmp (method, {"cubic", "natural", "v4"}))) + if (any (strcmp (method, {"cubic", "natural"}))) ## FIXME: implement missing interpolation methods. error ('griddata: "%s" interpolation not yet implemented', method); - elseif (strcmp (method, "nearest")) - ## Search index of nearest point. - idx = dsearch (x, y, tri, xi, yi); - valid = ! isnan (idx); - zi(valid) = z(idx(valid)); - elseif (strcmp (method, "linear")) ## Search for every point the enclosing triangle. tri_list = tsearch (x, y, tri, xi(:), yi(:)); @@ -167,6 +162,59 @@ ## Calculate zi by solving plane equation for xi, yi. zi(valid) = -(N(:,1).*xi(:)(valid) + N(:,2).*yi(:)(valid) + D) ./ N(:,3); + elseif (strcmp (method, "nearest")) + ## Search index of nearest point. + idx = dsearch (x, y, tri, xi, yi); + valid = ! isnan (idx); + zi(valid) = z(idx(valid)); + + elseif (strcmp (method, "v4")) + ## Use Biharmonic Spline Interpolation Green's Function method. + ## Compatible with Matlab v4 interpolation method, based on + ## D. Sandwell 1987 and Deng & Tang 2011. + + ## The free space Green Function which solves the two dimensional + ## Biharmonic PDE + ## + ## Delta(Delta(G(X))) = delta(X) + ## + ## for a point source yields + ## + ## G(X) = |X|^2 * (ln|X|-1) / (8 * pi) + ## + ## A N-point Biharmonic Interpolation at the point X is given by + ## + ## z(X) = sum_j_N (alpha_j * G(X-Xj)) + ## = sum_j_N (alpha_j * G(rj)) + ## + ## in which the coefficients alpha_j are the unknowns. rj is the + ## Euclidian distance between X and Xj. + ## From N datapoints {zi, Xi} an equation system can be formed: + ## + ## zi(Xi) = sum_j_N (alpha_j * G(Xi-Xj)) + ## = sum_j_N (alpha_j * G(rij)) + ## + ## Its inverse yields the unknowns alpha_j. + + ## Step1: Solve for weight coefficients alpha_j depending on the + ## Euclidian distances and the training data set {x,y,z} + r = sqrt ((x - x.').^2 + (y - y.').^2); # size N^2 + D = (r.^2) .* (log (r) - 1); + D(isnan (D)) = 0; # Fix Green Function for r=0 + alpha_j = D \ z; + + ## Step2 - Use alphas and Green's functions to get interpolated points. + ## Use dim3 projection for vectorized calculation to avoid loops. Memory + ## usage is proportional to Ni x N. + ## FIXME: if this approach is too memory intensive, revert portion to loop + x = permute (x, [3, 2, 1]); + y = permute (y, [3, 2, 1]); + alpha_j = permute (alpha_j, [3, 2, 1]); + r_i = sqrt ((xi - x).^2 + (yi - y).^2); # size Ni x N + Di = (r_i.^2) .* (log (r_i) - 1); + Di(isnan (Di)) = 0; # Fix Green's Function for r==0 + zi = sum (Di .* alpha_j, 3); + else error ('griddata: unknown interpolation METHOD: "%s"', method); endif @@ -220,7 +268,7 @@ %! 'method = "nearest neighbor"'}); %!testif HAVE_QHULL -%! [xx,yy] = meshgrid (linspace (-1, 1, 32)); +%! [xx, yy] = meshgrid (linspace (-1, 1, 32)); %! x = xx(:); %! x = x + 10*(2*round (rand (size (x))) - 1) * eps; %! y = yy(:); @@ -231,6 +279,30 @@ %! zz2(isnan (zz)) = NaN; %! assert (zz, zz2, 100*eps); +%!testif HAVE_QHULL +%! [xx, yy] = meshgrid (linspace (-1, 1, 5)); +%! x = xx(:); +%! x = x + 10*(2*round (rand (size (x))) - 1) * eps; +%! y = yy(:); +%! y = y + 10*(2*round (rand (size (y))) - 1) * eps; +%! z = 2*(x.^2 + y.^2); +%! zz = griddata (x,y,z,xx,yy, "v4"); +%! zz2 = 2*(xx.^2 + yy.^2); +%! zz2(isnan (zz)) = NaN; +%! assert (zz, zz2, 100*eps); + +%!testif HAVE_QHULL +%! [xx, yy] = meshgrid (linspace (-1, 1, 5)); +%! x = xx(:); +%! x = x + 10*(2*round (rand (size (x))) - 1) * eps; +%! y = yy(:); +%! y = y + 10*(2*round (rand (size (y))) - 1) * eps; +%! z = 2*(x.^2 + y.^2); +%! zz = griddata (x,y,z,xx,yy, "nearest"); +%! zz2 = 2*(xx.^2 + yy.^2); +%! zz2(isnan (zz)) = NaN; +%! assert (zz, zz2, 100*eps); + ## Test input validation %!error griddata () %!error griddata (1) @@ -248,5 +320,4 @@ %!error griddata (1:3, 1:3, 1:3, 1:3, 1:4) %!error <"cubic" .* not yet implemented> griddata (1,2,3,4,5, "cubic") %!error <"natural" .* not yet implemented> griddata (1,2,3,4,5, "natural") -%!error <"v4" .* not yet implemented> griddata (1,2,3,4,5, "v4") %!error griddata (1,2,3,4,5, "foobar") diff -r 63f5c721cfde -r bb929d5a34cb scripts/geometry/griddatan.m --- a/scripts/geometry/griddatan.m Sat Apr 11 16:26:19 2020 -0700 +++ b/scripts/geometry/griddatan.m Thu Apr 09 16:21:05 2020 -0400 @@ -89,8 +89,15 @@ yi(valid) = sum (y(tri(tri_list,:)) .* bary_list, 2); endif + elseif (any (strcmp (method, {"cubic", "v4"}))) + error ('griddata: "%s" METHOD only valid for 2-D inputs using "griddata"', method); + + elseif (strcmp (method, "natural")) + ## FIXME: implement missing interpolation method 'natural' for 3-D inputs. + error ('griddatan: "natural" interpolation METHOD not yet implemented'); + else - error ("griddatan: unknown interpolation METHOD"); + error ('griddatan: unknown interpolation METHOD: "%s"', method); endif endfunction @@ -121,3 +128,16 @@ %! y = [ 1; 2; 3; 4 ]; %! xi = [ .5, .5 ]; %! yi = griddatan (x, y, xi); + +## Test input validation +%!error griddatan () +%!error griddatan (1) +%!error griddatan (1,2) +%!error griddatan (1,2,3) +%!error griddatan (1,2,3,4,5) +%!error griddatan (1,2,3,4) +%!#error <"v4" METHOD only valid for 2-D inputs> +%! griddatan (ones(2,2,2), 2*ones(2,2,2), 3, "v4") +%!error <"cubic" METHOD only valid for> griddatan (1, 2, 3, "cubic") +%!error <"natural" .* not yet implemented> griddatan (1, 2, 3, "natural") +%!error griddatan (1, 2, 3, "foobar")