NDOF support added to trunk from ndof branch.
[blender.git] / extern / qhull / src / user_eg2.c
1 /*<html><pre>  -<a                             href="qh-qhull.htm"
2   >-------------------------------</a><a name="TOP">-</a>
3
4   user_eg2.c
5
6   sample code for calling qhull() from an application.
7
8   See user_eg.c for a simpler method using qh_new_qhull().
9   The method used here and in unix.c gives you additional
10   control over Qhull. 
11   
12   call with:
13
14      user_eg2 "triangulated cube/diamond options" "delaunay options" "halfspace options"
15
16   for example:
17
18      user_eg2                             # return summaries
19
20      user_eg2 "n" "o" "Fp"                # return normals, OFF, points
21
22      user_eg2 "QR0 p" "QR0 v p" "QR0 Fp"  # rotate input and return points
23                                          # 'v' returns Voronoi
24                                          # transform is rotated for halfspaces
25
26    main() makes three runs of qhull.
27
28      1) compute the convex hull of a cube, and incrementally add a diamond
29
30      2a) compute the Delaunay triangulation of random points, and add points.
31
32      2b) find the Delaunay triangle closest to a point.
33
34      3) compute the halfspace intersection of a diamond, and add a cube
35
36  notes:
37  
38    summaries are sent to stderr if other output formats are used
39
40    derived from unix.c and compiled by 'make user_eg2'
41
42    see qhull.h for data structures, macros, and user-callable functions.
43    
44    If you want to control all output to stdio and input to stdin,
45    set the #if below to "1" and delete all lines that contain "io.c".  
46    This prevents the loading of io.o.  Qhull will
47    still write to 'qh ferr' (stderr) for error reporting and tracing.
48
49    Defining #if 1, also prevents user.o from being loaded.
50 */
51
52 #include "qhull_a.h"
53
54 /*-------------------------------------------------
55 -internal function prototypes
56 */
57 void print_summary (void);
58 void makecube (coordT *points, int numpoints, int dim);
59 void adddiamond (coordT *points, int numpoints, int numnew, int dim);
60 void makeDelaunay (coordT *points, int numpoints, int dim);
61 void addDelaunay (coordT *points, int numpoints, int numnew, int dim);
62 void findDelaunay (int dim);
63 void makehalf (coordT *points, int numpoints, int dim);
64 void addhalf (coordT *points, int numpoints, int numnew, int dim, coordT *feasible);
65
66 /*-------------------------------------------------
67 -print_summary()
68 */
69 void print_summary (void) {
70   facetT *facet;
71   int k;
72
73   printf ("\n%d vertices and %d facets with normals:\n", 
74                  qh num_vertices, qh num_facets);
75   FORALLfacets {
76     for (k=0; k < qh hull_dim; k++) 
77       printf ("%6.2g ", facet->normal[k]);
78     printf ("\n");
79   }
80 }
81
82 /*--------------------------------------------------
83 -makecube- set points to vertices of cube
84   points is numpoints X dim
85 */
86 void makecube (coordT *points, int numpoints, int dim) {
87   int j,k;
88   coordT *point;
89
90   for (j=0; j<numpoints; j++) {
91     point= points + j*dim;
92     for (k=dim; k--; ) {
93       if (j & ( 1 << k))
94         point[k]= 1.0;
95       else
96         point[k]= -1.0;
97     }
98   }
99 } /*.makecube.*/
100
101 /*--------------------------------------------------
102 -adddiamond- add diamond to convex hull
103   points is numpoints+numnew X dim.
104   
105 notes:
106   qh_addpoint() does not make a copy of the point coordinates.
107
108   For inside points and some outside points, qh_findbestfacet performs 
109   an exhaustive search for a visible facet.  Algorithms that retain 
110   previously constructed hulls should be faster for on-line construction 
111   of the convex hull.
112 */
113 void adddiamond (coordT *points, int numpoints, int numnew, int dim) {
114   int j,k;
115   coordT *point;
116   facetT *facet;
117   boolT isoutside;
118   realT bestdist;
119
120   for (j= 0; j < numnew ; j++) {
121     point= points + (numpoints+j)*dim;
122     if (points == qh first_point)  /* in case of 'QRn' */
123       qh num_points= numpoints+j+1;
124     /* qh num_points sets the size of the points array.  You may
125        allocate the points elsewhere.  If so, qh_addpoint records
126        the point's address in qh other_points 
127     */
128     for (k=dim; k--; ) {
129       if (j/2 == k)
130         point[k]= (j & 1) ? 2.0 : -2.0;
131       else
132         point[k]= 0.0;
133     }
134     facet= qh_findbestfacet (point, !qh_ALL, &bestdist, &isoutside);
135     if (isoutside) {
136       if (!qh_addpoint (point, facet, False))
137         break;  /* user requested an early exit with 'TVn' or 'TCn' */
138     }
139     printf ("%d vertices and %d facets\n", 
140                  qh num_vertices, qh num_facets);
141     /* qh_produce_output(); */
142   }
143   if (qh DOcheckmax)
144     qh_check_maxout();
145   else if (qh KEEPnearinside)
146     qh_nearcoplanar();
147 } /*.adddiamond.*/
148
149 /*--------------------------------------------------
150 -makeDelaunay- set points for dim-1 Delaunay triangulation of random points
151   points is numpoints X dim.  Each point is projected to a paraboloid.
152 */
153 void makeDelaunay (coordT *points, int numpoints, int dim) {
154   int j,k, seed;
155   coordT *point, realr;
156
157   seed= time(NULL);
158   printf ("seed: %d\n", seed);
159   qh_RANDOMseed_( seed);
160   for (j=0; j<numpoints; j++) {
161     point= points + j*dim;
162     for (k= 0; k < dim-1; k++) {
163       realr= qh_RANDOMint;
164       point[k]= 2.0 * realr/(qh_RANDOMmax+1) - 1.0;
165     }
166   }
167   qh_setdelaunay (dim, numpoints, points);
168 } /*.makeDelaunay.*/
169
170 /*--------------------------------------------------
171 -addDelaunay- add points to dim-1 Delaunay triangulation
172   points is numpoints+numnew X dim.  Each point is projected to a paraboloid.
173 notes:
174   qh_addpoint() does not make a copy of the point coordinates.
175
176   Since qh_addpoint() is not given a visible facet, it performs a directed
177   search of all facets.  Algorithms that retain previously
178   constructed hulls may be faster.
179 */
180 void addDelaunay (coordT *points, int numpoints, int numnew, int dim) {
181   int j,k;
182   coordT *point, realr;
183   facetT *facet;
184   realT bestdist;
185   boolT isoutside;
186
187   for (j= 0; j < numnew ; j++) {
188     point= points + (numpoints+j)*dim;
189     if (points == qh first_point)  /* in case of 'QRn' */
190       qh num_points= numpoints+j+1;  
191     /* qh num_points sets the size of the points array.  You may
192        allocate the point elsewhere.  If so, qh_addpoint records
193        the point's address in qh other_points 
194     */
195     for (k= 0; k < dim-1; k++) {
196       realr= qh_RANDOMint;
197       point[k]= 2.0 * realr/(qh_RANDOMmax+1) - 1.0;
198     }
199     qh_setdelaunay (dim, 1, point);
200     facet= qh_findbestfacet (point, !qh_ALL, &bestdist, &isoutside);
201     if (isoutside) {
202       if (!qh_addpoint (point, facet, False))
203         break;  /* user requested an early exit with 'TVn' or 'TCn' */
204     }
205     qh_printpoint (stdout, "added point", point);
206     printf ("%d points, %d extra points, %d vertices, and %d facets in total\n", 
207                   qh num_points, qh_setsize (qh other_points),
208                   qh num_vertices, qh num_facets);
209     
210     /* qh_produce_output(); */
211   }
212   if (qh DOcheckmax)
213     qh_check_maxout();
214   else if (qh KEEPnearinside)
215     qh_nearcoplanar();
216 } /*.addDelaunay.*/
217
218 /*--------------------------------------------------
219 -findDelaunay- find Delaunay triangle for [0.5,0.5,...]
220   assumes dim < 100
221 */
222 void findDelaunay (int dim) {
223   int k;
224   coordT point[ 100];
225   boolT isoutside;
226   realT bestdist;
227   facetT *facet;
228   vertexT *vertex, **vertexp;
229
230   for (k= 0; k < dim-1; k++) 
231     point[k]= 0.5;
232   qh_setdelaunay (dim, 1, point);
233   facet= qh_findbestfacet (point, qh_ALL, &bestdist, &isoutside);
234   FOREACHvertex_(facet->vertices) {
235     for (k=0; k < dim-1; k++)
236       printf ("%5.2f ", vertex->point[k]);
237     printf ("\n");
238   }
239 } /*.findDelaunay.*/
240
241 /*--------------------------------------------------
242 -makehalf- set points to halfspaces for a (dim)-d diamond
243   points is numpoints X dim+1
244
245   each halfspace consists of dim coefficients followed by an offset
246 */
247 void makehalf (coordT *points, int numpoints, int dim) {
248   int j,k;
249   coordT *point;
250
251   for (j=0; j<numpoints; j++) {
252     point= points + j*(dim+1);
253     point[dim]= -1.0; /* offset */
254     for (k=dim; k--; ) {
255       if (j & ( 1 << k))
256         point[k]= 1.0;
257       else
258         point[k]= -1.0;
259     }
260   }
261 } /*.makehalf.*/
262
263 /*--------------------------------------------------
264 -addhalf- add halfspaces for a (dim)-d cube to the intersection
265   points is numpoints+numnew X dim+1
266 notes:
267   assumes dim < 100. 
268
269   For makehalf(), points is the initial set of halfspaces with offsets.
270   It is transformed by qh_sethalfspace_all into a
271   (dim)-d set of newpoints.  Qhull computed the convex hull of newpoints -
272   this is equivalent to the halfspace intersection of the
273   orginal halfspaces.
274
275   For addhalf(), the remainder of points stores the transforms of
276   the added halfspaces.  Qhull computes the convex hull of newpoints
277   and the added points.  qh_addpoint() does not make a copy of these points.
278
279   Since halfspace intersection is equivalent to a convex hull, 
280   qh_findbestfacet may perform an exhaustive search
281   for a visible facet.  Algorithms that retain previously constructed
282   intersections should be faster for on-line construction.
283 */
284 void addhalf (coordT *points, int numpoints, int numnew, int dim, coordT *feasible) {
285   int j,k;
286   coordT *point, normal[100], offset, *next;
287   facetT *facet;
288   boolT isoutside;
289   realT bestdist;
290
291   for (j= 0; j < numnew ; j++) {
292     offset= -1.0; 
293     for (k=dim; k--; ) {
294       if (j/2 == k) {
295         normal[k]= sqrt (dim);   /* to normalize as in makehalf */
296         if (j & 1)
297           normal[k]= -normal[k];
298       }else
299         normal[k]= 0.0;
300     }
301     point= points + (numpoints+j)* (dim+1);  /* does not use point[dim] */
302     qh_sethalfspace (dim, point, &next, normal, &offset, feasible);
303     facet= qh_findbestfacet (point, !qh_ALL, &bestdist, &isoutside);
304     if (isoutside) {
305       if (!qh_addpoint (point, facet, False))
306         break;  /* user requested an early exit with 'TVn' or 'TCn' */
307     }
308     qh_printpoint (stdout, "added offset -1 and normal", normal);
309     printf ("%d points, %d extra points, %d vertices, and %d facets in total\n", 
310                   qh num_points, qh_setsize (qh other_points),
311                   qh num_vertices, qh num_facets);
312     /* qh_produce_output(); */
313   }
314   if (qh DOcheckmax)
315     qh_check_maxout();
316   else if (qh KEEPnearinside)
317     qh_nearcoplanar();
318 } /*.addhalf.*/
319
320 #define DIM 3     /* dimension of points, must be < 31 for SIZEcube */
321 #define SIZEcube (1<<DIM)
322 #define SIZEdiamond (2*DIM)
323 #define TOTpoints (SIZEcube + SIZEdiamond)
324
325 /*--------------------------------------------------
326 -main- derived from call_qhull in user.c
327
328   see program header
329
330   this contains three runs of Qhull for convex hull, Delaunay
331   triangulation or Voronoi vertices, and halfspace intersection
332
333 */
334 int main (int argc, char *argv[]) {
335   boolT ismalloc;
336   int curlong, totlong, exitcode;
337   char options [2000];
338
339   printf ("This is the output from user_eg2.c\n\n\
340 It shows how qhull() may be called from an application.  It is not part\n\
341 of qhull itself.  If it appears accidently, please remove user_eg2.c from\n\
342 your project.\n\n");
343   ismalloc= False;      /* True if qh_freeqhull should 'free(array)' */
344   /*
345     Run 1: convex hull
346   */
347   qh_init_A (stdin, stdout, stderr, 0, NULL);
348   exitcode= setjmp (qh errexit);
349   if (!exitcode) {
350     coordT array[TOTpoints][DIM];
351
352     strcat (qh rbox_command, "user_eg cube");
353     sprintf (options, "qhull s Tcv Q11 %s ", argc >= 2 ? argv[1] : "");
354     qh_initflags (options);
355     printf( "\ncompute triangulated convex hull of cube after rotating input\n");
356     makecube (array[0], SIZEcube, DIM);
357     qh_init_B (array[0], SIZEcube, DIM, ismalloc);
358     qh_qhull();
359     qh_check_output();
360     qh_triangulate();  /* requires option 'Q11' if want to add points */ 
361     print_summary ();
362     if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
363       qh_check_points ();
364     printf( "\nadd points in a diamond\n");
365     adddiamond (array[0], SIZEcube, SIZEdiamond, DIM);
366     qh_check_output();
367     print_summary (); 
368     qh_produce_output();  /* delete this line to help avoid io.c */
369     if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
370       qh_check_points ();
371   }
372   qh NOerrexit= True;
373   qh_freeqhull (!qh_ALL);
374   qh_memfreeshort (&curlong, &totlong);
375   /*
376     Run 2: Delaunay triangulation
377   */
378   qh_init_A (stdin, stdout, stderr, 0, NULL);
379   exitcode= setjmp (qh errexit);
380   if (!exitcode) {
381     coordT array[TOTpoints][DIM];
382
383     strcat (qh rbox_command, "user_eg Delaunay");
384     sprintf (options, "qhull s d Tcv %s", argc >= 3 ? argv[2] : "");
385     qh_initflags (options);
386     printf( "\ncompute 2-d Delaunay triangulation\n");
387     makeDelaunay (array[0], SIZEcube, DIM);
388     /* Instead of makeDelaunay with qh_setdelaunay, you may
389        produce a 2-d array of points, set DIM to 2, and set 
390        qh PROJECTdelaunay to True.  qh_init_B will call 
391        qh_projectinput to project the points to the paraboloid
392        and add a point "at-infinity".
393     */
394     qh_init_B (array[0], SIZEcube, DIM, ismalloc);
395     qh_qhull();
396     /* If you want Voronoi ('v') without qh_produce_output(), call
397        qh_setvoronoi_all() after qh_qhull() */
398     qh_check_output();
399     print_summary ();
400     qh_produce_output();  /* delete this line to help avoid io.c */
401     if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
402       qh_check_points ();
403     printf( "\nadd points to triangulation\n");
404     addDelaunay (array[0], SIZEcube, SIZEdiamond, DIM); 
405     qh_check_output();
406     qh_produce_output();  /* delete this line to help avoid io.c */
407     if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
408       qh_check_points ();
409     printf( "\nfind Delaunay triangle closest to [0.5, 0.5, ...]\n");
410     findDelaunay (DIM);
411   }
412   qh NOerrexit= True;
413   qh_freeqhull (!qh_ALL);
414   qh_memfreeshort (&curlong, &totlong);
415   /*
416     Run 3: halfspace intersection
417   */
418   qh_init_A (stdin, stdout, stderr, 0, NULL);
419   exitcode= setjmp (qh errexit);
420   if (!exitcode) {
421     coordT array[TOTpoints][DIM+1];  /* +1 for halfspace offset */
422     pointT *points;
423
424     strcat (qh rbox_command, "user_eg halfspaces");
425     sprintf (options, "qhull H0 s Tcv %s", argc >= 4 ? argv[3] : "");
426     qh_initflags (options);
427     printf( "\ncompute halfspace intersection about the origin for a diamond\n");
428     makehalf (array[0], SIZEcube, DIM);
429     qh_setfeasible (DIM); /* from io.c, sets qh feasible_point from 'Hn,n' */
430     /* you may malloc and set qh feasible_point directly.  It is only used for
431        option 'Fp' */
432     points= qh_sethalfspace_all ( DIM+1, SIZEcube, array[0], qh feasible_point); 
433     qh_init_B (points, SIZEcube, DIM, True); /* qh_freeqhull frees points */
434     qh_qhull();
435     qh_check_output();
436     qh_produce_output();  /* delete this line to help avoid io.c */
437     if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
438       qh_check_points ();
439     printf( "\nadd halfspaces for cube to intersection\n");
440     addhalf (array[0], SIZEcube, SIZEdiamond, DIM, qh feasible_point); 
441     qh_check_output();
442     qh_produce_output();  /* delete this line to help avoid io.c */
443     if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
444       qh_check_points ();
445   }
446   qh NOerrexit= True;
447   qh NOerrexit= True;
448   qh_freeqhull (!qh_ALL);
449   qh_memfreeshort (&curlong, &totlong);
450   if (curlong || totlong)  /* could also check previous runs */
451     fprintf (stderr, "qhull internal warning (main): did not free %d bytes of long memory (%d pieces)\n",
452        totlong, curlong);
453   return exitcode;
454 } /* main */
455
456 #if 1    /* use 1 to prevent loading of io.o and user.o */
457 /*-------------------------------------------
458 -errexit- return exitcode to system after an error
459   assumes exitcode non-zero
460   prints useful information
461   see qh_errexit2() in qhull.c for 2 facets
462 */
463 void qh_errexit(int exitcode, facetT *facet, ridgeT *ridge) {
464
465   if (qh ERREXITcalled) {
466     fprintf (qh ferr, "qhull error while processing previous error.  Exit program\n");
467     exit(1);
468   }
469   qh ERREXITcalled= True;
470   if (!qh QHULLfinished)
471     qh hulltime= (unsigned)clock() - qh hulltime;
472   fprintf (qh ferr, "\nWhile executing: %s | %s\n", qh rbox_command, qh qhull_command);
473   fprintf(qh ferr, "Options selected:\n%s\n", qh qhull_options);
474   if (qh furthest_id >= 0) {
475     fprintf(qh ferr, "\nLast point added to hull was p%d", qh furthest_id);
476     if (zzval_(Ztotmerge))
477       fprintf(qh ferr, "  Last merge was #%d.", zzval_(Ztotmerge));
478     if (qh QHULLfinished)
479       fprintf(qh ferr, "\nQhull has finished constructing the hull.");
480     else if (qh POSTmerging)
481       fprintf(qh ferr, "\nQhull has started post-merging");
482     fprintf(qh ferr, "\n\n");
483   }
484   if (qh NOerrexit) {
485     fprintf (qh ferr, "qhull error while ending program.  Exit program\n");
486     exit(1);
487   }
488   if (!exitcode)
489     exitcode= qh_ERRqhull;
490   qh NOerrexit= True;
491   longjmp(qh errexit, exitcode);
492 } /* errexit */
493
494
495 /*-------------------------------------------
496 -errprint- prints out the information of the erroneous object
497     any parameter may be NULL, also prints neighbors and geomview output
498 */
499 void qh_errprint(char *string, facetT *atfacet, facetT *otherfacet, ridgeT *atridge, vertexT *atvertex) {
500
501   fprintf (qh ferr, "%s facets f%d f%d ridge r%d vertex v%d\n",
502            string, getid_(atfacet), getid_(otherfacet), getid_(atridge),
503            getid_(atvertex));
504 } /* errprint */
505
506
507 void qh_printfacetlist(facetT *facetlist, setT *facets, boolT printall) {
508   facetT *facet, **facetp;
509
510   /* remove these calls to help avoid io.c */
511   qh_printbegin (qh ferr, qh_PRINTfacets, facetlist, facets, printall);/*io.c*/
512   FORALLfacet_(facetlist)                                              /*io.c*/
513     qh_printafacet(qh ferr, qh_PRINTfacets, facet, printall);          /*io.c*/
514   FOREACHfacet_(facets)                                                /*io.c*/
515     qh_printafacet(qh ferr, qh_PRINTfacets, facet, printall);          /*io.c*/
516   qh_printend (qh ferr, qh_PRINTfacets, facetlist, facets, printall);  /*io.c*/
517
518   FORALLfacet_(facetlist)
519     fprintf( qh ferr, "facet f%d\n", facet->id);
520 } /* printfacetlist */
521
522
523
524 /*-----------------------------------------
525 -user_memsizes- allocate up to 10 additional, quick allocation sizes
526 */
527 void qh_user_memsizes (void) {
528
529   /* qh_memsize (size); */
530 } /* user_memsizes */
531
532 #endif