/*

Copyright (C) 2000  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 2, 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, write to the Free
Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301, USA.

*/

/*
  16. July 2000 - Kai Habel: first release

  25. September 2002 - Changes by Rafael Laboissiere <rafael@laboissiere.net>

  * Added Qbb option to normalize the input and avoid crashes in Octave.
  * delaunayn accepts now a second (optional) argument that must be a string
  containing extra options to the qhull command.
  * Fixed doc string.  The dimension of the result matrix is [m, dim+1], and
  not [n, dim-1].

  6. June 2006: Changes by Alexander Barth <abarth@marine.usf.edu>

  * triangulate non-simplicial facets 
  * allow options to be specified as cell array of strings
  * change the default options (for compatibility with matlab)
*/

#include <iostream>
#include <string>

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include "oct.h"
#include "ov-cell.h"

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

#ifdef NEED_QHULL_VERSION
char qh_version[] = "__delaunayn__.oct 2007-08-21";
#endif
FILE *outfile = stdout;
FILE *errfile = stderr;
char flags[250];
#endif

DEFUN_DLD (__delaunayn__, args, ,
	   "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {@var{T} =} __delaunayn__ (@var{P}[, @var{opt}])\n\
Internal function for delaunayn.\n\
@end deftypefn")

{
  octave_value_list retval;
#ifdef HAVE_QHULL
  retval(0) = 0.0;
  std::string options = "";

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

  Matrix p(args(0).matrix_value());
  const octave_idx_type dim = p.columns();
  const octave_idx_type n = p.rows();

  // default options
  if (dim <= 3)
    options = "Qt Qbb Qc";
  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
	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++)
	  {

	    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__: second argument must be a string, cell of stringsor empty");
	return retval;
      }
    } 

  //octave_stdout << "options " << options << std::endl;

  if (n > dim) 
    {
      p = p.transpose();
      double *pt_array = p.fortran_vec();
      boolT ismalloc = False;

      sprintf(flags,"qhull d %s",options.c_str());

      // If you want some debugging information replace the NULL
      // pointer with outfile.

      int exitcode =  qh_new_qhull (dim, n, pt_array, ismalloc,flags, 
				    NULL, errfile);

      if (exitcode)
	{
	  error("__delaunayn__: qhull failed.");
	  return retval;
	}

      // triangulate non-simplicial facets
      qh_triangulate(); 

      facetT *facet;
      vertexT *vertex, **vertexp;
      octave_idx_type nf = 0, i = 0;

      FORALLfacets
	{
	  if (!facet->upperdelaunay)
	    nf++;

	  // Double check
	  if (!facet->simplicial) 
	    {
	      error("__delaunayn__: Qhull returned non-simplicial facets.\n",
		    "Try delaunayn with different options.");
	      break;
	    }
	}

    if (!error_state) 
      {
	Matrix simpl(nf, dim+1);
	FORALLfacets
	  {
	    if (!facet->upperdelaunay) 
	      {
		octave_idx_type j = 0;
		FOREACHvertex_ (facet->vertices)
		  {
		    // if delaunayn crashes, enable this check
#if 0
		    if (j > dim)
		      {
			error("__delaunayn__: internal error. Qhull returned non-simplicial facets.");
			return retval;
		      }
#endif

		    simpl(i, j++) = 1 + qh_pointid(vertex->point);
		  }
		i++;
	      }
	  }
	retval(0) = simpl;
      }
    

    qh_freeqhull(!qh_ALL);
    //free long memory

    int curlong, totlong;
    qh_memfreeshort (&curlong, &totlong);
    //free short memory and memory allocator

    if (curlong || totlong)
      {
	warning("__delaunay__: did not free %d bytes of long memory (%d pieces)",
		totlong, curlong);
      }
    } 
  else if (n == dim + 1) 
    {
      // one should check if nx points span a simplex
      // I will look at this later.
      RowVector vec(n);
      for (octave_idx_type i = 0; i < n; i++) 
	{
	  vec(i) = i + 1.0;
	}
      retval(0) = vec;
    }
#else
  error ("__delaunayn__: not available in this version of Octave");
#endif
  return retval;
}