new round of warning fixes. we are now down to 24 with Xcode on blender
[blender.git] / extern / qhull / src / io.c
1 /*<html><pre>  -<a                             href="qh-io.htm"
2   >-------------------------------</a><a name="TOP">-</a>
3
4    io.c 
5    Input/Output routines of qhull application
6
7    see qh-io.htm and io.h
8
9    see user.c for qh_errprint and qh_printfacetlist
10
11    unix.c calls qh_readpoints and qh_produce_output
12
13    unix.c and user.c are the only callers of io.c functions
14    This allows the user to avoid loading io.o from qhull.a
15
16    copyright (c) 1993-2002 The Geometry Center        
17 */
18
19 #include "qhull_a.h"
20
21 /*========= -prototypes for internal functions ========= */
22
23 static int qh_compare_facetarea(const void *p1, const void *p2);
24 static int qh_compare_facetmerge(const void *p1, const void *p2);
25 static int qh_compare_facetvisit(const void *p1, const void *p2);
26 int qh_compare_vertexpoint(const void *p1, const void *p2); /* not used */
27
28 /*========= -functions in alphabetical order after qh_produce_output()  =====*/
29
30 /*-<a                             href="qh-io.htm#TOC"
31   >-------------------------------</a><a name="produce_output">-</a>
32   
33   qh_produce_output()
34     prints out the result of qhull in desired format
35     if qh.GETarea
36       computes and prints area and volume
37     qh.PRINTout[] is an array of output formats
38
39   notes:
40     prints output in qh.PRINTout order
41 */
42 void qh_produce_output(void) {
43   int i, tempsize= qh_setsize ((setT*)qhmem.tempstack), d_1;
44
45   if (qh VORONOI) {
46     qh_clearcenters (qh_ASvoronoi);
47     qh_vertexneighbors();
48   }
49   if (qh TRIangulate) {
50     qh_triangulate(); 
51     if (qh VERIFYoutput && !qh CHECKfrequently) 
52       qh_checkpolygon (qh facet_list);
53   }
54   qh_findgood_all (qh facet_list); 
55   if (qh GETarea)
56     qh_getarea(qh facet_list);
57   if (qh KEEParea || qh KEEPmerge || qh KEEPminArea < REALmax/2)
58     qh_markkeep (qh facet_list);
59   if (qh PRINTsummary)
60     qh_printsummary(qh ferr);
61   else if (qh PRINTout[0] == qh_PRINTnone)
62     qh_printsummary(qh fout);
63   for (i= 0; i < qh_PRINTEND; i++)
64     qh_printfacets (qh fout, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);
65   qh_allstatistics();
66   if (qh PRINTprecision && !qh MERGING && (qh JOGGLEmax > REALmax/2 || qh RERUN))
67     qh_printstats (qh ferr, qhstat precision, NULL);
68   if (qh VERIFYoutput && (zzval_(Zridge) > 0 || zzval_(Zridgemid) > 0)) 
69     qh_printstats (qh ferr, qhstat vridges, NULL);
70   if (qh PRINTstatistics) {
71     qh_collectstatistics();
72     qh_printstatistics(qh ferr, "");
73     qh_memstatistics (qh ferr);
74     d_1= sizeof(setT) + (qh hull_dim - 1) * SETelemsize;
75     fprintf(qh ferr, "\
76     size in bytes: merge %ld ridge %ld vertex %ld facet %ld\n\
77          normal %d ridge vertices %d facet vertices or neighbors %ld\n",
78             sizeof(mergeT), sizeof(ridgeT),
79             sizeof(vertexT), sizeof(facetT),
80             qh normal_size, d_1, d_1 + SETelemsize);
81   }
82   if (qh_setsize ((setT*)qhmem.tempstack) != tempsize) {
83     fprintf (qh ferr, "qhull internal error (qh_produce_output): temporary sets not empty (%d)\n",
84              qh_setsize ((setT*)qhmem.tempstack));
85     qh_errexit (qh_ERRqhull, NULL, NULL);
86   }
87 } /* produce_output */
88
89
90 /*-<a                             href="qh-io.htm#TOC"
91   >-------------------------------</a><a name="dfacet">-</a>
92   
93   dfacet( id )
94     print facet by id, for debugging
95
96 */
97 void dfacet (unsigned id) {
98   facetT *facet;
99
100   FORALLfacets {
101     if (facet->id == id) {
102       qh_printfacet (qh fout, facet);
103       break;
104     }
105   }
106 } /* dfacet */
107
108
109 /*-<a                             href="qh-io.htm#TOC"
110   >-------------------------------</a><a name="dvertex">-</a>
111   
112   dvertex( id )
113     print vertex by id, for debugging
114 */
115 void dvertex (unsigned id) {
116   vertexT *vertex;
117
118   FORALLvertices {
119     if (vertex->id == id) {
120       qh_printvertex (qh fout, vertex);
121       break;
122     }
123   }
124 } /* dvertex */
125
126
127 /*-<a                             href="qh-io.htm#TOC"
128   >-------------------------------</a><a name="compare_vertexpoint">-</a>
129   
130   qh_compare_vertexpoint( p1, p2 )
131     used by qsort() to order vertices by point id 
132 */
133 int qh_compare_vertexpoint(const void *p1, const void *p2) {
134   vertexT *a= *((vertexT **)p1), *b= *((vertexT **)p2);
135  
136   return ((qh_pointid(a->point) > qh_pointid(b->point)?1:-1));
137 } /* compare_vertexpoint */
138
139 /*-<a                             href="qh-io.htm#TOC"
140   >-------------------------------</a><a name="compare_facetarea">-</a>
141   
142   qh_compare_facetarea( p1, p2 )
143     used by qsort() to order facets by area
144 */
145 static int qh_compare_facetarea(const void *p1, const void *p2) {
146   facetT *a= *((facetT **)p1), *b= *((facetT **)p2);
147
148   if (!a->isarea)
149     return -1;
150   if (!b->isarea)
151     return 1; 
152   if (a->f.area > b->f.area)
153     return 1;
154   else if (a->f.area == b->f.area)
155     return 0;
156   return -1;
157 } /* compare_facetarea */
158
159 /*-<a                             href="qh-io.htm#TOC"
160   >-------------------------------</a><a name="compare_facetmerge">-</a>
161   
162   qh_compare_facetmerge( p1, p2 )
163     used by qsort() to order facets by number of merges
164 */
165 static int qh_compare_facetmerge(const void *p1, const void *p2) {
166   facetT *a= *((facetT **)p1), *b= *((facetT **)p2);
167  
168   return (a->nummerge - b->nummerge);
169 } /* compare_facetvisit */
170
171 /*-<a                             href="qh-io.htm#TOC"
172   >-------------------------------</a><a name="compare_facetvisit">-</a>
173   
174   qh_compare_facetvisit( p1, p2 )
175     used by qsort() to order facets by visit id or id
176 */
177 static int qh_compare_facetvisit(const void *p1, const void *p2) {
178   facetT *a= *((facetT **)p1), *b= *((facetT **)p2);
179   int i,j;
180
181   if (!(i= a->visitid))
182     i= - a->id; /* do not convert to int */
183   if (!(j= b->visitid))
184     j= - b->id;
185   return (i - j);
186 } /* compare_facetvisit */
187
188 /*-<a                             href="qh-io.htm#TOC"
189   >-------------------------------</a><a name="countfacets">-</a>
190   
191   qh_countfacets( facetlist, facets, printall, 
192           numfacets, numsimplicial, totneighbors, numridges, numcoplanar, numtricoplanars  )
193     count good facets for printing and set visitid
194     if allfacets, ignores qh_skipfacet()
195
196   notes:
197     qh_printsummary and qh_countfacets must match counts
198
199   returns:
200     numfacets, numsimplicial, total neighbors, numridges, coplanars
201     each facet with ->visitid indicating 1-relative position
202       ->visitid==0 indicates not good
203   
204   notes
205     numfacets >= numsimplicial
206     if qh.NEWfacets, 
207       does not count visible facets (matches qh_printafacet)
208
209   design:
210     for all facets on facetlist and in facets set
211       unless facet is skipped or visible (i.e., will be deleted)
212         mark facet->visitid
213         update counts
214 */
215 void qh_countfacets (facetT *facetlist, setT *facets, boolT printall,
216     int *numfacetsp, int *numsimplicialp, int *totneighborsp, int *numridgesp, int *numcoplanarsp, int *numtricoplanarsp) {
217   facetT *facet, **facetp;
218   int numfacets= 0, numsimplicial= 0, numridges= 0, totneighbors= 0, numcoplanars= 0, numtricoplanars= 0;
219
220   FORALLfacet_(facetlist) {
221     if ((facet->visible && qh NEWfacets)
222     || (!printall && qh_skipfacet(facet)))
223       facet->visitid= 0;
224     else {
225       facet->visitid= ++numfacets;
226       totneighbors += qh_setsize (facet->neighbors);
227       if (facet->simplicial) {
228         numsimplicial++;
229         if (facet->keepcentrum && facet->tricoplanar)
230           numtricoplanars++;
231       }else
232         numridges += qh_setsize (facet->ridges);
233       if (facet->coplanarset)
234         numcoplanars += qh_setsize (facet->coplanarset);
235     }
236   }
237   FOREACHfacet_(facets) {
238     if ((facet->visible && qh NEWfacets)
239     || (!printall && qh_skipfacet(facet)))
240       facet->visitid= 0;
241     else {
242       facet->visitid= ++numfacets;
243       totneighbors += qh_setsize (facet->neighbors);
244       if (facet->simplicial){
245         numsimplicial++;
246         if (facet->keepcentrum && facet->tricoplanar)
247           numtricoplanars++;
248       }else
249         numridges += qh_setsize (facet->ridges);
250       if (facet->coplanarset)
251         numcoplanars += qh_setsize (facet->coplanarset);
252     }
253   }
254   qh visit_id += numfacets+1;
255   *numfacetsp= numfacets;
256   *numsimplicialp= numsimplicial;
257   *totneighborsp= totneighbors;
258   *numridgesp= numridges;
259   *numcoplanarsp= numcoplanars;
260   *numtricoplanarsp= numtricoplanars;
261 } /* countfacets */
262
263 /*-<a                             href="qh-io.htm#TOC"
264   >-------------------------------</a><a name="detvnorm">-</a>
265   
266   qh_detvnorm( vertex, vertexA, centers, offset )
267     compute separating plane of the Voronoi diagram for a pair of input sites
268     centers= set of facets (i.e., Voronoi vertices)
269       facet->visitid= 0 iff vertex-at-infinity (i.e., unbounded)
270         
271   assumes:
272     qh_ASvoronoi and qh_vertexneighbors() already set
273   
274   returns:
275     norm
276       a pointer into qh.gm_matrix to qh.hull_dim-1 reals
277       copy the data before reusing qh.gm_matrix
278     offset
279       if 'QVn'
280         sign adjusted so that qh.GOODvertexp is inside
281       else
282         sign adjusted so that vertex is inside
283       
284     qh.gm_matrix= simplex of points from centers relative to first center
285     
286   notes:
287     in io.c so that code for 'v Tv' can be removed by removing io.c
288     returns pointer into qh.gm_matrix to avoid tracking of temporary memory
289   
290   design:
291     determine midpoint of input sites
292     build points as the set of Voronoi vertices
293     select a simplex from points (if necessary)
294       include midpoint if the Voronoi region is unbounded
295     relocate the first vertex of the simplex to the origin
296     compute the normalized hyperplane through the simplex
297     orient the hyperplane toward 'QVn' or 'vertex'
298     if 'Tv' or 'Ts'
299       if bounded
300         test that hyperplane is the perpendicular bisector of the input sites
301       test that Voronoi vertices not in the simplex are still on the hyperplane
302     free up temporary memory
303 */
304 pointT *qh_detvnorm (vertexT *vertex, vertexT *vertexA, setT *centers, realT *offsetp) {
305   facetT *facet, **facetp;
306   int  i, k, pointid, pointidA, point_i, point_n;
307   setT *simplex= NULL;
308   pointT *point, **pointp, *point0, *midpoint, *normal, *inpoint;
309   coordT *coord, *gmcoord, *normalp;
310   setT *points= qh_settemp (qh TEMPsize);
311   boolT nearzero= False;
312   boolT unbounded= False;
313   int numcenters= 0;
314   int dim= qh hull_dim - 1;
315   realT dist, offset, angle, zero= 0.0;
316
317   midpoint= qh gm_matrix + qh hull_dim * qh hull_dim;  /* last row */
318   for (k= 0; k < dim; k++)
319     midpoint[k]= (vertex->point[k] + vertexA->point[k])/2;
320   FOREACHfacet_(centers) {
321     numcenters++;
322     if (!facet->visitid)
323       unbounded= True;
324     else {
325       if (!facet->center)
326         facet->center= qh_facetcenter (facet->vertices);
327       qh_setappend (&points, facet->center);
328     }
329   }
330   if (numcenters > dim) {
331     simplex= qh_settemp (qh TEMPsize);
332     qh_setappend (&simplex, vertex->point);
333     if (unbounded)
334       qh_setappend (&simplex, midpoint);
335     qh_maxsimplex (dim, points, NULL, 0, &simplex);
336     qh_setdelnth (simplex, 0);
337   }else if (numcenters == dim) {
338     if (unbounded)
339       qh_setappend (&points, midpoint);
340     simplex= points; 
341   }else {
342     fprintf(qh ferr, "qh_detvnorm: too few points (%d) to compute separating plane\n", numcenters);
343     qh_errexit (qh_ERRqhull, NULL, NULL);
344   }
345   i= 0;
346   gmcoord= qh gm_matrix;
347   point0= SETfirstt_(simplex, pointT);
348   FOREACHpoint_(simplex) {
349     if (qh IStracing >= 4)
350       qh_printmatrix(qh ferr, "qh_detvnorm: Voronoi vertex or midpoint", 
351                               &point, 1, dim);
352     if (point != point0) {
353       qh gm_row[i++]= gmcoord;
354       coord= point0;
355       for (k= dim; k--; )
356         *(gmcoord++)= *point++ - *coord++;
357     }
358   }
359   qh gm_row[i]= gmcoord;  /* does not overlap midpoint, may be used later for qh_areasimplex */
360   normal= gmcoord;
361   qh_sethyperplane_gauss (dim, qh gm_row, point0, True,
362                 normal, &offset, &nearzero);
363   if (qh GOODvertexp == vertexA->point)
364     inpoint= vertexA->point;
365   else
366     inpoint= vertex->point;
367   zinc_(Zdistio);
368   dist= qh_distnorm (dim, inpoint, normal, &offset);
369   if (dist > 0) {
370     offset= -offset;
371     normalp= normal;
372     for (k= dim; k--; ) {
373       *normalp= -(*normalp);
374       normalp++;
375     }
376   }
377   if (qh VERIFYoutput || qh PRINTstatistics) {
378     pointid= qh_pointid (vertex->point);
379     pointidA= qh_pointid (vertexA->point);
380     if (!unbounded) {
381       zinc_(Zdiststat);
382       dist= qh_distnorm (dim, midpoint, normal, &offset);
383       if (dist < 0)
384         dist= -dist;
385       zzinc_(Zridgemid);
386       wwmax_(Wridgemidmax, dist);
387       wwadd_(Wridgemid, dist);
388       trace4((qh ferr, "qh_detvnorm: points %d %d midpoint dist %2.2g\n",
389                  pointid, pointidA, dist));
390       for (k= 0; k < dim; k++) 
391         midpoint[k]= vertexA->point[k] - vertex->point[k];  /* overwrites midpoint! */
392       qh_normalize (midpoint, dim, False);
393       angle= qh_distnorm (dim, midpoint, normal, &zero); /* qh_detangle uses dim+1 */
394       if (angle < 0.0)
395         angle= angle + 1.0;
396       else
397         angle= angle - 1.0;
398       if (angle < 0.0)
399         angle -= angle;
400       trace4((qh ferr, "qh_detvnorm: points %d %d angle %2.2g nearzero %d\n",
401                  pointid, pointidA, angle, nearzero));
402       if (nearzero) {
403         zzinc_(Zridge0);
404         wwmax_(Wridge0max, angle);
405         wwadd_(Wridge0, angle);
406       }else {
407         zzinc_(Zridgeok)
408         wwmax_(Wridgeokmax, angle);
409         wwadd_(Wridgeok, angle);
410       }
411     }
412     if (simplex != points) {
413       FOREACHpoint_i_(points) {
414         if (!qh_setin (simplex, point)) {
415           facet= SETelemt_(centers, point_i, facetT);
416           zinc_(Zdiststat);
417           dist= qh_distnorm (dim, point, normal, &offset);
418           if (dist < 0)
419             dist= -dist;
420           zzinc_(Zridge);
421           wwmax_(Wridgemax, dist);
422           wwadd_(Wridge, dist);
423           trace4((qh ferr, "qh_detvnorm: points %d %d Voronoi vertex %d dist %2.2g\n",
424                              pointid, pointidA, facet->visitid, dist));
425         }
426       }
427     }
428   }
429   *offsetp= offset;
430   if (simplex != points)
431     qh_settempfree (&simplex);
432   qh_settempfree (&points);
433   return normal;
434 } /* detvnorm */
435
436 /*-<a                             href="qh-io.htm#TOC"
437   >-------------------------------</a><a name="detvridge">-</a>
438
439   qh_detvridge( vertexA )
440     determine Voronoi ridge from 'seen' neighbors of vertexA
441     include one vertex-at-infinite if an !neighbor->visitid
442
443   returns:
444     temporary set of centers (facets, i.e., Voronoi vertices)
445     sorted by center id
446 */
447 setT *qh_detvridge (vertexT *vertex) {
448   setT *centers= qh_settemp (qh TEMPsize);
449   setT *tricenters= qh_settemp (qh TEMPsize);
450   facetT *neighbor, **neighborp;
451   boolT firstinf= True;
452   
453   FOREACHneighbor_(vertex) {
454     if (neighbor->seen) {
455       if (neighbor->visitid) {
456         if (!neighbor->tricoplanar || qh_setunique (&tricenters, neighbor->center)) 
457           qh_setappend (&centers, neighbor);
458       }else if (firstinf) {
459         firstinf= False;
460         qh_setappend (&centers, neighbor);
461       }
462     }
463   }
464   qsort (SETaddr_(centers, facetT), qh_setsize (centers),
465              sizeof (facetT *), qh_compare_facetvisit);
466   qh_settempfree (&tricenters);
467   return centers;
468 } /* detvridge */      
469
470 /*-<a                             href="qh-io.htm#TOC"
471   >-------------------------------</a><a name="detvridge3">-</a>
472
473   qh_detvridge3( atvertex, vertex )
474     determine 3-d Voronoi ridge from 'seen' neighbors of atvertex and vertex
475     include one vertex-at-infinite for !neighbor->visitid
476     assumes all facet->seen2= True
477
478   returns:
479     temporary set of centers (facets, i.e., Voronoi vertices)
480     listed in adjacency order (not oriented)
481     all facet->seen2= True
482
483   design:
484     mark all neighbors of atvertex
485     for each adjacent neighbor of both atvertex and vertex
486       if neighbor selected
487         add neighbor to set of Voronoi vertices
488 */
489 setT *qh_detvridge3 (vertexT *atvertex, vertexT *vertex) {
490   setT *centers= qh_settemp (qh TEMPsize);
491   setT *tricenters= qh_settemp (qh TEMPsize);
492   facetT *neighbor, **neighborp, *facet= NULL;
493   boolT firstinf= True;
494   
495   FOREACHneighbor_(atvertex)
496     neighbor->seen2= False;
497   FOREACHneighbor_(vertex) {
498     if (!neighbor->seen2) {
499       facet= neighbor;
500       break;
501     }
502   }
503   while (facet) { 
504     facet->seen2= True;
505     if (neighbor->seen) {
506       if (facet->visitid) {
507         if (!facet->tricoplanar || qh_setunique (&tricenters, facet->center)) 
508           qh_setappend (&centers, facet);
509       }else if (firstinf) {
510         firstinf= False;
511         qh_setappend (&centers, facet);
512       }
513     }
514     FOREACHneighbor_(facet) {
515       if (!neighbor->seen2) {
516         if (qh_setin (vertex->neighbors, neighbor))
517           break;
518         else
519           neighbor->seen2= True;
520       }
521     }
522     facet= neighbor;
523   }
524   if (qh CHECKfrequently) {
525     FOREACHneighbor_(vertex) {
526       if (!neighbor->seen2) {
527         fprintf (stderr, "qh_detvridge3: neigbors of vertex p%d are not connected at facet %d\n",
528                  qh_pointid (vertex->point), neighbor->id);
529         qh_errexit (qh_ERRqhull, neighbor, NULL);
530       }
531     }
532   }
533   FOREACHneighbor_(atvertex) 
534     neighbor->seen2= True;
535   qh_settempfree (&tricenters);
536   return centers;
537 } /* detvridge3 */      
538
539 /*-<a                             href="qh-io.htm#TOC"
540   >-------------------------------</a><a name="eachvoronoi">-</a>
541   
542   qh_eachvoronoi( fp, printvridge, vertex, visitall, innerouter, inorder )
543     if visitall,
544       visit all Voronoi ridges for vertex (i.e., an input site)
545     else
546       visit all unvisited Voronoi ridges for vertex
547       all vertex->seen= False if unvisited
548     assumes
549       all facet->seen= False
550       all facet->seen2= True (for qh_detvridge3)
551       all facet->visitid == 0 if vertex_at_infinity
552                          == index of Voronoi vertex 
553                          >= qh.num_facets if ignored
554     innerouter:
555       qh_RIDGEall--  both inner (bounded) and outer (unbounded) ridges
556       qh_RIDGEinner- only inner
557       qh_RIDGEouter- only outer
558       
559     if inorder
560       orders vertices for 3-d Voronoi diagrams
561   
562   returns:
563     number of visited ridges (does not include previously visited ridges)
564     
565     if printvridge,
566       calls printvridge( fp, vertex, vertexA, centers)
567         fp== any pointer (assumes FILE*)
568         vertex,vertexA= pair of input sites that define a Voronoi ridge
569         centers= set of facets (i.e., Voronoi vertices)
570                  ->visitid == index or 0 if vertex_at_infinity
571                  ordered for 3-d Voronoi diagram
572   notes:
573     uses qh.vertex_visit
574   
575   see:
576     qh_eachvoronoi_all()
577   
578   design:
579     mark selected neighbors of atvertex
580     for each selected neighbor (either Voronoi vertex or vertex-at-infinity)
581       for each unvisited vertex 
582         if atvertex and vertex share more than d-1 neighbors
583           bump totalcount
584           if printvridge defined
585             build the set of shared neighbors (i.e., Voronoi vertices)
586             call printvridge
587 */
588 int qh_eachvoronoi (FILE *fp, printvridgeT printvridge, vertexT *atvertex, boolT visitall, qh_RIDGE innerouter, boolT inorder) {
589   boolT unbounded;
590   int count;
591   facetT *neighbor, **neighborp, *neighborA, **neighborAp;
592   setT *centers;
593   setT *tricenters= qh_settemp (qh TEMPsize);
594
595   vertexT *vertex, **vertexp;
596   boolT firstinf;
597   unsigned int numfacets= (unsigned int)qh num_facets;
598   int totridges= 0;
599
600   qh vertex_visit++;
601   atvertex->seen= True;
602   if (visitall) {
603     FORALLvertices 
604       vertex->seen= False;
605   }
606   FOREACHneighbor_(atvertex) {
607     if (neighbor->visitid < numfacets) 
608       neighbor->seen= True;
609   }
610   FOREACHneighbor_(atvertex) {
611     if (neighbor->seen) {
612       FOREACHvertex_(neighbor->vertices) {
613         if (vertex->visitid != qh vertex_visit && !vertex->seen) {
614           vertex->visitid= qh vertex_visit;
615           count= 0;
616           firstinf= True;
617           qh_settruncate (tricenters, 0);
618           FOREACHneighborA_(vertex) {
619             if (neighborA->seen) {
620               if (neighborA->visitid) {
621                 if (!neighborA->tricoplanar || qh_setunique (&tricenters, neighborA->center))
622                   count++;
623               }else if (firstinf) {
624                 count++;
625                 firstinf= False;
626               }
627             }
628           }
629           if (count >= qh hull_dim - 1) {  /* e.g., 3 for 3-d Voronoi */
630             if (firstinf) {
631               if (innerouter == qh_RIDGEouter)
632                 continue;
633               unbounded= False;
634             }else {
635               if (innerouter == qh_RIDGEinner)
636                 continue;
637               unbounded= True;
638             }
639             totridges++;
640             trace4((qh ferr, "qh_eachvoronoi: Voronoi ridge of %d vertices between sites %d and %d\n",
641                   count, qh_pointid (atvertex->point), qh_pointid (vertex->point)));
642             if (printvridge) { 
643               if (inorder && qh hull_dim == 3+1) /* 3-d Voronoi diagram */
644                 centers= qh_detvridge3 (atvertex, vertex);
645               else
646                 centers= qh_detvridge (vertex);
647               (*printvridge) (fp, atvertex, vertex, centers, unbounded);
648               qh_settempfree (&centers);
649             }
650           }
651         }
652       }
653     }
654   }
655   FOREACHneighbor_(atvertex) 
656     neighbor->seen= False;
657   qh_settempfree (&tricenters);
658   return totridges;
659 } /* eachvoronoi */
660   
661
662 /*-<a                             href="qh-poly.htm#TOC"
663   >-------------------------------</a><a name="eachvoronoi_all">-</a>
664   
665   qh_eachvoronoi_all( fp, printvridge, isupper, innerouter, inorder )
666     visit all Voronoi ridges
667     
668     innerouter:
669       see qh_eachvoronoi()
670       
671     if inorder
672       orders vertices for 3-d Voronoi diagrams
673     
674   returns
675     total number of ridges 
676
677     if isupper == facet->upperdelaunay  (i.e., a Vornoi vertex)
678       facet->visitid= Voronoi vertex index (same as 'o' format)
679     else 
680       facet->visitid= 0
681
682     if printvridge,
683       calls printvridge( fp, vertex, vertexA, centers)
684       [see qh_eachvoronoi]
685       
686   notes:
687     Not used for qhull.exe
688     same effect as qh_printvdiagram but ridges not sorted by point id
689 */
690 int qh_eachvoronoi_all (FILE *fp, printvridgeT printvridge, boolT isupper, qh_RIDGE innerouter, boolT inorder) {
691   facetT *facet;
692   vertexT *vertex;
693   int numcenters= 1;  /* vertex 0 is vertex-at-infinity */
694   int totridges= 0;
695
696   qh_clearcenters (qh_ASvoronoi);
697   qh_vertexneighbors();
698   maximize_(qh visit_id, (unsigned) qh num_facets);
699   FORALLfacets {
700     facet->visitid= 0;
701     facet->seen= False;
702     facet->seen2= True;
703   }
704   FORALLfacets {
705     if (facet->upperdelaunay == isupper)
706       facet->visitid= numcenters++;
707   }
708   FORALLvertices 
709     vertex->seen= False;
710   FORALLvertices {
711     if (qh GOODvertex > 0 && qh_pointid(vertex->point)+1 != qh GOODvertex)
712       continue;
713     totridges += qh_eachvoronoi (fp, printvridge, vertex, 
714                    !qh_ALL, innerouter, inorder);
715   }
716   return totridges;
717 } /* eachvoronoi_all */
718       
719 /*-<a                             href="qh-io.htm#TOC"
720   >-------------------------------</a><a name="facet2point">-</a>
721   
722   qh_facet2point( facet, point0, point1, mindist )
723     return two projected temporary vertices for a 2-d facet
724     may be non-simplicial
725
726   returns:
727     point0 and point1 oriented and projected to the facet
728     returns mindist (maximum distance below plane)
729 */
730 void qh_facet2point(facetT *facet, pointT **point0, pointT **point1, realT *mindist) {
731   vertexT *vertex0, *vertex1;
732   realT dist;
733   
734   if (facet->toporient ^ qh_ORIENTclock) {
735     vertex0= SETfirstt_(facet->vertices, vertexT);
736     vertex1= SETsecondt_(facet->vertices, vertexT);
737   }else {
738     vertex1= SETfirstt_(facet->vertices, vertexT);
739     vertex0= SETsecondt_(facet->vertices, vertexT);
740   }
741   zadd_(Zdistio, 2);
742   qh_distplane(vertex0->point, facet, &dist);
743   *mindist= dist;
744   *point0= qh_projectpoint(vertex0->point, facet, dist);
745   qh_distplane(vertex1->point, facet, &dist);
746   minimize_(*mindist, dist);            
747   *point1= qh_projectpoint(vertex1->point, facet, dist);
748 } /* facet2point */
749
750
751 /*-<a                             href="qh-io.htm#TOC"
752   >-------------------------------</a><a name="facetvertices">-</a>
753   
754   qh_facetvertices( facetlist, facets, allfacets )
755     returns temporary set of vertices in a set and/or list of facets
756     if allfacets, ignores qh_skipfacet()
757
758   returns:
759     vertices with qh.vertex_visit
760     
761   notes:
762     optimized for allfacets of facet_list
763
764   design:
765     if allfacets of facet_list
766       create vertex set from vertex_list
767     else
768       for each selected facet in facets or facetlist
769         append unvisited vertices to vertex set
770 */
771 setT *qh_facetvertices (facetT *facetlist, setT *facets, boolT allfacets) {
772   setT *vertices;
773   facetT *facet, **facetp;
774   vertexT *vertex, **vertexp;
775
776   qh vertex_visit++;
777   if (facetlist == qh facet_list && allfacets && !facets) {
778     vertices= qh_settemp (qh num_vertices);
779     FORALLvertices {
780       vertex->visitid= qh vertex_visit; 
781       qh_setappend (&vertices, vertex);
782     }
783   }else {
784     vertices= qh_settemp (qh TEMPsize);
785     FORALLfacet_(facetlist) {
786       if (!allfacets && qh_skipfacet (facet))
787         continue;
788       FOREACHvertex_(facet->vertices) {
789         if (vertex->visitid != qh vertex_visit) {
790           vertex->visitid= qh vertex_visit;
791           qh_setappend (&vertices, vertex);
792         }
793       }
794     }
795   }
796   FOREACHfacet_(facets) {
797     if (!allfacets && qh_skipfacet (facet))
798       continue;
799     FOREACHvertex_(facet->vertices) {
800       if (vertex->visitid != qh vertex_visit) {
801         vertex->visitid= qh vertex_visit;
802         qh_setappend (&vertices, vertex);
803       }
804     }
805   }
806   return vertices;
807 } /* facetvertices */
808
809 /*-<a                             href="qh-geom.htm#TOC"
810   >-------------------------------</a><a name="geomplanes">-</a>
811   
812   qh_geomplanes( facet, outerplane, innerplane )
813     return outer and inner planes for Geomview 
814     qh.PRINTradius is size of vertices and points (includes qh.JOGGLEmax)
815
816   notes:
817     assume precise calculations in io.c with roundoff covered by qh_GEOMepsilon
818 */
819 void qh_geomplanes (facetT *facet, realT *outerplane, realT *innerplane) {
820   realT radius;
821
822   if (qh MERGING || qh JOGGLEmax < REALmax/2) {
823     qh_outerinner (facet, outerplane, innerplane);
824     radius= qh PRINTradius;
825     if (qh JOGGLEmax < REALmax/2)
826       radius -= qh JOGGLEmax * sqrt (qh hull_dim);  /* already accounted for in qh_outerinner() */
827     *outerplane += radius;
828     *innerplane -= radius;
829     if (qh PRINTcoplanar || qh PRINTspheres) {
830       *outerplane += qh MAXabs_coord * qh_GEOMepsilon;
831       *innerplane -= qh MAXabs_coord * qh_GEOMepsilon;
832     }
833   }else 
834     *innerplane= *outerplane= 0;
835 } /* geomplanes */
836
837
838 /*-<a                             href="qh-io.htm#TOC"
839   >-------------------------------</a><a name="markkeep">-</a>
840   
841   qh_markkeep( facetlist )
842     mark good facets that meet qh.KEEParea, qh.KEEPmerge, and qh.KEEPminArea
843     ignores visible facets (not part of convex hull)
844
845   returns:
846     may clear facet->good
847     recomputes qh.num_good
848
849   design:
850     get set of good facets
851     if qh.KEEParea
852       sort facets by area
853       clear facet->good for all but n largest facets
854     if qh.KEEPmerge
855       sort facets by merge count
856       clear facet->good for all but n most merged facets
857     if qh.KEEPminarea
858       clear facet->good if area too small
859     update qh.num_good    
860 */
861 void qh_markkeep (facetT *facetlist) {
862   facetT *facet, **facetp;
863   setT *facets= qh_settemp (qh num_facets);
864   int size, count;
865
866   trace2((qh ferr, "qh_markkeep: only keep %d largest and/or %d most merged facets and/or min area %.2g\n",
867           qh KEEParea, qh KEEPmerge, qh KEEPminArea));
868   FORALLfacet_(facetlist) {
869     if (!facet->visible && facet->good)
870       qh_setappend (&facets, facet);
871   }
872   size= qh_setsize (facets);
873   if (qh KEEParea) {
874     qsort (SETaddr_(facets, facetT), size,
875              sizeof (facetT *), qh_compare_facetarea);
876     if ((count= size - qh KEEParea) > 0) {
877       FOREACHfacet_(facets) {
878         facet->good= False;
879         if (--count == 0)
880           break;
881       }
882     }
883   }
884   if (qh KEEPmerge) {
885     qsort (SETaddr_(facets, facetT), size,
886              sizeof (facetT *), qh_compare_facetmerge);
887     if ((count= size - qh KEEPmerge) > 0) {
888       FOREACHfacet_(facets) {
889         facet->good= False;
890         if (--count == 0)
891           break;
892       }
893     }
894   }
895   if (qh KEEPminArea < REALmax/2) {
896     FOREACHfacet_(facets) {
897       if (!facet->isarea || facet->f.area < qh KEEPminArea)
898         facet->good= False;
899     }
900   }
901   qh_settempfree (&facets);
902   count= 0;
903   FORALLfacet_(facetlist) {
904     if (facet->good)
905       count++;
906   }
907   qh num_good= count;
908 } /* markkeep */
909
910
911 /*-<a                             href="qh-io.htm#TOC"
912   >-------------------------------</a><a name="markvoronoi">-</a>
913   
914   qh_markvoronoi( facetlist, facets, printall, islower, numcenters )
915     mark voronoi vertices for printing by site pairs
916   
917   returns:
918     temporary set of vertices indexed by pointid
919     islower set if printing lower hull (i.e., at least one facet is lower hull)
920     numcenters= total number of Voronoi vertices
921     bumps qh.printoutnum for vertex-at-infinity
922     clears all facet->seen and sets facet->seen2
923     
924     if selected
925       facet->visitid= Voronoi vertex id
926     else if upper hull (or 'Qu' and lower hull)
927       facet->visitid= 0
928     else
929       facet->visitid >= qh num_facets
930   
931   notes:
932     ignores qh.ATinfinity, if defined
933 */
934 setT *qh_markvoronoi (facetT *facetlist, setT *facets, boolT printall, boolT *islowerp, int *numcentersp) {
935   int numcenters=0;
936   facetT *facet, **facetp;
937   setT *vertices;
938   boolT islower= False;
939
940   qh printoutnum++;
941   qh_clearcenters (qh_ASvoronoi);  /* in case, qh_printvdiagram2 called by user */
942   qh_vertexneighbors();
943   vertices= qh_pointvertex();
944   if (qh ATinfinity) 
945     SETelem_(vertices, qh num_points-1)= NULL;
946   qh visit_id++;
947   maximize_(qh visit_id, (unsigned) qh num_facets);
948   FORALLfacet_(facetlist) { 
949     if (printall || !qh_skipfacet (facet)) {
950       if (!facet->upperdelaunay) {
951         islower= True;
952         break;
953       }
954     }
955   }
956   FOREACHfacet_(facets) {
957     if (printall || !qh_skipfacet (facet)) {
958       if (!facet->upperdelaunay) {
959         islower= True;
960         break;
961       }
962     }
963   }
964   FORALLfacets {
965     if (facet->normal && (facet->upperdelaunay == islower))
966       facet->visitid= 0;  /* facetlist or facets may overwrite */
967     else
968       facet->visitid= qh visit_id;
969     facet->seen= False;
970     facet->seen2= True;
971   }
972   numcenters++;  /* qh_INFINITE */
973   FORALLfacet_(facetlist) {
974     if (printall || !qh_skipfacet (facet))
975       facet->visitid= numcenters++;
976   }
977   FOREACHfacet_(facets) {
978     if (printall || !qh_skipfacet (facet))
979       facet->visitid= numcenters++;  
980   }
981   *islowerp= islower;
982   *numcentersp= numcenters;
983   trace2((qh ferr, "qh_markvoronoi: islower %d numcenters %d\n", islower, numcenters));
984   return vertices;
985 } /* markvoronoi */
986
987 /*-<a                             href="qh-io.htm#TOC"
988   >-------------------------------</a><a name="order_vertexneighbors">-</a>
989   
990   qh_order_vertexneighbors( vertex )
991     order facet neighbors of a 2-d or 3-d vertex by adjacency
992
993   notes:
994     does not orient the neighbors
995
996   design:
997     initialize a new neighbor set with the first facet in vertex->neighbors
998     while vertex->neighbors non-empty
999       select next neighbor in the previous facet's neighbor set
1000     set vertex->neighbors to the new neighbor set
1001 */
1002 void qh_order_vertexneighbors(vertexT *vertex) {
1003   setT *newset;
1004   facetT *facet, *neighbor, **neighborp;
1005
1006   trace4((qh ferr, "qh_order_vertexneighbors: order neighbors of v%d for 3-d\n", vertex->id));
1007   newset= qh_settemp (qh_setsize (vertex->neighbors));
1008   facet= (facetT*)qh_setdellast (vertex->neighbors);
1009   qh_setappend (&newset, facet);
1010   while (qh_setsize (vertex->neighbors)) {
1011     FOREACHneighbor_(vertex) {
1012       if (qh_setin (facet->neighbors, neighbor)) {
1013         qh_setdel(vertex->neighbors, neighbor);
1014         qh_setappend (&newset, neighbor);
1015         facet= neighbor;
1016         break;
1017       }
1018     }
1019     if (!neighbor) {
1020       fprintf (qh ferr, "qhull internal error (qh_order_vertexneighbors): no neighbor of v%d for f%d\n",
1021         vertex->id, facet->id);
1022       qh_errexit (qh_ERRqhull, facet, NULL);
1023     }
1024   }
1025   qh_setfree (&vertex->neighbors);
1026   qh_settemppop ();
1027   vertex->neighbors= newset;
1028 } /* order_vertexneighbors */
1029
1030 /*-<a                             href="qh-io.htm#TOC"
1031   >-------------------------------</a><a name="printafacet">-</a>
1032   
1033   qh_printafacet( fp, format, facet, printall )
1034     print facet to fp in given output format (see qh.PRINTout)
1035
1036   returns:
1037     nop if !printall and qh_skipfacet()
1038     nop if visible facet and NEWfacets and format != PRINTfacets
1039     must match qh_countfacets
1040
1041   notes
1042     preserves qh.visit_id
1043     facet->normal may be null if PREmerge/MERGEexact and STOPcone before merge
1044
1045   see
1046     qh_printbegin() and qh_printend()
1047
1048   design:
1049     test for printing facet
1050     call appropriate routine for format
1051     or output results directly
1052 */
1053 void qh_printafacet(FILE *fp, int format, facetT *facet, boolT printall) {
1054   realT color[4], offset, dist, outerplane, innerplane;
1055   boolT zerodiv;
1056   coordT *point, *normp, *coordp, **pointp, *feasiblep;
1057   int k;
1058   vertexT *vertex, **vertexp;
1059   facetT *neighbor, **neighborp;
1060
1061   if (!printall && qh_skipfacet (facet))
1062     return;
1063   if (facet->visible && qh NEWfacets && format != qh_PRINTfacets)
1064     return;
1065   qh printoutnum++;
1066   switch (format) {
1067   case qh_PRINTarea:
1068     if (facet->isarea) {
1069       fprintf (fp, qh_REAL_1, facet->f.area);
1070       fprintf (fp, "\n");
1071     }else
1072       fprintf (fp, "0\n");
1073     break;
1074   case qh_PRINTcoplanars:
1075     fprintf (fp, "%d", qh_setsize (facet->coplanarset));
1076     FOREACHpoint_(facet->coplanarset)
1077       fprintf (fp, " %d", qh_pointid (point));
1078     fprintf (fp, "\n");
1079     break;
1080   case qh_PRINTcentrums:
1081     qh_printcenter (fp, format, NULL, facet);
1082     break;
1083   case qh_PRINTfacets:
1084     qh_printfacet (fp, facet);
1085     break;
1086   case qh_PRINTfacets_xridge:
1087     qh_printfacetheader (fp, facet);
1088     break;
1089   case qh_PRINTgeom:  /* either 2 , 3, or 4-d by qh_printbegin */
1090     if (!facet->normal)
1091       break;
1092     for (k= qh hull_dim; k--; ) {
1093       color[k]= (facet->normal[k]+1.0)/2.0;
1094       maximize_(color[k], -1.0);
1095       minimize_(color[k], +1.0);
1096     }
1097     qh_projectdim3 (color, color);
1098     if (qh PRINTdim != qh hull_dim)
1099       qh_normalize2 (color, 3, True, NULL, NULL);
1100     if (qh hull_dim <= 2)
1101       qh_printfacet2geom (fp, facet, color);
1102     else if (qh hull_dim == 3) {
1103       if (facet->simplicial)
1104         qh_printfacet3geom_simplicial (fp, facet, color);
1105       else
1106         qh_printfacet3geom_nonsimplicial (fp, facet, color);
1107     }else {
1108       if (facet->simplicial)
1109         qh_printfacet4geom_simplicial (fp, facet, color);
1110       else
1111         qh_printfacet4geom_nonsimplicial (fp, facet, color);
1112     }
1113     break;
1114   case qh_PRINTids:
1115     fprintf (fp, "%d\n", facet->id);
1116     break;
1117   case qh_PRINTincidences:
1118   case qh_PRINToff:
1119   case qh_PRINTtriangles:
1120     if (qh hull_dim == 3 && format != qh_PRINTtriangles) 
1121       qh_printfacet3vertex (fp, facet, format);
1122     else if (facet->simplicial || qh hull_dim == 2 || format == qh_PRINToff)
1123       qh_printfacetNvertex_simplicial (fp, facet, format);
1124     else
1125       qh_printfacetNvertex_nonsimplicial (fp, facet, qh printoutvar++, format);
1126     break;
1127   case qh_PRINTinner:
1128     qh_outerinner (facet, NULL, &innerplane);
1129     offset= facet->offset - innerplane;
1130     goto LABELprintnorm;
1131     break; /* prevent warning */
1132   case qh_PRINTmerges:
1133     fprintf (fp, "%d\n", facet->nummerge);
1134     break;
1135   case qh_PRINTnormals:
1136     offset= facet->offset;
1137     goto LABELprintnorm;
1138     break; /* prevent warning */
1139   case qh_PRINTouter:
1140     qh_outerinner (facet, &outerplane, NULL);
1141     offset= facet->offset - outerplane;
1142   LABELprintnorm:
1143     if (!facet->normal) {
1144       fprintf (fp, "no normal for facet f%d\n", facet->id);
1145       break;
1146     }
1147     if (qh CDDoutput) {
1148       fprintf (fp, qh_REAL_1, -offset);
1149       for (k=0; k < qh hull_dim; k++) 
1150         fprintf (fp, qh_REAL_1, -facet->normal[k]);
1151     }else {
1152       for (k=0; k < qh hull_dim; k++) 
1153         fprintf (fp, qh_REAL_1, facet->normal[k]);
1154       fprintf (fp, qh_REAL_1, offset);
1155     }
1156     fprintf (fp, "\n");
1157     break;
1158   case qh_PRINTmathematica:  /* either 2 or 3-d by qh_printbegin */
1159     if (qh hull_dim == 2)
1160       qh_printfacet2math (fp, facet, qh printoutvar++);
1161     else 
1162       qh_printfacet3math (fp, facet, qh printoutvar++);
1163     break;
1164   case qh_PRINTneighbors:
1165     fprintf (fp, "%d", qh_setsize (facet->neighbors));
1166     FOREACHneighbor_(facet)
1167       fprintf (fp, " %d", 
1168                neighbor->visitid ? neighbor->visitid - 1: - neighbor->id);
1169     fprintf (fp, "\n");
1170     break;
1171   case qh_PRINTpointintersect:
1172     if (!qh feasible_point) {
1173       fprintf (fp, "qhull input error (qh_printafacet): option 'Fp' needs qh feasible_point\n");
1174       qh_errexit( qh_ERRinput, NULL, NULL);
1175     }
1176     if (facet->offset > 0)
1177       goto LABELprintinfinite;
1178     point= coordp= (coordT*)qh_memalloc (qh normal_size);
1179     normp= facet->normal;
1180     feasiblep= qh feasible_point;
1181     if (facet->offset < -qh MINdenom) {
1182       for (k= qh hull_dim; k--; )
1183         *(coordp++)= (*(normp++) / - facet->offset) + *(feasiblep++);
1184     }else {
1185       for (k= qh hull_dim; k--; ) {
1186         *(coordp++)= qh_divzero (*(normp++), facet->offset, qh MINdenom_1,
1187                                  &zerodiv) + *(feasiblep++);
1188         if (zerodiv) {
1189           qh_memfree (point, qh normal_size);
1190           goto LABELprintinfinite;
1191         }
1192       }
1193     }
1194     qh_printpoint (fp, NULL, point);
1195     qh_memfree (point, qh normal_size);
1196     break;
1197   LABELprintinfinite:
1198     for (k= qh hull_dim; k--; )
1199       fprintf (fp, qh_REAL_1, qh_INFINITE);
1200     fprintf (fp, "\n");   
1201     break;
1202   case qh_PRINTpointnearest:
1203     FOREACHpoint_(facet->coplanarset) {
1204       int id, id2;
1205       vertex= qh_nearvertex (facet, point, &dist);
1206       id= qh_pointid (vertex->point);
1207       id2= qh_pointid (point);
1208       fprintf (fp, "%d %d %d " qh_REAL_1 "\n", id, id2, facet->id, dist);
1209     }
1210     break;
1211   case qh_PRINTpoints:  /* VORONOI only by qh_printbegin */
1212     if (qh CDDoutput)
1213       fprintf (fp, "1 ");
1214     qh_printcenter (fp, format, NULL, facet);
1215     break;
1216   case qh_PRINTvertices:
1217     fprintf (fp, "%d", qh_setsize (facet->vertices));
1218     FOREACHvertex_(facet->vertices)
1219       fprintf (fp, " %d", qh_pointid (vertex->point));
1220     fprintf (fp, "\n");
1221     break;
1222   }
1223 } /* printafacet */
1224
1225 /*-<a                             href="qh-io.htm#TOC"
1226   >-------------------------------</a><a name="printbegin">-</a>
1227   
1228   qh_printbegin(  )
1229     prints header for all output formats
1230
1231   returns:
1232     checks for valid format
1233   
1234   notes:
1235     uses qh.visit_id for 3/4off
1236     changes qh.interior_point if printing centrums
1237     qh_countfacets clears facet->visitid for non-good facets
1238     
1239   see
1240     qh_printend() and qh_printafacet()
1241     
1242   design:
1243     count facets and related statistics
1244     print header for format
1245 */
1246 void qh_printbegin (FILE *fp, int format, facetT *facetlist, setT *facets, boolT printall) {
1247   int numfacets, numsimplicial, numridges, totneighbors, numcoplanars, numtricoplanars;
1248   int i, num;
1249   facetT *facet, **facetp;
1250   vertexT *vertex, **vertexp;
1251   setT *vertices;
1252   pointT *point, **pointp, *pointtemp;
1253
1254   qh printoutnum= 0;
1255   qh_countfacets (facetlist, facets, printall, &numfacets, &numsimplicial, 
1256       &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
1257   switch (format) {
1258   case qh_PRINTnone:
1259     break;
1260   case qh_PRINTarea:
1261     fprintf (fp, "%d\n", numfacets);
1262     break;
1263   case qh_PRINTcoplanars:
1264     fprintf (fp, "%d\n", numfacets);
1265     break;
1266   case qh_PRINTcentrums:
1267     if (qh CENTERtype == qh_ASnone)
1268       qh_clearcenters (qh_AScentrum);
1269     fprintf (fp, "%d\n%d\n", qh hull_dim, numfacets);
1270     break;
1271   case qh_PRINTfacets:
1272   case qh_PRINTfacets_xridge:
1273     if (facetlist)
1274       qh_printvertexlist (fp, "Vertices and facets:\n", facetlist, facets, printall);
1275     break;
1276   case qh_PRINTgeom: 
1277     if (qh hull_dim > 4)  /* qh_initqhull_globals also checks */
1278       goto LABELnoformat;
1279     if (qh VORONOI && qh hull_dim > 3)  /* PRINTdim == DROPdim == hull_dim-1 */
1280       goto LABELnoformat;
1281     if (qh hull_dim == 2 && (qh PRINTridges || qh DOintersections))
1282       fprintf (qh ferr, "qhull warning: output for ridges and intersections not implemented in 2-d\n");
1283     if (qh hull_dim == 4 && (qh PRINTinner || qh PRINTouter ||
1284                              (qh PRINTdim == 4 && qh PRINTcentrums)))
1285       fprintf (qh ferr, "qhull warning: output for outer/inner planes and centrums not implemented in 4-d\n");
1286     if (qh PRINTdim == 4 && (qh PRINTspheres))
1287       fprintf (qh ferr, "qhull warning: output for vertices not implemented in 4-d\n");
1288     if (qh PRINTdim == 4 && qh DOintersections && qh PRINTnoplanes)
1289       fprintf (qh ferr, "qhull warning: 'Gnh' generates no output in 4-d\n");
1290     if (qh PRINTdim == 2) {
1291       fprintf(fp, "{appearance {linewidth 3} LIST # %s | %s\n",
1292               qh rbox_command, qh qhull_command);
1293     }else if (qh PRINTdim == 3) {
1294       fprintf(fp, "{appearance {+edge -evert linewidth 2} LIST # %s | %s\n",
1295               qh rbox_command, qh qhull_command);
1296     }else if (qh PRINTdim == 4) {
1297       qh visit_id++;
1298       num= 0;
1299       FORALLfacet_(facetlist)    /* get number of ridges to be printed */
1300         qh_printend4geom (NULL, facet, &num, printall);
1301       FOREACHfacet_(facets)
1302         qh_printend4geom (NULL, facet, &num, printall);
1303       qh ridgeoutnum= num;
1304       qh printoutvar= 0;  /* counts number of ridges in output */
1305       fprintf (fp, "LIST # %s | %s\n", qh rbox_command, qh qhull_command);
1306     }
1307     if (qh PRINTdots) {
1308       qh printoutnum++;
1309       num= qh num_points + qh_setsize (qh other_points);
1310       if (qh DELAUNAY && qh ATinfinity)
1311         num--;
1312       if (qh PRINTdim == 4)
1313         fprintf (fp, "4VECT %d %d 1\n", num, num);
1314       else
1315         fprintf (fp, "VECT %d %d 1\n", num, num);
1316       for (i= num; i--; ) {
1317         if (i % 20 == 0)
1318           fprintf (fp, "\n");
1319         fprintf (fp, "1 ");
1320       }
1321       fprintf (fp, "# 1 point per line\n1 ");
1322       for (i= num-1; i--; ) {
1323         if (i % 20 == 0)
1324           fprintf (fp, "\n");
1325         fprintf (fp, "0 ");
1326       }
1327       fprintf (fp, "# 1 color for all\n");
1328       FORALLpoints {
1329         if (!qh DELAUNAY || !qh ATinfinity || qh_pointid(point) != qh num_points-1) {
1330           if (qh PRINTdim == 4)
1331             qh_printpoint (fp, NULL, point);
1332           else
1333             qh_printpoint3 (fp, point);
1334         }
1335       }
1336       FOREACHpoint_(qh other_points) {
1337         if (qh PRINTdim == 4)
1338           qh_printpoint (fp, NULL, point);
1339         else
1340           qh_printpoint3 (fp, point);
1341       }
1342       fprintf (fp, "0 1 1 1  # color of points\n");
1343     }
1344     if (qh PRINTdim == 4  && !qh PRINTnoplanes)
1345       /* 4dview loads up multiple 4OFF objects slowly */
1346       fprintf(fp, "4OFF %d %d 1\n", 3*qh ridgeoutnum, qh ridgeoutnum);
1347     qh PRINTcradius= 2 * qh DISTround;  /* include test DISTround */
1348     if (qh PREmerge) {
1349       maximize_(qh PRINTcradius, qh premerge_centrum + qh DISTround);
1350     }else if (qh POSTmerge)
1351       maximize_(qh PRINTcradius, qh postmerge_centrum + qh DISTround);
1352     qh PRINTradius= qh PRINTcradius;
1353     if (qh PRINTspheres + qh PRINTcoplanar)
1354       maximize_(qh PRINTradius, qh MAXabs_coord * qh_MINradius);
1355     if (qh premerge_cos < REALmax/2) {
1356       maximize_(qh PRINTradius, (1- qh premerge_cos) * qh MAXabs_coord);
1357     }else if (!qh PREmerge && qh POSTmerge && qh postmerge_cos < REALmax/2) {
1358       maximize_(qh PRINTradius, (1- qh postmerge_cos) * qh MAXabs_coord);
1359     }
1360     maximize_(qh PRINTradius, qh MINvisible); 
1361     if (qh JOGGLEmax < REALmax/2)
1362       qh PRINTradius += qh JOGGLEmax * sqrt (qh hull_dim);
1363     if (qh PRINTdim != 4 &&
1364         (qh PRINTcoplanar || qh PRINTspheres || qh PRINTcentrums)) {
1365       vertices= qh_facetvertices (facetlist, facets, printall);
1366       if (qh PRINTspheres && qh PRINTdim <= 3)
1367          qh_printspheres (fp, vertices, qh PRINTradius);
1368       if (qh PRINTcoplanar || qh PRINTcentrums) {
1369         qh firstcentrum= True;
1370         if (qh PRINTcoplanar&& !qh PRINTspheres) {
1371           FOREACHvertex_(vertices) 
1372             qh_printpointvect2 (fp, vertex->point, NULL,
1373                                 qh interior_point, qh PRINTradius);
1374         }
1375         FORALLfacet_(facetlist) {
1376           if (!printall && qh_skipfacet(facet))
1377             continue;
1378           if (!facet->normal)
1379             continue;
1380           if (qh PRINTcentrums && qh PRINTdim <= 3)
1381             qh_printcentrum (fp, facet, qh PRINTcradius);
1382           if (!qh PRINTcoplanar)
1383             continue;
1384           FOREACHpoint_(facet->coplanarset)
1385             qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius);
1386           FOREACHpoint_(facet->outsideset)
1387             qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius);
1388         }
1389         FOREACHfacet_(facets) {
1390           if (!printall && qh_skipfacet(facet))
1391             continue;
1392           if (!facet->normal)
1393             continue;
1394           if (qh PRINTcentrums && qh PRINTdim <= 3)
1395             qh_printcentrum (fp, facet, qh PRINTcradius);
1396           if (!qh PRINTcoplanar)
1397             continue;
1398           FOREACHpoint_(facet->coplanarset)
1399             qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius);
1400           FOREACHpoint_(facet->outsideset)
1401             qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius);
1402         }
1403       }
1404       qh_settempfree (&vertices);
1405     }
1406     qh visit_id++; /* for printing hyperplane intersections */
1407     break;
1408   case qh_PRINTids:
1409     fprintf (fp, "%d\n", numfacets);
1410     break;
1411   case qh_PRINTincidences:
1412     if (qh VORONOI && qh PRINTprecision)
1413       fprintf (qh ferr, "qhull warning: writing Delaunay.  Use 'p' or 'o' for Voronoi centers\n");
1414     qh printoutvar= qh vertex_id;  /* centrum id for non-simplicial facets */
1415     if (qh hull_dim <= 3)
1416       fprintf(fp, "%d\n", numfacets);
1417     else
1418       fprintf(fp, "%d\n", numsimplicial+numridges);
1419     break;
1420   case qh_PRINTinner:
1421   case qh_PRINTnormals:
1422   case qh_PRINTouter:
1423     if (qh CDDoutput)
1424       fprintf (fp, "%s | %s\nbegin\n    %d %d real\n", qh rbox_command, 
1425               qh qhull_command, numfacets, qh hull_dim+1);
1426     else
1427       fprintf (fp, "%d\n%d\n", qh hull_dim+1, numfacets);
1428     break;
1429   case qh_PRINTmathematica:  
1430     if (qh hull_dim > 3)  /* qh_initbuffers also checks */
1431       goto LABELnoformat;
1432     if (qh VORONOI)
1433       fprintf (qh ferr, "qhull warning: output is the Delaunay triangulation\n");
1434     fprintf(fp, "{\n");
1435     qh printoutvar= 0;   /* counts number of facets for notfirst */
1436     break;
1437   case qh_PRINTmerges:
1438     fprintf (fp, "%d\n", numfacets);
1439     break;
1440   case qh_PRINTpointintersect:
1441     fprintf (fp, "%d\n%d\n", qh hull_dim, numfacets);
1442     break;
1443   case qh_PRINTneighbors:
1444     fprintf (fp, "%d\n", numfacets);
1445     break;
1446   case qh_PRINToff:
1447   case qh_PRINTtriangles:
1448     if (qh VORONOI)
1449       goto LABELnoformat;
1450     num = qh hull_dim;
1451     if (format == qh_PRINToff || qh hull_dim == 2)
1452       fprintf (fp, "%d\n%d %d %d\n", num, 
1453         qh num_points+qh_setsize (qh other_points), numfacets, totneighbors/2);
1454     else { /* qh_PRINTtriangles */
1455       qh printoutvar= qh num_points+qh_setsize (qh other_points); /* first centrum */
1456       if (qh DELAUNAY)
1457         num--;  /* drop last dimension */
1458       fprintf (fp, "%d\n%d %d %d\n", num, qh printoutvar 
1459         + numfacets - numsimplicial, numsimplicial + numridges, totneighbors/2);
1460     }
1461     FORALLpoints
1462       qh_printpointid (qh fout, NULL, num, point, -1);
1463     FOREACHpoint_(qh other_points)
1464       qh_printpointid (qh fout, NULL, num, point, -1);
1465     if (format == qh_PRINTtriangles && qh hull_dim > 2) {
1466       FORALLfacets {
1467         if (!facet->simplicial && facet->visitid)
1468           qh_printcenter (qh fout, format, NULL, facet);
1469       }
1470     }
1471     break;
1472   case qh_PRINTpointnearest:
1473     fprintf (fp, "%d\n", numcoplanars);
1474     break;
1475   case qh_PRINTpoints:
1476     if (!qh VORONOI)
1477       goto LABELnoformat;
1478     if (qh CDDoutput)
1479       fprintf (fp, "%s | %s\nbegin\n%d %d real\n", qh rbox_command,
1480              qh qhull_command, numfacets, qh hull_dim);
1481     else
1482       fprintf (fp, "%d\n%d\n", qh hull_dim-1, numfacets);
1483     break;
1484   case qh_PRINTvertices:
1485     fprintf (fp, "%d\n", numfacets);
1486     break;
1487   case qh_PRINTsummary:
1488   default:
1489   LABELnoformat:
1490     fprintf (qh ferr, "qhull internal error (qh_printbegin): can not use this format for dimension %d\n",
1491          qh hull_dim);
1492     qh_errexit (qh_ERRqhull, NULL, NULL);
1493   }
1494 } /* printbegin */
1495
1496 /*-<a                             href="qh-io.htm#TOC"
1497   >-------------------------------</a><a name="printcenter">-</a>
1498   
1499   qh_printcenter( fp, string, facet )
1500     print facet->center as centrum or Voronoi center
1501     string may be NULL.  Don't include '%' codes.
1502     nop if qh CENTERtype neither CENTERvoronoi nor CENTERcentrum
1503     if upper envelope of Delaunay triangulation and point at-infinity
1504       prints qh_INFINITE instead;
1505
1506   notes:
1507     defines facet->center if needed
1508     if format=PRINTgeom, adds a 0 if would otherwise be 2-d
1509 */
1510 void qh_printcenter (FILE *fp, int format, char *string, facetT *facet) {
1511   int k, num;
1512
1513   if (qh CENTERtype != qh_ASvoronoi && qh CENTERtype != qh_AScentrum)
1514     return;
1515   if (string)
1516     fprintf (fp, string, facet->id);
1517   if (qh CENTERtype == qh_ASvoronoi) {
1518     num= qh hull_dim-1;
1519     if (!facet->normal || !facet->upperdelaunay || !qh ATinfinity) {
1520       if (!facet->center)
1521         facet->center= qh_facetcenter (facet->vertices);
1522       for (k=0; k < num; k++)
1523         fprintf (fp, qh_REAL_1, facet->center[k]);
1524     }else {
1525       for (k=0; k < num; k++)
1526         fprintf (fp, qh_REAL_1, qh_INFINITE);
1527     }
1528   }else /* qh CENTERtype == qh_AScentrum */ {
1529     num= qh hull_dim;
1530     if (format == qh_PRINTtriangles && qh DELAUNAY) 
1531       num--;
1532     if (!facet->center) 
1533       facet->center= qh_getcentrum (facet);
1534     for (k=0; k < num; k++)
1535       fprintf (fp, qh_REAL_1, facet->center[k]);
1536   }
1537   if (format == qh_PRINTgeom && num == 2)
1538     fprintf (fp, " 0\n");
1539   else
1540     fprintf (fp, "\n");
1541 } /* printcenter */
1542
1543 /*-<a                             href="qh-io.htm#TOC"
1544   >-------------------------------</a><a name="printcentrum">-</a>
1545   
1546   qh_printcentrum( fp, facet, radius )
1547     print centrum for a facet in OOGL format
1548     radius defines size of centrum
1549     2-d or 3-d only
1550
1551   returns:
1552     defines facet->center if needed
1553 */
1554 void qh_printcentrum (FILE *fp, facetT *facet, realT radius) {
1555   pointT *centrum, *projpt;
1556   boolT tempcentrum= False;
1557   realT xaxis[4], yaxis[4], normal[4], dist;
1558   realT green[3]={0, 1, 0};
1559   vertexT *apex;
1560   int k;
1561   
1562   if (qh CENTERtype == qh_AScentrum) {
1563     if (!facet->center)
1564       facet->center= qh_getcentrum (facet);
1565     centrum= facet->center;
1566   }else {
1567     centrum= qh_getcentrum (facet);
1568     tempcentrum= True;
1569   }
1570   fprintf (fp, "{appearance {-normal -edge normscale 0} ");
1571   if (qh firstcentrum) {
1572     qh firstcentrum= False;
1573     fprintf (fp, "{INST geom { define centrum CQUAD  # f%d\n\
1574 -0.3 -0.3 0.0001     0 0 1 1\n\
1575  0.3 -0.3 0.0001     0 0 1 1\n\
1576  0.3  0.3 0.0001     0 0 1 1\n\
1577 -0.3  0.3 0.0001     0 0 1 1 } transform { \n", facet->id);
1578   }else
1579     fprintf (fp, "{INST geom { : centrum } transform { # f%d\n", facet->id);
1580   apex= SETfirstt_(facet->vertices, vertexT);
1581   qh_distplane(apex->point, facet, &dist);
1582   projpt= qh_projectpoint(apex->point, facet, dist);
1583   for (k= qh hull_dim; k--; ) {
1584     xaxis[k]= projpt[k] - centrum[k];
1585     normal[k]= facet->normal[k];
1586   }
1587   if (qh hull_dim == 2) {
1588     xaxis[2]= 0;
1589     normal[2]= 0;
1590   }else if (qh hull_dim == 4) {
1591     qh_projectdim3 (xaxis, xaxis);
1592     qh_projectdim3 (normal, normal);
1593     qh_normalize2 (normal, qh PRINTdim, True, NULL, NULL);
1594   }
1595   qh_crossproduct (3, xaxis, normal, yaxis);
1596   fprintf (fp, "%8.4g %8.4g %8.4g 0\n", xaxis[0], xaxis[1], xaxis[2]);
1597   fprintf (fp, "%8.4g %8.4g %8.4g 0\n", yaxis[0], yaxis[1], yaxis[2]);
1598   fprintf (fp, "%8.4g %8.4g %8.4g 0\n", normal[0], normal[1], normal[2]);
1599   qh_printpoint3 (fp, centrum);
1600   fprintf (fp, "1 }}}\n"); 
1601   qh_memfree (projpt, qh normal_size);
1602   qh_printpointvect (fp, centrum, facet->normal, NULL, radius, green);
1603   if (tempcentrum)
1604     qh_memfree (centrum, qh normal_size);
1605 } /* printcentrum */
1606   
1607 /*-<a                             href="qh-io.htm#TOC"
1608   >-------------------------------</a><a name="printend">-</a>
1609   
1610   qh_printend( fp, format )
1611     prints trailer for all output formats
1612
1613   see:
1614     qh_printbegin() and qh_printafacet()
1615       
1616 */
1617 void qh_printend (FILE *fp, int format, facetT *facetlist, setT *facets, boolT printall) {
1618   int num;
1619   facetT *facet, **facetp;
1620
1621   if (!qh printoutnum)
1622     fprintf (qh ferr, "qhull warning: no facets printed\n");
1623   switch (format) {
1624   case qh_PRINTgeom:
1625     if (qh hull_dim == 4 && qh DROPdim < 0  && !qh PRINTnoplanes) {
1626       qh visit_id++;
1627       num= 0;
1628       FORALLfacet_(facetlist)
1629         qh_printend4geom (fp, facet,&num, printall);
1630       FOREACHfacet_(facets) 
1631         qh_printend4geom (fp, facet, &num, printall);
1632       if (num != qh ridgeoutnum || qh printoutvar != qh ridgeoutnum) {
1633         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);
1634         qh_errexit (qh_ERRqhull, NULL, NULL);
1635       }
1636     }else
1637       fprintf(fp, "}\n");
1638     break;
1639   case qh_PRINTinner:
1640   case qh_PRINTnormals:
1641   case qh_PRINTouter:
1642     if (qh CDDoutput) 
1643       fprintf (fp, "end\n");
1644     break;
1645   case qh_PRINTmathematica:
1646     fprintf(fp, "}\n");
1647     break;
1648   case qh_PRINTpoints:
1649     if (qh CDDoutput)
1650       fprintf (fp, "end\n");
1651     break;
1652   }
1653 } /* printend */
1654
1655 /*-<a                             href="qh-io.htm#TOC"
1656   >-------------------------------</a><a name="printend4geom">-</a>
1657   
1658   qh_printend4geom( fp, facet, numridges, printall )
1659     helper function for qh_printbegin/printend
1660
1661   returns:
1662     number of printed ridges
1663   
1664   notes:
1665     just counts printed ridges if fp=NULL
1666     uses facet->visitid
1667     must agree with qh_printfacet4geom...
1668
1669   design:
1670     computes color for facet from its normal
1671     prints each ridge of facet 
1672 */
1673 void qh_printend4geom (FILE *fp, facetT *facet, int *nump, boolT printall) {
1674   realT color[3];
1675   int i, num= *nump;
1676   facetT *neighbor, **neighborp;
1677   ridgeT *ridge, **ridgep;
1678   
1679   if (!printall && qh_skipfacet(facet))
1680     return;
1681   if (qh PRINTnoplanes || (facet->visible && qh NEWfacets))
1682     return;
1683   if (!facet->normal)
1684     return;
1685   if (fp) {
1686     for (i=0; i < 3; i++) {
1687       color[i]= (facet->normal[i]+1.0)/2.0;
1688       maximize_(color[i], -1.0);
1689       minimize_(color[i], +1.0);
1690     }
1691   }
1692   facet->visitid= qh visit_id;
1693   if (facet->simplicial) {
1694     FOREACHneighbor_(facet) {
1695       if (neighbor->visitid != qh visit_id) {
1696         if (fp)
1697           fprintf (fp, "3 %d %d %d %8.4g %8.4g %8.4g 1 # f%d f%d\n",
1698                  3*num, 3*num+1, 3*num+2, color[0], color[1], color[2],
1699                  facet->id, neighbor->id);
1700         num++;
1701       }
1702     }
1703   }else {
1704     FOREACHridge_(facet->ridges) {
1705       neighbor= otherfacet_(ridge, facet);
1706       if (neighbor->visitid != qh visit_id) {
1707         if (fp)
1708           fprintf (fp, "3 %d %d %d %8.4g %8.4g %8.4g 1 #r%d f%d f%d\n",
1709                  3*num, 3*num+1, 3*num+2, color[0], color[1], color[2],
1710                  ridge->id, facet->id, neighbor->id);
1711         num++;
1712       }
1713     }
1714   }
1715   *nump= num;
1716 } /* printend4geom */
1717
1718 /*-<a                             href="qh-io.htm#TOC"
1719   >-------------------------------</a><a name="printextremes">-</a>
1720   
1721   qh_printextremes( fp, facetlist, facets, printall )
1722     print extreme points for convex hulls or halfspace intersections
1723
1724   notes:
1725     #points, followed by ids, one per line
1726     
1727     sorted by id
1728     same order as qh_printpoints_out if no coplanar/interior points
1729 */
1730 void qh_printextremes (FILE *fp, facetT *facetlist, setT *facets, int printall) {
1731   setT *vertices, *points;
1732   pointT *point;
1733   vertexT *vertex, **vertexp;
1734   int id;
1735   int numpoints=0, point_i, point_n;
1736   int allpoints= qh num_points + qh_setsize (qh other_points);
1737
1738   points= qh_settemp (allpoints);
1739   qh_setzero (points, 0, allpoints);
1740   vertices= qh_facetvertices (facetlist, facets, printall);
1741   FOREACHvertex_(vertices) {
1742     id= qh_pointid (vertex->point);
1743     if (id >= 0) {
1744       SETelem_(points, id)= vertex->point;
1745       numpoints++;
1746     }
1747   }
1748   qh_settempfree (&vertices);
1749   fprintf (fp, "%d\n", numpoints);
1750   FOREACHpoint_i_(points) {
1751     if (point) 
1752       fprintf (fp, "%d\n", point_i);
1753   }
1754   qh_settempfree (&points);
1755 } /* printextremes */
1756
1757 /*-<a                             href="qh-io.htm#TOC"
1758   >-------------------------------</a><a name="printextremes_2d">-</a>
1759   
1760   qh_printextremes_2d( fp, facetlist, facets, printall )
1761     prints point ids for facets in qh_ORIENTclock order
1762
1763   notes:
1764     #points, followed by ids, one per line
1765     if facetlist/facets are disjoint than the output includes skips
1766     errors if facets form a loop
1767     does not print coplanar points
1768 */
1769 void qh_printextremes_2d (FILE *fp, facetT *facetlist, setT *facets, int printall) {
1770   int numfacets, numridges, totneighbors, numcoplanars, numsimplicial, numtricoplanars;
1771   setT *vertices;
1772   facetT *facet, *startfacet, *nextfacet;
1773   vertexT *vertexA, *vertexB;
1774
1775   qh_countfacets (facetlist, facets, printall, &numfacets, &numsimplicial, 
1776       &totneighbors, &numridges, &numcoplanars, &numtricoplanars); /* marks qh visit_id */
1777   vertices= qh_facetvertices (facetlist, facets, printall);
1778   fprintf(fp, "%d\n", qh_setsize (vertices));
1779   qh_settempfree (&vertices);
1780   if (!numfacets)
1781     return;
1782   facet= startfacet= facetlist ? facetlist : SETfirstt_(facets, facetT);
1783   qh vertex_visit++;
1784   qh visit_id++;
1785   do {
1786     if (facet->toporient ^ qh_ORIENTclock) {
1787       vertexA= SETfirstt_(facet->vertices, vertexT);
1788       vertexB= SETsecondt_(facet->vertices, vertexT);
1789       nextfacet= SETfirstt_(facet->neighbors, facetT);
1790     }else {
1791       vertexA= SETsecondt_(facet->vertices, vertexT);
1792       vertexB= SETfirstt_(facet->vertices, vertexT);
1793       nextfacet= SETsecondt_(facet->neighbors, facetT);
1794     }
1795     if (facet->visitid == qh visit_id) {
1796       fprintf(qh ferr, "qh_printextremes_2d: loop in facet list.  facet %d nextfacet %d\n",
1797                  facet->id, nextfacet->id);
1798       qh_errexit2 (qh_ERRqhull, facet, nextfacet);
1799     }
1800     if (facet->visitid) {
1801       if (vertexA->visitid != qh vertex_visit) {
1802         vertexA->visitid= qh vertex_visit;
1803         fprintf(fp, "%d\n", qh_pointid (vertexA->point));
1804       }
1805       if (vertexB->visitid != qh vertex_visit) {
1806         vertexB->visitid= qh vertex_visit;
1807         fprintf(fp, "%d\n", qh_pointid (vertexB->point));
1808       }
1809     }
1810     facet->visitid= qh visit_id;
1811     facet= nextfacet;
1812   }while (facet && facet != startfacet);
1813 } /* printextremes_2d */
1814
1815 /*-<a                             href="qh-io.htm#TOC"
1816   >-------------------------------</a><a name="printextremes_d">-</a>
1817   
1818   qh_printextremes_d( fp, facetlist, facets, printall )
1819     print extreme points of input sites for Delaunay triangulations
1820
1821   notes:
1822     #points, followed by ids, one per line
1823     
1824     unordered
1825 */
1826 void qh_printextremes_d (FILE *fp, facetT *facetlist, setT *facets, int printall) {
1827   setT *vertices;
1828   vertexT *vertex, **vertexp;
1829   boolT upperseen, lowerseen;
1830   facetT *neighbor, **neighborp;
1831   int numpoints=0;
1832
1833   vertices= qh_facetvertices (facetlist, facets, printall);
1834   qh_vertexneighbors();
1835   FOREACHvertex_(vertices) {
1836     upperseen= lowerseen= False;
1837     FOREACHneighbor_(vertex) {
1838       if (neighbor->upperdelaunay)
1839         upperseen= True;
1840       else
1841         lowerseen= True;
1842     }
1843     if (upperseen && lowerseen) {
1844       vertex->seen= True;
1845       numpoints++;
1846     }else
1847       vertex->seen= False;
1848   }
1849   fprintf (fp, "%d\n", numpoints);
1850   FOREACHvertex_(vertices) {
1851     if (vertex->seen)
1852       fprintf (fp, "%d\n", qh_pointid (vertex->point));
1853   }
1854   qh_settempfree (&vertices);
1855 } /* printextremes_d */
1856
1857 /*-<a                             href="qh-io.htm#TOC"
1858   >-------------------------------</a><a name="printfacet">-</a>
1859   
1860   qh_printfacet( fp, facet )
1861     prints all fields of a facet to fp
1862
1863   notes:
1864     ridges printed in neighbor order
1865 */
1866 void qh_printfacet(FILE *fp, facetT *facet) {
1867
1868   qh_printfacetheader (fp, facet);
1869   if (facet->ridges)
1870     qh_printfacetridges (fp, facet);
1871 } /* printfacet */
1872
1873
1874 /*-<a                             href="qh-io.htm#TOC"
1875   >-------------------------------</a><a name="printfacet2geom">-</a>
1876   
1877   qh_printfacet2geom( fp, facet, color )
1878     print facet as part of a 2-d VECT for Geomview
1879   
1880     notes:
1881       assume precise calculations in io.c with roundoff covered by qh_GEOMepsilon
1882       mindist is calculated within io.c.  maxoutside is calculated elsewhere
1883       so a DISTround error may have occured.
1884 */
1885 void qh_printfacet2geom(FILE *fp, facetT *facet, realT color[3]) {
1886   pointT *point0, *point1;
1887   realT mindist, innerplane, outerplane;
1888   int k;
1889
1890   qh_facet2point (facet, &point0, &point1, &mindist);
1891   qh_geomplanes (facet, &outerplane, &innerplane);
1892   if (qh PRINTouter || (!qh PRINTnoplanes && !qh PRINTinner))
1893     qh_printfacet2geom_points(fp, point0, point1, facet, outerplane, color);
1894   if (qh PRINTinner || (!qh PRINTnoplanes && !qh PRINTouter &&
1895                 outerplane - innerplane > 2 * qh MAXabs_coord * qh_GEOMepsilon)) {
1896     for(k= 3; k--; )
1897       color[k]= 1.0 - color[k];
1898     qh_printfacet2geom_points(fp, point0, point1, facet, innerplane, color);
1899   }
1900   qh_memfree (point1, qh normal_size);
1901   qh_memfree (point0, qh normal_size); 
1902 } /* printfacet2geom */
1903
1904 /*-<a                             href="qh-io.htm#TOC"
1905   >-------------------------------</a><a name="printfacet2geom_points">-</a>
1906   
1907   qh_printfacet2geom_points( fp, point1, point2, facet, offset, color )
1908     prints a 2-d facet as a VECT with 2 points at some offset.   
1909     The points are on the facet's plane.
1910 */
1911 void qh_printfacet2geom_points(FILE *fp, pointT *point1, pointT *point2,
1912                                facetT *facet, realT offset, realT color[3]) {
1913   pointT *p1= point1, *p2= point2;
1914
1915   fprintf(fp, "VECT 1 2 1 2 1 # f%d\n", facet->id);
1916   if (offset != 0.0) {
1917     p1= qh_projectpoint (p1, facet, -offset);
1918     p2= qh_projectpoint (p2, facet, -offset);
1919   }
1920   fprintf(fp, "%8.4g %8.4g %8.4g\n%8.4g %8.4g %8.4g\n",
1921            p1[0], p1[1], 0.0, p2[0], p2[1], 0.0);
1922   if (offset != 0.0) {
1923     qh_memfree (p1, qh normal_size);
1924     qh_memfree (p2, qh normal_size);
1925   }
1926   fprintf(fp, "%8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]);
1927 } /* printfacet2geom_points */
1928
1929
1930 /*-<a                             href="qh-io.htm#TOC"
1931   >-------------------------------</a><a name="printfacet2math">-</a>
1932   
1933   qh_printfacet2math( fp, facet, notfirst )
1934     print 2-d Mathematica output for a facet
1935     may be non-simplicial
1936
1937   notes:
1938     use %16.8f since Mathematica 2.2 does not handle exponential format
1939 */
1940 void qh_printfacet2math(FILE *fp, facetT *facet, int notfirst) {
1941   pointT *point0, *point1;
1942   realT mindist;
1943   
1944   qh_facet2point (facet, &point0, &point1, &mindist);
1945   if (notfirst)
1946     fprintf(fp, ",");
1947   fprintf(fp, "Line[{{%16.8f, %16.8f}, {%16.8f, %16.8f}}]\n",
1948           point0[0], point0[1], point1[0], point1[1]);
1949   qh_memfree (point1, qh normal_size);
1950   qh_memfree (point0, qh normal_size);
1951 } /* printfacet2math */
1952
1953
1954 /*-<a                             href="qh-io.htm#TOC"
1955   >-------------------------------</a><a name="printfacet3geom_nonsimplicial">-</a>
1956   
1957   qh_printfacet3geom_nonsimplicial( fp, facet, color )
1958     print Geomview OFF for a 3-d nonsimplicial facet.
1959     if DOintersections, prints ridges to unvisited neighbors (qh visit_id) 
1960
1961   notes
1962     uses facet->visitid for intersections and ridges
1963 */
1964 void qh_printfacet3geom_nonsimplicial(FILE *fp, facetT *facet, realT color[3]) {
1965   ridgeT *ridge, **ridgep;
1966   setT *projectedpoints, *vertices;
1967   vertexT *vertex, **vertexp, *vertexA, *vertexB;
1968   pointT *projpt, *point, **pointp;
1969   facetT *neighbor;
1970   realT dist, outerplane, innerplane;
1971   int cntvertices, k;
1972   realT black[3]={0, 0, 0}, green[3]={0, 1, 0};
1973
1974   qh_geomplanes (facet, &outerplane, &innerplane); 
1975   vertices= qh_facet3vertex (facet); /* oriented */
1976   cntvertices= qh_setsize(vertices);
1977   projectedpoints= qh_settemp(cntvertices);
1978   FOREACHvertex_(vertices) {
1979     zinc_(Zdistio);
1980     qh_distplane(vertex->point, facet, &dist);
1981     projpt= qh_projectpoint(vertex->point, facet, dist);
1982     qh_setappend (&projectedpoints, projpt);
1983   }
1984   if (qh PRINTouter || (!qh PRINTnoplanes && !qh PRINTinner))
1985     qh_printfacet3geom_points(fp, projectedpoints, facet, outerplane, color);
1986   if (qh PRINTinner || (!qh PRINTnoplanes && !qh PRINTouter &&
1987                 outerplane - innerplane > 2 * qh MAXabs_coord * qh_GEOMepsilon)) {
1988     for (k=3; k--; )
1989       color[k]= 1.0 - color[k];
1990     qh_printfacet3geom_points(fp, projectedpoints, facet, innerplane, color);
1991   }
1992   FOREACHpoint_(projectedpoints)
1993     qh_memfree (point, qh normal_size);
1994   qh_settempfree(&projectedpoints);
1995   qh_settempfree(&vertices);
1996   if ((qh DOintersections || qh PRINTridges)
1997   && (!facet->visible || !qh NEWfacets)) {
1998     facet->visitid= qh visit_id;
1999     FOREACHridge_(facet->ridges) {
2000       neighbor= otherfacet_(ridge, facet);
2001       if (neighbor->visitid != qh visit_id) {
2002         if (qh DOintersections)
2003           qh_printhyperplaneintersection(fp, facet, neighbor, ridge->vertices, black);
2004         if (qh PRINTridges) {
2005           vertexA= SETfirstt_(ridge->vertices, vertexT);
2006           vertexB= SETsecondt_(ridge->vertices, vertexT);
2007           qh_printline3geom (fp, vertexA->point, vertexB->point, green);
2008         }
2009       }
2010     }
2011   }
2012 } /* printfacet3geom_nonsimplicial */
2013
2014 /*-<a                             href="qh-io.htm#TOC"
2015   >-------------------------------</a><a name="printfacet3geom_points">-</a>
2016   
2017   qh_printfacet3geom_points( fp, points, facet, offset )
2018     prints a 3-d facet as OFF Geomview object. 
2019     offset is relative to the facet's hyperplane
2020     Facet is determined as a list of points
2021 */
2022 void qh_printfacet3geom_points(FILE *fp, setT *points, facetT *facet, realT offset, realT color[3]) {
2023   int k, n= qh_setsize(points), i;
2024   pointT *point, **pointp;
2025   setT *printpoints;
2026
2027   fprintf(fp, "{ OFF %d 1 1 # f%d\n", n, facet->id);
2028   if (offset != 0.0) {
2029     printpoints= qh_settemp (n);
2030     FOREACHpoint_(points) 
2031       qh_setappend (&printpoints, qh_projectpoint(point, facet, -offset));
2032   }else
2033     printpoints= points;
2034   FOREACHpoint_(printpoints) {
2035     for (k=0; k < qh hull_dim; k++) {
2036       if (k == qh DROPdim)
2037         fprintf(fp, "0 ");
2038       else
2039         fprintf(fp, "%8.4g ", point[k]);
2040     }
2041     if (printpoints != points)
2042       qh_memfree (point, qh normal_size);
2043     fprintf (fp, "\n");
2044   }
2045   if (printpoints != points)
2046     qh_settempfree (&printpoints);
2047   fprintf(fp, "%d ", n);
2048   for(i= 0; i < n; i++)
2049     fprintf(fp, "%d ", i);
2050   fprintf(fp, "%8.4g %8.4g %8.4g 1.0 }\n", color[0], color[1], color[2]);
2051 } /* printfacet3geom_points */
2052
2053
2054 /*-<a                             href="qh-io.htm#TOC"
2055   >-------------------------------</a><a name="printfacet3geom_simplicial">-</a>
2056   
2057   qh_printfacet3geom_simplicial(  )
2058     print Geomview OFF for a 3-d simplicial facet.
2059
2060   notes:
2061     may flip color
2062     uses facet->visitid for intersections and ridges
2063
2064     assume precise calculations in io.c with roundoff covered by qh_GEOMepsilon
2065     innerplane may be off by qh DISTround.  Maxoutside is calculated elsewhere
2066     so a DISTround error may have occured.
2067 */
2068 void qh_printfacet3geom_simplicial(FILE *fp, facetT *facet, realT color[3]) {
2069   setT *points, *vertices;
2070   vertexT *vertex, **vertexp, *vertexA, *vertexB;
2071   facetT *neighbor, **neighborp;
2072   realT outerplane, innerplane;
2073   realT black[3]={0, 0, 0}, green[3]={0, 1, 0};
2074   int k;
2075
2076   qh_geomplanes (facet, &outerplane, &innerplane); 
2077   vertices= qh_facet3vertex (facet);
2078   points= qh_settemp (qh TEMPsize);
2079   FOREACHvertex_(vertices)
2080     qh_setappend(&points, vertex->point);
2081   if (qh PRINTouter || (!qh PRINTnoplanes && !qh PRINTinner))
2082     qh_printfacet3geom_points(fp, points, facet, outerplane, color);
2083   if (qh PRINTinner || (!qh PRINTnoplanes && !qh PRINTouter &&
2084               outerplane - innerplane > 2 * qh MAXabs_coord * qh_GEOMepsilon)) {
2085     for (k= 3; k--; )
2086       color[k]= 1.0 - color[k];
2087     qh_printfacet3geom_points(fp, points, facet, innerplane, color);
2088   }
2089   qh_settempfree(&points);
2090   qh_settempfree(&vertices);
2091   if ((qh DOintersections || qh PRINTridges)
2092   && (!facet->visible || !qh NEWfacets)) {
2093     facet->visitid= qh visit_id;
2094     FOREACHneighbor_(facet) {
2095       if (neighbor->visitid != qh visit_id) {
2096         vertices= qh_setnew_delnthsorted (facet->vertices, qh hull_dim,
2097                           SETindex_(facet->neighbors, neighbor), 0);
2098         if (qh DOintersections)
2099            qh_printhyperplaneintersection(fp, facet, neighbor, vertices, black); 
2100         if (qh PRINTridges) {
2101           vertexA= SETfirstt_(vertices, vertexT);
2102           vertexB= SETsecondt_(vertices, vertexT);
2103           qh_printline3geom (fp, vertexA->point, vertexB->point, green);
2104         }
2105         qh_setfree(&vertices);
2106       }
2107     }
2108   }
2109 } /* printfacet3geom_simplicial */
2110
2111 /*-<a                             href="qh-io.htm#TOC"
2112   >-------------------------------</a><a name="printfacet3math">-</a>
2113   
2114   qh_printfacet3math( fp, facet, notfirst )
2115     print 3-d Mathematica output for a facet
2116
2117   notes:
2118     may be non-simplicial
2119     use %16.8f since Mathematica 2.2 does not handle exponential format
2120 */
2121 void qh_printfacet3math (FILE *fp, facetT *facet, int notfirst) {
2122   vertexT *vertex, **vertexp;
2123   setT *points, *vertices;
2124   pointT *point, **pointp;
2125   boolT firstpoint= True;
2126   realT dist;
2127   
2128   if (notfirst)
2129     fprintf(fp, ",\n");
2130   vertices= qh_facet3vertex (facet);
2131   points= qh_settemp (qh_setsize (vertices));
2132   FOREACHvertex_(vertices) {
2133     zinc_(Zdistio);
2134     qh_distplane(vertex->point, facet, &dist);
2135     point= qh_projectpoint(vertex->point, facet, dist);
2136     qh_setappend (&points, point);
2137   }
2138   fprintf(fp, "Polygon[{");
2139   FOREACHpoint_(points) {
2140     if (firstpoint)
2141       firstpoint= False;
2142     else
2143       fprintf(fp, ",\n");
2144     fprintf(fp, "{%16.8f, %16.8f, %16.8f}", point[0], point[1], point[2]);
2145   }
2146   FOREACHpoint_(points)
2147     qh_memfree (point, qh normal_size);
2148   qh_settempfree(&points);
2149   qh_settempfree(&vertices);
2150   fprintf(fp, "}]");
2151 } /* printfacet3math */
2152
2153
2154 /*-<a                             href="qh-io.htm#TOC"
2155   >-------------------------------</a><a name="printfacet3vertex">-</a>
2156   
2157   qh_printfacet3vertex( fp, facet, format )
2158     print vertices in a 3-d facet as point ids
2159
2160   notes:
2161     prints number of vertices first if format == qh_PRINToff
2162     the facet may be non-simplicial
2163 */
2164 void qh_printfacet3vertex(FILE *fp, facetT *facet, int format) {
2165   vertexT *vertex, **vertexp;
2166   setT *vertices;
2167
2168   vertices= qh_facet3vertex (facet);
2169   if (format == qh_PRINToff)
2170     fprintf (fp, "%d ", qh_setsize (vertices));
2171   FOREACHvertex_(vertices) 
2172     fprintf (fp, "%d ", qh_pointid(vertex->point));
2173   fprintf (fp, "\n");
2174   qh_settempfree(&vertices);
2175 } /* printfacet3vertex */
2176
2177
2178 /*-<a                             href="qh-io.htm#TOC"
2179   >-------------------------------</a><a name="printfacet4geom_nonsimplicial">-</a>
2180   
2181   qh_printfacet4geom_nonsimplicial(  )
2182     print Geomview 4OFF file for a 4d nonsimplicial facet
2183     prints all ridges to unvisited neighbors (qh.visit_id)
2184     if qh.DROPdim
2185       prints in OFF format
2186   
2187   notes:
2188     must agree with printend4geom()
2189 */
2190 void qh_printfacet4geom_nonsimplicial(FILE *fp, facetT *facet, realT color[3]) {
2191   facetT *neighbor;
2192   ridgeT *ridge, **ridgep;
2193   vertexT *vertex, **vertexp;
2194   pointT *point;
2195   int k;
2196   realT dist;
2197   
2198   facet->visitid= qh visit_id;
2199   if (qh PRINTnoplanes || (facet->visible && qh NEWfacets))
2200     return;
2201   FOREACHridge_(facet->ridges) {
2202     neighbor= otherfacet_(ridge, facet);
2203     if (neighbor->visitid == qh visit_id) 
2204       continue;
2205     if (qh PRINTtransparent && !neighbor->good)
2206       continue;  
2207     if (qh DOintersections)
2208       qh_printhyperplaneintersection(fp, facet, neighbor, ridge->vertices, color);
2209     else {
2210       if (qh DROPdim >= 0) 
2211         fprintf(fp, "OFF 3 1 1 # f%d\n", facet->id);
2212       else {
2213         qh printoutvar++;
2214         fprintf (fp, "# r%d between f%d f%d\n", ridge->id, facet->id, neighbor->id);
2215       }
2216       FOREACHvertex_(ridge->vertices) {
2217         zinc_(Zdistio);
2218         qh_distplane(vertex->point,facet, &dist);
2219         point=qh_projectpoint(vertex->point,facet, dist);
2220         for(k= 0; k < qh hull_dim; k++) {
2221           if (k != qh DROPdim)
2222             fprintf(fp, "%8.4g ", point[k]);
2223         }
2224         fprintf (fp, "\n");
2225         qh_memfree (point, qh normal_size);
2226       }
2227       if (qh DROPdim >= 0)
2228         fprintf(fp, "3 0 1 2 %8.4g %8.4g %8.4g\n", color[0], color[1], color[2]);
2229     }
2230   }
2231 } /* printfacet4geom_nonsimplicial */
2232
2233
2234 /*-<a                             href="qh-io.htm#TOC"
2235   >-------------------------------</a><a name="printfacet4geom_simplicial">-</a>
2236   
2237   qh_printfacet4geom_simplicial( fp, facet, color )
2238     print Geomview 4OFF file for a 4d simplicial facet
2239     prints triangles for unvisited neighbors (qh.visit_id)
2240
2241   notes:
2242     must agree with printend4geom()
2243 */
2244 void qh_printfacet4geom_simplicial(FILE *fp, facetT *facet, realT color[3]) {
2245   setT *vertices;
2246   facetT *neighbor, **neighborp;
2247   vertexT *vertex, **vertexp;
2248   int k;
2249   
2250   facet->visitid= qh visit_id;
2251   if (qh PRINTnoplanes || (facet->visible && qh NEWfacets))
2252     return;
2253   FOREACHneighbor_(facet) {
2254     if (neighbor->visitid == qh visit_id)
2255       continue;
2256     if (qh PRINTtransparent && !neighbor->good)
2257       continue;  
2258     vertices= qh_setnew_delnthsorted (facet->vertices, qh hull_dim,
2259                           SETindex_(facet->neighbors, neighbor), 0);
2260     if (qh DOintersections)
2261       qh_printhyperplaneintersection(fp, facet, neighbor, vertices, color);
2262     else {
2263       if (qh DROPdim >= 0) 
2264         fprintf(fp, "OFF 3 1 1 # ridge between f%d f%d\n",
2265                 facet->id, neighbor->id);
2266       else {
2267         qh printoutvar++;
2268         fprintf (fp, "# ridge between f%d f%d\n", facet->id, neighbor->id);
2269       }
2270       FOREACHvertex_(vertices) {
2271         for(k= 0; k < qh hull_dim; k++) {
2272           if (k != qh DROPdim)
2273             fprintf(fp, "%8.4g ", vertex->point[k]);
2274         }
2275         fprintf (fp, "\n");
2276       }
2277       if (qh DROPdim >= 0) 
2278         fprintf(fp, "3 0 1 2 %8.4g %8.4g %8.4g\n", color[0], color[1], color[2]);
2279     }
2280     qh_setfree(&vertices);
2281   }
2282 } /* printfacet4geom_simplicial */
2283
2284
2285 /*-<a                             href="qh-io.htm#TOC"
2286   >-------------------------------</a><a name="printfacetNvertex_nonsimplicial">-</a>
2287   
2288   qh_printfacetNvertex_nonsimplicial( fp, facet, id, format )
2289     print vertices for an N-d non-simplicial facet
2290     triangulates each ridge to the id
2291 */
2292 void qh_printfacetNvertex_nonsimplicial(FILE *fp, facetT *facet, int id, int format) {
2293   vertexT *vertex, **vertexp;
2294   ridgeT *ridge, **ridgep;
2295
2296   if (facet->visible && qh NEWfacets)
2297     return;
2298   FOREACHridge_(facet->ridges) {
2299     if (format == qh_PRINTtriangles)
2300       fprintf(fp, "%d ", qh hull_dim);
2301     fprintf(fp, "%d ", id);
2302     if ((ridge->top == facet) ^ qh_ORIENTclock) {
2303       FOREACHvertex_(ridge->vertices)
2304         fprintf(fp, "%d ", qh_pointid(vertex->point));
2305     }else {
2306       FOREACHvertexreverse12_(ridge->vertices)
2307         fprintf(fp, "%d ", qh_pointid(vertex->point));
2308     }
2309     fprintf(fp, "\n");
2310   }
2311 } /* printfacetNvertex_nonsimplicial */
2312
2313
2314 /*-<a                             href="qh-io.htm#TOC"
2315   >-------------------------------</a><a name="printfacetNvertex_simplicial">-</a>
2316   
2317   qh_printfacetNvertex_simplicial( fp, facet, format )
2318     print vertices for an N-d simplicial facet
2319     prints vertices for non-simplicial facets
2320       2-d facets (orientation preserved by qh_mergefacet2d)
2321       PRINToff ('o') for 4-d and higher
2322 */
2323 void qh_printfacetNvertex_simplicial(FILE *fp, facetT *facet, int format) {
2324   vertexT *vertex, **vertexp;
2325
2326   if (format == qh_PRINToff || format == qh_PRINTtriangles)
2327     fprintf (fp, "%d ", qh_setsize (facet->vertices));
2328   if ((facet->toporient ^ qh_ORIENTclock) 
2329   || (qh hull_dim > 2 && !facet->simplicial)) {
2330     FOREACHvertex_(facet->vertices)
2331       fprintf(fp, "%d ", qh_pointid(vertex->point));
2332   }else {
2333     FOREACHvertexreverse12_(facet->vertices)
2334       fprintf(fp, "%d ", qh_pointid(vertex->point));
2335   }
2336   fprintf(fp, "\n");
2337 } /* printfacetNvertex_simplicial */
2338
2339
2340 /*-<a                             href="qh-io.htm#TOC"
2341   >-------------------------------</a><a name="printfacetheader">-</a>
2342   
2343   qh_printfacetheader( fp, facet )
2344     prints header fields of a facet to fp
2345     
2346   notes:
2347     for 'f' output and debugging
2348 */
2349 void qh_printfacetheader(FILE *fp, facetT *facet) {
2350   pointT *point, **pointp, *furthest;
2351   facetT *neighbor, **neighborp;
2352   realT dist;
2353
2354   if (facet == qh_MERGEridge) {
2355     fprintf (fp, " MERGEridge\n");
2356     return;
2357   }else if (facet == qh_DUPLICATEridge) {
2358     fprintf (fp, " DUPLICATEridge\n");
2359     return;
2360   }else if (!facet) {
2361     fprintf (fp, " NULLfacet\n");
2362     return;
2363   }
2364   qh old_randomdist= qh RANDOMdist;
2365   qh RANDOMdist= False;
2366   fprintf(fp, "- f%d\n", facet->id);
2367   fprintf(fp, "    - flags:");
2368   if (facet->toporient) 
2369     fprintf(fp, " top");
2370   else
2371     fprintf(fp, " bottom");
2372   if (facet->simplicial)
2373     fprintf(fp, " simplicial");
2374   if (facet->tricoplanar)
2375     fprintf(fp, " tricoplanar");
2376   if (facet->upperdelaunay)
2377     fprintf(fp, " upperDelaunay");
2378   if (facet->visible)
2379     fprintf(fp, " visible");
2380   if (facet->newfacet)
2381     fprintf(fp, " new");
2382   if (facet->tested)
2383     fprintf(fp, " tested");
2384   if (!facet->good)
2385     fprintf(fp, " notG");
2386   if (facet->seen)
2387     fprintf(fp, " seen");
2388   if (facet->coplanar)
2389     fprintf(fp, " coplanar");
2390   if (facet->mergehorizon)
2391     fprintf(fp, " mergehorizon");
2392   if (facet->keepcentrum)
2393     fprintf(fp, " keepcentrum");
2394   if (facet->dupridge)
2395     fprintf(fp, " dupridge");
2396   if (facet->mergeridge && !facet->mergeridge2)
2397     fprintf(fp, " mergeridge1");
2398   if (facet->mergeridge2)
2399     fprintf(fp, " mergeridge2");
2400   if (facet->newmerge)
2401     fprintf(fp, " newmerge");
2402   if (facet->flipped) 
2403     fprintf(fp, " flipped");
2404   if (facet->notfurthest) 
2405     fprintf(fp, " notfurthest");
2406   if (facet->degenerate)
2407     fprintf(fp, " degenerate");
2408   if (facet->redundant)
2409     fprintf(fp, " redundant");
2410   fprintf(fp, "\n");
2411   if (facet->isarea)
2412     fprintf(fp, "    - area: %2.2g\n", facet->f.area);
2413   else if (qh NEWfacets && facet->visible && facet->f.replace)
2414     fprintf(fp, "    - replacement: f%d\n", facet->f.replace->id);
2415   else if (facet->newfacet) {
2416     if (facet->f.samecycle && facet->f.samecycle != facet)
2417       fprintf(fp, "    - shares same visible/horizon as f%d\n", facet->f.samecycle->id);
2418   }else if (facet->tricoplanar /* !isarea */) {
2419     if (facet->f.triowner)
2420       fprintf(fp, "    - owner of normal & centrum is facet f%d\n", facet->f.triowner->id);
2421   }else if (facet->f.newcycle)
2422     fprintf(fp, "    - was horizon to f%d\n", facet->f.newcycle->id);
2423   if (facet->nummerge)
2424     fprintf(fp, "    - merges: %d\n", facet->nummerge);
2425   qh_printpointid(fp, "    - normal: ", qh hull_dim, facet->normal, -1);
2426   fprintf(fp, "    - offset: %10.7g\n", facet->offset);
2427   if (qh CENTERtype == qh_ASvoronoi || facet->center)
2428     qh_printcenter (fp, qh_PRINTfacets, "    - center: ", facet);
2429 #if qh_MAXoutside
2430   if (facet->maxoutside > qh DISTround)
2431     fprintf(fp, "    - maxoutside: %10.7g\n", facet->maxoutside);
2432 #endif
2433   if (!SETempty_(facet->outsideset)) {
2434     furthest= (pointT*)qh_setlast(facet->outsideset);
2435     if (qh_setsize (facet->outsideset) < 6) {
2436       fprintf(fp, "    - outside set (furthest p%d):\n", qh_pointid(furthest));
2437       FOREACHpoint_(facet->outsideset)
2438         qh_printpoint(fp, "     ", point);
2439     }else if (qh_setsize (facet->outsideset) < 21) {
2440       qh_printpoints(fp, "    - outside set:", facet->outsideset);
2441     }else {
2442       fprintf(fp, "    - outside set:  %d points.", qh_setsize(facet->outsideset));
2443       qh_printpoint(fp, "  Furthest", furthest);
2444     }
2445 #if !qh_COMPUTEfurthest
2446     fprintf(fp, "    - furthest distance= %2.2g\n", facet->furthestdist);
2447 #endif
2448   }
2449   if (!SETempty_(facet->coplanarset)) {
2450     furthest= (pointT*)qh_setlast(facet->coplanarset);
2451     if (qh_setsize (facet->coplanarset) < 6) {
2452       fprintf(fp, "    - coplanar set (furthest p%d):\n", qh_pointid(furthest));
2453       FOREACHpoint_(facet->coplanarset)
2454         qh_printpoint(fp, "     ", point);
2455     }else if (qh_setsize (facet->coplanarset) < 21) {
2456       qh_printpoints(fp, "    - coplanar set:", facet->coplanarset);
2457     }else {
2458       fprintf(fp, "    - coplanar set:  %d points.", qh_setsize(facet->coplanarset));
2459       qh_printpoint(fp, "  Furthest", furthest);
2460     }
2461     zinc_(Zdistio);
2462     qh_distplane (furthest, facet, &dist);
2463     fprintf(fp, "      furthest distance= %2.2g\n", dist);
2464   }
2465   qh_printvertices (fp, "    - vertices:", facet->vertices);
2466   fprintf(fp, "    - neighboring facets: ");
2467   FOREACHneighbor_(facet) {
2468     if (neighbor == qh_MERGEridge)
2469       fprintf(fp, " MERGE");
2470     else if (neighbor == qh_DUPLICATEridge)
2471       fprintf(fp, " DUP");
2472     else
2473       fprintf(fp, " f%d", neighbor->id);
2474   }
2475   fprintf(fp, "\n");
2476   qh RANDOMdist= qh old_randomdist;
2477 } /* printfacetheader */
2478
2479
2480 /*-<a                             href="qh-io.htm#TOC"
2481   >-------------------------------</a><a name="printfacetridges">-</a>
2482   
2483   qh_printfacetridges( fp, facet )
2484     prints ridges of a facet to fp
2485
2486   notes:
2487     ridges printed in neighbor order
2488     assumes the ridges exist
2489     for 'f' output
2490 */
2491 void qh_printfacetridges(FILE *fp, facetT *facet) {
2492   facetT *neighbor, **neighborp;
2493   ridgeT *ridge, **ridgep;
2494   int numridges= 0;
2495
2496
2497   if (facet->visible && qh NEWfacets) {
2498     fprintf(fp, "    - ridges (ids may be garbage):");
2499     FOREACHridge_(facet->ridges)
2500       fprintf(fp, " r%d", ridge->id);
2501     fprintf(fp, "\n");
2502   }else {
2503     fprintf(fp, "    - ridges:\n");
2504     FOREACHridge_(facet->ridges)
2505       ridge->seen= False;
2506     if (qh hull_dim == 3) {
2507       ridge= SETfirstt_(facet->ridges, ridgeT);
2508       while (ridge && !ridge->seen) {
2509         ridge->seen= True;
2510         qh_printridge(fp, ridge);
2511         numridges++;
2512         ridge= qh_nextridge3d (ridge, facet, NULL);
2513         }
2514     }else {
2515       FOREACHneighbor_(facet) {
2516         FOREACHridge_(facet->ridges) {
2517           if (otherfacet_(ridge,facet) == neighbor) {
2518             ridge->seen= True;
2519             qh_printridge(fp, ridge);
2520             numridges++;
2521           }
2522         }
2523       }
2524     }
2525     if (numridges != qh_setsize (facet->ridges)) {
2526       fprintf (fp, "     - all ridges:");
2527       FOREACHridge_(facet->ridges) 
2528         fprintf (fp, " r%d", ridge->id);
2529         fprintf (fp, "\n");
2530     }
2531     FOREACHridge_(facet->ridges) {
2532       if (!ridge->seen) 
2533         qh_printridge(fp, ridge);
2534     }
2535   }
2536 } /* printfacetridges */
2537
2538 /*-<a                             href="qh-io.htm#TOC"
2539   >-------------------------------</a><a name="printfacets">-</a>
2540   
2541   qh_printfacets( fp, format, facetlist, facets, printall )
2542     prints facetlist and/or facet set in output format
2543   
2544   notes:
2545     also used for specialized formats ('FO' and summary)
2546     turns off 'Rn' option since want actual numbers
2547 */
2548 void qh_printfacets(FILE *fp, int format, facetT *facetlist, setT *facets, boolT printall) {
2549   int numfacets, numsimplicial, numridges, totneighbors, numcoplanars, numtricoplanars;
2550   facetT *facet, **facetp;
2551   setT *vertices;
2552   coordT *center;
2553   realT outerplane, innerplane;
2554
2555   qh old_randomdist= qh RANDOMdist;
2556   qh RANDOMdist= False;
2557   if (qh CDDoutput && (format == qh_PRINTcentrums || format == qh_PRINTpointintersect || format == qh_PRINToff))
2558     fprintf (qh ferr, "qhull warning: CDD format is not available for centrums, halfspace\nintersections, and OFF file format.\n");
2559   if (format == qh_PRINTnone)
2560     ; /* print nothing */
2561   else if (format == qh_PRINTaverage) {
2562     vertices= qh_facetvertices (facetlist, facets, printall);
2563     center= qh_getcenter (vertices);
2564     fprintf (fp, "%d 1\n", qh hull_dim);
2565     qh_printpointid (fp, NULL, qh hull_dim, center, -1);
2566     qh_memfree (center, qh normal_size);
2567     qh_settempfree (&vertices);
2568   }else if (format == qh_PRINTextremes) {
2569     if (qh DELAUNAY)
2570       qh_printextremes_d (fp, facetlist, facets, printall);
2571     else if (qh hull_dim == 2)
2572       qh_printextremes_2d (fp, facetlist, facets, printall);
2573     else 
2574       qh_printextremes (fp, facetlist, facets, printall);
2575   }else if (format == qh_PRINToptions)
2576     fprintf(fp, "Options selected for Qhull %s:\n%s\n", qh_VERSION, qh qhull_options);
2577   else if (format == qh_PRINTpoints && !qh VORONOI)
2578     qh_printpoints_out (fp, facetlist, facets, printall);
2579   else if (format == qh_PRINTqhull)
2580     fprintf (fp, "%s | %s\n", qh rbox_command, qh qhull_command);
2581   else if (format == qh_PRINTsize) {
2582     fprintf (fp, "0\n2 ");
2583     fprintf (fp, qh_REAL_1, qh totarea);
2584     fprintf (fp, qh_REAL_1, qh totvol);
2585     fprintf (fp, "\n");
2586   }else if (format == qh_PRINTsummary) {
2587     qh_countfacets (facetlist, facets, printall, &numfacets, &numsimplicial, 
2588       &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
2589     vertices= qh_facetvertices (facetlist, facets, printall); 
2590     fprintf (fp, "10 %d %d %d %d %d %d %d %d %d %d\n2 ", qh hull_dim, 
2591                 qh num_points + qh_setsize (qh other_points),
2592                 qh num_vertices, qh num_facets - qh num_visible,
2593                 qh_setsize (vertices), numfacets, numcoplanars, 
2594                 numfacets - numsimplicial, zzval_(Zdelvertextot), 
2595                 numtricoplanars);
2596     qh_settempfree (&vertices);
2597     qh_outerinner (NULL, &outerplane, &innerplane);
2598     fprintf (fp, qh_REAL_2n, outerplane, innerplane);
2599   }else if (format == qh_PRINTvneighbors)
2600     qh_printvneighbors (fp, facetlist, facets, printall);
2601   else if (qh VORONOI && format == qh_PRINToff)
2602     qh_printvoronoi (fp, format, facetlist, facets, printall);
2603   else if (qh VORONOI && format == qh_PRINTgeom) {
2604     qh_printbegin (fp, format, facetlist, facets, printall);
2605     qh_printvoronoi (fp, format, facetlist, facets, printall);
2606     qh_printend (fp, format, facetlist, facets, printall);
2607   }else if (qh VORONOI 
2608   && (format == qh_PRINTvertices || format == qh_PRINTinner || format == qh_PRINTouter))
2609     qh_printvdiagram (fp, format, facetlist, facets, printall);
2610   else {
2611     qh_printbegin (fp, format, facetlist, facets, printall);
2612     FORALLfacet_(facetlist)
2613       qh_printafacet (fp, format, facet, printall);
2614     FOREACHfacet_(facets) 
2615       qh_printafacet (fp, format, facet, printall);
2616     qh_printend (fp, format, facetlist, facets, printall);
2617   }
2618   qh RANDOMdist= qh old_randomdist;
2619 } /* printfacets */
2620
2621
2622 /*-<a                             href="qh-io.htm#TOC"
2623   >-------------------------------</a><a name="printhelp_degenerate">-</a>
2624   
2625   qh_printhelp_degenerate( fp )
2626     prints descriptive message for precision error
2627
2628   notes:
2629     no message if qh_QUICKhelp
2630 */
2631 void qh_printhelp_degenerate(FILE *fp) {
2632   
2633   if (qh MERGEexact || qh PREmerge || qh JOGGLEmax < REALmax/2) 
2634     fprintf(fp, "\n\
2635 A Qhull error has occurred.  Qhull should have corrected the above\n\
2636 precision error.  Please send the input and all of the output to\n\
2637 qhull_bug@geom.umn.edu\n");
2638   else if (!qh_QUICKhelp) {
2639     fprintf(fp, "\n\
2640 Precision problems were detected during construction of the convex hull.\n\
2641 This occurs because convex hull algorithms assume that calculations are\n\
2642 exact, but floating-point arithmetic has roundoff errors.\n\
2643 \n\
2644 To correct for precision problems, do not use 'Q0'.  By default, Qhull\n\
2645 selects 'C-0' or 'Qx' and merges non-convex facets.  With option 'QJ',\n\
2646 Qhull joggles the input to prevent precision problems.  See \"Imprecision\n\
2647 in Qhull\" (qh-impre.htm).\n\
2648 \n\
2649 If you use 'Q0', the output may include\n\
2650 coplanar ridges, concave ridges, and flipped facets.  In 4-d and higher,\n\
2651 Qhull may produce a ridge with four neighbors or two facets with the same \n\
2652 vertices.  Qhull reports these events when they occur.  It stops when a\n\
2653 concave ridge, flipped facet, or duplicate facet occurs.\n");
2654 #if REALfloat
2655     fprintf (fp, "\
2656 \n\
2657 Qhull is currently using single precision arithmetic.  The following\n\
2658 will probably remove the precision problems:\n\
2659   - recompile qhull for double precision (#define REALfloat 0 in user.h).\n");
2660 #endif
2661     if (qh DELAUNAY && !qh SCALElast && qh MAXabs_coord > 1e4)
2662       fprintf( fp, "\
2663 \n\
2664 When computing the Delaunay triangulation of coordinates > 1.0,\n\
2665   - use 'Qbb' to scale the last coordinate to [0,m] (max previous coordinate)\n");
2666     if (qh DELAUNAY && !qh ATinfinity) 
2667       fprintf( fp, "\
2668 When computing the Delaunay triangulation:\n\
2669   - use 'Qz' to add a point at-infinity.  This reduces precision problems.\n");
2670  
2671     fprintf(fp, "\
2672 \n\
2673 If you need triangular output:\n\
2674   - use option 'Qt' to triangulate the output\n\
2675   - use option 'QJ' to joggle the input points and remove precision errors\n\
2676   - use option 'Ft'.  It triangulates non-simplicial facets with added points.\n\
2677 \n\
2678 If you must use 'Q0',\n\
2679 try one or more of the following options.  They can not guarantee an output.\n\
2680   - use 'QbB' to scale the input to a cube.\n\
2681   - use 'Po' to produce output and prevent partitioning for flipped facets\n\
2682   - use 'V0' to set min. distance to visible facet as 0 instead of roundoff\n\
2683   - use 'En' to specify a maximum roundoff error less than %2.2g.\n\
2684   - options 'Qf', 'Qbb', and 'QR0' may also help\n",
2685                qh DISTround);
2686     fprintf(fp, "\
2687 \n\
2688 To guarantee simplicial output:\n\
2689   - use option 'Qt' to triangulate the output\n\
2690   - use option 'QJ' to joggle the input points and remove precision errors\n\
2691   - use option 'Ft' to triangulate the output by adding points\n\
2692   - use exact arithmetic (see \"Imprecision in Qhull\", qh-impre.htm)\n\
2693 ");
2694   }
2695 } /* printhelp_degenerate */
2696
2697
2698 /*-<a                             href="qh-io.htm#TOC"
2699   >-------------------------------</a><a name="printhelp_singular">-</a>
2700   
2701   qh_printhelp_singular( fp )
2702     prints descriptive message for singular input
2703 */
2704 void qh_printhelp_singular(FILE *fp) {
2705   facetT *facet;
2706   vertexT *vertex, **vertexp;
2707   realT min, max, *coord, dist;
2708   int i,k;
2709   
2710   fprintf(fp, "\n\
2711 The input to qhull appears to be less than %d dimensional, or a\n\
2712 computation has overflowed.\n\n\
2713 Qhull could not construct a clearly convex simplex from points:\n",
2714            qh hull_dim);
2715   qh_printvertexlist (fp, "", qh facet_list, NULL, qh_ALL);
2716   if (!qh_QUICKhelp)
2717     fprintf(fp, "\n\
2718 The center point is coplanar with a facet, or a vertex is coplanar\n\
2719 with a neighboring facet.  The maximum round off error for\n\
2720 computing distances is %2.2g.  The center point, facets and distances\n\
2721 to the center point are as follows:\n\n", qh DISTround);
2722   qh_printpointid (fp, "center point", qh hull_dim, qh interior_point, -1);
2723   fprintf (fp, "\n");
2724   FORALLfacets {
2725     fprintf (fp, "facet");
2726     FOREACHvertex_(facet->vertices)
2727       fprintf (fp, " p%d", qh_pointid(vertex->point));
2728     zinc_(Zdistio);
2729     qh_distplane(qh interior_point, facet, &dist);
2730     fprintf (fp, " distance= %4.2g\n", dist);
2731   }
2732   if (!qh_QUICKhelp) {
2733     if (qh HALFspace) 
2734       fprintf (fp, "\n\
2735 These points are the dual of the given halfspaces.  They indicate that\n\
2736 the intersection is degenerate.\n");
2737     fprintf (fp,"\n\
2738 These points either have a maximum or minimum x-coordinate, or\n\
2739 they maximize the determinant for k coordinates.  Trial points\n\
2740 are first selected from points that maximize a coordinate.\n");
2741     if (qh hull_dim >= qh_INITIALmax)
2742       fprintf (fp, "\n\
2743 Because of the high dimension, the min x-coordinate and max-coordinate\n\
2744 points are used if the determinant is non-zero.  Option 'Qs' will\n\
2745 do a better, though much slower, job.  Instead of 'Qs', you can change\n\
2746 the points by randomly rotating the input with 'QR0'.\n");
2747   }
2748   fprintf (fp, "\nThe min and max coordinates for each dimension are:\n");
2749   for (k=0; k < qh hull_dim; k++) {
2750     min= REALmax;
2751     max= -REALmin;
2752     for (i=qh num_points, coord= qh first_point+k; i--; coord += qh hull_dim) {
2753       maximize_(max, *coord);
2754       minimize_(min, *coord);
2755     }
2756     fprintf (fp, "  %d:  %8.4g  %8.4g  difference= %4.4g\n", k, min, max, max-min);
2757   }
2758   if (!qh_QUICKhelp) {
2759     fprintf (fp, "\n\
2760 If the input should be full dimensional, you have several options that\n\
2761 may determine an initial simplex:\n\
2762   - use 'QJ'  to joggle the input and make it full dimensional\n\
2763   - use 'QbB' to scale the points to the unit cube\n\
2764   - use 'QR0' to randomly rotate the input for different maximum points\n\
2765   - use 'Qs'  to search all points for the initial simplex\n\
2766   - use 'En'  to specify a maximum roundoff error less than %2.2g.\n\
2767   - trace execution with 'T3' to see the determinant for each point.\n",
2768                      qh DISTround);
2769 #if REALfloat
2770     fprintf (fp, "\
2771   - recompile qhull for double precision (#define REALfloat 0 in qhull.h).\n");
2772 #endif
2773     fprintf (fp, "\n\
2774 If the input is lower dimensional:\n\
2775   - use 'QJ' to joggle the input and make it full dimensional\n\
2776   - use 'Qbk:0Bk:0' to delete coordinate k from the input.  You should\n\
2777     pick the coordinate with the least range.  The hull will have the\n\
2778     correct topology.\n\
2779   - determine the flat containing the points, rotate the points\n\
2780     into a coordinate plane, and delete the other coordinates.\n\
2781   - add one or more points to make the input full dimensional.\n\
2782 ");
2783     if (qh DELAUNAY && !qh ATinfinity)
2784       fprintf (fp, "\n\n\
2785 This is a Delaunay triangulation and the input is co-circular or co-spherical:\n\
2786   - use 'Qz' to add a point \"at infinity\" (i.e., above the paraboloid)\n\
2787   - or use 'QJ' to joggle the input and avoid co-circular data\n");
2788   }
2789 } /* printhelp_singular */
2790
2791 /*-<a                             href="qh-io.htm#TOC"
2792   >-------------------------------</a><a name="printhyperplaneintersection">-</a>
2793   
2794   qh_printhyperplaneintersection( fp, facet1, facet2, vertices, color )
2795     print Geomview OFF or 4OFF for the intersection of two hyperplanes in 3-d or 4-d
2796 */
2797 void qh_printhyperplaneintersection(FILE *fp, facetT *facet1, facetT *facet2,
2798                    setT *vertices, realT color[3]) {
2799   realT costheta, denominator, dist1, dist2, s, t, mindenom, p[4];
2800   vertexT *vertex, **vertexp;
2801   int i, k;
2802   boolT nearzero1, nearzero2;
2803   
2804   costheta= qh_getangle(facet1->normal, facet2->normal);
2805   denominator= 1 - costheta * costheta;
2806   i= qh_setsize(vertices);
2807   if (qh hull_dim == 3)
2808     fprintf(fp, "VECT 1 %d 1 %d 1 ", i, i);
2809   else if (qh hull_dim == 4 && qh DROPdim >= 0)
2810     fprintf(fp, "OFF 3 1 1 ");
2811   else
2812     qh printoutvar++;
2813   fprintf (fp, "# intersect f%d f%d\n", facet1->id, facet2->id);
2814   mindenom= 1 / (10.0 * qh MAXabs_coord);
2815   FOREACHvertex_(vertices) {
2816     zadd_(Zdistio, 2);
2817     qh_distplane(vertex->point, facet1, &dist1);
2818     qh_distplane(vertex->point, facet2, &dist2);
2819     s= qh_divzero (-dist1 + costheta * dist2, denominator,mindenom,&nearzero1);
2820     t= qh_divzero (-dist2 + costheta * dist1, denominator,mindenom,&nearzero2);
2821     if (nearzero1 || nearzero2)
2822       s= t= 0.0;
2823     for(k= qh hull_dim; k--; )
2824       p[k]= vertex->point[k] + facet1->normal[k] * s + facet2->normal[k] * t;
2825     if (qh PRINTdim <= 3) {
2826       qh_projectdim3 (p, p);
2827       fprintf(fp, "%8.4g %8.4g %8.4g # ", p[0], p[1], p[2]);
2828     }else 
2829       fprintf(fp, "%8.4g %8.4g %8.4g %8.4g # ", p[0], p[1], p[2], p[3]);
2830     if (nearzero1+nearzero2)
2831       fprintf (fp, "p%d (coplanar facets)\n", qh_pointid (vertex->point));
2832     else
2833       fprintf (fp, "projected p%d\n", qh_pointid (vertex->point));
2834   }
2835   if (qh hull_dim == 3)
2836     fprintf(fp, "%8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]); 
2837   else if (qh hull_dim == 4 && qh DROPdim >= 0)  
2838     fprintf(fp, "3 0 1 2 %8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]);
2839 } /* printhyperplaneintersection */
2840
2841 /*-<a                             href="qh-io.htm#TOC"
2842   >-------------------------------</a><a name="printline3geom">-</a>
2843   
2844   qh_printline3geom( fp, pointA, pointB, color )
2845     prints a line as a VECT
2846     prints 0's for qh.DROPdim
2847   
2848   notes:
2849     if pointA == pointB, 
2850       it's a 1 point VECT
2851 */
2852 void qh_printline3geom (FILE *fp, pointT *pointA, pointT *pointB, realT color[3]) {
2853   int k;
2854   realT pA[4], pB[4];
2855
2856   qh_projectdim3(pointA, pA);
2857   qh_projectdim3(pointB, pB);
2858   if ((fabs(pA[0] - pB[0]) > 1e-3) || 
2859       (fabs(pA[1] - pB[1]) > 1e-3) || 
2860       (fabs(pA[2] - pB[2]) > 1e-3)) {
2861     fprintf (fp, "VECT 1 2 1 2 1\n");
2862     for (k= 0; k < 3; k++)
2863        fprintf (fp, "%8.4g ", pB[k]);
2864     fprintf (fp, " # p%d\n", qh_pointid (pointB));
2865   }else
2866     fprintf (fp, "VECT 1 1 1 1 1\n");
2867   for (k=0; k < 3; k++)
2868     fprintf (fp, "%8.4g ", pA[k]);
2869   fprintf (fp, " # p%d\n", qh_pointid (pointA));
2870   fprintf (fp, "%8.4g %8.4g %8.4g 1\n", color[0], color[1], color[2]);
2871 }
2872
2873 /*-<a                             href="qh-io.htm#TOC"
2874   >-------------------------------</a><a name="printneighborhood">-</a>
2875   
2876   qh_printneighborhood( fp, format, facetA, facetB, printall )
2877     print neighborhood of one or two facets
2878
2879   notes:
2880     calls qh_findgood_all() 
2881     bumps qh.visit_id
2882 */
2883 void qh_printneighborhood (FILE *fp, int format, facetT *facetA, facetT *facetB, boolT printall) {
2884   facetT *neighbor, **neighborp, *facet;
2885   setT *facets;
2886
2887   if (format == qh_PRINTnone)
2888     return;
2889   qh_findgood_all (qh facet_list);
2890   if (facetA == facetB)
2891     facetB= NULL;
2892   facets= qh_settemp (2*(qh_setsize (facetA->neighbors)+1));
2893   qh visit_id++;
2894   for (facet= facetA; facet; facet= ((facet == facetA) ? facetB : NULL)) {
2895     if (facet->visitid != qh visit_id) {
2896       facet->visitid= qh visit_id;
2897       qh_setappend (&facets, facet);
2898     }
2899     FOREACHneighbor_(facet) {
2900       if (neighbor->visitid == qh visit_id)
2901         continue;
2902       neighbor->visitid= qh visit_id;
2903       if (printall || !qh_skipfacet (neighbor))
2904         qh_setappend (&facets, neighbor);
2905     }
2906   }
2907   qh_printfacets (fp, format, NULL, facets, printall);
2908   qh_settempfree (&facets);
2909 } /* printneighborhood */
2910
2911 /*-<a                             href="qh-io.htm#TOC"
2912   >-------------------------------</a><a name="printpoint">-</a>
2913   
2914   qh_printpoint( fp, string, point )
2915   qh_printpointid( fp, string, dim, point, id )
2916     prints the coordinates of a point
2917
2918   returns:
2919     if string is defined
2920       prints 'string p%d' (skips p%d if id=-1)
2921
2922   notes:
2923     nop if point is NULL
2924     prints id unless it is undefined (-1)
2925 */
2926 void qh_printpoint(FILE *fp, char *string, pointT *point) {
2927   int id= qh_pointid( point);
2928
2929   qh_printpointid( fp, string, qh hull_dim, point, id);
2930 } /* printpoint */
2931
2932 void qh_printpointid(FILE *fp, char *string, int dim, pointT *point, int id) {
2933   int k;
2934   realT r; /*bug fix*/
2935   
2936   if (!point)
2937     return;
2938   if (string) {
2939     fputs (string, fp);
2940    if (id != -1)
2941       fprintf(fp, " p%d: ", id);
2942   }
2943   for(k= dim; k--; ) {
2944     r= *point++;
2945     if (string)
2946       fprintf(fp, " %8.4g", r);
2947     else
2948       fprintf(fp, qh_REAL_1, r);
2949   }
2950   fprintf(fp, "\n");
2951 } /* printpointid */
2952
2953 /*-<a                             href="qh-io.htm#TOC"
2954   >-------------------------------</a><a name="printpoint3">-</a>
2955   
2956   qh_printpoint3( fp, point )
2957     prints 2-d, 3-d, or 4-d point as Geomview 3-d coordinates
2958 */
2959 void qh_printpoint3 (FILE *fp, pointT *point) {
2960   int k;
2961   realT p[4];
2962   
2963   qh_projectdim3 (point, p);
2964   for (k=0; k < 3; k++)
2965     fprintf (fp, "%8.4g ", p[k]);
2966   fprintf (fp, " # p%d\n", qh_pointid (point));
2967 } /* printpoint3 */
2968
2969 /*----------------------------------------
2970 -printpoints- print pointids for a set of points starting at index 
2971    see geom.c
2972 */
2973
2974 /*-<a                             href="qh-io.htm#TOC"
2975   >-------------------------------</a><a name="printpoints_out">-</a>
2976   
2977   qh_printpoints_out( fp, facetlist, facets, printall )
2978     prints vertices, coplanar/inside points, for facets by their point coordinates
2979     allows qh.CDDoutput
2980
2981   notes:
2982     same format as qhull input
2983     if no coplanar/interior points,
2984       same order as qh_printextremes
2985 */
2986 void qh_printpoints_out (FILE *fp, facetT *facetlist, setT *facets, int printall) {
2987   int allpoints= qh num_points + qh_setsize (qh other_points);
2988   int numpoints=0, point_i, point_n;
2989   setT *vertices, *points;
2990   facetT *facet, **facetp;
2991   pointT *point, **pointp;
2992   vertexT *vertex, **vertexp;
2993   int id;
2994
2995   points= qh_settemp (allpoints);
2996   qh_setzero (points, 0, allpoints);
2997   vertices= qh_facetvertices (facetlist, facets, printall);
2998   FOREACHvertex_(vertices) {
2999     id= qh_pointid (vertex->point);
3000     if (id >= 0)
3001       SETelem_(points, id)= vertex->point;
3002   }
3003   if (qh KEEPinside || qh KEEPcoplanar || qh KEEPnearinside) {
3004     FORALLfacet_(facetlist) {
3005       if (!printall && qh_skipfacet(facet))
3006         continue;
3007       FOREACHpoint_(facet->coplanarset) {
3008         id= qh_pointid (point);
3009         if (id >= 0)
3010           SETelem_(points, id)= point;
3011       }
3012     }
3013     FOREACHfacet_(facets) {
3014       if (!printall && qh_skipfacet(facet))
3015         continue;
3016       FOREACHpoint_(facet->coplanarset) {
3017         id= qh_pointid (point);
3018         if (id >= 0)
3019           SETelem_(points, id)= point;
3020       }
3021     }
3022   }
3023   qh_settempfree (&vertices);
3024   FOREACHpoint_i_(points) {
3025     if (point)
3026       numpoints++;
3027   }
3028   if (qh CDDoutput)
3029     fprintf (fp, "%s | %s\nbegin\n%d %d real\n", qh rbox_command,
3030              qh qhull_command, numpoints, qh hull_dim + 1);
3031   else
3032     fprintf (fp, "%d\n%d\n", qh hull_dim, numpoints);
3033   FOREACHpoint_i_(points) {
3034     if (point) {
3035       if (qh CDDoutput)
3036         fprintf (fp, "1 ");
3037       qh_printpoint (fp, NULL, point);
3038     }
3039   }
3040   if (qh CDDoutput)
3041     fprintf (fp, "end\n");
3042   qh_settempfree (&points);
3043 } /* printpoints_out */
3044   
3045
3046 /*-<a                             href="qh-io.htm#TOC"
3047   >-------------------------------</a><a name="printpointvect">-</a>
3048   
3049   qh_printpointvect( fp, point, normal, center, radius, color )
3050     prints a 2-d, 3-d, or 4-d point as 3-d VECT's relative to normal or to center point
3051 */
3052 void qh_printpointvect (FILE *fp, pointT *point, coordT *normal, pointT *center, realT radius, realT color[3]) {
3053   realT diff[4], pointA[4];
3054   int k;
3055   
3056   for (k= qh hull_dim; k--; ) {
3057     if (center)
3058       diff[k]= point[k]-center[k];
3059     else if (normal) 
3060       diff[k]= normal[k];
3061     else
3062       diff[k]= 0;
3063   }
3064   if (center)
3065     qh_normalize2 (diff, qh hull_dim, True, NULL, NULL);
3066   for (k= qh hull_dim; k--; ) 
3067     pointA[k]= point[k]+diff[k] * radius;
3068   qh_printline3geom (fp, point, pointA, color);
3069 } /* printpointvect */  
3070
3071 /*-<a                             href="qh-io.htm#TOC"
3072   >-------------------------------</a><a name="printpointvect2">-</a>
3073   
3074   qh_printpointvect2( fp, point, normal, center, radius )
3075     prints a 2-d, 3-d, or 4-d point as 2 3-d VECT's for an imprecise point
3076 */
3077 void qh_printpointvect2 (FILE *fp, pointT *point, coordT *normal, pointT *center, realT radius) {
3078   realT red[3]={1, 0, 0}, yellow[3]={1, 1, 0};
3079
3080   qh_printpointvect (fp, point, normal, center, radius, red);
3081   qh_printpointvect (fp, point, normal, center, -radius, yellow);
3082 } /* printpointvect2 */
3083
3084 /*-<a                             href="qh-io.htm#TOC"
3085   >-------------------------------</a><a name="printridge">-</a>
3086   
3087   qh_printridge( fp, ridge )
3088     prints the information in a ridge
3089
3090   notes:
3091     for qh_printfacetridges()
3092 */
3093 void qh_printridge(FILE *fp, ridgeT *ridge) {
3094   
3095   fprintf(fp, "     - r%d", ridge->id);
3096   if (ridge->tested)
3097     fprintf (fp, " tested");
3098   if (ridge->nonconvex)
3099     fprintf (fp, " nonconvex");
3100   fprintf (fp, "\n");
3101   qh_printvertices (fp, "           vertices:", ridge->vertices);
3102   if (ridge->top && ridge->bottom)
3103     fprintf(fp, "           between f%d and f%d\n",
3104             ridge->top->id, ridge->bottom->id);
3105 } /* printridge */
3106
3107 /*-<a                             href="qh-io.htm#TOC"
3108   >-------------------------------</a><a name="printspheres">-</a>
3109   
3110   qh_printspheres( fp, vertices, radius )
3111     prints 3-d vertices as OFF spheres
3112
3113   notes:
3114     inflated octahedron from Stuart Levy earth/mksphere2
3115 */
3116 void qh_printspheres(FILE *fp, setT *vertices, realT radius) {
3117   vertexT *vertex, **vertexp;
3118
3119   qh printoutnum++;
3120   fprintf (fp, "{appearance {-edge -normal normscale 0} {\n\
3121 INST geom {define vsphere OFF\n\
3122 18 32 48\n\
3123 \n\
3124 0 0 1\n\
3125 1 0 0\n\
3126 0 1 0\n\
3127 -1 0 0\n\
3128 0 -1 0\n\
3129 0 0 -1\n\
3130 0.707107 0 0.707107\n\
3131 0 -0.707107 0.707107\n\
3132 0.707107 -0.707107 0\n\
3133 -0.707107 0 0.707107\n\
3134 -0.707107 -0.707107 0\n\
3135 0 0.707107 0.707107\n\
3136 -0.707107 0.707107 0\n\
3137 0.707107 0.707107 0\n\
3138 0.707107 0 -0.707107\n\
3139 0 0.707107 -0.707107\n\
3140 -0.707107 0 -0.707107\n\
3141 0 -0.707107 -0.707107\n\
3142 \n\
3143 3 0 6 11\n\
3144 3 0 7 6 \n\
3145 3 0 9 7 \n\
3146 3 0 11 9\n\
3147 3 1 6 8 \n\
3148 3 1 8 14\n\
3149 3 1 13 6\n\
3150 3 1 14 13\n\
3151 3 2 11 13\n\
3152 3 2 12 11\n\
3153 3 2 13 15\n\
3154 3 2 15 12\n\
3155 3 3 9 12\n\
3156 3 3 10 9\n\
3157 3 3 12 16\n\
3158 3 3 16 10\n\
3159 3 4 7 10\n\
3160 3 4 8 7\n\
3161 3 4 10 17\n\
3162 3 4 17 8\n\
3163 3 5 14 17\n\
3164 3 5 15 14\n\
3165 3 5 16 15\n\
3166 3 5 17 16\n\
3167 3 6 13 11\n\
3168 3 7 8 6\n\
3169 3 9 10 7\n\
3170 3 11 12 9\n\
3171 3 14 8 17\n\
3172 3 15 13 14\n\
3173 3 16 12 15\n\
3174 3 17 10 16\n} transforms { TLIST\n");
3175   FOREACHvertex_(vertices) {
3176     fprintf(fp, "%8.4g 0 0 0 # v%d\n 0 %8.4g 0 0\n0 0 %8.4g 0\n",
3177       radius, vertex->id, radius, radius);
3178     qh_printpoint3 (fp, vertex->point);
3179     fprintf (fp, "1\n");
3180   }
3181   fprintf (fp, "}}}\n");
3182 } /* printspheres */
3183
3184
3185 /*----------------------------------------------
3186 -printsummary-
3187                 see qhull.c
3188 */
3189
3190 /*-<a                             href="qh-io.htm#TOC"
3191   >-------------------------------</a><a name="printvdiagram">-</a>
3192   
3193   qh_printvdiagram( fp, format, facetlist, facets, printall )
3194     print voronoi diagram
3195       # of pairs of input sites
3196       #indices site1 site2 vertex1 ...
3197     
3198     sites indexed by input point id
3199       point 0 is the first input point
3200     vertices indexed by 'o' and 'p' order
3201       vertex 0 is the 'vertex-at-infinity'
3202       vertex 1 is the first Voronoi vertex
3203
3204   see:
3205     qh_printvoronoi()
3206     qh_eachvoronoi_all()
3207
3208   notes:
3209     if all facets are upperdelaunay, 
3210       prints upper hull (furthest-site Voronoi diagram)
3211 */
3212 void qh_printvdiagram (FILE *fp, int format, facetT *facetlist, setT *facets, boolT printall) {
3213   setT *vertices;
3214   int totcount, numcenters;
3215   boolT islower;
3216   qh_RIDGE innerouter= qh_RIDGEall;
3217   printvridgeT printvridge= NULL;
3218
3219   if (format == qh_PRINTvertices) {
3220     innerouter= qh_RIDGEall;
3221     printvridge= qh_printvridge;
3222   }else if (format == qh_PRINTinner) {
3223     innerouter= qh_RIDGEinner;
3224     printvridge= qh_printvnorm;
3225   }else if (format == qh_PRINTouter) {
3226     innerouter= qh_RIDGEouter;
3227     printvridge= qh_printvnorm;
3228   }else {
3229     fprintf(qh ferr, "qh_printvdiagram: unknown print format %d.\n", format);
3230     qh_errexit (qh_ERRinput, NULL, NULL);
3231   }
3232   vertices= qh_markvoronoi (facetlist, facets, printall, &islower, &numcenters);
3233   totcount= qh_printvdiagram2 (NULL, NULL, vertices, innerouter, False);
3234   fprintf (fp, "%d\n", totcount);
3235   totcount= qh_printvdiagram2 (fp, printvridge, vertices, innerouter, True /* inorder*/);
3236   qh_settempfree (&vertices);
3237 #if 0  /* for testing qh_eachvoronoi_all */
3238   fprintf (fp, "\n");
3239   totcount= qh_eachvoronoi_all(fp, printvridge, qh UPPERdelaunay, innerouter, True /* inorder*/);
3240   fprintf (fp, "%d\n", totcount);
3241 #endif
3242 } /* printvdiagram */
3243   
3244 /*-<a                             href="qh-io.htm#TOC"
3245   >-------------------------------</a><a name="printvdiagram2">-</a>
3246   
3247   qh_printvdiagram2( fp, printvridge, vertices, innerouter, inorder )
3248     visit all pairs of input sites (vertices) for selected Voronoi vertices
3249     vertices may include NULLs
3250   
3251   innerouter:
3252     qh_RIDGEall   print inner ridges (bounded) and outer ridges (unbounded)
3253     qh_RIDGEinner print only inner ridges
3254     qh_RIDGEouter print only outer ridges
3255   
3256   inorder:
3257     print 3-d Voronoi vertices in order
3258   
3259   assumes:
3260     qh_markvoronoi marked facet->visitid for Voronoi vertices
3261     all facet->seen= False
3262     all facet->seen2= True
3263   
3264   returns:
3265     total number of Voronoi ridges 
3266     if printvridge,
3267       calls printvridge( fp, vertex, vertexA, centers) for each ridge
3268       [see qh_eachvoronoi()]
3269   
3270   see:
3271     qh_eachvoronoi_all()
3272 */
3273 int qh_printvdiagram2 (FILE *fp, printvridgeT printvridge, setT *vertices, qh_RIDGE innerouter, boolT inorder) {
3274   int totcount= 0;
3275   int vertex_i, vertex_n;
3276   vertexT *vertex;
3277
3278   FORALLvertices 
3279     vertex->seen= False;
3280   FOREACHvertex_i_(vertices) {
3281     if (vertex) {
3282       if (qh GOODvertex > 0 && qh_pointid(vertex->point)+1 != qh GOODvertex)
3283         continue;
3284       totcount += qh_eachvoronoi (fp, printvridge, vertex, !qh_ALL, innerouter, inorder);
3285     }
3286   }
3287   return totcount;
3288 } /* printvdiagram2 */
3289   
3290 /*-<a                             href="qh-io.htm#TOC"
3291   >-------------------------------</a><a name="printvertex">-</a>
3292   
3293   qh_printvertex( fp, vertex )
3294     prints the information in a vertex
3295 */
3296 void qh_printvertex(FILE *fp, vertexT *vertex) {
3297   pointT *point;
3298   int k, count= 0;
3299   facetT *neighbor, **neighborp;
3300   realT r; /*bug fix*/
3301
3302   if (!vertex) {
3303     fprintf (fp, "  NULLvertex\n");
3304     return;
3305   }
3306   fprintf(fp, "- p%d (v%d):", qh_pointid(vertex->point), vertex->id);
3307   point= vertex->point;
3308   if (point) {
3309     for(k= qh hull_dim; k--; ) {
3310       r= *point++;
3311       fprintf(fp, " %5.2g", r);
3312     }
3313   }
3314   if (vertex->deleted)
3315     fprintf(fp, " deleted");
3316   if (vertex->delridge)
3317     fprintf (fp, " ridgedeleted");
3318   fprintf(fp, "\n");
3319   if (vertex->neighbors) {
3320     fprintf(fp, "  neighbors:");
3321     FOREACHneighbor_(vertex) {
3322       if (++count % 100 == 0)
3323         fprintf (fp, "\n     ");
3324       fprintf(fp, " f%d", neighbor->id);
3325     }
3326     fprintf(fp, "\n");
3327   }
3328 } /* printvertex */
3329
3330
3331 /*-<a                             href="qh-io.htm#TOC"
3332   >-------------------------------</a><a name="printvertexlist">-</a>
3333   
3334   qh_printvertexlist( fp, string, facetlist, facets, printall )
3335     prints vertices used by a facetlist or facet set
3336     tests qh_skipfacet() if !printall
3337 */
3338 void qh_printvertexlist (FILE *fp, char* string, facetT *facetlist, 
3339                          setT *facets, boolT printall) {
3340   vertexT *vertex, **vertexp;
3341   setT *vertices;
3342   
3343   vertices= qh_facetvertices (facetlist, facets, printall);
3344   fputs (string, fp);
3345   FOREACHvertex_(vertices)
3346     qh_printvertex(fp, vertex);
3347   qh_settempfree (&vertices);
3348 } /* printvertexlist */
3349
3350
3351 /*-<a                             href="qh-io.htm#TOC"
3352   >-------------------------------</a><a name="printvertices">-</a>
3353   
3354   qh_printvertices( fp, string, vertices )
3355     prints vertices in a set
3356 */
3357 void qh_printvertices(FILE *fp, char* string, setT *vertices) {
3358   vertexT *vertex, **vertexp;
3359   
3360   fputs (string, fp);
3361   FOREACHvertex_(vertices) 
3362     fprintf (fp, " p%d (v%d)", qh_pointid(vertex->point), vertex->id);
3363   fprintf(fp, "\n");
3364 } /* printvertices */
3365
3366 /*-<a                             href="qh-io.htm#TOC"
3367   >-------------------------------</a><a name="printvneighbors">-</a>
3368   
3369   qh_printvneighbors( fp, facetlist, facets, printall )
3370     print vertex neighbors of vertices in facetlist and facets ('FN')
3371
3372   notes:
3373     qh_countfacets clears facet->visitid for non-printed facets
3374
3375   design:
3376     collect facet count and related statistics
3377     if necessary, build neighbor sets for each vertex
3378     collect vertices in facetlist and facets
3379     build a point array for point->vertex and point->coplanar facet
3380     for each point
3381       list vertex neighbors or coplanar facet
3382 */
3383 void qh_printvneighbors (FILE *fp, facetT* facetlist, setT *facets, boolT printall) {
3384   int numfacets, numsimplicial, numridges, totneighbors, numneighbors, numcoplanars, numtricoplanars;
3385   setT *vertices, *vertex_points, *coplanar_points;
3386   int numpoints= qh num_points + qh_setsize (qh other_points);
3387   vertexT *vertex, **vertexp;
3388   int vertex_i, vertex_n;
3389   facetT *facet, **facetp, *neighbor, **neighborp;
3390   pointT *point, **pointp;
3391
3392   qh_countfacets (facetlist, facets, printall, &numfacets, &numsimplicial, 
3393       &totneighbors, &numridges, &numcoplanars, &numtricoplanars);  /* sets facet->visitid */
3394   fprintf (fp, "%d\n", numpoints);
3395   qh_vertexneighbors();
3396   vertices= qh_facetvertices (facetlist, facets, printall);
3397   vertex_points= qh_settemp (numpoints);
3398   coplanar_points= qh_settemp (numpoints);
3399   qh_setzero (vertex_points, 0, numpoints);
3400   qh_setzero (coplanar_points, 0, numpoints);
3401   FOREACHvertex_(vertices)
3402     qh_point_add (vertex_points, vertex->point, vertex);
3403   FORALLfacet_(facetlist) {
3404     FOREACHpoint_(facet->coplanarset)
3405       qh_point_add (coplanar_points, point, facet);
3406   }
3407   FOREACHfacet_(facets) {
3408     FOREACHpoint_(facet->coplanarset)
3409       qh_point_add (coplanar_points, point, facet);
3410   }
3411   FOREACHvertex_i_(vertex_points) {
3412     if (vertex) { 
3413       numneighbors= qh_setsize (vertex->neighbors);
3414       fprintf (fp, "%d", numneighbors);
3415       if (qh hull_dim == 3)
3416         qh_order_vertexneighbors (vertex);
3417       else if (qh hull_dim >= 4)
3418         qsort (SETaddr_(vertex->neighbors, facetT), numneighbors,
3419              sizeof (facetT *), qh_compare_facetvisit);
3420       FOREACHneighbor_(vertex) 
3421         fprintf (fp, " %d", 
3422                  neighbor->visitid ? neighbor->visitid - 1 : - neighbor->id);
3423       fprintf (fp, "\n");
3424     }else if ((facet= SETelemt_(coplanar_points, vertex_i, facetT)))
3425       fprintf (fp, "1 %d\n",
3426                   facet->visitid ? facet->visitid - 1 : - facet->id);
3427     else
3428       fprintf (fp, "0\n");
3429   }
3430   qh_settempfree (&coplanar_points);
3431   qh_settempfree (&vertex_points);
3432   qh_settempfree (&vertices);
3433 } /* printvneighbors */
3434
3435 /*-<a                             href="qh-io.htm#TOC"
3436   >-------------------------------</a><a name="printvoronoi">-</a>
3437   
3438   qh_printvoronoi( fp, format, facetlist, facets, printall )
3439     print voronoi diagram in 'o' or 'G' format
3440     for 'o' format
3441       prints voronoi centers for each facet and for infinity
3442       for each vertex, lists ids of printed facets or infinity
3443       assumes facetlist and facets are disjoint
3444     for 'G' format
3445       prints an OFF object
3446       adds a 0 coordinate to center
3447       prints infinity but does not list in vertices
3448
3449   see:
3450     qh_printvdiagram()
3451
3452   notes:
3453     if 'o', 
3454       prints a line for each point except "at-infinity"
3455     if all facets are upperdelaunay, 
3456       reverses lower and upper hull
3457 */
3458 void qh_printvoronoi (FILE *fp, int format, facetT *facetlist, setT *facets, boolT printall) {
3459   int k, numcenters, numvertices= 0, numneighbors, numinf, vid=1, vertex_i, vertex_n;
3460   facetT *facet, **facetp, *neighbor, **neighborp;
3461   setT *vertices;
3462   vertexT *vertex;
3463   boolT islower;
3464   unsigned int numfacets= (unsigned int) qh num_facets;
3465
3466   vertices= qh_markvoronoi (facetlist, facets, printall, &islower, &numcenters);
3467   FOREACHvertex_i_(vertices) {
3468     if (vertex) {
3469       numvertices++;
3470       numneighbors = numinf = 0;
3471       FOREACHneighbor_(vertex) {
3472         if (neighbor->visitid == 0)
3473           numinf= 1;
3474         else if (neighbor->visitid < numfacets)
3475           numneighbors++;
3476       }
3477       if (numinf && !numneighbors) {
3478         SETelem_(vertices, vertex_i)= NULL;
3479         numvertices--;
3480       }
3481     }
3482   }
3483   if (format == qh_PRINTgeom) 
3484     fprintf (fp, "{appearance {+edge -face} OFF %d %d 1 # Voronoi centers and cells\n", 
3485                 numcenters, numvertices);
3486   else
3487     fprintf (fp, "%d\n%d %d 1\n", qh hull_dim-1, numcenters, qh_setsize(vertices));
3488   if (format == qh_PRINTgeom) {
3489     for (k= qh hull_dim-1; k--; )
3490       fprintf (fp, qh_REAL_1, 0.0);
3491     fprintf (fp, " 0 # infinity not used\n");
3492   }else {
3493     for (k= qh hull_dim-1; k--; )
3494       fprintf (fp, qh_REAL_1, qh_INFINITE);
3495     fprintf (fp, "\n");
3496   }
3497   FORALLfacet_(facetlist) {
3498     if (facet->visitid && facet->visitid < numfacets) {
3499       if (format == qh_PRINTgeom)
3500         fprintf (fp, "# %d f%d\n", vid++, facet->id);
3501       qh_printcenter (fp, format, NULL, facet);
3502     }
3503   }
3504   FOREACHfacet_(facets) {
3505     if (facet->visitid && facet->visitid < numfacets) {
3506       if (format == qh_PRINTgeom)
3507         fprintf (fp, "# %d f%d\n", vid++, facet->id);
3508       qh_printcenter (fp, format, NULL, facet);
3509     }
3510   }
3511   FOREACHvertex_i_(vertices) {
3512     numneighbors= 0;
3513     numinf=0;
3514     if (vertex) {
3515       if (qh hull_dim == 3)
3516         qh_order_vertexneighbors(vertex);
3517       else if (qh hull_dim >= 4)
3518         qsort (SETaddr_(vertex->neighbors, vertexT), 
3519              qh_setsize (vertex->neighbors),
3520              sizeof (facetT *), qh_compare_facetvisit);
3521       FOREACHneighbor_(vertex) {
3522         if (neighbor->visitid == 0)
3523           numinf= 1;
3524         else if (neighbor->visitid < numfacets)
3525           numneighbors++;
3526       }
3527     }
3528     if (format == qh_PRINTgeom) {
3529       if (vertex) {
3530         fprintf (fp, "%d", numneighbors);
3531         if (vertex) {
3532           FOREACHneighbor_(vertex) {
3533             if (neighbor->visitid && neighbor->visitid < numfacets)
3534               fprintf (fp, " %d", neighbor->visitid);
3535           }
3536         }
3537         fprintf (fp, " # p%d (v%d)\n", vertex_i, vertex->id);
3538       }else
3539         fprintf (fp, " # p%d is coplanar or isolated\n", vertex_i);
3540     }else {
3541       if (numinf)
3542         numneighbors++;
3543       fprintf (fp, "%d", numneighbors);
3544       if (vertex) {
3545         FOREACHneighbor_(vertex) {
3546           if (neighbor->visitid == 0) {
3547             if (numinf) {
3548               numinf= 0;
3549               fprintf (fp, " %d", neighbor->visitid);
3550             }
3551           }else if (neighbor->visitid < numfacets)
3552             fprintf (fp, " %d", neighbor->visitid);
3553         }
3554       }
3555       fprintf (fp, "\n");
3556     }
3557   }
3558   if (format == qh_PRINTgeom)
3559     fprintf (fp, "}\n");
3560   qh_settempfree (&vertices);
3561 } /* printvoronoi */
3562   
3563 /*-<a                             href="qh-io.htm#TOC"
3564   >-------------------------------</a><a name="printvnorm">-</a>
3565   
3566   qh_printvnorm( fp, vertex, vertexA, centers, unbounded )
3567     print one separating plane of the Voronoi diagram for a pair of input sites
3568     unbounded==True if centers includes vertex-at-infinity
3569   
3570   assumes:
3571     qh_ASvoronoi and qh_vertexneighbors() already set
3572     
3573   see:
3574     qh_printvdiagram()
3575     qh_eachvoronoi()
3576 */
3577 void qh_printvnorm (FILE *fp, vertexT *vertex, vertexT *vertexA, setT *centers, boolT unbounded) {
3578   pointT *normal;
3579   realT offset;
3580   int k;
3581   
3582   normal= qh_detvnorm (vertex, vertexA, centers, &offset);
3583   fprintf (fp, "%d %d %d ", 
3584       2+qh hull_dim, qh_pointid (vertex->point), qh_pointid (vertexA->point));
3585   for (k= 0; k< qh hull_dim-1; k++)
3586     fprintf (fp, qh_REAL_1, normal[k]);
3587   fprintf (fp, qh_REAL_1, offset);
3588   fprintf (fp, "\n");
3589 } /* printvnorm */
3590
3591 /*-<a                             href="qh-io.htm#TOC"
3592   >-------------------------------</a><a name="printvridge">-</a>
3593   
3594   qh_printvridge( fp, vertex, vertexA, centers, unbounded )
3595     print one ridge of the Voronoi diagram for a pair of input sites
3596     unbounded==True if centers includes vertex-at-infinity
3597   
3598   see:
3599     qh_printvdiagram()
3600   
3601   notes:
3602     the user may use a different function
3603 */
3604 void qh_printvridge (FILE *fp, vertexT *vertex, vertexT *vertexA, setT *centers, boolT unbounded) {
3605   facetT *facet, **facetp;
3606
3607   fprintf (fp, "%d %d %d", qh_setsize (centers)+2, 
3608        qh_pointid (vertex->point), qh_pointid (vertexA->point));
3609   FOREACHfacet_(centers) 
3610     fprintf (fp, " %d", facet->visitid);
3611   fprintf (fp, "\n");
3612 } /* printvridge */
3613
3614 /*-<a                             href="qh-io.htm#TOC"
3615   >-------------------------------</a><a name="projectdim3">-</a>
3616   
3617   qh_projectdim3( source, destination )
3618     project 2-d 3-d or 4-d point to a 3-d point
3619     uses qh.DROPdim and qh.hull_dim
3620     source and destination may be the same
3621     
3622   notes:
3623     allocate 4 elements to destination just in case
3624 */
3625 void qh_projectdim3 (pointT *source, pointT *destination) {
3626   int i,k;
3627
3628   for (k= 0, i=0; k < qh hull_dim; k++) {
3629     if (qh hull_dim == 4) {
3630       if (k != qh DROPdim)
3631         destination[i++]= source[k];
3632     }else if (k == qh DROPdim)
3633       destination[i++]= 0;
3634     else
3635       destination[i++]= source[k];
3636   }
3637   while (i < 3)
3638     destination[i++]= 0.0;
3639 } /* projectdim3 */
3640
3641 /*-<a                             href="qh-io.htm#TOC"
3642   >-------------------------------</a><a name="readfeasible">-</a>
3643   
3644   qh_readfeasible( dim, remainder )
3645     read feasible point from remainder string and qh.fin
3646
3647   returns:
3648     number of lines read from qh.fin
3649     sets qh.FEASIBLEpoint with malloc'd coordinates
3650
3651   notes:
3652     checks for qh.HALFspace
3653     assumes dim > 1
3654
3655   see:
3656     qh_setfeasible
3657 */
3658 int qh_readfeasible (int dim, char *remainder) {
3659   boolT isfirst= True;
3660   int linecount= 0, tokcount= 0;
3661   char *s, *t, firstline[qh_MAXfirst+1];
3662   coordT *coords, value;
3663
3664   if (!qh HALFspace) {
3665     fprintf  (qh ferr, "qhull input error: feasible point (dim 1 coords) is only valid for halfspace intersection\n");
3666     qh_errexit (qh_ERRinput, NULL, NULL);
3667   }  
3668   if (qh feasible_string)
3669     fprintf  (qh ferr, "qhull input warning: feasible point (dim 1 coords) overrides 'Hn,n,n' feasible point for halfspace intersection\n");
3670   if (!(qh feasible_point= (coordT*)malloc (dim* sizeof(coordT)))) {
3671     fprintf(qh ferr, "qhull error: insufficient memory for feasible point\n");
3672     qh_errexit(qh_ERRmem, NULL, NULL);
3673   }
3674   coords= qh feasible_point;
3675   while ((s= (isfirst ?  remainder : fgets(firstline, qh_MAXfirst, qh fin)))) {
3676     if (isfirst)
3677       isfirst= False;
3678     else
3679       linecount++;
3680     while (*s) {
3681       while (isspace(*s))
3682         s++;
3683       value= qh_strtod (s, &t);
3684       if (s == t)
3685         break;
3686       s= t;
3687       *(coords++)= value;
3688       if (++tokcount == dim) {
3689         while (isspace (*s))
3690           s++;
3691         qh_strtod (s, &t);
3692         if (s != t) {
3693           fprintf (qh ferr, "qhull input error: coordinates for feasible point do not finish out the line: %s\n",
3694                s);
3695           qh_errexit (qh_ERRinput, NULL, NULL);
3696         }
3697         return linecount;
3698       }
3699     }
3700   }
3701   fprintf (qh ferr, "qhull input error: only %d coordinates.  Could not read %d-d feasible point.\n",
3702            tokcount, dim);
3703   qh_errexit (qh_ERRinput, NULL, NULL);
3704   return 0;
3705 } /* readfeasible */
3706
3707 /*-<a                             href="qh-io.htm#TOC"
3708   >-------------------------------</a><a name="readpoints">-</a>
3709   
3710   qh_readpoints( numpoints, dimension, ismalloc )
3711     read points from qh.fin into qh.first_point, qh.num_points
3712     qh.fin is lines of coordinates, one per vertex, first line number of points
3713     if 'rbox D4',
3714       gives message
3715     if qh.ATinfinity,
3716       adds point-at-infinity for Delaunay triangulations
3717
3718   returns:
3719     number of points, array of point coordinates, dimension, ismalloc True
3720     if qh.DELAUNAY & !qh.PROJECTinput, projects points to paraboloid
3721         and clears qh.PROJECTdelaunay
3722     if qh.HALFspace, reads optional feasible point, reads halfspaces,
3723         converts to dual.
3724
3725   for feasible point in "cdd format" in 3-d:
3726     3 1
3727     coordinates
3728     comments
3729     begin
3730     n 4 real/integer
3731     ...
3732     end
3733
3734   notes:
3735     dimension will change in qh_initqhull_globals if qh.PROJECTinput
3736     uses malloc() since qh_mem not initialized
3737     FIXUP: this routine needs rewriting
3738 */
3739 coordT *qh_readpoints(int *numpoints, int *dimension, boolT *ismalloc) {
3740   coordT *points, *coords, *infinity= NULL;
3741   realT paraboloid, maxboloid= -REALmax, value;
3742   realT *coordp= NULL, *offsetp= NULL, *normalp= NULL;
3743   char *s, *t, firstline[qh_MAXfirst+1];
3744   int diminput=0, numinput=0, dimfeasible= 0, newnum, k, tempi;
3745   int firsttext=0, firstshort=0, firstlong=0, firstpoint=0;
3746   int tokcount= 0, linecount=0, maxcount, coordcount=0;
3747   boolT islong, isfirst= True, wasbegin= False;
3748   boolT isdelaunay= qh DELAUNAY && !qh PROJECTinput;
3749
3750   if (qh CDDinput) {
3751     while ((s= fgets(firstline, qh_MAXfirst, qh fin))) {
3752       linecount++;
3753       if (qh HALFspace && linecount == 1 && isdigit(*s)) {
3754         dimfeasible= qh_strtol (s, &s); 
3755         while (isspace(*s))
3756           s++;
3757         if (qh_strtol (s, &s) == 1)
3758           linecount += qh_readfeasible (dimfeasible, s);
3759         else
3760           dimfeasible= 0;
3761       }else if (!memcmp (firstline, "begin", 5) || !memcmp (firstline, "BEGIN", 5))
3762         break;
3763       else if (!*qh rbox_command)
3764         strncat(qh rbox_command, s, sizeof (qh rbox_command)-1);
3765     }
3766     if (!s) {
3767       fprintf (qh ferr, "qhull input error: missing \"begin\" for cdd-formated input\n");
3768       qh_errexit (qh_ERRinput, NULL, NULL);
3769     }
3770   }
3771   while(!numinput && (s= fgets(firstline, qh_MAXfirst, qh fin))) {
3772     linecount++;
3773     if (!memcmp (s, "begin", 5) || !memcmp (s, "BEGIN", 5))
3774       wasbegin= True;
3775     while (*s) {
3776       while (isspace(*s))
3777         s++;
3778       if (!*s)
3779         break;
3780       if (!isdigit(*s)) {
3781         if (!*qh rbox_command) {
3782           strncat(qh rbox_command, s, sizeof (qh rbox_command)-1);
3783           firsttext= linecount;
3784         }
3785         break;
3786       }
3787       if (!diminput) 
3788         diminput= qh_strtol (s, &s);
3789       else {
3790         numinput= qh_strtol (s, &s);
3791         if (numinput == 1 && diminput >= 2 && qh HALFspace && !qh CDDinput) {
3792           linecount += qh_readfeasible (diminput, s); /* checks if ok */
3793           dimfeasible= diminput;
3794           diminput= numinput= 0;
3795         }else 
3796           break;
3797       }
3798     }
3799   }
3800   if (!s) {
3801     fprintf(qh ferr, "qhull input error: short input file.  Did not find dimension and number of points\n");
3802     qh_errexit(qh_ERRinput, NULL, NULL);
3803   }
3804   if (diminput > numinput) {
3805     tempi= diminput;    /* exchange dim and n, e.g., for cdd input format */
3806     diminput= numinput;
3807     numinput= tempi;
3808   }
3809   if (diminput < 2) {
3810     fprintf(qh ferr,"qhull input error: dimension %d (first number) should be at least 2\n",
3811             diminput);
3812     qh_errexit(qh_ERRinput, NULL, NULL);
3813   }
3814   if (isdelaunay) {
3815     qh PROJECTdelaunay= False;
3816     if (qh CDDinput)
3817       *dimension= diminput;
3818     else
3819       *dimension= diminput+1;
3820     *numpoints= numinput;
3821     if (qh ATinfinity)
3822       (*numpoints)++;
3823   }else if (qh HALFspace) {
3824     *dimension= diminput - 1;
3825     *numpoints= numinput;
3826     if (diminput < 3) {
3827       fprintf(qh ferr,"qhull input error: dimension %d (first number, includes offset) should be at least 3 for halfspaces\n",
3828             diminput);
3829       qh_errexit(qh_ERRinput, NULL, NULL);
3830     }
3831     if (dimfeasible) {
3832       if (dimfeasible != *dimension) {
3833         fprintf(qh ferr,"qhull input error: dimension %d of feasible point is not one less than dimension %d for halfspaces\n",
3834           dimfeasible, diminput);
3835         qh_errexit(qh_ERRinput, NULL, NULL);
3836       }
3837     }else 
3838       qh_setfeasible (*dimension);
3839   }else {
3840     if (qh CDDinput) 
3841       *dimension= diminput-1;
3842     else
3843       *dimension= diminput;
3844     *numpoints= numinput;
3845   }
3846   qh normal_size= *dimension * sizeof(coordT); /* for tracing with qh_printpoint */
3847   if (qh HALFspace) {
3848     qh half_space= coordp= (coordT*) malloc (qh normal_size + sizeof(coordT));
3849     if (qh CDDinput) {
3850       offsetp= qh half_space;
3851       normalp= offsetp + 1;
3852     }else {
3853       normalp= qh half_space;
3854       offsetp= normalp + *dimension;
3855     }
3856   } 
3857   qh maxline= diminput * (qh_REALdigits + 5);
3858   maximize_(qh maxline, 500);
3859   qh line= (char*)malloc ((qh maxline+1) * sizeof (char));
3860   *ismalloc= True;  /* use malloc since memory not setup */
3861   coords= points= qh temp_malloc= 
3862         (coordT*)malloc((*numpoints)*(*dimension)*sizeof(coordT));
3863   if (!coords || !qh line || (qh HALFspace && !qh half_space)) {
3864     fprintf(qh ferr, "qhull error: insufficient memory to read %d points\n",
3865             numinput);
3866     qh_errexit(qh_ERRmem, NULL, NULL);
3867   }
3868   if (isdelaunay && qh ATinfinity) {
3869     infinity= points + numinput * (*dimension);
3870     for (k= (*dimension) - 1; k--; )
3871       infinity[k]= 0.0;
3872   }
3873   maxcount= numinput * diminput;
3874   paraboloid= 0.0;
3875   while ((s= (isfirst ?  s : fgets(qh line, qh maxline, qh fin)))) {
3876     if (!isfirst) {
3877       linecount++;
3878       if (*s == 'e' || *s == 'E') {
3879         if (!memcmp (s, "end", 3) || !memcmp (s, "END", 3)) {
3880           if (qh CDDinput )
3881             break;
3882           else if (wasbegin) 
3883             fprintf (qh ferr, "qhull input warning: the input appears to be in cdd format.  If so, use 'Fd'\n");
3884         }
3885       }
3886     }
3887     islong= False;
3888     while (*s) {
3889       while (isspace(*s))
3890         s++;
3891       value= qh_strtod (s, &t);
3892       if (s == t) {
3893         if (!*qh rbox_command)
3894          strncat(qh rbox_command, s, sizeof (qh rbox_command)-1);
3895         if (*s && !firsttext) 
3896           firsttext= linecount;
3897         if (!islong && !firstshort && coordcount)
3898           firstshort= linecount;
3899         break;
3900       }
3901       if (!firstpoint)
3902         firstpoint= linecount;
3903       s= t;
3904       if (++tokcount > maxcount)
3905         continue;
3906       if (qh HALFspace) {
3907         if (qh CDDinput) 
3908           *(coordp++)= -value; /* both coefficients and offset */