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