#include <stdio.h> 
#include <stdlib.h> 
#include <libqhull_r/libqhull_r.h>

qhT* new_qhull_handler()
{
  qhT qh_qh;
  qhT *qh;
  FILE *errfile= stderr;

  QHULL_LIB_CHECK
  qh = (struct qhT*) malloc (sizeof (struct qhT));
  qh_zero(qh, errfile);
  return qh;
}

int delaunay_init_and_compute(qhT *qh, int dim, int numpoints, coordT *points, int* numcells, const char* flags)
{

  boolT ismalloc = False;
  FILE *errfile= stderr;
  char flagsarr[250];
  int exitcode;
  int err;
  facetT *facet;

  /* Set options */
  if (!flags)
    {
      if (dim <= 3)
        flags = "qhull d Qt Qbb Qc Qz";
      else
        flags = "qhull d Qt Qbb Qc Qx";
    }
  snprintf(flagsarr, sizeof flagsarr, "%s", flags);

  /* Run qhull*/
  exitcode = qh_new_qhull(qh, dim, numpoints, points, ismalloc, flagsarr, NULL, errfile);
  if (exitcode) return exitcode;
  /* Triangulate all the non-simplex cells generated by qhull*/
  qh_triangulate(qh);

  /* Count number of generated cells and double check that they are simplex*/
  (*numcells) = 0;
  FORALLfacets
    {
      if (! facet->upperdelaunay) (*numcells)++;
      if (! facet->simplicial) return -1;
    }

  return 0;
}

int delaunay_fill_cells(qhT *qh, int dim, int num_cells, int *cells)
{

  int icell = 0;
  facetT *facet;
  vertexT *vertex, **vertexp;
  int ivertex;

  /* Fill-in the cells connectivities */
  FORALLfacets
  {
    if (! facet->upperdelaunay)
      {
        ivertex = 0;
        FOREACHvertex_ (facet->vertices)
          {
            cells[(dim+1)*icell+ivertex] = 1 + qh_pointid(qh, vertex->point);
            ivertex++;
          }
        icell++;
      }
  }

  return 0;
}

int delaunay_free(qhT *qh)
{

  int curlong, totlong;

  /* Free memory from qhull and check that no memory is leaked */
  qh_freeqhull (qh, ! qh_ALL);
  qh_memfreeshort (qh, &curlong, &totlong);
  free(qh);
  if (curlong || totlong) return -2;

  return 0;
}

