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