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