1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
|
/*<html><pre> -<a href="qh-qhull.htm"
>-------------------------------</a><a name="TOP">-</a>
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; j<numpoints; j++) {
point= points + j*dim;
for (k=dim; k--; ) {
if (j & ( 1 << k))
point[k]= 1.0;
else
point[k]= -1.0;
}
}
} /*.makecube.*/
/*--------------------------------------------------
-adddiamond- add diamond to convex hull
points is numpoints+numnew X dim.
notes:
qh_addpoint() does not make a copy of the point coordinates.
For inside points and some outside points, qh_findbestfacet performs
an exhaustive search for a visible facet. Algorithms that retain
previously constructed hulls should be faster for on-line construction
of the convex hull.
*/
void adddiamond (coordT *points, int numpoints, int numnew, int dim) {
int j,k;
coordT *point;
facetT *facet;
boolT isoutside;
realT bestdist;
for (j= 0; j < numnew ; j++) {
point= points + (numpoints+j)*dim;
if (points == qh first_point) /* in case of 'QRn' */
qh num_points= numpoints+j+1;
/* qh num_points sets the size of the points array. You may
allocate the points elsewhere. If so, qh_addpoint records
the point's address in qh other_points
*/
for (k=dim; k--; ) {
if (j/2 == k)
point[k]= (j & 1) ? 2.0 : -2.0;
else
point[k]= 0.0;
}
facet= qh_findbestfacet (point, !qh_ALL, &bestdist, &isoutside);
if (isoutside) {
if (!qh_addpoint (point, facet, False))
break; /* user requested an early exit with 'TVn' or 'TCn' */
}
printf ("%d vertices and %d facets\n",
qh num_vertices, qh num_facets);
/* qh_produce_output(); */
}
if (qh DOcheckmax)
qh_check_maxout();
else if (qh KEEPnearinside)
qh_nearcoplanar();
} /*.adddiamond.*/
/*--------------------------------------------------
-makeDelaunay- set points for dim-1 Delaunay triangulation of random points
points is numpoints X dim. Each point is projected to a paraboloid.
*/
void makeDelaunay (coordT *points, int numpoints, int dim) {
int j,k, seed;
coordT *point, realr;
seed= time(NULL);
printf ("seed: %d\n", seed);
qh_RANDOMseed_( seed);
for (j=0; j<numpoints; j++) {
point= points + j*dim;
for (k= 0; k < dim-1; k++) {
realr= qh_RANDOMint;
point[k]= 2.0 * realr/(qh_RANDOMmax+1) - 1.0;
}
}
qh_setdelaunay (dim, numpoints, points);
} /*.makeDelaunay.*/
/*--------------------------------------------------
-addDelaunay- add points to dim-1 Delaunay triangulation
points is numpoints+numnew X dim. Each point is projected to a paraboloid.
notes:
qh_addpoint() does not make a copy of the point coordinates.
Since qh_addpoint() is not given a visible facet, it performs a directed
search of all facets. Algorithms that retain previously
constructed hulls may be faster.
*/
void addDelaunay (coordT *points, int numpoints, int numnew, int dim) {
int j,k;
coordT *point, realr;
facetT *facet;
realT bestdist;
boolT isoutside;
for (j= 0; j < numnew ; j++) {
point= points + (numpoints+j)*dim;
if (points == qh first_point) /* in case of 'QRn' */
qh num_points= numpoints+j+1;
/* qh num_points sets the size of the points array. You may
allocate the point elsewhere. If so, qh_addpoint records
the point's address in qh other_points
*/
for (k= 0; k < dim-1; k++) {
realr= qh_RANDOMint;
point[k]= 2.0 * realr/(qh_RANDOMmax+1) - 1.0;
}
qh_setdelaunay (dim, 1, point);
facet= qh_findbestfacet (point, !qh_ALL, &bestdist, &isoutside);
if (isoutside) {
if (!qh_addpoint (point, facet, False))
break; /* user requested an early exit with 'TVn' or 'TCn' */
}
qh_printpoint (stdout, "added point", point);
printf ("%d points, %d extra points, %d vertices, and %d facets in total\n",
qh num_points, qh_setsize (qh other_points),
qh num_vertices, qh num_facets);
/* qh_produce_output(); */
}
if (qh DOcheckmax)
qh_check_maxout();
else if (qh KEEPnearinside)
qh_nearcoplanar();
} /*.addDelaunay.*/
/*--------------------------------------------------
-findDelaunay- find Delaunay triangle for [0.5,0.5,...]
assumes dim < 100
*/
void findDelaunay (int dim) {
int k;
coordT point[ 100];
boolT isoutside;
realT bestdist;
facetT *facet;
vertexT *vertex, **vertexp;
for (k= 0; k < dim-1; k++)
point[k]= 0.5;
qh_setdelaunay (dim, 1, point);
facet= qh_findbestfacet (point, qh_ALL, &bestdist, &isoutside);
FOREACHvertex_(facet->vertices) {
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<numpoints; j++) {
point= points + j*(dim+1);
point[dim]= -1.0; /* offset */
for (k=dim; k--; ) {
if (j & ( 1 << k))
point[k]= 1.0;
else
point[k]= -1.0;
}
}
} /*.makehalf.*/
/*--------------------------------------------------
-addhalf- add halfspaces for a (dim)-d cube to the intersection
points is numpoints+numnew X dim+1
notes:
assumes dim < 100.
For makehalf(), points is the initial set of halfspaces with offsets.
It is transformed by qh_sethalfspace_all into a
(dim)-d set of newpoints. Qhull computed the convex hull of newpoints -
this is equivalent to the halfspace intersection of the
orginal halfspaces.
For addhalf(), the remainder of points stores the transforms of
the added halfspaces. Qhull computes the convex hull of newpoints
and the added points. qh_addpoint() does not make a copy of these points.
Since halfspace intersection is equivalent to a convex hull,
qh_findbestfacet may perform an exhaustive search
for a visible facet. Algorithms that retain previously constructed
intersections should be faster for on-line construction.
*/
void addhalf (coordT *points, int numpoints, int numnew, int dim, coordT *feasible) {
int j,k;
coordT *point, normal[100], offset, *next;
facetT *facet;
boolT isoutside;
realT bestdist;
for (j= 0; j < numnew ; j++) {
offset= -1.0;
for (k=dim; k--; ) {
if (j/2 == k) {
normal[k]= sqrt (dim); /* to normalize as in makehalf */
if (j & 1)
normal[k]= -normal[k];
}else
normal[k]= 0.0;
}
point= points + (numpoints+j)* (dim+1); /* does not use point[dim] */
qh_sethalfspace (dim, point, &next, normal, &offset, feasible);
facet= qh_findbestfacet (point, !qh_ALL, &bestdist, &isoutside);
if (isoutside) {
if (!qh_addpoint (point, facet, False))
break; /* user requested an early exit with 'TVn' or 'TCn' */
}
qh_printpoint (stdout, "added offset -1 and normal", normal);
printf ("%d points, %d extra points, %d vertices, and %d facets in total\n",
qh num_points, qh_setsize (qh other_points),
qh num_vertices, qh num_facets);
/* qh_produce_output(); */
}
if (qh DOcheckmax)
qh_check_maxout();
else if (qh KEEPnearinside)
qh_nearcoplanar();
} /*.addhalf.*/
#define DIM 3 /* dimension of points, must be < 31 for SIZEcube */
#define SIZEcube (1<<DIM)
#define SIZEdiamond (2*DIM)
#define TOTpoints (SIZEcube + SIZEdiamond)
/*--------------------------------------------------
-main- derived from call_qhull in user.c
see program header
this contains three runs of Qhull for convex hull, Delaunay
triangulation or Voronoi vertices, and halfspace intersection
*/
int main (int argc, char *argv[]) {
boolT ismalloc;
int curlong, totlong, exitcode;
char options [2000];
printf ("This is the output from user_eg2.c\n\n\
It shows how qhull() may be called from an application. It is not part\n\
of qhull itself. If it appears accidently, please remove user_eg2.c from\n\
your project.\n\n");
ismalloc= False; /* True if qh_freeqhull should 'free(array)' */
/*
Run 1: convex hull
*/
qh_init_A (stdin, stdout, stderr, 0, NULL);
exitcode= setjmp (qh errexit);
if (!exitcode) {
coordT array[TOTpoints][DIM];
strcat (qh rbox_command, "user_eg cube");
sprintf (options, "qhull s Tcv Q11 %s ", argc >= 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
|
