libqhull-ours: io.c Source File

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