# HG changeset patch # User Nicholas R. Jankowski # Date 1581907033 18000 # Node ID 37ee7c5cc6b2584053c2ada4d235908866bda1f3 # Parent 38a9f6444eb015714ca13ef043e657e403bb0c24 griddata.m: Add support for 3D inputs (bug #57323). * NEWS: Announce change in Matlab compatibility section. * griddata.m: Rewrite docstring to add new calling forms and document them. Detect 3-D input in input validation and call griddata3. Update error message strings. Update and add new BIST tests for input validation. diff -r 38a9f6444eb0 -r 37ee7c5cc6b2 NEWS --- a/NEWS Tue Mar 10 07:51:55 2020 -0400 +++ b/NEWS Sun Feb 16 21:37:13 2020 -0500 @@ -41,6 +41,8 @@ ### Matlab compatibility +- The function `griddata` now accepts 3-D inputs by passing data +directly to `griddata3`. ### Alphabetical list of new functions added in Octave 7 diff -r 38a9f6444eb0 -r 37ee7c5cc6b2 scripts/geometry/griddata.m --- a/scripts/geometry/griddata.m Tue Mar 10 07:51:55 2020 -0400 +++ b/scripts/geometry/griddata.m Sun Feb 16 21:37:13 2020 -0500 @@ -27,123 +27,158 @@ ## @deftypefn {} {@var{zi} =} griddata (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}) ## @deftypefnx {} {@var{zi} =} griddata (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{method}) ## @deftypefnx {} {[@var{xi}, @var{yi}, @var{zi}] =} griddata (@dots{}) -## -## Generate a regular mesh from irregular data using interpolation. +## @deftypefnx {} {@var{vi} =} griddata (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi}) +## @deftypefnx {} {@var{vi} =} griddata (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi}, @var{method}) +## @deftypefnx {} {@var{vi} =} griddata (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi}, @var{method}, @var{options}) ## -## The function is defined by @code{@var{z} = f (@var{x}, @var{y})}. Inputs +## Generate a regular mesh from irregular data using 2-D or 3-D interpolation. +## +## For 2-D interpolation, the function is defined by +## @code{@var{z} = f (@var{x}, @var{y})}. Inputs ## @code{@var{x}, @var{y}, @var{z}} are vectors of the same length or -## @code{@var{x}, @var{y}} are vectors and @code{@var{z}} is matrix. +## @code{@var{x}, @var{y}} are vectors and @code{@var{z}} is a matrix. ## ## The interpolation points are all @code{(@var{xi}, @var{yi})}. If @var{xi}, ## @var{yi} are vectors then they are made into a 2-D mesh. ## -## The interpolation method can be @qcode{"nearest"}, @qcode{"cubic"} or -## @qcode{"linear"}. If method is omitted it defaults to @qcode{"linear"}. +## For 3-D interpolation, the function is defined by +## @code{@var{v} = f (@var{x}, @var{y}, @var{z})}, and the interpolation points +## 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"}. +## +## For 3-D interpolation, the optional argument @var{options} is passed +## directly to Qhull when computing the Delaunay triangulation used for +## interpolation. See @code{delaunayn} for information on the defaults and +## how to pass different values. ## @seealso{griddata3, griddatan, delaunay} ## @end deftypefn ## Algorithm: xi and yi are not "meshgridded" if both are vectors ## of the same size (for compatibility) -function [rx, ry, rz] = griddata (x, y, z, xi, yi, method = "linear") +function [rx, ry, rz] = griddata (x, y, z, varargin) - if (nargin < 5 || nargin > 7) + if (nargin < 5) print_usage (); endif - if (ischar (method)) - method = tolower (method); - endif - - ## Meshgrid if x and y are vectors but z is matrix - if (isvector (x) && isvector (y) && all ([numel(y), numel(x)] == size (z))) - [x, y] = meshgrid (x, y); - endif - - if (isvector (x) && isvector (y) && isvector (z)) - if (! isequal (length (x), length (y), length (z))) - error ("griddata: X, Y, and Z must be vectors of the same length"); - endif - elseif (! size_equal (x, y, z)) - error ("griddata: lengths of X, Y must match the columns and rows of Z"); - endif - - ## Meshgrid xi and yi if they are a row and column vector. - if (rows (xi) == 1 && columns (yi) == 1) - [xi, yi] = meshgrid (xi, yi); - elseif (isvector (xi) && isvector (yi)) - ## Otherwise, convert to column vectors - xi = xi(:); - yi = yi(:); - endif - - if (! size_equal (xi, yi)) - error ("griddata: XI and YI must be vectors or matrices of same size"); - endif - - x = x(:); - y = y(:); - z = z(:); - - ## Triangulate data. - tri = delaunay (x, y); - zi = NaN (size (xi)); - - if (strcmp (method, "cubic")) - error ("griddata: cubic interpolation not yet implemented"); + if (nargin > 6) + ## Current 2D implementation has nargin max = 6, since no triangulation + ## options are passed to the 2D algorithm. 3D algorithm requires nargin >=7 - 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(:)); - - ## Only keep the points within triangles. - valid = ! isnan (tri_list); - tri_list = tri_list(valid); - nr_t = rows (tri_list); - - tri = tri(tri_list,:); - - ## Assign x,y,z for each point of triangle. - x1 = x(tri(:,1)); - x2 = x(tri(:,2)); - x3 = x(tri(:,3)); - - y1 = y(tri(:,1)); - y2 = y(tri(:,2)); - y3 = y(tri(:,3)); - - z1 = z(tri(:,1)); - z2 = z(tri(:,2)); - z3 = z(tri(:,3)); - - ## Calculate norm vector. - N = cross ([x2-x1, y2-y1, z2-z1], [x3-x1, y3-y1, z3-z1]); - ## Normalize. - N = diag (norm (N, "rows")) \ N; - - ## Calculate D of plane equation - ## Ax+By+Cz+D = 0; - D = -(N(:,1) .* x1 + N(:,2) .* y1 + N(:,3) .* z1); - - ## Calculate zi by solving plane equation for xi, yi. - zi(valid) = -(N(:,1).*xi(:)(valid) + N(:,2).*yi(:)(valid) + D) ./ N(:,3); + if (nargout > 1) + error ("griddata: only one output argument valid for 3D interpolation"); + endif + rx = griddata3 (x, y, z, varargin{:}); else - error ("griddata: unknown interpolation METHOD"); - endif + ## for nargin 5 or 6, assign varargin terms to variables for 2D algorithm + xi = varargin{1}; + yi = varargin{2}; + + if (nargin == 6) + if (! ischar (varargin{3})) + error ("griddata: unknown interpolation METHOD"); + endif + method = tolower (varargin{3}); + else + method = "linear"; + endif + + ## Meshgrid if x and y are vectors but z is matrix + if (isvector (x) && isvector (y) && all ([numel(y), numel(x)] == size (z))) + [x, y] = meshgrid (x, y); + endif + + if (isvector (x) && isvector (y) && isvector (z)) + if (! isequal (length (x), length (y), length (z))) + error ("griddata: X, Y, and Z must be vectors of the same length"); + endif + elseif (! size_equal (x, y, z)) + error ("griddata: lengths of X, Y must match the columns and rows of Z"); + endif + + ## Meshgrid xi and yi if they are a row and column vector. + if (rows (xi) == 1 && columns (yi) == 1) + [xi, yi] = meshgrid (xi, yi); + elseif (isvector (xi) && isvector (yi)) + ## Otherwise, convert to column vectors + xi = xi(:); + yi = yi(:); + endif + + if (! size_equal (xi, yi)) + error ("griddata: XI and YI must be vectors or matrices of same size"); + endif + + x = x(:); + y = y(:); + z = z(:); + + ## Triangulate data. + tri = delaunay (x, y); + zi = NaN (size (xi)); + + if (any (strcmp (method, {"cubic", "natural", "v4"}))) + ## FIXME: implement missing interpolation methods. + error ('griddata: "%s" interpolation not yet implemented', method); - if (nargout > 1) - rx = xi; - ry = yi; - rz = zi; - else - rx = zi; + 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(:)); + + ## Only keep the points within triangles. + valid = ! isnan (tri_list); + tri_list = tri_list(valid); + nr_t = rows (tri_list); + + tri = tri(tri_list,:); + + ## Assign x,y,z for each point of triangle. + x1 = x(tri(:,1)); + x2 = x(tri(:,2)); + x3 = x(tri(:,3)); + + y1 = y(tri(:,1)); + y2 = y(tri(:,2)); + y3 = y(tri(:,3)); + + z1 = z(tri(:,1)); + z2 = z(tri(:,2)); + z3 = z(tri(:,3)); + + ## Calculate norm vector. + N = cross ([x2-x1, y2-y1, z2-z1], [x3-x1, y3-y1, z3-z1]); + ## Normalize. + N = diag (norm (N, "rows")) \ N; + + ## Calculate D of plane equation: Ax+By+Cz+D = 0 + D = -(N(:,1) .* x1 + N(:,2) .* y1 + N(:,3) .* z1); + + ## Calculate zi by solving plane equation for xi, yi. + zi(valid) = -(N(:,1).*xi(:)(valid) + N(:,2).*yi(:)(valid) + D) ./ N(:,3); + + else + error ('griddata: unknown interpolation METHOD: "%s"', method); + endif + + if (nargout > 1) + rx = xi; + ry = yi; + rz = zi; + else + rx = zi; + endif + endif endfunction @@ -155,10 +190,10 @@ %! x = 2*rand (100,1) - 1; %! y = 2*rand (size (x)) - 1; %! z = sin (2*(x.^2 + y.^2)); -%! [xx,yy] = meshgrid (linspace (-1,1,32)); +%! [xx,yy] = meshgrid (linspace (-1, 1, 32)); %! zz = griddata (x,y,z,xx,yy); %! mesh (xx, yy, zz); -%! title ("nonuniform grid sampled at 100 points"); +%! title ("non-uniform grid sampled at 100 points"); %!demo %! clf; @@ -166,10 +201,11 @@ %! x = 2*rand (1000,1) - 1; %! y = 2*rand (size (x)) - 1; %! z = sin (2*(x.^2 + y.^2)); -%! [xx,yy] = meshgrid (linspace (-1,1,32)); +%! [xx,yy] = meshgrid (linspace (-1, 1, 32)); %! zz = griddata (x,y,z,xx,yy); %! mesh (xx, yy, zz); -%! title ("nonuniform grid sampled at 1000 points"); +%! title ({"non-uniform grid sampled at 1,000 points", +%! 'method = "linear"'}); %!demo %! clf; @@ -177,19 +213,20 @@ %! x = 2*rand (1000,1) - 1; %! y = 2*rand (size (x)) - 1; %! z = sin (2*(x.^2 + y.^2)); -%! [xx,yy] = meshgrid (linspace (-1,1,32)); +%! [xx,yy] = meshgrid (linspace (-1, 1, 32)); %! zz = griddata (x,y,z,xx,yy,"nearest"); %! mesh (xx, yy, zz); -%! title ("nonuniform grid sampled at 1000 points with nearest neighbor"); +%! title ({"non-uniform grid sampled at 1,000 points", +%! '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(:); %! y = y + 10*(2*round (rand (size (y))) - 1) * eps; %! z = sin (2*(x.^2 + y.^2)); -%! zz = griddata (x,y,z,xx,yy,"linear"); +%! zz = griddata (x,y,z,xx,yy, "linear"); %! zz2 = sin (2*(xx.^2 + yy.^2)); %! zz2(isnan (zz)) = NaN; %! assert (zz, zz2, 100*eps); @@ -200,11 +237,16 @@ %!error griddata (1,2) %!error griddata (1,2,3) %!error griddata (1,2,3,4) -%!error griddata (1,2,3,4,5,6,7) +%!error [xi,yi] = griddata (1,2,3,4,5,6,7) +%!error griddata (1,2,3,4,5, {"linear"}) +%!error griddata (1:4, 1:3, 1:3, 1:3, 1:3) +%!error griddata (1:3, 1:4, 1:3, 1:3, 1:3) %!error griddata (1:3, 1:3, 1:4, 1:3, 1:3) -%!error griddata (1:3, 1:4, 1:3, 1:3, 1:3) -%!error griddata (1:4, 1:3, 1:3, 1:3, 1:3) %!error griddata (1:4, 1:3, ones (4,4), 1:3, 1:3) %!error griddata (1:4, 1:3, ones (3,5), 1:3, 1:3) -%!error griddata (1:3, 1:3, 1:3, 1:4, 1:3) -%!error griddata (1:3, 1:3, 1:3, 1:3, 1:4) +%!error griddata (1:3, 1:3, 1:3, 1:4, 1:3) +%!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")