00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "qhull_a.h"
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 void qh_produce_output(void) {
00040 int tempsize= qh_setsize(qhmem.tempstack);
00041
00042 qh_prepare_output();
00043 qh_produce_output2();
00044 if (qh_setsize(qhmem.tempstack) != tempsize) {
00045 qh_fprintf(qh ferr, 6206, "qhull internal error (qh_produce_output): temporary sets not empty(%d)\n",
00046 qh_setsize(qhmem.tempstack));
00047 qh_errexit(qh_ERRqhull, NULL, NULL);
00048 }
00049 }
00050
00051
00052 void qh_produce_output2(void) {
00053 int i, tempsize= qh_setsize(qhmem.tempstack), d_1;
00054
00055 if (qh PRINTsummary)
00056 qh_printsummary(qh ferr);
00057 else if (qh PRINTout[0] == qh_PRINTnone)
00058 qh_printsummary(qh fout);
00059 for (i=0; i < qh_PRINTEND; i++)
00060 qh_printfacets(qh fout, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
00061 qh_allstatistics();
00062 if (qh PRINTprecision && !qh MERGING && (qh JOGGLEmax > REALmax/2 || qh RERUN))
00063 qh_printstats(qh ferr, qhstat precision, NULL);
00064 if (qh VERIFYoutput && (zzval_(Zridge) > 0 || zzval_(Zridgemid) > 0))
00065 qh_printstats(qh ferr, qhstat vridges, NULL);
00066 if (qh PRINTstatistics) {
00067 qh_printstatistics(qh ferr, "");
00068 qh_memstatistics(qh ferr);
00069 d_1= sizeof(setT) + (qh hull_dim - 1) * SETelemsize;
00070 qh_fprintf(qh ferr, 8040, "\
00071 size in bytes: merge %d ridge %d vertex %d facet %d\n\
00072 normal %d ridge vertices %d facet vertices or neighbors %d\n",
00073 (int)sizeof(mergeT), (int)sizeof(ridgeT),
00074 (int)sizeof(vertexT), (int)sizeof(facetT),
00075 qh normal_size, d_1, d_1 + SETelemsize);
00076 }
00077 if (qh_setsize(qhmem.tempstack) != tempsize) {
00078 qh_fprintf(qh ferr, 6065, "qhull internal error (qh_produce_output2): temporary sets not empty(%d)\n",
00079 qh_setsize(qhmem.tempstack));
00080 qh_errexit(qh_ERRqhull, NULL, NULL);
00081 }
00082 }
00083
00084
00085
00086
00087
00088
00089
00090
00091 void dfacet(unsigned id) {
00092 facetT *facet;
00093
00094 FORALLfacets {
00095 if (facet->id == id) {
00096 qh_printfacet(qh fout, facet);
00097 break;
00098 }
00099 }
00100 }
00101
00102
00103
00104
00105
00106
00107
00108
00109 void dvertex(unsigned id) {
00110 vertexT *vertex;
00111
00112 FORALLvertices {
00113 if (vertex->id == id) {
00114 qh_printvertex(qh fout, vertex);
00115 break;
00116 }
00117 }
00118 }
00119
00120
00121
00122
00123
00124
00125
00126
00127 int qh_compare_vertexpoint(const void *p1, const void *p2) {
00128 const vertexT *a= *((vertexT *const*)p1), *b= *((vertexT *const*)p2);
00129
00130 return((qh_pointid(a->point) > qh_pointid(b->point)?1:-1));
00131 }
00132
00133
00134
00135
00136
00137
00138
00139 int qh_compare_facetarea(const void *p1, const void *p2) {
00140 const facetT *a= *((facetT *const*)p1), *b= *((facetT *const*)p2);
00141
00142 if (!a->isarea)
00143 return -1;
00144 if (!b->isarea)
00145 return 1;
00146 if (a->f.area > b->f.area)
00147 return 1;
00148 else if (a->f.area == b->f.area)
00149 return 0;
00150 return -1;
00151 }
00152
00153
00154
00155
00156
00157
00158
00159 int qh_compare_facetmerge(const void *p1, const void *p2) {
00160 const facetT *a= *((facetT *const*)p1), *b= *((facetT *const*)p2);
00161
00162 return(a->nummerge - b->nummerge);
00163 }
00164
00165
00166
00167
00168
00169
00170
00171 int qh_compare_facetvisit(const void *p1, const void *p2) {
00172 const facetT *a= *((facetT *const*)p1), *b= *((facetT *const*)p2);
00173 int i,j;
00174
00175 if (!(i= a->visitid))
00176 i= 0 - a->id;
00177 if (!(j= b->visitid))
00178 j= 0 - b->id;
00179 return(i - j);
00180 }
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191 void qh_copyfilename(char *filename, int size, const char* source, int length) {
00192 char c= *source;
00193
00194 if (length > size + 1) {
00195 qh_fprintf(qh ferr, 6040, "qhull error: filename is more than %d characters, %s\n", size-1, source);
00196 qh_errexit(qh_ERRinput, NULL, NULL);
00197 }
00198 strncpy(filename, source, length);
00199 filename[length]= '\0';
00200 if (c == '\'' || c == '"') {
00201 char *s= filename + 1;
00202 char *t= filename;
00203 while (*s) {
00204 if (*s == c) {
00205 if (s[-1] == '\\')
00206 t[-1]= c;
00207 }else
00208 *t++= *s;
00209 s++;
00210 }
00211 *t= '\0';
00212 }
00213 }
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242 void qh_countfacets(facetT *facetlist, setT *facets, boolT printall,
00243 int *numfacetsp, int *numsimplicialp, int *totneighborsp, int *numridgesp, int *numcoplanarsp, int *numtricoplanarsp) {
00244 facetT *facet, **facetp;
00245 int numfacets= 0, numsimplicial= 0, numridges= 0, totneighbors= 0, numcoplanars= 0, numtricoplanars= 0;
00246
00247 FORALLfacet_(facetlist) {
00248 if ((facet->visible && qh NEWfacets)
00249 || (!printall && qh_skipfacet(facet)))
00250 facet->visitid= 0;
00251 else {
00252 facet->visitid= ++numfacets;
00253 totneighbors += qh_setsize(facet->neighbors);
00254 if (facet->simplicial) {
00255 numsimplicial++;
00256 if (facet->keepcentrum && facet->tricoplanar)
00257 numtricoplanars++;
00258 }else
00259 numridges += qh_setsize(facet->ridges);
00260 if (facet->coplanarset)
00261 numcoplanars += qh_setsize(facet->coplanarset);
00262 }
00263 }
00264
00265 FOREACHfacet_(facets) {
00266 if ((facet->visible && qh NEWfacets)
00267 || (!printall && qh_skipfacet(facet)))
00268 facet->visitid= 0;
00269 else {
00270 facet->visitid= ++numfacets;
00271 totneighbors += qh_setsize(facet->neighbors);
00272 if (facet->simplicial){
00273 numsimplicial++;
00274 if (facet->keepcentrum && facet->tricoplanar)
00275 numtricoplanars++;
00276 }else
00277 numridges += qh_setsize(facet->ridges);
00278 if (facet->coplanarset)
00279 numcoplanars += qh_setsize(facet->coplanarset);
00280 }
00281 }
00282 qh visit_id += numfacets+1;
00283 *numfacetsp= numfacets;
00284 *numsimplicialp= numsimplicial;
00285 *totneighborsp= totneighbors;
00286 *numridgesp= numridges;
00287 *numcoplanarsp= numcoplanars;
00288 *numtricoplanarsp= numtricoplanars;
00289 }
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332 pointT *qh_detvnorm(vertexT *vertex, vertexT *vertexA, setT *centers, realT *offsetp) {
00333 facetT *facet, **facetp;
00334 int i, k, pointid, pointidA, point_i, point_n;
00335 setT *simplex= NULL;
00336 pointT *point, **pointp, *point0, *midpoint, *normal, *inpoint;
00337 coordT *coord, *gmcoord, *normalp;
00338 setT *points= qh_settemp(qh TEMPsize);
00339 boolT nearzero= False;
00340 boolT unbounded= False;
00341 int numcenters= 0;
00342 int dim= qh hull_dim - 1;
00343 realT dist, offset, angle, zero= 0.0;
00344
00345 midpoint= qh gm_matrix + qh hull_dim * qh hull_dim;
00346 for (k=0; k < dim; k++)
00347 midpoint[k]= (vertex->point[k] + vertexA->point[k])/2;
00348 FOREACHfacet_(centers) {
00349 numcenters++;
00350 if (!facet->visitid)
00351 unbounded= True;
00352 else {
00353 if (!facet->center)
00354 facet->center= qh_facetcenter(facet->vertices);
00355 qh_setappend(&points, facet->center);
00356 }
00357 }
00358 if (numcenters > dim) {
00359 simplex= qh_settemp(qh TEMPsize);
00360 qh_setappend(&simplex, vertex->point);
00361 if (unbounded)
00362 qh_setappend(&simplex, midpoint);
00363 qh_maxsimplex(dim, points, NULL, 0, &simplex);
00364 qh_setdelnth(simplex, 0);
00365 }else if (numcenters == dim) {
00366 if (unbounded)
00367 qh_setappend(&points, midpoint);
00368 simplex= points;
00369 }else {
00370 qh_fprintf(qh ferr, 6216, "qhull internal error (qh_detvnorm): too few points(%d) to compute separating plane\n", numcenters);
00371 qh_errexit(qh_ERRqhull, NULL, NULL);
00372 }
00373 i= 0;
00374 gmcoord= qh gm_matrix;
00375 point0= SETfirstt_(simplex, pointT);
00376 FOREACHpoint_(simplex) {
00377 if (qh IStracing >= 4)
00378 qh_printmatrix(qh ferr, "qh_detvnorm: Voronoi vertex or midpoint",
00379 &point, 1, dim);
00380 if (point != point0) {
00381 qh gm_row[i++]= gmcoord;
00382 coord= point0;
00383 for (k=dim; k--; )
00384 *(gmcoord++)= *point++ - *coord++;
00385 }
00386 }
00387 qh gm_row[i]= gmcoord;
00388 normal= gmcoord;
00389 qh_sethyperplane_gauss(dim, qh gm_row, point0, True,
00390 normal, &offset, &nearzero);
00391 if (qh GOODvertexp == vertexA->point)
00392 inpoint= vertexA->point;
00393 else
00394 inpoint= vertex->point;
00395 zinc_(Zdistio);
00396 dist= qh_distnorm(dim, inpoint, normal, &offset);
00397 if (dist > 0) {
00398 offset= -offset;
00399 normalp= normal;
00400 for (k=dim; k--; ) {
00401 *normalp= -(*normalp);
00402 normalp++;
00403 }
00404 }
00405 if (qh VERIFYoutput || qh PRINTstatistics) {
00406 pointid= qh_pointid(vertex->point);
00407 pointidA= qh_pointid(vertexA->point);
00408 if (!unbounded) {
00409 zinc_(Zdiststat);
00410 dist= qh_distnorm(dim, midpoint, normal, &offset);
00411 if (dist < 0)
00412 dist= -dist;
00413 zzinc_(Zridgemid);
00414 wwmax_(Wridgemidmax, dist);
00415 wwadd_(Wridgemid, dist);
00416 trace4((qh ferr, 4014, "qh_detvnorm: points %d %d midpoint dist %2.2g\n",
00417 pointid, pointidA, dist));
00418 for (k=0; k < dim; k++)
00419 midpoint[k]= vertexA->point[k] - vertex->point[k];
00420 qh_normalize(midpoint, dim, False);
00421 angle= qh_distnorm(dim, midpoint, normal, &zero);
00422 if (angle < 0.0)
00423 angle= angle + 1.0;
00424 else
00425 angle= angle - 1.0;
00426 if (angle < 0.0)
00427 angle -= angle;
00428 trace4((qh ferr, 4015, "qh_detvnorm: points %d %d angle %2.2g nearzero %d\n",
00429 pointid, pointidA, angle, nearzero));
00430 if (nearzero) {
00431 zzinc_(Zridge0);
00432 wwmax_(Wridge0max, angle);
00433 wwadd_(Wridge0, angle);
00434 }else {
00435 zzinc_(Zridgeok)
00436 wwmax_(Wridgeokmax, angle);
00437 wwadd_(Wridgeok, angle);
00438 }
00439 }
00440 if (simplex != points) {
00441 FOREACHpoint_i_(points) {
00442 if (!qh_setin(simplex, point)) {
00443 facet= SETelemt_(centers, point_i, facetT);
00444 zinc_(Zdiststat);
00445 dist= qh_distnorm(dim, point, normal, &offset);
00446 if (dist < 0)
00447 dist= -dist;
00448 zzinc_(Zridge);
00449 wwmax_(Wridgemax, dist);
00450 wwadd_(Wridge, dist);
00451 trace4((qh ferr, 4016, "qh_detvnorm: points %d %d Voronoi vertex %d dist %2.2g\n",
00452 pointid, pointidA, facet->visitid, dist));
00453 }
00454 }
00455 }
00456 }
00457 *offsetp= offset;
00458 if (simplex != points)
00459 qh_settempfree(&simplex);
00460 qh_settempfree(&points);
00461 return normal;
00462 }
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475 setT *qh_detvridge(vertexT *vertex) {
00476 setT *centers= qh_settemp(qh TEMPsize);
00477 setT *tricenters= qh_settemp(qh TEMPsize);
00478 facetT *neighbor, **neighborp;
00479 boolT firstinf= True;
00480
00481 FOREACHneighbor_(vertex) {
00482 if (neighbor->seen) {
00483 if (neighbor->visitid) {
00484 if (!neighbor->tricoplanar || qh_setunique(&tricenters, neighbor->center))
00485 qh_setappend(¢ers, neighbor);
00486 }else if (firstinf) {
00487 firstinf= False;
00488 qh_setappend(¢ers, neighbor);
00489 }
00490 }
00491 }
00492 qsort(SETaddr_(centers, facetT), (size_t)qh_setsize(centers),
00493 sizeof(facetT *), qh_compare_facetvisit);
00494 qh_settempfree(&tricenters);
00495 return centers;
00496 }
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517 setT *qh_detvridge3 (vertexT *atvertex, vertexT *vertex) {
00518 setT *centers= qh_settemp(qh TEMPsize);
00519 setT *tricenters= qh_settemp(qh TEMPsize);
00520 facetT *neighbor, **neighborp, *facet= NULL;
00521 boolT firstinf= True;
00522
00523 FOREACHneighbor_(atvertex)
00524 neighbor->seen2= False;
00525 FOREACHneighbor_(vertex) {
00526 if (!neighbor->seen2) {
00527 facet= neighbor;
00528 break;
00529 }
00530 }
00531 while (facet) {
00532 facet->seen2= True;
00533 if (neighbor->seen) {
00534 if (facet->visitid) {
00535 if (!facet->tricoplanar || qh_setunique(&tricenters, facet->center))
00536 qh_setappend(¢ers, facet);
00537 }else if (firstinf) {
00538 firstinf= False;
00539 qh_setappend(¢ers, facet);
00540 }
00541 }
00542 FOREACHneighbor_(facet) {
00543 if (!neighbor->seen2) {
00544 if (qh_setin(vertex->neighbors, neighbor))
00545 break;
00546 else
00547 neighbor->seen2= True;
00548 }
00549 }
00550 facet= neighbor;
00551 }
00552 if (qh CHECKfrequently) {
00553 FOREACHneighbor_(vertex) {
00554 if (!neighbor->seen2) {
00555 qh_fprintf(qh ferr, 6217, "qhull internal error (qh_detvridge3): neighbors of vertex p%d are not connected at facet %d\n",
00556 qh_pointid(vertex->point), neighbor->id);
00557 qh_errexit(qh_ERRqhull, neighbor, NULL);
00558 }
00559 }
00560 }
00561 FOREACHneighbor_(atvertex)
00562 neighbor->seen2= True;
00563 qh_settempfree(&tricenters);
00564 return centers;
00565 }
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616 int qh_eachvoronoi(FILE *fp, printvridgeT printvridge, vertexT *atvertex, boolT visitall, qh_RIDGE innerouter, boolT inorder) {
00617 boolT unbounded;
00618 int count;
00619 facetT *neighbor, **neighborp, *neighborA, **neighborAp;
00620 setT *centers;
00621 setT *tricenters= qh_settemp(qh TEMPsize);
00622
00623 vertexT *vertex, **vertexp;
00624 boolT firstinf;
00625 unsigned int numfacets= (unsigned int)qh num_facets;
00626 int totridges= 0;
00627
00628 qh vertex_visit++;
00629 atvertex->seen= True;
00630 if (visitall) {
00631 FORALLvertices
00632 vertex->seen= False;
00633 }
00634 FOREACHneighbor_(atvertex) {
00635 if (neighbor->visitid < numfacets)
00636 neighbor->seen= True;
00637 }
00638 FOREACHneighbor_(atvertex) {
00639 if (neighbor->seen) {
00640 FOREACHvertex_(neighbor->vertices) {
00641 if (vertex->visitid != qh vertex_visit && !vertex->seen) {
00642 vertex->visitid= qh vertex_visit;
00643 count= 0;
00644 firstinf= True;
00645 qh_settruncate(tricenters, 0);
00646 FOREACHneighborA_(vertex) {
00647 if (neighborA->seen) {
00648 if (neighborA->visitid) {
00649 if (!neighborA->tricoplanar || qh_setunique(&tricenters, neighborA->center))
00650 count++;
00651 }else if (firstinf) {
00652 count++;
00653 firstinf= False;
00654 }
00655 }
00656 }
00657 if (count >= qh hull_dim - 1) {
00658 if (firstinf) {
00659 if (innerouter == qh_RIDGEouter)
00660 continue;
00661 unbounded= False;
00662 }else {
00663 if (innerouter == qh_RIDGEinner)
00664 continue;
00665 unbounded= True;
00666 }
00667 totridges++;
00668 trace4((qh ferr, 4017, "qh_eachvoronoi: Voronoi ridge of %d vertices between sites %d and %d\n",
00669 count, qh_pointid(atvertex->point), qh_pointid(vertex->point)));
00670 if (printvridge && fp) {
00671 if (inorder && qh hull_dim == 3+1)
00672 centers= qh_detvridge3 (atvertex, vertex);
00673 else
00674 centers= qh_detvridge(vertex);
00675 (*printvridge) (fp, atvertex, vertex, centers, unbounded);
00676 qh_settempfree(¢ers);
00677 }
00678 }
00679 }
00680 }
00681 }
00682 }
00683 FOREACHneighbor_(atvertex)
00684 neighbor->seen= False;
00685 qh_settempfree(&tricenters);
00686 return totridges;
00687 }
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718 int qh_eachvoronoi_all(FILE *fp, printvridgeT printvridge, boolT isUpper, qh_RIDGE innerouter, boolT inorder) {
00719 facetT *facet;
00720 vertexT *vertex;
00721 int numcenters= 1;
00722 int totridges= 0;
00723
00724 qh_clearcenters(qh_ASvoronoi);
00725 qh_vertexneighbors();
00726 maximize_(qh visit_id, (unsigned) qh num_facets);
00727 FORALLfacets {
00728 facet->visitid= 0;
00729 facet->seen= False;
00730 facet->seen2= True;
00731 }
00732 FORALLfacets {
00733 if (facet->upperdelaunay == isUpper)
00734 facet->visitid= numcenters++;
00735 }
00736 FORALLvertices
00737 vertex->seen= False;
00738 FORALLvertices {
00739 if (qh GOODvertex > 0 && qh_pointid(vertex->point)+1 != qh GOODvertex)
00740 continue;
00741 totridges += qh_eachvoronoi(fp, printvridge, vertex,
00742 !qh_ALL, innerouter, inorder);
00743 }
00744 return totridges;
00745 }
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758 void qh_facet2point(facetT *facet, pointT **point0, pointT **point1, realT *mindist) {
00759 vertexT *vertex0, *vertex1;
00760 realT dist;
00761
00762 if (facet->toporient ^ qh_ORIENTclock) {
00763 vertex0= SETfirstt_(facet->vertices, vertexT);
00764 vertex1= SETsecondt_(facet->vertices, vertexT);
00765 }else {
00766 vertex1= SETfirstt_(facet->vertices, vertexT);
00767 vertex0= SETsecondt_(facet->vertices, vertexT);
00768 }
00769 zadd_(Zdistio, 2);
00770 qh_distplane(vertex0->point, facet, &dist);
00771 *mindist= dist;
00772 *point0= qh_projectpoint(vertex0->point, facet, dist);
00773 qh_distplane(vertex1->point, facet, &dist);
00774 minimize_(*mindist, dist);
00775 *point1= qh_projectpoint(vertex1->point, facet, dist);
00776 }
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799 setT *qh_facetvertices(facetT *facetlist, setT *facets, boolT allfacets) {
00800 setT *vertices;
00801 facetT *facet, **facetp;
00802 vertexT *vertex, **vertexp;
00803
00804 qh vertex_visit++;
00805 if (facetlist == qh facet_list && allfacets && !facets) {
00806 vertices= qh_settemp(qh num_vertices);
00807 FORALLvertices {
00808 vertex->visitid= qh vertex_visit;
00809 qh_setappend(&vertices, vertex);
00810 }
00811 }else {
00812 vertices= qh_settemp(qh TEMPsize);
00813 FORALLfacet_(facetlist) {
00814 if (!allfacets && qh_skipfacet(facet))
00815 continue;
00816 FOREACHvertex_(facet->vertices) {
00817 if (vertex->visitid != qh vertex_visit) {
00818 vertex->visitid= qh vertex_visit;
00819 qh_setappend(&vertices, vertex);
00820 }
00821 }
00822 }
00823 }
00824 FOREACHfacet_(facets) {
00825 if (!allfacets && qh_skipfacet(facet))
00826 continue;
00827 FOREACHvertex_(facet->vertices) {
00828 if (vertex->visitid != qh vertex_visit) {
00829 vertex->visitid= qh vertex_visit;
00830 qh_setappend(&vertices, vertex);
00831 }
00832 }
00833 }
00834 return vertices;
00835 }
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847 void qh_geomplanes(facetT *facet, realT *outerplane, realT *innerplane) {
00848 realT radius;
00849
00850 if (qh MERGING || qh JOGGLEmax < REALmax/2) {
00851 qh_outerinner(facet, outerplane, innerplane);
00852 radius= qh PRINTradius;
00853 if (qh JOGGLEmax < REALmax/2)
00854 radius -= qh JOGGLEmax * sqrt((realT)qh hull_dim);
00855 *outerplane += radius;
00856 *innerplane -= radius;
00857 if (qh PRINTcoplanar || qh PRINTspheres) {
00858 *outerplane += qh MAXabs_coord * qh_GEOMepsilon;
00859 *innerplane -= qh MAXabs_coord * qh_GEOMepsilon;
00860 }
00861 }else
00862 *innerplane= *outerplane= 0;
00863 }
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889 void qh_markkeep(facetT *facetlist) {
00890 facetT *facet, **facetp;
00891 setT *facets= qh_settemp(qh num_facets);
00892 int size, count;
00893
00894 trace2((qh ferr, 2006, "qh_markkeep: only keep %d largest and/or %d most merged facets and/or min area %.2g\n",
00895 qh KEEParea, qh KEEPmerge, qh KEEPminArea));
00896 FORALLfacet_(facetlist) {
00897 if (!facet->visible && facet->good)
00898 qh_setappend(&facets, facet);
00899 }
00900 size= qh_setsize(facets);
00901 if (qh KEEParea) {
00902 qsort(SETaddr_(facets, facetT), (size_t)size,
00903 sizeof(facetT *), qh_compare_facetarea);
00904 if ((count= size - qh KEEParea) > 0) {
00905 FOREACHfacet_(facets) {
00906 facet->good= False;
00907 if (--count == 0)
00908 break;
00909 }
00910 }
00911 }
00912 if (qh KEEPmerge) {
00913 qsort(SETaddr_(facets, facetT), (size_t)size,
00914 sizeof(facetT *), qh_compare_facetmerge);
00915 if ((count= size - qh KEEPmerge) > 0) {
00916 FOREACHfacet_(facets) {
00917 facet->good= False;
00918 if (--count == 0)
00919 break;
00920 }
00921 }
00922 }
00923 if (qh KEEPminArea < REALmax/2) {
00924 FOREACHfacet_(facets) {
00925 if (!facet->isarea || facet->f.area < qh KEEPminArea)
00926 facet->good= False;
00927 }
00928 }
00929 qh_settempfree(&facets);
00930 count= 0;
00931 FORALLfacet_(facetlist) {
00932 if (facet->good)
00933 count++;
00934 }
00935 qh num_good= count;
00936 }
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962 setT *qh_markvoronoi(facetT *facetlist, setT *facets, boolT printall, boolT *isLowerp, int *numcentersp) {
00963 int numcenters=0;
00964 facetT *facet, **facetp;
00965 setT *vertices;
00966 boolT isLower= False;
00967
00968 qh printoutnum++;
00969 qh_clearcenters(qh_ASvoronoi);
00970 qh_vertexneighbors();
00971 vertices= qh_pointvertex();
00972 if (qh ATinfinity)
00973 SETelem_(vertices, qh num_points-1)= NULL;
00974 qh visit_id++;
00975 maximize_(qh visit_id, (unsigned) qh num_facets);
00976 FORALLfacet_(facetlist) {
00977 if (printall || !qh_skipfacet(facet)) {
00978 if (!facet->upperdelaunay) {
00979 isLower= True;
00980 break;
00981 }
00982 }
00983 }
00984 FOREACHfacet_(facets) {
00985 if (printall || !qh_skipfacet(facet)) {
00986 if (!facet->upperdelaunay) {
00987 isLower= True;
00988 break;
00989 }
00990 }
00991 }
00992 FORALLfacets {
00993 if (facet->normal && (facet->upperdelaunay == isLower))
00994 facet->visitid= 0;
00995 else
00996 facet->visitid= qh visit_id;
00997 facet->seen= False;
00998 facet->seen2= True;
00999 }
01000 numcenters++;
01001 FORALLfacet_(facetlist) {
01002 if (printall || !qh_skipfacet(facet))
01003 facet->visitid= numcenters++;
01004 }
01005 FOREACHfacet_(facets) {
01006 if (printall || !qh_skipfacet(facet))
01007 facet->visitid= numcenters++;
01008 }
01009 *isLowerp= isLower;
01010 *numcentersp= numcenters;
01011 trace2((qh ferr, 2007, "qh_markvoronoi: isLower %d numcenters %d\n", isLower, numcenters));
01012 return vertices;
01013 }
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030 void qh_order_vertexneighbors(vertexT *vertex) {
01031 setT *newset;
01032 facetT *facet, *neighbor, **neighborp;
01033
01034 trace4((qh ferr, 4018, "qh_order_vertexneighbors: order neighbors of v%d for 3-d\n", vertex->id));
01035 newset= qh_settemp(qh_setsize(vertex->neighbors));
01036 facet= (facetT*)qh_setdellast(vertex->neighbors);
01037 qh_setappend(&newset, facet);
01038 while (qh_setsize(vertex->neighbors)) {
01039 FOREACHneighbor_(vertex) {
01040 if (qh_setin(facet->neighbors, neighbor)) {
01041 qh_setdel(vertex->neighbors, neighbor);
01042 qh_setappend(&newset, neighbor);
01043 facet= neighbor;
01044 break;
01045 }
01046 }
01047 if (!neighbor) {
01048 qh_fprintf(qh ferr, 6066, "qhull internal error (qh_order_vertexneighbors): no neighbor of v%d for f%d\n",
01049 vertex->id, facet->id);
01050 qh_errexit(qh_ERRqhull, facet, NULL);
01051 }
01052 }
01053 qh_setfree(&vertex->neighbors);
01054 qh_settemppop();
01055 vertex->neighbors= newset;
01056 }
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069 void qh_prepare_output(void) {
01070 if (qh VORONOI) {
01071 qh_clearcenters (qh_ASvoronoi);
01072 qh_vertexneighbors();
01073 }
01074 if (qh TRIangulate && !qh hasTriangulation) {
01075 qh_triangulate();
01076 if (qh VERIFYoutput && !qh CHECKfrequently)
01077 qh_checkpolygon (qh facet_list);
01078 }
01079 qh_findgood_all (qh facet_list);
01080 if (qh GETarea)
01081 qh_getarea(qh facet_list);
01082 if (qh KEEParea || qh KEEPmerge || qh KEEPminArea < REALmax/2)
01083 qh_markkeep (qh facet_list);
01084 if (qh PRINTstatistics)
01085 qh_collectstatistics();
01086 }
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111 void qh_printafacet(FILE *fp, qh_PRINT format, facetT *facet, boolT printall) {
01112 realT color[4], offset, dist, outerplane, innerplane;
01113 boolT zerodiv;
01114 coordT *point, *normp, *coordp, **pointp, *feasiblep;
01115 int k;
01116 vertexT *vertex, **vertexp;
01117 facetT *neighbor, **neighborp;
01118
01119 if (!printall && qh_skipfacet(facet))
01120 return;
01121 if (facet->visible && qh NEWfacets && format != qh_PRINTfacets)
01122 return;
01123 qh printoutnum++;
01124 switch (format) {
01125 case qh_PRINTarea:
01126 if (facet->isarea) {
01127 qh_fprintf(fp, 9009, qh_REAL_1, facet->f.area);
01128 qh_fprintf(fp, 9010, "\n");
01129 }else
01130 qh_fprintf(fp, 9011, "0\n");
01131 break;
01132 case qh_PRINTcoplanars:
01133 qh_fprintf(fp, 9012, "%d", qh_setsize(facet->coplanarset));
01134 FOREACHpoint_(facet->coplanarset)
01135 qh_fprintf(fp, 9013, " %d", qh_pointid(point));
01136 qh_fprintf(fp, 9014, "\n");
01137 break;
01138 case qh_PRINTcentrums:
01139 qh_printcenter(fp, format, NULL, facet);
01140 break;
01141 case qh_PRINTfacets:
01142 qh_printfacet(fp, facet);
01143 break;
01144 case qh_PRINTfacets_xridge:
01145 qh_printfacetheader(fp, facet);
01146 break;
01147 case qh_PRINTgeom:
01148 if (!facet->normal)
01149 break;
01150 for (k=qh hull_dim; k--; ) {
01151 color[k]= (facet->normal[k]+1.0)/2.0;
01152 maximize_(color[k], -1.0);
01153 minimize_(color[k], +1.0);
01154 }
01155 qh_projectdim3 (color, color);
01156 if (qh PRINTdim != qh hull_dim)
01157 qh_normalize2 (color, 3, True, NULL, NULL);
01158 if (qh hull_dim <= 2)
01159 qh_printfacet2geom(fp, facet, color);
01160 else if (qh hull_dim == 3) {
01161 if (facet->simplicial)
01162 qh_printfacet3geom_simplicial(fp, facet, color);
01163 else
01164 qh_printfacet3geom_nonsimplicial(fp, facet, color);
01165 }else {
01166 if (facet->simplicial)
01167 qh_printfacet4geom_simplicial(fp, facet, color);
01168 else
01169 qh_printfacet4geom_nonsimplicial(fp, facet, color);
01170 }
01171 break;
01172 case qh_PRINTids:
01173 qh_fprintf(fp, 9015, "%d\n", facet->id);
01174 break;
01175 case qh_PRINTincidences:
01176 case qh_PRINToff:
01177 case qh_PRINTtriangles:
01178 if (qh hull_dim == 3 && format != qh_PRINTtriangles)
01179 qh_printfacet3vertex(fp, facet, format);
01180 else if (facet->simplicial || qh hull_dim == 2 || format == qh_PRINToff)
01181 qh_printfacetNvertex_simplicial(fp, facet, format);
01182 else
01183 qh_printfacetNvertex_nonsimplicial(fp, facet, qh printoutvar++, format);
01184 break;
01185 case qh_PRINTinner:
01186 qh_outerinner(facet, NULL, &innerplane);
01187 offset= facet->offset - innerplane;
01188 goto LABELprintnorm;
01189 break;
01190 case qh_PRINTmerges:
01191 qh_fprintf(fp, 9016, "%d\n", facet->nummerge);
01192 break;
01193 case qh_PRINTnormals:
01194 offset= facet->offset;
01195 goto LABELprintnorm;
01196 break;
01197 case qh_PRINTouter:
01198 qh_outerinner(facet, &outerplane, NULL);
01199 offset= facet->offset - outerplane;
01200 LABELprintnorm:
01201 if (!facet->normal) {
01202 qh_fprintf(fp, 9017, "no normal for facet f%d\n", facet->id);
01203 break;
01204 }
01205 if (qh CDDoutput) {
01206 qh_fprintf(fp, 9018, qh_REAL_1, -offset);
01207 for (k=0; k < qh hull_dim; k++)
01208 qh_fprintf(fp, 9019, qh_REAL_1, -facet->normal[k]);
01209 }else {
01210 for (k=0; k < qh hull_dim; k++)
01211 qh_fprintf(fp, 9020, qh_REAL_1, facet->normal[k]);
01212 qh_fprintf(fp, 9021, qh_REAL_1, offset);
01213 }
01214 qh_fprintf(fp, 9022, "\n");
01215 break;
01216 case qh_PRINTmathematica:
01217 case qh_PRINTmaple:
01218 if (qh hull_dim == 2)
01219 qh_printfacet2math(fp, facet, format, qh printoutvar++);
01220 else
01221 qh_printfacet3math(fp, facet, format, qh printoutvar++);
01222 break;
01223 case qh_PRINTneighbors:
01224 qh_fprintf(fp, 9023, "%d", qh_setsize(facet->neighbors));
01225 FOREACHneighbor_(facet)
01226 qh_fprintf(fp, 9024, " %d",
01227 neighbor->visitid ? neighbor->visitid - 1: 0 - neighbor->id);
01228 qh_fprintf(fp, 9025, "\n");
01229 break;
01230 case qh_PRINTpointintersect:
01231 if (!qh feasible_point) {
01232 qh_fprintf(qh ferr, 6067, "qhull input error (qh_printafacet): option 'Fp' needs qh feasible_point\n");
01233 qh_errexit( qh_ERRinput, NULL, NULL);
01234 }
01235 if (facet->offset > 0)
01236 goto LABELprintinfinite;
01237 point= coordp= (coordT*)qh_memalloc(qh normal_size);
01238 normp= facet->normal;
01239 feasiblep= qh feasible_point;
01240 if (facet->offset < -qh MINdenom) {
01241 for (k=qh hull_dim; k--; )
01242 *(coordp++)= (*(normp++) / - facet->offset) + *(feasiblep++);
01243 }else {
01244 for (k=qh hull_dim; k--; ) {
01245 *(coordp++)= qh_divzero(*(normp++), facet->offset, qh MINdenom_1,
01246 &zerodiv) + *(feasiblep++);
01247 if (zerodiv) {
01248 qh_memfree(point, qh normal_size);
01249 goto LABELprintinfinite;
01250 }
01251 }
01252 }
01253 qh_printpoint(fp, NULL, point);
01254 qh_memfree(point, qh normal_size);
01255 break;
01256 LABELprintinfinite:
01257 for (k=qh hull_dim; k--; )
01258 qh_fprintf(fp, 9026, qh_REAL_1, qh_INFINITE);
01259 qh_fprintf(fp, 9027, "\n");
01260 break;
01261 case qh_PRINTpointnearest:
01262 FOREACHpoint_(facet->coplanarset) {
01263 int id, id2;
01264 vertex= qh_nearvertex(facet, point, &dist);
01265 id= qh_pointid(vertex->point);
01266 id2= qh_pointid(point);
01267 qh_fprintf(fp, 9028, "%d %d %d " qh_REAL_1 "\n", id, id2, facet->id, dist);
01268 }
01269 break;
01270 case qh_PRINTpoints:
01271 if (qh CDDoutput)
01272 qh_fprintf(fp, 9029, "1 ");
01273 qh_printcenter(fp, format, NULL, facet);
01274 break;
01275 case qh_PRINTvertices:
01276 qh_fprintf(fp, 9030, "%d", qh_setsize(facet->vertices));
01277 FOREACHvertex_(facet->vertices)
01278 qh_fprintf(fp, 9031, " %d", qh_pointid(vertex->point));
01279 qh_fprintf(fp, 9032, "\n");
01280 break;
01281 default:
01282 break;
01283 }
01284 }
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307 void qh_printbegin(FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall) {
01308 int numfacets, numsimplicial, numridges, totneighbors, numcoplanars, numtricoplanars;
01309 int i, num;
01310 facetT *facet, **facetp;
01311 vertexT *vertex, **vertexp;
01312 setT *vertices;
01313 pointT *point, **pointp, *pointtemp;
01314
01315 qh printoutnum= 0;
01316 qh_countfacets(facetlist, facets, printall, &numfacets, &numsimplicial,
01317 &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
01318 switch (format) {
01319 case qh_PRINTnone:
01320 break;
01321 case qh_PRINTarea:
01322 qh_fprintf(fp, 9033, "%d\n", numfacets);
01323 break;
01324 case qh_PRINTcoplanars:
01325 qh_fprintf(fp, 9034, "%d\n", numfacets);
01326 break;
01327 case qh_PRINTcentrums:
01328 if (qh CENTERtype == qh_ASnone)
01329 qh_clearcenters(qh_AScentrum);
01330 qh_fprintf(fp, 9035, "%d\n%d\n", qh hull_dim, numfacets);
01331 break;
01332 case qh_PRINTfacets:
01333 case qh_PRINTfacets_xridge:
01334 if (facetlist)
01335 qh_printvertexlist(fp, "Vertices and facets:\n", facetlist, facets, printall);
01336 break;
01337 case qh_PRINTgeom:
01338 if (qh hull_dim > 4)
01339 goto LABELnoformat;
01340 if (qh VORONOI && qh hull_dim > 3)
01341 goto LABELnoformat;
01342 if (qh hull_dim == 2 && (qh PRINTridges || qh DOintersections))
01343 qh_fprintf(qh ferr, 7049, "qhull warning: output for ridges and intersections not implemented in 2-d\n");
01344 if (qh hull_dim == 4 && (qh PRINTinner || qh PRINTouter ||
01345 (qh PRINTdim == 4 && qh PRINTcentrums)))
01346 qh_fprintf(qh ferr, 7050, "qhull warning: output for outer/inner planes and centrums not implemented in 4-d\n");
01347 if (qh PRINTdim == 4 && (qh PRINTspheres))
01348 qh_fprintf(qh ferr, 7051, "qhull warning: output for vertices not implemented in 4-d\n");
01349 if (qh PRINTdim == 4 && qh DOintersections && qh PRINTnoplanes)
01350 qh_fprintf(qh ferr, 7052, "qhull warning: 'Gnh' generates no output in 4-d\n");
01351 if (qh PRINTdim == 2) {
01352 qh_fprintf(fp, 9036, "{appearance {linewidth 3} LIST # %s | %s\n",
01353 qh rbox_command, qh qhull_command);
01354 }else if (qh PRINTdim == 3) {
01355 qh_fprintf(fp, 9037, "{appearance {+edge -evert linewidth 2} LIST # %s | %s\n",
01356 qh rbox_command, qh qhull_command);
01357 }else if (qh PRINTdim == 4) {
01358 qh visit_id++;
01359 num= 0;
01360 FORALLfacet_(facetlist)
01361 qh_printend4geom(NULL, facet, &num, printall);
01362 FOREACHfacet_(facets)
01363 qh_printend4geom(NULL, facet, &num, printall);
01364 qh ridgeoutnum= num;
01365 qh printoutvar= 0;
01366 qh_fprintf(fp, 9038, "LIST # %s | %s\n", qh rbox_command, qh qhull_command);
01367 }
01368
01369 if (qh PRINTdots) {
01370 qh printoutnum++;
01371 num= qh num_points + qh_setsize(qh other_points);
01372 if (qh DELAUNAY && qh ATinfinity)
01373 num--;
01374 if (qh PRINTdim == 4)
01375 qh_fprintf(fp, 9039, "4VECT %d %d 1\n", num, num);
01376 else
01377 qh_fprintf(fp, 9040, "VECT %d %d 1\n", num, num);
01378
01379 for (i=num; i--; ) {
01380 if (i % 20 == 0)
01381 qh_fprintf(fp, 9041, "\n");
01382 qh_fprintf(fp, 9042, "1 ");
01383 }
01384 qh_fprintf(fp, 9043, "# 1 point per line\n1 ");
01385 for (i=num-1; i--; ) {
01386 if (i % 20 == 0)
01387 qh_fprintf(fp, 9044, "\n");
01388 qh_fprintf(fp, 9045, "0 ");
01389 }
01390 qh_fprintf(fp, 9046, "# 1 color for all\n");
01391 FORALLpoints {
01392 if (!qh DELAUNAY || !qh ATinfinity || qh_pointid(point) != qh num_points-1) {
01393 if (qh PRINTdim == 4)
01394 qh_printpoint(fp, NULL, point);
01395 else
01396 qh_printpoint3 (fp, point);
01397 }
01398 }
01399 FOREACHpoint_(qh other_points) {
01400 if (qh PRINTdim == 4)
01401 qh_printpoint(fp, NULL, point);
01402 else
01403 qh_printpoint3 (fp, point);
01404 }
01405 qh_fprintf(fp, 9047, "0 1 1 1 # color of points\n");
01406 }
01407
01408 if (qh PRINTdim == 4 && !qh PRINTnoplanes)
01409
01410 qh_fprintf(fp, 9048, "4OFF %d %d 1\n", 3*qh ridgeoutnum, qh ridgeoutnum);
01411 qh PRINTcradius= 2 * qh DISTround;
01412 if (qh PREmerge) {
01413 maximize_(qh PRINTcradius, qh premerge_centrum + qh DISTround);
01414 }else if (qh POSTmerge)
01415 maximize_(qh PRINTcradius, qh postmerge_centrum + qh DISTround);
01416 qh PRINTradius= qh PRINTcradius;
01417 if (qh PRINTspheres + qh PRINTcoplanar)
01418 maximize_(qh PRINTradius, qh MAXabs_coord * qh_MINradius);
01419 if (qh premerge_cos < REALmax/2) {
01420 maximize_(qh PRINTradius, (1- qh premerge_cos) * qh MAXabs_coord);
01421 }else if (!qh PREmerge && qh POSTmerge && qh postmerge_cos < REALmax/2) {
01422 maximize_(qh PRINTradius, (1- qh postmerge_cos) * qh MAXabs_coord);
01423 }
01424 maximize_(qh PRINTradius, qh MINvisible);
01425 if (qh JOGGLEmax < REALmax/2)
01426 qh PRINTradius += qh JOGGLEmax * sqrt((realT)qh hull_dim);
01427 if (qh PRINTdim != 4 &&
01428 (qh PRINTcoplanar || qh PRINTspheres || qh PRINTcentrums)) {
01429 vertices= qh_facetvertices(facetlist, facets, printall);
01430 if (qh PRINTspheres && qh PRINTdim <= 3)
01431 qh_printspheres(fp, vertices, qh PRINTradius);
01432 if (qh PRINTcoplanar || qh PRINTcentrums) {
01433 qh firstcentrum= True;
01434 if (qh PRINTcoplanar&& !qh PRINTspheres) {
01435 FOREACHvertex_(vertices)
01436 qh_printpointvect2 (fp, vertex->point, NULL, qh interior_point, qh PRINTradius);
01437 }
01438 FORALLfacet_(facetlist) {
01439 if (!printall && qh_skipfacet(facet))
01440 continue;
01441 if (!facet->normal)
01442 continue;
01443 if (qh PRINTcentrums && qh PRINTdim <= 3)
01444 qh_printcentrum(fp, facet, qh PRINTcradius);
01445 if (!qh PRINTcoplanar)
01446 continue;
01447 FOREACHpoint_(facet->coplanarset)
01448 qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius);
01449 FOREACHpoint_(facet->outsideset)
01450 qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius);
01451 }
01452 FOREACHfacet_(facets) {
01453 if (!printall && qh_skipfacet(facet))
01454 continue;
01455 if (!facet->normal)
01456 continue;
01457 if (qh PRINTcentrums && qh PRINTdim <= 3)
01458 qh_printcentrum(fp, facet, qh PRINTcradius);
01459 if (!qh PRINTcoplanar)
01460 continue;
01461 FOREACHpoint_(facet->coplanarset)
01462 qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius);
01463 FOREACHpoint_(facet->outsideset)
01464 qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius);
01465 }
01466 }
01467 qh_settempfree(&vertices);
01468 }
01469 qh visit_id++;
01470 break;
01471 case qh_PRINTids:
01472 qh_fprintf(fp, 9049, "%d\n", numfacets);
01473 break;
01474 case qh_PRINTincidences:
01475 if (qh VORONOI && qh PRINTprecision)
01476 qh_fprintf(qh ferr, 7053, "qhull warning: writing Delaunay. Use 'p' or 'o' for Voronoi centers\n");
01477 qh printoutvar= qh vertex_id;
01478 if (qh hull_dim <= 3)
01479 qh_fprintf(fp, 9050, "%d\n", numfacets);
01480 else
01481 qh_fprintf(fp, 9051, "%d\n", numsimplicial+numridges);
01482 break;
01483 case qh_PRINTinner:
01484 case qh_PRINTnormals:
01485 case qh_PRINTouter:
01486 if (qh CDDoutput)
01487 qh_fprintf(fp, 9052, "%s | %s\nbegin\n %d %d real\n", qh rbox_command,
01488 qh qhull_command, numfacets, qh hull_dim+1);
01489 else
01490 qh_fprintf(fp, 9053, "%d\n%d\n", qh hull_dim+1, numfacets);
01491 break;
01492 case qh_PRINTmathematica:
01493 case qh_PRINTmaple:
01494 if (qh hull_dim > 3)
01495 goto LABELnoformat;
01496 if (qh VORONOI)
01497 qh_fprintf(qh ferr, 7054, "qhull warning: output is the Delaunay triangulation\n");
01498 if (format == qh_PRINTmaple) {
01499 if (qh hull_dim == 2)
01500 qh_fprintf(fp, 9054, "PLOT(CURVES(\n");
01501 else
01502 qh_fprintf(fp, 9055, "PLOT3D(POLYGONS(\n");
01503 }else
01504 qh_fprintf(fp, 9056, "{\n");
01505 qh printoutvar= 0;
01506 break;
01507 case qh_PRINTmerges:
01508 qh_fprintf(fp, 9057, "%d\n", numfacets);
01509 break;
01510 case qh_PRINTpointintersect:
01511 qh_fprintf(fp, 9058, "%d\n%d\n", qh hull_dim, numfacets);
01512 break;
01513 case qh_PRINTneighbors:
01514 qh_fprintf(fp, 9059, "%d\n", numfacets);
01515 break;
01516 case qh_PRINToff:
01517 case qh_PRINTtriangles:
01518 if (qh VORONOI)
01519 goto LABELnoformat;
01520 num = qh hull_dim;
01521 if (format == qh_PRINToff || qh hull_dim == 2)
01522 qh_fprintf(fp, 9060, "%d\n%d %d %d\n", num,
01523 qh num_points+qh_setsize(qh other_points), numfacets, totneighbors/2);
01524 else {
01525 qh printoutvar= qh num_points+qh_setsize(qh other_points);
01526 if (qh DELAUNAY)
01527 num--;
01528 qh_fprintf(fp, 9061, "%d\n%d %d %d\n", num, qh printoutvar
01529 + numfacets - numsimplicial, numsimplicial + numridges, totneighbors/2);
01530 }
01531 FORALLpoints
01532 qh_printpointid(qh fout, NULL, num, point, -1);
01533 FOREACHpoint_(qh other_points)
01534 qh_printpointid(qh fout, NULL, num, point, -1);
01535 if (format == qh_PRINTtriangles && qh hull_dim > 2) {
01536 FORALLfacets {
01537 if (!facet->simplicial && facet->visitid)
01538 qh_printcenter(qh fout, format, NULL, facet);
01539 }
01540 }
01541 break;
01542 case qh_PRINTpointnearest:
01543 qh_fprintf(fp, 9062, "%d\n", numcoplanars);
01544 break;
01545 case qh_PRINTpoints:
01546 if (!qh VORONOI)
01547 goto LABELnoformat;
01548 if (qh CDDoutput)
01549 qh_fprintf(fp, 9063, "%s | %s\nbegin\n%d %d real\n", qh rbox_command,
01550 qh qhull_command, numfacets, qh hull_dim);
01551 else
01552 qh_fprintf(fp, 9064, "%d\n%d\n", qh hull_dim-1, numfacets);
01553 break;
01554 case qh_PRINTvertices:
01555 qh_fprintf(fp, 9065, "%d\n", numfacets);
01556 break;
01557 case qh_PRINTsummary:
01558 default:
01559 LABELnoformat:
01560 qh_fprintf(qh ferr, 6068, "qhull internal error (qh_printbegin): can not use this format for dimension %d\n",
01561 qh hull_dim);
01562 qh_errexit(qh_ERRqhull, NULL, NULL);
01563 }
01564 }
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575
01576
01577
01578
01579
01580
01581 void qh_printcenter(FILE *fp, qh_PRINT format, const char *string, facetT *facet) {
01582 int k, num;
01583
01584 if (qh CENTERtype != qh_ASvoronoi && qh CENTERtype != qh_AScentrum)
01585 return;
01586 if (string)
01587 qh_fprintf(fp, 9066, string);
01588 if (qh CENTERtype == qh_ASvoronoi) {
01589 num= qh hull_dim-1;
01590 if (!facet->normal || !facet->upperdelaunay || !qh ATinfinity) {
01591 if (!facet->center)
01592 facet->center= qh_facetcenter(facet->vertices);
01593 for (k=0; k < num; k++)
01594 qh_fprintf(fp, 9067, qh_REAL_1, facet->center[k]);
01595 }else {
01596 for (k=0; k < num; k++)
01597 qh_fprintf(fp, 9068, qh_REAL_1, qh_INFINITE);
01598 }
01599 }else {
01600 num= qh hull_dim;
01601 if (format == qh_PRINTtriangles && qh DELAUNAY)
01602 num--;
01603 if (!facet->center)
01604 facet->center= qh_getcentrum(facet);
01605 for (k=0; k < num; k++)
01606 qh_fprintf(fp, 9069, qh_REAL_1, facet->center[k]);
01607 }
01608 if (format == qh_PRINTgeom && num == 2)
01609 qh_fprintf(fp, 9070, " 0\n");
01610 else
01611 qh_fprintf(fp, 9071, "\n");
01612 }
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625 void qh_printcentrum(FILE *fp, facetT *facet, realT radius) {
01626 pointT *centrum, *projpt;
01627 boolT tempcentrum= False;
01628 realT xaxis[4], yaxis[4], normal[4], dist;
01629 realT green[3]={0, 1, 0};
01630 vertexT *apex;
01631 int k;
01632
01633 if (qh CENTERtype == qh_AScentrum) {
01634 if (!facet->center)
01635 facet->center= qh_getcentrum(facet);
01636 centrum= facet->center;
01637 }else {
01638 centrum= qh_getcentrum(facet);
01639 tempcentrum= True;
01640 }
01641 qh_fprintf(fp, 9072, "{appearance {-normal -edge normscale 0} ");
01642 if (qh firstcentrum) {
01643 qh firstcentrum= False;
01644 qh_fprintf(fp, 9073, "{INST geom { define centrum CQUAD # f%d\n\
01645 -0.3 -0.3 0.0001 0 0 1 1\n\
01646 0.3 -0.3 0.0001 0 0 1 1\n\
01647 0.3 0.3 0.0001 0 0 1 1\n\
01648 -0.3 0.3 0.0001 0 0 1 1 } transform { \n", facet->id);
01649 }else
01650 qh_fprintf(fp, 9074, "{INST geom { : centrum } transform { # f%d\n", facet->id);
01651 apex= SETfirstt_(facet->vertices, vertexT);
01652 qh_distplane(apex->point, facet, &dist);
01653 projpt= qh_projectpoint(apex->point, facet, dist);
01654 for (k=qh hull_dim; k--; ) {
01655 xaxis[k]= projpt[k] - centrum[k];
01656 normal[k]= facet->normal[k];
01657 }
01658 if (qh hull_dim == 2) {
01659 xaxis[2]= 0;
01660 normal[2]= 0;
01661 }else if (qh hull_dim == 4) {
01662 qh_projectdim3 (xaxis, xaxis);
01663 qh_projectdim3 (normal, normal);
01664 qh_normalize2 (normal, qh PRINTdim, True, NULL, NULL);
01665 }
01666 qh_crossproduct(3, xaxis, normal, yaxis);
01667 qh_fprintf(fp, 9075, "%8.4g %8.4g %8.4g 0\n", xaxis[0], xaxis[1], xaxis[2]);
01668 qh_fprintf(fp, 9076, "%8.4g %8.4g %8.4g 0\n", yaxis[0], yaxis[1], yaxis[2]);
01669 qh_fprintf(fp, 9077, "%8.4g %8.4g %8.4g 0\n", normal[0], normal[1], normal[2]);
01670 qh_printpoint3 (fp, centrum);
01671 qh_fprintf(fp, 9078, "1 }}}\n");
01672 qh_memfree(projpt, qh normal_size);
01673 qh_printpointvect(fp, centrum, facet->normal, NULL, radius, green);
01674 if (tempcentrum)
01675 qh_memfree(centrum, qh normal_size);
01676 }
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688 void qh_printend(FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall) {
01689 int num;
01690 facetT *facet, **facetp;
01691
01692 if (!qh printoutnum)
01693 qh_fprintf(qh ferr, 7055, "qhull warning: no facets printed\n");
01694 switch (format) {
01695 case qh_PRINTgeom:
01696 if (qh hull_dim == 4 && qh DROPdim < 0 && !qh PRINTnoplanes) {
01697 qh visit_id++;
01698 num= 0;
01699 FORALLfacet_(facetlist)
01700 qh_printend4geom(fp, facet,&num, printall);
01701 FOREACHfacet_(facets)
01702 qh_printend4geom(fp, facet, &num, printall);
01703 if (num != qh ridgeoutnum || qh printoutvar != qh ridgeoutnum) {
01704 qh_fprintf(qh ferr, 6069, "qhull internal error (qh_printend): number of ridges %d != number printed %d and at end %d\n", qh ridgeoutnum, qh printoutvar, num);
01705 qh_errexit(qh_ERRqhull, NULL, NULL);
01706 }
01707 }else
01708 qh_fprintf(fp, 9079, "}\n");
01709 break;
01710 case qh_PRINTinner:
01711 case qh_PRINTnormals:
01712 case qh_PRINTouter:
01713 if (qh CDDoutput)
01714 qh_fprintf(fp, 9080, "end\n");
01715 break;
01716 case qh_PRINTmaple:
01717 qh_fprintf(fp, 9081, "));\n");
01718 break;
01719 case qh_PRINTmathematica:
01720 qh_fprintf(fp, 9082, "}\n");
01721 break;
01722 case qh_PRINTpoints:
01723 if (qh CDDoutput)
01724 qh_fprintf(fp, 9083, "end\n");
01725 break;
01726 default:
01727 break;
01728 }
01729 }
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749 void qh_printend4geom(FILE *fp, facetT *facet, int *nump, boolT printall) {
01750 realT color[3];
01751 int i, num= *nump;
01752 facetT *neighbor, **neighborp;
01753 ridgeT *ridge, **ridgep;
01754
01755 if (!printall && qh_skipfacet(facet))
01756 return;
01757 if (qh PRINTnoplanes || (facet->visible && qh NEWfacets))
01758 return;
01759 if (!facet->normal)
01760 return;
01761 if (fp) {
01762 for (i=0; i < 3; i++) {
01763 color[i]= (facet->normal[i]+1.0)/2.0;
01764 maximize_(color[i], -1.0);
01765 minimize_(color[i], +1.0);
01766 }
01767 }
01768 facet->visitid= qh visit_id;
01769 if (facet->simplicial) {
01770 FOREACHneighbor_(facet) {
01771 if (neighbor->visitid != qh visit_id) {
01772 if (fp)
01773 qh_fprintf(fp, 9084, "3 %d %d %d %8.4g %8.4g %8.4g 1 # f%d f%d\n",
01774 3*num, 3*num+1, 3*num+2, color[0], color[1], color[2],
01775 facet->id, neighbor->id);
01776 num++;
01777 }
01778 }
01779 }else {
01780 FOREACHridge_(facet->ridges) {
01781 neighbor= otherfacet_(ridge, facet);
01782 if (neighbor->visitid != qh visit_id) {
01783 if (fp)
01784 qh_fprintf(fp, 9085, "3 %d %d %d %8.4g %8.4g %8.4g 1 #r%d f%d f%d\n",
01785 3*num, 3*num+1, 3*num+2, color[0], color[1], color[2],
01786 ridge->id, facet->id, neighbor->id);
01787 num++;
01788 }
01789 }
01790 }
01791 *nump= num;
01792 }
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806 void qh_printextremes(FILE *fp, facetT *facetlist, setT *facets, boolT printall) {
01807 setT *vertices, *points;
01808 pointT *point;
01809 vertexT *vertex, **vertexp;
01810 int id;
01811 int numpoints=0, point_i, point_n;
01812 int allpoints= qh num_points + qh_setsize(qh other_points);
01813
01814 points= qh_settemp(allpoints);
01815 qh_setzero(points, 0, allpoints);
01816 vertices= qh_facetvertices(facetlist, facets, printall);
01817 FOREACHvertex_(vertices) {
01818 id= qh_pointid(vertex->point);
01819 if (id >= 0) {
01820 SETelem_(points, id)= vertex->point;
01821 numpoints++;
01822 }
01823 }
01824 qh_settempfree(&vertices);
01825 qh_fprintf(fp, 9086, "%d\n", numpoints);
01826 FOREACHpoint_i_(points) {
01827 if (point)
01828 qh_fprintf(fp, 9087, "%d\n", point_i);
01829 }
01830 qh_settempfree(&points);
01831 }
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845 void qh_printextremes_2d(FILE *fp, facetT *facetlist, setT *facets, boolT printall) {
01846 int numfacets, numridges, totneighbors, numcoplanars, numsimplicial, numtricoplanars;
01847 setT *vertices;
01848 facetT *facet, *startfacet, *nextfacet;
01849 vertexT *vertexA, *vertexB;
01850
01851 qh_countfacets(facetlist, facets, printall, &numfacets, &numsimplicial,
01852 &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
01853 vertices= qh_facetvertices(facetlist, facets, printall);
01854 qh_fprintf(fp, 9088, "%d\n", qh_setsize(vertices));
01855 qh_settempfree(&vertices);
01856 if (!numfacets)
01857 return;
01858 facet= startfacet= facetlist ? facetlist : SETfirstt_(facets, facetT);
01859 qh vertex_visit++;
01860 qh visit_id++;
01861 do {
01862 if (facet->toporient ^ qh_ORIENTclock) {
01863 vertexA= SETfirstt_(facet->vertices, vertexT);
01864 vertexB= SETsecondt_(facet->vertices, vertexT);
01865 nextfacet= SETfirstt_(facet->neighbors, facetT);
01866 }else {
01867 vertexA= SETsecondt_(facet->vertices, vertexT);
01868 vertexB= SETfirstt_(facet->vertices, vertexT);
01869 nextfacet= SETsecondt_(facet->neighbors, facetT);
01870 }
01871 if (facet->visitid == qh visit_id) {
01872 qh_fprintf(qh ferr, 6218, "Qhull internal error (qh_printextremes_2d): loop in facet list. facet %d nextfacet %d\n",
01873 facet->id, nextfacet->id);
01874 qh_errexit2 (qh_ERRqhull, facet, nextfacet);
01875 }
01876 if (facet->visitid) {
01877 if (vertexA->visitid != qh vertex_visit) {
01878 vertexA->visitid= qh vertex_visit;
01879 qh_fprintf(fp, 9089, "%d\n", qh_pointid(vertexA->point));
01880 }
01881 if (vertexB->visitid != qh vertex_visit) {
01882 vertexB->visitid= qh vertex_visit;
01883 qh_fprintf(fp, 9090, "%d\n", qh_pointid(vertexB->point));
01884 }
01885 }
01886 facet->visitid= qh visit_id;
01887 facet= nextfacet;
01888 }while (facet && facet != startfacet);
01889 }
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902 void qh_printextremes_d(FILE *fp, facetT *facetlist, setT *facets, boolT printall) {
01903 setT *vertices;
01904 vertexT *vertex, **vertexp;
01905 boolT upperseen, lowerseen;
01906 facetT *neighbor, **neighborp;
01907 int numpoints=0;
01908
01909 vertices= qh_facetvertices(facetlist, facets, printall);
01910 qh_vertexneighbors();
01911 FOREACHvertex_(vertices) {
01912 upperseen= lowerseen= False;
01913 FOREACHneighbor_(vertex) {
01914 if (neighbor->upperdelaunay)
01915 upperseen= True;
01916 else
01917 lowerseen= True;
01918 }
01919 if (upperseen && lowerseen) {
01920 vertex->seen= True;
01921 numpoints++;
01922 }else
01923 vertex->seen= False;
01924 }
01925 qh_fprintf(fp, 9091, "%d\n", numpoints);
01926 FOREACHvertex_(vertices) {
01927 if (vertex->seen)
01928 qh_fprintf(fp, 9092, "%d\n", qh_pointid(vertex->point));
01929 }
01930 qh_settempfree(&vertices);
01931 }
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942 void qh_printfacet(FILE *fp, facetT *facet) {
01943
01944 qh_printfacetheader(fp, facet);
01945 if (facet->ridges)
01946 qh_printfacetridges(fp, facet);
01947 }
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961 void qh_printfacet2geom(FILE *fp, facetT *facet, realT color[3]) {
01962 pointT *point0, *point1;
01963 realT mindist, innerplane, outerplane;
01964 int k;
01965
01966 qh_facet2point(facet, &point0, &point1, &mindist);
01967 qh_geomplanes(facet, &outerplane, &innerplane);
01968 if (qh PRINTouter || (!qh PRINTnoplanes && !qh PRINTinner))
01969 qh_printfacet2geom_points(fp, point0, point1, facet, outerplane, color);
01970 if (qh PRINTinner || (!qh PRINTnoplanes && !qh PRINTouter &&
01971 outerplane - innerplane > 2 * qh MAXabs_coord * qh_GEOMepsilon)) {
01972 for (k=3; k--; )
01973 color[k]= 1.0 - color[k];
01974 qh_printfacet2geom_points(fp, point0, point1, facet, innerplane, color);
01975 }
01976 qh_memfree(point1, qh normal_size);
01977 qh_memfree(point0, qh normal_size);
01978 }
01979
01980
01981
01982
01983
01984
01985
01986
01987 void qh_printfacet2geom_points(FILE *fp, pointT *point1, pointT *point2,
01988 facetT *facet, realT offset, realT color[3]) {
01989 pointT *p1= point1, *p2= point2;
01990
01991 qh_fprintf(fp, 9093, "VECT 1 2 1 2 1 # f%d\n", facet->id);
01992 if (offset != 0.0) {
01993 p1= qh_projectpoint(p1, facet, -offset);
01994 p2= qh_projectpoint(p2, facet, -offset);
01995 }
01996 qh_fprintf(fp, 9094, "%8.4g %8.4g %8.4g\n%8.4g %8.4g %8.4g\n",
01997 p1[0], p1[1], 0.0, p2[0], p2[1], 0.0);
01998 if (offset != 0.0) {
01999 qh_memfree(p1, qh normal_size);
02000 qh_memfree(p2, qh normal_size);
02001 }
02002 qh_fprintf(fp, 9095, "%8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]);
02003 }
02004
02005
02006
02007
02008
02009
02010
02011
02012
02013
02014
02015
02016
02017 void qh_printfacet2math(FILE *fp, facetT *facet, qh_PRINT format, int notfirst) {
02018 pointT *point0, *point1;
02019 realT mindist;
02020 const char *pointfmt;
02021
02022 qh_facet2point(facet, &point0, &point1, &mindist);
02023 if (notfirst)
02024 qh_fprintf(fp, 9096, ",");
02025 if (format == qh_PRINTmaple)
02026 pointfmt= "[[%16.8f, %16.8f], [%16.8f, %16.8f]]\n";
02027 else
02028 pointfmt= "Line[{{%16.8f, %16.8f}, {%16.8f, %16.8f}}]\n";
02029 qh_fprintf(fp, 9097, pointfmt, point0[0], point0[1], point1[0], point1[1]);
02030 qh_memfree(point1, qh normal_size);
02031 qh_memfree(point0, qh normal_size);
02032 }
02033
02034
02035
02036
02037
02038
02039
02040
02041
02042
02043
02044
02045 void qh_printfacet3geom_nonsimplicial(FILE *fp, facetT *facet, realT color[3]) {
02046 ridgeT *ridge, **ridgep;
02047 setT *projectedpoints, *vertices;
02048 vertexT *vertex, **vertexp, *vertexA, *vertexB;
02049 pointT *projpt, *point, **pointp;
02050 facetT *neighbor;
02051 realT dist, outerplane, innerplane;
02052 int cntvertices, k;
02053 realT black[3]={0, 0, 0}, green[3]={0, 1, 0};
02054
02055 qh_geomplanes(facet, &outerplane, &innerplane);
02056 vertices= qh_facet3vertex(facet);
02057 cntvertices= qh_setsize(vertices);
02058 projectedpoints= qh_settemp(cntvertices);
02059 FOREACHvertex_(vertices) {
02060 zinc_(Zdistio);
02061 qh_distplane(vertex->point, facet, &dist);
02062 projpt= qh_projectpoint(vertex->point, facet, dist);
02063 qh_setappend(&projectedpoints, projpt);
02064 }
02065 if (qh PRINTouter || (!qh PRINTnoplanes && !qh PRINTinner))
02066 qh_printfacet3geom_points(fp, projectedpoints, facet, outerplane, color);
02067 if (qh PRINTinner || (!qh PRINTnoplanes && !qh PRINTouter &&
02068 outerplane - innerplane > 2 * qh MAXabs_coord * qh_GEOMepsilon)) {
02069 for (k=3; k--; )
02070 color[k]= 1.0 - color[k];
02071 qh_printfacet3geom_points(fp, projectedpoints, facet, innerplane, color);
02072 }
02073 FOREACHpoint_(projectedpoints)
02074 qh_memfree(point, qh normal_size);
02075 qh_settempfree(&projectedpoints);
02076 qh_settempfree(&vertices);
02077 if ((qh DOintersections || qh PRINTridges)
02078 && (!facet->visible || !qh NEWfacets)) {
02079 facet->visitid= qh visit_id;
02080 FOREACHridge_(facet->ridges) {
02081 neighbor= otherfacet_(ridge, facet);
02082 if (neighbor->visitid != qh visit_id) {
02083 if (qh DOintersections)
02084 qh_printhyperplaneintersection(fp, facet, neighbor, ridge->vertices, black);
02085 if (qh PRINTridges) {
02086 vertexA= SETfirstt_(ridge->vertices, vertexT);
02087 vertexB= SETsecondt_(ridge->vertices, vertexT);
02088 qh_printline3geom(fp, vertexA->point, vertexB->point, green);
02089 }
02090 }
02091 }
02092 }
02093 }
02094
02095
02096
02097
02098
02099
02100
02101
02102
02103 void qh_printfacet3geom_points(FILE *fp, setT *points, facetT *facet, realT offset, realT color[3]) {
02104 int k, n= qh_setsize(points), i;
02105 pointT *point, **pointp;
02106 setT *printpoints;
02107
02108 qh_fprintf(fp, 9098, "{ OFF %d 1 1 # f%d\n", n, facet->id);
02109 if (offset != 0.0) {
02110 printpoints= qh_settemp(n);
02111 FOREACHpoint_(points)
02112 qh_setappend(&printpoints, qh_projectpoint(point, facet, -offset));
02113 }else
02114 printpoints= points;
02115 FOREACHpoint_(printpoints) {
02116 for (k=0; k < qh hull_dim; k++) {
02117 if (k == qh DROPdim)
02118 qh_fprintf(fp, 9099, "0 ");
02119 else
02120 qh_fprintf(fp, 9100, "%8.4g ", point[k]);
02121 }
02122 if (printpoints != points)
02123 qh_memfree(point, qh normal_size);
02124 qh_fprintf(fp, 9101, "\n");
02125 }
02126 if (printpoints != points)
02127 qh_settempfree(&printpoints);
02128 qh_fprintf(fp, 9102, "%d ", n);
02129 for (i=0; i < n; i++)
02130 qh_fprintf(fp, 9103, "%d ", i);
02131 qh_fprintf(fp, 9104, "%8.4g %8.4g %8.4g 1.0 }\n", color[0], color[1], color[2]);
02132 }
02133
02134
02135
02136
02137
02138
02139
02140
02141
02142
02143
02144
02145
02146
02147
02148
02149 void qh_printfacet3geom_simplicial(FILE *fp, facetT *facet, realT color[3]) {
02150 setT *points, *vertices;
02151 vertexT *vertex, **vertexp, *vertexA, *vertexB;
02152 facetT *neighbor, **neighborp;
02153 realT outerplane, innerplane;
02154 realT black[3]={0, 0, 0}, green[3]={0, 1, 0};
02155 int k;
02156
02157 qh_geomplanes(facet, &outerplane, &innerplane);
02158 vertices= qh_facet3vertex(facet);
02159 points= qh_settemp(qh TEMPsize);
02160 FOREACHvertex_(vertices)
02161 qh_setappend(&points, vertex->point);
02162 if (qh PRINTouter || (!qh PRINTnoplanes && !qh PRINTinner))
02163 qh_printfacet3geom_points(fp, points, facet, outerplane, color);
02164 if (qh PRINTinner || (!qh PRINTnoplanes && !qh PRINTouter &&
02165 outerplane - innerplane > 2 * qh MAXabs_coord * qh_GEOMepsilon)) {
02166 for (k=3; k--; )
02167 color[k]= 1.0 - color[k];
02168 qh_printfacet3geom_points(fp, points, facet, innerplane, color);
02169 }
02170 qh_settempfree(&points);
02171 qh_settempfree(&vertices);
02172 if ((qh DOintersections || qh PRINTridges)
02173 && (!facet->visible || !qh NEWfacets)) {
02174 facet->visitid= qh visit_id;
02175 FOREACHneighbor_(facet) {
02176 if (neighbor->visitid != qh visit_id) {
02177 vertices= qh_setnew_delnthsorted(facet->vertices, qh hull_dim,
02178 SETindex_(facet->neighbors, neighbor), 0);
02179 if (qh DOintersections)
02180 qh_printhyperplaneintersection(fp, facet, neighbor, vertices, black);
02181 if (qh PRINTridges) {
02182 vertexA= SETfirstt_(vertices, vertexT);
02183 vertexB= SETsecondt_(vertices, vertexT);
02184 qh_printline3geom(fp, vertexA->point, vertexB->point, green);
02185 }
02186 qh_setfree(&vertices);
02187 }
02188 }
02189 }
02190 }
02191
02192
02193
02194
02195
02196
02197
02198
02199
02200
02201
02202
02203 void qh_printfacet3math(FILE *fp, facetT *facet, qh_PRINT format, int notfirst) {
02204 vertexT *vertex, **vertexp;
02205 setT *points, *vertices;
02206 pointT *point, **pointp;
02207 boolT firstpoint= True;
02208 realT dist;
02209 const char *pointfmt, *endfmt;
02210
02211 if (notfirst)
02212 qh_fprintf(fp, 9105, ",\n");
02213 vertices= qh_facet3vertex(facet);
02214 points= qh_settemp(qh_setsize(vertices));
02215 FOREACHvertex_(vertices) {
02216 zinc_(Zdistio);
02217 qh_distplane(vertex->point, facet, &dist);
02218 point= qh_projectpoint(vertex->point, facet, dist);
02219 qh_setappend(&points, point);
02220 }
02221 if (format == qh_PRINTmaple) {
02222 qh_fprintf(fp, 9106, "[");
02223 pointfmt= "[%16.8f, %16.8f, %16.8f]";
02224 endfmt= "]";
02225 }else {
02226 qh_fprintf(fp, 9107, "Polygon[{");
02227 pointfmt= "{%16.8f, %16.8f, %16.8f}";
02228 endfmt= "}]";
02229 }
02230 FOREACHpoint_(points) {
02231 if (firstpoint)
02232 firstpoint= False;
02233 else
02234 qh_fprintf(fp, 9108, ",\n");
02235 qh_fprintf(fp, 9109, pointfmt, point[0], point[1], point[2]);
02236 }
02237 FOREACHpoint_(points)
02238 qh_memfree(point, qh normal_size);
02239 qh_settempfree(&points);
02240 qh_settempfree(&vertices);
02241 qh_fprintf(fp, 9110, endfmt);
02242 }
02243
02244
02245
02246
02247
02248
02249
02250
02251
02252
02253
02254
02255 void qh_printfacet3vertex(FILE *fp, facetT *facet, qh_PRINT format) {
02256 vertexT *vertex, **vertexp;
02257 setT *vertices;
02258
02259 vertices= qh_facet3vertex(facet);
02260 if (format == qh_PRINToff)
02261 qh_fprintf(fp, 9111, "%d ", qh_setsize(vertices));
02262 FOREACHvertex_(vertices)
02263 qh_fprintf(fp, 9112, "%d ", qh_pointid(vertex->point));
02264 qh_fprintf(fp, 9113, "\n");
02265 qh_settempfree(&vertices);
02266 }
02267
02268
02269
02270
02271
02272
02273
02274
02275
02276
02277
02278
02279
02280
02281 void qh_printfacet4geom_nonsimplicial(FILE *fp, facetT *facet, realT color[3]) {
02282 facetT *neighbor;
02283 ridgeT *ridge, **ridgep;
02284 vertexT *vertex, **vertexp;
02285 pointT *point;
02286 int k;
02287 realT dist;
02288
02289 facet->visitid= qh visit_id;
02290 if (qh PRINTnoplanes || (facet->visible && qh NEWfacets))
02291 return;
02292 FOREACHridge_(facet->ridges) {
02293 neighbor= otherfacet_(ridge, facet);
02294 if (neighbor->visitid == qh visit_id)
02295 continue;
02296 if (qh PRINTtransparent && !neighbor->good)
02297 continue;
02298 if (qh DOintersections)
02299 qh_printhyperplaneintersection(fp, facet, neighbor, ridge->vertices, color);
02300 else {
02301 if (qh DROPdim >= 0)
02302 qh_fprintf(fp, 9114, "OFF 3 1 1 # f%d\n", facet->id);
02303 else {
02304 qh printoutvar++;
02305 qh_fprintf(fp, 9115, "# r%d between f%d f%d\n", ridge->id, facet->id, neighbor->id);
02306 }
02307 FOREACHvertex_(ridge->vertices) {
02308 zinc_(Zdistio);
02309 qh_distplane(vertex->point,facet, &dist);
02310 point=qh_projectpoint(vertex->point,facet, dist);
02311 for (k=0; k < qh hull_dim; k++) {
02312 if (k != qh DROPdim)
02313 qh_fprintf(fp, 9116, "%8.4g ", point[k]);
02314 }
02315 qh_fprintf(fp, 9117, "\n");
02316 qh_memfree(point, qh normal_size);
02317 }
02318 if (qh DROPdim >= 0)
02319 qh_fprintf(fp, 9118, "3 0 1 2 %8.4g %8.4g %8.4g\n", color[0], color[1], color[2]);
02320 }
02321 }
02322 }
02323
02324
02325
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335 void qh_printfacet4geom_simplicial(FILE *fp, facetT *facet, realT color[3]) {
02336 setT *vertices;
02337 facetT *neighbor, **neighborp;
02338 vertexT *vertex, **vertexp;
02339 int k;
02340
02341 facet->visitid= qh visit_id;
02342 if (qh PRINTnoplanes || (facet->visible && qh NEWfacets))
02343 return;
02344 FOREACHneighbor_(facet) {
02345 if (neighbor->visitid == qh visit_id)
02346 continue;
02347 if (qh PRINTtransparent && !neighbor->good)
02348 continue;
02349 vertices= qh_setnew_delnthsorted(facet->vertices, qh hull_dim,
02350 SETindex_(facet->neighbors, neighbor), 0);
02351 if (qh DOintersections)
02352 qh_printhyperplaneintersection(fp, facet, neighbor, vertices, color);
02353 else {
02354 if (qh DROPdim >= 0)
02355 qh_fprintf(fp, 9119, "OFF 3 1 1 # ridge between f%d f%d\n",
02356 facet->id, neighbor->id);
02357 else {
02358 qh printoutvar++;
02359 qh_fprintf(fp, 9120, "# ridge between f%d f%d\n", facet->id, neighbor->id);
02360 }
02361 FOREACHvertex_(vertices) {
02362 for (k=0; k < qh hull_dim; k++) {
02363 if (k != qh DROPdim)
02364 qh_fprintf(fp, 9121, "%8.4g ", vertex->point[k]);
02365 }
02366 qh_fprintf(fp, 9122, "\n");
02367 }
02368 if (qh DROPdim >= 0)
02369 qh_fprintf(fp, 9123, "3 0 1 2 %8.4g %8.4g %8.4g\n", color[0], color[1], color[2]);
02370 }
02371 qh_setfree(&vertices);
02372 }
02373 }
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383 void qh_printfacetNvertex_nonsimplicial(FILE *fp, facetT *facet, int id, qh_PRINT format) {
02384 vertexT *vertex, **vertexp;
02385 ridgeT *ridge, **ridgep;
02386
02387 if (facet->visible && qh NEWfacets)
02388 return;
02389 FOREACHridge_(facet->ridges) {
02390 if (format == qh_PRINTtriangles)
02391 qh_fprintf(fp, 9124, "%d ", qh hull_dim);
02392 qh_fprintf(fp, 9125, "%d ", id);
02393 if ((ridge->top == facet) ^ qh_ORIENTclock) {
02394 FOREACHvertex_(ridge->vertices)
02395 qh_fprintf(fp, 9126, "%d ", qh_pointid(vertex->point));
02396 }else {
02397 FOREACHvertexreverse12_(ridge->vertices)
02398 qh_fprintf(fp, 9127, "%d ", qh_pointid(vertex->point));
02399 }
02400 qh_fprintf(fp, 9128, "\n");
02401 }
02402 }
02403
02404
02405
02406
02407
02408
02409
02410
02411
02412
02413
02414 void qh_printfacetNvertex_simplicial(FILE *fp, facetT *facet, qh_PRINT format) {
02415 vertexT *vertex, **vertexp;
02416
02417 if (format == qh_PRINToff || format == qh_PRINTtriangles)
02418 qh_fprintf(fp, 9129, "%d ", qh_setsize(facet->vertices));
02419 if ((facet->toporient ^ qh_ORIENTclock)
02420 || (qh hull_dim > 2 && !facet->simplicial)) {
02421 FOREACHvertex_(facet->vertices)
02422 qh_fprintf(fp, 9130, "%d ", qh_pointid(vertex->point));
02423 }else {
02424 FOREACHvertexreverse12_(facet->vertices)
02425 qh_fprintf(fp, 9131, "%d ", qh_pointid(vertex->point));
02426 }
02427 qh_fprintf(fp, 9132, "\n");
02428 }
02429
02430
02431
02432
02433
02434
02435
02436
02437
02438
02439
02440
02441 void qh_printfacetheader(FILE *fp, facetT *facet) {
02442 pointT *point, **pointp, *furthest;
02443 facetT *neighbor, **neighborp;
02444 realT dist;
02445
02446 if (facet == qh_MERGEridge) {
02447 qh_fprintf(fp, 9133, " MERGEridge\n");
02448 return;
02449 }else if (facet == qh_DUPLICATEridge) {
02450 qh_fprintf(fp, 9134, " DUPLICATEridge\n");
02451 return;
02452 }else if (!facet) {
02453 qh_fprintf(fp, 9135, " NULLfacet\n");
02454 return;
02455 }
02456 qh old_randomdist= qh RANDOMdist;
02457 qh RANDOMdist= False;
02458 qh_fprintf(fp, 9136, "- f%d\n", facet->id);
02459 qh_fprintf(fp, 9137, " - flags:");
02460 if (facet->toporient)
02461 qh_fprintf(fp, 9138, " top");
02462 else
02463 qh_fprintf(fp, 9139, " bottom");
02464 if (facet->simplicial)
02465 qh_fprintf(fp, 9140, " simplicial");
02466 if (facet->tricoplanar)
02467 qh_fprintf(fp, 9141, " tricoplanar");
02468 if (facet->upperdelaunay)
02469 qh_fprintf(fp, 9142, " upperDelaunay");
02470 if (facet->visible)
02471 qh_fprintf(fp, 9143, " visible");
02472 if (facet->newfacet)
02473 qh_fprintf(fp, 9144, " new");
02474 if (facet->tested)
02475 qh_fprintf(fp, 9145, " tested");
02476 if (!facet->good)
02477 qh_fprintf(fp, 9146, " notG");
02478 if (facet->seen)
02479 qh_fprintf(fp, 9147, " seen");
02480 if (facet->coplanar)
02481 qh_fprintf(fp, 9148, " coplanar");
02482 if (facet->mergehorizon)
02483 qh_fprintf(fp, 9149, " mergehorizon");
02484 if (facet->keepcentrum)
02485 qh_fprintf(fp, 9150, " keepcentrum");
02486 if (facet->dupridge)
02487 qh_fprintf(fp, 9151, " dupridge");
02488 if (facet->mergeridge && !facet->mergeridge2)
02489 qh_fprintf(fp, 9152, " mergeridge1");
02490 if (facet->mergeridge2)
02491 qh_fprintf(fp, 9153, " mergeridge2");
02492 if (facet->newmerge)
02493 qh_fprintf(fp, 9154, " newmerge");
02494 if (facet->flipped)
02495 qh_fprintf(fp, 9155, " flipped");
02496 if (facet->notfurthest)
02497 qh_fprintf(fp, 9156, " notfurthest");
02498 if (facet->degenerate)
02499 qh_fprintf(fp, 9157, " degenerate");
02500 if (facet->redundant)
02501 qh_fprintf(fp, 9158, " redundant");
02502 qh_fprintf(fp, 9159, "\n");
02503 if (facet->isarea)
02504 qh_fprintf(fp, 9160, " - area: %2.2g\n", facet->f.area);
02505 else if (qh NEWfacets && facet->visible && facet->f.replace)
02506 qh_fprintf(fp, 9161, " - replacement: f%d\n", facet->f.replace->id);
02507 else if (facet->newfacet) {
02508 if (facet->f.samecycle && facet->f.samecycle != facet)
02509 qh_fprintf(fp, 9162, " - shares same visible/horizon as f%d\n", facet->f.samecycle->id);
02510 }else if (facet->tricoplanar ) {
02511 if (facet->f.triowner)
02512 qh_fprintf(fp, 9163, " - owner of normal & centrum is facet f%d\n", facet->f.triowner->id);
02513 }else if (facet->f.newcycle)
02514 qh_fprintf(fp, 9164, " - was horizon to f%d\n", facet->f.newcycle->id);
02515 if (facet->nummerge)
02516 qh_fprintf(fp, 9165, " - merges: %d\n", facet->nummerge);
02517 qh_printpointid(fp, " - normal: ", qh hull_dim, facet->normal, -1);
02518 qh_fprintf(fp, 9166, " - offset: %10.7g\n", facet->offset);
02519 if (qh CENTERtype == qh_ASvoronoi || facet->center)
02520 qh_printcenter(fp, qh_PRINTfacets, " - center: ", facet);
02521 #if qh_MAXoutside
02522 if (facet->maxoutside > qh DISTround)
02523 qh_fprintf(fp, 9167, " - maxoutside: %10.7g\n", facet->maxoutside);
02524 #endif
02525 if (!SETempty_(facet->outsideset)) {
02526 furthest= (pointT*)qh_setlast(facet->outsideset);
02527 if (qh_setsize(facet->outsideset) < 6) {
02528 qh_fprintf(fp, 9168, " - outside set(furthest p%d):\n", qh_pointid(furthest));
02529 FOREACHpoint_(facet->outsideset)
02530 qh_printpoint(fp, " ", point);
02531 }else if (qh_setsize(facet->outsideset) < 21) {
02532 qh_printpoints(fp, " - outside set:", facet->outsideset);
02533 }else {
02534 qh_fprintf(fp, 9169, " - outside set: %d points.", qh_setsize(facet->outsideset));
02535 qh_printpoint(fp, " Furthest", furthest);
02536 }
02537 #if !qh_COMPUTEfurthest
02538 qh_fprintf(fp, 9170, " - furthest distance= %2.2g\n", facet->furthestdist);
02539 #endif
02540 }
02541 if (!SETempty_(facet->coplanarset)) {
02542 furthest= (pointT*)qh_setlast(facet->coplanarset);
02543 if (qh_setsize(facet->coplanarset) < 6) {
02544 qh_fprintf(fp, 9171, " - coplanar set(furthest p%d):\n", qh_pointid(furthest));
02545 FOREACHpoint_(facet->coplanarset)
02546 qh_printpoint(fp, " ", point);
02547 }else if (qh_setsize(facet->coplanarset) < 21) {
02548 qh_printpoints(fp, " - coplanar set:", facet->coplanarset);
02549 }else {
02550 qh_fprintf(fp, 9172, " - coplanar set: %d points.", qh_setsize(facet->coplanarset));
02551 qh_printpoint(fp, " Furthest", furthest);
02552 }
02553 zinc_(Zdistio);
02554 qh_distplane(furthest, facet, &dist);
02555 qh_fprintf(fp, 9173, " furthest distance= %2.2g\n", dist);
02556 }
02557 qh_printvertices(fp, " - vertices:", facet->vertices);
02558 qh_fprintf(fp, 9174, " - neighboring facets:");
02559 FOREACHneighbor_(facet) {
02560 if (neighbor == qh_MERGEridge)
02561 qh_fprintf(fp, 9175, " MERGE");
02562 else if (neighbor == qh_DUPLICATEridge)
02563 qh_fprintf(fp, 9176, " DUP");
02564 else
02565 qh_fprintf(fp, 9177, " f%d", neighbor->id);
02566 }
02567 qh_fprintf(fp, 9178, "\n");
02568 qh RANDOMdist= qh old_randomdist;
02569 }
02570
02571
02572
02573
02574
02575
02576
02577
02578
02579
02580
02581
02582
02583
02584 void qh_printfacetridges(FILE *fp, facetT *facet) {
02585 facetT *neighbor, **neighborp;
02586 ridgeT *ridge, **ridgep;
02587 int numridges= 0;
02588
02589
02590 if (facet->visible && qh NEWfacets) {
02591 qh_fprintf(fp, 9179, " - ridges(ids may be garbage):");
02592 FOREACHridge_(facet->ridges)
02593 qh_fprintf(fp, 9180, " r%d", ridge->id);
02594 qh_fprintf(fp, 9181, "\n");
02595 }else {
02596 qh_fprintf(fp, 9182, " - ridges:\n");
02597 FOREACHridge_(facet->ridges)
02598 ridge->seen= False;
02599 if (qh hull_dim == 3) {
02600 ridge= SETfirstt_(facet->ridges, ridgeT);
02601 while (ridge && !ridge->seen) {
02602 ridge->seen= True;
02603 qh_printridge(fp, ridge);
02604 numridges++;
02605 ridge= qh_nextridge3d(ridge, facet, NULL);
02606 }
02607 }else {
02608 FOREACHneighbor_(facet) {
02609 FOREACHridge_(facet->ridges) {
02610 if (otherfacet_(ridge,facet) == neighbor) {
02611 ridge->seen= True;
02612 qh_printridge(fp, ridge);
02613 numridges++;
02614 }
02615 }
02616 }
02617 }
02618 if (numridges != qh_setsize(facet->ridges)) {
02619 qh_fprintf(fp, 9183, " - all ridges:");
02620 FOREACHridge_(facet->ridges)
02621 qh_fprintf(fp, 9184, " r%d", ridge->id);
02622 qh_fprintf(fp, 9185, "\n");
02623 }
02624 FOREACHridge_(facet->ridges) {
02625 if (!ridge->seen)
02626 qh_printridge(fp, ridge);
02627 }
02628 }
02629 }
02630
02631
02632
02633
02634
02635
02636
02637
02638
02639
02640
02641 void qh_printfacets(FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall) {
02642 int numfacets, numsimplicial, numridges, totneighbors, numcoplanars, numtricoplanars;
02643 facetT *facet, **facetp;
02644 setT *vertices;
02645 coordT *center;
02646 realT outerplane, innerplane;
02647
02648 qh old_randomdist= qh RANDOMdist;
02649 qh RANDOMdist= False;
02650 if (qh CDDoutput && (format == qh_PRINTcentrums || format == qh_PRINTpointintersect || format == qh_PRINToff))
02651 qh_fprintf(qh ferr, 7056, "qhull warning: CDD format is not available for centrums, halfspace\nintersections, and OFF file format.\n");
02652 if (format == qh_PRINTnone)
02653 ;
02654 else if (format == qh_PRINTaverage) {
02655 vertices= qh_facetvertices(facetlist, facets, printall);
02656 center= qh_getcenter(vertices);
02657 qh_fprintf(fp, 9186, "%d 1\n", qh hull_dim);
02658 qh_printpointid(fp, NULL, qh hull_dim, center, -1);
02659 qh_memfree(center, qh normal_size);
02660 qh_settempfree(&vertices);
02661 }else if (format == qh_PRINTextremes) {
02662 if (qh DELAUNAY)
02663 qh_printextremes_d(fp, facetlist, facets, printall);
02664 else if (qh hull_dim == 2)
02665 qh_printextremes_2d(fp, facetlist, facets, printall);
02666 else
02667 qh_printextremes(fp, facetlist, facets, printall);
02668 }else if (format == qh_PRINToptions)
02669 qh_fprintf(fp, 9187, "Options selected for Qhull %s:\n%s\n", qh_version, qh qhull_options);
02670 else if (format == qh_PRINTpoints && !qh VORONOI)
02671 qh_printpoints_out(fp, facetlist, facets, printall);
02672 else if (format == qh_PRINTqhull)
02673 qh_fprintf(fp, 9188, "%s | %s\n", qh rbox_command, qh qhull_command);
02674 else if (format == qh_PRINTsize) {
02675 qh_fprintf(fp, 9189, "0\n2 ");
02676 qh_fprintf(fp, 9190, qh_REAL_1, qh totarea);
02677 qh_fprintf(fp, 9191, qh_REAL_1, qh totvol);
02678 qh_fprintf(fp, 9192, "\n");
02679 }else if (format == qh_PRINTsummary) {
02680 qh_countfacets(facetlist, facets, printall, &numfacets, &numsimplicial,
02681 &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
02682 vertices= qh_facetvertices(facetlist, facets, printall);
02683 qh_fprintf(fp, 9193, "10 %d %d %d %d %d %d %d %d %d %d\n2 ", qh hull_dim,
02684 qh num_points + qh_setsize(qh other_points),
02685 qh num_vertices, qh num_facets - qh num_visible,
02686 qh_setsize(vertices), numfacets, numcoplanars,
02687 numfacets - numsimplicial, zzval_(Zdelvertextot),
02688 numtricoplanars);
02689 qh_settempfree(&vertices);
02690 qh_outerinner(NULL, &outerplane, &innerplane);
02691 qh_fprintf(fp, 9194, qh_REAL_2n, outerplane, innerplane);
02692 }else if (format == qh_PRINTvneighbors)
02693 qh_printvneighbors(fp, facetlist, facets, printall);
02694 else if (qh VORONOI && format == qh_PRINToff)
02695 qh_printvoronoi(fp, format, facetlist, facets, printall);
02696 else if (qh VORONOI && format == qh_PRINTgeom) {
02697 qh_printbegin(fp, format, facetlist, facets, printall);
02698 qh_printvoronoi(fp, format, facetlist, facets, printall);
02699 qh_printend(fp, format, facetlist, facets, printall);
02700 }else if (qh VORONOI
02701 && (format == qh_PRINTvertices || format == qh_PRINTinner || format == qh_PRINTouter))
02702 qh_printvdiagram(fp, format, facetlist, facets, printall);
02703 else {
02704 qh_printbegin(fp, format, facetlist, facets, printall);
02705 FORALLfacet_(facetlist)
02706 qh_printafacet(fp, format, facet, printall);
02707 FOREACHfacet_(facets)
02708 qh_printafacet(fp, format, facet, printall);
02709 qh_printend(fp, format, facetlist, facets, printall);
02710 }
02711 qh RANDOMdist= qh old_randomdist;
02712 }
02713
02714
02715
02716
02717
02718
02719
02720
02721 void qh_printhyperplaneintersection(FILE *fp, facetT *facet1, facetT *facet2,
02722 setT *vertices, realT color[3]) {
02723 realT costheta, denominator, dist1, dist2, s, t, mindenom, p[4];
02724 vertexT *vertex, **vertexp;
02725 int i, k;
02726 boolT nearzero1, nearzero2;
02727
02728 costheta= qh_getangle(facet1->normal, facet2->normal);
02729 denominator= 1 - costheta * costheta;
02730 i= qh_setsize(vertices);
02731 if (qh hull_dim == 3)
02732 qh_fprintf(fp, 9195, "VECT 1 %d 1 %d 1 ", i, i);
02733 else if (qh hull_dim == 4 && qh DROPdim >= 0)
02734 qh_fprintf(fp, 9196, "OFF 3 1 1 ");
02735 else
02736 qh printoutvar++;
02737 qh_fprintf(fp, 9197, "# intersect f%d f%d\n", facet1->id, facet2->id);
02738 mindenom= 1 / (10.0 * qh MAXabs_coord);
02739 FOREACHvertex_(vertices) {
02740 zadd_(Zdistio, 2);
02741 qh_distplane(vertex->point, facet1, &dist1);
02742 qh_distplane(vertex->point, facet2, &dist2);
02743 s= qh_divzero(-dist1 + costheta * dist2, denominator,mindenom,&nearzero1);
02744 t= qh_divzero(-dist2 + costheta * dist1, denominator,mindenom,&nearzero2);
02745 if (nearzero1 || nearzero2)
02746 s= t= 0.0;
02747 for (k=qh hull_dim; k--; )
02748 p[k]= vertex->point[k] + facet1->normal[k] * s + facet2->normal[k] * t;
02749 if (qh PRINTdim <= 3) {
02750 qh_projectdim3 (p, p);
02751 qh_fprintf(fp, 9198, "%8.4g %8.4g %8.4g # ", p[0], p[1], p[2]);
02752 }else
02753 qh_fprintf(fp, 9199, "%8.4g %8.4g %8.4g %8.4g # ", p[0], p[1], p[2], p[3]);
02754 if (nearzero1+nearzero2)
02755 qh_fprintf(fp, 9200, "p%d(coplanar facets)\n", qh_pointid(vertex->point));
02756 else
02757 qh_fprintf(fp, 9201, "projected p%d\n", qh_pointid(vertex->point));
02758 }
02759 if (qh hull_dim == 3)
02760 qh_fprintf(fp, 9202, "%8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]);
02761 else if (qh hull_dim == 4 && qh DROPdim >= 0)
02762 qh_fprintf(fp, 9203, "3 0 1 2 %8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]);
02763 }
02764
02765
02766
02767
02768
02769
02770
02771
02772
02773
02774
02775
02776 void qh_printline3geom(FILE *fp, pointT *pointA, pointT *pointB, realT color[3]) {
02777 int k;
02778 realT pA[4], pB[4];
02779
02780 qh_projectdim3(pointA, pA);
02781 qh_projectdim3(pointB, pB);
02782 if ((fabs(pA[0] - pB[0]) > 1e-3) ||
02783 (fabs(pA[1] - pB[1]) > 1e-3) ||
02784 (fabs(pA[2] - pB[2]) > 1e-3)) {
02785 qh_fprintf(fp, 9204, "VECT 1 2 1 2 1\n");
02786 for (k=0; k < 3; k++)
02787 qh_fprintf(fp, 9205, "%8.4g ", pB[k]);
02788 qh_fprintf(fp, 9206, " # p%d\n", qh_pointid(pointB));
02789 }else
02790 qh_fprintf(fp, 9207, "VECT 1 1 1 1 1\n");
02791 for (k=0; k < 3; k++)
02792 qh_fprintf(fp, 9208, "%8.4g ", pA[k]);
02793 qh_fprintf(fp, 9209, " # p%d\n", qh_pointid(pointA));
02794 qh_fprintf(fp, 9210, "%8.4g %8.4g %8.4g 1\n", color[0], color[1], color[2]);
02795 }
02796
02797
02798
02799
02800
02801
02802
02803
02804
02805
02806
02807 void qh_printneighborhood(FILE *fp, qh_PRINT format, facetT *facetA, facetT *facetB, boolT printall) {
02808 facetT *neighbor, **neighborp, *facet;
02809 setT *facets;
02810
02811 if (format == qh_PRINTnone)
02812 return;
02813 qh_findgood_all(qh facet_list);
02814 if (facetA == facetB)
02815 facetB= NULL;
02816 facets= qh_settemp(2*(qh_setsize(facetA->neighbors)+1));
02817 qh visit_id++;
02818 for (facet= facetA; facet; facet= ((facet == facetA) ? facetB : NULL)) {
02819 if (facet->visitid != qh visit_id) {
02820 facet->visitid= qh visit_id;
02821 qh_setappend(&facets, facet);
02822 }
02823 FOREACHneighbor_(facet) {
02824 if (neighbor->visitid == qh visit_id)
02825 continue;
02826 neighbor->visitid= qh visit_id;
02827 if (printall || !qh_skipfacet(neighbor))
02828 qh_setappend(&facets, neighbor);
02829 }
02830 }
02831 qh_printfacets(fp, format, NULL, facets, printall);
02832 qh_settempfree(&facets);
02833 }
02834
02835
02836
02837
02838
02839
02840
02841
02842
02843
02844
02845
02846
02847
02848
02849
02850
02851 void qh_printpoint(FILE *fp, const char *string, pointT *point) {
02852 int id= qh_pointid( point);
02853
02854 qh_printpointid( fp, string, qh hull_dim, point, id);
02855 }
02856
02857 void qh_printpointid(FILE *fp, const char *string, int dim, pointT *point, int id) {
02858 int k;
02859 realT r;
02860
02861 if (!point)
02862 return;
02863 if (string) {
02864 qh_fprintf(fp, 9211, "%s", string);
02865 if (id != -1)
02866 qh_fprintf(fp, 9212, " p%d: ", id);
02867 }
02868 for (k=dim; k--; ) {
02869 r= *point++;
02870 if (string)
02871 qh_fprintf(fp, 9213, " %8.4g", r);
02872 else
02873 qh_fprintf(fp, 9214, qh_REAL_1, r);
02874 }
02875 qh_fprintf(fp, 9215, "\n");
02876 }
02877
02878
02879
02880
02881
02882
02883
02884 void qh_printpoint3 (FILE *fp, pointT *point) {
02885 int k;
02886 realT p[4];
02887
02888 qh_projectdim3 (point, p);
02889 for (k=0; k < 3; k++)
02890 qh_fprintf(fp, 9216, "%8.4g ", p[k]);
02891 qh_fprintf(fp, 9217, " # p%d\n", qh_pointid(point));
02892 }
02893
02894
02895
02896
02897
02898
02899
02900
02901
02902
02903
02904
02905
02906
02907
02908
02909
02910
02911 void qh_printpoints_out(FILE *fp, facetT *facetlist, setT *facets, boolT printall) {
02912 int allpoints= qh num_points + qh_setsize(qh other_points);
02913 int numpoints=0, point_i, point_n;
02914 setT *vertices, *points;
02915 facetT *facet, **facetp;
02916 pointT *point, **pointp;
02917 vertexT *vertex, **vertexp;
02918 int id;
02919
02920 points= qh_settemp(allpoints);
02921 qh_setzero(points, 0, allpoints);
02922 vertices= qh_facetvertices(facetlist, facets, printall);
02923 FOREACHvertex_(vertices) {
02924 id= qh_pointid(vertex->point);
02925 if (id >= 0)
02926 SETelem_(points, id)= vertex->point;
02927 }
02928 if (qh KEEPinside || qh KEEPcoplanar || qh KEEPnearinside) {
02929 FORALLfacet_(facetlist) {
02930 if (!printall && qh_skipfacet(facet))
02931 continue;
02932 FOREACHpoint_(facet->coplanarset) {
02933 id= qh_pointid(point);
02934 if (id >= 0)
02935 SETelem_(points, id)= point;
02936 }
02937 }
02938 FOREACHfacet_(facets) {
02939 if (!printall && qh_skipfacet(facet))
02940 continue;
02941 FOREACHpoint_(facet->coplanarset) {
02942 id= qh_pointid(point);
02943 if (id >= 0)
02944 SETelem_(points, id)= point;
02945 }
02946 }
02947 }
02948 qh_settempfree(&vertices);
02949 FOREACHpoint_i_(points) {
02950 if (point)
02951 numpoints++;
02952 }
02953 if (qh CDDoutput)
02954 qh_fprintf(fp, 9218, "%s | %s\nbegin\n%d %d real\n", qh rbox_command,
02955 qh qhull_command, numpoints, qh hull_dim + 1);
02956 else
02957 qh_fprintf(fp, 9219, "%d\n%d\n", qh hull_dim, numpoints);
02958 FOREACHpoint_i_(points) {
02959 if (point) {
02960 if (qh CDDoutput)
02961 qh_fprintf(fp, 9220, "1 ");
02962 qh_printpoint(fp, NULL, point);
02963 }
02964 }
02965 if (qh CDDoutput)
02966 qh_fprintf(fp, 9221, "end\n");
02967 qh_settempfree(&points);
02968 }
02969
02970
02971
02972
02973
02974
02975
02976
02977 void qh_printpointvect(FILE *fp, pointT *point, coordT *normal, pointT *center, realT radius, realT color[3]) {
02978 realT diff[4], pointA[4];
02979 int k;
02980
02981 for (k=qh hull_dim; k--; ) {
02982 if (center)
02983 diff[k]= point[k]-center[k];
02984 else if (normal)
02985 diff[k]= normal[k];
02986 else
02987 diff[k]= 0;
02988 }
02989 if (center)
02990 qh_normalize2 (diff, qh hull_dim, True, NULL, NULL);
02991 for (k=qh hull_dim; k--; )
02992 pointA[k]= point[k]+diff[k] * radius;
02993 qh_printline3geom(fp, point, pointA, color);
02994 }
02995
02996
02997
02998
02999
03000
03001
03002 void qh_printpointvect2 (FILE *fp, pointT *point, coordT *normal, pointT *center, realT radius) {
03003 realT red[3]={1, 0, 0}, yellow[3]={1, 1, 0};
03004
03005 qh_printpointvect(fp, point, normal, center, radius, red);
03006 qh_printpointvect(fp, point, normal, center, -radius, yellow);
03007 }
03008
03009
03010
03011
03012
03013
03014
03015
03016
03017
03018
03019 void qh_printridge(FILE *fp, ridgeT *ridge) {
03020
03021 qh_fprintf(fp, 9222, " - r%d", ridge->id);
03022 if (ridge->tested)
03023 qh_fprintf(fp, 9223, " tested");
03024 if (ridge->nonconvex)
03025 qh_fprintf(fp, 9224, " nonconvex");
03026 qh_fprintf(fp, 9225, "\n");
03027 qh_printvertices(fp, " vertices:", ridge->vertices);
03028 if (ridge->top && ridge->bottom)
03029 qh_fprintf(fp, 9226, " between f%d and f%d\n",
03030 ridge->top->id, ridge->bottom->id);
03031 }
03032
03033
03034
03035
03036
03037
03038
03039
03040
03041
03042 void qh_printspheres(FILE *fp, setT *vertices, realT radius) {
03043 vertexT *vertex, **vertexp;
03044
03045 qh printoutnum++;
03046 qh_fprintf(fp, 9227, "{appearance {-edge -normal normscale 0} {\n\
03047 INST geom {define vsphere OFF\n\
03048 18 32 48\n\
03049 \n\
03050 0 0 1\n\
03051 1 0 0\n\
03052 0 1 0\n\
03053 -1 0 0\n\
03054 0 -1 0\n\
03055 0 0 -1\n\
03056 0.707107 0 0.707107\n\
03057 0 -0.707107 0.707107\n\
03058 0.707107 -0.707107 0\n\
03059 -0.707107 0 0.707107\n\
03060 -0.707107 -0.707107 0\n\
03061 0 0.707107 0.707107\n\
03062 -0.707107 0.707107 0\n\
03063 0.707107 0.707107 0\n\
03064 0.707107 0 -0.707107\n\
03065 0 0.707107 -0.707107\n\
03066 -0.707107 0 -0.707107\n\
03067 0 -0.707107 -0.707107\n\
03068 \n\
03069 3 0 6 11\n\
03070 3 0 7 6 \n\
03071 3 0 9 7 \n\
03072 3 0 11 9\n\
03073 3 1 6 8 \n\
03074 3 1 8 14\n\
03075 3 1 13 6\n\
03076 3 1 14 13\n\
03077 3 2 11 13\n\
03078 3 2 12 11\n\
03079 3 2 13 15\n\
03080 3 2 15 12\n\
03081 3 3 9 12\n\
03082 3 3 10 9\n\
03083 3 3 12 16\n\
03084 3 3 16 10\n\
03085 3 4 7 10\n\
03086 3 4 8 7\n\
03087 3 4 10 17\n\
03088 3 4 17 8\n\
03089 3 5 14 17\n\
03090 3 5 15 14\n\
03091 3 5 16 15\n\
03092 3 5 17 16\n\
03093 3 6 13 11\n\
03094 3 7 8 6\n\
03095 3 9 10 7\n\
03096 3 11 12 9\n\
03097 3 14 8 17\n\
03098 3 15 13 14\n\
03099 3 16 12 15\n\
03100 3 17 10 16\n} transforms { TLIST\n");
03101 FOREACHvertex_(vertices) {
03102 qh_fprintf(fp, 9228, "%8.4g 0 0 0 # v%d\n 0 %8.4g 0 0\n0 0 %8.4g 0\n",
03103 radius, vertex->id, radius, radius);
03104 qh_printpoint3 (fp, vertex->point);
03105 qh_fprintf(fp, 9229, "1\n");
03106 }
03107 qh_fprintf(fp, 9230, "}}}\n");
03108 }
03109
03110
03111
03112
03113
03114
03115
03116
03117
03118
03119
03120
03121
03122
03123
03124
03125
03126
03127
03128
03129
03130
03131
03132
03133
03134
03135
03136
03137
03138 void qh_printvdiagram(FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall) {
03139 setT *vertices;
03140 int totcount, numcenters;
03141 boolT isLower;
03142 qh_RIDGE innerouter= qh_RIDGEall;
03143 printvridgeT printvridge= NULL;
03144
03145 if (format == qh_PRINTvertices) {
03146 innerouter= qh_RIDGEall;
03147 printvridge= qh_printvridge;
03148 }else if (format == qh_PRINTinner) {
03149 innerouter= qh_RIDGEinner;
03150 printvridge= qh_printvnorm;
03151 }else if (format == qh_PRINTouter) {
03152 innerouter= qh_RIDGEouter;
03153 printvridge= qh_printvnorm;
03154 }else {
03155 qh_fprintf(qh ferr, 6219, "Qhull internal error (qh_printvdiagram): unknown print format %d.\n", format);
03156 qh_errexit(qh_ERRinput, NULL, NULL);
03157 }
03158 vertices= qh_markvoronoi(facetlist, facets, printall, &isLower, &numcenters);
03159 totcount= qh_printvdiagram2 (NULL, NULL, vertices, innerouter, False);
03160 qh_fprintf(fp, 9231, "%d\n", totcount);
03161 totcount= qh_printvdiagram2 (fp, printvridge, vertices, innerouter, True );
03162 qh_settempfree(&vertices);
03163 #if 0
03164 qh_fprintf(fp, 9232, "\n");
03165 totcount= qh_eachvoronoi_all(fp, printvridge, qh UPPERdelaunay, innerouter, True );
03166 qh_fprintf(fp, 9233, "%d\n", totcount);
03167 #endif
03168 }
03169
03170
03171
03172
03173
03174
03175
03176
03177
03178
03179
03180
03181
03182
03183
03184
03185
03186
03187
03188
03189
03190
03191
03192
03193
03194
03195
03196
03197
03198
03199 int qh_printvdiagram2 (FILE *fp, printvridgeT printvridge, setT *vertices, qh_RIDGE innerouter, boolT inorder) {
03200 int totcount= 0;
03201 int vertex_i, vertex_n;
03202 vertexT *vertex;
03203
03204 FORALLvertices
03205 vertex->seen= False;
03206 FOREACHvertex_i_(vertices) {
03207 if (vertex) {
03208 if (qh GOODvertex > 0 && qh_pointid(vertex->point)+1 != qh GOODvertex)
03209 continue;
03210 totcount += qh_eachvoronoi(fp, printvridge, vertex, !qh_ALL, innerouter, inorder);
03211 }
03212 }
03213 return totcount;
03214 }
03215
03216
03217
03218
03219
03220
03221
03222
03223 void qh_printvertex(FILE *fp, vertexT *vertex) {
03224 pointT *point;
03225 int k, count= 0;
03226 facetT *neighbor, **neighborp;
03227 realT r;
03228
03229 if (!vertex) {
03230 qh_fprintf(fp, 9234, " NULLvertex\n");
03231 return;
03232 }
03233 qh_fprintf(fp, 9235, "- p%d(v%d):", qh_pointid(vertex->point), vertex->id);
03234 point= vertex->point;
03235 if (point) {
03236 for (k=qh hull_dim; k--; ) {
03237 r= *point++;
03238 qh_fprintf(fp, 9236, " %5.2g", r);
03239 }
03240 }
03241 if (vertex->deleted)
03242 qh_fprintf(fp, 9237, " deleted");
03243 if (vertex->delridge)
03244 qh_fprintf(fp, 9238, " ridgedeleted");
03245 qh_fprintf(fp, 9239, "\n");
03246 if (vertex->neighbors) {
03247 qh_fprintf(fp, 9240, " neighbors:");
03248 FOREACHneighbor_(vertex) {
03249 if (++count % 100 == 0)
03250 qh_fprintf(fp, 9241, "\n ");
03251 qh_fprintf(fp, 9242, " f%d", neighbor->id);
03252 }
03253 qh_fprintf(fp, 9243, "\n");
03254 }
03255 }
03256
03257
03258
03259
03260
03261
03262
03263
03264
03265 void qh_printvertexlist(FILE *fp, const char* string, facetT *facetlist,
03266 setT *facets, boolT printall) {
03267 vertexT *vertex, **vertexp;
03268 setT *vertices;
03269
03270 vertices= qh_facetvertices(facetlist, facets, printall);
03271 qh_fprintf(fp, 9244, "%s", string);
03272 FOREACHvertex_(vertices)
03273 qh_printvertex(fp, vertex);
03274 qh_settempfree(&vertices);
03275 }
03276
03277
03278
03279
03280
03281
03282
03283
03284
03285 void qh_printvertices(FILE *fp, const char* string, setT *vertices) {
03286 vertexT *vertex, **vertexp;
03287
03288 qh_fprintf(fp, 9245, "%s", string);
03289 FOREACHvertex_(vertices)
03290 qh_fprintf(fp, 9246, " p%d(v%d)", qh_pointid(vertex->point), vertex->id);
03291 qh_fprintf(fp, 9247, "\n");
03292 }
03293
03294
03295
03296
03297
03298
03299
03300
03301
03302
03303
03304
03305
03306
03307
03308
03309
03310
03311 void qh_printvneighbors(FILE *fp, facetT* facetlist, setT *facets, boolT printall) {
03312 int numfacets, numsimplicial, numridges, totneighbors, numneighbors, numcoplanars, numtricoplanars;
03313 setT *vertices, *vertex_points, *coplanar_points;
03314 int numpoints= qh num_points + qh_setsize(qh other_points);
03315 vertexT *vertex, **vertexp;
03316 int vertex_i, vertex_n;
03317 facetT *facet, **facetp, *neighbor, **neighborp;
03318 pointT *point, **pointp;
03319
03320 qh_countfacets(facetlist, facets, printall, &numfacets, &numsimplicial,
03321 &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
03322 qh_fprintf(fp, 9248, "%d\n", numpoints);
03323 qh_vertexneighbors();
03324 vertices= qh_facetvertices(facetlist, facets, printall);
03325 vertex_points= qh_settemp(numpoints);
03326 coplanar_points= qh_settemp(numpoints);
03327 qh_setzero(vertex_points, 0, numpoints);
03328 qh_setzero(coplanar_points, 0, numpoints);
03329 FOREACHvertex_(vertices)
03330 qh_point_add(vertex_points, vertex->point, vertex);
03331 FORALLfacet_(facetlist) {
03332 FOREACHpoint_(facet->coplanarset)
03333 qh_point_add(coplanar_points, point, facet);
03334 }
03335 FOREACHfacet_(facets) {
03336 FOREACHpoint_(facet->coplanarset)
03337 qh_point_add(coplanar_points, point, facet);
03338 }
03339 FOREACHvertex_i_(vertex_points) {
03340 if (vertex) {
03341 numneighbors= qh_setsize(vertex->neighbors);
03342 qh_fprintf(fp, 9249, "%d", numneighbors);
03343 if (qh hull_dim == 3)
03344 qh_order_vertexneighbors(vertex);
03345 else if (qh hull_dim >= 4)
03346 qsort(SETaddr_(vertex->neighbors, facetT), (size_t)numneighbors,
03347 sizeof(facetT *), qh_compare_facetvisit);
03348 FOREACHneighbor_(vertex)
03349 qh_fprintf(fp, 9250, " %d",
03350 neighbor->visitid ? neighbor->visitid - 1 : 0 - neighbor->id);
03351 qh_fprintf(fp, 9251, "\n");
03352 }else if ((facet= SETelemt_(coplanar_points, vertex_i, facetT)))
03353 qh_fprintf(fp, 9252, "1 %d\n",
03354 facet->visitid ? facet->visitid - 1 : 0 - facet->id);
03355 else
03356 qh_fprintf(fp, 9253, "0\n");
03357 }
03358 qh_settempfree(&coplanar_points);
03359 qh_settempfree(&vertex_points);
03360 qh_settempfree(&vertices);
03361 }
03362
03363
03364
03365
03366
03367
03368
03369
03370
03371
03372
03373
03374
03375
03376
03377
03378
03379
03380
03381
03382
03383
03384
03385
03386 void qh_printvoronoi(FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall) {
03387 int k, numcenters, numvertices= 0, numneighbors, numinf, vid=1, vertex_i, vertex_n;
03388 facetT *facet, **facetp, *neighbor, **neighborp;
03389 setT *vertices;
03390 vertexT *vertex;
03391 boolT isLower;
03392 unsigned int numfacets= (unsigned int) qh num_facets;
03393
03394 vertices= qh_markvoronoi(facetlist, facets, printall, &isLower, &numcenters);
03395 FOREACHvertex_i_(vertices) {
03396 if (vertex) {
03397 numvertices++;
03398 numneighbors = numinf = 0;
03399 FOREACHneighbor_(vertex) {
03400 if (neighbor->visitid == 0)
03401 numinf= 1;
03402 else if (neighbor->visitid < numfacets)
03403 numneighbors++;
03404 }
03405 if (numinf && !numneighbors) {
03406 SETelem_(vertices, vertex_i)= NULL;
03407 numvertices--;
03408 }
03409 }
03410 }
03411 if (format == qh_PRINTgeom)
03412 qh_fprintf(fp, 9254, "{appearance {+edge -face} OFF %d %d 1 # Voronoi centers and cells\n",
03413 numcenters, numvertices);
03414 else
03415 qh_fprintf(fp, 9255, "%d\n%d %d 1\n", qh hull_dim-1, numcenters, qh_setsize(vertices));
03416 if (format == qh_PRINTgeom) {
03417 for (k=qh hull_dim-1; k--; )
03418 qh_fprintf(fp, 9256, qh_REAL_1, 0.0);
03419 qh_fprintf(fp, 9257, " 0 # infinity not used\n");
03420 }else {
03421 for (k=qh hull_dim-1; k--; )
03422 qh_fprintf(fp, 9258, qh_REAL_1, qh_INFINITE);
03423 qh_fprintf(fp, 9259, "\n");
03424 }
03425 FORALLfacet_(facetlist) {
03426 if (facet->visitid && facet->visitid < numfacets) {
03427 if (format == qh_PRINTgeom)
03428 qh_fprintf(fp, 9260, "# %d f%d\n", vid++, facet->id);
03429 qh_printcenter(fp, format, NULL, facet);
03430 }
03431 }
03432 FOREACHfacet_(facets) {
03433 if (facet->visitid && facet->visitid < numfacets) {
03434 if (format == qh_PRINTgeom)
03435 qh_fprintf(fp, 9261, "# %d f%d\n", vid++, facet->id);
03436 qh_printcenter(fp, format, NULL, facet);
03437 }
03438 }
03439 FOREACHvertex_i_(vertices) {
03440 numneighbors= 0;
03441 numinf=0;
03442 if (vertex) {
03443 if (qh hull_dim == 3)
03444 qh_order_vertexneighbors(vertex);
03445 else if (qh hull_dim >= 4)
03446 qsort(SETaddr_(vertex->neighbors, facetT),
03447 (size_t)qh_setsize(vertex->neighbors),
03448 sizeof(facetT *), qh_compare_facetvisit);
03449 FOREACHneighbor_(vertex) {
03450 if (neighbor->visitid == 0)
03451 numinf= 1;
03452 else if (neighbor->visitid < numfacets)
03453 numneighbors++;
03454 }
03455 }
03456 if (format == qh_PRINTgeom) {
03457 if (vertex) {
03458 qh_fprintf(fp, 9262, "%d", numneighbors);
03459 FOREACHneighbor_(vertex) {
03460 if (neighbor->visitid && neighbor->visitid < numfacets)
03461 qh_fprintf(fp, 9263, " %d", neighbor->visitid);
03462 }
03463 qh_fprintf(fp, 9264, " # p%d(v%d)\n", vertex_i, vertex->id);
03464 }else
03465 qh_fprintf(fp, 9265, " # p%d is coplanar or isolated\n", vertex_i);
03466 }else {
03467 if (numinf)
03468 numneighbors++;
03469 qh_fprintf(fp, 9266, "%d", numneighbors);
03470 if (vertex) {
03471 FOREACHneighbor_(vertex) {
03472 if (neighbor->visitid == 0) {
03473 if (numinf) {
03474 numinf= 0;
03475 qh_fprintf(fp, 9267, " %d", neighbor->visitid);
03476 }
03477 }else if (neighbor->visitid < numfacets)
03478 qh_fprintf(fp, 9268, " %d", neighbor->visitid);
03479 }
03480 }
03481 qh_fprintf(fp, 9269, "\n");
03482 }
03483 }
03484 if (format == qh_PRINTgeom)
03485 qh_fprintf(fp, 9270, "}\n");
03486 qh_settempfree(&vertices);
03487 }
03488
03489
03490
03491
03492
03493
03494
03495
03496
03497
03498
03499
03500
03501
03502
03503
03504
03505
03506 void qh_printvnorm(FILE *fp, vertexT *vertex, vertexT *vertexA, setT *centers, boolT unbounded) {
03507 pointT *normal;
03508 realT offset;
03509 int k;
03510 QHULL_UNUSED(unbounded);
03511
03512 normal= qh_detvnorm(vertex, vertexA, centers, &offset);
03513 qh_fprintf(fp, 9271, "%d %d %d ",
03514 2+qh hull_dim, qh_pointid(vertex->point), qh_pointid(vertexA->point));
03515 for (k=0; k< qh hull_dim-1; k++)
03516 qh_fprintf(fp, 9272, qh_REAL_1, normal[k]);
03517 qh_fprintf(fp, 9273, qh_REAL_1, offset);
03518 qh_fprintf(fp, 9274, "\n");
03519 }
03520
03521
03522
03523
03524
03525
03526
03527
03528
03529
03530
03531
03532
03533
03534
03535 void qh_printvridge(FILE *fp, vertexT *vertex, vertexT *vertexA, setT *centers, boolT unbounded) {
03536 facetT *facet, **facetp;
03537 QHULL_UNUSED(unbounded);
03538
03539 qh_fprintf(fp, 9275, "%d %d %d", qh_setsize(centers)+2,
03540 qh_pointid(vertex->point), qh_pointid(vertexA->point));
03541 FOREACHfacet_(centers)
03542 qh_fprintf(fp, 9276, " %d", facet->visitid);
03543 qh_fprintf(fp, 9277, "\n");
03544 }
03545
03546
03547
03548
03549
03550
03551
03552
03553
03554
03555
03556
03557 void qh_projectdim3 (pointT *source, pointT *destination) {
03558 int i,k;
03559
03560 for (k=0, i=0; k < qh hull_dim; k++) {
03561 if (qh hull_dim == 4) {
03562 if (k != qh DROPdim)
03563 destination[i++]= source[k];
03564 }else if (k == qh DROPdim)
03565 destination[i++]= 0;
03566 else
03567 destination[i++]= source[k];
03568 }
03569 while (i < 3)
03570 destination[i++]= 0.0;
03571 }
03572
03573
03574
03575
03576
03577
03578
03579
03580
03581
03582
03583
03584
03585
03586
03587
03588
03589
03590 int qh_readfeasible(int dim, const char *curline) {
03591 boolT isfirst= True;
03592 int linecount= 0, tokcount= 0;
03593 const char *s;
03594 char *t, firstline[qh_MAXfirst+1];
03595 coordT *coords, value;
03596
03597 if (!qh HALFspace) {
03598 qh_fprintf(qh ferr, 6070, "qhull input error: feasible point(dim 1 coords) is only valid for halfspace intersection\n");
03599 qh_errexit(qh_ERRinput, NULL, NULL);
03600 }
03601 if (qh feasible_string)
03602 qh_fprintf(qh ferr, 7057, "qhull input warning: feasible point(dim 1 coords) overrides 'Hn,n,n' feasible point for halfspace intersection\n");
03603 if (!(qh feasible_point= (coordT*)qh_malloc(dim* sizeof(coordT)))) {
03604 qh_fprintf(qh ferr, 6071, "qhull error: insufficient memory for feasible point\n");
03605 qh_errexit(qh_ERRmem, NULL, NULL);
03606 }
03607 coords= qh feasible_point;
03608 while ((s= (isfirst ? curline : fgets(firstline, qh_MAXfirst, qh fin)))) {
03609 if (isfirst)
03610 isfirst= False;
03611 else
03612 linecount++;
03613 while (*s) {
03614 while (isspace(*s))
03615 s++;
03616 value= qh_strtod(s, &t);
03617 if (s == t)
03618 break;
03619 s= t;
03620 *(coords++)= value;
03621 if (++tokcount == dim) {
03622 while (isspace(*s))
03623 s++;
03624 qh_strtod(s, &t);
03625 if (s != t) {
03626 qh_fprintf(qh ferr, 6072, "qhull input error: coordinates for feasible point do not finish out the line: %s\n",
03627 s);
03628 qh_errexit(qh_ERRinput, NULL, NULL);
03629 }
03630 return linecount;
03631 }
03632 }
03633 }
03634 qh_fprintf(qh ferr, 6073, "qhull input error: only %d coordinates. Could not read %d-d feasible point.\n",
03635 tokcount, dim);
03636 qh_errexit(qh_ERRinput, NULL, NULL);
03637 return 0;
03638 }
03639
03640
03641
03642
03643
03644
03645
03646
03647
03648
03649
03650
03651
03652
03653
03654
03655
03656
03657
03658
03659
03660
03661
03662
03663
03664
03665
03666
03667
03668
03669
03670
03671
03672 coordT *qh_readpoints(int *numpoints, int *dimension, boolT *ismalloc) {
03673 coordT *points, *coords, *infinity= NULL;
03674 realT paraboloid, maxboloid= -REALmax, value;
03675 realT *coordp= NULL, *offsetp= NULL, *normalp= NULL;
03676 char *s= 0, *t, firstline[qh_MAXfirst+1];
03677 int diminput=0, numinput=0, dimfeasible= 0, newnum, k, tempi;
03678 int firsttext=0, firstshort=0, firstlong=0, firstpoint=0;
03679 int tokcount= 0, linecount=0, maxcount, coordcount=0;
03680 boolT islong, isfirst= True, wasbegin= False;
03681 boolT isdelaunay= qh DELAUNAY && !qh PROJECTinput;
03682
03683 if (qh CDDinput) {
03684 while ((s= fgets(firstline, qh_MAXfirst, qh fin))) {
03685 linecount++;
03686 if (qh HALFspace && linecount == 1 && isdigit(*s)) {
03687 dimfeasible= qh_strtol(s, &s);
03688 while (isspace(*s))
03689 s++;
03690 if (qh_strtol(s, &s) == 1)
03691 linecount += qh_readfeasible(dimfeasible, s);
03692 else
03693 dimfeasible= 0;
03694 }else if (!memcmp(firstline, "begin", (size_t)5) || !memcmp(firstline, "BEGIN", (size_t)5))
03695 break;
03696 else if (!*qh rbox_command)
03697 strncat(qh rbox_command, s, sizeof(qh rbox_command)-1);
03698 }
03699 if (!s) {
03700 qh_fprintf(qh ferr, 6074, "qhull input error: missing \"begin\" for cdd-formated input\n");
03701 qh_errexit(qh_ERRinput, NULL, NULL);
03702 }
03703 }
03704 while (!numinput && (s= fgets(firstline, qh_MAXfirst, qh fin))) {
03705 linecount++;
03706 if (!memcmp(s, "begin", (size_t)5) || !memcmp(s, "BEGIN", (size_t)5))
03707 wasbegin= True;
03708 while (*s) {
03709 while (isspace(*s))
03710 s++;
03711 if (!*s)
03712 break;
03713 if (!isdigit(*s)) {
03714 if (!*qh rbox_command) {
03715 strncat(qh rbox_command, s, sizeof(qh rbox_command)-1);
03716 firsttext= linecount;
03717 }
03718 break;
03719 }
03720 if (!diminput)
03721 diminput= qh_strtol(s, &s);
03722 else {
03723 numinput= qh_strtol(s, &s);
03724 if (numinput == 1 && diminput >= 2 && qh HALFspace && !qh CDDinput) {
03725 linecount += qh_readfeasible(diminput, s);
03726 dimfeasible= diminput;
03727 diminput= numinput= 0;
03728 }else
03729 break;
03730 }
03731 }
03732 }
03733 if (!s) {
03734 qh_fprintf(qh ferr, 6075, "qhull input error: short input file. Did not find dimension and number of points\n");
03735 qh_errexit(qh_ERRinput, NULL, NULL);
03736 }
03737 if (diminput > numinput) {
03738 tempi= diminput;
03739 diminput= numinput;
03740 numinput= tempi;
03741 }
03742 if (diminput < 2) {
03743 qh_fprintf(qh ferr, 6220,"qhull input error: dimension %d(first number) should be at least 2\n",
03744 diminput);
03745 qh_errexit(qh_ERRinput, NULL, NULL);
03746 }
03747 if (isdelaunay) {
03748 qh PROJECTdelaunay= False;
03749 if (qh CDDinput)
03750 *dimension= diminput;
03751 else
03752 *dimension= diminput+1;
03753 *numpoints= numinput;
03754 if (qh ATinfinity)
03755 (*numpoints)++;
03756 }else if (qh HALFspace) {
03757 *dimension= diminput - 1;
03758 *numpoints= numinput;
03759 if (diminput < 3) {
03760 qh_fprintf(qh ferr, 6221,"qhull input error: dimension %d(first number, includes offset) should be at least 3 for halfspaces\n",
03761 diminput);
03762 qh_errexit(qh_ERRinput, NULL, NULL);
03763 }
03764 if (dimfeasible) {
03765 if (dimfeasible != *dimension) {
03766 qh_fprintf(qh ferr, 6222,"qhull input error: dimension %d of feasible point is not one less than dimension %d for halfspaces\n",
03767 dimfeasible, diminput);
03768 qh_errexit(qh_ERRinput, NULL, NULL);
03769 }
03770 }else
03771 qh_setfeasible(*dimension);
03772 }else {
03773 if (qh CDDinput)
03774 *dimension= diminput-1;
03775 else
03776 *dimension= diminput;
03777 *numpoints= numinput;
03778 }
03779 qh normal_size= *dimension * sizeof(coordT);
03780 if (qh HALFspace) {
03781 qh half_space= coordp= (coordT*)qh_malloc(qh normal_size + sizeof(coordT));
03782 if (qh CDDinput) {
03783 offsetp= qh half_space;
03784 normalp= offsetp + 1;
03785 }else {
03786 normalp= qh half_space;
03787 offsetp= normalp + *dimension;
03788 }
03789 }
03790 qh maxline= diminput * (qh_REALdigits + 5);
03791 maximize_(qh maxline, 500);
03792 qh line= (char*)qh_malloc((qh maxline+1) * sizeof(char));
03793 *ismalloc= True;
03794 coords= points= qh temp_malloc=
03795 (coordT*)qh_malloc((*numpoints)*(*dimension)*sizeof(coordT));
03796 if (!coords || !qh line || (qh HALFspace && !qh half_space)) {
03797 qh_fprintf(qh ferr, 6076, "qhull error: insufficient memory to read %d points\n",
03798 numinput);
03799 qh_errexit(qh_ERRmem, NULL, NULL);
03800 }
03801 if (isdelaunay && qh ATinfinity) {
03802 infinity= points + numinput * (*dimension);
03803 for (k= (*dimension) - 1; k--; )
03804 infinity[k]= 0.0;
03805 }
03806 maxcount= numinput * diminput;
03807 paraboloid= 0.0;
03808 while ((s= (isfirst ? s : fgets(qh line, qh maxline, qh fin)))) {
03809 if (!isfirst) {
03810 linecount++;
03811 if (*s == 'e' || *s == 'E') {
03812 if (!memcmp(s, "end", (size_t)3) || !memcmp(s, "END", (size_t)3)) {
03813 if (qh CDDinput )
03814 break;
03815 else if (wasbegin)
03816 qh_fprintf(qh ferr, 7058, "qhull input warning: the input appears to be in cdd format. If so, use 'Fd'\n");
03817 }
03818 }
03819 }
03820 islong= False;
03821 while (*s) {
03822 while (isspace(*s))
03823 s++;
03824 value= qh_strtod(s, &t);
03825 if (s == t) {
03826 if (!*qh rbox_command)
03827 strncat(qh rbox_command, s, sizeof(qh rbox_command)-1);
03828 if (*s && !firsttext)
03829 firsttext= linecount;
03830 if (!islong && !firstshort && coordcount)
03831 firstshort= linecount;
03832 break;
03833 }
03834 if (!firstpoint)
03835 firstpoint= linecount;
03836 s= t;
03837 if (++tokcount > maxcount)
03838 continue;
03839 if (qh HALFspace) {
03840 if (qh CDDinput)
03841 *(coordp++)= -value;
03842 else
03843 *(coordp++)= value;
03844 }else {
03845 *(coords++)= value;
03846 if (qh CDDinput && !coordcount) {
03847 if (value != 1.0) {
03848 qh_fprintf(qh ferr, 6077, "qhull input error: for cdd format, point at line %d does not start with '1'\n",
03849 linecount);
03850 qh_errexit(qh_ERRinput, NULL, NULL);
03851 }
03852 coords--;
03853 }else if (isdelaunay) {
03854 paraboloid += value * value;
03855 if (qh ATinfinity) {
03856 if (qh CDDinput)
03857 infinity[coordcount-1] += value;
03858 else
03859 infinity[coordcount] += value;
03860 }
03861 }
03862 }
03863 if (++coordcount == diminput) {
03864 coordcount= 0;
03865 if (isdelaunay) {
03866 *(coords++)= paraboloid;
03867 maximize_(maxboloid, paraboloid);
03868 paraboloid= 0.0;
03869 }else if (qh HALFspace) {
03870 if (!qh_sethalfspace(*dimension, coords, &coords, normalp, offsetp, qh feasible_point)) {
03871 qh_fprintf(qh ferr, 8048, "The halfspace was on line %d\n", linecount);
03872 if (wasbegin)
03873 qh_fprintf(qh ferr, 8049, "The input appears to be in cdd format. If so, you should use option 'Fd'\n");
03874 qh_errexit(qh_ERRinput, NULL, NULL);
03875 }
03876 coordp= qh half_space;
03877 }
03878 while (isspace(*s))
03879 s++;
03880 if (*s) {
03881 islong= True;
03882 if (!firstlong)
03883 firstlong= linecount;
03884 }
03885 }
03886 }
03887 if (!islong && !firstshort && coordcount)
03888 firstshort= linecount;
03889 if (!isfirst && s - qh line >= qh maxline) {
03890 qh_fprintf(qh ferr, 6078, "qhull input error: line %d contained more than %d characters\n",
03891 linecount, (int) (s - qh line));
03892 qh_errexit(qh_ERRinput, NULL, NULL);
03893 }
03894 isfirst= False;
03895 }
03896 if (tokcount != maxcount) {
03897 newnum= fmin_(numinput, tokcount/diminput);
03898 qh_fprintf(qh ferr, 7073,"\
03899 qhull warning: instead of %d %d-dimensional points, input contains\n\
03900 %d points and %d extra coordinates. Line %d is the first\npoint",
03901 numinput, diminput, tokcount/diminput, tokcount % diminput, firstpoint);
03902 if (firsttext)
03903 qh_fprintf(qh ferr, 8051, ", line %d is the first comment", firsttext);
03904 if (firstshort)
03905 qh_fprintf(qh ferr, 8052, ", line %d is the first short\nline", firstshort);
03906 if (firstlong)
03907 qh_fprintf(qh ferr, 8053, ", line %d is the first long line", firstlong);
03908 qh_fprintf(qh ferr, 8054, ". Continue with %d points.\n", newnum);
03909 numinput= newnum;
03910 if (isdelaunay && qh ATinfinity) {
03911 for (k= tokcount % diminput; k--; )
03912 infinity[k] -= *(--coords);
03913 *numpoints= newnum+1;
03914 }else {
03915 coords -= tokcount % diminput;
03916 *numpoints= newnum;
03917 }
03918 }
03919 if (isdelaunay && qh ATinfinity) {
03920 for (k= (*dimension) -1; k--; )
03921 infinity[k] /= numinput;
03922 if (coords == infinity)
03923 coords += (*dimension) -1;
03924 else {
03925 for (k=0; k < (*dimension) -1; k++)
03926 *(coords++)= infinity[k];
03927 }
03928 *(coords++)= maxboloid * 1.1;
03929 }
03930 if (qh rbox_command[0]) {
03931 qh rbox_command[strlen(qh rbox_command)-1]= '\0';
03932 if (!strcmp(qh rbox_command, "./rbox D4"))
03933 qh_fprintf(qh ferr, 8055, "\n\
03934 This is the qhull test case. If any errors or core dumps occur,\n\
03935 recompile qhull with 'make new'. If errors still occur, there is\n\
03936 an incompatibility. You should try a different compiler. You can also\n\
03937 change the choices in user.h. If you discover the source of the problem,\n\
03938 please send mail to qhull_bug@qhull.org.\n\
03939 \n\
03940 Type 'qhull' for a short list of options.\n");
03941 }
03942 qh_free(qh line);
03943 qh line= NULL;
03944 if (qh half_space) {
03945 qh_free(qh half_space);
03946 qh half_space= NULL;
03947 }
03948 qh temp_malloc= NULL;
03949 trace1((qh ferr, 1008,"qh_readpoints: read in %d %d-dimensional points\n",
03950 numinput, diminput));
03951 return(points);
03952 }
03953
03954
03955
03956
03957
03958
03959
03960
03961
03962
03963
03964
03965 void qh_setfeasible(int dim) {
03966 int tokcount= 0;
03967 char *s;
03968 coordT *coords, value;
03969
03970 if (!(s= qh feasible_string)) {
03971 qh_fprintf(qh ferr, 6223, "\
03972 qhull input error: halfspace intersection needs a feasible point.\n\
03973 Either prepend the input with 1 point or use 'Hn,n,n'. See manual.\n");
03974 qh_errexit(qh_ERRinput, NULL, NULL);
03975 }
03976 if (!(qh feasible_point= (pointT*)qh_malloc(dim * sizeof(coordT)))) {
03977 qh_fprintf(qh ferr, 6079, "qhull error: insufficient memory for 'Hn,n,n'\n");
03978 qh_errexit(qh_ERRmem, NULL, NULL);
03979 }
03980 coords= qh feasible_point;
03981 while (*s) {
03982 value= qh_strtod(s, &s);
03983 if (++tokcount > dim) {
03984 qh_fprintf(qh ferr, 7059, "qhull input warning: more coordinates for 'H%s' than dimension %d\n",
03985 qh feasible_string, dim);
03986 break;
03987 }
03988 *(coords++)= value;
03989 if (*s)
03990 s++;
03991 }
03992 while (++tokcount <= dim)
03993 *(coords++)= 0.0;
03994 }
03995
03996
03997
03998
03999
04000
04001
04002
04003
04004
04005 boolT qh_skipfacet(facetT *facet) {
04006 facetT *neighbor, **neighborp;
04007
04008 if (qh PRINTneighbors) {
04009 if (facet->good)
04010 return !qh PRINTgood;
04011 FOREACHneighbor_(facet) {
04012 if (neighbor->good)
04013 return False;
04014 }
04015 return True;
04016 }else if (qh PRINTgood)
04017 return !facet->good;
04018 else if (!facet->normal)
04019 return True;
04020 return(!qh_inthresholds(facet->normal, NULL));
04021 }
04022
04023
04024
04025
04026
04027
04028
04029
04030
04031
04032
04033
04034
04035 char *qh_skipfilename(char *filename) {
04036 char *s= filename;
04037 char c;
04038
04039 while (*s && isspace(*s))
04040 s++;
04041 c= *s++;
04042 if (c == '\0') {
04043 qh_fprintf(qh ferr, 6204, "qhull input error: filename expected, none found.\n");
04044 qh_errexit(qh_ERRinput, NULL, NULL);
04045 }
04046 if (c == '\'' || c == '"') {
04047 while (*s !=c || s[-1] == '\\') {
04048 if (!*s) {
04049 qh_fprintf(qh ferr, 6203, "qhull input error: missing quote after filename -- %s\n", filename);
04050 qh_errexit(qh_ERRinput, NULL, NULL);
04051 }
04052 s++;
04053 }
04054 s++;
04055 }
04056 else while (*s && !isspace(*s))
04057 s++;
04058 return s;
04059 }
04060
