/* ---------------------------------
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 */