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