diff --git a/Cdir/convexhull.h b/Cdir/convexhull.h new file mode 100644 --- /dev/null +++ b/Cdir/convexhull.h @@ -0,0 +1,59 @@ +typedef struct Vertex { + unsigned id; + double* point; +} VertexT; + +typedef struct FullVertex { + unsigned id; + double* point; + unsigned* neighfacets; + unsigned nneighfacets; + unsigned* neighvertices; + unsigned nneighsvertices; + unsigned* neighridges; + unsigned nneighridges; +} FullVertexT; + +typedef struct SetOfVertices { + VertexT* vertices; + unsigned nvertices; +} SetOfVerticesT; + +typedef struct Ridge { + VertexT* vertices; + unsigned ridgeOf1; + unsigned ridgeOf2; + unsigned nvertices; + unsigned id; + unsigned** edges; + unsigned nedges; +} RidgeT; + +typedef struct Face { + VertexT* vertices; + unsigned nvertices; + RidgeT* ridges; + unsigned* ridgesids; + unsigned nridges; + double* center; + double* normal; + double offset; + double area; + unsigned* neighbors; + unsigned neighborsize; + int family; + unsigned** edges; + unsigned nedges; +} FaceT; + +typedef struct ConvexHull { + unsigned dim; + FullVertexT* vertices; + unsigned nvertices; + FaceT* faces; + unsigned nfaces; + RidgeT* ridges; + unsigned nridges; + unsigned** edges; + unsigned nedges; +} ConvexHullT; diff --git a/Cdir/delaunay.h b/Cdir/delaunay.h new file mode 100644 --- /dev/null +++ b/Cdir/delaunay.h @@ -0,0 +1,50 @@ +typedef struct Site { + unsigned id; + unsigned* neighsites; + unsigned nneighsites; + unsigned* neighridgesids; + unsigned nneighridges; + unsigned* neightiles; + unsigned nneightiles; +} SiteT; + +typedef struct Simplex { + unsigned* sitesids; + double* center; + double radius; + double volume; +} SimplexT; + +typedef struct SubTile { + unsigned id; + SimplexT simplex; + unsigned ridgeOf1; + int ridgeOf2; + double* normal; + double offset; + unsigned flag; +} SubTileT; + +#define INIT_SUBTILE(X) SubTileT X = {.flag = 0} + +typedef struct Tile { + unsigned id; + SimplexT simplex; + unsigned* neighbors; + unsigned nneighbors; + unsigned* ridgesids; + unsigned nridges; // = dim+1 + int family; + int orientation; +} TileT; + +typedef struct Tesselation { + SiteT* sites; + TileT* tiles; + unsigned ntiles; + SubTileT* subtiles; + unsigned nsubtiles; +} TesselationT; + +TesselationT* tesselation(double*, unsigned, unsigned, unsigned, unsigned, double, unsigned*); +void testdel2(); diff --git a/Cdir/geom_r.h b/Cdir/geom_r.h new file mode 100644 --- /dev/null +++ b/Cdir/geom_r.h @@ -0,0 +1,184 @@ +/*
  ---------------------------------
+
+  geom_r.h
+    header file for geometric routines
+
+   see qh-geom_r.htm and geom_r.c
+
+   Copyright (c) 1993-2015 The Geometry Center.
+   $Id: //main/2015/qhull/src/libqhull_r/geom_r.h#3 $$Change: 2079 $
+   $DateTime: 2016/02/07 17:43:34 $$Author: bbarber $
+*/
+
+#ifndef qhDEFgeom
+#define qhDEFgeom 1
+
+#include "libqhull_r.h"
+
+/* ============ -macros- ======================== */
+
+/*----------------------------------
+
+  fabs_(a)
+    returns the absolute value of a
+*/
+#define fabs_( a ) ((( a ) < 0 ) ? -( a ):( a ))
+
+/*----------------------------------
+
+  fmax_(a,b)
+    returns the maximum value of a and b
+*/
+#define fmax_( a,b )  ( ( a ) < ( b ) ? ( b ) : ( a ) )
+
+/*----------------------------------
+
+  fmin_(a,b)
+    returns the minimum value of a and b
+*/
+#define fmin_( a,b )  ( ( a ) > ( b ) ? ( b ) : ( a ) )
+
+/*----------------------------------
+
+  maximize_(maxval, val)
+    set maxval to val if val is greater than maxval
+*/
+#define maximize_( maxval, val ) { if (( maxval ) < ( val )) ( maxval )= ( val ); }
+
+/*----------------------------------
+
+  minimize_(minval, val)
+    set minval to val if val is less than minval
+*/
+#define minimize_( minval, val ) { if (( minval ) > ( val )) ( minval )= ( val ); }
+
+/*----------------------------------
+
+  det2_(a1, a2,
+        b1, b2)
+
+    compute a 2-d determinate
+*/
+#define det2_( a1,a2,b1,b2 ) (( a1 )*( b2 ) - ( a2 )*( b1 ))
+
+/*----------------------------------
+
+  det3_(a1, a2, a3,
+       b1, b2, b3,
+       c1, c2, c3)
+
+    compute a 3-d determinate
+*/
+#define det3_( a1,a2,a3,b1,b2,b3,c1,c2,c3 ) ( ( a1 )*det2_( b2,b3,c2,c3 ) \
+                - ( b1 )*det2_( a2,a3,c2,c3 ) + ( c1 )*det2_( a2,a3,b2,b3 ) )
+
+/*----------------------------------
+
+  dX( p1, p2 )
+  dY( p1, p2 )
+  dZ( p1, p2 )
+
+    given two indices into rows[],
+
+    compute the difference between X, Y, or Z coordinates
+*/
+#define dX( p1,p2 )  ( *( rows[p1] ) - *( rows[p2] ))
+#define dY( p1,p2 )  ( *( rows[p1]+1 ) - *( rows[p2]+1 ))
+#define dZ( p1,p2 )  ( *( rows[p1]+2 ) - *( rows[p2]+2 ))
+#define dW( p1,p2 )  ( *( rows[p1]+3 ) - *( rows[p2]+3 ))
+
+/*============= prototypes in alphabetical order, infrequent at end ======= */
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void    qh_backnormal(qhT *qh, realT **rows, int numrow, int numcol, boolT sign, coordT *normal, boolT *nearzero);
+void    qh_distplane(qhT *qh, pointT *point, facetT *facet, realT *dist);
+facetT *qh_findbest(qhT *qh, pointT *point, facetT *startfacet,
+                     boolT bestoutside, boolT isnewfacets, boolT noupper,
+                     realT *dist, boolT *isoutside, int *numpart);
+facetT *qh_findbesthorizon(qhT *qh, boolT ischeckmax, pointT *point,
+                     facetT *startfacet, boolT noupper, realT *bestdist, int *numpart);
+facetT *qh_findbestnew(qhT *qh, pointT *point, facetT *startfacet, realT *dist,
+                     boolT bestoutside, boolT *isoutside, int *numpart);
+void    qh_gausselim(qhT *qh, realT **rows, int numrow, int numcol, boolT *sign, boolT *nearzero);
+realT   qh_getangle(qhT *qh, pointT *vect1, pointT *vect2);
+pointT *qh_getcenter(qhT *qh, setT *vertices);
+pointT *qh_getcentrum(qhT *qh, facetT *facet);
+realT   qh_getdistance(qhT *qh, facetT *facet, facetT *neighbor, realT *mindist, realT *maxdist);
+void    qh_normalize(qhT *qh, coordT *normal, int dim, boolT toporient);
+void    qh_normalize2(qhT *qh, coordT *normal, int dim, boolT toporient,
+            realT *minnorm, boolT *ismin);
+pointT *qh_projectpoint(qhT *qh, pointT *point, facetT *facet, realT dist);
+
+void    qh_setfacetplane(qhT *qh, facetT *newfacets);
+void    qh_sethyperplane_det(qhT *qh, int dim, coordT **rows, coordT *point0,
+              boolT toporient, coordT *normal, realT *offset, boolT *nearzero);
+void    qh_sethyperplane_gauss(qhT *qh, int dim, coordT **rows, pointT *point0,
+             boolT toporient, coordT *normal, coordT *offset, boolT *nearzero);
+boolT   qh_sharpnewfacets(qhT *qh);
+
+/*========= infrequently used code in geom2_r.c =============*/
+
+coordT *qh_copypoints(qhT *qh, coordT *points, int numpoints, int dimension);
+void    qh_crossproduct(int dim, realT vecA[3], realT vecB[3], realT vecC[3]);
+realT   qh_determinant(qhT *qh, realT **rows, int dim, boolT *nearzero);
+realT   qh_detjoggle(qhT *qh, pointT *points, int numpoints, int dimension);
+void    qh_detroundoff(qhT *qh);
+realT   qh_detsimplex(qhT *qh, pointT *apex, setT *points, int dim, boolT *nearzero);
+realT   qh_distnorm(int dim, pointT *point, pointT *normal, realT *offsetp);
+realT   qh_distround(qhT *qh, int dimension, realT maxabs, realT maxsumabs);
+realT   qh_divzero(realT numer, realT denom, realT mindenom1, boolT *zerodiv);
+realT   qh_facetarea(qhT *qh, facetT *facet);
+realT   qh_facetarea_simplex(qhT *qh, int dim, coordT *apex, setT *vertices,
+          vertexT *notvertex,  boolT toporient, coordT *normal, realT *offset);
+pointT *qh_facetcenter(qhT *qh, setT *vertices);
+facetT *qh_findgooddist(qhT *qh, pointT *point, facetT *facetA, realT *distp, facetT **facetlist);
+void    qh_getarea(qhT *qh, facetT *facetlist);
+boolT   qh_gram_schmidt(qhT *qh, int dim, realT **rows);
+boolT   qh_inthresholds(qhT *qh, coordT *normal, realT *angle);
+void    qh_joggleinput(qhT *qh);
+realT  *qh_maxabsval(realT *normal, int dim);
+setT   *qh_maxmin(qhT *qh, pointT *points, int numpoints, int dimension);
+realT   qh_maxouter(qhT *qh);
+void    qh_maxsimplex(qhT *qh, int dim, setT *maxpoints, pointT *points, int numpoints, setT **simplex);
+realT   qh_minabsval(realT *normal, int dim);
+int     qh_mindiff(realT *vecA, realT *vecB, int dim);
+boolT   qh_orientoutside(qhT *qh, facetT *facet);
+void    qh_outerinner(qhT *qh, facetT *facet, realT *outerplane, realT *innerplane);
+coordT  qh_pointdist(pointT *point1, pointT *point2, int dim);
+void    qh_printmatrix(qhT *qh, FILE *fp, const char *string, realT **rows, int numrow, int numcol);
+void    qh_printpoints(qhT *qh, FILE *fp, const char *string, setT *points);
+void    qh_projectinput(qhT *qh);
+void    qh_projectpoints(qhT *qh, signed char *project, int n, realT *points,
+             int numpoints, int dim, realT *newpoints, int newdim);
+void    qh_rotateinput(qhT *qh, realT **rows);
+void    qh_rotatepoints(qhT *qh, realT *points, int numpoints, int dim, realT **rows);
+void    qh_scaleinput(qhT *qh);
+void    qh_scalelast(qhT *qh, coordT *points, int numpoints, int dim, coordT low,
+                   coordT high, coordT newhigh);
+void    qh_scalepoints(qhT *qh, pointT *points, int numpoints, int dim,
+                realT *newlows, realT *newhighs);
+boolT   qh_sethalfspace(qhT *qh, int dim, coordT *coords, coordT **nextp,
+              coordT *normal, coordT *offset, coordT *feasible);
+coordT *qh_sethalfspace_all(qhT *qh, int dim, int count, coordT *halfspaces, pointT *feasible);
+pointT *qh_voronoi_center(qhT *qh, int dim, setT *points);
+
+#ifdef __cplusplus
+} /* extern "C"*/
+#endif
+
+#endif /* qhDEFgeom */
+
+
+
diff --git a/Cdir/io_r.h b/Cdir/io_r.h
new file mode 100644
--- /dev/null
+++ b/Cdir/io_r.h
@@ -0,0 +1,167 @@
+/*
  ---------------------------------
+
+   io_r.h
+   declarations of Input/Output functions
+
+   see README, libqhull_r.h and io_r.c
+
+   Copyright (c) 1993-2015 The Geometry Center.
+   $Id: //main/2015/qhull/src/libqhull_r/io_r.h#3 $$Change: 2079 $
+   $DateTime: 2016/02/07 17:43:34 $$Author: bbarber $
+*/
+
+#ifndef qhDEFio
+#define qhDEFio 1
+
+#include "libqhull_r.h"
+
+/*============ constants and flags ==================*/
+
+/*----------------------------------
+
+  qh_MAXfirst
+    maximum length of first two lines of stdin
+*/
+#define qh_MAXfirst  200
+
+/*----------------------------------
+
+  qh_MINradius
+    min radius for Gp and Gv, fraction of maxcoord
+*/
+#define qh_MINradius 0.02
+
+/*----------------------------------
+
+  qh_GEOMepsilon
+    adjust outer planes for 'lines closer' and geomview roundoff.
+    This prevents bleed through.
+*/
+#define qh_GEOMepsilon 2e-3
+
+/*----------------------------------
+
+  qh_WHITESPACE
+    possible values of white space
+*/
+#define qh_WHITESPACE " \n\t\v\r\f"
+
+
+/*----------------------------------
+
+  qh_RIDGE
+    to select which ridges to print in qh_eachvoronoi
+*/
+typedef enum
+{
+    qh_RIDGEall = 0, qh_RIDGEinner, qh_RIDGEouter
+}
+qh_RIDGE;
+
+/*----------------------------------
+
+  printvridgeT
+    prints results of qh_printvdiagram
+
+  see:
+    qh_printvridge for an example
+*/
+typedef void (*printvridgeT)(qhT *qh, FILE *fp, vertexT *vertex, vertexT *vertexA, setT *centers, boolT unbounded);
+
+/*============== -prototypes in alphabetical order =========*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void    qh_dfacet(qhT *qh, unsigned id);
+void    qh_dvertex(qhT *qh, unsigned id);
+int     qh_compare_facetarea(const void *p1, const void *p2);
+int     qh_compare_facetmerge(const void *p1, const void *p2);
+int     qh_compare_facetvisit(const void *p1, const void *p2);
+/* int  qh_compare_vertexpoint(const void *p1, const void *p2); Not useable since it depends on qh */
+void    qh_copyfilename(qhT *qh, char *filename, int size, const char* source, int length);
+void    qh_countfacets(qhT *qh, facetT *facetlist, setT *facets, boolT printall,
+              int *numfacetsp, int *numsimplicialp, int *totneighborsp,
+              int *numridgesp, int *numcoplanarsp, int *numnumtricoplanarsp);
+pointT *qh_detvnorm(qhT *qh, vertexT *vertex, vertexT *vertexA, setT *centers, realT *offsetp);
+setT   *qh_detvridge(qhT *qh, vertexT *vertex);
+setT   *qh_detvridge3(qhT *qh, vertexT *atvertex, vertexT *vertex);
+int     qh_eachvoronoi(qhT *qh, FILE *fp, printvridgeT printvridge, vertexT *atvertex, boolT visitall, qh_RIDGE innerouter, boolT inorder);
+int     qh_eachvoronoi_all(qhT *qh, FILE *fp, printvridgeT printvridge, boolT isUpper, qh_RIDGE innerouter, boolT inorder);
+void    qh_facet2point(qhT *qh, facetT *facet, pointT **point0, pointT **point1, realT *mindist);
+setT   *qh_facetvertices(qhT *qh, facetT *facetlist, setT *facets, boolT allfacets);
+void    qh_geomplanes(qhT *qh, facetT *facet, realT *outerplane, realT *innerplane);
+void    qh_markkeep(qhT *qh, facetT *facetlist);
+setT   *qh_markvoronoi(qhT *qh, facetT *facetlist, setT *facets, boolT printall, boolT *isLowerp, int *numcentersp);
+void    qh_order_vertexneighbors(qhT *qh, vertexT *vertex);
+void    qh_prepare_output(qhT *qh);
+void    qh_printafacet(qhT *qh, FILE *fp, qh_PRINT format, facetT *facet, boolT printall);
+void    qh_printbegin(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall);
+void    qh_printcenter(qhT *qh, FILE *fp, qh_PRINT format, const char *string, facetT *facet);
+void    qh_printcentrum(qhT *qh, FILE *fp, facetT *facet, realT radius);
+void    qh_printend(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall);
+void    qh_printend4geom(qhT *qh, FILE *fp, facetT *facet, int *num, boolT printall);
+void    qh_printextremes(qhT *qh, FILE *fp, facetT *facetlist, setT *facets, boolT printall);
+void    qh_printextremes_2d(qhT *qh, FILE *fp, facetT *facetlist, setT *facets, boolT printall);
+void    qh_printextremes_d(qhT *qh, FILE *fp, facetT *facetlist, setT *facets, boolT printall);
+void    qh_printfacet(qhT *qh, FILE *fp, facetT *facet);
+void    qh_printfacet2math(qhT *qh, FILE *fp, facetT *facet, qh_PRINT format, int notfirst);
+void    qh_printfacet2geom(qhT *qh, FILE *fp, facetT *facet, realT color[3]);
+void    qh_printfacet2geom_points(qhT *qh, FILE *fp, pointT *point1, pointT *point2,
+                               facetT *facet, realT offset, realT color[3]);
+void    qh_printfacet3math(qhT *qh, FILE *fp, facetT *facet, qh_PRINT format, int notfirst);
+void    qh_printfacet3geom_nonsimplicial(qhT *qh, FILE *fp, facetT *facet, realT color[3]);
+void    qh_printfacet3geom_points(qhT *qh, FILE *fp, setT *points, facetT *facet, realT offset, realT color[3]);
+void    qh_printfacet3geom_simplicial(qhT *qh, FILE *fp, facetT *facet, realT color[3]);
+void    qh_printfacet3vertex(qhT *qh, FILE *fp, facetT *facet, qh_PRINT format);
+void    qh_printfacet4geom_nonsimplicial(qhT *qh, FILE *fp, facetT *facet, realT color[3]);
+void    qh_printfacet4geom_simplicial(qhT *qh, FILE *fp, facetT *facet, realT color[3]);
+void    qh_printfacetNvertex_nonsimplicial(qhT *qh, FILE *fp, facetT *facet, int id, qh_PRINT format);
+void    qh_printfacetNvertex_simplicial(qhT *qh, FILE *fp, facetT *facet, qh_PRINT format);
+void    qh_printfacetheader(qhT *qh, FILE *fp, facetT *facet);
+void    qh_printfacetridges(qhT *qh, FILE *fp, facetT *facet);
+void    qh_printfacets(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall);
+void    qh_printhyperplaneintersection(qhT *qh, FILE *fp, facetT *facet1, facetT *facet2,
+                   setT *vertices, realT color[3]);
+void    qh_printneighborhood(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetA, facetT *facetB, boolT printall);
+void    qh_printline3geom(qhT *qh, FILE *fp, pointT *pointA, pointT *pointB, realT color[3]);
+void    qh_printpoint(qhT *qh, FILE *fp, const char *string, pointT *point);
+void    qh_printpointid(qhT *qh, FILE *fp, const char *string, int dim, pointT *point, int id);
+void    qh_printpoint3(qhT *qh, FILE *fp, pointT *point);
+void    qh_printpoints_out(qhT *qh, FILE *fp, facetT *facetlist, setT *facets, boolT printall);
+void    qh_printpointvect(qhT *qh, FILE *fp, pointT *point, coordT *normal, pointT *center, realT radius, realT color[3]);
+void    qh_printpointvect2(qhT *qh, FILE *fp, pointT *point, coordT *normal, pointT *center, realT radius);
+void    qh_printridge(qhT *qh, FILE *fp, ridgeT *ridge);
+void    qh_printspheres(qhT *qh, FILE *fp, setT *vertices, realT radius);
+void    qh_printvdiagram(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall);
+int     qh_printvdiagram2(qhT *qh, FILE *fp, printvridgeT printvridge, setT *vertices, qh_RIDGE innerouter, boolT inorder);
+void    qh_printvertex(qhT *qh, FILE *fp, vertexT *vertex);
+void    qh_printvertexlist(qhT *qh, FILE *fp, const char* string, facetT *facetlist,
+                         setT *facets, boolT printall);
+void    qh_printvertices(qhT *qh, FILE *fp, const char* string, setT *vertices);
+void    qh_printvneighbors(qhT *qh, FILE *fp, facetT* facetlist, setT *facets, boolT printall);
+void    qh_printvoronoi(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall);
+void    qh_printvnorm(qhT *qh, FILE *fp, vertexT *vertex, vertexT *vertexA, setT *centers, boolT unbounded);
+void    qh_printvridge(qhT *qh, FILE *fp, vertexT *vertex, vertexT *vertexA, setT *centers, boolT unbounded);
+void    qh_produce_output(qhT *qh);
+void    qh_produce_output2(qhT *qh);
+void    qh_projectdim3(qhT *qh, pointT *source, pointT *destination);
+int     qh_readfeasible(qhT *qh, int dim, const char *curline);
+coordT *qh_readpoints(qhT *qh, int *numpoints, int *dimension, boolT *ismalloc);
+void    qh_setfeasible(qhT *qh, int dim);
+boolT   qh_skipfacet(qhT *qh, facetT *facet);
+char   *qh_skipfilename(qhT *qh, char *filename);
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /* qhDEFio */
diff --git a/Cdir/libqhull_r.h b/Cdir/libqhull_r.h
new file mode 100644
--- /dev/null
+++ b/Cdir/libqhull_r.h
@@ -0,0 +1,1134 @@
+/*
  ---------------------------------
+
+   libqhull_r.h
+   user-level header file for using qhull.a library
+
+   see qh-qhull_r.htm, qhull_ra.h
+
+   Copyright (c) 1993-2015 The Geometry Center.
+   $Id: //main/2015/qhull/src/libqhull_r/libqhull_r.h#8 $$Change: 2079 $
+   $DateTime: 2016/02/07 17:43:34 $$Author: bbarber $
+
+   includes function prototypes for libqhull_r.c, geom_r.c, global_r.c, io_r.c, user.c
+
+   use mem_r.h for mem_r.c
+   use qset_r.h for qset_r.c
+
+   see unix_r.c for an example of using libqhull_r.h
+
+   recompile qhull if you change this file
+*/
+
+#ifndef qhDEFlibqhull
+#define qhDEFlibqhull 1
+
+/*=========================== -included files ==============*/
+
+/* user_r.h first for QHULL_CRTDBG */
+#include "user_r.h"      /* user definable constants (e.g., realT). */
+
+#include "mem_r.h"   /* Needed for qhT in libqhull_r.h */
+#include "qset_r.h"   /* Needed for QHULL_LIB_CHECK */
+/* include stat_r.h after defining boolT.  statT needed for qhT in libqhull_r.h */
+
+#include 
+#include 
+#include 
+#include 
+
+#ifndef __STDC__
+#ifndef __cplusplus
+#if     !_MSC_VER
+#error  Neither __STDC__ nor __cplusplus is defined.  Please use strict ANSI C or C++ to compile
+#error  Qhull.  You may need to turn off compiler extensions in your project configuration.  If
+#error  your compiler is a standard C compiler, you can delete this warning from libqhull_r.h
+#endif
+#endif
+#endif
+
+/*============ constants and basic types ====================*/
+
+extern const char qh_version[]; /* defined in global_r.c */
+extern const char qh_version2[]; /* defined in global_r.c */
+
+/*----------------------------------
+
+  coordT
+    coordinates and coefficients are stored as realT (i.e., double)
+
+  notes:
+    Qhull works well if realT is 'float'.  If so joggle (QJ) is not effective.
+
+    Could use 'float' for data and 'double' for calculations (realT vs. coordT)
+      This requires many type casts, and adjusted error bounds.
+      Also C compilers may do expressions in double anyway.
+*/
+#define coordT realT
+
+/*----------------------------------
+
+  pointT
+    a point is an array of coordinates, usually qh.hull_dim
+    qh_pointid returns
+      qh_IDnone if point==0 or qh is undefined
+      qh_IDinterior for qh.interior_point
+      qh_IDunknown if point is neither in qh.first_point... nor qh.other_points
+
+  notes:
+    qh.STOPcone and qh.STOPpoint assume that qh_IDunknown==-1 (other negative numbers indicate points)
+    qh_IDunknown is also returned by getid_() for unknown facet, ridge, or vertex
+*/
+#define pointT coordT
+typedef enum
+{
+    qh_IDnone = -3, qh_IDinterior = -2, qh_IDunknown = -1
+}
+qh_pointT;
+
+/*----------------------------------
+
+  flagT
+    Boolean flag as a bit
+*/
+#define flagT unsigned int
+
+/*----------------------------------
+
+  boolT
+    boolean value, either True or False
+
+  notes:
+    needed for portability
+    Use qh_False/qh_True as synonyms
+*/
+#define boolT unsigned int
+#ifdef False
+#undef False
+#endif
+#ifdef True
+#undef True
+#endif
+#define False 0
+#define True 1
+#define qh_False 0
+#define qh_True 1
+
+#include "stat_r.h"  /* needs boolT */
+
+/*----------------------------------
+
+  qh_CENTER
+    to distinguish facet->center
+*/
+typedef enum
+{
+    qh_ASnone = 0,   /* If not MERGING and not VORONOI */
+    qh_ASvoronoi,    /* Set by qh_clearcenters on qh_prepare_output, or if not MERGING and VORONOI */
+    qh_AScentrum     /* If MERGING (assumed during merging) */
+}
+qh_CENTER;
+
+/*----------------------------------
+
+  qh_PRINT
+    output formats for printing (qh.PRINTout).
+    'Fa' 'FV' 'Fc' 'FC'
+
+
+   notes:
+   some of these names are similar to qhT names.  The similar names are only
+   used in switch statements in qh_printbegin() etc.
+*/
+typedef enum {qh_PRINTnone= 0,
+  qh_PRINTarea, qh_PRINTaverage,           /* 'Fa' 'FV' 'Fc' 'FC' */
+  qh_PRINTcoplanars, qh_PRINTcentrums,
+  qh_PRINTfacets, qh_PRINTfacets_xridge,   /* 'f' 'FF' 'G' 'FI' 'Fi' 'Fn' */
+  qh_PRINTgeom, qh_PRINTids, qh_PRINTinner, qh_PRINTneighbors,
+  qh_PRINTnormals, qh_PRINTouter, qh_PRINTmaple, /* 'n' 'Fo' 'i' 'm' 'Fm' 'FM', 'o' */
+  qh_PRINTincidences, qh_PRINTmathematica, qh_PRINTmerges, qh_PRINToff,
+  qh_PRINToptions, qh_PRINTpointintersect, /* 'FO' 'Fp' 'FP' 'p' 'FQ' 'FS' */
+  qh_PRINTpointnearest, qh_PRINTpoints, qh_PRINTqhull, qh_PRINTsize,
+  qh_PRINTsummary, qh_PRINTtriangles,      /* 'Fs' 'Ft' 'Fv' 'FN' 'Fx' */
+  qh_PRINTvertices, qh_PRINTvneighbors, qh_PRINTextremes,
+  qh_PRINTEND} qh_PRINT;
+
+/*----------------------------------
+
+  qh_ALL
+    argument flag for selecting everything
+*/
+#define qh_ALL      True
+#define qh_NOupper  True     /* argument for qh_findbest */
+#define qh_IScheckmax  True     /* argument for qh_findbesthorizon */
+#define qh_ISnewfacets  True     /* argument for qh_findbest */
+#define qh_RESETvisible  True     /* argument for qh_resetlists */
+
+/*----------------------------------
+
+  qh_ERR
+    Qhull exit codes, for indicating errors
+    See: MSG_ERROR and MSG_WARNING [user.h]
+*/
+#define qh_ERRnone  0    /* no error occurred during qhull */
+#define qh_ERRinput 1    /* input inconsistency */
+#define qh_ERRsingular 2 /* singular input data */
+#define qh_ERRprec  3    /* precision error */
+#define qh_ERRmem   4    /* insufficient memory, matches mem_r.h */
+#define qh_ERRqhull 5    /* internal error detected, matches mem_r.h */
+
+/*----------------------------------
+
+qh_FILEstderr
+Fake stderr to distinguish error output from normal output
+For C++ interface.  Must redefine qh_fprintf_qhull
+*/
+#define qh_FILEstderr ((FILE*)1)
+
+/* ============ -structures- ====================
+   each of the following structures is defined by a typedef
+   all realT and coordT fields occur at the beginning of a structure
+        (otherwise space may be wasted due to alignment)
+   define all flags together and pack into 32-bit number
+
+   DEFqhT and DEFsetT are likewise defined in
+   mem_r.h, qset_r.h, and stat_r.h.
+
+*/
+
+typedef struct vertexT vertexT;
+typedef struct ridgeT ridgeT;
+typedef struct facetT facetT;
+
+#ifndef DEFqhT
+#define DEFqhT 1
+typedef struct qhT qhT;          /* defined below */
+#endif
+
+#ifndef DEFsetT
+#define DEFsetT 1
+typedef struct setT setT;          /* defined in qset_r.h */
+#endif
+
+/*----------------------------------
+
+  facetT
+    defines a facet
+
+  notes:
+   qhull() generates the hull as a list of facets.
+
+  topological information:
+    f.previous,next     doubly-linked list of facets
+    f.vertices          set of vertices
+    f.ridges            set of ridges
+    f.neighbors         set of neighbors
+    f.toporient         True if facet has top-orientation (else bottom)
+
+  geometric information:
+    f.offset,normal     hyperplane equation
+    f.maxoutside        offset to outer plane -- all points inside
+    f.center            centrum for testing convexity or Voronoi center for output
+    f.simplicial        True if facet is simplicial
+    f.flipped           True if facet does not include qh.interior_point
+
+  for constructing hull:
+    f.visible           True if facet on list of visible facets (will be deleted)
+    f.newfacet          True if facet on list of newly created facets
+    f.coplanarset       set of points coplanar with this facet
+                        (includes near-inside points for later testing)
+    f.outsideset        set of points outside of this facet
+    f.furthestdist      distance to furthest point of outside set
+    f.visitid           marks visited facets during a loop
+    f.replace           replacement facet for to-be-deleted, visible facets
+    f.samecycle,newcycle cycle of facets for merging into horizon facet
+
+  see below for other flags and fields
+*/
+struct facetT {
+#if !qh_COMPUTEfurthest
+  coordT   furthestdist;/* distance to furthest point of outsideset */
+#endif
+#if qh_MAXoutside
+  coordT   maxoutside;  /* max computed distance of point to facet
+                        Before QHULLfinished this is an approximation
+                        since maxdist not always set for mergefacet
+                        Actual outer plane is +DISTround and
+                        computed outer plane is +2*DISTround */
+#endif
+  coordT   offset;      /* exact offset of hyperplane from origin */
+  coordT  *normal;      /* normal of hyperplane, hull_dim coefficients */
+                        /*   if ->tricoplanar, shared with a neighbor */
+  union {               /* in order of testing */
+   realT   area;        /* area of facet, only in io_r.c if  ->isarea */
+   facetT *replace;     /*  replacement facet if ->visible and NEWfacets
+                             is NULL only if qh_mergedegen_redundant or interior */
+   facetT *samecycle;   /*  cycle of facets from the same visible/horizon intersection,
+                             if ->newfacet */
+   facetT *newcycle;    /*  in horizon facet, current samecycle of new facets */
+   facetT *trivisible;  /* visible facet for ->tricoplanar facets during qh_triangulate() */
+   facetT *triowner;    /* owner facet for ->tricoplanar, !isarea facets w/ ->keepcentrum */
+  }f;
+  coordT  *center;      /* set according to qh.CENTERtype */
+                        /*   qh_ASnone:    no center (not MERGING) */
+                        /*   qh_AScentrum: centrum for testing convexity (qh_getcentrum) */
+                        /*                 assumed qh_AScentrum while merging */
+                        /*   qh_ASvoronoi: Voronoi center (qh_facetcenter) */
+                        /* after constructing the hull, it may be changed (qh_clearcenter) */
+                        /* if tricoplanar and !keepcentrum, shared with a neighbor */
+  facetT  *previous;    /* previous facet in the facet_list */
+  facetT  *next;        /* next facet in the facet_list */
+  setT    *vertices;    /* vertices for this facet, inverse sorted by ID
+                           if simplicial, 1st vertex was apex/furthest */
+  setT    *ridges;      /* explicit ridges for nonsimplicial facets.
+                           for simplicial facets, neighbors define the ridges */
+  setT    *neighbors;   /* neighbors of the facet.  If simplicial, the kth
+                           neighbor is opposite the kth vertex, and the first
+                           neighbor is the horizon facet for the first vertex*/
+  setT    *outsideset;  /* set of points outside this facet
+                           if non-empty, last point is furthest
+                           if NARROWhull, includes coplanars for partitioning*/
+  setT    *coplanarset; /* set of points coplanar with this facet
+                           > qh.min_vertex and <= facet->max_outside
+                           a point is assigned to the furthest facet
+                           if non-empty, last point is furthest away */
+  unsigned visitid;     /* visit_id, for visiting all neighbors,
+                           all uses are independent */
+  unsigned id;          /* unique identifier from qh.facet_id */
+  unsigned nummerge:9;  /* number of merges */
+#define qh_MAXnummerge 511 /*     2^9-1, 32 flags total, see "flags:" in io_r.c */
+  flagT    tricoplanar:1; /* True if TRIangulate and simplicial and coplanar with a neighbor */
+                          /*   all tricoplanars share the same apex */
+                          /*   all tricoplanars share the same ->center, ->normal, ->offset, ->maxoutside */
+                          /*     ->keepcentrum is true for the owner.  It has the ->coplanareset */
+                          /*   if ->degenerate, does not span facet (one logical ridge) */
+                          /*   during qh_triangulate, f.trivisible points to original facet */
+  flagT    newfacet:1;  /* True if facet on qh.newfacet_list (new or merged) */
+  flagT    visible:1;   /* True if visible facet (will be deleted) */
+  flagT    toporient:1; /* True if created with top orientation
+                           after merging, use ridge orientation */
+  flagT    simplicial:1;/* True if simplicial facet, ->ridges may be implicit */
+  flagT    seen:1;      /* used to perform operations only once, like visitid */
+  flagT    seen2:1;     /* used to perform operations only once, like visitid */
+  flagT    flipped:1;   /* True if facet is flipped */
+  flagT    upperdelaunay:1; /* True if facet is upper envelope of Delaunay triangulation */
+  flagT    notfurthest:1; /* True if last point of outsideset is not furthest*/
+
+/*-------- flags primarily for output ---------*/
+  flagT    good:1;      /* True if a facet marked good for output */
+  flagT    isarea:1;    /* True if facet->f.area is defined */
+
+/*-------- flags for merging ------------------*/
+  flagT    dupridge:1;  /* True if duplicate ridge in facet */
+  flagT    mergeridge:1; /* True if facet or neighbor contains a qh_MERGEridge
+                            ->normal defined (also defined for mergeridge2) */
+  flagT    mergeridge2:1; /* True if neighbor contains a qh_MERGEridge (qhT *qh, mark_dupridges */
+  flagT    coplanar:1;  /* True if horizon facet is coplanar at last use */
+  flagT     mergehorizon:1; /* True if will merge into horizon (->coplanar) */
+  flagT     cycledone:1;/* True if mergecycle_all already done */
+  flagT    tested:1;    /* True if facet convexity has been tested (false after merge */
+  flagT    keepcentrum:1; /* True if keep old centrum after a merge, or marks owner for ->tricoplanar */
+  flagT    newmerge:1;  /* True if facet is newly merged for reducevertices */
+  flagT    degenerate:1; /* True if facet is degenerate (degen_mergeset or ->tricoplanar) */
+  flagT    redundant:1;  /* True if facet is redundant (degen_mergeset) */
+};
+
+
+/*----------------------------------
+
+  ridgeT
+    defines a ridge
+
+  notes:
+  a ridge is hull_dim-1 simplex between two neighboring facets.  If the
+  facets are non-simplicial, there may be more than one ridge between
+  two facets.  E.G. a 4-d hypercube has two triangles between each pair
+  of neighboring facets.
+
+  topological information:
+    vertices            a set of vertices
+    top,bottom          neighboring facets with orientation
+
+  geometric information:
+    tested              True if ridge is clearly convex
+    nonconvex           True if ridge is non-convex
+*/
+struct ridgeT {
+  setT    *vertices;    /* vertices belonging to this ridge, inverse sorted by ID
+                           NULL if a degen ridge (matchsame) */
+  facetT  *top;         /* top facet this ridge is part of */
+  facetT  *bottom;      /* bottom facet this ridge is part of */
+  unsigned id;          /* unique identifier.  Same size as vertex_id and ridge_id */
+  flagT    seen:1;      /* used to perform operations only once */
+  flagT    tested:1;    /* True when ridge is tested for convexity */
+  flagT    nonconvex:1; /* True if getmergeset detected a non-convex neighbor
+                           only one ridge between neighbors may have nonconvex */
+};
+
+/*----------------------------------
+
+  vertexT
+     defines a vertex
+
+  topological information:
+    next,previous       doubly-linked list of all vertices
+    neighbors           set of adjacent facets (only if qh.VERTEXneighbors)
+
+  geometric information:
+    point               array of DIM3 coordinates
+*/
+struct vertexT {
+  vertexT *next;        /* next vertex in vertex_list */
+  vertexT *previous;    /* previous vertex in vertex_list */
+  pointT  *point;       /* hull_dim coordinates (coordT) */
+  setT    *neighbors;   /* neighboring facets of vertex, qh_vertexneighbors()
+                           inits in io_r.c or after first merge */
+  unsigned id;          /* unique identifier.  Same size as qh.vertex_id and qh.ridge_id */
+  unsigned visitid;    /* for use with qh.vertex_visit, size must match */
+  flagT    seen:1;      /* used to perform operations only once */
+  flagT    seen2:1;     /* another seen flag */
+  flagT    delridge:1;  /* vertex was part of a deleted ridge */
+  flagT    deleted:1;   /* true if vertex on qh.del_vertices */
+  flagT    newlist:1;   /* true if vertex on qh.newvertex_list */
+};
+
+/*======= -global variables -qh ============================*/
+
+/*----------------------------------
+
+  qhT
+   All global variables for qhull are in qhT.  It includes qhmemT, qhstatT, and rbox globals
+
+   This version of Qhull is reentrant, but it is not thread-safe.
+
+   Do not run separate threads on the same instance of qhT.
+
+   QHULL_LIB_CHECK checks that a program and the corresponding
+   qhull library were built with the same type of header files.
+*/
+
+#define QHULL_NON_REENTRANT 0
+#define QHULL_QH_POINTER 1
+#define QHULL_REENTRANT 2
+
+#define QHULL_LIB_TYPE QHULL_REENTRANT
+
+#define QHULL_LIB_CHECK qh_lib_check(QHULL_LIB_TYPE, sizeof(qhT), sizeof(vertexT), sizeof(ridgeT), sizeof(facetT), sizeof(setT), sizeof(qhmemT));
+#define QHULL_LIB_CHECK_RBOX qh_lib_check(QHULL_LIB_TYPE, sizeof(qhT), sizeof(vertexT), sizeof(ridgeT), sizeof(facetT), 0, 0);
+
+struct qhT {
+
+/*----------------------------------
+
+  qh constants
+    configuration flags and constants for Qhull
+
+  notes:
+    The user configures Qhull by defining flags.  They are
+    copied into qh by qh_setflags().  qh-quick_r.htm#options defines the flags.
+*/
+  boolT ALLpoints;        /* true 'Qs' if search all points for initial simplex */
+  boolT ANGLEmerge;       /* true 'Qa' if sort potential merges by angle */
+  boolT APPROXhull;       /* true 'Wn' if MINoutside set */
+  realT   MINoutside;     /*   'Wn' min. distance for an outside point */
+  boolT ANNOTATEoutput;   /* true 'Ta' if annotate output with message codes */
+  boolT ATinfinity;       /* true 'Qz' if point num_points-1 is "at-infinity"
+                             for improving precision in Delaunay triangulations */
+  boolT AVOIDold;         /* true 'Q4' if avoid old->new merges */
+  boolT BESToutside;      /* true 'Qf' if partition points into best outsideset */
+  boolT CDDinput;         /* true 'Pc' if input uses CDD format (1.0/offset first) */
+  boolT CDDoutput;        /* true 'PC' if print normals in CDD format (offset first) */
+  boolT CHECKfrequently;  /* true 'Tc' if checking frequently */
+  realT premerge_cos;     /*   'A-n'   cos_max when pre merging */
+  realT postmerge_cos;    /*   'An'    cos_max when post merging */
+  boolT DELAUNAY;         /* true 'd' if computing DELAUNAY triangulation */
+  boolT DOintersections;  /* true 'Gh' if print hyperplane intersections */
+  int   DROPdim;          /* drops dim 'GDn' for 4-d -> 3-d output */
+  boolT FORCEoutput;      /* true 'Po' if forcing output despite degeneracies */
+  int   GOODpoint;        /* 1+n for 'QGn', good facet if visible/not(-) from point n*/
+  pointT *GOODpointp;     /*   the actual point */
+  boolT GOODthreshold;    /* true if qh.lower_threshold/upper_threshold defined
+                             false if qh.SPLITthreshold */
+  int   GOODvertex;       /* 1+n, good facet if vertex for point n */
+  pointT *GOODvertexp;     /*   the actual point */
+  boolT HALFspace;        /* true 'Hn,n,n' if halfspace intersection */
+  boolT ISqhullQh;        /* Set by Qhull.cpp on initialization */
+  int   IStracing;        /* trace execution, 0=none, 1=least, 4=most, -1=events */
+  int   KEEParea;         /* 'PAn' number of largest facets to keep */
+  boolT KEEPcoplanar;     /* true 'Qc' if keeping nearest facet for coplanar points */
+  boolT KEEPinside;       /* true 'Qi' if keeping nearest facet for inside points
+                              set automatically if 'd Qc' */
+  int   KEEPmerge;        /* 'PMn' number of facets to keep with most merges */
+  realT KEEPminArea;      /* 'PFn' minimum facet area to keep */
+  realT MAXcoplanar;      /* 'Un' max distance below a facet to be coplanar*/
+  boolT MERGEexact;       /* true 'Qx' if exact merges (coplanar, degen, dupridge, flipped) */
+  boolT MERGEindependent; /* true 'Q2' if merging independent sets */
+  boolT MERGING;          /* true if exact-, pre- or post-merging, with angle and centrum tests */
+  realT   premerge_centrum;  /*   'C-n' centrum_radius when pre merging.  Default is round-off */
+  realT   postmerge_centrum; /*   'Cn' centrum_radius when post merging.  Default is round-off */
+  boolT MERGEvertices;    /* true 'Q3' if merging redundant vertices */
+  realT MINvisible;       /* 'Vn' min. distance for a facet to be visible */
+  boolT NOnarrow;         /* true 'Q10' if no special processing for narrow distributions */
+  boolT NOnearinside;     /* true 'Q8' if ignore near-inside points when partitioning */
+  boolT NOpremerge;       /* true 'Q0' if no defaults for C-0 or Qx */
+  boolT NOwide;           /* true 'Q12' if no error on wide merge due to duplicate ridge */
+  boolT ONLYgood;         /* true 'Qg' if process points with good visible or horizon facets */
+  boolT ONLYmax;          /* true 'Qm' if only process points that increase max_outside */
+  boolT PICKfurthest;     /* true 'Q9' if process furthest of furthest points*/
+  boolT POSTmerge;        /* true if merging after buildhull (Cn or An) */
+  boolT PREmerge;         /* true if merging during buildhull (C-n or A-n) */
+                        /* NOTE: some of these names are similar to qh_PRINT names */
+  boolT PRINTcentrums;    /* true 'Gc' if printing centrums */
+  boolT PRINTcoplanar;    /* true 'Gp' if printing coplanar points */
+  int   PRINTdim;         /* print dimension for Geomview output */
+  boolT PRINTdots;        /* true 'Ga' if printing all points as dots */
+  boolT PRINTgood;        /* true 'Pg' if printing good facets */
+  boolT PRINTinner;       /* true 'Gi' if printing inner planes */
+  boolT PRINTneighbors;   /* true 'PG' if printing neighbors of good facets */
+  boolT PRINTnoplanes;    /* true 'Gn' if printing no planes */
+  boolT PRINToptions1st;  /* true 'FO' if printing options to stderr */
+  boolT PRINTouter;       /* true 'Go' if printing outer planes */
+  boolT PRINTprecision;   /* false 'Pp' if not reporting precision problems */
+  qh_PRINT PRINTout[qh_PRINTEND]; /* list of output formats to print */
+  boolT PRINTridges;      /* true 'Gr' if print ridges */
+  boolT PRINTspheres;     /* true 'Gv' if print vertices as spheres */
+  boolT PRINTstatistics;  /* true 'Ts' if printing statistics to stderr */
+  boolT PRINTsummary;     /* true 's' if printing summary to stderr */
+  boolT PRINTtransparent; /* true 'Gt' if print transparent outer ridges */
+  boolT PROJECTdelaunay;  /* true if DELAUNAY, no readpoints() and
+                             need projectinput() for Delaunay in qh_init_B */
+  int   PROJECTinput;     /* number of projected dimensions 'bn:0Bn:0' */
+  boolT QUICKhelp;        /* true if quick help message for degen input */
+  boolT RANDOMdist;       /* true if randomly change distplane and setfacetplane */
+  realT RANDOMfactor;     /*    maximum random perturbation */
+  realT RANDOMa;          /*    qh_randomfactor is randr * RANDOMa + RANDOMb */
+  realT RANDOMb;
+  boolT RANDOMoutside;    /* true if select a random outside point */
+  int   REPORTfreq;       /* buildtracing reports every n facets */
+  int   REPORTfreq2;      /* tracemerging reports every REPORTfreq/2 facets */
+  int   RERUN;            /* 'TRn' rerun qhull n times (qh.build_cnt) */
+  int   ROTATErandom;     /* 'QRn' seed, 0 time, >= rotate input */
+  boolT SCALEinput;       /* true 'Qbk' if scaling input */
+  boolT SCALElast;        /* true 'Qbb' if scale last coord to max prev coord */
+  boolT SETroundoff;      /* true 'E' if qh.DISTround is predefined */
+  boolT SKIPcheckmax;     /* true 'Q5' if skip qh_check_maxout */
+  boolT SKIPconvex;       /* true 'Q6' if skip convexity testing during pre-merge */
+  boolT SPLITthresholds;  /* true if upper_/lower_threshold defines a region
+                               used only for printing (!for qh.ONLYgood) */
+  int   STOPcone;         /* 'TCn' 1+n for stopping after cone for point n */
+                          /*       also used by qh_build_withresart for err exit*/
+  int   STOPpoint;        /* 'TVn' 'TV-n' 1+n for stopping after/before(-)
+                                        adding point n */
+  int   TESTpoints;       /* 'QTn' num of test points after qh.num_points.  Test points always coplanar. */
+  boolT TESTvneighbors;   /*  true 'Qv' if test vertex neighbors at end */
+  int   TRACElevel;       /* 'Tn' conditional IStracing level */
+  int   TRACElastrun;     /*  qh.TRACElevel applies to last qh.RERUN */
+  int   TRACEpoint;       /* 'TPn' start tracing when point n is a vertex */
+  realT TRACEdist;        /* 'TWn' start tracing when merge distance too big */
+  int   TRACEmerge;       /* 'TMn' start tracing before this merge */
+  boolT TRIangulate;      /* true 'Qt' if triangulate non-simplicial facets */
+  boolT TRInormals;       /* true 'Q11' if triangulate duplicates ->normal and ->center (sets Qt) */
+  boolT UPPERdelaunay;    /* true 'Qu' if computing furthest-site Delaunay */
+  boolT USEstdout;        /* true 'Tz' if using stdout instead of stderr */
+  boolT VERIFYoutput;     /* true 'Tv' if verify output at end of qhull */
+  boolT VIRTUALmemory;    /* true 'Q7' if depth-first processing in buildhull */
+  boolT VORONOI;          /* true 'v' if computing Voronoi diagram */
+
+  /*--------input constants ---------*/
+  realT AREAfactor;       /* 1/(hull_dim-1)! for converting det's to area */
+  boolT DOcheckmax;       /* true if calling qh_check_maxout (qhT *qh, qh_initqhull_globals) */
+  char  *feasible_string;  /* feasible point 'Hn,n,n' for halfspace intersection */
+  coordT *feasible_point;  /*    as coordinates, both malloc'd */
+  boolT GETarea;          /* true 'Fa', 'FA', 'FS', 'PAn', 'PFn' if compute facet area/Voronoi volume in io_r.c */
+  boolT KEEPnearinside;   /* true if near-inside points in coplanarset */
+  int   hull_dim;         /* dimension of hull, set by initbuffers */
+  int   input_dim;        /* dimension of input, set by initbuffers */
+  int   num_points;       /* number of input points */
+  pointT *first_point;    /* array of input points, see POINTSmalloc */
+  boolT POINTSmalloc;     /*   true if qh.first_point/num_points allocated */
+  pointT *input_points;   /* copy of original qh.first_point for input points for qh_joggleinput */
+  boolT input_malloc;     /* true if qh.input_points malloc'd */
+  char  qhull_command[256];/* command line that invoked this program */
+  int   qhull_commandsiz2; /*    size of qhull_command at qh_clear_outputflags */
+  char  rbox_command[256]; /* command line that produced the input points */
+  char  qhull_options[512];/* descriptive list of options */
+  int   qhull_optionlen;  /*    length of last line */
+  int   qhull_optionsiz;  /*    size of qhull_options at qh_build_withrestart */
+  int   qhull_optionsiz2; /*    size of qhull_options at qh_clear_outputflags */
+  int   run_id;           /* non-zero, random identifier for this instance of qhull */
+  boolT VERTEXneighbors;  /* true if maintaining vertex neighbors */
+  boolT ZEROcentrum;      /* true if 'C-0' or 'C-0 Qx'.  sets ZEROall_ok */
+  realT *upper_threshold; /* don't print if facet->normal[k]>=upper_threshold[k]
+                             must set either GOODthreshold or SPLITthreshold
+                             if Delaunay, default is 0.0 for upper envelope */
+  realT *lower_threshold; /* don't print if facet->normal[k] <=lower_threshold[k] */
+  realT *upper_bound;     /* scale point[k] to new upper bound */
+  realT *lower_bound;     /* scale point[k] to new lower bound
+                             project if both upper_ and lower_bound == 0 */
+
+/*----------------------------------
+
+  qh precision constants
+    precision constants for Qhull
+
+  notes:
+    qh_detroundoff(qh) computes the maximum roundoff error for distance
+    and other computations.  It also sets default values for the
+    qh constants above.
+*/
+  realT ANGLEround;       /* max round off error for angles */
+  realT centrum_radius;   /* max centrum radius for convexity (roundoff added) */
+  realT cos_max;          /* max cosine for convexity (roundoff added) */
+  realT DISTround;        /* max round off error for distances, 'E' overrides qh_distround() */
+  realT MAXabs_coord;     /* max absolute coordinate */
+  realT MAXlastcoord;     /* max last coordinate for qh_scalelast */
+  realT MAXsumcoord;      /* max sum of coordinates */
+  realT MAXwidth;         /* max rectilinear width of point coordinates */
+  realT MINdenom_1;       /* min. abs. value for 1/x */
+  realT MINdenom;         /*    use divzero if denominator < MINdenom */
+  realT MINdenom_1_2;     /* min. abs. val for 1/x that allows normalization */
+  realT MINdenom_2;       /*    use divzero if denominator < MINdenom_2 */
+  realT MINlastcoord;     /* min. last coordinate for qh_scalelast */
+  boolT NARROWhull;       /* set in qh_initialhull if angle < qh_MAXnarrow */
+  realT *NEARzero;        /* hull_dim array for near zero in gausselim */
+  realT NEARinside;       /* keep points for qh_check_maxout if close to facet */
+  realT ONEmerge;         /* max distance for merging simplicial facets */
+  realT outside_err;      /* application's epsilon for coplanar points
+                             qh_check_bestdist() qh_check_points() reports error if point outside */
+  realT WIDEfacet;        /* size of wide facet for skipping ridge in
+                             area computation and locking centrum */
+
+/*----------------------------------
+
+  qh internal constants
+    internal constants for Qhull
+*/
+  char qhull[sizeof("qhull")]; /* "qhull" for checking ownership while debugging */
+  jmp_buf errexit;        /* exit label for qh_errexit, defined by setjmp() and NOerrexit */
+  char jmpXtra[40];       /* extra bytes in case jmp_buf is defined wrong by compiler */
+  jmp_buf restartexit;    /* restart label for qh_errexit, defined by setjmp() and ALLOWrestart */
+  char jmpXtra2[40];      /* extra bytes in case jmp_buf is defined wrong by compiler*/
+  FILE *fin;              /* pointer to input file, init by qh_initqhull_start */
+  FILE *fout;             /* pointer to output file */
+  FILE *ferr;             /* pointer to error file */
+  pointT *interior_point; /* center point of the initial simplex*/
+  int normal_size;     /* size in bytes for facet normals and point coords*/
+  int center_size;     /* size in bytes for Voronoi centers */
+  int   TEMPsize;         /* size for small, temporary sets (in quick mem) */
+
+/*----------------------------------
+
+  qh facet and vertex lists
+    defines lists of facets, new facets, visible facets, vertices, and
+    new vertices.  Includes counts, next ids, and trace ids.
+  see:
+    qh_resetlists()
+*/
+  facetT *facet_list;     /* first facet */
+  facetT  *facet_tail;     /* end of facet_list (dummy facet) */
+  facetT *facet_next;     /* next facet for buildhull()
+                             previous facets do not have outside sets
+                             NARROWhull: previous facets may have coplanar outside sets for qh_outcoplanar */
+  facetT *newfacet_list;  /* list of new facets to end of facet_list */
+  facetT *visible_list;   /* list of visible facets preceding newfacet_list,
+                             facet->visible set */
+  int       num_visible;  /* current number of visible facets */
+  unsigned tracefacet_id;  /* set at init, then can print whenever */
+  facetT *tracefacet;     /*   set in newfacet/mergefacet, undone in delfacet*/
+  unsigned tracevertex_id;  /* set at buildtracing, can print whenever */
+  vertexT *tracevertex;     /*   set in newvertex, undone in delvertex*/
+  vertexT *vertex_list;     /* list of all vertices, to vertex_tail */
+  vertexT  *vertex_tail;    /*      end of vertex_list (dummy vertex) */
+  vertexT *newvertex_list; /* list of vertices in newfacet_list, to vertex_tail
+                             all vertices have 'newlist' set */
+  int   num_facets;       /* number of facets in facet_list
+                             includes visible faces (num_visible) */
+  int   num_vertices;     /* number of vertices in facet_list */
+  int   num_outside;      /* number of points in outsidesets (for tracing and RANDOMoutside)
+                               includes coplanar outsideset points for NARROWhull/qh_outcoplanar() */
+  int   num_good;         /* number of good facets (after findgood_all) */
+  unsigned facet_id;      /* ID of next, new facet from newfacet() */
+  unsigned ridge_id;      /* ID of next, new ridge from newridge() */
+  unsigned vertex_id;     /* ID of next, new vertex from newvertex() */
+
+/*----------------------------------
+
+  qh global variables
+    defines minimum and maximum distances, next visit ids, several flags,
+    and other global variables.
+    initialize in qh_initbuild or qh_maxmin if used in qh_buildhull
+*/
+  unsigned long hulltime; /* ignore time to set up input and randomize */
+                          /*   use unsigned to avoid wrap-around errors */
+  boolT ALLOWrestart;     /* true if qh_precision can use qh.restartexit */
+  int   build_cnt;        /* number of calls to qh_initbuild */
+  qh_CENTER CENTERtype;   /* current type of facet->center, qh_CENTER */
+  int   furthest_id;      /* pointid of furthest point, for tracing */
+  facetT *GOODclosest;    /* closest facet to GOODthreshold in qh_findgood */
+  boolT hasAreaVolume;    /* true if totarea, totvol was defined by qh_getarea */
+  boolT hasTriangulation; /* true if triangulation created by qh_triangulate */
+  realT JOGGLEmax;        /* set 'QJn' if randomly joggle input */
+  boolT maxoutdone;       /* set qh_check_maxout(), cleared by qh_addpoint() */
+  realT max_outside;      /* maximum distance from a point to a facet,
+                               before roundoff, not simplicial vertices
+                               actual outer plane is +DISTround and
+                               computed outer plane is +2*DISTround */
+  realT max_vertex;       /* maximum distance (>0) from vertex to a facet,
+                               before roundoff, due to a merge */
+  realT min_vertex;       /* minimum distance (<0) from vertex to a facet,
+                               before roundoff, due to a merge
+                               if qh.JOGGLEmax, qh_makenewplanes sets it
+                               recomputed if qh.DOcheckmax, default -qh.DISTround */
+  boolT NEWfacets;        /* true while visible facets invalid due to new or merge
+                              from makecone/attachnewfacets to deletevisible */
+  boolT findbestnew;      /* true if partitioning calls qh_findbestnew */
+  boolT findbest_notsharp; /* true if new facets are at least 90 degrees */
+  boolT NOerrexit;        /* true if qh.errexit is not available, cleared after setjmp */
+  realT PRINTcradius;     /* radius for printing centrums */
+  realT PRINTradius;      /* radius for printing vertex spheres and points */
+  boolT POSTmerging;      /* true when post merging */
+  int   printoutvar;      /* temporary variable for qh_printbegin, etc. */
+  int   printoutnum;      /* number of facets printed */
+  boolT QHULLfinished;    /* True after qhull() is finished */
+  realT totarea;          /* 'FA': total facet area computed by qh_getarea, hasAreaVolume */
+  realT totvol;           /* 'FA': total volume computed by qh_getarea, hasAreaVolume */
+  unsigned int visit_id;  /* unique ID for searching neighborhoods, */
+  unsigned int vertex_visit; /* unique ID for searching vertices, reset with qh_buildtracing */
+  boolT ZEROall_ok;       /* True if qh_checkzero always succeeds */
+  boolT WAScoplanar;      /* True if qh_partitioncoplanar (qhT *qh, qh_check_maxout) */
+
+/*----------------------------------
+
+  qh global sets
+    defines sets for merging, initial simplex, hashing, extra input points,
+    and deleted vertices
+*/
+  setT *facet_mergeset;   /* temporary set of merges to be done */
+  setT *degen_mergeset;   /* temporary set of degenerate and redundant merges */
+  setT *hash_table;       /* hash table for matching ridges in qh_matchfacets
+                             size is setsize() */
+  setT *other_points;     /* additional points */
+  setT *del_vertices;     /* vertices to partition and delete with visible
+                             facets.  Have deleted set for checkfacet */
+
+/*----------------------------------
+
+  qh global buffers
+    defines buffers for maxtrix operations, input, and error messages
+*/
+  coordT *gm_matrix;      /* (dim+1)Xdim matrix for geom_r.c */
+  coordT **gm_row;        /* array of gm_matrix rows */
+  char* line;             /* malloc'd input line of maxline+1 chars */
+  int maxline;
+  coordT *half_space;     /* malloc'd input array for halfspace (qh.normal_size+coordT) */
+  coordT *temp_malloc;    /* malloc'd input array for points */
+
+/*----------------------------------
+
+  qh static variables
+    defines static variables for individual functions
+
+  notes:
+    do not use 'static' within a function.  Multiple instances of qhull
+    may exist.
+
+    do not assume zero initialization, 'QPn' may cause a restart
+*/
+  boolT ERREXITcalled;    /* true during qh_errexit (qhT *qh, prevents duplicate calls */
+  boolT firstcentrum;     /* for qh_printcentrum */
+  boolT old_randomdist;   /* save RANDOMdist flag during io, tracing, or statistics */
+  setT *coplanarfacetset;  /* set of coplanar facets for searching qh_findbesthorizon() */
+  realT last_low;         /* qh_scalelast parameters for qh_setdelaunay */
+  realT last_high;
+  realT last_newhigh;
+  unsigned lastreport;    /* for qh_buildtracing */
+  int mergereport;        /* for qh_tracemerging */
+  setT *old_tempstack;    /* for saving qh->qhmem.tempstack in save_qhull */
+  int   ridgeoutnum;      /* number of ridges for 4OFF output (qh_printbegin,etc) */
+
+/*----------------------------------
+
+  qh memory management, rbox globals, and statistics
+
+  Replaces global data structures defined for libqhull
+*/
+  int     last_random;    /* Last random number from qh_rand (random_r.c) */
+  jmp_buf rbox_errexit;   /* errexit from rboxlib_r.c, defined by qh_rboxpoints() only */
+  char    jmpXtra3[40];   /* extra bytes in case jmp_buf is defined wrong by compiler */
+  int     rbox_isinteger;
+  double  rbox_out_offset;
+  void *  cpp_object;     /* C++ pointer.  Currently used by RboxPoints.qh_fprintf_rbox */
+
+  /* Last, otherwise zero'd by qh_initqhull_start2 (global_r.c */
+  qhmemT  qhmem;          /* Qhull managed memory (mem_r.h) */
+  /* After qhmem because its size depends on the number of statistics */
+  qhstatT qhstat;         /* Qhull statistics (stat_r.h) */
+};
+
+/*=========== -macros- =========================*/
+
+/*----------------------------------
+
+  otherfacet_(ridge, facet)
+    return neighboring facet for a ridge in facet
+*/
+#define otherfacet_(ridge, facet) \
+                        (((ridge)->top == (facet)) ? (ridge)->bottom : (ridge)->top)
+
+/*----------------------------------
+
+  getid_(p)
+    return int ID for facet, ridge, or vertex
+    return qh_IDunknown(-1) if NULL
+*/
+#define getid_(p)       ((p) ? (int)((p)->id) : qh_IDunknown)
+
+/*============== FORALL macros ===================*/
+
+/*----------------------------------
+
+  FORALLfacets { ... }
+    assign 'facet' to each facet in qh.facet_list
+
+  notes:
+    uses 'facetT *facet;'
+    assumes last facet is a sentinel
+    assumes qh defined
+
+  see:
+    FORALLfacet_( facetlist )
+*/
+#define FORALLfacets for (facet=qh->facet_list;facet && facet->next;facet=facet->next)
+
+/*----------------------------------
+
+  FORALLpoints { ... }
+    assign 'point' to each point in qh.first_point, qh.num_points
+
+  notes:
+    assumes qh defined
+
+  declare:
+    coordT *point, *pointtemp;
+*/
+#define FORALLpoints FORALLpoint_(qh, qh->first_point, qh->num_points)
+
+/*----------------------------------
+
+  FORALLpoint_( qh, points, num) { ... }
+    assign 'point' to each point in points array of num points
+
+  declare:
+    coordT *point, *pointtemp;
+*/
+#define FORALLpoint_(qh, points, num) for (point= (points), \
+      pointtemp= (points)+qh->hull_dim*(num); point < pointtemp; point += qh->hull_dim)
+
+/*----------------------------------
+
+  FORALLvertices { ... }
+    assign 'vertex' to each vertex in qh.vertex_list
+
+  declare:
+    vertexT *vertex;
+
+  notes:
+    assumes qh.vertex_list terminated with a sentinel
+    assumes qh defined
+*/
+#define FORALLvertices for (vertex=qh->vertex_list;vertex && vertex->next;vertex= vertex->next)
+
+/*----------------------------------
+
+  FOREACHfacet_( facets ) { ... }
+    assign 'facet' to each facet in facets
+
+  declare:
+    facetT *facet, **facetp;
+
+  see:
+    FOREACHsetelement_
+*/
+#define FOREACHfacet_(facets)    FOREACHsetelement_(facetT, facets, facet)
+
+/*----------------------------------
+
+  FOREACHneighbor_( facet ) { ... }
+    assign 'neighbor' to each neighbor in facet->neighbors
+
+  FOREACHneighbor_( vertex ) { ... }
+    assign 'neighbor' to each neighbor in vertex->neighbors
+
+  declare:
+    facetT *neighbor, **neighborp;
+
+  see:
+    FOREACHsetelement_
+*/
+#define FOREACHneighbor_(facet)  FOREACHsetelement_(facetT, facet->neighbors, neighbor)
+
+/*----------------------------------
+
+  FOREACHpoint_( points ) { ... }
+    assign 'point' to each point in points set
+
+  declare:
+    pointT *point, **pointp;
+
+  see:
+    FOREACHsetelement_
+*/
+#define FOREACHpoint_(points)    FOREACHsetelement_(pointT, points, point)
+
+/*----------------------------------
+
+  FOREACHridge_( ridges ) { ... }
+    assign 'ridge' to each ridge in ridges set
+
+  declare:
+    ridgeT *ridge, **ridgep;
+
+  see:
+    FOREACHsetelement_
+*/
+#define FOREACHridge_(ridges)    FOREACHsetelement_(ridgeT, ridges, ridge)
+
+/*----------------------------------
+
+  FOREACHvertex_( vertices ) { ... }
+    assign 'vertex' to each vertex in vertices set
+
+  declare:
+    vertexT *vertex, **vertexp;
+
+  see:
+    FOREACHsetelement_
+*/
+#define FOREACHvertex_(vertices) FOREACHsetelement_(vertexT, vertices,vertex)
+
+/*----------------------------------
+
+  FOREACHfacet_i_( qh, facets ) { ... }
+    assign 'facet' and 'facet_i' for each facet in facets set
+
+  declare:
+    facetT *facet;
+    int     facet_n, facet_i;
+
+  see:
+    FOREACHsetelement_i_
+*/
+#define FOREACHfacet_i_(qh, facets)    FOREACHsetelement_i_(qh, facetT, facets, facet)
+
+/*----------------------------------
+
+  FOREACHneighbor_i_( qh, facet ) { ... }
+    assign 'neighbor' and 'neighbor_i' for each neighbor in facet->neighbors
+
+  FOREACHneighbor_i_( qh, vertex ) { ... }
+    assign 'neighbor' and 'neighbor_i' for each neighbor in vertex->neighbors
+
+  declare:
+    facetT *neighbor;
+    int     neighbor_n, neighbor_i;
+
+  see:
+    FOREACHsetelement_i_
+*/
+#define FOREACHneighbor_i_(qh, facet)  FOREACHsetelement_i_(qh, facetT, facet->neighbors, neighbor)
+
+/*----------------------------------
+
+  FOREACHpoint_i_( qh, points ) { ... }
+    assign 'point' and 'point_i' for each point in points set
+
+  declare:
+    pointT *point;
+    int     point_n, point_i;
+
+  see:
+    FOREACHsetelement_i_
+*/
+#define FOREACHpoint_i_(qh, points)    FOREACHsetelement_i_(qh, pointT, points, point)
+
+/*----------------------------------
+
+  FOREACHridge_i_( qh, ridges ) { ... }
+    assign 'ridge' and 'ridge_i' for each ridge in ridges set
+
+  declare:
+    ridgeT *ridge;
+    int     ridge_n, ridge_i;
+
+  see:
+    FOREACHsetelement_i_
+*/
+#define FOREACHridge_i_(qh, ridges)    FOREACHsetelement_i_(qh, ridgeT, ridges, ridge)
+
+/*----------------------------------
+
+  FOREACHvertex_i_( qh, vertices ) { ... }
+    assign 'vertex' and 'vertex_i' for each vertex in vertices set
+
+  declare:
+    vertexT *vertex;
+    int     vertex_n, vertex_i;
+
+  see:
+    FOREACHsetelement_i_
+*/
+#define FOREACHvertex_i_(qh, vertices) FOREACHsetelement_i_(qh, vertexT, vertices,vertex)
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/********* -libqhull_r.c prototypes (duplicated from qhull_ra.h) **********************/
+
+void    qh_qhull(qhT *qh);
+boolT   qh_addpoint(qhT *qh, pointT *furthest, facetT *facet, boolT checkdist);
+void    qh_printsummary(qhT *qh, FILE *fp);
+
+/********* -user.c prototypes (alphabetical) **********************/
+
+void    qh_errexit(qhT *qh, int exitcode, facetT *facet, ridgeT *ridge);
+void    qh_errprint(qhT *qh, const char* string, facetT *atfacet, facetT *otherfacet, ridgeT *atridge, vertexT *atvertex);
+int     qh_new_qhull(qhT *qh, int dim, int numpoints, coordT *points, boolT ismalloc,
+                char *qhull_cmd, FILE *outfile, FILE *errfile);
+void    qh_printfacetlist(qhT *qh, facetT *facetlist, setT *facets, boolT printall);
+void    qh_printhelp_degenerate(qhT *qh, FILE *fp);
+void    qh_printhelp_narrowhull(qhT *qh, FILE *fp, realT minangle);
+void    qh_printhelp_singular(qhT *qh, FILE *fp);
+void    qh_user_memsizes(qhT *qh);
+
+/********* -usermem_r.c prototypes (alphabetical) **********************/
+void    qh_exit(int exitcode);
+void    qh_fprintf_stderr(int msgcode, const char *fmt, ... );
+void    qh_free(void *mem);
+void   *qh_malloc(size_t size);
+
+/********* -userprintf_r.c and userprintf_rbox_r.c prototypes **********************/
+void    qh_fprintf(qhT *qh, FILE *fp, int msgcode, const char *fmt, ... );
+void    qh_fprintf_rbox(qhT *qh, FILE *fp, int msgcode, const char *fmt, ... );
+
+/***** -geom_r.c/geom2_r.c/random_r.c prototypes (duplicated from geom_r.h, random_r.h) ****************/
+
+facetT *qh_findbest(qhT *qh, pointT *point, facetT *startfacet,
+                     boolT bestoutside, boolT newfacets, boolT noupper,
+                     realT *dist, boolT *isoutside, int *numpart);
+facetT *qh_findbestnew(qhT *qh, pointT *point, facetT *startfacet,
+                     realT *dist, boolT bestoutside, boolT *isoutside, int *numpart);
+boolT   qh_gram_schmidt(qhT *qh, int dim, realT **rows);
+void    qh_outerinner(qhT *qh, facetT *facet, realT *outerplane, realT *innerplane);
+void    qh_printsummary(qhT *qh, FILE *fp);
+void    qh_projectinput(qhT *qh);
+void    qh_randommatrix(qhT *qh, realT *buffer, int dim, realT **row);
+void    qh_rotateinput(qhT *qh, realT **rows);
+void    qh_scaleinput(qhT *qh);
+void    qh_setdelaunay(qhT *qh, int dim, int count, pointT *points);
+coordT  *qh_sethalfspace_all(qhT *qh, int dim, int count, coordT *halfspaces, pointT *feasible);
+
+/***** -global_r.c prototypes (alphabetical) ***********************/
+
+unsigned long qh_clock(qhT *qh);
+void    qh_checkflags(qhT *qh, char *command, char *hiddenflags);
+void    qh_clear_outputflags(qhT *qh);
+void    qh_freebuffers(qhT *qh);
+void    qh_freeqhull(qhT *qh, boolT allmem);
+void    qh_init_A(qhT *qh, FILE *infile, FILE *outfile, FILE *errfile, int argc, char *argv[]);
+void    qh_init_B(qhT *qh, coordT *points, int numpoints, int dim, boolT ismalloc);
+void    qh_init_qhull_command(qhT *qh, int argc, char *argv[]);
+void    qh_initbuffers(qhT *qh, coordT *points, int numpoints, int dim, boolT ismalloc);
+void    qh_initflags(qhT *qh, char *command);
+void    qh_initqhull_buffers(qhT *qh);
+void    qh_initqhull_globals(qhT *qh, coordT *points, int numpoints, int dim, boolT ismalloc);
+void    qh_initqhull_mem(qhT *qh);
+void    qh_initqhull_outputflags(qhT *qh);
+void    qh_initqhull_start(qhT *qh, FILE *infile, FILE *outfile, FILE *errfile);
+void    qh_initqhull_start2(qhT *qh, FILE *infile, FILE *outfile, FILE *errfile);
+void    qh_initthresholds(qhT *qh, char *command);
+void    qh_lib_check(int qhullLibraryType, int qhTsize, int vertexTsize, int ridgeTsize, int facetTsize, int setTsize, int qhmemTsize);
+void    qh_option(qhT *qh, const char *option, int *i, realT *r);
+void    qh_zero(qhT *qh, FILE *errfile);
+
+/***** -io_r.c prototypes (duplicated from io_r.h) ***********************/
+
+void    qh_dfacet(qhT *qh, unsigned id);
+void    qh_dvertex(qhT *qh, unsigned id);
+void    qh_printneighborhood(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetA, facetT *facetB, boolT printall);
+void    qh_produce_output(qhT *qh);
+coordT *qh_readpoints(qhT *qh, int *numpoints, int *dimension, boolT *ismalloc);
+
+
+/********* -mem_r.c prototypes (duplicated from mem_r.h) **********************/
+
+void qh_meminit(qhT *qh, FILE *ferr);
+void qh_memfreeshort(qhT *qh, int *curlong, int *totlong);
+
+/********* -poly_r.c/poly2_r.c prototypes (duplicated from poly_r.h) **********************/
+
+void    qh_check_output(qhT *qh);
+void    qh_check_points(qhT *qh);
+setT   *qh_facetvertices(qhT *qh, facetT *facetlist, setT *facets, boolT allfacets);
+facetT *qh_findbestfacet(qhT *qh, pointT *point, boolT bestoutside,
+           realT *bestdist, boolT *isoutside);
+vertexT *qh_nearvertex(qhT *qh, facetT *facet, pointT *point, realT *bestdistp);
+pointT *qh_point(qhT *qh, int id);
+setT   *qh_pointfacet(qhT *qh /*qh.facet_list*/);
+int     qh_pointid(qhT *qh, pointT *point);
+setT   *qh_pointvertex(qhT *qh /*qh.facet_list*/);
+void    qh_setvoronoi_all(qhT *qh);
+void    qh_triangulate(qhT *qh /*qh.facet_list*/);
+
+/********* -rboxpoints_r.c prototypes **********************/
+int     qh_rboxpoints(qhT *qh, char* rbox_command);
+void    qh_errexit_rbox(qhT *qh, int exitcode);
+
+/********* -stat_r.c prototypes (duplicated from stat_r.h) **********************/
+
+void    qh_collectstatistics(qhT *qh);
+void    qh_printallstatistics(qhT *qh, FILE *fp, const char *string);
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /* qhDEFlibqhull */
diff --git a/Cdir/mem_r.h b/Cdir/mem_r.h
new file mode 100644
--- /dev/null
+++ b/Cdir/mem_r.h
@@ -0,0 +1,234 @@
+/*
  ---------------------------------
+
+   mem_r.h
+     prototypes for memory management functions
+
+   see qh-mem_r.htm, mem_r.c and qset_r.h
+
+   for error handling, writes message and calls
+     qh_errexit(qhT *qh, qhmem_ERRmem, NULL, NULL) if insufficient memory
+       and
+     qh_errexit(qhT *qh, qhmem_ERRqhull, NULL, NULL) otherwise
+
+   Copyright (c) 1993-2015 The Geometry Center.
+   $Id: //main/2015/qhull/src/libqhull_r/mem_r.h#4 $$Change: 2079 $
+   $DateTime: 2016/02/07 17:43:34 $$Author: bbarber $
+*/
+
+#ifndef qhDEFmem
+#define qhDEFmem 1
+
+#include 
+
+#ifndef DEFsetT
+#define DEFsetT 1
+typedef struct setT setT;          /* defined in qset_r.h */
+#endif
+
+#ifndef DEFqhT
+#define DEFqhT 1
+typedef struct qhT qhT;          /* defined in libqhull_r.h */
+#endif
+
+/*---------------------------------
+
+  qh_NOmem
+    turn off quick-fit memory allocation
+
+  notes:
+    mem_r.c implements Quickfit memory allocation for about 20% time
+    savings.  If it fails on your machine, try to locate the
+    problem, and send the answer to qhull@qhull.org.  If this can
+    not be done, define qh_NOmem to use malloc/free instead.
+
+   #define qh_NOmem
+*/
+
+/*---------------------------------
+
+qh_TRACEshort
+Trace short and quick memory allocations at T5
+
+*/
+#define qh_TRACEshort
+
+/*-------------------------------------------
+    to avoid bus errors, memory allocation must consider alignment requirements.
+    malloc() automatically takes care of alignment.   Since mem_r.c manages
+    its own memory, we need to explicitly specify alignment in
+    qh_meminitbuffers().
+
+    A safe choice is sizeof(double).  sizeof(float) may be used if doubles
+    do not occur in data structures and pointers are the same size.  Be careful
+    of machines (e.g., DEC Alpha) with large pointers.  If gcc is available,
+    use __alignof__(double) or fmax_(__alignof__(float), __alignof__(void *)).
+
+   see qh_MEMalign in user.h for qhull's alignment
+*/
+
+#define qhmem_ERRmem 4    /* matches qh_ERRmem in libqhull_r.h */
+#define qhmem_ERRqhull 5  /* matches qh_ERRqhull in libqhull_r.h */
+
+/*----------------------------------
+
+  ptr_intT
+    for casting a void * to an integer-type that holds a pointer
+    Used for integer expressions (e.g., computing qh_gethash() in poly_r.c)
+
+  notes:
+    WARN64 -- these notes indicate 64-bit issues
+    On 64-bit machines, a pointer may be larger than an 'int'.
+    qh_meminit()/mem_r.c checks that 'ptr_intT' holds a 'void*'
+    ptr_intT is typically a signed value, but not necessarily so
+    size_t is typically unsigned, but should match the parameter type
+    Qhull uses int instead of size_t except for system calls such as malloc, qsort, qh_malloc, etc.
+    This matches Qt convention and is easier to work with.
+*/
+#if (defined(__MINGW64__)) && defined(_WIN64)
+typedef long long ptr_intT;
+#elif (_MSC_VER) && defined(_WIN64)
+typedef long long ptr_intT;
+#else
+typedef long ptr_intT;
+#endif
+
+/*----------------------------------
+
+  qhmemT
+    global memory structure for mem_r.c
+
+ notes:
+   users should ignore qhmem except for writing extensions
+   qhmem is allocated in mem_r.c
+
+   qhmem could be swapable like qh and qhstat, but then
+   multiple qh's and qhmem's would need to keep in synch.
+   A swapable qhmem would also waste memory buffers.  As long
+   as memory operations are atomic, there is no problem with
+   multiple qh structures being active at the same time.
+   If you need separate address spaces, you can swap the
+   contents of qh->qhmem.
+*/
+typedef struct qhmemT qhmemT;
+
+/* Update qhmem in mem_r.c if add or remove fields */
+struct qhmemT {               /* global memory management variables */
+  int      BUFsize;           /* size of memory allocation buffer */
+  int      BUFinit;           /* initial size of memory allocation buffer */
+  int      TABLEsize;         /* actual number of sizes in free list table */
+  int      NUMsizes;          /* maximum number of sizes in free list table */
+  int      LASTsize;          /* last size in free list table */
+  int      ALIGNmask;         /* worst-case alignment, must be 2^n-1 */
+  void   **freelists;          /* free list table, linked by offset 0 */
+  int     *sizetable;         /* size of each freelist */
+  int     *indextable;        /* size->index table */
+  void    *curbuffer;         /* current buffer, linked by offset 0 */
+  void    *freemem;           /*   free memory in curbuffer */
+  int      freesize;          /*   size of freemem in bytes */
+  setT    *tempstack;         /* stack of temporary memory, managed by users */
+  FILE    *ferr;              /* file for reporting errors when 'qh' may be undefined */
+  int      IStracing;         /* =5 if tracing memory allocations */
+  int      cntquick;          /* count of quick allocations */
+                              /* Note: removing statistics doesn't effect speed */
+  int      cntshort;          /* count of short allocations */
+  int      cntlong;           /* count of long allocations */
+  int      freeshort;         /* count of short memfrees */
+  int      freelong;          /* count of long memfrees */
+  int      totbuffer;         /* total short memory buffers minus buffer links */
+  int      totdropped;        /* total dropped memory at end of short memory buffers (e.g., freesize) */
+  int      totfree;           /* total size of free, short memory on freelists */
+  int      totlong;           /* total size of long memory in use */
+  int      maxlong;           /*   maximum totlong */
+  int      totshort;          /* total size of short memory in use */
+  int      totunused;         /* total unused short memory (estimated, short size - request size of first allocations) */
+  int      cntlarger;         /* count of setlarger's */
+  int      totlarger;         /* total copied by setlarger */
+};
+
+
+/*==================== -macros ====================*/
+
+/*----------------------------------
+
+  qh_memalloc_(qh, insize, freelistp, object, type)
+    returns object of size bytes
+        assumes size<=qh->qhmem.LASTsize and void **freelistp is a temp
+*/
+
+#if defined qh_NOmem
+#define qh_memalloc_(qh, insize, freelistp, object, type) {\
+  object= (type*)qh_memalloc(qh, insize); }
+#elif defined qh_TRACEshort
+#define qh_memalloc_(qh, insize, freelistp, object, type) {\
+    freelistp= NULL; /* Avoid warnings */ \
+    object= (type*)qh_memalloc(qh, insize); }
+#else /* !qh_NOmem */
+
+#define qh_memalloc_(qh, insize, freelistp, object, type) {\
+  freelistp= qh->qhmem.freelists + qh->qhmem.indextable[insize];\
+  if ((object= (type*)*freelistp)) {\
+    qh->qhmem.totshort += qh->qhmem.sizetable[qh->qhmem.indextable[insize]]; \
+    qh->qhmem.totfree -= qh->qhmem.sizetable[qh->qhmem.indextable[insize]]; \
+    qh->qhmem.cntquick++;  \
+    *freelistp= *((void **)*freelistp);\
+  }else object= (type*)qh_memalloc(qh, insize);}
+#endif
+
+/*----------------------------------
+
+  qh_memfree_(qh, object, insize, freelistp)
+    free up an object
+
+  notes:
+    object may be NULL
+    assumes size<=qh->qhmem.LASTsize and void **freelistp is a temp
+*/
+#if defined qh_NOmem
+#define qh_memfree_(qh, object, insize, freelistp) {\
+  qh_memfree(qh, object, insize); }
+#elif defined qh_TRACEshort
+#define qh_memfree_(qh, object, insize, freelistp) {\
+    freelistp= NULL; /* Avoid warnings */ \
+    qh_memfree(qh, object, insize); }
+#else /* !qh_NOmem */
+
+#define qh_memfree_(qh, object, insize, freelistp) {\
+  if (object) { \
+    qh->qhmem.freeshort++;\
+    freelistp= qh->qhmem.freelists + qh->qhmem.indextable[insize];\
+    qh->qhmem.totshort -= qh->qhmem.sizetable[qh->qhmem.indextable[insize]]; \
+    qh->qhmem.totfree += qh->qhmem.sizetable[qh->qhmem.indextable[insize]]; \
+    *((void **)object)= *freelistp;\
+    *freelistp= object;}}
+#endif
+
+/*=============== prototypes in alphabetical order ============*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void *qh_memalloc(qhT *qh, int insize);
+void qh_memcheck(qhT *qh);
+void qh_memfree(qhT *qh, void *object, int insize);
+void qh_memfreeshort(qhT *qh, int *curlong, int *totlong);
+void qh_meminit(qhT *qh, FILE *ferr);
+void qh_meminitbuffers(qhT *qh, int tracelevel, int alignment, int numsizes,
+                        int bufsize, int bufinit);
+void qh_memsetup(qhT *qh);
+void qh_memsize(qhT *qh, int size);
+void qh_memstatistics(qhT *qh, FILE *fp);
+void qh_memtotal(qhT *qh, int *totlong, int *curlong, int *totshort, int *curshort, int *maxlong, int *totbuffer);
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /* qhDEFmem */
diff --git a/Cdir/merge_r.h b/Cdir/merge_r.h
new file mode 100644
--- /dev/null
+++ b/Cdir/merge_r.h
@@ -0,0 +1,186 @@
+/*
  ---------------------------------
+
+   merge_r.h
+   header file for merge_r.c
+
+   see qh-merge_r.htm and merge_r.c
+
+   Copyright (c) 1993-2015 C.B. Barber.
+   $Id: //main/2015/qhull/src/libqhull_r/merge_r.h#3 $$Change: 2079 $
+   $DateTime: 2016/02/07 17:43:34 $$Author: bbarber $
+*/
+
+#ifndef qhDEFmerge
+#define qhDEFmerge 1
+
+#include "libqhull_r.h"
+
+
+/*============ -constants- ==============*/
+
+/*----------------------------------
+
+  qh_ANGLEredundant
+    indicates redundant merge in mergeT->angle
+*/
+#define qh_ANGLEredundant 6.0
+
+/*----------------------------------
+
+  qh_ANGLEdegen
+    indicates degenerate facet in mergeT->angle
+*/
+#define qh_ANGLEdegen     5.0
+
+/*----------------------------------
+
+  qh_ANGLEconcave
+    offset to indicate concave facets in mergeT->angle
+
+  notes:
+    concave facets are assigned the range of [2,4] in mergeT->angle
+    roundoff error may make the angle less than 2
+*/
+#define qh_ANGLEconcave  1.5
+
+/*----------------------------------
+
+  MRG... (mergeType)
+    indicates the type of a merge (mergeT->type)
+*/
+typedef enum {  /* in sort order for facet_mergeset */
+  MRGnone= 0,
+  MRGcoplanar,          /* centrum coplanar */
+  MRGanglecoplanar,     /* angle coplanar */
+                        /* could detect half concave ridges */
+  MRGconcave,           /* concave ridge */
+  MRGflip,              /* flipped facet. facet1 == facet2 */
+  MRGridge,             /* duplicate ridge (qh_MERGEridge) */
+                        /* degen and redundant go onto degen_mergeset */
+  MRGdegen,             /* degenerate facet (!enough neighbors) facet1 == facet2 */
+  MRGredundant,         /* redundant facet (vertex subset) */
+                        /* merge_degenredundant assumes degen < redundant */
+  MRGmirror,            /* mirror facet from qh_triangulate */
+  ENDmrg
+} mergeType;
+
+/*----------------------------------
+
+  qh_MERGEapex
+    flag for qh_mergefacet() to indicate an apex merge
+*/
+#define qh_MERGEapex     True
+
+/*============ -structures- ====================*/
+
+/*----------------------------------
+
+  mergeT
+    structure used to merge facets
+*/
+
+typedef struct mergeT mergeT;
+struct mergeT {         /* initialize in qh_appendmergeset */
+  realT   angle;        /* angle between normals of facet1 and facet2 */
+  facetT *facet1;       /* will merge facet1 into facet2 */
+  facetT *facet2;
+  mergeType type;
+};
+
+
+/*=========== -macros- =========================*/
+
+/*----------------------------------
+
+  FOREACHmerge_( merges ) {...}
+    assign 'merge' to each merge in merges
+
+  notes:
+    uses 'mergeT *merge, **mergep;'
+    if qh_mergefacet(),
+      restart since qh.facet_mergeset may change
+    see FOREACHsetelement_
+*/
+#define FOREACHmerge_( merges ) FOREACHsetelement_(mergeT, merges, merge)
+
+/*============ prototypes in alphabetical order after pre/postmerge =======*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void    qh_premerge(qhT *qh, vertexT *apex, realT maxcentrum, realT maxangle);
+void    qh_postmerge(qhT *qh, const char *reason, realT maxcentrum, realT maxangle,
+             boolT vneighbors);
+void    qh_all_merges(qhT *qh, boolT othermerge, boolT vneighbors);
+void    qh_appendmergeset(qhT *qh, facetT *facet, facetT *neighbor, mergeType mergetype, realT *angle);
+setT   *qh_basevertices(qhT *qh, facetT *samecycle);
+void    qh_checkconnect(qhT *qh /* qh.new_facets */);
+boolT   qh_checkzero(qhT *qh, boolT testall);
+int     qh_compareangle(const void *p1, const void *p2);
+int     qh_comparemerge(const void *p1, const void *p2);
+int     qh_comparevisit(const void *p1, const void *p2);
+void    qh_copynonconvex(qhT *qh, ridgeT *atridge);
+void    qh_degen_redundant_facet(qhT *qh, facetT *facet);
+void    qh_degen_redundant_neighbors(qhT *qh, facetT *facet, facetT *delfacet);
+vertexT *qh_find_newvertex(qhT *qh, vertexT *oldvertex, setT *vertices, setT *ridges);
+void    qh_findbest_test(qhT *qh, boolT testcentrum, facetT *facet, facetT *neighbor,
+           facetT **bestfacet, realT *distp, realT *mindistp, realT *maxdistp);
+facetT *qh_findbestneighbor(qhT *qh, facetT *facet, realT *distp, realT *mindistp, realT *maxdistp);
+void    qh_flippedmerges(qhT *qh, facetT *facetlist, boolT *wasmerge);
+void    qh_forcedmerges(qhT *qh, boolT *wasmerge);
+void    qh_getmergeset(qhT *qh, facetT *facetlist);
+void    qh_getmergeset_initial(qhT *qh, facetT *facetlist);
+void    qh_hashridge(qhT *qh, setT *hashtable, int hashsize, ridgeT *ridge, vertexT *oldvertex);
+ridgeT *qh_hashridge_find(qhT *qh, setT *hashtable, int hashsize, ridgeT *ridge,
+              vertexT *vertex, vertexT *oldvertex, int *hashslot);
+void    qh_makeridges(qhT *qh, facetT *facet);
+void    qh_mark_dupridges(qhT *qh, facetT *facetlist);
+void    qh_maydropneighbor(qhT *qh, facetT *facet);
+int     qh_merge_degenredundant(qhT *qh);
+void    qh_merge_nonconvex(qhT *qh, facetT *facet1, facetT *facet2, mergeType mergetype);
+void    qh_mergecycle(qhT *qh, facetT *samecycle, facetT *newfacet);
+void    qh_mergecycle_all(qhT *qh, facetT *facetlist, boolT *wasmerge);
+void    qh_mergecycle_facets(qhT *qh, facetT *samecycle, facetT *newfacet);
+void    qh_mergecycle_neighbors(qhT *qh, facetT *samecycle, facetT *newfacet);
+void    qh_mergecycle_ridges(qhT *qh, facetT *samecycle, facetT *newfacet);
+void    qh_mergecycle_vneighbors(qhT *qh, facetT *samecycle, facetT *newfacet);
+void    qh_mergefacet(qhT *qh, facetT *facet1, facetT *facet2, realT *mindist, realT *maxdist, boolT mergeapex);
+void    qh_mergefacet2d(qhT *qh, facetT *facet1, facetT *facet2);
+void    qh_mergeneighbors(qhT *qh, facetT *facet1, facetT *facet2);
+void    qh_mergeridges(qhT *qh, facetT *facet1, facetT *facet2);
+void    qh_mergesimplex(qhT *qh, facetT *facet1, facetT *facet2, boolT mergeapex);
+void    qh_mergevertex_del(qhT *qh, vertexT *vertex, facetT *facet1, facetT *facet2);
+void    qh_mergevertex_neighbors(qhT *qh, facetT *facet1, facetT *facet2);
+void    qh_mergevertices(qhT *qh, setT *vertices1, setT **vertices);
+setT   *qh_neighbor_intersections(qhT *qh, vertexT *vertex);
+void    qh_newvertices(qhT *qh, setT *vertices);
+boolT   qh_reducevertices(qhT *qh);
+vertexT *qh_redundant_vertex(qhT *qh, vertexT *vertex);
+boolT   qh_remove_extravertices(qhT *qh, facetT *facet);
+vertexT *qh_rename_sharedvertex(qhT *qh, vertexT *vertex, facetT *facet);
+void    qh_renameridgevertex(qhT *qh, ridgeT *ridge, vertexT *oldvertex, vertexT *newvertex);
+void    qh_renamevertex(qhT *qh, vertexT *oldvertex, vertexT *newvertex, setT *ridges,
+                        facetT *oldfacet, facetT *neighborA);
+boolT   qh_test_appendmerge(qhT *qh, facetT *facet, facetT *neighbor);
+boolT   qh_test_vneighbors(qhT *qh /* qh.newfacet_list */);
+void    qh_tracemerge(qhT *qh, facetT *facet1, facetT *facet2);
+void    qh_tracemerging(qhT *qh);
+void    qh_updatetested(qhT *qh, facetT *facet1, facetT *facet2);
+setT   *qh_vertexridges(qhT *qh, vertexT *vertex);
+void    qh_vertexridges_facet(qhT *qh, vertexT *vertex, facetT *facet, setT **ridges);
+void    qh_willdelete(qhT *qh, facetT *facet, facetT *replace);
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /* qhDEFmerge */
diff --git a/Cdir/poly_r.h b/Cdir/poly_r.h
new file mode 100644
--- /dev/null
+++ b/Cdir/poly_r.h
@@ -0,0 +1,303 @@
+/*
  ---------------------------------
+
+   poly_r.h
+   header file for poly_r.c and poly2_r.c
+
+   see qh-poly_r.htm, libqhull_r.h and poly_r.c
+
+   Copyright (c) 1993-2015 The Geometry Center.
+   $Id: //main/2015/qhull/src/libqhull_r/poly_r.h#5 $$Change: 2079 $
+   $DateTime: 2016/02/07 17:43:34 $$Author: bbarber $
+*/
+
+#ifndef qhDEFpoly
+#define qhDEFpoly 1
+
+#include "libqhull_r.h"
+
+/*===============   constants ========================== */
+
+/*----------------------------------
+
+  ALGORITHMfault
+    use as argument to checkconvex() to report errors during buildhull
+*/
+#define qh_ALGORITHMfault 0
+
+/*----------------------------------
+
+  DATAfault
+    use as argument to checkconvex() to report errors during initialhull
+*/
+#define qh_DATAfault 1
+
+/*----------------------------------
+
+  DUPLICATEridge
+    special value for facet->neighbor to indicate a duplicate ridge
+
+  notes:
+    set by matchneighbor, used by matchmatch and mark_dupridge
+*/
+#define qh_DUPLICATEridge (facetT *)1L
+
+/*----------------------------------
+
+  MERGEridge       flag in facet
+    special value for facet->neighbor to indicate a merged ridge
+
+  notes:
+    set by matchneighbor, used by matchmatch and mark_dupridge
+*/
+#define qh_MERGEridge (facetT *)2L
+
+
+/*============ -structures- ====================*/
+
+/*=========== -macros- =========================*/
+
+/*----------------------------------
+
+  FORALLfacet_( facetlist ) { ... }
+    assign 'facet' to each facet in facetlist
+
+  notes:
+    uses 'facetT *facet;'
+    assumes last facet is a sentinel
+
+  see:
+    FORALLfacets
+*/
+#define FORALLfacet_( facetlist ) if (facetlist ) for ( facet=( facetlist ); facet && facet->next; facet= facet->next )
+
+/*----------------------------------
+
+  FORALLnew_facets { ... }
+    assign 'newfacet' to each facet in qh.newfacet_list
+
+  notes:
+    uses 'facetT *newfacet;'
+    at exit, newfacet==NULL
+*/
+#define FORALLnew_facets for ( newfacet=qh->newfacet_list;newfacet && newfacet->next;newfacet=newfacet->next )
+
+/*----------------------------------
+
+  FORALLvertex_( vertexlist ) { ... }
+    assign 'vertex' to each vertex in vertexlist
+
+  notes:
+    uses 'vertexT *vertex;'
+    at exit, vertex==NULL
+*/
+#define FORALLvertex_( vertexlist ) for (vertex=( vertexlist );vertex && vertex->next;vertex= vertex->next )
+
+/*----------------------------------
+
+  FORALLvisible_facets { ... }
+    assign 'visible' to each visible facet in qh.visible_list
+
+  notes:
+    uses 'vacetT *visible;'
+    at exit, visible==NULL
+*/
+#define FORALLvisible_facets for (visible=qh->visible_list; visible && visible->visible; visible= visible->next)
+
+/*----------------------------------
+
+  FORALLsame_( newfacet ) { ... }
+    assign 'same' to each facet in newfacet->f.samecycle
+
+  notes:
+    uses 'facetT *same;'
+    stops when it returns to newfacet
+*/
+#define FORALLsame_(newfacet) for (same= newfacet->f.samecycle; same != newfacet; same= same->f.samecycle)
+
+/*----------------------------------
+
+  FORALLsame_cycle_( newfacet ) { ... }
+    assign 'same' to each facet in newfacet->f.samecycle
+
+  notes:
+    uses 'facetT *same;'
+    at exit, same == NULL
+*/
+#define FORALLsame_cycle_(newfacet) \
+     for (same= newfacet->f.samecycle; \
+         same; same= (same == newfacet ?  NULL : same->f.samecycle))
+
+/*----------------------------------
+
+  FOREACHneighborA_( facet ) { ... }
+    assign 'neighborA' to each neighbor in facet->neighbors
+
+  FOREACHneighborA_( vertex ) { ... }
+    assign 'neighborA' to each neighbor in vertex->neighbors
+
+  declare:
+    facetT *neighborA, **neighborAp;
+
+  see:
+    FOREACHsetelement_
+*/
+#define FOREACHneighborA_(facet)  FOREACHsetelement_(facetT, facet->neighbors, neighborA)
+
+/*----------------------------------
+
+  FOREACHvisible_( facets ) { ... }
+    assign 'visible' to each facet in facets
+
+  notes:
+    uses 'facetT *facet, *facetp;'
+    see FOREACHsetelement_
+*/
+#define FOREACHvisible_(facets) FOREACHsetelement_(facetT, facets, visible)
+
+/*----------------------------------
+
+  FOREACHnewfacet_( facets ) { ... }
+    assign 'newfacet' to each facet in facets
+
+  notes:
+    uses 'facetT *newfacet, *newfacetp;'
+    see FOREACHsetelement_
+*/
+#define FOREACHnewfacet_(facets) FOREACHsetelement_(facetT, facets, newfacet)
+
+/*----------------------------------
+
+  FOREACHvertexA_( vertices ) { ... }
+    assign 'vertexA' to each vertex in vertices
+
+  notes:
+    uses 'vertexT *vertexA, *vertexAp;'
+    see FOREACHsetelement_
+*/
+#define FOREACHvertexA_(vertices) FOREACHsetelement_(vertexT, vertices, vertexA)
+
+/*----------------------------------
+
+  FOREACHvertexreverse12_( vertices ) { ... }
+    assign 'vertex' to each vertex in vertices
+    reverse order of first two vertices
+
+  notes:
+    uses 'vertexT *vertex, *vertexp;'
+    see FOREACHsetelement_
+*/
+#define FOREACHvertexreverse12_(vertices) FOREACHsetelementreverse12_(vertexT, vertices, vertex)
+
+
+/*=============== prototypes poly_r.c in alphabetical order ================*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void    qh_appendfacet(qhT *qh, facetT *facet);
+void    qh_appendvertex(qhT *qh, vertexT *vertex);
+void    qh_attachnewfacets(qhT *qh /* qh.visible_list, qh.newfacet_list */);
+boolT   qh_checkflipped(qhT *qh, facetT *facet, realT *dist, boolT allerror);
+void    qh_delfacet(qhT *qh, facetT *facet);
+void    qh_deletevisible(qhT *qh /* qh.visible_list, qh.horizon_list */);
+setT   *qh_facetintersect(qhT *qh, facetT *facetA, facetT *facetB, int *skipAp,int *skipBp, int extra);
+int     qh_gethash(qhT *qh, int hashsize, setT *set, int size, int firstindex, void *skipelem);
+facetT *qh_makenewfacet(qhT *qh, setT *vertices, boolT toporient, facetT *facet);
+void    qh_makenewplanes(qhT *qh /* qh.newfacet_list */);
+facetT *qh_makenew_nonsimplicial(qhT *qh, facetT *visible, vertexT *apex, int *numnew);
+facetT *qh_makenew_simplicial(qhT *qh, facetT *visible, vertexT *apex, int *numnew);
+void    qh_matchneighbor(qhT *qh, facetT *newfacet, int newskip, int hashsize,
+                          int *hashcount);
+void    qh_matchnewfacets(qhT *qh);
+boolT   qh_matchvertices(qhT *qh, int firstindex, setT *verticesA, int skipA,
+                          setT *verticesB, int *skipB, boolT *same);
+facetT *qh_newfacet(qhT *qh);
+ridgeT *qh_newridge(qhT *qh);
+int     qh_pointid(qhT *qh, pointT *point);
+void    qh_removefacet(qhT *qh, facetT *facet);
+void    qh_removevertex(qhT *qh, vertexT *vertex);
+void    qh_updatevertices(qhT *qh);
+
+
+/*========== -prototypes poly2_r.c in alphabetical order ===========*/
+
+void    qh_addhash(void *newelem, setT *hashtable, int hashsize, int hash);
+void    qh_check_bestdist(qhT *qh);
+void    qh_check_dupridge(qhT *qh, facetT *facet1, realT dist1, facetT *facet2, realT dist2);
+void    qh_check_maxout(qhT *qh);
+void    qh_check_output(qhT *qh);
+void    qh_check_point(qhT *qh, pointT *point, facetT *facet, realT *maxoutside, realT *maxdist, facetT **errfacet1, facetT **errfacet2);
+void    qh_check_points(qhT *qh);
+void    qh_checkconvex(qhT *qh, facetT *facetlist, int fault);
+void    qh_checkfacet(qhT *qh, facetT *facet, boolT newmerge, boolT *waserrorp);
+void    qh_checkflipped_all(qhT *qh, facetT *facetlist);
+void    qh_checkpolygon(qhT *qh, facetT *facetlist);
+void    qh_checkvertex(qhT *qh, vertexT *vertex);
+void    qh_clearcenters(qhT *qh, qh_CENTER type);
+void    qh_createsimplex(qhT *qh, setT *vertices);
+void    qh_delridge(qhT *qh, ridgeT *ridge);
+void    qh_delvertex(qhT *qh, vertexT *vertex);
+setT   *qh_facet3vertex(qhT *qh, facetT *facet);
+facetT *qh_findbestfacet(qhT *qh, pointT *point, boolT bestoutside,
+           realT *bestdist, boolT *isoutside);
+facetT *qh_findbestlower(qhT *qh, facetT *upperfacet, pointT *point, realT *bestdistp, int *numpart);
+facetT *qh_findfacet_all(qhT *qh, pointT *point, realT *bestdist, boolT *isoutside,
+                          int *numpart);
+int     qh_findgood(qhT *qh, facetT *facetlist, int goodhorizon);
+void    qh_findgood_all(qhT *qh, facetT *facetlist);
+void    qh_furthestnext(qhT *qh /* qh.facet_list */);
+void    qh_furthestout(qhT *qh, facetT *facet);
+void    qh_infiniteloop(qhT *qh, facetT *facet);
+void    qh_initbuild(qhT *qh);
+void    qh_initialhull(qhT *qh, setT *vertices);
+setT   *qh_initialvertices(qhT *qh, int dim, setT *maxpoints, pointT *points, int numpoints);
+vertexT *qh_isvertex(pointT *point, setT *vertices);
+vertexT *qh_makenewfacets(qhT *qh, pointT *point /*horizon_list, visible_list*/);
+void    qh_matchduplicates(qhT *qh, facetT *atfacet, int atskip, int hashsize, int *hashcount);
+void    qh_nearcoplanar(qhT *qh /* qh.facet_list */);
+vertexT *qh_nearvertex(qhT *qh, facetT *facet, pointT *point, realT *bestdistp);
+int     qh_newhashtable(qhT *qh, int newsize);
+vertexT *qh_newvertex(qhT *qh, pointT *point);
+ridgeT *qh_nextridge3d(ridgeT *atridge, facetT *facet, vertexT **vertexp);
+void    qh_outcoplanar(qhT *qh /* qh.facet_list */);
+pointT *qh_point(qhT *qh, int id);
+void    qh_point_add(qhT *qh, setT *set, pointT *point, void *elem);
+setT   *qh_pointfacet(qhT *qh /*qh.facet_list*/);
+setT   *qh_pointvertex(qhT *qh /*qh.facet_list*/);
+void    qh_prependfacet(qhT *qh, facetT *facet, facetT **facetlist);
+void    qh_printhashtable(qhT *qh, FILE *fp);
+void    qh_printlists(qhT *qh);
+void    qh_resetlists(qhT *qh, boolT stats, boolT resetVisible /*qh.newvertex_list qh.newfacet_list qh.visible_list*/);
+void    qh_setvoronoi_all(qhT *qh);
+void    qh_triangulate(qhT *qh /*qh.facet_list*/);
+void    qh_triangulate_facet(qhT *qh, facetT *facetA, vertexT **first_vertex);
+void    qh_triangulate_link(qhT *qh, facetT *oldfacetA, facetT *facetA, facetT *oldfacetB, facetT *facetB);
+void    qh_triangulate_mirror(qhT *qh, facetT *facetA, facetT *facetB);
+void    qh_triangulate_null(qhT *qh, facetT *facetA);
+void    qh_vertexintersect(qhT *qh, setT **vertexsetA,setT *vertexsetB);
+setT   *qh_vertexintersect_new(qhT *qh, setT *vertexsetA,setT *vertexsetB);
+void    qh_vertexneighbors(qhT *qh /*qh.facet_list*/);
+boolT   qh_vertexsubset(setT *vertexsetA, setT *vertexsetB);
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /* qhDEFpoly */
diff --git a/Cdir/qhull_ra.h b/Cdir/qhull_ra.h
new file mode 100644
--- /dev/null
+++ b/Cdir/qhull_ra.h
@@ -0,0 +1,158 @@
+/*
  ---------------------------------
+
+   qhull_ra.h
+   all header files for compiling qhull with reentrant code
+   included before C++ headers for user_r.h:QHULL_CRTDBG
+
+   see qh-qhull.htm
+
+   see libqhull_r.h for user-level definitions
+
+   see user_r.h for user-definable constants
+
+   defines internal functions for libqhull_r.c global_r.c
+
+   Copyright (c) 1993-2015 The Geometry Center.
+   $Id: //main/2015/qhull/src/libqhull_r/qhull_ra.h#6 $$Change: 2079 $
+   $DateTime: 2016/02/07 17:43:34 $$Author: bbarber $
+
+   Notes:  grep for ((" and (" to catch fprintf("lkasdjf");
+           full parens around (x?y:z)
+           use '#include "libqhull_r/qhull_ra.h"' to avoid name clashes
+*/
+
+#ifndef qhDEFqhulla
+#define qhDEFqhulla 1
+
+#include "libqhull_r.h"  /* Includes user_r.h and data types */
+
+#include "stat_r.h"
+#include "random_r.h"
+#include "mem_r.h"
+#include "qset_r.h"
+#include "geom_r.h"
+#include "merge_r.h"
+#include "poly_r.h"
+#include "io_r.h"
+
+#include 
+#include 
+#include 
+#include     /* some compilers will not need float.h */
+#include 
+#include 
+#include 
+#include 
+#include 
+/*** uncomment here and qset_r.c
+     if string.h does not define memcpy()
+#include 
+*/
+
+#if qh_CLOCKtype == 2  /* defined in user_r.h from libqhull_r.h */
+#include 
+#include 
+#include 
+#endif
+
+#ifdef _MSC_VER  /* Microsoft Visual C++ -- warning level 4 */
+#pragma warning( disable : 4100)  /* unreferenced formal parameter */
+#pragma warning( disable : 4127)  /* conditional expression is constant */
+#pragma warning( disable : 4706)  /* assignment within conditional function */
+#pragma warning( disable : 4996)  /* function was declared deprecated(strcpy, localtime, etc.) */
+#endif
+
+/* ======= -macros- =========== */
+
+/*----------------------------------
+
+  traceN((qh, qh->ferr, 0Nnnn, "format\n", vars));
+    calls qh_fprintf if qh.IStracing >= N
+
+    Add debugging traps to the end of qh_fprintf
+
+  notes:
+    removing tracing reduces code size but doesn't change execution speed
+*/
+#ifndef qh_NOtrace
+#define trace0(args) {if (qh->IStracing) qh_fprintf args;}
+#define trace1(args) {if (qh->IStracing >= 1) qh_fprintf args;}
+#define trace2(args) {if (qh->IStracing >= 2) qh_fprintf args;}
+#define trace3(args) {if (qh->IStracing >= 3) qh_fprintf args;}
+#define trace4(args) {if (qh->IStracing >= 4) qh_fprintf args;}
+#define trace5(args) {if (qh->IStracing >= 5) qh_fprintf args;}
+#else /* qh_NOtrace */
+#define trace0(args) {}
+#define trace1(args) {}
+#define trace2(args) {}
+#define trace3(args) {}
+#define trace4(args) {}
+#define trace5(args) {}
+#endif /* qh_NOtrace */
+
+/*----------------------------------
+
+  Define an unused variable to avoid compiler warnings
+
+  Derived from Qt's corelib/global/qglobal.h
+
+*/
+
+#if defined(__cplusplus) && defined(__INTEL_COMPILER) && !defined(QHULL_OS_WIN)
+template 
+inline void qhullUnused(T &x) { (void)x; }
+#  define QHULL_UNUSED(x) qhullUnused(x);
+#else
+#  define QHULL_UNUSED(x) (void)x;
+#endif
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/***** -libqhull_r.c prototypes (alphabetical after qhull) ********************/
+
+void    qh_qhull(qhT *qh);
+boolT   qh_addpoint(qhT *qh, pointT *furthest, facetT *facet, boolT checkdist);
+void    qh_buildhull(qhT *qh);
+void    qh_buildtracing(qhT *qh, pointT *furthest, facetT *facet);
+void    qh_build_withrestart(qhT *qh);
+void    qh_errexit2(qhT *qh, int exitcode, facetT *facet, facetT *otherfacet);
+void    qh_findhorizon(qhT *qh, pointT *point, facetT *facet, int *goodvisible,int *goodhorizon);
+pointT *qh_nextfurthest(qhT *qh, facetT **visible);
+void    qh_partitionall(qhT *qh, setT *vertices, pointT *points,int npoints);
+void    qh_partitioncoplanar(qhT *qh, pointT *point, facetT *facet, realT *dist);
+void    qh_partitionpoint(qhT *qh, pointT *point, facetT *facet);
+void    qh_partitionvisible(qhT *qh, boolT allpoints, int *numpoints);
+void    qh_precision(qhT *qh, const char *reason);
+void    qh_printsummary(qhT *qh, FILE *fp);
+
+/***** -global_r.c internal prototypes (alphabetical) ***********************/
+
+void    qh_appendprint(qhT *qh, qh_PRINT format);
+void    qh_freebuild(qhT *qh, boolT allmem);
+void    qh_freebuffers(qhT *qh);
+void    qh_initbuffers(qhT *qh, coordT *points, int numpoints, int dim, boolT ismalloc);
+
+/***** -stat_r.c internal prototypes (alphabetical) ***********************/
+
+void    qh_allstatA(qhT *qh);
+void    qh_allstatB(qhT *qh);
+void    qh_allstatC(qhT *qh);
+void    qh_allstatD(qhT *qh);
+void    qh_allstatE(qhT *qh);
+void    qh_allstatE2(qhT *qh);
+void    qh_allstatF(qhT *qh);
+void    qh_allstatG(qhT *qh);
+void    qh_allstatH(qhT *qh);
+void    qh_freebuffers(qhT *qh);
+void    qh_initbuffers(qhT *qh, coordT *points, int numpoints, int dim, boolT ismalloc);
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /* qhDEFqhulla */
diff --git a/Cdir/qset_r.h b/Cdir/qset_r.h
new file mode 100644
--- /dev/null
+++ b/Cdir/qset_r.h
@@ -0,0 +1,502 @@
+/*
  ---------------------------------
+
+   qset_r.h
+     header file for qset_r.c that implements set
+
+   see qh-set_r.htm and qset_r.c
+
+   only uses mem_r.c, malloc/free
+
+   for error handling, writes message and calls
+      qh_errexit(qhT *qh, qhmem_ERRqhull, NULL, NULL);
+
+   set operations satisfy the following properties:
+    - sets have a max size, the actual size (if different) is stored at the end
+    - every set is NULL terminated
+    - sets may be sorted or unsorted, the caller must distinguish this
+
+   Copyright (c) 1993-2015 The Geometry Center.
+   $Id: //main/2015/qhull/src/libqhull_r/qset_r.h#4 $$Change: 2079 $
+   $DateTime: 2016/02/07 17:43:34 $$Author: bbarber $
+*/
+
+#ifndef qhDEFset
+#define qhDEFset 1
+
+#include 
+
+/*================= -structures- ===============*/
+
+#ifndef DEFsetT
+#define DEFsetT 1
+typedef struct setT setT;   /* a set is a sorted or unsorted array of pointers */
+#endif
+
+#ifndef DEFqhT
+#define DEFqhT 1
+typedef struct qhT qhT;          /* defined in libqhull_r.h */
+#endif
+
+/* [jan'15] Decided not to use countT.  Most sets are small.  The code uses signed tests */
+
+/*------------------------------------------
+
+setT
+  a set or list of pointers with maximum size and actual size.
+
+variations:
+  unsorted, unique   -- a list of unique pointers with NULL terminator
+                           user guarantees uniqueness
+  sorted             -- a sorted list of unique pointers with NULL terminator
+                           qset_r.c guarantees uniqueness
+  unsorted           -- a list of pointers terminated with NULL
+  indexed            -- an array of pointers with NULL elements
+
+structure for set of n elements:
+
+        --------------
+        |  maxsize
+        --------------
+        |  e[0] - a pointer, may be NULL for indexed sets
+        --------------
+        |  e[1]
+
+        --------------
+        |  ...
+        --------------
+        |  e[n-1]
+        --------------
+        |  e[n] = NULL
+        --------------
+        |  ...
+        --------------
+        |  e[maxsize] - n+1 or NULL (determines actual size of set)
+        --------------
+
+*/
+
+/*-- setelemT -- internal type to allow both pointers and indices
+*/
+typedef union setelemT setelemT;
+union setelemT {
+  void    *p;
+  int   i;         /* integer used for e[maxSize] */
+};
+
+struct setT {
+  int maxsize;          /* maximum number of elements (except NULL) */
+  setelemT e[1];        /* array of pointers, tail is NULL */
+                        /* last slot (unless NULL) is actual size+1
+                           e[maxsize]==NULL or e[e[maxsize]-1]==NULL */
+                        /* this may generate a warning since e[] contains
+                           maxsize elements */
+};
+
+/*=========== -constants- =========================*/
+
+/*-------------------------------------
+
+  SETelemsize
+    size of a set element in bytes
+*/
+#define SETelemsize ((int)sizeof(setelemT))
+
+
+/*=========== -macros- =========================*/
+
+/*-------------------------------------
+
+   FOREACHsetelement_(type, set, variable)
+     define FOREACH iterator
+
+   declare:
+     assumes *variable and **variablep are declared
+     no space in "variable)" [DEC Alpha cc compiler]
+
+   each iteration:
+     variable is set element
+     variablep is one beyond variable.
+
+   to repeat an element:
+     variablep--; / *repeat* /
+
+   at exit:
+     variable is NULL at end of loop
+
+   example:
+     #define FOREACHfacet_( facets ) FOREACHsetelement_( facetT, facets, facet )
+
+   notes:
+     use FOREACHsetelement_i_() if need index or include NULLs
+
+   WARNING:
+     nested loops can't use the same variable (define another FOREACH)
+
+     needs braces if nested inside another FOREACH
+     this includes intervening blocks, e.g. FOREACH...{ if () FOREACH...} )
+*/
+#define FOREACHsetelement_(type, set, variable) \
+        if (((variable= NULL), set)) for (\
+          variable##p= (type **)&((set)->e[0].p); \
+          (variable= *variable##p++);)
+
+/*------------------------------------------
+
+   FOREACHsetelement_i_(qh, type, set, variable)
+     define indexed FOREACH iterator
+
+   declare:
+     type *variable, variable_n, variable_i;
+
+   each iteration:
+     variable is set element, may be NULL
+     variable_i is index, variable_n is qh_setsize()
+
+   to repeat an element:
+     variable_i--; variable_n-- repeats for deleted element
+
+   at exit:
+     variable==NULL and variable_i==variable_n
+
+   example:
+     #define FOREACHfacet_i_( qh, facets ) FOREACHsetelement_i_( qh, facetT, facets, facet )
+
+   WARNING:
+     nested loops can't use the same variable (define another FOREACH)
+
+     needs braces if nested inside another FOREACH
+     this includes intervening blocks, e.g. FOREACH...{ if () FOREACH...} )
+*/
+#define FOREACHsetelement_i_(qh, type, set, variable) \
+        if (((variable= NULL), set)) for (\
+          variable##_i= 0, variable= (type *)((set)->e[0].p), \
+                   variable##_n= qh_setsize(qh, set);\
+          variable##_i < variable##_n;\
+          variable= (type *)((set)->e[++variable##_i].p) )
+
+/*----------------------------------------
+
+   FOREACHsetelementreverse_(qh, type, set, variable)-
+     define FOREACH iterator in reverse order
+
+   declare:
+     assumes *variable and **variablep are declared
+     also declare 'int variabletemp'
+
+   each iteration:
+     variable is set element
+
+   to repeat an element:
+     variabletemp++; / *repeat* /
+
+   at exit:
+     variable is NULL
+
+   example:
+     #define FOREACHvertexreverse_( vertices ) FOREACHsetelementreverse_( vertexT, vertices, vertex )
+
+   notes:
+     use FOREACHsetelementreverse12_() to reverse first two elements
+     WARNING: needs braces if nested inside another FOREACH
+*/
+#define FOREACHsetelementreverse_(qh, type, set, variable) \
+        if (((variable= NULL), set)) for (\
+           variable##temp= qh_setsize(qh, set)-1, variable= qh_setlast(qh, set);\
+           variable; variable= \
+           ((--variable##temp >= 0) ? SETelemt_(set, variable##temp, type) : NULL))
+
+/*-------------------------------------
+
+   FOREACHsetelementreverse12_(type, set, variable)-
+     define FOREACH iterator with e[1] and e[0] reversed
+
+   declare:
+     assumes *variable and **variablep are declared
+
+   each iteration:
+     variable is set element
+     variablep is one after variable.
+
+   to repeat an element:
+     variablep--; / *repeat* /
+
+   at exit:
+     variable is NULL at end of loop
+
+   example
+     #define FOREACHvertexreverse12_( vertices ) FOREACHsetelementreverse12_( vertexT, vertices, vertex )
+
+   notes:
+     WARNING: needs braces if nested inside another FOREACH
+*/
+#define FOREACHsetelementreverse12_(type, set, variable) \
+        if (((variable= NULL), set)) for (\
+          variable##p= (type **)&((set)->e[1].p); \
+          (variable= *variable##p); \
+          variable##p == ((type **)&((set)->e[0].p))?variable##p += 2: \
+              (variable##p == ((type **)&((set)->e[1].p))?variable##p--:variable##p++))
+
+/*-------------------------------------
+
+   FOREACHelem_( set )-
+     iterate elements in a set
+
+   declare:
+     void *elem, *elemp;
+
+   each iteration:
+     elem is set element
+     elemp is one beyond
+
+   to repeat an element:
+     elemp--; / *repeat* /
+
+   at exit:
+     elem == NULL at end of loop
+
+   example:
+     FOREACHelem_(set) {
+
+   notes:
+     WARNING: needs braces if nested inside another FOREACH
+*/
+#define FOREACHelem_(set) FOREACHsetelement_(void, set, elem)
+
+/*-------------------------------------
+
+   FOREACHset_( set )-
+     iterate a set of sets
+
+   declare:
+     setT *set, **setp;
+
+   each iteration:
+     set is set element
+     setp is one beyond
+
+   to repeat an element:
+     setp--; / *repeat* /
+
+   at exit:
+     set == NULL at end of loop
+
+   example
+     FOREACHset_(sets) {
+
+   notes:
+     WARNING: needs braces if nested inside another FOREACH
+*/
+#define FOREACHset_(sets) FOREACHsetelement_(setT, sets, set)
+
+/*-------------------------------------------
+
+   SETindex_( set, elem )
+     return index of elem in set
+
+   notes:
+     for use with FOREACH iteration
+     WARN64 -- Maximum set size is 2G
+
+   example:
+     i= SETindex_(ridges, ridge)
+*/
+#define SETindex_(set, elem) ((int)((void **)elem##p - (void **)&(set)->e[1].p))
+
+/*-----------------------------------------
+
+   SETref_( elem )
+     l.h.s. for modifying the current element in a FOREACH iteration
+
+   example:
+     SETref_(ridge)= anotherridge;
+*/
+#define SETref_(elem) (elem##p[-1])
+
+/*-----------------------------------------
+
+   SETelem_(set, n)
+     return the n'th element of set
+
+   notes:
+      assumes that n is valid [0..size] and that set is defined
+      use SETelemt_() for type cast
+*/
+#define SETelem_(set, n)           ((set)->e[n].p)
+
+/*-----------------------------------------
+
+   SETelemt_(set, n, type)
+     return the n'th element of set as a type
+
+   notes:
+      assumes that n is valid [0..size] and that set is defined
+*/
+#define SETelemt_(set, n, type)    ((type*)((set)->e[n].p))
+
+/*-----------------------------------------
+
+   SETelemaddr_(set, n, type)
+     return address of the n'th element of a set
+
+   notes:
+      assumes that n is valid [0..size] and set is defined
+*/
+#define SETelemaddr_(set, n, type) ((type **)(&((set)->e[n].p)))
+
+/*-----------------------------------------
+
+   SETfirst_(set)
+     return first element of set
+
+*/
+#define SETfirst_(set)             ((set)->e[0].p)
+
+/*-----------------------------------------
+
+   SETfirstt_(set, type)
+     return first element of set as a type
+
+*/
+#define SETfirstt_(set, type)      ((type*)((set)->e[0].p))
+
+/*-----------------------------------------
+
+   SETsecond_(set)
+     return second element of set
+
+*/
+#define SETsecond_(set)            ((set)->e[1].p)
+
+/*-----------------------------------------
+
+   SETsecondt_(set, type)
+     return second element of set as a type
+*/
+#define SETsecondt_(set, type)     ((type*)((set)->e[1].p))
+
+/*-----------------------------------------
+
+   SETaddr_(set, type)
+       return address of set's elements
+*/
+#define SETaddr_(set,type)         ((type **)(&((set)->e[0].p)))
+
+/*-----------------------------------------
+
+   SETreturnsize_(set, size)
+     return size of a set
+
+   notes:
+      set must be defined
+      use qh_setsize(qhT *qh, set) unless speed is critical
+*/
+#define SETreturnsize_(set, size) (((size)= ((set)->e[(set)->maxsize].i))?(--(size)):((size)= (set)->maxsize))
+
+/*-----------------------------------------
+
+   SETempty_(set)
+     return true(1) if set is empty
+
+   notes:
+      set may be NULL
+*/
+#define SETempty_(set)            (!set || (SETfirst_(set) ? 0 : 1))
+
+/*---------------------------------
+
+  SETsizeaddr_(set)
+    return pointer to 'actual size+1' of set (set CANNOT be NULL!!)
+    Its type is setelemT* for strict aliasing
+    All SETelemaddr_ must be cast to setelemT
+
+
+  notes:
+    *SETsizeaddr==NULL or e[*SETsizeaddr-1].p==NULL
+*/
+#define SETsizeaddr_(set) (&((set)->e[(set)->maxsize]))
+
+/*-----------------------------------------
+
+   SETtruncate_(set, size)
+     truncate set to size
+
+   see:
+     qh_settruncate()
+
+*/
+#define SETtruncate_(set, size) {set->e[set->maxsize].i= size+1; /* maybe overwritten */ \
+      set->e[size].p= NULL;}
+
+/*======= prototypes in alphabetical order ============*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void  qh_setaddsorted(qhT *qh, setT **setp, void *elem);
+void  qh_setaddnth(qhT *qh, setT **setp, int nth, void *newelem);
+void  qh_setappend(qhT *qh, setT **setp, void *elem);
+void  qh_setappend_set(qhT *qh, setT **setp, setT *setA);
+void  qh_setappend2ndlast(qhT *qh, setT **setp, void *elem);
+void  qh_setcheck(qhT *qh, setT *set, const char *tname, unsigned id);
+void  qh_setcompact(qhT *qh, setT *set);
+setT *qh_setcopy(qhT *qh, setT *set, int extra);
+void *qh_setdel(setT *set, void *elem);
+void *qh_setdellast(setT *set);
+void *qh_setdelnth(qhT *qh, setT *set, int nth);
+void *qh_setdelnthsorted(qhT *qh, setT *set, int nth);
+void *qh_setdelsorted(setT *set, void *newelem);
+setT *qh_setduplicate(qhT *qh, setT *set, int elemsize);
+void **qh_setendpointer(setT *set);
+int   qh_setequal(setT *setA, setT *setB);
+int   qh_setequal_except(setT *setA, void *skipelemA, setT *setB, void *skipelemB);
+int   qh_setequal_skip(setT *setA, int skipA, setT *setB, int skipB);
+void  qh_setfree(qhT *qh, setT **set);
+void  qh_setfree2(qhT *qh, setT **setp, int elemsize);
+void  qh_setfreelong(qhT *qh, setT **set);
+int   qh_setin(setT *set, void *setelem);
+int qh_setindex(setT *set, void *setelem);
+void  qh_setlarger(qhT *qh, setT **setp);
+void *qh_setlast(setT *set);
+setT *qh_setnew(qhT *qh, int size);
+setT *qh_setnew_delnthsorted(qhT *qh, setT *set, int size, int nth, int prepend);
+void  qh_setprint(qhT *qh, FILE *fp, const char* string, setT *set);
+void  qh_setreplace(qhT *qh, setT *set, void *oldelem, void *newelem);
+int qh_setsize(qhT *qh, setT *set);
+setT *qh_settemp(qhT *qh, int setsize);
+void  qh_settempfree(qhT *qh, setT **set);
+void  qh_settempfree_all(qhT *qh);
+setT *qh_settemppop(qhT *qh);
+void  qh_settemppush(qhT *qh, setT *set);
+void  qh_settruncate(qhT *qh, setT *set, int size);
+int   qh_setunique(qhT *qh, setT **set, void *elem);
+void  qh_setzero(qhT *qh, setT *set, int idx, int size);
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /* qhDEFset */
diff --git a/Cdir/random_r.h b/Cdir/random_r.h
new file mode 100644
--- /dev/null
+++ b/Cdir/random_r.h
@@ -0,0 +1,41 @@
+/*
  ---------------------------------
+
+  random.h
+    header file for random and utility routines
+
+   see qh-geom_r.htm and random_r.c
+
+   Copyright (c) 1993-2015 The Geometry Center.
+   $Id: //main/2015/qhull/src/libqhull_r/random_r.h#4 $$Change: 2079 $
+   $DateTime: 2016/02/07 17:43:34 $$Author: bbarber $
+*/
+
+#ifndef qhDEFrandom
+#define qhDEFrandom 1
+
+#include "libqhull_r.h"
+
+/*============= prototypes in alphabetical order ======= */
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int     qh_argv_to_command(int argc, char *argv[], char* command, int max_size);
+int     qh_argv_to_command_size(int argc, char *argv[]);
+int     qh_rand(qhT *qh);
+void    qh_srand(qhT *qh, int seed);
+realT   qh_randomfactor(qhT *qh, realT scale, realT offset);
+void    qh_randommatrix(qhT *qh, realT *buffer, int dim, realT **row);
+int     qh_strtol(const char *s, char **endp);
+double  qh_strtod(const char *s, char **endp);
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /* qhDEFrandom */
+
+
+
diff --git a/Cdir/stat_r.h b/Cdir/stat_r.h
new file mode 100644
--- /dev/null
+++ b/Cdir/stat_r.h
@@ -0,0 +1,533 @@
+/*
  ---------------------------------
+
+   stat_r.h
+     contains all statistics that are collected for qhull
+
+   see qh-stat_r.htm and stat_r.c
+
+   Copyright (c) 1993-2015 The Geometry Center.
+   $Id: //main/2015/qhull/src/libqhull_r/stat_r.h#5 $$Change: 2079 $
+   $DateTime: 2016/02/07 17:43:34 $$Author: bbarber $
+
+   recompile qhull if you change this file
+
+   Integer statistics are Z* while real statistics are W*.
+
+   define MAYdebugx to call a routine at every statistic event
+
+*/
+
+#ifndef qhDEFstat
+#define qhDEFstat 1
+
+/* Depends on realT.  Do not include libqhull_r to avoid circular dependency */
+
+#ifndef DEFqhT
+#define DEFqhT 1
+typedef struct qhT qhT;         /* Defined by libqhull_r.h */
+#endif
+
+#ifndef DEFqhstatT
+#define DEFqhstatT 1
+typedef struct qhstatT qhstatT; /* Defined here */
+#endif
+
+/*---------------------------------
+
+  qh_KEEPstatistics
+    0 turns off statistic gathering (except zzdef/zzinc/zzadd/zzval/wwval)
+*/
+#ifndef qh_KEEPstatistics
+#define qh_KEEPstatistics 1
+#endif
+
+/*---------------------------------
+
+  Zxxx for integers, Wxxx for reals
+
+  notes:
+    be sure that all statistics are defined in stat_r.c
+      otherwise initialization may core dump
+    can pick up all statistics by:
+      grep '[zw].*_[(][ZW]' *.c >z.x
+    remove trailers with query">-
+    remove leaders with  query-replace-regexp [ ^I]+  (
+*/
+#if qh_KEEPstatistics
+enum qh_statistics {     /* alphabetical after Z/W */
+    Zacoplanar,
+    Wacoplanarmax,
+    Wacoplanartot,
+    Zangle,
+    Wangle,
+    Wanglemax,
+    Wanglemin,
+    Zangletests,
+    Wareatot,
+    Wareamax,
+    Wareamin,
+    Zavoidold,
+    Wavoidoldmax,
+    Wavoidoldtot,
+    Zback0,
+    Zbestcentrum,
+    Zbestdist,
+    Zbestlower,
+    Zbestlowerall,
+    Zbestloweralln,
+    Zbestlowerv,
+    Zcentrumtests,
+    Zcheckpart,
+    Zcomputefurthest,
+    Zconcave,
+    Wconcavemax,
+    Wconcavetot,
+    Zconcaveridges,
+    Zconcaveridge,
+    Zcoplanar,
+    Wcoplanarmax,
+    Wcoplanartot,
+    Zcoplanarangle,
+    Zcoplanarcentrum,
+    Zcoplanarhorizon,
+    Zcoplanarinside,
+    Zcoplanarpart,
+    Zcoplanarridges,
+    Wcpu,
+    Zcyclefacetmax,
+    Zcyclefacettot,
+    Zcyclehorizon,
+    Zcyclevertex,
+    Zdegen,
+    Wdegenmax,
+    Wdegentot,
+    Zdegenvertex,
+    Zdelfacetdup,
+    Zdelridge,
+    Zdelvertextot,
+    Zdelvertexmax,
+    Zdetsimplex,
+    Zdistcheck,
+    Zdistconvex,
+    Zdistgood,
+    Zdistio,
+    Zdistplane,
+    Zdiststat,
+    Zdistvertex,
+    Zdistzero,
+    Zdoc1,
+    Zdoc2,
+    Zdoc3,
+    Zdoc4,
+    Zdoc5,
+    Zdoc6,
+    Zdoc7,
+    Zdoc8,
+    Zdoc9,
+    Zdoc10,
+    Zdoc11,
+    Zdoc12,
+    Zdropdegen,
+    Zdropneighbor,
+    Zdupflip,
+    Zduplicate,
+    Wduplicatemax,
+    Wduplicatetot,
+    Zdupridge,
+    Zdupsame,
+    Zflipped,
+    Wflippedmax,
+    Wflippedtot,
+    Zflippedfacets,
+    Zfindbest,
+    Zfindbestmax,
+    Zfindbesttot,
+    Zfindcoplanar,
+    Zfindfail,
+    Zfindhorizon,
+    Zfindhorizonmax,
+    Zfindhorizontot,
+    Zfindjump,
+    Zfindnew,
+    Zfindnewmax,
+    Zfindnewtot,
+    Zfindnewjump,
+    Zfindnewsharp,
+    Zgauss0,
+    Zgoodfacet,
+    Zhashlookup,
+    Zhashridge,
+    Zhashridgetest,
+    Zhashtests,
+    Zinsidevisible,
+    Zintersect,
+    Zintersectfail,
+    Zintersectmax,
+    Zintersectnum,
+    Zintersecttot,
+    Zmaxneighbors,
+    Wmaxout,
+    Wmaxoutside,
+    Zmaxridges,
+    Zmaxvertex,
+    Zmaxvertices,
+    Zmaxvneighbors,
+    Zmemfacets,
+    Zmempoints,
+    Zmemridges,
+    Zmemvertices,
+    Zmergeflipdup,
+    Zmergehorizon,
+    Zmergeinittot,
+    Zmergeinitmax,
+    Zmergeinittot2,
+    Zmergeintohorizon,
+    Zmergenew,
+    Zmergesettot,
+    Zmergesetmax,
+    Zmergesettot2,
+    Zmergesimplex,
+    Zmergevertex,
+    Wmindenom,
+    Wminvertex,
+    Zminnorm,
+    Zmultiridge,
+    Znearlysingular,
+    Zneighbor,
+    Wnewbalance,
+    Wnewbalance2,
+    Znewfacettot,
+    Znewfacetmax,
+    Znewvertex,
+    Wnewvertex,
+    Wnewvertexmax,
+    Znoarea,
+    Znonsimplicial,
+    Znowsimplicial,
+    Znotgood,
+    Znotgoodnew,
+    Znotmax,
+    Znumfacets,
+    Znummergemax,
+    Znummergetot,
+    Znumneighbors,
+    Znumridges,
+    Znumvertices,
+    Znumvisibility,
+    Znumvneighbors,
+    Zonehorizon,
+    Zpartangle,
+    Zpartcoplanar,
+    Zpartflip,
+    Zparthorizon,
+    Zpartinside,
+    Zpartition,
+    Zpartitionall,
+    Zpartnear,
+    Zpbalance,
+    Wpbalance,
+    Wpbalance2,
+    Zpostfacets,
+    Zpremergetot,
+    Zprocessed,
+    Zremvertex,
+    Zremvertexdel,
+    Zrenameall,
+    Zrenamepinch,
+    Zrenameshare,
+    Zretry,
+    Wretrymax,
+    Zridge,
+    Wridge,
+    Wridgemax,
+    Zridge0,
+    Wridge0,
+    Wridge0max,
+    Zridgemid,
+    Wridgemid,
+    Wridgemidmax,
+    Zridgeok,
+    Wridgeok,
+    Wridgeokmax,
+    Zsearchpoints,
+    Zsetplane,
+    Ztestvneighbor,
+    Ztotcheck,
+    Ztothorizon,
+    Ztotmerge,
+    Ztotpartcoplanar,
+    Ztotpartition,
+    Ztotridges,
+    Ztotvertices,
+    Ztotvisible,
+    Ztricoplanar,
+    Ztricoplanarmax,
+    Ztricoplanartot,
+    Ztridegen,
+    Ztrimirror,
+    Ztrinull,
+    Wvertexmax,
+    Wvertexmin,
+    Zvertexridge,
+    Zvertexridgetot,
+    Zvertexridgemax,
+    Zvertices,
+    Zvisfacettot,
+    Zvisfacetmax,
+    Zvisit,
+    Zvisit2max,
+    Zvisvertextot,
+    Zvisvertexmax,
+    Zvvisit,
+    Zvvisit2max,
+    Zwidefacet,
+    Zwidevertices,
+    ZEND};
+
+/*---------------------------------
+
+  Zxxx/Wxxx statistics that remain defined if qh_KEEPstatistics=0
+
+  notes:
+    be sure to use zzdef, zzinc, etc. with these statistics (no double checking!)
+*/
+#else
+enum qh_statistics {     /* for zzdef etc. macros */
+  Zback0,
+  Zbestdist,
+  Zcentrumtests,
+  Zcheckpart,
+  Zconcaveridges,
+  Zcoplanarhorizon,
+  Zcoplanarpart,
+  Zcoplanarridges,
+  Zcyclefacettot,
+  Zcyclehorizon,
+  Zdelvertextot,
+  Zdistcheck,
+  Zdistconvex,
+  Zdistzero,
+  Zdoc1,
+  Zdoc2,
+  Zdoc3,
+  Zdoc11,
+  Zflippedfacets,
+  Zgauss0,
+  Zminnorm,
+  Zmultiridge,
+  Znearlysingular,
+  Wnewvertexmax,
+  Znumvisibility,
+  Zpartcoplanar,
+  Zpartition,
+  Zpartitionall,
+  Zprocessed,
+  Zretry,
+  Zridge,
+  Wridge,
+  Wridgemax,
+  Zridge0,
+  Wridge0,
+  Wridge0max,
+  Zridgemid,
+  Wridgemid,
+  Wridgemidmax,
+  Zridgeok,
+  Wridgeok,
+  Wridgeokmax,
+  Zsetplane,
+  Ztotcheck,
+  Ztotmerge,
+    ZEND};
+#endif
+
+/*---------------------------------
+
+  ztype
+    the type of a statistic sets its initial value.
+
+  notes:
+    The type should be the same as the macro for collecting the statistic
+*/
+enum ztypes {zdoc,zinc,zadd,zmax,zmin,ZTYPEreal,wadd,wmax,wmin,ZTYPEend};
+
+/*========== macros and constants =============*/
+
+/*----------------------------------
+
+  MAYdebugx
+    define as maydebug() to be called frequently for error trapping
+*/
+#define MAYdebugx
+
+/*----------------------------------
+
+  zzdef_, zdef_( type, name, doc, -1)
+    define a statistic (assumes 'qhstat.next= 0;')
+
+  zdef_( type, name, doc, count)
+    define an averaged statistic
+    printed as name/count
+*/
+#define zzdef_(stype,name,string,cnt) qh->qhstat.id[qh->qhstat.next++]=name; \
+   qh->qhstat.doc[name]= string; qh->qhstat.count[name]= cnt; qh->qhstat.type[name]= stype
+#if qh_KEEPstatistics
+#define zdef_(stype,name,string,cnt) qh->qhstat.id[qh->qhstat.next++]=name; \
+   qh->qhstat.doc[name]= string; qh->qhstat.count[name]= cnt; qh->qhstat.type[name]= stype
+#else
+#define zdef_(type,name,doc,count)
+#endif
+
+/*----------------------------------
+
+  zzinc_( name ), zinc_( name)
+    increment an integer statistic
+*/
+#define zzinc_(id) {MAYdebugx; qh->qhstat.stats[id].i++;}
+#if qh_KEEPstatistics
+#define zinc_(id) {MAYdebugx; qh->qhstat.stats[id].i++;}
+#else
+#define zinc_(id) {}
+#endif
+
+/*----------------------------------
+
+  zzadd_( name, value ), zadd_( name, value ), wadd_( name, value )
+    add value to an integer or real statistic
+*/
+#define zzadd_(id, val) {MAYdebugx; qh->qhstat.stats[id].i += (val);}
+#define wwadd_(id, val) {MAYdebugx; qh->qhstat.stats[id].r += (val);}
+#if qh_KEEPstatistics
+#define zadd_(id, val) {MAYdebugx; qh->qhstat.stats[id].i += (val);}
+#define wadd_(id, val) {MAYdebugx; qh->qhstat.stats[id].r += (val);}
+#else
+#define zadd_(id, val) {}
+#define wadd_(id, val) {}
+#endif
+
+/*----------------------------------
+
+  zzval_( name ), zval_( name ), wwval_( name )
+    set or return value of a statistic
+*/
+#define zzval_(id) ((qh->qhstat.stats[id]).i)
+#define wwval_(id) ((qh->qhstat.stats[id]).r)
+#if qh_KEEPstatistics
+#define zval_(id) ((qh->qhstat.stats[id]).i)
+#define wval_(id) ((qh->qhstat.stats[id]).r)
+#else
+#define zval_(id) qh->qhstat.tempi
+#define wval_(id) qh->qhstat.tempr
+#endif
+
+/*----------------------------------
+
+  zmax_( id, val ), wmax_( id, value )
+    maximize id with val
+*/
+#define wwmax_(id, val) {MAYdebugx; maximize_(qh->qhstat.stats[id].r,(val));}
+#if qh_KEEPstatistics
+#define zmax_(id, val) {MAYdebugx; maximize_(qh->qhstat.stats[id].i,(val));}
+#define wmax_(id, val) {MAYdebugx; maximize_(qh->qhstat.stats[id].r,(val));}
+#else
+#define zmax_(id, val) {}
+#define wmax_(id, val) {}
+#endif
+
+/*----------------------------------
+
+  zmin_( id, val ), wmin_( id, value )
+    minimize id with val
+*/
+#if qh_KEEPstatistics
+#define zmin_(id, val) {MAYdebugx; minimize_(qh->qhstat.stats[id].i,(val));}
+#define wmin_(id, val) {MAYdebugx; minimize_(qh->qhstat.stats[id].r,(val));}
+#else
+#define zmin_(id, val) {}
+#define wmin_(id, val) {}
+#endif
+
+/*================== stat_r.h types ==============*/
+
+
+/*----------------------------------
+
+  intrealT
+    union of integer and real, used for statistics
+*/
+typedef union intrealT intrealT;    /* union of int and realT */
+union intrealT {
+    int i;
+    realT r;
+};
+
+/*----------------------------------
+
+  qhstat
+    Data structure for statistics, similar to qh and qhrbox
+
+    Allocated as part of qhT (libqhull_r.h)
+*/
+
+struct qhstatT {
+  intrealT   stats[ZEND];     /* integer and real statistics */
+  unsigned   char id[ZEND+10]; /* id's in print order */
+  const char *doc[ZEND];       /* array of documentation strings */
+  short int  count[ZEND];     /* -1 if none, else index of count to use */
+  char       type[ZEND];      /* type, see ztypes above */
+  char       printed[ZEND];   /* true, if statistic has been printed */
+  intrealT   init[ZTYPEend];  /* initial values by types, set initstatistics */
+
+  int        next;            /* next index for zdef_ */
+  int        precision;       /* index for precision problems */
+  int        vridges;         /* index for Voronoi ridges */
+  int        tempi;
+  realT      tempr;
+};
+
+/*========== function prototypes ===========*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void    qh_allstatA(qhT *qh);
+void    qh_allstatB(qhT *qh);
+void    qh_allstatC(qhT *qh);
+void    qh_allstatD(qhT *qh);
+void    qh_allstatE(qhT *qh);
+void    qh_allstatE2(qhT *qh);
+void    qh_allstatF(qhT *qh);
+void    qh_allstatG(qhT *qh);
+void    qh_allstatH(qhT *qh);
+void    qh_allstatI(qhT *qh);
+void    qh_allstatistics(qhT *qh);
+void    qh_collectstatistics(qhT *qh);
+void    qh_initstatistics(qhT *qh);
+boolT   qh_newstats(qhT *qh, int idx, int *nextindex);
+boolT   qh_nostatistic(qhT *qh, int i);
+void    qh_printallstatistics(qhT *qh, FILE *fp, const char *string);
+void    qh_printstatistics(qhT *qh, FILE *fp, const char *string);
+void    qh_printstatlevel(qhT *qh, FILE *fp, int id);
+void    qh_printstats(qhT *qh, FILE *fp, int idx, int *nextindex);
+realT   qh_stddev(int num, realT tot, realT tot2, realT *ave);
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif
+
+#endif   /* qhDEFstat */
diff --git a/Cdir/user_r.h b/Cdir/user_r.h
new file mode 100644
--- /dev/null
+++ b/Cdir/user_r.h
@@ -0,0 +1,882 @@
+/*
  ---------------------------------
+
+   user.h
+   user redefinable constants
+
+   for each source file, user_r.h is included first
+
+   see qh-user_r.htm.  see COPYING for copyright information.
+
+   See user_r.c for sample code.
+
+   before reading any code, review libqhull_r.h for data structure definitions
+
+Sections:
+   ============= qhull library constants ======================
+   ============= data types and configuration macros ==========
+   ============= performance related constants ================
+   ============= memory constants =============================
+   ============= joggle constants =============================
+   ============= conditional compilation ======================
+   ============= -merge constants- ============================
+
+Code flags --
+  NOerrors -- the code does not call qh_errexit()
+  WARN64 -- the code may be incompatible with 64-bit pointers
+
+*/
+
+#include 
+
+#ifndef qhDEFuser
+#define qhDEFuser 1
+
+/* Derived from Qt's corelib/global/qglobal.h */
+#if !defined(SAG_COM) && !defined(__CYGWIN__) && (defined(WIN64) || defined(_WIN64) || defined(__WIN64__) || defined(WIN32) || defined(_WIN32) || defined(__WIN32__) || defined(__NT__))
+#   define QHULL_OS_WIN
+#elif defined(__MWERKS__) && defined(__INTEL__) /* Metrowerks discontinued before the release of Intel Macs */
+#   define QHULL_OS_WIN
+#endif
+
+/*============================================================*/
+/*============= qhull library constants ======================*/
+/*============================================================*/
+
+/*----------------------------------
+
+  FILENAMElen -- max length for TI and TO filenames
+
+*/
+
+#define qh_FILENAMElen 500
+
+/*----------------------------------
+
+  msgcode -- Unique message codes for qh_fprintf
+
+  If add new messages, assign these values and increment in user.h and user_r.h
+  See QhullError.h for 10000 errors.
+
+  def counters =  [27, 1048, 2059, 3026, 4068, 5003,
+     6273, 7081, 8147, 9411, 10000, 11029]
+
+  See: qh_ERR* [libqhull_r.h]
+*/
+
+#define MSG_TRACE0 0
+#define MSG_TRACE1 1000
+#define MSG_TRACE2 2000
+#define MSG_TRACE3 3000
+#define MSG_TRACE4 4000
+#define MSG_TRACE5 5000
+#define MSG_ERROR  6000   /* errors written to qh.ferr */
+#define MSG_WARNING 7000
+#define MSG_STDERR  8000  /* log messages Written to qh.ferr */
+#define MSG_OUTPUT  9000
+#define MSG_QHULL_ERROR 10000 /* errors thrown by QhullError.cpp (QHULLlastError is in QhullError.h) */
+#define MSG_FIXUP  11000  /* FIXUP QH11... */
+#define MSG_MAXLEN  3000 /* qh_printhelp_degenerate() in user.c */
+
+
+/*----------------------------------
+
+  qh_OPTIONline -- max length of an option line 'FO'
+*/
+#define qh_OPTIONline 80
+
+/*============================================================*/
+/*============= data types and configuration macros ==========*/
+/*============================================================*/
+
+/*----------------------------------
+
+  realT
+    set the size of floating point numbers
+
+  qh_REALdigits
+    maximimum number of significant digits
+
+  qh_REAL_1, qh_REAL_2n, qh_REAL_3n
+    format strings for printf
+
+  qh_REALmax, qh_REALmin
+    maximum and minimum (near zero) values
+
+  qh_REALepsilon
+    machine roundoff.  Maximum roundoff error for addition and multiplication.
+
+  notes:
+   Select whether to store floating point numbers in single precision (float)
+   or double precision (double).
+
+   Use 'float' to save about 8% in time and 25% in space.  This is particularly
+   helpful if high-d where convex hulls are space limited.  Using 'float' also
+   reduces the printed size of Qhull's output since numbers have 8 digits of
+   precision.
+
+   Use 'double' when greater arithmetic precision is needed.  This is needed
+   for Delaunay triangulations and Voronoi diagrams when you are not merging
+   facets.
+
+   If 'double' gives insufficient precision, your data probably includes
+   degeneracies.  If so you should use facet merging (done by default)
+   or exact arithmetic (see imprecision section of manual, qh-impre.htm).
+   You may also use option 'Po' to force output despite precision errors.
+
+   You may use 'long double', but many format statements need to be changed
+   and you may need a 'long double' square root routine.  S. Grundmann
+   (sg@eeiwzb.et.tu-dresden.de) has done this.  He reports that the code runs
+   much slower with little gain in precision.
+
+   WARNING: on some machines,    int f(){realT a= REALmax;return (a == REALmax);}
+      returns False.  Use (a > REALmax/2) instead of (a == REALmax).
+
+   REALfloat =   1      all numbers are 'float' type
+             =   0      all numbers are 'double' type
+*/
+#define REALfloat 0
+
+#if (REALfloat == 1)
+#define realT float
+#define REALmax FLT_MAX
+#define REALmin FLT_MIN
+#define REALepsilon FLT_EPSILON
+#define qh_REALdigits 8   /* maximum number of significant digits */
+#define qh_REAL_1 "%6.8g "
+#define qh_REAL_2n "%6.8g %6.8g\n"
+#define qh_REAL_3n "%6.8g %6.8g %6.8g\n"
+
+#elif (REALfloat == 0)
+#define realT double
+#define REALmax DBL_MAX
+#define REALmin DBL_MIN
+#define REALepsilon DBL_EPSILON
+#define qh_REALdigits 16    /* maximum number of significant digits */
+#define qh_REAL_1 "%6.16g "
+#define qh_REAL_2n "%6.16g %6.16g\n"
+#define qh_REAL_3n "%6.16g %6.16g %6.16g\n"
+
+#else
+#error unknown float option
+#endif
+
+/*----------------------------------
+
+  countT
+    The type for counts and identifiers (e.g., the number of points, vertex identifiers)
+    Currently used by C++ code-only.  Decided against using it for setT because most sets are small.
+
+    Defined as 'int' for C-code compatibility and QH11026
+
+    FIXUP QH11026 countT may be defined as a unsigned value, but several code issues need to be solved first.  See countT in Changes.txt
+*/
+
+#ifndef DEFcountT
+#define DEFcountT 1
+typedef int countT;
+#endif
+#define COUNTmax 0x7fffffff
+
+
+/*----------------------------------
+
+  qh_CPUclock
+    define the clock() function for reporting the total time spent by Qhull
+    returns CPU ticks as a 'long int'
+    qh_CPUclock is only used for reporting the total time spent by Qhull
+
+  qh_SECticks
+    the number of clock ticks per second
+
+  notes:
+    looks for CLOCKS_PER_SEC, CLOCKS_PER_SECOND, or assumes microseconds
+    to define a custom clock, set qh_CLOCKtype to 0
+
+    if your system does not use clock() to return CPU ticks, replace
+    qh_CPUclock with the corresponding function.  It is converted
+    to 'unsigned long' to prevent wrap-around during long runs.  By default,
+     defines clock_t as 'long'
+
+   Set qh_CLOCKtype to
+
+     1          for CLOCKS_PER_SEC, CLOCKS_PER_SECOND, or microsecond
+                Note:  may fail if more than 1 hour elapsed time
+
+     2          use qh_clock() with POSIX times() (see global_r.c)
+*/
+#define qh_CLOCKtype 1  /* change to the desired number */
+
+#if (qh_CLOCKtype == 1)
+
+#if defined(CLOCKS_PER_SECOND)
+#define qh_CPUclock    ((unsigned long)clock())  /* return CPU clock */
+#define qh_SECticks CLOCKS_PER_SECOND
+
+#elif defined(CLOCKS_PER_SEC)
+#define qh_CPUclock    ((unsigned long)clock())  /* return CPU clock */
+#define qh_SECticks CLOCKS_PER_SEC
+
+#elif defined(CLK_TCK)
+#define qh_CPUclock    ((unsigned long)clock())  /* return CPU clock */
+#define qh_SECticks CLK_TCK
+
+#else
+#define qh_CPUclock    ((unsigned long)clock())  /* return CPU clock */
+#define qh_SECticks 1E6
+#endif
+
+#elif (qh_CLOCKtype == 2)
+#define qh_CPUclock    qh_clock()  /* return CPU clock */
+#define qh_SECticks 100
+
+#else /* qh_CLOCKtype == ? */
+#error unknown clock option
+#endif
+
+/*----------------------------------
+
+  qh_RANDOMtype, qh_RANDOMmax, qh_RANDOMseed
+    define random number generator
+
+    qh_RANDOMint generates a random integer between 0 and qh_RANDOMmax.
+    qh_RANDOMseed sets the random number seed for qh_RANDOMint
+
+  Set qh_RANDOMtype (default 5) to:
+    1       for random() with 31 bits (UCB)
+    2       for rand() with RAND_MAX or 15 bits (system 5)
+    3       for rand() with 31 bits (Sun)
+    4       for lrand48() with 31 bits (Solaris)
+    5       for qh_rand(qh) with 31 bits (included with Qhull, requires 'qh')
+
+  notes:
+    Random numbers are used by rbox to generate point sets.  Random
+    numbers are used by Qhull to rotate the input ('QRn' option),
+    simulate a randomized algorithm ('Qr' option), and to simulate
+    roundoff errors ('Rn' option).
+
+    Random number generators differ between systems.  Most systems provide
+    rand() but the period varies.  The period of rand() is not critical
+    since qhull does not normally use random numbers.
+
+    The default generator is Park & Miller's minimal standard random
+    number generator [CACM 31:1195 '88].  It is included with Qhull.
+
+    If qh_RANDOMmax is wrong, qhull will report a warning and Geomview
+    output will likely be invisible.
+*/
+#define qh_RANDOMtype 5   /* *** change to the desired number *** */
+
+#if (qh_RANDOMtype == 1)
+#define qh_RANDOMmax ((realT)0x7fffffffUL)  /* 31 bits, random()/MAX */
+#define qh_RANDOMint random()
+#define qh_RANDOMseed_(qh, seed) srandom(seed);
+
+#elif (qh_RANDOMtype == 2)
+#ifdef RAND_MAX
+#define qh_RANDOMmax ((realT)RAND_MAX)
+#else
+#define qh_RANDOMmax ((realT)32767)   /* 15 bits (System 5) */
+#endif
+#define qh_RANDOMint  rand()
+#define qh_RANDOMseed_(qh, seed) srand((unsigned)seed);
+
+#elif (qh_RANDOMtype == 3)
+#define qh_RANDOMmax ((realT)0x7fffffffUL)  /* 31 bits, Sun */
+#define qh_RANDOMint  rand()
+#define qh_RANDOMseed_(qh, seed) srand((unsigned)seed);
+
+#elif (qh_RANDOMtype == 4)
+#define qh_RANDOMmax ((realT)0x7fffffffUL)  /* 31 bits, lrand38()/MAX */
+#define qh_RANDOMint lrand48()
+#define qh_RANDOMseed_(qh, seed) srand48(seed);
+
+#elif (qh_RANDOMtype == 5)  /* 'qh' is an implicit parameter */
+#define qh_RANDOMmax ((realT)2147483646UL)  /* 31 bits, qh_rand/MAX */
+#define qh_RANDOMint qh_rand(qh)
+#define qh_RANDOMseed_(qh, seed) qh_srand(qh, seed);
+/* unlike rand(), never returns 0 */
+
+#else
+#error: unknown random option
+#endif
+
+/*----------------------------------
+
+  qh_ORIENTclock
+    0 for inward pointing normals by Geomview convention
+*/
+#define qh_ORIENTclock 0
+
+
+/*============================================================*/
+/*============= joggle constants =============================*/
+/*============================================================*/
+
+/*----------------------------------
+
+qh_JOGGLEdefault
+default qh.JOGGLEmax is qh.DISTround * qh_JOGGLEdefault
+
+notes:
+rbox s r 100 | qhull QJ1e-15 QR0 generates 90% faults at distround 7e-16
+rbox s r 100 | qhull QJ1e-14 QR0 generates 70% faults
+rbox s r 100 | qhull QJ1e-13 QR0 generates 35% faults
+rbox s r 100 | qhull QJ1e-12 QR0 generates 8% faults
+rbox s r 100 | qhull QJ1e-11 QR0 generates 1% faults
+rbox s r 100 | qhull QJ1e-10 QR0 generates 0% faults
+rbox 1000 W0 | qhull QJ1e-12 QR0 generates 86% faults
+rbox 1000 W0 | qhull QJ1e-11 QR0 generates 20% faults
+rbox 1000 W0 | qhull QJ1e-10 QR0 generates 2% faults
+the later have about 20 points per facet, each of which may interfere
+
+pick a value large enough to avoid retries on most inputs
+*/
+#define qh_JOGGLEdefault 30000.0
+
+/*----------------------------------
+
+qh_JOGGLEincrease
+factor to increase qh.JOGGLEmax on qh_JOGGLEretry or qh_JOGGLEagain
+*/
+#define qh_JOGGLEincrease 10.0
+
+/*----------------------------------
+
+qh_JOGGLEretry
+if ZZretry = qh_JOGGLEretry, increase qh.JOGGLEmax
+
+notes:
+try twice at the original value in case of bad luck the first time
+*/
+#define qh_JOGGLEretry 2
+
+/*----------------------------------
+
+qh_JOGGLEagain
+every following qh_JOGGLEagain, increase qh.JOGGLEmax
+
+notes:
+1 is OK since it's already failed qh_JOGGLEretry times
+*/
+#define qh_JOGGLEagain 1
+
+/*----------------------------------
+
+qh_JOGGLEmaxincrease
+maximum qh.JOGGLEmax due to qh_JOGGLEincrease
+relative to qh.MAXwidth
+
+notes:
+qh.joggleinput will retry at this value until qh_JOGGLEmaxretry
+*/
+#define qh_JOGGLEmaxincrease 1e-2
+
+/*----------------------------------
+
+qh_JOGGLEmaxretry
+stop after qh_JOGGLEmaxretry attempts
+*/
+#define qh_JOGGLEmaxretry 100
+
+/*============================================================*/
+/*============= performance related constants ================*/
+/*============================================================*/
+
+/*----------------------------------
+
+  qh_HASHfactor
+    total hash slots / used hash slots.  Must be at least 1.1.
+
+  notes:
+    =2 for at worst 50% occupancy for qh.hash_table and normally 25% occupancy
+*/
+#define qh_HASHfactor 2
+
+/*----------------------------------
+
+  qh_VERIFYdirect
+    with 'Tv' verify all points against all facets if op count is smaller
+
+  notes:
+    if greater, calls qh_check_bestdist() instead
+*/
+#define qh_VERIFYdirect 1000000
+
+/*----------------------------------
+
+  qh_INITIALsearch
+     if qh_INITIALmax, search points up to this dimension
+*/
+#define qh_INITIALsearch 6
+
+/*----------------------------------
+
+  qh_INITIALmax
+    if dim >= qh_INITIALmax, use min/max coordinate points for initial simplex
+
+  notes:
+    from points with non-zero determinants
+    use option 'Qs' to override (much slower)
+*/
+#define qh_INITIALmax 8
+
+/*============================================================*/
+/*============= memory constants =============================*/
+/*============================================================*/
+
+/*----------------------------------
+
+  qh_MEMalign
+    memory alignment for qh_meminitbuffers() in global_r.c
+
+  notes:
+    to avoid bus errors, memory allocation must consider alignment requirements.
+    malloc() automatically takes care of alignment.   Since mem_r.c manages
+    its own memory, we need to explicitly specify alignment in
+    qh_meminitbuffers().
+
+    A safe choice is sizeof(double).  sizeof(float) may be used if doubles
+    do not occur in data structures and pointers are the same size.  Be careful
+    of machines (e.g., DEC Alpha) with large pointers.
+
+    If using gcc, best alignment is [fmax_() is defined in geom_r.h]
+              #define qh_MEMalign fmax_(__alignof__(realT),__alignof__(void *))
+*/
+#define qh_MEMalign ((int)(fmax_(sizeof(realT), sizeof(void *))))
+
+/*----------------------------------
+
+  qh_MEMbufsize
+    size of additional memory buffers
+
+  notes:
+    used for qh_meminitbuffers() in global_r.c
+*/
+#define qh_MEMbufsize 0x10000       /* allocate 64K memory buffers */
+
+/*----------------------------------
+
+  qh_MEMinitbuf
+    size of initial memory buffer
+
+  notes:
+    use for qh_meminitbuffers() in global_r.c
+*/
+#define qh_MEMinitbuf 0x20000      /* initially allocate 128K buffer */
+
+/*----------------------------------
+
+  qh_INFINITE
+    on output, indicates Voronoi center at infinity
+*/
+#define qh_INFINITE  -10.101
+
+/*----------------------------------
+
+  qh_DEFAULTbox
+    default box size (Geomview expects 0.5)
+
+  qh_DEFAULTbox
+    default box size for integer coorindate (rbox only)
+*/
+#define qh_DEFAULTbox 0.5
+#define qh_DEFAULTzbox 1e6
+
+/*============================================================*/
+/*============= conditional compilation ======================*/
+/*============================================================*/
+
+/*----------------------------------
+
+  __cplusplus
+    defined by C++ compilers
+
+  __MSC_VER
+    defined by Microsoft Visual C++
+
+  __MWERKS__ && __INTEL__
+    defined by Metrowerks when compiling for Windows (not Intel-based Macintosh)
+
+  __MWERKS__ && __POWERPC__
+    defined by Metrowerks when compiling for PowerPC-based Macintosh
+
+  __STDC__
+    defined for strict ANSI C
+*/
+
+/*----------------------------------
+
+  qh_COMPUTEfurthest
+    compute furthest distance to an outside point instead of storing it with the facet
+    =1 to compute furthest
+
+  notes:
+    computing furthest saves memory but costs time
+      about 40% more distance tests for partitioning
+      removes facet->furthestdist
+*/
+#define qh_COMPUTEfurthest 0
+
+/*----------------------------------
+
+  qh_KEEPstatistics
+    =0 removes most of statistic gathering and reporting
+
+  notes:
+    if 0, code size is reduced by about 4%.
+*/
+#define qh_KEEPstatistics 1
+
+/*----------------------------------
+
+  qh_MAXoutside
+    record outer plane for each facet
+    =1 to record facet->maxoutside
+
+  notes:
+    this takes a realT per facet and slightly slows down qhull
+    it produces better outer planes for geomview output
+*/
+#define qh_MAXoutside 1
+
+/*----------------------------------
+
+  qh_NOmerge
+    disables facet merging if defined
+
+  notes:
+    This saves about 10% space.
+
+    Unless 'Q0'
+      qh_NOmerge sets 'QJ' to avoid precision errors
+
+    #define qh_NOmerge
+
+  see:
+    qh_NOmem in mem_r.c
+
+    see user_r.c/user_eg.c for removing io_r.o
+*/
+
+/*----------------------------------
+
+  qh_NOtrace
+    no tracing if defined
+
+  notes:
+    This saves about 5% space.
+
+    #define qh_NOtrace
+*/
+
+#if 0  /* sample code */
+    exitcode= qh_new_qhull(qhT *qh, dim, numpoints, points, ismalloc,
+                      flags, outfile, errfile);
+    qh_freeqhull(qhT *qh, !qh_ALL); /* frees long memory used by second call */
+    qh_memfreeshort(qhT *qh, &curlong, &totlong);  /* frees short memory and memory allocator */
+#endif
+
+/*----------------------------------
+
+  qh_QUICKhelp
+    =1 to use abbreviated help messages, e.g., for degenerate inputs
+*/
+#define qh_QUICKhelp    0
+
+/*============================================================*/
+/*============= -merge constants- ============================*/
+/*============================================================*/
+/*
+   These constants effect facet merging.  You probably will not need
+   to modify them.  They effect the performance of facet merging.
+*/
+
+/*----------------------------------
+
+  qh_DIMmergeVertex
+    max dimension for vertex merging (it is not effective in high-d)
+*/
+#define qh_DIMmergeVertex 6
+
+/*----------------------------------
+
+  qh_DIMreduceBuild
+     max dimension for vertex reduction during build (slow in high-d)
+*/
+#define qh_DIMreduceBuild 5
+
+/*----------------------------------
+
+  qh_BESTcentrum
+     if > 2*dim+n vertices, qh_findbestneighbor() tests centrums (faster)
+     else, qh_findbestneighbor() tests all vertices (much better merges)
+
+  qh_BESTcentrum2
+     if qh_BESTcentrum2 * DIM3 + BESTcentrum < #vertices tests centrums
+*/
+#define qh_BESTcentrum 20
+#define qh_BESTcentrum2 2
+
+/*----------------------------------
+
+  qh_BESTnonconvex
+    if > dim+n neighbors, qh_findbestneighbor() tests nonconvex ridges.
+
+  notes:
+    It is needed because qh_findbestneighbor is slow for large facets
+*/
+#define qh_BESTnonconvex 15
+
+/*----------------------------------
+
+  qh_MAXnewmerges
+    if >n newmerges, qh_merge_nonconvex() calls qh_reducevertices_centrums.
+
+  notes:
+    It is needed because postmerge can merge many facets at once
+*/
+#define qh_MAXnewmerges 2
+
+/*----------------------------------
+
+  qh_MAXnewcentrum
+    if <= dim+n vertices (n approximates the number of merges),
+      reset the centrum in qh_updatetested() and qh_mergecycle_facets()
+
+  notes:
+    needed to reduce cost and because centrums may move too much if
+    many vertices in high-d
+*/
+#define qh_MAXnewcentrum 5
+
+/*----------------------------------
+
+  qh_COPLANARratio
+    for 3-d+ merging, qh.MINvisible is n*premerge_centrum
+
+  notes:
+    for non-merging, it's DISTround
+*/
+#define qh_COPLANARratio 3
+
+/*----------------------------------
+
+  qh_DISToutside
+    When is a point clearly outside of a facet?
+    Stops search in qh_findbestnew or qh_partitionall
+    qh_findbest uses qh.MINoutside since since it is only called if no merges.
+
+  notes:
+    'Qf' always searches for best facet
+    if !qh.MERGING, same as qh.MINoutside.
+    if qh_USEfindbestnew, increase value since neighboring facets may be ill-behaved
+      [Note: Zdelvertextot occurs normally with interior points]
+            RBOX 1000 s Z1 G1e-13 t1001188774 | QHULL Tv
+    When there is a sharp edge, need to move points to a
+    clearly good facet; otherwise may be lost in another partitioning.
+    if too big then O(n^2) behavior for partitioning in cone
+    if very small then important points not processed
+    Needed in qh_partitionall for
+      RBOX 1000 s Z1 G1e-13 t1001032651 | QHULL Tv
+    Needed in qh_findbestnew for many instances of
+      RBOX 1000 s Z1 G1e-13 t | QHULL Tv
+
+  See:
+    qh_DISToutside -- when is a point clearly outside of a facet
+    qh_SEARCHdist -- when is facet coplanar with the best facet?
+    qh_USEfindbestnew -- when to use qh_findbestnew for qh_partitionpoint()
+*/
+#define qh_DISToutside ((qh_USEfindbestnew ? 2 : 1) * \
+     fmax_((qh->MERGING ? 2 : 1)*qh->MINoutside, qh->max_outside))
+
+/*----------------------------------
+
+  qh_RATIOnearinside
+    ratio of qh.NEARinside to qh.ONEmerge for retaining inside points for
+    qh_check_maxout().
+
+  notes:
+    This is overkill since do not know the correct value.
+    It effects whether 'Qc' reports all coplanar points
+    Not used for 'd' since non-extreme points are coplanar
+*/
+#define qh_RATIOnearinside 5
+
+/*----------------------------------
+
+  qh_SEARCHdist
+    When is a facet coplanar with the best facet?
+    qh_findbesthorizon: all coplanar facets of the best facet need to be searched.
+
+  See:
+    qh_DISToutside -- when is a point clearly outside of a facet
+    qh_SEARCHdist -- when is facet coplanar with the best facet?
+    qh_USEfindbestnew -- when to use qh_findbestnew for qh_partitionpoint()
+*/
+#define qh_SEARCHdist ((qh_USEfindbestnew ? 2 : 1) * \
+      (qh->max_outside + 2 * qh->DISTround + fmax_( qh->MINvisible, qh->MAXcoplanar)));
+
+/*----------------------------------
+
+  qh_USEfindbestnew
+     Always use qh_findbestnew for qh_partitionpoint, otherwise use
+     qh_findbestnew if merged new facet or sharpnewfacets.
+
+  See:
+    qh_DISToutside -- when is a point clearly outside of a facet
+    qh_SEARCHdist -- when is facet coplanar with the best facet?
+    qh_USEfindbestnew -- when to use qh_findbestnew for qh_partitionpoint()
+*/
+#define qh_USEfindbestnew (zzval_(Ztotmerge) > 50)
+
+/*----------------------------------
+
+  qh_WIDEcoplanar
+    n*MAXcoplanar or n*MINvisible for a WIDEfacet
+
+    if vertex is further than qh.WIDEfacet from the hyperplane
+    then its ridges are not counted in computing the area, and
+    the facet's centrum is frozen.
+
+  notes:
+   qh.WIDEfacet= max(qh.MAXoutside,qh_WIDEcoplanar*qh.MAXcoplanar,
+      qh_WIDEcoplanar * qh.MINvisible);
+*/
+#define qh_WIDEcoplanar 6
+
+/*----------------------------------
+
+  qh_WIDEduplicate
+    Merge ratio for errexit from qh_forcedmerges due to duplicate ridge
+    Override with option Q12 no-wide-duplicate
+
+    Notes:
+      Merging a duplicate ridge can lead to very wide facets.
+      A future release of qhull will avoid duplicate ridges by removing duplicate sub-ridges from the horizon
+*/
+#define qh_WIDEduplicate 100
+
+/*----------------------------------
+
+  qh_MAXnarrow
+    max. cosine in initial hull that sets qh.NARROWhull
+
+  notes:
+    If qh.NARROWhull, the initial partition does not make
+    coplanar points.  If narrow, a coplanar point can be
+    coplanar to two facets of opposite orientations and
+    distant from the exact convex hull.
+
+    Conservative estimate.  Don't actually see problems until it is -1.0
+*/
+#define qh_MAXnarrow -0.99999999
+
+/*----------------------------------
+
+  qh_WARNnarrow
+    max. cosine in initial hull to warn about qh.NARROWhull
+
+  notes:
+    this is a conservative estimate.
+    Don't actually see problems until it is -1.0.  See qh-impre.htm
+*/
+#define qh_WARNnarrow -0.999999999999999
+
+/*----------------------------------
+
+  qh_ZEROdelaunay
+    a zero Delaunay facet occurs for input sites coplanar with their convex hull
+    the last normal coefficient of a zero Delaunay facet is within
+        qh_ZEROdelaunay * qh.ANGLEround of 0
+
+  notes:
+    qh_ZEROdelaunay does not allow for joggled input ('QJ').
+
+    You can avoid zero Delaunay facets by surrounding the input with a box.
+
+    Use option 'PDk:-n' to explicitly define zero Delaunay facets
+      k= dimension of input sites (e.g., 3 for 3-d Delaunay triangulation)
+      n= the cutoff for zero Delaunay facets (e.g., 'PD3:-1e-12')
+*/
+#define qh_ZEROdelaunay 2
+
+/*============================================================*/
+/*============= Microsoft DevStudio ==========================*/
+/*============================================================*/
+
+/*
+   Finding Memory Leaks Using the CRT Library
+   https://msdn.microsoft.com/en-us/library/x98tx3cf(v=vs.100).aspx
+
+   Reports enabled in qh_lib_check for Debug window and stderr
+
+   From 2005=>msvcr80d, 2010=>msvcr100d, 2012=>msvcr110d
+
+   Watch: {,,msvcr80d.dll}_crtBreakAlloc  Value from {n} in the leak report
+   _CrtSetBreakAlloc(689); // qh_lib_check() [global_r.c]
+
+   Examples
+     http://free-cad.sourceforge.net/SrcDocu/d2/d7f/MemDebug_8cpp_source.html
+     https://github.com/illlust/Game/blob/master/library/MemoryLeak.cpp
+*/
+#if 0   /* off (0) by default for QHULL_CRTDBG */
+#define QHULL_CRTDBG
+#endif
+
+#if defined(_MSC_VER) && defined(_DEBUG) && defined(QHULL_CRTDBG)
+#define _CRTDBG_MAP_ALLOC
+#include 
+#include 
+#endif
+
+#endif /* qh_DEFuser */
+
+
+
diff --git a/Cdir/utils.h b/Cdir/utils.h
new file mode 100644
--- /dev/null
+++ b/Cdir/utils.h
@@ -0,0 +1,20 @@
+
+int cmpfunc(const void*, const void*);
+int cmpfuncdbl(const void*, const void*);
+void qsortu(unsigned*, unsigned);
+
+double square(double);
+
+void appendu(unsigned, unsigned**, unsigned, unsigned*);
+
+double* nanvector(int);
+
+double* middle(double*, double*, unsigned);
+
+double* getpoint(double*, unsigned, unsigned);
+
+double dotproduct(double*, double*, unsigned);
+
+unsigned* uzeros(unsigned);
+
+double squaredDistance(double*, double*, unsigned);
diff --git a/ChangeLog.md b/ChangeLog.md
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -13,3 +13,5 @@
     there was a hint on stackoverflow that single char dirs do not work well in cabal
     works with stack once, after stack install uniform-geometry can use it.
 -     cannot go to hpack, trying stack init - not working either
+0.1.0.4 
+the package failed to build on hackage; it builds with stack, but fails with cabal
diff --git a/qhull.cabal b/qhull.cabal
--- a/qhull.cabal
+++ b/qhull.cabal
@@ -1,5 +1,5 @@
 name:                qhull
-version:             0.1.0.3
+version:             0.1.0.4 
 synopsis:       Delaunay triangulation, Voronoi diagrams and convex hulls. 
 description:    Based on the qhull C library. 
 
@@ -112,6 +112,22 @@
                      , ./Cdir/convexhull.c
                      , ./Cdir/utils.c
                      , ./Cdir/halfspaces.c
+-- additions for 0.1.0.4
+  install-includes:   ./Cdir/qhull_ra.h
+                    , ./Cdir/libqhull_r.h
+                    , ./Cdir/user_r.h
+                    , ./Cdir/convexhull.h
+                    , ./Cdir/delaunay.h
+                    , ./Cdir/geom_r.h
+                    , ./Cdir/io_r.h
+                    , ./Cdir/mem_r.h
+                    , ./Cdir/merge_r.h
+                    , ./Cdir/poly_r.h
+                    , ./Cdir/qset_r.h
+                    , ./Cdir/random_r.h
+                    , ./Cdir/stat_r.h
+                    , ./Cdir/user_r.h
+                    , ./Cdir/utils.h
 
   ghc-options:         -O0 -Wall