--- a/src/DLD-FUNCTIONS/convhulln.cc
+++ b/src/DLD-FUNCTIONS/convhulln.cc
@@ -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);
 */