graspit: /opt/ros/diamondback/stacks/graspit_simulator/graspit/graspit

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