libqhull: user_eg.c Source File

Go to the documentation of this file.
00001 /*<html><pre>  -<a                             href="../libqhull/qh-user.htm"
00002   >-------------------------------</a><a name="TOP">-</a>
00003 
00004   user_eg.c
00005   sample code for calling qhull() from an application
00006 
00007   call with:
00008 
00009      user_eg "cube/diamond options" "delaunay options" "halfspace options"
00010 
00011   for example:
00012 
00013      user_eg                             # return summaries
00014 
00015      user_eg "n" "o" "Fp"                # return normals, OFF, points
00016 
00017      user_eg "n Qt" "o" "Fp"             # triangulated cube
00018 
00019      user_eg "QR0 p" "QR0 v p" "QR0 Fp"  # rotate input and return points
00020                                          # 'v' returns Voronoi
00021                                          # transform is rotated for halfspaces
00022 
00023    main() makes three runs of qhull.
00024 
00025      1) compute the convex hull of a cube
00026 
00027      2a) compute the Delaunay triangulation of random points
00028 
00029      2b) find the Delaunay triangle closest to a point.
00030 
00031      3) compute the halfspace intersection of a diamond
00032 
00033  notes:
00034 
00035    For another example, see main() in unix.c and user_eg2.c.
00036    These examples, call qh_qhull() directly.  They allow
00037    tighter control on the code loaded with Qhull.
00038 
00039    For a C++ example, see user_eg3.cpp
00040 
00041    Summaries are sent to stderr if other output formats are used
00042 
00043    compiled by 'make user_eg'
00044 
00045    see libqhull.h for data structures, macros, and user-callable functions.
00046 */
00047 
00048 #define qh_QHimport
00049 #include "qhull_a.h"
00050 
00051 /*-------------------------------------------------
00052 -internal function prototypes
00053 */
00054 void print_summary (void);
00055 void makecube (coordT *points, int numpoints, int dim);
00056 void makeDelaunay (coordT *points, int numpoints, int dim, int seed);
00057 void findDelaunay (int dim);
00058 void makehalf (coordT *points, int numpoints, int dim);
00059 
00060 /*-------------------------------------------------
00061 -print_summary()
00062 */
00063 void print_summary (void) {
00064   facetT *facet;
00065   int k;
00066 
00067   printf ("\n%d vertices and %d facets with normals:\n",
00068                  qh num_vertices, qh num_facets);
00069   FORALLfacets {
00070     for (k=0; k < qh hull_dim; k++)
00071       printf ("%6.2g ", facet->normal[k]);
00072     printf ("\n");
00073   }
00074 }
00075 
00076 /*--------------------------------------------------
00077 -makecube- set points to vertices of cube
00078   points is numpoints X dim
00079 */
00080 void makecube (coordT *points, int numpoints, int dim) {
00081   int j,k;
00082   coordT *point;
00083 
00084   for (j=0; j<numpoints; j++) {
00085     point= points + j*dim;
00086     for (k=dim; k--; ) {
00087       if (j & ( 1 << k))
00088         point[k]= 1.0;
00089       else
00090         point[k]= -1.0;
00091     }
00092   }
00093 } /*.makecube.*/
00094 
00095 /*--------------------------------------------------
00096 -makeDelaunay- set points for dim Delaunay triangulation of random points
00097   points is numpoints X dim.
00098 notes:
00099   makeDelaunay() in user_eg2.c uses qh_setdelaunay() to project points in place.
00100 */
00101 void makeDelaunay (coordT *points, int numpoints, int dim, int seed) {
00102   int j,k;
00103   coordT *point, realr;
00104 
00105   printf ("seed: %d\n", seed);
00106   qh_RANDOMseed_( seed);
00107   for (j=0; j<numpoints; j++) {
00108     point= points + j*dim;
00109     for (k= 0; k < dim; k++) {
00110       realr= qh_RANDOMint;
00111       point[k]= 2.0 * realr/(qh_RANDOMmax+1) - 1.0;
00112     }
00113   }
00114 } /*.makeDelaunay.*/
00115 
00116 /*--------------------------------------------------
00117 -findDelaunay- find Delaunay triangle for [0.5,0.5,...]
00118   assumes dim < 100
00119 notes:
00120   calls qh_setdelaunay() to project the point to a parabaloid
00121 warning:
00122   This is not implemented for tricoplanar facets ('Qt'),
00123   See <a href="../html/qh-code.htm#findfacet">locate a facet with qh_findbestfacet()</a>
00124 */
00125 void findDelaunay (int dim) {
00126   int k;
00127   coordT point[ 100];
00128   boolT isoutside;
00129   realT bestdist;
00130   facetT *facet;
00131   vertexT *vertex, **vertexp;
00132 
00133   for (k= 0; k < dim; k++)
00134     point[k]= 0.5;
00135   qh_setdelaunay (dim+1, 1, point);
00136   facet= qh_findbestfacet (point, qh_ALL, &bestdist, &isoutside);
00137   if (facet->tricoplanar) {
00138     fprintf(stderr, "findDelaunay: not implemented for triangulated, non-simplicial Delaunay regions (tricoplanar facet, f%d).\n",
00139        facet->id);
00140     qh_errexit (qh_ERRqhull, facet, NULL);
00141   }
00142   FOREACHvertex_(facet->vertices) {
00143     for (k=0; k < dim; k++)
00144       printf ("%5.2f ", vertex->point[k]);
00145     printf ("\n");
00146   }
00147 } /*.findDelaunay.*/
00148 
00149 /*--------------------------------------------------
00150 -makehalf- set points to halfspaces for a (dim)-dimensional diamond
00151   points is numpoints X dim+1
00152 
00153   each halfspace consists of dim coefficients followed by an offset
00154 */
00155 void makehalf (coordT *points, int numpoints, int dim) {
00156   int j,k;
00157   coordT *point;
00158 
00159   for (j=0; j<numpoints; j++) {
00160     point= points + j*(dim+1);
00161     point[dim]= -1.0; /* offset */
00162     for (k=dim; k--; ) {
00163       if (j & ( 1 << k))
00164         point[k]= 1.0;
00165       else
00166         point[k]= -1.0;
00167     }
00168   }
00169 } /*.makehalf.*/
00170 
00171 #define DIM 3     /* dimension of points, must be < 31 for SIZEcube */
00172 #define SIZEcube (1<<DIM)
00173 #define SIZEdiamond (2*DIM)
00174 #define TOTpoints (SIZEcube + SIZEdiamond)
00175 
00176 /*--------------------------------------------------
00177 -main- derived from call_qhull in user.c
00178 
00179   see program header
00180 
00181   this contains three runs of Qhull for convex hull, Delaunay
00182   triangulation or Voronoi vertices, and halfspace intersection
00183 
00184 */
00185 int main (int argc, char *argv[]) {
00186   int dim= DIM;             /* dimension of points */
00187   int numpoints;            /* number of points */
00188   coordT points[(DIM+1)*TOTpoints]; /* array of coordinates for each point */
00189   coordT *rows[TOTpoints];
00190   boolT ismalloc= False;    /* True if qhull should free points in qh_freeqhull() or reallocation */
00191   char flags[250];          /* option flags for qhull, see qh-quick.htm */
00192   FILE *outfile= stdout;    /* output from qh_produce_output()
00193                                use NULL to skip qh_produce_output() */
00194   FILE *errfile= stderr;    /* error messages from qhull code */
00195   int exitcode;             /* 0 if no error from qhull */
00196   facetT *facet;            /* set by FORALLfacets */
00197   int curlong, totlong;     /* memory remaining after qh_memfreeshort */
00198   int i;
00199 
00200   printf ("This is the output from user_eg.c\n\n\
00201 It shows how qhull() may be called from an application using the qhull\n\
00202 shared library.  It is not part of qhull itself.  If it appears accidently,\n\
00203 please remove user_eg.c from your project.\n\n");
00204 
00205 #if qh_QHpointer  /* see user.h */
00206   if (qh_qh){
00207       printf ("QH6233: Qhull link error.  The global variable qh_qh was not initialized\n\
00208 to NULL by global.c.  Please compile user_eg.c with -Dqh_QHpointer_dllimport\n\
00209 as well as -Dqh_QHpointer, or use libqhullstatic, or use a different tool chain.\n\n");
00210       return -1;
00211   }
00212 #endif
00213 
00214   /*
00215     Run 1: convex hull
00216   */
00217   printf( "\ncompute convex hull of cube after rotating input\n");
00218   sprintf (flags, "qhull s Tcv %s", argc >= 2 ? argv[1] : "");
00219   numpoints= SIZEcube;
00220   makecube (points, numpoints, DIM);
00221   for (i=numpoints; i--; )
00222     rows[i]= points+dim*i;
00223   qh_printmatrix (outfile, "input", rows, numpoints, dim);
00224   exitcode= qh_new_qhull (dim, numpoints, points, ismalloc,
00225                       flags, outfile, errfile);
00226   if (!exitcode) {                  /* if no error */
00227     /* 'qh facet_list' contains the convex hull */
00228     print_summary();
00229     FORALLfacets {
00230        /* ... your code ... */
00231     }
00232   }
00233   qh_freeqhull(!qh_ALL);                   /* free long memory  */
00234   qh_memfreeshort (&curlong, &totlong);    /* free short memory and memory allocator */
00235   if (curlong || totlong)
00236     fprintf (errfile, "qhull internal warning (user_eg, #1): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
00237 
00238   /*
00239     Run 2: Delaunay triangulation
00240   */
00241 
00242   printf( "\ncompute %d-d Delaunay triangulation\n", dim);
00243   sprintf (flags, "qhull s d Tcv %s", argc >= 3 ? argv[2] : "");
00244   numpoints= SIZEcube;
00245   makeDelaunay (points, numpoints, dim, (int)time(NULL));
00246   for (i=numpoints; i--; )
00247     rows[i]= points+dim*i;
00248   qh_printmatrix (outfile, "input", rows, numpoints, dim);
00249   exitcode= qh_new_qhull (dim, numpoints, points, ismalloc,
00250                       flags, outfile, errfile);
00251   if (!exitcode) {                  /* if no error */
00252     /* 'qh facet_list' contains the convex hull */
00253     /* If you want a Voronoi diagram ('v') and do not request output (i.e., outfile=NULL),
00254        call qh_setvoronoi_all() after qh_new_qhull(). */
00255     print_summary();
00256     FORALLfacets {
00257        /* ... your code ... */
00258     }
00259     printf( "\nfind %d-d Delaunay triangle closest to [0.5, 0.5, ...]\n", dim);
00260     exitcode= setjmp (qh errexit);
00261     if (!exitcode) {
00262       /* Trap Qhull errors in findDelaunay().  Without the setjmp(), Qhull
00263          will exit() after reporting an error */
00264       qh NOerrexit= False;
00265       findDelaunay (DIM);
00266     }
00267     qh NOerrexit= True;
00268   }
00269 #if qh_QHpointer  /* see user.h */
00270   {
00271     qhT *oldqhA, *oldqhB;
00272     coordT pointsB[DIM*TOTpoints]; /* array of coordinates for each point */
00273 
00274 
00275     printf( "\nsave first triangulation and compute a new triangulation\n");
00276     oldqhA= qh_save_qhull();
00277     sprintf (flags, "qhull s d Tcv %s", argc >= 3 ? argv[2] : "");
00278     numpoints= SIZEcube;
00279     makeDelaunay (pointsB, numpoints, dim, (int)time(NULL)+1);
00280     for (i=numpoints; i--; )
00281       rows[i]= pointsB+dim*i;
00282     qh_printmatrix (outfile, "input", rows, numpoints, dim);
00283     exitcode= qh_new_qhull (dim, numpoints, pointsB, ismalloc,
00284                       flags, outfile, errfile);
00285     if (!exitcode)
00286       print_summary();
00287     printf( "\nsave second triangulation and restore first one\n");
00288     oldqhB= qh_save_qhull();
00289     qh_restore_qhull (&oldqhA);
00290     print_summary();
00291     printf( "\nfree first triangulation and restore second one.\n");
00292     qh_freeqhull (qh_ALL);               /* free short and long memory used by first call */
00293                                          /* do not use qh_memfreeshort */
00294     qh_restore_qhull (&oldqhB);
00295     print_summary();
00296   }
00297 #endif
00298   qh_freeqhull(!qh_ALL);                 /* free long memory */
00299   qh_memfreeshort (&curlong, &totlong);  /* free short memory and memory allocator */
00300   if (curlong || totlong)
00301     fprintf (errfile, "qhull internal warning (user_eg, #2): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
00302 
00303   /*
00304     Run 3: halfspace intersection about the origin
00305   */
00306   printf( "\ncompute halfspace intersection about the origin for a diamond\n");
00307   sprintf (flags, "qhull H0 s Tcv %s", argc >= 4 ? argv[3] : "Fp");
00308   numpoints= SIZEcube;
00309   makehalf (points, numpoints, dim);
00310   for (i=numpoints; i--; )
00311     rows[i]= points+(dim+1)*i;
00312   qh_printmatrix (outfile, "input as halfspace coefficients + offsets", rows, numpoints, dim+1);
00313   /* use qh_sethalfspace_all to transform the halfspaces yourself.
00314      If so, set 'qh feasible_point and do not use option 'Hn,...' [it would retransform the halfspaces]
00315   */
00316   exitcode= qh_new_qhull (dim+1, numpoints, points, ismalloc,
00317                       flags, outfile, errfile);
00318   if (!exitcode)
00319     print_summary();
00320   qh_freeqhull (!qh_ALL);
00321   qh_memfreeshort (&curlong, &totlong);
00322   if (curlong || totlong)  /* could also check previous runs */
00323     fprintf (stderr, "qhull internal warning (user_eg, #3): did not free %d bytes of long memory (%d pieces)\n",
00324        totlong, curlong);
00325   return exitcode;
00326 } /* main */
00327