# HG changeset patch
# User Nicholas R. Jankowski <jankowskin@asme.org>
# Date 1581907033 18000
#      Sun Feb 16 21:37:13 2020 -0500
# Node ID f271c34441c09ad137aef044886c792dd84891f4
# Parent  a81fda0cb5f7800f80610a1599304467b8a2f2ef
griddata compatibility fixed to allow 3D inputs (bug #57323)

* /scripts/general/griddata.m: modified to accept 3D input as passthrough for griddata3
* /NEWS: added griddata 3D change to Matlab compatibility section

diff -r a81fda0cb5f7 -r f271c34441c0 NEWS
--- a/NEWS	Sun Feb 16 21:08:35 2020 -0500
+++ b/NEWS	Sun Feb 16 21:37:13 2020 -0500
@@ -178,6 +178,9 @@
 
 - Several validator functions of type `mustBe*` have been added.  See
   the list of new functions below.

+- The function `griddata` now accepts 2D and 3D inputs by passing data directly
+  to `griddata3`.  
 
 ### Alphabetical list of new functions added in Octave 6
 
diff -r a81fda0cb5f7 -r f271c34441c0 scripts/geometry/griddata.m
--- a/scripts/geometry/griddata.m	Sun Feb 16 21:08:35 2020 -0500
+++ b/scripts/geometry/griddata.m	Sun Feb 16 21:37:13 2020 -0500
@@ -27,123 +27,164 @@
 ## @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 interpolation.
+## Generate a regular mesh from irregular data using 2-D or 3-D interpolation.
 ##
-## The function is defined by @code{@var{z} = f (@var{x}, @var{y})}.  Inputs
+## 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.
 ##
-## 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
+  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
 
-  ## 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");
+    if (nargout > 1)
+      error("griddata: only one output argument valid for 3D interpolation")
+    else
+      rx = griddata3 (x, y, z, varargin{:});
     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");
-
-  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
+    ## for nargin 5 or 6, assign varargin terms to variables for 2D algorithm
+    xi = varargin{1};
+    yi = varargin{2};
 
-  if (nargout > 1)
-    rx = xi;
-    ry = yi;
-    rz = zi;
-  else
-    rx = zi;
+    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
