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