--- a/NEWS	Mon Oct 24 16:16:50 2011 -0400
+++ b/NEWS	Tue Oct 25 10:17:23 2011 -0700
@@ -46,6 +46,16 @@
               in pattern and in string.  Note that Matlab documents this 
               behavior but the implementation does *not* always follow it.
 
+ ** Geometry functions derived from Qhull (convhull, delaunay, voronoi)
+    have been revamped.  The options passed to the underlying qhull command
+    have been changed for better results or for Matlab compatibility.
+
+    convhull : Default options are "Qt" for 2D, 3D, 4D inputs
+               Default options are "Qt Qx" for 5D and higher
+    delaunay : Default options are "Qt Qbb Qc Qz" for 2D and 3D inputs
+               Default options are "Qt Qbb Qc Qx" for 4D and higher
+    voronoi  : No default arguments
+
  ** Matlab-compatible preference functions:
 
       addpref  getpref  ispref  rmpref  setpref

--- a/scripts/geometry/convhull.m	Mon Oct 24 16:16:50 2011 -0400
+++ b/scripts/geometry/convhull.m	Tue Oct 25 10:17:23 2011 -0700
@@ -18,35 +18,44 @@
 
 ## -*- texinfo -*-
 ## @deftypefn  {Function File} {@var{H} =} convhull (@var{x}, @var{y})
-## @deftypefnx {Function File} {@var{H} =} convhull (@var{x}, @var{y}, @var{opt})
-## Return the index vector to the points of the enclosing convex hull.  The
-## data points are defined by the x and y vectors.
+## @deftypefnx {Function File} {@var{H} =} convhull (@var{x}, @var{y}, @var{options})
+## Compute the convex hull of the set of points defined by the
+## vectors @var{x} and @var{y}.  The hull @var{H} is an index vector into
+## the set of points and specifies which points form the enclosing hull.
 ##
-## A third optional argument, which must be a string, contains extra options
-## passed to the underlying qhull command.  See the documentation for the
-## Qhull library for details.
+## An optional third argument, which must be a string or cell array of strings,
+## contains options passed to the underlying qhull command.
+## See the documentation for the Qhull library for details
+## @url{http://www.qhull.org/html/qh-quick.htm#options}.
+## The default option is @code{@{"Qt"@}}.
 ##
-## @seealso{delaunay, convhulln}
+## If @var{options} is not present or @code{[]} then the default arguments are
+## used.  Otherwise, @var{options} replaces the default argument list. 
+## 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.
+##
+## @seealso{convhulln, delaunay, voronoi}
 ## @end deftypefn
 
 ## Author: Kai Habel <kai.habel@gmx.de>
 
-function H = convhull (x, y, opt)
+function H = convhull (x, y, options)
 
   if (nargin != 2 && nargin != 3)
     print_usage ();
   endif
 
-  if (isvector (x) && isvector (y) && length (x) == length (y))
-    if (nargin == 2)
-      i = convhulln ([x(:), y(:)]);
-    elseif (ischar (opt) || iscell (opt))
-      i = convhulln ([x(:), y(:)], opt);
-    else
-      error ("convhull: third argument must be a string or cell array of strings");
-    endif
+  if (! (isvector (x) && isvector (y) && length (x) == length (y))
+      && ! size_equal (x, y))
+    error ("convhull: X and Y must be the same size");
+  elseif (nargin == 3 && ! (ischar (options) || iscellstr (options)))
+    error ("convhull: OPTIONS must be a string or cell array of strings");
+  endif
+
+  if (nargin == 2)
+    i = convhulln ([x(:), y(:)]);
   else
-    error ("convhull: first two input arguments must be vectors of same size");
+    i = convhulln ([x(:), y(:)], options);
   endif
 
   n = rows (i);
@@ -71,16 +80,21 @@
   endfor
 
   H(n + 1) = H(1);
+
 endfunction
 
-%!testif HAVE_QHULL
-%! x = -3:0.5:3;
-%! y = abs (sin (x));
-%! assert (convhull (x, y, {"s","Qci","Tcv","Pp"}), [1;7;13;12;11;10;4;3;2;1])
 
 %!demo
 %! x = -3:0.05:3;
 %! y = abs (sin (x));
 %! k = convhull (x, y);
-%! plot (x(k),y(k),'r-',x,y,'b+');
+%! plot (x(k),y(k),"r-;convex hull;", x,y,"b+;points;");
 %! axis ([-3.05, 3.05, -0.05, 1.05]);
+
+%!testif HAVE_QHULL
+%! x = -3:0.5:3;
+%! y = abs (sin (x));
+%! assert (convhull (x, y), [1;7;13;12;11;10;4;3;2;1])
+
+%% FIXME: Need input validation tests
+

--- a/scripts/geometry/delaunay.m	Mon Oct 24 16:16:50 2011 -0400
+++ b/scripts/geometry/delaunay.m	Tue Oct 25 10:17:23 2011 -0700
@@ -17,74 +17,103 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} {@var{tri} =} delaunay (@var{x}, @var{y})
-## @deftypefnx {Function File} {@var{tri} =} delaunay (@var{x}, @var{y}, @var{opt})
-## The return matrix of size [n, 3] contains a set triangles which are
-## described by the indices to the data point x and y vector.
-## The triangulation satisfies the Delaunay circum-circle criterion.
-## No other data point is in the circum-circle of the defining triangle.
+## @deftypefn  {Function File} {} delaunay (@var{x}, @var{y})
+## @deftypefnx {Function File} {@var{tri} =} delaunay (@var{x}, @var{y})
+## @deftypefnx {Function File} {@var{tri} =} delaunay (@var{x}, @var{y}, @var{options})
+## 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 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.
 ##
-## A third optional argument, which must be a string, contains extra options
-## passed to the underlying qhull command.  See the documentation for the
-## Qhull library for details.
+## An optional third argument, which must be a string or cell array of strings,
+## contains options passed to the underlying qhull command.
+## See the documentation for the Qhull library for details
+## @url{http://www.qhull.org/html/qh-quick.htm#options}.
+## The default options are @code{@{"Qt", "Qbb", "Qc", "Qz"@}}.
+##
+## If @var{options} is not present or @code{[]} then the default arguments are
+## used.  Otherwise, @var{options} replaces the default argument list. 
+## 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 (size (x));
+## y = rand (1, 10);
 ## T = delaunay (x, y);
-## X = [x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1))];
-## Y = [y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1))];
+## 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)) ];
 ## axis ([0,1,0,1]);
-## plot (X, Y, "b", x, y, "r*");
+## plot (VX, VY, "b", x, y, "r*");
 ## @end group
 ## @end example
-## @seealso{voronoi, delaunay3, delaunayn}
+## @seealso{delaunay3, delaunayn, convhull, voronoi}
 ## @end deftypefn
 
 ## Author: Kai Habel <kai.habel@gmx.de>
 
-function ret = delaunay (x, y, opt)
+function tri = delaunay (x, y, options)
 
   if (nargin != 2 && nargin != 3)
     print_usage ();
   endif
 
-  if ((isvector (x) && isvector (y) && length (x) == length (y))
-      || size_equal (x, y))
-    if (nargin == 2)
-      tri = delaunayn ([x(:), y(:)]);
-    elseif (ischar (opt) || iscellstr (opt))
-      tri = delaunayn ([x(:), y(:)], opt);
-    else
-      error ("delaunay: third argument must be a string");
-    endif
+  if (! (isvector (x) && isvector (y) && length (x) == length (y))
+      && ! size_equal (x, y))
+    error ("delaunay: X and Y must be the same size");
+  elseif (nargin == 3 && ! (ischar (options) || iscellstr (options)))
+    error ("delaunay: OPTIONS must be a string or cell array of strings");
+  endif
+
+  if (nargin == 2)
+    T = delaunayn ([x(:), y(:)]);
   else
-    error ("delaunay: first two input arguments must be matrices of same size");
+    T = delaunayn ([x(:), y(:)], options);
   endif
 
   if (nargout == 0)
     x = x(:).';
     y = y(:).';
-    X = [x(tri(:,1)); x(tri(:,2)); x(tri(:,3)); x(tri(:,1))];
-    Y = [y(tri(:,1)); y(tri(:,2)); y(tri(:,3)); y(tri(:,1))];
-    plot(X, Y, 'b', x, y, 'r*');
+    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*");
   else
-    ret = tri;
+    tri = T;
   endif
+
 endfunction
 
+
+%!demo
+%! old_state = rand ("state");
+%! restore_state = onCleanup (@() rand ("state", old_state));
+%! 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)) ];
+%! axis ([0,1,0,1]);
+%! plot (VX,VY,"b", x,y,"r*");
+
+%!testif HAVE_QHULL
+%! x = [-1, 0, 1, 0];
+%! y = [0, 1, 0, -1];
+%! assert (sortrows (sort (delaunay (x, y), 2)), [1,2,4;2,3,4]);
+
 %!testif HAVE_QHULL
 %! x = [-1, 0, 1, 0, 0];
 %! y = [0, 1, 0, -1, 0];
-%! assert (sortrows (sort (delaunay (x, y), 2)), [1,2,5;1,4,5;2,3,5;3,4,5])
+%! assert (sortrows (sort (delaunay (x, y), 2)), [1,2,5;1,4,5;2,3,5;3,4,5]);
 
-%!demo
-%! rand ('state', 1);
-%! x = rand(1,10);
-%! y = rand(size(x));
-%! T = delaunay(x,y);
-%! X = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ];
-%! Y = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ];
-%! axis([0,1,0,1]);
-%! plot(X,Y,'b',x,y,'r*');
+%% FIXME: Need input validation tests
+

--- a/scripts/geometry/delaunay3.m	Mon Oct 24 16:16:50 2011 -0400
+++ b/scripts/geometry/delaunay3.m	Tue Oct 25 10:17:23 2011 -0700
@@ -17,42 +17,61 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} {@var{T} =} delaunay3 (@var{x}, @var{y}, @var{z})
-## @deftypefnx {Function File} {@var{T} =} delaunay3 (@var{x}, @var{y}, @var{z}, @var{opt})
-## A matrix of size [n, 4] is returned.  Each row contains a
-## set of tetrahedron which are
-## described by the indices to the data point vectors (x,y,z).
+## @deftypefn  {Function File} {@var{tetr} =} delaunay3 (@var{x}, @var{y}, @var{z})
+## @deftypefnx {Function File} {@var{tetr} =} delaunay3 (@var{x}, @var{y}, @var{z}, @var{options})
+## Compute the Delaunay triangulation for a 3-D set of points.
+## 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.
 ##
-## A fourth optional argument, which must be a string or cell array of strings,
-## contains extra options passed to the underlying qhull command.  See the
-## documentation for the Qhull library for details.
-## @seealso{delaunay, delaunayn}
+## The set of tetrahedrons @var{tetr} 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.
+##
+## An optional fourth argument, which must be a string or cell array of strings,
+## contains options passed to the underlying qhull command.
+## See the documentation for the Qhull library for details
+## @url{http://www.qhull.org/html/qh-quick.htm#options}.
+## The default options are @code{@{"Qt", "Qbb", "Qc", "Qz"@}}.
+##
+## If @var{options} is not present or @code{[]} then the default arguments are
+## used.  Otherwise, @var{options} replaces the default argument list. 
+## 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.
+##
+## @seealso{delaunay, delaunayn, convhull, voronoi}
 ## @end deftypefn
 
 ## Author: Kai Habel <kai.habel@gmx.de>
 
-function tetr = delaunay3 (x, y, z, opt)
+function tetr = delaunay3 (x, y, z, options)
 
-  if (nargin != 3 && nargin != 4)
+  if (nargin < 3 || nargin > 4)
     print_usage ();
   endif
 
-  if (isvector (x) && isvector (y) &&isvector (z)
-      && length (x) == length (y) && length(x) == length (z))
-    if (nargin == 3)
-      tetr = delaunayn ([x(:), y(:), z(:)]);
-    elseif (ischar (opt) || iscell (opt))
-      tetr = delaunayn ([x(:), y(:), z(:)], opt);
-    else
-      error ("delaunay3: fourth argument must be a string or cell array of strings");
-    endif
+  if (! (isvector (x) && isvector (y) && isvector (z)
+         && length (x) == length (y) && length(x) == length (z)))
+    error ("delaunay: X, Y, and Z must be the same size");
+  elseif (nargin == 4 && ! (ischar (options) || iscellstr (options)))
+    error ("delaunay3: OPTIONS must be a string or cell array of strings");
+  endif
+
+  if (nargin == 3)
+    tetr = delaunayn ([x(:), y(:), z(:)]);
   else
-    error ("delaunay3: first three input arguments must be vectors of same size");
+    tetr = delaunayn ([x(:), y(:), z(:)], options);
   endif
 
 endfunction
 
+
 %!testif HAVE_QHULL
 %! x = [-1, -1, 1, 0, -1]; y = [-1, 1, 1, 0, -1]; z = [0, 0, 0, 1, 1];
 %! assert (sortrows (sort (delaunay3 (x, y, z), 2)), [1,2,3,4;1,2,4,5])
 
+%% FIXME: Need input validation tests
+

--- a/scripts/geometry/delaunayn.m	Mon Oct 24 16:16:50 2011 -0400
+++ b/scripts/geometry/delaunayn.m	Tue Oct 25 10:17:23 2011 -0700
@@ -17,66 +17,77 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} {@var{t} =} delaunayn (@var{p})
-## @deftypefnx {Function File} {@var{t} =} delaunayn (@var{p}, @var{opt})
-## Form the Delaunay triangulation for a set of points.
-## The Delaunay triangulation is a tessellation of the convex hull of the
-## points such that no n-sphere defined by the n-triangles contains
+## @deftypefn  {Function File} {@var{T} =} delaunayn (@var{pts})
+## @deftypefnx {Function File} {@var{T} =} delaunayn (@var{pts}, @var{options})
+## Compute the Delaunay triangulation for an N-dimensional set of points.
+## The Delaunay triangulation is a tessellation of the convex hull of a set
+## of points such that no N-sphere defined by the N-triangles contains
 ## any other points from the set.
-## The input matrix @var{p} of size @code{[n, dim]} contains @math{n}
-## points in a space of dimension dim.  The return matrix @var{t} has the
-## size @code{[m, dim+1]}.  It contains for each row a set of indices to
-## the points, which describes a simplex of dimension dim.  For example,
-## a 2-D simplex is a triangle and 3-D simplex is a tetrahedron.
+## 
+## The input matrix @var{pts} of size [n, dim] contains n points in a space of
+## dimension dim.  The return matrix @var{T} has size [m, dim+1].  Each row
+## of @var{T} contains a set of indices back into the original set of points
+## @var{pts} which describes a simplex of dimension dim.  For example, a 2-D
+## simplex is a triangle and 3-D simplex is a tetrahedron.
 ##
-## Extra options for the underlying Qhull command can be specified by the
-## second argument.  This argument is a cell array of strings.  The default
-## options depend on the dimension of the input:
+## An optional second argument, which must be a string or cell array of strings,
+## contains options passed to the underlying qhull command.
+## See the documentation for the Qhull library for details
+## @url{http://www.qhull.org/html/qh-quick.htm#options}.
+## The default options depend on the dimension of the input:
 ##
 ## @itemize
-## @item 2D and 3D: @var{opt} = @code{@{"Qt", "Qbb", "Qc"@}}
+## @item 2-D and 3-D: @var{options} = @code{@{"Qt", "Qbb", "Qc", "Qz"@}}
 ##
-## @item 4D and higher: @var{opt} = @code{@{"Qt", "Qbb", "Qc", "Qz"@}}
+## @item 4-D and higher: @var{options} = @code{@{"Qt", "Qbb", "Qc", "Qx"@}}
 ## @end itemize
 ##
-## If @var{opt} is [], then the default arguments are used.  If @var{opt}
-## is @code{@{"@w{}"@}}, then none of the default arguments are used by Qhull.
-## See the Qhull documentation for the available options.
+## If @var{options} is not present or @code{[]} then the default arguments are
+## used.  Otherwise, @var{options} replaces the default argument list. 
+## 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.
 ##
-## All options can also be specified as single string, for example
-## @code{"Qt Qbb Qc Qz"}.
-##
+## @seealso{delaunay, delaunay3, convhulln, voronoin}
 ## @end deftypefn
 
-function t = delaunayn (p, varargin)
+function T = delaunayn (pts, varargin)
+
   if (nargin < 1)
     print_usage ();
   endif
 
-  t = __delaunayn__ (p, varargin{:});
+  T = __delaunayn__ (pts, varargin{:});
 
-  if (isa (p, "single"))
+  if (isa (pts, "single"))
     myeps = eps ("single");
   else
     myeps = eps;
   endif
 
-  ## Try to remove the zero volume simplices. The volume of the i-th simplex is
-  ## given by abs(det(p(t(i,1:end-1),:)-p(t(i,2:end),:)))/prod(1:n)
-  ## (reference http://en.wikipedia.org/wiki/Simplex). Any simplex with a
-  ## relative volume less than some arbitrary criteria is rejected. The
+  ## Try to remove the zero volume simplices.  The volume of the i-th simplex is
+  ## given by abs(det(pts(T(i,1:end-1),:)-pts(T(i,2:end),:)))/prod(1:n)
+  ## (reference http://en.wikipedia.org/wiki/Simplex).  Any simplex with a
+  ## relative volume less than some arbitrary criteria is rejected.  The
   ## criteria we use is the volume of the simplex corresponding to an
   ## orthogonal simplex is equal edge length all equal to the edge length of
-  ## the original simplex. If the relative volume is 1e3*eps then the simplex
-  ## is rejected. Note division of the two volumes means that the factor
+  ## the original simplex.  If the relative volume is 1e3*eps then the simplex
+  ## is rejected.  Note division of the two volumes means that the factor
   ## prod(1:n) is dropped.
   idx = [];
-  [nt, n] = size (t);
+  [nt, n] = size (T);
+  ## FIXME: Vectorize this for loop or convert to delaunayn to .oct function
   for i = 1:nt
-    X = p(t(i,1:end-1),:) - p(t(i,2:end),:);
-    if (abs (det (X)) /  sqrt (sum (X .^ 2, 2)) < 1e3 * myeps)
-     idx = [idx, i];
+    X = pts(T(i,1:end-1),:) - pts(T(i,2:end),:);
+    if (abs (det (X)) / sqrt (sum (X .^ 2, 2)) < 1e3 * myeps)
+      idx(end+1) = i;
     endif
   endfor
-  t(idx,:) = [];
+  T(idx,:) = [];
+
 endfunction
+
+
+%% FIXME: Need tests for delaunayn
+
+%% FIXME: Need input validation tests
+

--- a/scripts/geometry/voronoi.m	Mon Oct 24 16:16:50 2011 -0400
+++ b/scripts/geometry/voronoi.m	Tue Oct 25 10:17:23 2011 -0700
@@ -18,18 +18,27 @@
 
 ## -*- texinfo -*-
 ## @deftypefn  {Function File} {} voronoi (@var{x}, @var{y})
-## @deftypefnx {Function File} {} voronoi (@var{x}, @var{y}, "plotstyle")
-## @deftypefnx {Function File} {} voronoi (@var{x}, @var{y}, "plotstyle", @var{options})
+## @deftypefnx {Function File} {} voronoi (@var{x}, @var{y}, @var{options})
+## @deftypefnx {Function File} {} voronoi (@dots{}, "linespec")
+## @deftypefnx {Function File} {} voronoi (@var{hax}, @dots{})
+## @deftypefnx {Function File} {@var{h} =} voronoi (@dots{})
 ## @deftypefnx {Function File} {[@var{vx}, @var{vy}] =} voronoi (@dots{})
-## Plot Voronoi diagram of points @code{(@var{x}, @var{y})}.
+## Plot the Voronoi diagram of points @code{(@var{x}, @var{y})}.
 ## The Voronoi facets with points at infinity are not drawn.
-## [@var{vx}, @var{vy}] = voronoi(@dots{}) returns the vertices instead of
-## plotting the
-## diagram. plot (@var{vx}, @var{vy}) shows the Voronoi diagram.
+## 
+## If "linespec" is given it is used to set the color and line style of the
+## plot.  If an axis graphics handle @var{hax} is supplied then the Voronoi
+## diagram is drawn on the specified axis rather than in a new figure.
 ##
-## A fourth optional argument, which must be a string, contains extra options
-## passed to the underlying qhull command.  See the documentation for the
-## Qhull library for details.
+## The @var{options} argument, which must be a string or cell array of strings,
+## contains options passed to the underlying qhull command.
+## See the documentation for the Qhull library for details
+## @url{http://www.qhull.org/html/qh-quick.htm#options}.
+##
+## If a single output argument is requested then the Voronoi diagram will be
+## plotted and a graphics handle to the plot is returned.
+## [@var{vx}, @var{vy}] = voronoi(@dots{}) returns the Voronoi vertices
+## instead of plotting the diagram.
 ##
 ## @example
 ## @group
@@ -57,23 +66,23 @@
 ## Added optional fourth argument to pass options to the underlying
 ## qhull command
 
-function [vvx, vvy] = voronoi (varargin)
+function [vx, vy] = voronoi (varargin)
 
   if (nargin < 1)
     print_usage ();
+  elseif (nargout > 2)
+    error ("voronoi: No more than two output arguments supported");
   endif
 
   narg = 1;
   if (isscalar (varargin{1}) && ishandle (varargin{1}))
     handl = varargin{1};
-    narg++;
     if (! strcmp (get (handl, "type"), "axes"))
       error ("voronoi: expecting first argument to be an axes object");
     endif
-  else
-    if (nargout < 2)
-      handl = gca ();
-    endif
+    narg++;
+  elseif (nargout < 2)
+    handl = gca ();
   endif
 
   if (nargin < 1 + narg || nargin > 3 + narg)
@@ -87,24 +96,22 @@
   if (narg <= nargin)
     if (iscell (varargin{narg}))
       opts = varargin(narg++);
-    elseif (ismatrix (varargin{narg}))
-      ## Accept but ignore the triangulation
+    elseif (isnumeric (varargin{narg}))
+      ## Accept, but ignore, the triangulation
       narg++;
     endif
   endif
 
   linespec = {"b"};
-  if (narg <= nargin)
-    if (ischar (varargin{narg}))
-      linespec = varargin(narg);
-    endif
+  if (narg <= nargin && ischar (varargin{narg}))
+    linespec = varargin(narg);
   endif
 
   lx = length (x);
   ly = length (y);
 
   if (lx != ly)
-    error ("voronoi: arguments must be vectors of same length");
+    error ("voronoi: X and Y must be vectors of the same length");
   endif
 
   ## Add box to approximate rays to infinity. For Voronoi diagrams the
@@ -126,12 +133,12 @@
 
   [p, c, infi] = __voronoi__ ([[x(:) ; xbox(:)], [y(:); ybox(:)]], opts{:});
 
-  idx = find (!infi);
+  idx = find (! infi);
   ll = length (idx);
   c = c(idx).';
   k = sum (cellfun ("length", c));
-  edges = cell2mat(cellfun (@(x) [x ; [x(end), x(1:end-1)]], c,
-                            "uniformoutput", false));
+  edges = cell2mat (cellfun (@(x) [x ; [x(end), x(1:end-1)]], c,
+                             "uniformoutput", false));
 
   ## Identify the unique edges of the Voronoi diagram
   edges = sortrows (sort (edges).').';
@@ -147,32 +154,34 @@
   edges (:, edgeoutside) = [];
 
   ## Get points of the diagram
-  vx = reshape (p (edges, 1), size(edges));
-  vy = reshape (p (edges, 2), size(edges));
+  Vvx = reshape (p(edges, 1), size (edges));
+  Vvy = reshape (p(edges, 2), size (edges));
 
   if (nargout < 2)
     lim = [xmin, xmax, ymin, ymax];
-    h = plot (handl, vx, vy, linespec{:}, x, y, '+');
+    h = plot (handl, Vvx, Vvy, linespec{:}, x, y, '+');
     axis (lim + 0.1 * [[-1, 1] * (lim (2) - lim (1)), ...
                        [-1, 1] * (lim (4) - lim (3))]);
     if (nargout == 1)
-      vxx = h;
+      vx = h;
     endif
-  elseif (nargout == 2)
-    vvx = vx;
-    vvy = vy;
   else
-    error ("voronoi: only two or zero output arguments supported");
+    vx = Vvx;
+    vy = Vvy;
   endif
 
 endfunction
 
-%!testif HAVE_QHULL
-%! phi=linspace(-pi,3/4*pi,8);
-%! [x,y]=pol2cart(phi,1);
-%! [vx,vy]=voronoi(x,y);
-%! assert(vx(2,:),zeros(1,size(vx,2)),eps);
-%! assert(vy(2,:),zeros(1,size(vy,2)),eps);
 
 %!demo
 %! voronoi (rand(10,1), rand(10,1));
+
+%!testif HAVE_QHULL
+%! phi = linspace (-pi, 3/4*pi, 8);
+%! [x,y] = pol2cart (phi, 1);
+%! [vx,vy] = voronoi (x,y);
+%! assert(vx(2,:), zeros (1, columns (vx)), eps);
+%! assert(vy(2,:), zeros (1, columns (vy)), eps);
+
+%% FIXME: Need input validation tests
+

--- a/scripts/geometry/voronoin.m	Mon Oct 24 16:16:50 2011 -0400
+++ b/scripts/geometry/voronoin.m	Tue Oct 25 10:17:23 2011 -0700
@@ -20,14 +20,15 @@
 ## @deftypefn  {Function File} {[@var{C}, @var{F}] =} voronoin (@var{pts})
 ## @deftypefnx {Function File} {[@var{C}, @var{F}] =} voronoin (@var{pts}, @var{options})
 ## Compute N-dimensional Voronoi facets.  The input matrix @var{pts}
-## of size [n, dim] contains n points of dimension dim.
+## of size [n, dim] contains n points in a space of dimension dim.
 ## @var{C} contains the points of the Voronoi facets.  The list @var{F}
-## contains for each facet the indices of the Voronoi points.
+## contains, for each facet, the indices of the Voronoi points.
 ##
-## A second optional argument, which must be a string, contains extra options
-## passed to the underlying qhull command.  See the documentation for the
-## Qhull library for details.
-## @seealso{voronoin, delaunay, convhull}
+## An optional second argument, which must be a string or cell array of strings,
+## contains options passed to the underlying qhull command.
+## See the documentation for the Qhull library for details
+## @url{http://www.qhull.org/html/qh-quick.htm#options}.
+## @seealso{voronoi, convhulln, delaunayn}
 ## @end deftypefn
 
 ## Author: Kai Habel <kai.habel@gmx.de>
@@ -43,17 +44,24 @@
     print_usage ();
   endif
 
-  [np, dims] = size (pts);
-  if (np > dims)
-    if (nargin == 1)
-      [C, F, infi] = __voronoi__ (pts);
-    elseif (ischar (options) || iscellstr (options))
-      [C, F, infi] = __voronoi__ (pts, options);
-    else
-      error ("voronoin: second argument must be a string or cell array of strings");
-    endif
+  [np, dim] = size (pts);
+
+  if (np <= dim)
+    error ("voronoin: number of points must be greater than their dimension");
+  elseif (nargin == 2 && ! (ischar (options) || iscellstr (options)))
+    error ("voronoin: OPTIONS argument must be a string or cell array of strings");
+  endif
 
+  if (nargin == 1)
+    [C, F] = __voronoi__ (pts);
   else
-    error ("voronoin: number of points must be greater than their dimension");
+    [C, F] = __voronoi__ (pts, options);
   endif
+
 endfunction
+
+
+%% FIXME: Need functional tests
+
+%% FIXME: Need input validation tests
+

--- a/src/DLD-FUNCTIONS/__delaunayn__.cc	Mon Oct 24 16:16:50 2011 -0400
+++ b/src/DLD-FUNCTIONS/__delaunayn__.cc	Tue Oct 25 10:17:23 2011 -0700
@@ -62,8 +62,8 @@
 
 DEFUN_DLD (__delaunayn__, args, ,
            "-*- texinfo -*-\n\
-@deftypefn  {Loadable Function} {@var{T} =} __delaunayn__ (@var{P})\n\
-@deftypefnx {Loadable Function} {@var{T} =} __delaunayn__ (@var{P}, @var{opt})\n\
+@deftypefn  {Loadable Function} {@var{T} =} __delaunayn__ (@var{pts})\n\
+@deftypefnx {Loadable Function} {@var{T} =} __delaunayn__ (@var{pts}, @var{options})\n\
 Undocumented internal function.\n\
 @end deftypefn")
 
@@ -73,7 +73,6 @@
 #ifdef HAVE_QHULL
 
   retval(0) = 0.0;
-  std::string options = "";
 
   int nargin = args.length ();
   if (nargin < 1 || nargin > 2)
@@ -86,67 +85,52 @@
   const octave_idx_type dim = p.columns ();
   const octave_idx_type n = p.rows ();
 
-  // default options
+  // Default options
+  std::string options;
   if (dim <= 3)
-    options = "Qt Qbb Qc";
+    options = "Qt Qbb Qc Qz";
   else
     options = "Qt Qbb Qc Qx";
 
-
   if (nargin == 2)
     {
-    if (args(1).is_empty ())
-      {
-        // keep default options
-      }
-    else if (args(1).is_string ())
-      {
-        // option string is directly provided
+      if (args(1).is_string ())
         options = args(1).string_value ();
-      }
-    else if (args(1).is_cell ())
-      {
-        options = "";
-
-        Cell c = args(1).cell_value ();
-        for (octave_idx_type i = 0; i < c.numel (); i++)
-          {
+      else if (args(1).is_empty ())
+        ;  // Use default options
+      else if (args(1).is_cellstr ())
+        {
+          options = "";
+          Array<std::string> tmp = args(1).cellstr_value ();
 
-            if (! c.elem(i).is_string ())
-              {
-                error ("__delaunayn__: all options must be strings");
-                return retval;
-              }
-
-            options = options + c.elem(i).string_value () + " ";
-          }
-      }
-    else
-      {
-        error ("__delaunayn__: OPT argument must be a string, cell array of strings, or empty");
-        return retval;
-      }
+          for (octave_idx_type i = 0; i < tmp.numel (); i++)
+            options += tmp(i) + " ";
+        }
+      else
+        {
+          error ("__delaunayn__: OPTIONS argument must be a string, cell array of strings, or empty");
+          return retval;
+        }
     }
 
-  //octave_stdout << "options " << options << std::endl;
-
   if (n > dim + 1)
     {
       p = p.transpose ();
       double *pt_array = p.fortran_vec ();
       boolT ismalloc = false;
 
-      OCTAVE_LOCAL_BUFFER (char, flags, 250);
+      // Qhull flags argument is not const char*
+      OCTAVE_LOCAL_BUFFER (char, flags, 9 + options.length());
 
       sprintf (flags, "qhull d %s", options.c_str ());
 
-      // If you want some debugging information replace the 0 pointer
-      // with stdout or some other file open for writing.
-
+      // Replace the 0 pointer with stdout for debugging information
       FILE *outfile = 0;
       FILE *errfile = stderr;
-
-      if (! qh_new_qhull (dim, n, pt_array, ismalloc, flags, outfile, errfile))
+      
+      int exitcode = qh_new_qhull (dim, n, pt_array, 
+                                   ismalloc, flags, outfile, errfile);
+      if (! exitcode)
         {
           // triangulate non-simplicial facets
           qh_triangulate ();
@@ -160,54 +144,48 @@
               if (! facet->upperdelaunay)
                 nf++;
 
-              // Double check
+              // Double check.  Non-simplicial facets will cause segfault below
               if (! facet->simplicial)
                 {
                   error ("__delaunayn__: Qhull returned non-simplicial facets -- try delaunayn with different options");
+                  exitcode = 1;
                   break;
                 }
             }
 
-          Matrix simpl (nf, dim+1);
-
-          FORALLfacets
+          if (! exitcode)
             {
-              if (! facet->upperdelaunay)
-                {
-                  octave_idx_type j = 0;
+              Matrix simpl (nf, dim+1);
 
-                  FOREACHvertex_ (facet->vertices)
+              FORALLfacets
+                {
+                  if (! facet->upperdelaunay)
                     {
-                      // if delaunayn crashes, enable this check
-#if 0
-                      if (j > dim)
-                        {
-                          error ("__delaunayn__: internal error. Qhull returned non-simplicial facets");
-                          return retval;
-                        }
-#endif
+                      octave_idx_type j = 0;
 
-                      simpl(i, j++) = 1 + qh_pointid(vertex->point);
+                      FOREACHvertex_ (facet->vertices)
+                        {
+                          simpl(i, j++) = 1 + qh_pointid(vertex->point);
+                        }
+                      i++;
                     }
-                  i++;
                 }
-            }
-
-          retval(0) = simpl;
 
-          // free long memory
-          qh_freeqhull (! qh_ALL);
-
-          // free short memory and memory allocator
-          int curlong, totlong;
-          qh_memfreeshort (&curlong, &totlong);
-
-          if (curlong || totlong)
-            warning ("__delaunay__: did not free %d bytes of long memory (%d pieces)",
-                     totlong, curlong);
+              retval(0) = simpl;
+            }
         }
       else
         error ("__delaunayn__: qhull failed");
+
+      // Free memory from Qhull
+      qh_freeqhull (! qh_ALL);
+
+      int curlong, totlong;
+      qh_memfreeshort (&curlong, &totlong);
+
+      if (curlong || totlong)
+        warning ("__delaunay__: did not free %d bytes of long memory (%d pieces)",
+                 totlong, curlong);
     }
   else if (n == dim + 1)
     {

--- a/src/DLD-FUNCTIONS/__voronoi__.cc	Mon Oct 24 16:16:50 2011 -0400
+++ b/src/DLD-FUNCTIONS/__voronoi__.cc	Tue Oct 25 10:17:23 2011 -0700
@@ -53,10 +53,11 @@
 #endif
 #endif
 
-DEFUN_DLD (__voronoi__, args, ,
+DEFUN_DLD (__voronoi__, args, nargout,
         "-*- texinfo -*-\n\
-@deftypefn  {Loadable Function} {@var{tri} =} __voronoi__ (@var{point})\n\
-@deftypefnx {Loadable Function} {@var{tri} =} __voronoi__ (@var{point}, @var{options})\n\
+@deftypefn  {Loadable Function} {@var{C}, @var{F} =} __voronoi__ (@var{pts})\n\
+@deftypefnx {Loadable Function} {@var{C}, @var{F} =} __voronoi__ (@var{pts}, @var{options})\n\
+@deftypefnx {Loadable Function} {@var{C}, @var{F}, @var{Inf_Pts} =} __voronoi__ (@dots{})\n\
 Undocumented internal function.\n\
 @end deftypefn")
 {
@@ -73,57 +74,50 @@
       return retval;
     }
 
-  std::string options = "qhull v Fv T0";
+  std::string options = "";
 
   if (nargin == 2)
     {
-      if (args(1).is_cellstr ())
+      if (args(1).is_string ())
+        options = args(1).string_value ();
+      else if (args(1).is_empty ())
+        ;  // Use default options
+      else if (args(1).is_cellstr ())
         {
+          options = "";
           Array<std::string> tmp = args(1).cellstr_value ();
 
           for (octave_idx_type i = 0; i < tmp.numel (); i++)
-            options += " " + tmp(i);
+            options += tmp(i) + " ";
         }
-      else if (args(1).is_string ())
-        options += " " + args(1).string_value ();
       else
         {
-          error ("__voronoi__: OPTIONS argument must be a string or cellstr");
+          error ("__voronoi__: OPTIONS argument must be a string, cell array of strings, or empty");
           return retval;
         }
     }
 
   Matrix p (args(0).matrix_value ());
-
   const octave_idx_type dim = p.columns ();
   const octave_idx_type np = p.rows ();
+
   p = p.transpose ();
-
   double *pt_array = p.fortran_vec ();
 
-  //double  pt_array[dim * np];
-  //for (int i = 0; i < np; i++)
-  //  {
-  //    for (int j = 0; j < dim; j++)
-  //      {
-  //        pt_array[j+i*dim] = p(i,j);
-  //      }
-  //  }
-
   boolT ismalloc = false;
 
-  // If you want some debugging information replace the 0 pointer
-  // with stdout or some other file open for writing.
-
+  // Replace the 0 pointer with stdout for debugging information
   FILE *outfile = 0;
   FILE *errfile = stderr;
 
-  // Qhull flags argument is not const char*...
-  OCTAVE_LOCAL_BUFFER (char, flags, options.length () + 1);
+  // Qhull flags argument is not const char*
+  OCTAVE_LOCAL_BUFFER (char, flags, 12 + options.length ());
 
-  strcpy (flags, options.c_str ());
+  sprintf (flags, "qhull v FV %s", options.c_str ());
 
-  if (! qh_new_qhull (dim, np, pt_array, ismalloc, flags, outfile, errfile))
+  int exitcode = qh_new_qhull (dim, np, pt_array, 
+                               ismalloc, flags, outfile, errfile);
+  if (! exitcode) 
     {
       facetT *facet;
       vertexT *vertex;
@@ -134,6 +128,7 @@
       for (i = 0; i < np; i++)
         ni[i] = 0;
       qh_setvoronoi_all ();
+
       bool infinity_seen = false;
       facetT *neighbor, **neighborp;
       coordT *voronoi_vertex;
@@ -147,6 +142,7 @@
         {
           if (qh hull_dim == 3)
             qh_order_vertexneighbors (vertex);
+          
           infinity_seen = false;
 
           FOREACHneighbor_ (vertex)
@@ -173,9 +169,7 @@
       for (octave_idx_type d = 0; d < dim; d++)
         v(0,d) = octave_Inf;
 
-      boolMatrix AtInf (np, 1);
-      for (i = 0; i < np; i++)
-        AtInf(i) = false;
+      boolMatrix AtInf (np, 1, false);
       octave_value_list F (np, octave_value ());
       k = 0;
       i = 0;
@@ -189,6 +183,7 @@
         {
           if (qh hull_dim == 3)
             qh_order_vertexneighbors(vertex);
+
           infinity_seen = false;
           RowVector facet_list (ni[k++]);
           m = 0;
@@ -229,24 +224,27 @@
       for (i = 0; i < r; i++)
         C.elem (i) = F(i);
 
-      retval(0) = v;
+      if (nargout == 3)
+        {
+          AtInf.resize (r, 1);
+          retval(2) = AtInf;
+        }
       retval(1) = C;
-      AtInf.resize (r, 1);
-      retval(2) = AtInf;
-
-      // free long memory
-      qh_freeqhull (! qh_ALL);
-
-      // free short memory and memory allocator
-      int curlong, totlong;
-      qh_memfreeshort (&curlong, &totlong);
-
-      if (curlong || totlong)
-        warning ("__voronoi__: did not free %d bytes of long memory (%d pieces)", totlong, curlong);
+      retval(0) = v;
     }
   else
     error ("__voronoi__: qhull failed");
 
+  // free memory from Qhull
+  qh_freeqhull (! qh_ALL);
+
+  int curlong, totlong;
+  qh_memfreeshort (&curlong, &totlong);
+
+  if (curlong || totlong)
+    warning ("__voronoi__: did not free %d bytes of long memory (%d pieces)",
+             totlong, curlong);
+
 #else
   error ("__voronoi__: not available in this version of Octave");
 #endif

--- a/src/DLD-FUNCTIONS/convhulln.cc	Mon Oct 24 16:16:50 2011 -0400
+++ b/src/DLD-FUNCTIONS/convhulln.cc	Tue Oct 25 10:17:23 2011 -0700
@@ -53,24 +53,39 @@
 
 DEFUN_DLD (convhulln, args, nargout,
   "-*- texinfo -*-\n\
-@deftypefn  {Loadable Function} {@var{h} =} convhulln (@var{p})\n\
-@deftypefnx {Loadable Function} {@var{h} =} convhulln (@var{p}, @var{opt})\n\
+@deftypefn  {Loadable Function} {@var{h} =} convhulln (@var{pts})\n\
+@deftypefnx {Loadable Function} {@var{h} =} convhulln (@var{pts}, @var{options})\n\
 @deftypefnx {Loadable Function} {[@var{h}, @var{v}] =} convhulln (@dots{})\n\
-Return an index vector to the points of the enclosing convex hull.\n\
-The input matrix of size [n, dim] contains n points of dimension dim.\n\n\
-If a second optional argument is given, it must be a string or cell array\n\
-of strings containing options for the underlying qhull command.  (See\n\
-the Qhull documentation for the available options.)  The default options\n\
-are \"s Qci Tcv\".\n\
-If the second output @var{v} is requested the volume of the convex hull is\n\
-calculated.\n\n\
-@seealso{convhull, delaunayn}\n\
+Compute the convex hull of the set of points @var{pts} which is a matrix\n\
+of size [n, dim] containing n points in a space of dimension dim.\n\
+The hull @var{h} is an index vector into the set of points and specifies\n\
+which points form the enclosing hull.\n\
+\n\
+An optional second argument, which must be a string or cell array of strings,\n\
+contains options passed to the underlying qhull command.\n\
+See the documentation for the Qhull library for details\n\
+@url{http://www.qhull.org/html/qh-quick.htm#options}.\n\
+The default options depend on the dimension of the input:\n\
+\n\
+@itemize\n\
+@item 2D, 3D, 4D: @var{options} = @code{@{\"Qt\"@}}\n\
+\n\
+@item 5D and higher: @var{options} = @code{@{\"Qt\", \"Qx\"@}}\n\
+@end itemize\n\
+\n\
+If @var{options} is not present or @code{[]} then the default arguments are\n\
+used.  Otherwise, @var{options} replaces the default argument list.\n\
+To append user options to the defaults it is necessary to repeat the\n\
+default arguments in @var{options}.  Use a null string to pass no arguments.\n\
+\n\
+If the second output @var{v} is requested the volume of the enclosing\n\
+convex hull is calculated.\n\n\
+@seealso{convhull, delaunayn, voronoin}\n\
 @end deftypefn")
 {
   octave_value_list retval;
 
 #ifdef HAVE_QHULL
-  std::string options;
 
   int nargin = args.length ();
   if (nargin < 1 || nargin > 2)
@@ -79,85 +94,88 @@
       return retval;
     }
 
+  Matrix p (args(0).matrix_value ());
+  const octave_idx_type dim = p.columns ();
+  const octave_idx_type n = p.rows ();
+
+  // Default options
+  std::string options;
+  if (dim <= 4)
+    options = "Qt";
+  else
+    options = "Qt Qx";
+
   if (nargin == 2)
     {
-      if (args (1).is_string ())
+      if (args(1).is_string ())
         options = args(1).string_value ();
-      else if (args(1).is_cell ())
+      else if (args(1).is_empty ())
+        ; // Use default options
+      else if (args(1).is_cellstr ())
         {
-          Cell c = args(1).cell_value ();
           options = "";
-          for (octave_idx_type i = 0; i < c.numel (); i++)
-            {
-              if (! c.elem(i).is_string ())
-                {
-                  error ("convhulln: OPT must be a string or cell array of strings");
-                  return retval;
-                }
+          Array<std::string> tmp = args(1).cellstr_value ();
 
-              options = options + c.elem(i).string_value() + " ";
-            }
+          for (octave_idx_type i = 0; i < tmp.numel (); i++)
+            options += tmp(i) + " ";
         }
       else
         {
-          error ("convhulln: OPT must be a string or cell array of strings");
+          error ("convhulln: OPTIONS must be a string, cell array of strings, or empty");
           return retval;
         }
-    }
-  else
-    // turn on some consistency checks
-    options = "s Qci Tcv";
+     }
 
-  Matrix p (args(0).matrix_value ());
 
-  const octave_idx_type dim = p.columns ();
-  const octave_idx_type n = p.rows ();
   p = p.transpose ();
-
   double *pt_array = p.fortran_vec ();
-
-  boolT ismalloc = False;
-
-  std::ostringstream buf;
+  boolT ismalloc = false;
 
-  buf << "qhull QJ " << options;
-
-  std::string buf_string = buf.str ();
+  // FIXME: we can't just pass options.c_str () to qh_new_qhull
+  // because the argument is not declared const.  Ugh.  Unless qh_new_qhull
+  // really needs to modify this argument, someone should fix QHULL.
+  OCTAVE_LOCAL_BUFFER (char, flags, 7 + options.length ());
 
-  // FIXME -- we can't just pass buf_string.c_str () to qh_new_qhull
-  // because the argument is not declared const.  Ugh.  Unless
-  // qh_new_qhull really needs to modify this argument, someone should
-  // fix QHULL.
+  sprintf (flags, "qhull %s", options.c_str ());
 
-  OCTAVE_LOCAL_BUFFER (char, flags, buf_string.length () + 1);
-
-  strcpy (flags, buf_string.c_str ());
-
-  if (! qh_new_qhull (dim, n, pt_array, ismalloc, flags, 0, stderr))
+  // Replace the 0 pointer with stdout for debugging information
+  FILE *outfile = 0;
+  FILE *errfile = stderr;
+      
+  int exitcode = qh_new_qhull (dim, n, pt_array, 
+                               ismalloc, flags, outfile, errfile);
+  if (! exitcode)
     {
-      // If you want some debugging information replace the NULL
-      // pointer with stdout
-
       vertexT *vertex, **vertexp;
       facetT *facet;
       setT *vertices;
-      unsigned int nf = qh num_facets;
+      bool nonsimp_seen = false;
+      octave_idx_type nf = qh num_facets;
 
-      Matrix idx (nf, dim);
+      Matrix idx (nf, dim + 1);
 
-      octave_idx_type j, i = 0;
+      octave_idx_type i = 0, j;
       FORALLfacets
         {
           j = 0;
-          if (! facet->simplicial)
-            // should never happen with QJ
-            error ("convhulln: non-simplicial facet");
+
+          if (! nonsimp_seen && ! facet->simplicial)
+            {
+              nonsimp_seen = true;
 
+              if (options.find ("QJ") != std::string::npos)
+                {
+                  // should never happen with QJ
+                  error ("convhulln: qhull failed.  Option 'QJ' returned non-simplicial facet");
+                  break;
+                }
+            }
           if (dim == 3)
             {
               vertices = qh_facet3vertex (facet);
               FOREACHvertex_ (vertices)
                 idx(i, j++) = 1 + qh_pointid(vertex->point);
+
               qh_settempfree (&vertices);
             }
           else
@@ -174,11 +192,15 @@
                 }
             }
           if (j < dim)
-            // likewise but less fatal
             warning ("facet %d only has %d vertices", i, j);
+
           i++;
         }
 
+      // Remove extra dimension if all facets were simplicial
+      if (! nonsimp_seen)
+        idx.resize (nf, dim, 0.0);
+
       if (nargout == 2)
         // calculate volume of convex hull
         // taken from qhull src/geom2.c
@@ -213,19 +235,21 @@
           retval(1) = octave_value (qh totvol);
         }
 
-      retval(0) = octave_value (idx);
+      retval(0) = idx;
     }
+  else
+    error ("convhulln: qhull failed");
 
-  // free long memory
+  // free memory from Qhull
   qh_freeqhull (! qh_ALL);
 
-  // free short memory and memory allocator
   int curlong, totlong;
   qh_memfreeshort (&curlong, &totlong);
 
   if (curlong || totlong)
     warning ("convhulln: did not free %d bytes of long memory (%d pieces)",
-            totlong, curlong);
+             totlong, curlong);
+
 #else
   error ("convhulln: not available in this version of Octave");
 #endif
@@ -236,10 +260,23 @@
 /*
 %!testif HAVE_QHULL
 %! cube = [0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1];
-%! [h, v] = convhulln(cube,'Pp');
+%! [h, v] = convhulln (cube);
+%! assert (size (h), [6 4]); 
+%! h = sortrows (sort (h, 2), [1:4]);
+%! assert (h, [1 2 3 4; 1 2 5 6; 1 4 5 8; 2 3 6 7; 3 4 7 8; 5 6 7 8]);
+%! assert (v, 1, 10*eps);
+
+%!testif HAVE_QHULL
+%! cube = [0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1];
+%! [h, v] = convhulln (cube, "QJ");
+%! assert (size (h), [12 3]); 
+%! assert (sortrows (sort (h, 2), [1:3]), [1 2 4; 1 2 5; 1 4 5; 2 3 4; 2 3 6; 2 5 6; 3 4 8; 3 6 7; 3 7 8; 4 5 8; 5 6 8; 6 7 8]);
 %! assert (v, 1.0, 1e6*eps);
+
 %!testif HAVE_QHULL
 %! tetrahedron = [1 1 1;-1 -1 1;-1 1 -1;1 -1 -1];
-%! [h, v] = convhulln(tetrahedron,'Pp');
-%! assert (v, 8/3, 1e6*eps);
+%! [h, v] = convhulln (tetrahedron);
+%! h = sortrows (sort (h, 2), [1 2 3]);
+%! assert (h, [1 2 3;1 2 4; 1 3 4; 2 3 4]);
+%! assert (v, 8/3, 10*eps);
 */