--- 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
+