Free cookie consent management tool by TermsFeed Policy Generator

source: branches/HeuristicLab.Analysis.AlgorithmBehavior/qhull-2012.1/src/libqhull/io.c @ 10635

Last change on this file since 10635 was 10207, checked in by ascheibe, 11 years ago

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