/* ---------------------------------
user_eg2.c
sample code for calling qhull() from an application.
See user_eg.c for a simpler method using qh_new_qhull().
The method used here and in unix.c gives you additional
control over Qhull.
call with:
user_eg2 "triangulated cube/diamond options" "delaunay options" "halfspace options"
for example:
user_eg2 # return summaries
user_eg2 "n" "o" "Fp" # return normals, OFF, points
user_eg2 "QR0 p" "QR0 v p" "QR0 Fp" # rotate input and return points
# 'v' returns Voronoi
# transform is rotated for halfspaces
main() makes three runs of qhull.
1) compute the convex hull of a cube, and incrementally add a diamond
2a) compute the Delaunay triangulation of random points, and add points.
2b) find the Delaunay triangle closest to a point.
3) compute the halfspace intersection of a diamond, and add a cube
notes:
summaries are sent to stderr if other output formats are used
derived from unix.c and compiled by 'make user_eg2'
see qhull.h for data structures, macros, and user-callable functions.
If you want to control all output to stdio and input to stdin,
set the #if below to "1" and delete all lines that contain "io.c".
This prevents the loading of io.o. Qhull will
still write to 'qh ferr' (stderr) for error reporting and tracing.
Defining #if 1, also prevents user.o from being loaded.
*/
#include "qhull_a.h"
/*-------------------------------------------------
-internal function prototypes
*/
void print_summary (void);
void makecube (coordT *points, int numpoints, int dim);
void adddiamond (coordT *points, int numpoints, int numnew, int dim);
void makeDelaunay (coordT *points, int numpoints, int dim);
void addDelaunay (coordT *points, int numpoints, int numnew, int dim);
void findDelaunay (int dim);
void makehalf (coordT *points, int numpoints, int dim);
void addhalf (coordT *points, int numpoints, int numnew, int dim, coordT *feasible);
/*-------------------------------------------------
-print_summary()
*/
void print_summary (void) {
facetT *facet;
int k;
printf ("\n%d vertices and %d facets with normals:\n",
qh num_vertices, qh num_facets);
FORALLfacets {
for (k=0; k < qh hull_dim; k++)
printf ("%6.2g ", facet->normal[k]);
printf ("\n");
}
}
/*--------------------------------------------------
-makecube- set points to vertices of cube
points is numpoints X dim
*/
void makecube (coordT *points, int numpoints, int dim) {
int j,k;
coordT *point;
for (j=0; jlocate a facet with qh_findbestfacet()
*/
void findDelaunay (int dim) {
int k;
coordT point[ 100];
boolT isoutside;
realT bestdist;
facetT *facet;
vertexT *vertex, **vertexp;
for (k= 0; k < dim-1; k++)
point[k]= 0.5;
qh_setdelaunay (dim, 1, point);
facet= qh_findbestfacet (point, qh_ALL, &bestdist, &isoutside);
if (facet->tricoplanar) {
fprintf(stderr, "findDelaunay: not implemented for triangulated, non-simplicial Delaunay regions (tricoplanar facet, f%d).\n",
facet->id);
qh_errexit (qh_ERRqhull, facet, NULL);
}
FOREACHvertex_(facet->vertices) {
for (k=0; k < dim-1; k++)
printf ("%5.2f ", vertex->point[k]);
printf ("\n");
}
} /*.findDelaunay.*/
/*--------------------------------------------------
-makehalf- set points to halfspaces for a (dim)-d diamond
points is numpoints X dim+1
each halfspace consists of dim coefficients followed by an offset
*/
void makehalf (coordT *points, int numpoints, int dim) {
int j,k;
coordT *point;
for (j=0; j= 2 ? argv[1] : "");
qh_initflags (options);
printf( "\ncompute triangulated convex hull of cube after rotating input\n");
makecube (array[0], SIZEcube, DIM);
qh_init_B (array[0], SIZEcube, DIM, ismalloc);
qh_qhull();
qh_check_output();
qh_triangulate(); /* requires option 'Q11' if want to add points */
print_summary ();
if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
qh_check_points ();
printf( "\nadd points in a diamond\n");
adddiamond (array[0], SIZEcube, SIZEdiamond, DIM);
qh_check_output();
print_summary ();
qh_produce_output(); /* delete this line to help avoid io.c */
if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
qh_check_points ();
}
qh NOerrexit= True;
qh_freeqhull (!qh_ALL);
qh_memfreeshort (&curlong, &totlong);
/*
Run 2: Delaunay triangulation
*/
qh_init_A (stdin, stdout, stderr, 0, NULL);
exitcode= setjmp (qh errexit);
if (!exitcode) {
coordT array[TOTpoints][DIM];
strcat (qh rbox_command, "user_eg Delaunay");
sprintf (options, "qhull s d Tcv %s", argc >= 3 ? argv[2] : "");
qh_initflags (options);
printf( "\ncompute %d-d Delaunay triangulation\n", DIM-1);
makeDelaunay (array[0], SIZEcube, DIM);
/* Instead of makeDelaunay with qh_setdelaunay, you may
produce a 2-d array of points, set DIM to 2, and set
qh PROJECTdelaunay to True. qh_init_B will call
qh_projectinput to project the points to the paraboloid
and add a point "at-infinity".
*/
qh_init_B (array[0], SIZEcube, DIM, ismalloc);
qh_qhull();
/* If you want Voronoi ('v') without qh_produce_output(), call
qh_setvoronoi_all() after qh_qhull() */
qh_check_output();
print_summary ();
qh_produce_output(); /* delete this line to help avoid io.c */
if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
qh_check_points ();
printf( "\nadd points to triangulation\n");
addDelaunay (array[0], SIZEcube, SIZEdiamond, DIM);
qh_check_output();
qh_produce_output(); /* delete this line to help avoid io.c */
if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
qh_check_points ();
printf( "\nfind Delaunay triangle closest to [0.5, 0.5, ...]\n");
findDelaunay (DIM);
}
qh NOerrexit= True;
qh_freeqhull (!qh_ALL);
qh_memfreeshort (&curlong, &totlong);
/*
Run 3: halfspace intersection
*/
qh_init_A (stdin, stdout, stderr, 0, NULL);
exitcode= setjmp (qh errexit);
if (!exitcode) {
coordT array[TOTpoints][DIM+1]; /* +1 for halfspace offset */
pointT *points;
strcat (qh rbox_command, "user_eg halfspaces");
sprintf (options, "qhull H0 s Tcv %s", argc >= 4 ? argv[3] : "");
qh_initflags (options);
printf( "\ncompute halfspace intersection about the origin for a diamond\n");
makehalf (array[0], SIZEcube, DIM);
qh_setfeasible (DIM); /* from io.c, sets qh feasible_point from 'Hn,n' */
/* you may malloc and set qh feasible_point directly. It is only used for
option 'Fp' */
points= qh_sethalfspace_all ( DIM+1, SIZEcube, array[0], qh feasible_point);
qh_init_B (points, SIZEcube, DIM, True); /* qh_freeqhull frees points */
qh_qhull();
qh_check_output();
qh_produce_output(); /* delete this line to help avoid io.c */
if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
qh_check_points ();
printf( "\nadd halfspaces for cube to intersection\n");
addhalf (array[0], SIZEcube, SIZEdiamond, DIM, qh feasible_point);
qh_check_output();
qh_produce_output(); /* delete this line to help avoid io.c */
if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
qh_check_points ();
}
qh NOerrexit= True;
qh NOerrexit= True;
qh_freeqhull (!qh_ALL);
qh_memfreeshort (&curlong, &totlong);
if (curlong || totlong) /* could also check previous runs */
fprintf (stderr, "qhull internal warning (main): did not free %d bytes of long memory (%d pieces)\n",
totlong, curlong);
return exitcode;
} /* main */
#if 1 /* use 1 to prevent loading of io.o and user.o */
/*-------------------------------------------
-errexit- return exitcode to system after an error
assumes exitcode non-zero
prints useful information
see qh_errexit2() in qhull.c for 2 facets
*/
void qh_errexit(int exitcode, facetT *facet, ridgeT *ridge) {
if (qh ERREXITcalled) {
fprintf (qh ferr, "qhull error while processing previous error. Exit program\n");
exit(1);
}
qh ERREXITcalled= True;
if (!qh QHULLfinished)
qh hulltime= (unsigned)clock() - qh hulltime;
fprintf (qh ferr, "\nWhile executing: %s | %s\n", qh rbox_command, qh qhull_command);
fprintf(qh ferr, "Options selected:\n%s\n", qh qhull_options);
if (qh furthest_id >= 0) {
fprintf(qh ferr, "\nLast point added to hull was p%d", qh furthest_id);
if (zzval_(Ztotmerge))
fprintf(qh ferr, " Last merge was #%d.", zzval_(Ztotmerge));
if (qh QHULLfinished)
fprintf(qh ferr, "\nQhull has finished constructing the hull.");
else if (qh POSTmerging)
fprintf(qh ferr, "\nQhull has started post-merging");
fprintf(qh ferr, "\n\n");
}
if (qh NOerrexit) {
fprintf (qh ferr, "qhull error while ending program. Exit program\n");
exit(1);
}
if (!exitcode)
exitcode= qh_ERRqhull;
qh NOerrexit= True;
longjmp(qh errexit, exitcode);
} /* errexit */
/*-------------------------------------------
-errprint- prints out the information of the erroneous object
any parameter may be NULL, also prints neighbors and geomview output
*/
void qh_errprint(char *string, facetT *atfacet, facetT *otherfacet, ridgeT *atridge, vertexT *atvertex) {
fprintf (qh ferr, "%s facets f%d f%d ridge r%d vertex v%d\n",
string, getid_(atfacet), getid_(otherfacet), getid_(atridge),
getid_(atvertex));
} /* errprint */
void qh_printfacetlist(facetT *facetlist, setT *facets, boolT printall) {
facetT *facet, **facetp;
/* remove these calls to help avoid io.c */
qh_printbegin (qh ferr, qh_PRINTfacets, facetlist, facets, printall);/*io.c*/
FORALLfacet_(facetlist) /*io.c*/
qh_printafacet(qh ferr, qh_PRINTfacets, facet, printall); /*io.c*/
FOREACHfacet_(facets) /*io.c*/
qh_printafacet(qh ferr, qh_PRINTfacets, facet, printall); /*io.c*/
qh_printend (qh ferr, qh_PRINTfacets, facetlist, facets, printall); /*io.c*/
FORALLfacet_(facetlist)
fprintf( qh ferr, "facet f%d\n", facet->id);
} /* printfacetlist */
/*-----------------------------------------
-user_memsizes- allocate up to 10 additional, quick allocation sizes
*/
void qh_user_memsizes (void) {
/* qh_memsize (size); */
} /* user_memsizes */
#endif