00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047 #include "qhull_a.h"
00048
00049 #ifdef USE_DMALLOC
00050 #include "dmalloc.h"
00051 #endif
00052
00053 char qh_version[] = "user_eg 97/3/31";
00054
00055
00056
00057
00058 void print_summary (void);
00059 void makecube (coordT *points, int numpoints, int dim);
00060 void adddiamond (coordT *points, int numpoints, int numnew, int dim);
00061 void makeDelaunay (coordT *points, int numpoints, int dim);
00062 void addDelaunay (coordT *points, int numpoints, int numnew, int dim);
00063 void findDelaunay (int dim);
00064 void makehalf (coordT *points, int numpoints, int dim);
00065 void addhalf (coordT *points, int numpoints, int numnew, int dim, coordT *feasible);
00066
00067
00068
00069
00070 void print_summary (void) {
00071 facetT *facet;
00072 int k;
00073
00074 printf ("\n%d vertices and %d facets with normals:\n",
00075 qh num_vertices, qh num_facets);
00076 FORALLfacets {
00077 for (k=0; k < qh hull_dim; k++)
00078 printf ("%6.2g ", facet->normal[k]);
00079 printf ("\n");
00080 }
00081 }
00082
00083
00084
00085
00086
00087 void makecube (coordT *points, int numpoints, int dim) {
00088 int j,k;
00089 coordT *point;
00090
00091 for (j=0; j<numpoints; j++) {
00092 point= points + j*dim;
00093 for (k=dim; k--; ) {
00094 if (j & ( 1 << k))
00095 point[k]= 1.0;
00096 else
00097 point[k]= -1.0;
00098 }
00099 }
00100 }
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114 void adddiamond (coordT *points, int numpoints, int numnew, int dim) {
00115 int j,k;
00116 coordT *point;
00117 facetT *facet;
00118 boolT isoutside;
00119 realT bestdist;
00120
00121 for (j= 0; j < numnew ; j++) {
00122 point= points + (numpoints+j)*dim;
00123 if (points == qh first_point)
00124 qh num_points= numpoints+j+1;
00125
00126
00127
00128
00129 for (k=dim; k--; ) {
00130 if (j/2 == k)
00131 point[k]= (j & 1) ? 2.0 : -2.0;
00132 else
00133 point[k]= 0.0;
00134 }
00135 facet= qh_findbestfacet (point, !qh_ALL, &bestdist, &isoutside);
00136 if (isoutside) {
00137 if (!qh_addpoint (point, facet, False))
00138 break;
00139 }
00140 printf ("%d vertices and %d facets\n",
00141 qh num_vertices, qh num_facets);
00142
00143 }
00144 if (qh DOcheckmax)
00145 qh_check_maxout();
00146 else if (qh KEEPnearinside)
00147 qh_nearcoplanar();
00148 }
00149
00150
00151
00152
00153
00154 void makeDelaunay (coordT *points, int numpoints, int dim) {
00155 int j,k, seed;
00156 coordT *point, realr;
00157
00158 seed= time(NULL);
00159 printf ("seed: %d\n", seed);
00160 qh_RANDOMseed_( seed);
00161 for (j=0; j<numpoints; j++) {
00162 point= points + j*dim;
00163 for (k= 0; k < dim-1; k++) {
00164 realr= qh_RANDOMint;
00165 point[k]= 2.0 * realr/(qh_RANDOMmax+1) - 1.0;
00166 }
00167 }
00168 qh_setdelaunay (dim, numpoints, points);
00169 }
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181 void addDelaunay (coordT *points, int numpoints, int numnew, int dim) {
00182 int j,k;
00183 coordT *point, realr;
00184 facetT *facet;
00185 realT bestdist;
00186 boolT isoutside;
00187
00188 for (j= 0; j < numnew ; j++) {
00189 point= points + (numpoints+j)*dim;
00190 if (points == qh first_point)
00191 qh num_points= numpoints+j+1;
00192
00193
00194
00195
00196 for (k= 0; k < dim-1; k++) {
00197 realr= qh_RANDOMint;
00198 point[k]= 2.0 * realr/(qh_RANDOMmax+1) - 1.0;
00199 }
00200 qh_setdelaunay (dim, 1, point);
00201 facet= qh_findbestfacet (point, !qh_ALL, &bestdist, &isoutside);
00202 if (isoutside) {
00203 if (!qh_addpoint (point, facet, False))
00204 break;
00205 }
00206 qh_printpoint (stdout, "added point", point);
00207 printf ("%d points, %d extra points, %d vertices, and %d facets in total\n",
00208 qh num_points, qh_setsize (qh other_points),
00209 qh num_vertices, qh num_facets);
00210
00211
00212 }
00213 if (qh DOcheckmax)
00214 qh_check_maxout();
00215 else if (qh KEEPnearinside)
00216 qh_nearcoplanar();
00217 }
00218
00219
00220
00221
00222
00223 void findDelaunay (int dim) {
00224 int k;
00225 coordT point[ 100];
00226 boolT isoutside;
00227 realT bestdist;
00228 facetT *facet;
00229 vertexT *vertex, **vertexp;
00230
00231 for (k= 0; k < dim-1; k++)
00232 point[k]= 0.5;
00233 qh_setdelaunay (dim, 1, point);
00234 facet= qh_findbestfacet (point, qh_ALL, &bestdist, &isoutside);
00235 FOREACHvertex_(facet->vertices) {
00236 for (k=0; k < dim-1; k++)
00237 printf ("%5.2f ", vertex->point[k]);
00238 printf ("\n");
00239 }
00240 }
00241
00242
00243
00244
00245
00246
00247
00248 void makehalf (coordT *points, int numpoints, int dim) {
00249 int j,k;
00250 coordT *point;
00251
00252 for (j=0; j<numpoints; j++) {
00253 point= points + j*(dim+1);
00254 point[dim]= -1.0;
00255 for (k=dim; k--; ) {
00256 if (j & ( 1 << k))
00257 point[k]= 1.0;
00258 else
00259 point[k]= -1.0;
00260 }
00261 }
00262 }
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285 void addhalf (coordT *points, int numpoints, int numnew, int dim, coordT *feasible) {
00286 int j,k;
00287 coordT *point, normal[100], offset, *next;
00288 facetT *facet;
00289 boolT isoutside;
00290 realT bestdist;
00291
00292 for (j= 0; j < numnew ; j++) {
00293 offset= -1.0;
00294 for (k=dim; k--; ) {
00295 if (j/2 == k) {
00296 normal[k]= sqrt (dim);
00297 if (j & 1)
00298 normal[k]= -normal[k];
00299 }else
00300 normal[k]= 0.0;
00301 }
00302 point= points + (numpoints+j)* (dim+1);
00303 qh_sethalfspace (dim, point, &next, normal, &offset, feasible);
00304 facet= qh_findbestfacet (point, !qh_ALL, &bestdist, &isoutside);
00305 if (isoutside) {
00306 if (!qh_addpoint (point, facet, False))
00307 break;
00308 }
00309 qh_printpoint (stdout, "added offset -1 and normal", normal);
00310 printf ("%d points, %d extra points, %d vertices, and %d facets in total\n",
00311 qh num_points, qh_setsize (qh other_points),
00312 qh num_vertices, qh num_facets);
00313
00314 }
00315 if (qh DOcheckmax)
00316 qh_check_maxout();
00317 else if (qh KEEPnearinside)
00318 qh_nearcoplanar();
00319 }
00320
00321 #define DIM 3
00322 #define SIZEcube (1<<DIM)
00323 #define SIZEdiamond (2*DIM)
00324 #define TOTpoints (SIZEcube + SIZEdiamond)
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335 int main (int argc, char *argv[]) {
00336 boolT ismalloc;
00337 int curlong, totlong, exitcode;
00338 char options [2000];
00339
00340 printf ("This is the output from user_eg.c\n\n\
00341 It shows how qhull() may be called from an application. It is not part\n\
00342 of qhull itself. If it appears accidently, please remove user_eg.c from\n\
00343 your project.\n\n");
00344 ismalloc= False;
00345
00346
00347
00348 qh_init_A (stdin, stdout, stderr, 0, NULL);
00349 exitcode= setjmp (qh errexit);
00350 if (!exitcode) {
00351 coordT array[TOTpoints][DIM];
00352
00353 strcat (qh rbox_command, "user_eg cube");
00354 sprintf (options, "qhull s Tcv %s", argc >= 2 ? argv[1] : "");
00355 qh_initflags (options);
00356 printf( "\ncompute convex hull of cube\n");
00357 makecube (array[0], SIZEcube, DIM);
00358 qh_init_B (array[0], SIZEcube, DIM, ismalloc);
00359 qh_qhull();
00360 qh_check_output();
00361 print_summary ();
00362 if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
00363 qh_check_points ();
00364 printf( "\nadd points in a diamond\n");
00365 adddiamond (array[0], SIZEcube, SIZEdiamond, DIM);
00366 qh_check_output();
00367 print_summary ();
00368 qh_produce_output();
00369 if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
00370 qh_check_points ();
00371 }
00372 qh NOerrexit= True;
00373 qh_freeqhull (!qh_ALL);
00374 qh_memfreeshort (&curlong, &totlong);
00375
00376
00377
00378 qh_init_A (stdin, stdout, stderr, 0, NULL);
00379 exitcode= setjmp (qh errexit);
00380 if (!exitcode) {
00381 coordT array[TOTpoints][DIM];
00382
00383 strcat (qh rbox_command, "user_eg Delaunay");
00384 sprintf (options, "qhull s d Tcv %s", argc >= 3 ? argv[2] : "");
00385 qh_initflags (options);
00386 printf( "\ncompute 2-d Delaunay triangulation\n");
00387 makeDelaunay (array[0], SIZEcube, DIM);
00388
00389
00390
00391
00392
00393
00394 qh_init_B (array[0], SIZEcube, DIM, ismalloc);
00395 qh_qhull();
00396
00397
00398 qh_check_output();
00399 print_summary ();
00400 qh_produce_output();
00401 if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
00402 qh_check_points ();
00403 printf( "\nadd points to triangulation\n");
00404 addDelaunay (array[0], SIZEcube, SIZEdiamond, DIM);
00405 qh_check_output();
00406 qh_produce_output();
00407 if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
00408 qh_check_points ();
00409 printf( "\nfind Delaunay triangle closest to [0.5, 0.5, ...]\n");
00410 findDelaunay (DIM);
00411 }
00412 qh NOerrexit= True;
00413 qh_freeqhull (!qh_ALL);
00414 qh_memfreeshort (&curlong, &totlong);
00415
00416
00417
00418 qh_init_A (stdin, stdout, stderr, 0, NULL);
00419 exitcode= setjmp (qh errexit);
00420 if (!exitcode) {
00421 coordT array[TOTpoints][DIM+1];
00422 pointT *points;
00423
00424 strcat (qh rbox_command, "user_eg halfspaces");
00425 sprintf (options, "qhull H0 s Tcv %s", argc >= 4 ? argv[3] : "");
00426 qh_initflags (options);
00427 printf( "\ncompute halfspace intersection about the origin for a diamond\n");
00428 makehalf (array[0], SIZEcube, DIM);
00429 qh_setfeasible (DIM);
00430
00431
00432 points= qh_sethalfspace_all ( DIM+1, SIZEcube, array[0], qh feasible_point);
00433 qh_init_B (points, SIZEcube, DIM, True);
00434 qh_qhull();
00435 qh_check_output();
00436 qh_produce_output();
00437 if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
00438 qh_check_points ();
00439 printf( "\nadd halfspaces for cube to intersection\n");
00440 addhalf (array[0], SIZEcube, SIZEdiamond, DIM, qh feasible_point);
00441 qh_check_output();
00442 qh_produce_output();
00443 if (qh VERIFYoutput && !qh STOPpoint && !qh STOPcone)
00444 qh_check_points ();
00445 }
00446 qh NOerrexit= True;
00447 qh NOerrexit= True;
00448 qh_freeqhull (!qh_ALL);
00449 qh_memfreeshort (&curlong, &totlong);
00450 if (curlong || totlong)
00451 fprintf (stderr, "qhull internal warning (main): did not free %d bytes of long memory (%d pieces)\n",
00452 totlong, curlong);
00453 return exitcode;
00454 }
00455
00456 #if 1
00457
00458
00459
00460
00461
00462
00463 void qh_errexit(int exitcode, facetT *facet, ridgeT *ridge) {
00464
00465 if (qh ERREXITcalled) {
00466 fprintf (qh ferr, "qhull error while processing previous error. Exit program\n");
00467 exit(1);
00468 }
00469 qh ERREXITcalled= True;
00470 if (!qh QHULLfinished)
00471 qh hulltime= (unsigned)clock() - qh hulltime;
00472 fprintf (qh ferr, "\nWhile executing: %s | %s\n", qh rbox_command, qh qhull_command);
00473 fprintf(qh ferr, "Options selected:\n%s\n", qh qhull_options);
00474 if (qh furthest_id >= 0) {
00475 fprintf(qh ferr, "\nLast point added to hull was p%d", qh furthest_id);
00476 if (zzval_(Ztotmerge))
00477 fprintf(qh ferr, " Last merge was #%d.", zzval_(Ztotmerge));
00478 if (qh QHULLfinished)
00479 fprintf(qh ferr, "\nQhull has finished constructing the hull.");
00480 else if (qh POSTmerging)
00481 fprintf(qh ferr, "\nQhull has started post-merging");
00482 fprintf(qh ferr, "\n\n");
00483 }
00484 if (qh NOerrexit) {
00485 fprintf (qh ferr, "qhull error while ending program. Exit program\n");
00486 exit(1);
00487 }
00488 if (!exitcode)
00489 exitcode= qh_ERRqhull;
00490 qh NOerrexit= True;
00491 longjmp(qh errexit, exitcode);
00492 }
00493
00494
00495
00496
00497
00498
00499 void qh_errprint(char *string, facetT *atfacet, facetT *otherfacet, ridgeT *atridge, vertexT *atvertex) {
00500
00501 fprintf (qh ferr, "%s facets f%d f%d ridge r%d vertex v%d\n",
00502 string, getid_(atfacet), getid_(otherfacet), getid_(atridge),
00503 getid_(atvertex));
00504 }
00505
00506
00507 void qh_printfacetlist(facetT *facetlist, setT *facets, boolT printall) {
00508 facetT *facet, **facetp;
00509
00510
00511 qh_printbegin (qh ferr, qh_PRINTfacets, facetlist, facets, printall);
00512 FORALLfacet_(facetlist)
00513 qh_printafacet(qh ferr, qh_PRINTfacets, facet, printall);
00514 FOREACHfacet_(facets)
00515 qh_printafacet(qh ferr, qh_PRINTfacets, facet, printall);
00516 qh_printend (qh ferr, qh_PRINTfacets, facetlist, facets, printall);
00517
00518 FORALLfacet_(facetlist)
00519 fprintf( qh ferr, "facet f%d\n", facet->id);
00520 }
00521
00522
00523
00524
00525
00526
00527 void qh_user_memsizes (void) {
00528
00529
00530 }
00531
00532 #endif