--- a/scripts/geometry/delaunay.m	Fri Sep 26 08:54:40 2014 -0700
+++ b/scripts/geometry/delaunay.m	Fri Sep 26 09:02:53 2014 -0700
@@ -17,22 +17,32 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} {} delaunay (@var{x}, @var{y})
-## @deftypefnx {Function File} {} delaunay (@var{x})
-## @deftypefnx {Function File} {} delaunay (@dots{}, @var{options})
-## @deftypefnx {Function File} {@var{tri} =} delaunay (@dots{})
-## Compute the Delaunay triangulation for a 2-D set of points.
-## The return value @var{tri} is a set of triangles which satisfies the
-## Delaunay circum-circle criterion, i.e., only a single data point from
-## [@var{x}, @var{y}] is within the circum-circle of the defining triangle.
-## The input @var{x} may also be a matrix with two columns where the first
-## column contains x-data and the second y-data.
+## @deftypefn  {Function File} {@var{tri} =} delaunay (@var{x}, @var{y})
+## @deftypefnx {Function File} {@var{tetr} =} delaunay (@var{x}, @var{y}, @var{z})
+## @deftypefnx {Function File} {@var{tri} =} delaunay (@var{x})
+## @deftypefnx {Function File} {@var{tri} =} delaunay (@dots{}, @var{options})
+## Compute the Delaunay triangulation for a 2-D or 3-D set of points.
 ##
-## The set of triangles @var{tri} is a matrix of size [n, 3].  Each
-## row defines a triangle and the three columns are the three vertices
-## of the triangle.  The value of @code{@var{tri}(i,j)} is an index into
-## @var{x} and @var{y} for the location of the j-th vertex of the i-th
-## triangle.
+## For 2-D sets, the return value @var{tri} is a set of triangles which
+## satisfies the Delaunay circum-circle criterion, i.e., only a single data
+## point from [@var{x}, @var{y}] is within the circum-circle of the defining
+## triangle.  The set of triangles @var{tri} is a matrix of size [n, 3].  Each
+## row defines a triangle and the three columns are the three vertices of the
+## triangle.  The value of @code{@var{tri}(i,j)} is an index into @var{x} and
+## @var{y} for the location of the j-th vertex of the i-th triangle.
+##
+## For 3-D sets, the return value @var{tetr} is a set of tetrahedrons which
+## satisfies the Delaunay circum-circle criterion, i.e., only a single data
+## point from [@var{x}, @var{y}, @var{z}] is within the circum-circle of the
+## defining tetrahedron.  The set of tetrahedrons is a matrix of size [n, 4].
+## Each row defines a tetrahedron and the four columns are the four vertices of
+## the tetrahedron.  The value of @code{@var{tetr}(i,j)} is an index into
+## @var{x}, @var{y}, @var{z} for the location of the j-th vertex of the i-th
+## tetrahedron.
+##
+## The input @var{x} may also be a matrix with two or three columns where the
+## first column contains x-data, the second y-data, and the optional third
+## column contains z-data.
 ##
 ## The optional last argument, which must be a string or cell array of strings,
 ## contains options passed to the underlying qhull command.
@@ -45,63 +55,84 @@
 ## To append user options to the defaults it is necessary to repeat the 
 ## default arguments in @var{options}.  Use a null string to pass no arguments.
 ##
-## If no output argument is specified the resulting Delaunay triangulation 
-## is plotted along with the original set of points.
-##
 ## @example
 ## @group
 ## x = rand (1, 10);
 ## y = rand (1, 10);
-## T = delaunay (x, y);
-## VX = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ];
-## VY = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ];
+## tri = delaunay (x, y);
+## triplot (tri, x, y);
+## hold on;
+## plot (x, y, "r*");
 ## axis ([0,1,0,1]);
-## plot (VX, VY, "b", x, y, "r*");
 ## @end group
 ## @end example
-## @seealso{delaunay3, delaunayn, convhull, voronoi, triplot, trimesh, trisurf}
+## @seealso{delaunayn, convhull, voronoi, triplot, trimesh, tetramesh, trisurf}
 ## @end deftypefn
 
 ## Author: Kai Habel <kai.habel@gmx.de>
 
 function tri = delaunay (varargin)
 
-  if (nargin < 1 || nargin > 3)
+  if (nargin < 1 || nargin > 4)
     print_usage ();
   endif
 
+  z = [];
   options = [];
 
   switch (nargin)
 
     case 1
-      if (! ismatrix (varargin{1}) || columns (varargin{1}) != 2)
-          error ("delaunay: X must be a matrix with 2 columns");
+      if (! ismatrix (varargin{1})
+          || (columns (varargin{1}) != 2 && columns (varargin{1}) != 3))
+          error ("delaunay: X must be a matrix with 2 or 3 columns");
       else
         x = varargin{1}(:,1);
         y = varargin{1}(:,2);
+        if (columns (varargin{1}) == 3)
+          z = varargin{1}(:,3);
+        endif
       endif
     
     case 2
       if (isnumeric (varargin{2}))
         x = varargin{1};
         y = varargin{2};
-      elseif (ischar (varargin{2}) || iscellstr (varargin{2}))
+      elseif (! (ischar (varargin{2}) || iscellstr (varargin{2})))
+        error ("delaunay: OPTIONS must be a string or cell array of strings");
+      else
         options = varargin{2};
-        if (! ismatrix (varargin{1}) && columns (varargin{1}) != 2)
-            error ("delaunay: X must be a matrix with 2 columns");
+        ncols = columns (varargin{1});
+
+        if (! ismatrix (varargin{1}) || (ncols != 2 && ncols != 3))
+          error ("delaunay: X must be a matrix with 2 or 3 columns");
         else
           x = varargin{1}(:,1);
           y = varargin{1}(:,2);
+          if (ncols == 3)
+            z = varargin{1}(:,3);
+          endif
         endif
-      else
-        error ("delaunay: OPTIONS must be a string or cell array of strings");
       endif
 
     case 3
+      if (isnumeric (varargin{3}))
+        x = varargin{1};
+        y = varargin{2};
+        z = varargin{3};
+      elseif (! (ischar (varargin{3}) || iscellstr (varargin{3})))
+        error ("delaunay: OPTIONS must be a string or cell array of strings");
+      else
+        x = varargin{1};
+        y = varargin{2};
+        options = varargin{3};
+      endif
+
+    case 4
       x = varargin{1};
       y = varargin{2};
-      options = varargin{3};
+      z = varargin{3};
+      options = varargin{4};
 
       if (! (ischar (options) || iscellstr (options)))
         error ("delaunay: OPTIONS must be a string or cell array of strings");
@@ -109,20 +140,16 @@
 
   endswitch
 
-  if (! (isequal(size(x),size(y))))
-    error ("delaunay: X and Y must be the same size");
-  endif
-
-  T = delaunayn ([x(:), y(:)], options);
-
-  if (nargout == 0)
-    x = x(:).';
-    y = y(:).';
-    VX = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ];
-    VY = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ];
-    plot (VX, VY, "b", x, y, "r*");
+  if (isempty (z))
+    if (! size_equal (x, y))
+      error ("delaunay: X and Y must be the same size");
+    endif
+    tri = delaunayn ([x(:), y(:)], options);
   else
-    tri = T;
+    if (! size_equal (x, y, z))
+      error ("delaunay: X, Y, and Z must be the same size");
+    endif
+    tri = delaunayn ([x(:), y(:), z(:)], options);
   endif
 
 endfunction
@@ -134,11 +161,11 @@
 %! rand ("state", 1);
 %! x = rand (1,10);
 %! y = rand (1,10);
-%! T = delaunay (x,y);
-%! VX = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ];
-%! VY = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ];
+%! tri = delaunay (x,y);
 %! clf;
-%! plot (VX,VY,"b", x,y,"r*");
+%! triplot (tri, x, y);
+%! hold on;
+%! plot (x, y, "r*");
 %! axis ([0,1,0,1]);
 
 %!testif HAVE_QHULL
@@ -165,5 +192,19 @@
 %! y = [5 7 8; 1 2 3];
 %! assert (sortrows (sort (delaunay (x, y), 2)), [1,2,4;1,3,4;1,3,5;3,4,6]);
 
-%% FIXME: Need input validation tests
+## Test 3-D input
+%!testif HAVE_QHULL
+%! x = [-1, -1, 1, 0, -1]; y = [-1, 1, 1, 0, -1]; z = [0, 0, 0, 1, 1];
+%! assert (sortrows (sort (delaunay (x, y, z), 2)), [1,2,3,4;1,2,4,5])
 
+%% Input validation tests
+%!error delaunay ()
+%!error delaunay (1,2,3,4,5)
+%!error <X must be a matrix with 2 or 3 columns> delaunay (ones (2,4))
+%!error <OPTIONS must be a string or cell array> delaunay (ones (2,2), struct())
+%!error <X must be a matrix with 2 or 3 columns> delaunay (ones (2,4), "")
+%!error <OPTIONS must be a string or cell array> delaunay (ones (2,2), ones (2,2), struct())
+%!error <OPTIONS must be a string or cell array> delaunay (ones (2,2), ones (2,2), ones (2,2), struct())
+%!error <X and Y must be the same size> delaunay (1, [1 2])
+%!error <X, Y, and Z must be the same size> delaunay (1, [1 2], [1 2])
+