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