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