# HG changeset patch
# User Nicholas R. Jankowski <jankowskin@asme.org>
# Date 1575007242 18000
#      Fri Nov 29 01:00:42 2019 -0500
# Node ID f3a534a92690a764996bc9e9db53146e49b2f47d
# Parent  dcf6569cd9738f26288565a35504f52680c69f66
griddata 3D compatibility fix, and griddata3 method input fix (bug #57323)
* /scripts/geometry/griddata.m: enabled griddata to accept griddata3 input for Matlab compatibility
* /scripts/geometry/griddata3.m: fixed problem where it failed to pass method to griddatan

diff -r dcf6569cd973 -r f3a534a92690 scripts/geometry/griddata.m
--- a/scripts/geometry/griddata.m	Thu Nov 28 22:25:07 2019 +0100
+++ b/scripts/geometry/griddata.m	Fri Nov 29 01:00:42 2019 -0500
@@ -20,18 +20,32 @@
 ## @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} =} griddata3 (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi})
+## @deftypefnx {} {@var{vi} =} griddata3 (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi}, @var{method})
+## @deftypefnx {} {@var{vi} =} griddata3 (@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
 
@@ -39,106 +53,135 @@
 ## Adapted-by:  Alexander Barth <barth.alexander@gmail.com>
 ##              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
+    
+    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
+    ## 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 (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
-  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
+    ## 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
+    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(:);
+    x = x(:);
+    y = y(:);
+    z = z(:);
 
-  ## Triangulate data.
-  tri = delaunay (x, y);
-  zi = NaN (size (xi));
+    ## Triangulate data.
+    tri = delaunay (x, y);
+    zi = NaN (size (xi));
 
-  if (strcmp (method, "cubic"))
-    error ("griddata: cubic interpolation not yet implemented");
+    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, "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(:));
+    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);
+      ## Only keep the points within triangles.
+      valid = ! isnan (tri_list);
+      tri_list = tri_list(valid);
+      nr_t = rows (tri_list);
 
-    tri = tri(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));
+      ## 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));
+      y1 = y(tri(:,1));
+      y2 = y(tri(:,2));
+      y3 = y(tri(:,3));
 
-    z1 = z(tri(:,1));
-    z2 = z(tri(:,2));
-    z3 = z(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 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 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);
+      ## 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
+    else
+      error ("griddata: unknown interpolation METHOD");
+    endif
 
-  if (nargout > 1)
-    rx = xi;
-    ry = yi;
-    rz = zi;
-  else
-    rx = zi;
+    if (nargout > 1)
+      rx = xi;
+      ry = yi;
+      rz = zi;
+    else
+      rx = zi;
+    endif
+
   endif
 
 endfunction
@@ -189,6 +232,33 @@
 %! zz2(isnan (zz)) = NaN;
 %! assert (zz, zz2, 100*eps);
 
+## Verify gridtata3 tests pass when called through griddata
+%!testif HAVE_QHULL
+%! old_state = rand ("state");
+%! restore_state = onCleanup (@() rand ("state", old_state));
+%! rand ("state", 0);
+%! x = 2 * rand (1000, 1) - 1;
+%! y = 2 * rand (1000, 1) - 1;
+%! z = 2 * rand (1000, 1) - 1;
+%! v = x.^2 + y.^2 + z.^2;
+%! [xi, yi, zi] = meshgrid (-0.8:0.2:0.8);
+%! vi = griddata (x, y, z, v, xi, yi, zi, "linear");
+%! vv = vi - xi.^2 - yi.^2 - zi.^2;
+%! assert (max (abs (vv(:))), 0, 0.1);
+
+%!testif HAVE_QHULL
+%! old_state = rand ("state");
+%! restore_state = onCleanup (@() rand ("state", old_state));
+%! rand ("state", 0);
+%! x = 2 * rand (1000, 1) - 1;
+%! y = 2 * rand (1000, 1) - 1;
+%! z = 2 * rand (1000, 1) - 1;
+%! v = x.^2 + y.^2 + z.^2;
+%! [xi, yi, zi] = meshgrid (-0.8:0.2:0.8);
+%! vi = griddata (x, y, z, v, xi, yi, zi, "nearest");
+%! vv = vi - xi.^2 - yi.^2 - zi.^2;
+%! assert (max (abs (vv(:))), 0.385, 0.1);
+
 ## Test input validation
 %!error griddata ()
 %!error griddata (1)
@@ -203,3 +273,4 @@
 %!error <the columns and rows of Z> griddata (1:4, 1:3, ones (3,5), 1:3, 1:3)
 %!error <matrices of same size> griddata (1:3, 1:3, 1:3, 1:4, 1:3)
 %!error <matrices of same size> griddata (1:3, 1:3, 1:3, 1:3, 1:4)
+
diff -r dcf6569cd973 -r f3a534a92690 scripts/geometry/griddata3.m
--- a/scripts/geometry/griddata3.m	Thu Nov 28 22:25:07 2019 +0100
+++ b/scripts/geometry/griddata3.m	Fri Nov 29 01:00:42 2019 -0500
@@ -38,7 +38,7 @@
 
 ## Author: David Bateman <dbateman@free.fr>
 
-function vi = griddata3 (x, y, z, v, xi, yi, zi, method, varargin)
+function vi = griddata3 (x, y, z, v, xi, yi, zi, method = 'linear', varargin)
 
   if (nargin < 7)
     print_usage ();
@@ -69,7 +69,7 @@
     error ("griddata3: XI, YI, and ZI must be vectors or matrices of the same size");
   endif
 
-  vi = griddatan ([x(:), y(:), z(:)], v(:), [xi(:), yi(:), zi(:)], varargin{:});
+  vi = griddatan ([x(:), y(:), z(:)], v(:), [xi(:), yi(:), zi(:)], method, varargin{:});
   vi = reshape (vi, size (xi));
 
 endfunction
@@ -99,4 +99,4 @@
 %! [xi, yi, zi] = meshgrid (-0.8:0.2:0.8);
 %! vi = griddata3 (x, y, z, v, xi, yi, zi, "nearest");
 %! vv = vi - xi.^2 - yi.^2 - zi.^2;
-%! assert (max (abs (vv(:))), 0, 0.1);
+%! assert (max (abs (vv(:))), 0.385, 0.1);
