--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/dldfcn/__voronoi__.cc	Sat Aug 18 16:23:39 2012 -0400
@@ -0,0 +1,334 @@
+/*
+
+Copyright (C) 2000-2012 Kai Habel
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+/*
+20. Augiust 2000 - Kai Habel: first release
+*/
+
+/*
+2003-12-14 Rafael Laboissiere <rafael@laboissiere.net>
+Added optional second argument to pass options to the underlying
+qhull command
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <cstdio>
+
+#include <list>
+
+#include "lo-ieee.h"
+
+#include "Cell.h"
+#include "defun-dld.h"
+#include "error.h"
+#include "oct-obj.h"
+#include "unwind-prot.h"
+
+#if defined (HAVE_QHULL)
+# include "oct-qhull.h"
+# if defined (NEED_QHULL_VERSION)
+char qh_version[] = "__voronoi__.oct 2007-07-24";
+# endif
+#endif
+
+static void
+close_fcn (FILE *f)
+{
+  gnulib::fclose (f);
+}
+
+DEFUN_DLD (__voronoi__, args, ,
+        "-*- texinfo -*-\n\
+@deftypefn  {Loadable Function} {@var{C}, @var{F} =} __voronoi__ (@var{caller}, @var{pts})\n\
+@deftypefnx {Loadable Function} {@var{C}, @var{F} =} __voronoi__ (@var{caller}, @var{pts}, @var{options})\n\
+@deftypefnx {Loadable Function} {@var{C}, @var{F}, @var{Inf_Pts} =} __voronoi__ (@dots{})\n\
+Undocumented internal function.\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+
+  std::string caller = args(0).string_value ();
+
+#if defined (HAVE_QHULL)
+
+  retval(0) = 0.0;
+
+  int nargin = args.length ();
+  if (nargin < 2 || nargin > 3)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  Matrix points = args(1).matrix_value ();
+  const octave_idx_type dim = points.columns ();
+  const octave_idx_type num_points = points.rows ();
+
+  points = points.transpose ();
+
+  std::string options;
+
+  if (dim <= 4)
+    options = " Qbb";
+  else
+    options = " Qbb Qx";
+
+  if (nargin == 3)
+    {
+      octave_value opt_arg = args(2);
+
+      if (opt_arg.is_string ())
+        options = " " + opt_arg.string_value ();
+      else if (opt_arg.is_empty ())
+        ; // Use default options.
+      else if (opt_arg.is_cellstr ())
+        {
+          options = "";
+
+          Array<std::string> tmp = opt_arg.cellstr_value ();
+
+          for (octave_idx_type i = 0; i < tmp.numel (); i++)
+            options += " " + tmp(i);
+        }
+      else
+        {
+          error ("%s: OPTIONS must be a string, cell array of strings, or empty",
+                 caller.c_str ());
+          return retval;
+        }
+    }
+
+  boolT ismalloc = false;
+
+  unwind_protect frame;
+
+  // Replace the outfile pointer with stdout for debugging information.
+#if defined (OCTAVE_HAVE_WINDOWS_FILESYSTEM) && ! defined (OCTAVE_HAVE_POSIX_FILESYSTEM)
+  FILE *outfile = gnulib::fopen ("NUL", "w");
+#else
+  FILE *outfile = gnulib::fopen ("/dev/null", "w");
+#endif
+  FILE *errfile = stderr;
+
+  if (outfile)
+    frame.add_fcn (close_fcn, outfile);
+  else
+    {
+      error ("__voronoi__: unable to create temporary file for output");
+      return retval;
+    }
+
+  // qh_new_qhull command and points arguments are not const...
+
+  std::string cmd = "qhull v" + options;
+
+  OCTAVE_LOCAL_BUFFER (char, cmd_str, cmd.length () + 1);
+
+  strcpy (cmd_str, cmd.c_str ());
+
+  int exitcode = qh_new_qhull (dim, num_points, points.fortran_vec (),
+                               ismalloc, cmd_str, outfile, errfile);
+  if (! exitcode) 
+    {
+      // Calling findgood_all provides the number of Voronoi vertices
+      // (sets qh num_good).
+
+      qh_findgood_all (qh facet_list);
+
+      octave_idx_type num_voronoi_regions
+        = qh num_vertices - qh_setsize (qh del_vertices);
+
+      octave_idx_type num_voronoi_vertices = qh num_good;
+
+      // Find the voronoi centers for all facets.
+
+      qh_setvoronoi_all ();
+
+      facetT *facet;
+      vertexT *vertex;
+      octave_idx_type k;
+
+      // Find the number of Voronoi vertices for each Voronoi cell and
+      // store them in NI so we can use them later to set the dimensions
+      // of the RowVector objects used to collect them.
+
+      FORALLfacets
+        {
+          facet->seen = false;
+        }
+      
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, ni, num_voronoi_regions);
+      for (octave_idx_type i = 0; i < num_voronoi_regions; i++)
+        ni[i] = 0;
+
+      k = 0;
+
+      FORALLvertices
+        {
+          if (qh hull_dim == 3)
+            qh_order_vertexneighbors (vertex);
+          
+          bool infinity_seen = false;
+
+          facetT *neighbor, **neighborp;
+
+          FOREACHneighbor_ (vertex)
+            {
+              if (neighbor->upperdelaunay)
+                {
+                  if (! infinity_seen)
+                    {
+                      infinity_seen = true;
+                      ni[k]++;
+                    }
+                }
+              else
+                {
+                  neighbor->seen = true;
+                  ni[k]++;
+                }
+            }
+
+          k++;
+        }
+
+      // If Qhull finds fewer regions than points, we will pad the end
+      // of the at_inf and C arrays so that they always contain at least
+      // as many elements as the given points array.
+
+      // FIXME -- is it possible (or does it make sense) for
+      // num_voronoi_regions to ever be larger than num_points?
+
+      octave_idx_type nr = (num_points > num_voronoi_regions
+                            ? num_points : num_voronoi_regions);
+
+      boolMatrix at_inf (nr, 1, false);
+
+      // The list of Voronoi vertices.  The first element is always
+      // Inf.
+      Matrix F (num_voronoi_vertices+1, dim);
+
+      for (octave_idx_type d = 0; d < dim; d++)
+        F(0,d) = octave_Inf;
+
+      // The cell array of vectors of indices into F that represent the
+      // vertices of the Voronoi regions (cells).
+
+      Cell C (nr, 1);
+
+      // Now loop through the list of vertices again and store the
+      // coordinates of the Voronoi vertices and the lists of indices
+      // for the cells.
+
+      FORALLfacets
+        {
+          facet->seen = false;
+        }
+
+      octave_idx_type i = 0;
+      k = 0;
+
+      FORALLvertices
+        {
+          if (qh hull_dim == 3)
+            qh_order_vertexneighbors (vertex);
+
+          bool infinity_seen = false;
+
+          octave_idx_type idx = qh_pointid (vertex->point);
+
+          octave_idx_type num_vertices = ni[k++];
+
+          // Qhull seems to sometimes produces regions with a single
+          // vertex.  Is that a bug?  How can a region have just one
+          // vertex?  Let's skip it.
+
+          if (num_vertices == 1)
+            continue;
+
+          RowVector facet_list (num_vertices);
+
+          octave_idx_type m = 0;
+
+          facetT *neighbor, **neighborp;
+
+          FOREACHneighbor_(vertex)
+            {
+              if (neighbor->upperdelaunay)
+                {
+                  if (! infinity_seen)
+                    {
+                      infinity_seen = true;
+                      facet_list(m++) = 1;
+                      at_inf(idx) = true;
+                    }
+                }
+              else
+                {
+                  if (! neighbor->seen)
+                    {
+                      i++;
+                      for (octave_idx_type d = 0; d < dim; d++)
+                        F(i,d) = neighbor->center[d];
+
+                      neighbor->seen = true;
+                      neighbor->visitid = i;
+                    }
+
+                  facet_list(m++) = neighbor->visitid + 1;
+                }
+            }
+
+          C(idx) = facet_list;
+        }
+
+      retval(2) = at_inf;
+      retval(1) = C;
+      retval(0) = F;
+    }
+  else
+    error ("%s: qhull failed", caller.c_str ());
+
+  // Free memory from Qhull
+  qh_freeqhull (! qh_ALL);
+
+  int curlong, totlong;
+  qh_memfreeshort (&curlong, &totlong);
+
+  if (curlong || totlong)
+    warning ("%s: qhull did not free %d bytes of long memory (%d pieces)",
+             caller.c_str (), totlong, curlong);
+
+#else
+  error ("%s: not available in this version of Octave", caller.c_str ());
+#endif
+
+  return retval;
+}
+
+/*
+## No test needed for internal helper function.
+%!assert (1)
+*/