/*
--------------------------------- 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; jvertices) { 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 2-d Delaunay triangulation\n"); 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