######################################################################## ## ## Copyright (C) 1999-2020 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or . ## ## This file is part of Octave. ## ## Octave is free software: you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## . ## ######################################################################## ## -*- texinfo -*- ## @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{}) ## @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}) ## ## 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. ## ## 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. ## ## 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, varargin) if (nargin < 5) print_usage (); endif if (nargin > 6) ## current 2D implementaiton has nargin max = 6, since no triangulation ## options are passed to the 2D algorithm. 3D algorithm requires nargin >=7 if (nargout > 1) error("griddata: only one output argument valid for 3D interpolation") else rx = griddata3 (x, y, z, varargin{:}); endif else ## 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})) method = varargin{3}; else error ("griddata: unknown interpolation METHOD") endif else method = "linear"; 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 any (strcmp (method, {"cubic","natural","v4"})) ## TODO: implement compatible implementations of all 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(:)); ## 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"); endif if (nargout > 1) rx = xi; ry = yi; rz = zi; else rx = zi; endif endif endfunction %!demo %! clf; %! colormap ("default"); %! 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)); %! zz = griddata (x,y,z,xx,yy); %! mesh (xx, yy, zz); %! title ("nonuniform grid sampled at 100 points"); %!demo %! clf; %! colormap ("default"); %! 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)); %! zz = griddata (x,y,z,xx,yy); %! mesh (xx, yy, zz); %! title ("nonuniform grid sampled at 1000 points"); %!demo %! clf; %! colormap ("default"); %! 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)); %! zz = griddata (x,y,z,xx,yy,"nearest"); %! mesh (xx, yy, zz); %! title ("nonuniform grid sampled at 1000 points with nearest neighbor"); %!testif HAVE_QHULL %! [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"); %! zz2 = sin (2*(xx.^2 + yy.^2)); %! zz2(isnan (zz)) = NaN; %! assert (zz, zz2, 100*eps); ## Test input validation %!error griddata () %!error griddata (1) %!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 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)