/*

Copyright (C) 2000, 2007 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/>.

*/

/*
29. July 2000 - Kai Habel: first release
2002-04-22 Paul Kienzle
* Use warning(...) function rather than writing to cerr
2006-05-01 Tom Holroyd
* add support for consistent winding in all dimensions; output is
* guaranteed to be simplicial.
*/

#include <sstream>

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include "oct.h"
#include "Cell.h"

#ifdef HAVE_QHULL
#if defined(HAVE__SNPRINTF) && !defined(HAVE_SNPRINTF)
#define snprintf _snprintf
#endif

extern "C" {
#include <qhull/qhull_a.h>
}

#ifdef NEED_QHULL_VERSION
char qh_version[] = "convhulln.oct 2007-07-24";
#endif
#endif

DEFUN_DLD (convhulln, args, ,
	"-*- texinfo -*-\n\
@deftypefn {Loadable Function} {@var{H} =} convhulln (@var{p})\n\
@deftypefnx {Loadable Function} {@var{H} =} convhulln (@var{p}, @var{opt})\n\
Returns 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\n\
@seealso{convhull, delaunayn}\n\
@end deftypefn")
{
  octave_value_list retval;

#ifdef HAVE_QHULL
  std::string options;

  int nargin = args.length ();
  if (nargin < 1 || nargin > 2) 
    {
      print_usage ();
      return retval;
    }

  if (nargin == 2) 
    {
      if (args (1).is_string ()) 
	options = args(1).string_value ();
      else if (args(1).is_cell ())
	{
	  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: second argument must be a string or cell array of strings");
		  return retval;
		}

	      options = options + c.elem(i).string_value() + " ";
	    }
	}
      else
	{
	  error ("convhulln: second argument must be a string or cell array of strings");
	  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;

  OCTAVE_LOCAL_BUFFER (char, flags, 250);

  // hmm, lots of options for qhull here
  // QJ guarantees that the output will be triangles
  snprintf (flags, 250, "qhull QJ %s", options.c_str ());

  if (! qh_new_qhull (dim, n, pt_array, ismalloc, flags, NULL, stderr)) 
    {
      // 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;

      Matrix idx (nf, dim);

      octave_idx_type j, i = 0;
      FORALLfacets 
	{
	  j = 0;
	  if (! facet->simplicial)
	    // should never happen with QJ
	    error ("convhulln: non-simplicial facet");

	  if (dim == 3) 
	    {
	      vertices = qh_facet3vertex (facet);
	      FOREACHvertex_ (vertices)
		idx(i, j++) = 1 + qh_pointid(vertex->point);
	      qh_settempfree (&vertices);
	    } 
	  else 
	    {
	      if (facet->toporient ^ qh_ORIENTclock) 
		{
		  FOREACHvertex_ (facet->vertices)
		    idx(i, j++) = 1 + qh_pointid(vertex->point);
		} 
	      else 
		{
		  FOREACHvertexreverse12_ (facet->vertices)
		    idx(i, j++) = 1 + qh_pointid(vertex->point);
		}
	    }
	  if (j < dim)
	    // likewise but less fatal
	    warning ("facet %d only has %d vertices", i, j);
	  i++;
	}
      retval(0) = octave_value (idx);
    }

  // free long memory
  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);
#else
  error ("convhulln: not available in this version of Octave");
#endif

  return retval;
}