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


hpp-fcl
Author(s):
autogenerated on Fri Jun 2 2023 02:39:01