/*
--------------------------------- user_eg_r.c sample code for calling qhull() from an application. Uses reentrant libqhull_r call with: user_eg "cube/diamond options" "delaunay options" "halfspace options" for example: user_eg # return summaries user_eg "n" "o" "Fp" # return normals, OFF, points user_eg "n Qt" "o" "Fp" # triangulated cube user_eg "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 2a) compute the Delaunay triangulation of random points 2b) find the Delaunay triangle closest to a point. 3) compute the halfspace intersection of a diamond notes: For another example, see main() in unix_r.c and user_eg2_r.c. These examples, call qh_qhull() directly. They allow tighter control on the code loaded with Qhull. For a C++ example, see user_eg3/user_eg3_r.cpp Summaries are sent to stderr if other output formats are used compiled by 'make bin/user_eg' see libqhull_r.h for data structures, macros, and user-callable functions. */ #include "libqhull_r/qhull_ra.h" /*------------------------------------------------- -internal function prototypes */ void print_summary(qhT *qh); void makecube(coordT *points, int numpoints, int dim); void makeDelaunay(qhT *qh, coordT *points, int numpoints, int dim, int seed); void findDelaunay(qhT *qh, int dim); void makehalf(coordT *points, int numpoints, int dim); /*------------------------------------------------- -print_summary(qh) */ void print_summary(qhT *qh) { 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() calls qh_setdelaunay() to project the point to a parabaloid warning: Errors if it finds a tricoplanar facet ('Qt'). The corresponding Delaunay triangle is in the set of tricoplanar facets or one of their neighbors. This search is not implemented here. */ void findDelaunay(qhT *qh, int dim) { int k; coordT point[ 100]; boolT isoutside; realT bestdist; facetT *facet; vertexT *vertex, **vertexp; for (k=0; k < dim; k++) point[k]= 0.5; qh_setdelaunay(qh, dim+1, 1, point); facet= qh_findbestfacet(qh, point, qh_ALL, &bestdist, &isoutside); if (facet->tricoplanar) { fprintf(stderr, "findDelaunay: search not implemented for triangulated, non-simplicial Delaunay regions (tricoplanar facet, f%d).\n", facet->id); qh_errexit(qh, qh_ERRqhull, facet, NULL); } FOREACHvertex_(facet->vertices) { for (k=0; k < dim; k++) printf("%5.2f ", vertex->point[k]); printf("\n"); } } /*.findDelaunay.*/ /*-------------------------------------------------- -makehalf- set points to halfspaces for a (dim)-dimensional 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] : ""); numpoints= SIZEcube; makecube(points, numpoints, DIM); for (i=numpoints; i--; ) rows[i]= points+dim*i; qh_printmatrix(qh, outfile, "input", rows, numpoints, dim); fflush(NULL); exitcode= qh_new_qhull(qh, dim, numpoints, points, ismalloc, flags, outfile, errfile); fflush(NULL); if (!exitcode) { /* if no error */ /* 'qh->facet_list' contains the convex hull */ print_summary(qh); FORALLfacets { /* ... your code ... */ } } #ifdef qh_NOmem qh_freeqhull(qh, qh_ALL); #else qh_freeqhull(qh, !qh_ALL); /* free long memory */ qh_memfreeshort(qh, &curlong, &totlong); /* free short memory and memory allocator */ if (curlong || totlong) fprintf(errfile, "qhull internal warning (user_eg, #1): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong); #endif /* Run 2: Delaunay triangulation, reusing the previous qh/qh_qh */ printf( "\n========\ncompute %d-d Delaunay triangulation\n", dim); sprintf(flags, "qhull s d Tcv %s", argc >= 3 ? argv[2] : ""); numpoints= SIZEcube; makeDelaunay(qh, points, numpoints, dim, (int)time(NULL)); for (i=numpoints; i--; ) rows[i]= points+dim*i; qh_printmatrix(qh, outfile, "input", rows, numpoints, dim); fflush(NULL); exitcode= qh_new_qhull(qh, dim, numpoints, points, ismalloc, flags, outfile, errfile); fflush(NULL); if (!exitcode) { /* if no error */ /* 'qh->facet_list' contains the convex hull */ /* If you want a Voronoi diagram ('v') and do not request output (i.e., outfile=NULL), call qh_setvoronoi_all() after qh_new_qhull(). */ print_summary(qh); FORALLfacets { /* ... your code ... */ } printf( "\nfind %d-d Delaunay triangle or adjacent triangle closest to [0.5, 0.5, ...]\n", dim); exitcode= setjmp(qh->errexit); if (!exitcode) { /* Trap Qhull errors from findDelaunay(). Without the setjmp(), Qhull will exit() after reporting an error */ qh->NOerrexit= False; findDelaunay(qh, DIM); } qh->NOerrexit= True; } { coordT pointsB[DIM*TOTpoints]; /* array of coordinates for each point */ qhT qh_qhB; /* Create a new instance of Qhull (qhB) */ qhT *qhB= &qh_qhB; qh_zero(qhB, errfile); printf( "\n========\nCompute a new triangulation as a separate instance of Qhull\n"); sprintf(flags, "qhull s d Tcv %s", argc >= 3 ? argv[2] : ""); numpoints= SIZEcube; makeDelaunay(qhB, pointsB, numpoints, dim, (int)time(NULL)+1); for (i=numpoints; i--; ) rows[i]= pointsB+dim*i; qh_printmatrix(qhB, outfile, "input", rows, numpoints, dim); fflush(NULL); exitcode= qh_new_qhull(qhB, dim, numpoints, pointsB, ismalloc, flags, outfile, errfile); fflush(NULL); if (!exitcode) print_summary(qhB); printf( "\n========\nFree memory allocated by the new instance of Qhull, and redisplay the old results.\n"); #ifdef qh_NOmem qh_freeqhull(qh, qh_ALL); #else qh_freeqhull(qhB, !qh_ALL); /* free long memory */ qh_memfreeshort(qhB, &curlong, &totlong); /* free short memory and memory allocator */ if (curlong || totlong) fprintf(errfile, "qhull internal warning (user_eg, #4): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong); #endif printf( "\n\n"); print_summary(qh); /* The other instance is unchanged */ /* Exiting the block frees qh_qhB */ } #ifdef qh_NOmem qh_freeqhull(qh, qh_ALL); #else qh_freeqhull(qh, !qh_ALL); /* free long memory */ qh_memfreeshort(qh, &curlong, &totlong); /* free short memory and memory allocator */ if (curlong || totlong) fprintf(errfile, "qhull internal warning (user_eg, #2): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong); #endif /* Run 3: halfspace intersection about the origin */ printf( "\n========\ncompute halfspace intersection about the origin for a diamond\n"); sprintf(flags, "qhull H0 s Tcv %s", argc >= 4 ? argv[3] : "Fp"); numpoints= SIZEcube; makehalf(points, numpoints, dim); for (i=numpoints; i--; ) rows[i]= points+(dim+1)*i; qh_printmatrix(qh, outfile, "input as halfspace coefficients + offsets", rows, numpoints, dim+1); fflush(NULL); /* use qh_sethalfspace_all to transform the halfspaces yourself. If so, set 'qh->feasible_point and do not use option 'Hn,...' [it would retransform the halfspaces] */ exitcode= qh_new_qhull(qh, dim+1, numpoints, points, ismalloc, flags, outfile, errfile); fflush(NULL); if (!exitcode) print_summary(qh); #ifdef qh_NOmem qh_freeqhull(qh, qh_ALL); #else qh_freeqhull(qh, !qh_ALL); qh_memfreeshort(qh, &curlong, &totlong); if (curlong || totlong) /* could also check previous runs */ fprintf(stderr, "qhull internal warning (user_eg, #3): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong); #endif return exitcode; } /* main */