1
/*<html><pre> -<a href="qh-poly.htm"
2
>-------------------------------</a><a name="TOP">-</a>
5
implements polygons and simplices
7
see qh-poly.htm, poly.h and qhull.h
9
frequently used code is in poly.c
11
copyright (c) 1993-2002, The Geometry Center
16
/*======== functions in alphabetical order ==========*/
18
/*-<a href="qh-poly.htm#TOC"
19
>-------------------------------</a><a name="addhash">-</a>
21
qh_addhash( newelem, hashtable, hashsize, hash )
22
add newelem to linear hash table at hash if not already there
24
void qh_addhash (void* newelem, setT *hashtable, int hashsize, unsigned hash) {
28
for (scan= (int)hash; (elem= SETelem_(hashtable, scan));
29
scan= (++scan >= hashsize ? 0 : scan)) {
33
/* loop terminates because qh_HASHfactor >= 1.1 by qh_initbuffers */
35
SETelem_(hashtable, scan)= newelem;
38
/*-<a href="qh-poly.htm#TOC"
39
>-------------------------------</a><a name="check_bestdist">-</a>
42
check that all points are within max_outside of the nearest facet
47
qh_check_maxout(), qh_outerinner()
50
only called from qh_check_points()
51
seldom used since qh.MERGING is almost always set
52
if notverified>0 at end of routine
53
some points were well inside the hull. If the hull contains
54
a lens-shaped component, these points were not verified. Use
55
options 'Qi Tv' to verify all points. (Exhaustive check also verifies)
58
determine facet for each point (if any)
60
start with the assigned facet or with the first facet
61
find the best facet for the point and check all coplanar facets
62
error if point is outside of facet
64
void qh_check_bestdist (void) {
65
boolT waserror= False, unassigned;
66
facetT *facet, *bestfacet, *errfacet1= NULL, *errfacet2= NULL;
68
realT dist, maxoutside, maxdist= -REALmax;
70
int numpart= 0, facet_i, facet_n, notgood= 0, notverified= 0;
73
trace1((qh ferr, "qh_check_bestdist: check points below nearest facet. Facet_list f%d\n",
75
maxoutside= qh_maxouter();
76
maxoutside += qh DISTround;
77
/* one more qh.DISTround for check computation */
78
trace1((qh ferr, "qh_check_bestdist: check that all points are within %2.2g of best facet\n", maxoutside));
79
facets= qh_pointfacet (/*qh facet_list*/);
80
if (!qh_QUICKhelp && qh PRINTprecision)
81
fprintf (qh ferr, "\n\
82
qhull output completed. Verifying that %d points are\n\
83
below %2.2g of the nearest %sfacet.\n",
84
qh_setsize(facets), maxoutside, (qh ONLYgood ? "good " : ""));
85
FOREACHfacet_i_(facets) { /* for each point with facet assignment */
92
point= qh_point(facet_i);
93
if (point == qh GOODpointp)
95
qh_distplane(point, facet, &dist);
97
bestfacet= qh_findbesthorizon (!qh_IScheckmax, point, facet, qh_NOupper, &dist, &numpart);
98
/* occurs after statistics reported */
99
maximize_(maxdist, dist);
100
if (dist > maxoutside) {
101
if (qh ONLYgood && !bestfacet->good
102
&& !((bestfacet= qh_findgooddist (point, bestfacet, &dist, &facetlist))
103
&& dist > maxoutside))
107
fprintf(qh ferr, "qhull precision error: point p%d is outside facet f%d, distance= %6.8g maxoutside= %6.8g\n",
108
facet_i, bestfacet->id, dist, maxoutside);
109
if (errfacet1 != bestfacet) {
110
errfacet2= errfacet1;
111
errfacet1= bestfacet;
114
}else if (unassigned && dist < -qh MAXcoplanar)
117
qh_settempfree (&facets);
118
if (notverified && !qh DELAUNAY && !qh_QUICKhelp && qh PRINTprecision)
119
fprintf(qh ferr, "\n%d points were well inside the hull. If the hull contains\n\
120
a lens-shaped component, these points were not verified. Use\n\
121
options 'Qci Tv' to verify all points.\n", notverified);
122
if (maxdist > qh outside_err) {
123
fprintf( qh ferr, "qhull precision error (qh_check_bestdist): a coplanar point is %6.2g from convex hull. The maximum value (qh.outside_err) is %6.2g\n",
124
maxdist, qh outside_err);
125
qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);
126
}else if (waserror && qh outside_err > REALmax/2)
127
qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);
129
; /* the error was logged to qh.ferr but does not effect the output */
130
trace0((qh ferr, "qh_check_bestdist: max distance outside %2.2g\n", maxdist));
131
} /* check_bestdist */
133
/*-<a href="qh-poly.htm#TOC"
134
>-------------------------------</a><a name="check_maxout">-</a>
137
updates qh.max_outside by checking all points against bestfacet
138
if qh.ONLYgood, ignores !good facets
141
updates facet->maxoutside via qh_findbesthorizon()
143
if printing qh.min_vertex (qh_outerinner),
144
it is updated to the current vertices
145
removes inside/coplanar points from coplanarset as needed
148
defines coplanar as min_vertex instead of MAXcoplanar
149
may not need to check near-inside points because of qh.MAXcoplanar
150
and qh.KEEPnearinside (before it was -DISTround)
156
if qh.min_vertex is needed
157
for all neighbors of all vertices
158
test distance from vertex to neighbor
159
determine facet for each point (if any)
160
for each point with an assigned facet
161
find the best facet for the point and check all coplanar facets
162
(updates outer planes)
163
remove near-inside points from coplanar sets
166
void qh_check_maxout (void) {
167
facetT *facet, *bestfacet, *neighbor, **neighborp, *facetlist;
168
realT dist, maxoutside, minvertex, old_maxoutside;
170
int numpart= 0, facet_i, facet_n, notgood= 0;
171
setT *facets, *vertices;
174
trace1((qh ferr, "qh_check_maxout: check and update maxoutside for each facet.\n"));
175
maxoutside= minvertex= 0;
176
if (qh VERTEXneighbors
177
&& (qh PRINTsummary || qh KEEPinside || qh KEEPcoplanar
178
|| qh TRACElevel || qh PRINTstatistics
179
|| qh PRINTout[0] == qh_PRINTsummary || qh PRINTout[0] == qh_PRINTnone)) {
180
trace1((qh ferr, "qh_check_maxout: determine actual maxoutside and minvertex\n"));
181
vertices= qh_pointvertex (/*qh facet_list*/);
183
FOREACHneighbor_(vertex) {
184
zinc_(Zdistvertex); /* distance also computed by main loop below */
185
qh_distplane (vertex->point, neighbor, &dist);
186
minimize_(minvertex, dist);
187
if (-dist > qh TRACEdist || dist > qh TRACEdist
188
|| neighbor == qh tracefacet || vertex == qh tracevertex)
189
fprintf (qh ferr, "qh_check_maxout: p%d (v%d) is %.2g from f%d\n",
190
qh_pointid (vertex->point), vertex->id, dist, neighbor->id);
194
wmin_(Wminvertex, qh min_vertex);
196
qh min_vertex= minvertex;
197
qh_settempfree (&vertices);
199
facets= qh_pointfacet (/*qh facet_list*/);
201
old_maxoutside= fmax_(qh max_outside, maxoutside);
202
FOREACHfacet_i_(facets) { /* for each point with facet assignment */
204
point= qh_point(facet_i);
205
if (point == qh GOODpointp)
208
qh_distplane(point, facet, &dist);
210
bestfacet= qh_findbesthorizon (qh_IScheckmax, point, facet, !qh_NOupper, &dist, &numpart);
211
if (bestfacet && dist > maxoutside) {
212
if (qh ONLYgood && !bestfacet->good
213
&& !((bestfacet= qh_findgooddist (point, bestfacet, &dist, &facetlist))
214
&& dist > maxoutside))
219
if (dist > qh TRACEdist || (bestfacet && bestfacet == qh tracefacet))
220
fprintf (qh ferr, "qh_check_maxout: p%d is %.2g above f%d\n",
221
qh_pointid (point), dist, bestfacet->id);
225
(maxoutside > 2*old_maxoutside);
226
/* if qh.maxoutside increases substantially, qh_SEARCHdist is not valid
227
e.g., RBOX 5000 s Z1 G1e-13 t1001200614 | qhull */
228
zzadd_(Zcheckpart, numpart);
229
qh_settempfree (&facets);
230
wval_(Wmaxout)= maxoutside - qh max_outside;
231
wmax_(Wmaxoutside, qh max_outside);
232
qh max_outside= maxoutside;
233
qh_nearcoplanar (/*qh.facet_list*/);
235
trace1((qh ferr, "qh_check_maxout: maxoutside %2.2g, min_vertex %2.2g, outside of not good %d\n",
236
maxoutside, qh min_vertex, notgood));
238
#else /* qh_NOmerge */
239
void qh_check_maxout (void) {
243
/*-<a href="qh-poly.htm#TOC"
244
>-------------------------------</a><a name="check_output">-</a>
247
performs the checks at the end of qhull algorithm
248
Maybe called after voronoi output. Will recompute otherwise centrums are Voronoi centers instead
250
void qh_check_output (void) {
255
if (qh VERIFYoutput | qh IStracing | qh CHECKfrequently) {
256
qh_checkpolygon (qh facet_list);
257
qh_checkflipped_all (qh facet_list);
258
qh_checkconvex (qh facet_list, qh_ALGORITHMfault);
259
}else if (!qh MERGING && qh_newstats (qhstat precision, &i)) {
260
qh_checkflipped_all (qh facet_list);
261
qh_checkconvex (qh facet_list, qh_ALGORITHMfault);
267
/*-<a href="qh-poly.htm#TOC"
268
>-------------------------------</a><a name="check_point">-</a>
270
qh_check_point( point, facet, maxoutside, maxdist, errfacet1, errfacet2 )
271
check that point is less than maxoutside from facet
273
void qh_check_point (pointT *point, facetT *facet, realT *maxoutside, realT *maxdist, facetT **errfacet1, facetT **errfacet2) {
276
/* occurs after statistics reported */
277
qh_distplane(point, facet, &dist);
278
if (dist > *maxoutside) {
279
if (*errfacet1 != facet) {
280
*errfacet2= *errfacet1;
283
fprintf(qh ferr, "qhull precision error: point p%d is outside facet f%d, distance= %6.8g maxoutside= %6.8g\n",
284
qh_pointid(point), facet->id, dist, *maxoutside);
286
maximize_(*maxdist, dist);
287
} /* qh_check_point */
290
/*-<a href="qh-poly.htm#TOC"
291
>-------------------------------</a><a name="check_points">-</a>
294
checks that all points are inside all facets
297
if many points and qh_check_maxout not called (i.e., !qh.MERGING),
298
calls qh_findbesthorizon (seldom done).
299
ignores flipped facets
300
maxoutside includes 2 qh.DISTrounds
301
one qh.DISTround for the computed distances in qh_check_points
302
qh_printafacet and qh_printsummary needs only one qh.DISTround
303
the computation for qh.VERIFYdirect does not account for qh.other_points
307
use qh_check_bestdist()
311
check that point is inside facet
313
void qh_check_points (void) {
314
facetT *facet, *errfacet1= NULL, *errfacet2= NULL;
315
realT total, maxoutside, maxdist= -REALmax;
316
pointT *point, **pointp, *pointtemp;
319
maxoutside= qh_maxouter();
320
maxoutside += qh DISTround;
321
/* one more qh.DISTround for check computation */
322
trace1((qh ferr, "qh_check_points: check all points below %2.2g of all facet planes\n",
324
if (qh num_good) /* miss counts other_points and !good facets */
325
total= (float) qh num_good * qh num_points;
327
total= (float) qh num_facets * qh num_points;
328
if (total >= qh_VERIFYdirect && !qh maxoutdone) {
329
if (!qh_QUICKhelp && qh SKIPcheckmax && qh MERGING)
330
fprintf (qh ferr, "\n\
331
qhull input warning: merging without checking outer planes ('Q5' or 'Po').\n\
332
Verify may report that a point is outside of a facet.\n");
335
if (qh_MAXoutside && qh maxoutdone)
341
fprintf (qh ferr, "\n\
342
qhull input warning: exact merge ('Qx'). Verify may report that a point\n\
343
is outside of a facet. See qh-optq.htm#Qx\n");
344
else if (qh SKIPcheckmax || qh NOnearinside)
345
fprintf (qh ferr, "\n\
346
qhull input warning: no outer plane check ('Q5') or no processing of\n\
347
near-inside points ('Q8'). Verify may report that a point is outside\n\
350
if (qh PRINTprecision) {
352
fprintf (qh ferr, "\n\
353
Output completed. Verifying that all points are below outer planes of\n\
354
all %sfacets. Will make %2.0f distance computations.\n",
355
(qh ONLYgood ? "good " : ""), total);
357
fprintf (qh ferr, "\n\
358
Output completed. Verifying that all points are below %2.2g of\n\
359
all %sfacets. Will make %2.0f distance computations.\n",
360
maxoutside, (qh ONLYgood ? "good " : ""), total);
363
if (!facet->good && qh ONLYgood)
367
if (!facet->normal) {
368
fprintf( qh ferr, "qhull warning (qh_check_points): missing normal for facet f%d\n", facet->id);
373
maxoutside= facet->maxoutside + 2* qh DISTround;
374
/* one DISTround to actual point and another to computed point */
378
if (point != qh GOODpointp)
379
qh_check_point (point, facet, &maxoutside, &maxdist, &errfacet1, &errfacet2);
381
FOREACHpoint_(qh other_points) {
382
if (point != qh GOODpointp)
383
qh_check_point (point, facet, &maxoutside, &maxdist, &errfacet1, &errfacet2);
386
if (maxdist > qh outside_err) {
387
fprintf( qh ferr, "qhull precision error (qh_check_points): a coplanar point is %6.2g from convex hull. The maximum value (qh.outside_err) is %6.2g\n",
388
maxdist, qh outside_err );
389
qh_errexit2( qh_ERRprec, errfacet1, errfacet2 );
390
}else if (errfacet1 && qh outside_err > REALmax/2)
391
qh_errexit2( qh_ERRprec, errfacet1, errfacet2 );
393
; /* the error was logged to qh.ferr but does not effect the output */
394
trace0((qh ferr, "qh_check_points: max distance outside %2.2g\n", maxdist));
399
/*-<a href="qh-poly.htm#TOC"
400
>-------------------------------</a><a name="checkconvex">-</a>
402
qh_checkconvex( facetlist, fault )
403
check that each ridge in facetlist is convex
404
fault = qh_DATAfault if reporting errors
405
= qh_ALGORITHMfault otherwise
408
counts Zconcaveridges and Zcoplanarridges
409
errors if concaveridge or if merging an coplanar ridge
413
tests vertices for neighboring simplicial facets
415
tests vertices for neighboring simplicial facets
417
tests centrums of neighboring facets
421
report flipped facets
422
if ZEROcentrum and simplicial neighbors
423
test vertices for neighboring simplicial facets
425
test centrum against all neighbors
427
void qh_checkconvex(facetT *facetlist, int fault) {
428
facetT *facet, *neighbor, **neighborp, *errfacet1=NULL, *errfacet2=NULL;
432
boolT waserror= False, centrum_warning= False, tempcentrum= False, allsimplicial;
435
trace1((qh ferr, "qh_checkconvex: check all ridges are convex\n"));
437
zzval_(Zconcaveridges)= 0;
438
zzval_(Zcoplanarridges)= 0;
440
FORALLfacet_(facetlist) {
441
if (facet->flipped) {
442
qh_precision ("flipped facet");
443
fprintf (qh ferr, "qhull precision error: f%d is flipped (interior point is outside)\n",
449
if (qh MERGING && (!qh ZEROcentrum || !facet->simplicial || facet->tricoplanar))
450
allsimplicial= False;
454
FOREACHneighbor_(facet) {
455
vertex= SETelemt_(facet->vertices, neighbor_i++, vertexT);
456
if (!neighbor->simplicial || neighbor->tricoplanar) {
457
allsimplicial= False;
460
qh_distplane (vertex->point, neighbor, &dist);
461
if (dist > -qh DISTround) {
462
if (fault == qh_DATAfault) {
463
qh_precision ("coplanar or concave ridge");
464
fprintf (qh ferr, "qhull precision error: initial simplex is not convex. Distance=%.2g\n", dist);
465
qh_errexit(qh_ERRsingular, NULL, NULL);
467
if (dist > qh DISTround) {
468
zzinc_(Zconcaveridges);
469
qh_precision ("concave ridge");
470
fprintf (qh ferr, "qhull precision error: f%d is concave to f%d, since p%d (v%d) is %6.4g above\n",
471
facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist);
475
}else if (qh ZEROcentrum) {
476
if (dist > 0) { /* qh_checkzero checks that dist < - qh DISTround */
477
zzinc_(Zcoplanarridges);
478
qh_precision ("coplanar ridge");
479
fprintf (qh ferr, "qhull precision error: f%d is clearly not convex to f%d, since p%d (v%d) is %6.4g above\n",
480
facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist);
486
zzinc_(Zcoplanarridges);
487
qh_precision ("coplanar ridge");
488
trace0((qh ferr, "qhull precision error: f%d may be coplanar to f%d, since p%d (v%d) is within %6.4g during p%d\n",
489
facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist, qh furthest_id));
494
if (!allsimplicial) {
495
if (qh CENTERtype == qh_AScentrum) {
497
facet->center= qh_getcentrum (facet);
498
centrum= facet->center;
500
if (!centrum_warning && (!facet->simplicial || facet->tricoplanar)) {
501
centrum_warning= True;
502
fprintf (qh ferr, "qhull note: recomputing centrums for convexity test. This may lead to false, precision errors.\n");
504
centrum= qh_getcentrum(facet);
507
FOREACHneighbor_(facet) {
508
if (qh ZEROcentrum && facet->simplicial && neighbor->simplicial)
510
if (facet->tricoplanar || neighbor->tricoplanar)
513
qh_distplane (centrum, neighbor, &dist);
514
if (dist > qh DISTround) {
515
zzinc_(Zconcaveridges);
516
qh_precision ("concave ridge");
517
fprintf (qh ferr, "qhull precision error: f%d is concave to f%d. Centrum of f%d is %6.4g above f%d\n",
518
facet->id, neighbor->id, facet->id, dist, neighbor->id);
522
}else if (dist >= 0.0) { /* if arithmetic always rounds the same,
523
can test against centrum radius instead */
524
zzinc_(Zcoplanarridges);
525
qh_precision ("coplanar ridge");
526
fprintf (qh ferr, "qhull precision error: f%d is coplanar or concave to f%d. Centrum of f%d is %6.4g above f%d\n",
527
facet->id, neighbor->id, facet->id, dist, neighbor->id);
534
qh_memfree(centrum, qh normal_size);
537
if (waserror && !qh FORCEoutput)
538
qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);
542
/*-<a href="qh-poly.htm#TOC"
543
>-------------------------------</a><a name="checkfacet">-</a>
545
qh_checkfacet( facet, newmerge, waserror )
546
checks for consistency errors in facet
547
newmerge set if from merge.c
550
sets waserror if any error occurs
553
vertex ids are inverse sorted
554
unless newmerge, at least hull_dim neighbors and vertices (exactly if simplicial)
555
if non-simplicial, at least as many ridges as neighbors
556
neighbors are not duplicated
557
ridges are not duplicated
558
in 3-d, ridges=verticies
559
(qh.hull_dim-1) ridge vertices
560
neighbors are reciprocated
561
ridge neighbors are facet neighbors and a ridge for every neighbor
562
simplicial neighbors match facetintersect
563
vertex intersection matches vertices of common ridges
564
vertex neighbors and facet vertices agree
565
all ridges have distinct vertex sets
573
check sizes of neighbors and vertices
574
check for qh_MERGEridge and qh_DUPLICATEridge flags
577
check ridges, neighbors, and vertices
579
void qh_checkfacet(facetT *facet, boolT newmerge, boolT *waserrorp) {
580
facetT *neighbor, **neighborp, *errother=NULL;
581
ridgeT *ridge, **ridgep, *errridge= NULL, *ridge2;
582
vertexT *vertex, **vertexp;
583
unsigned previousid= INT_MAX;
584
int numneighbors, numvertices, numridges=0, numRvertices=0;
585
boolT waserror= False;
586
int skipA, skipB, ridge_i, ridge_n, i;
589
if (facet->visible) {
590
fprintf (qh ferr, "qhull internal error (qh_checkfacet): facet f%d is on the visible_list\n",
592
qh_errexit (qh_ERRqhull, facet, NULL);
594
if (!facet->normal) {
595
fprintf (qh ferr, "qhull internal error (qh_checkfacet): facet f%d does not have a normal\n",
599
qh_setcheck (facet->vertices, "vertices for f", facet->id);
600
qh_setcheck (facet->ridges, "ridges for f", facet->id);
601
qh_setcheck (facet->outsideset, "outsideset for f", facet->id);
602
qh_setcheck (facet->coplanarset, "coplanarset for f", facet->id);
603
qh_setcheck (facet->neighbors, "neighbors for f", facet->id);
604
FOREACHvertex_(facet->vertices) {
605
if (vertex->deleted) {
606
fprintf(qh ferr, "qhull internal error (qh_checkfacet): deleted vertex v%d in f%d\n", vertex->id, facet->id);
607
qh_errprint ("ERRONEOUS", NULL, NULL, NULL, vertex);
610
if (vertex->id >= previousid) {
611
fprintf(qh ferr, "qhull internal error (qh_checkfacet): vertices of f%d are not in descending id order at v%d\n", facet->id, vertex->id);
615
previousid= vertex->id;
617
numneighbors= qh_setsize(facet->neighbors);
618
numvertices= qh_setsize(facet->vertices);
619
numridges= qh_setsize(facet->ridges);
620
if (facet->simplicial) {
621
if (numvertices+numneighbors != 2*qh hull_dim
622
&& !facet->degenerate && !facet->redundant) {
623
fprintf(qh ferr, "qhull internal error (qh_checkfacet): for simplicial facet f%d, #vertices %d + #neighbors %d != 2*qh hull_dim\n",
624
facet->id, numvertices, numneighbors);
625
qh_setprint (qh ferr, "", facet->neighbors);
628
}else { /* non-simplicial */
630
&&(numvertices < qh hull_dim || numneighbors < qh hull_dim)
631
&& !facet->degenerate && !facet->redundant) {
632
fprintf(qh ferr, "qhull internal error (qh_checkfacet): for facet f%d, #vertices %d or #neighbors %d < qh hull_dim\n",
633
facet->id, numvertices, numneighbors);
636
/* in 3-d, can get a vertex twice in an edge list, e.g., RBOX 1000 s W1e-13 t995849315 D2 | QHULL d Tc Tv TP624 TW1e-13 T4 */
637
if (numridges < numneighbors
638
||(qh hull_dim == 3 && numvertices > numridges && !qh NEWfacets)
639
||(qh hull_dim == 2 && numridges + numvertices + numneighbors != 6)) {
640
if (!facet->degenerate && !facet->redundant) {
641
fprintf(qh ferr, "qhull internal error (qh_checkfacet): for facet f%d, #ridges %d < #neighbors %d or (3-d) > #vertices %d or (2-d) not all 2\n",
642
facet->id, numridges, numneighbors, numvertices);
647
FOREACHneighbor_(facet) {
648
if (neighbor == qh_MERGEridge || neighbor == qh_DUPLICATEridge) {
649
fprintf(qh ferr, "qhull internal error (qh_checkfacet): facet f%d still has a MERGE or DUP neighbor\n", facet->id);
650
qh_errexit (qh_ERRqhull, facet, NULL);
652
neighbor->seen= True;
654
FOREACHneighbor_(facet) {
655
if (!qh_setin(neighbor->neighbors, facet)) {
656
fprintf(qh ferr, "qhull internal error (qh_checkfacet): facet f%d has neighbor f%d, but f%d does not have neighbor f%d\n",
657
facet->id, neighbor->id, neighbor->id, facet->id);
661
if (!neighbor->seen) {
662
fprintf(qh ferr, "qhull internal error (qh_checkfacet): facet f%d has a duplicate neighbor f%d\n",
663
facet->id, neighbor->id);
667
neighbor->seen= False;
669
FOREACHridge_(facet->ridges) {
670
qh_setcheck (ridge->vertices, "vertices for r", ridge->id);
673
FOREACHridge_(facet->ridges) {
675
fprintf(qh ferr, "qhull internal error (qh_checkfacet): facet f%d has a duplicate ridge r%d\n",
676
facet->id, ridge->id);
681
numRvertices= qh_setsize(ridge->vertices);
682
if (numRvertices != qh hull_dim - 1) {
683
fprintf(qh ferr, "qhull internal error (qh_checkfacet): ridge between f%d and f%d has %d vertices\n",
684
ridge->top->id, ridge->bottom->id, numRvertices);
688
neighbor= otherfacet_(ridge, facet);
689
neighbor->seen= True;
690
if (!qh_setin(facet->neighbors, neighbor)) {
691
fprintf(qh ferr, "qhull internal error (qh_checkfacet): for facet f%d, neighbor f%d of ridge r%d not in facet\n",
692
facet->id, neighbor->id, ridge->id);
697
if (!facet->simplicial) {
698
FOREACHneighbor_(facet) {
699
if (!neighbor->seen) {
700
fprintf(qh ferr, "qhull internal error (qh_checkfacet): facet f%d does not have a ridge for neighbor f%d\n",
701
facet->id, neighbor->id);
705
intersection= qh_vertexintersect_new(facet->vertices, neighbor->vertices);
706
qh_settemppush (intersection);
707
FOREACHvertex_(facet->vertices) {
709
vertex->seen2= False;
711
FOREACHvertex_(intersection)
713
FOREACHridge_(facet->ridges) {
714
if (neighbor != otherfacet_(ridge, facet))
716
FOREACHvertex_(ridge->vertices) {
718
fprintf (qh ferr, "qhull internal error (qh_checkfacet): vertex v%d in r%d not in f%d intersect f%d\n",
719
vertex->id, ridge->id, facet->id, neighbor->id);
720
qh_errexit (qh_ERRqhull, facet, ridge);
726
FOREACHvertex_(intersection) {
727
if (!vertex->seen2) {
728
if (qh IStracing >=3 || !qh MERGING) {
729
fprintf (qh ferr, "qhull precision error (qh_checkfacet): vertex v%d in f%d intersect f%d but\n\
730
not in a ridge. This is ok under merging. Last point was p%d\n",
731
vertex->id, facet->id, neighbor->id, qh furthest_id);
732
if (!qh FORCEoutput && !qh MERGING) {
733
qh_errprint ("ERRONEOUS", facet, neighbor, NULL, vertex);
735
qh_errexit (qh_ERRqhull, NULL, NULL);
741
qh_settempfree (&intersection);
743
}else { /* simplicial */
744
FOREACHneighbor_(facet) {
745
if (neighbor->simplicial) {
746
skipA= SETindex_(facet->neighbors, neighbor);
747
skipB= qh_setindex (neighbor->neighbors, facet);
748
if (!qh_setequal_skip (facet->vertices, skipA, neighbor->vertices, skipB)) {
749
fprintf (qh ferr, "qhull internal error (qh_checkfacet): facet f%d skip %d and neighbor f%d skip %d do not match \n",
750
facet->id, skipA, neighbor->id, skipB);
757
if (qh hull_dim < 5 && (qh IStracing > 2 || qh CHECKfrequently)) {
758
FOREACHridge_i_(facet->ridges) { /* expensive */
759
for (i= ridge_i+1; i < ridge_n; i++) {
760
ridge2= SETelemt_(facet->ridges, i, ridgeT);
761
if (qh_setequal (ridge->vertices, ridge2->vertices)) {
762
fprintf (qh ferr, "qh_checkfacet: ridges r%d and r%d have the same vertices\n",
763
ridge->id, ridge2->id);
771
qh_errprint("ERRONEOUS", facet, errother, errridge, NULL);
777
/*-<a href="qh-poly.htm#TOC"
778
>-------------------------------</a><a name="checkflipped_all">-</a>
780
qh_checkflipped_all( facetlist )
781
checks orientation of facets in list against interior point
783
void qh_checkflipped_all (facetT *facetlist) {
785
boolT waserror= False;
788
if (facetlist == qh facet_list)
789
zzval_(Zflippedfacets)= 0;
790
FORALLfacet_(facetlist) {
791
if (facet->normal && !qh_checkflipped (facet, &dist, !qh_ALL)) {
792
fprintf(qh ferr, "qhull precision error: facet f%d is flipped, distance= %6.12g\n",
794
if (!qh FORCEoutput) {
795
qh_errprint("ERRONEOUS", facet, NULL, NULL, NULL);
801
fprintf (qh ferr, "\n\
802
A flipped facet occurs when its distance to the interior point is\n\
803
greater than %2.2g, the maximum roundoff error.\n", -qh DISTround);
804
qh_errexit(qh_ERRprec, NULL, NULL);
806
} /* checkflipped_all */
808
/*-<a href="qh-poly.htm#TOC"
809
>-------------------------------</a><a name="checkpolygon">-</a>
811
qh_checkpolygon( facetlist )
812
checks the correctness of the structure
815
call with either qh.facet_list or qh.newfacet_list
816
checks num_facets and num_vertices if qh.facet_list
820
checks facet and outside set
821
initializes vertexlist
824
if checking all facets (qh.facetlist)
826
if qh.VERTEXneighbors
827
check vertex neighbors and count
830
void qh_checkpolygon(facetT *facetlist) {
832
vertexT *vertex, **vertexp, *vertexlist;
833
int numfacets= 0, numvertices= 0, numridges= 0;
834
int totvneighbors= 0, totvertices= 0;
835
boolT waserror= False, nextseen= False, visibleseen= False;
837
trace1((qh ferr, "qh_checkpolygon: check all facets from f%d\n", facetlist->id));
838
if (facetlist != qh facet_list || qh ONLYgood)
840
FORALLfacet_(facetlist) {
841
if (facet == qh visible_list)
843
if (!facet->visible) {
845
if (facet == qh facet_next)
847
else if (qh_setsize (facet->outsideset)) {
849
#if !qh_COMPUTEfurthest
850
|| facet->furthestdist >= qh MINoutside
853
fprintf (qh ferr, "qhull internal error (qh_checkpolygon): f%d has outside points before qh facet_next\n",
855
qh_errexit (qh_ERRqhull, facet, NULL);
860
qh_checkfacet(facet, False, &waserror);
863
if (qh visible_list && !visibleseen && facetlist == qh facet_list) {
864
fprintf (qh ferr, "qhull internal error (qh_checkpolygon): visible list f%d no longer on facet list\n", qh visible_list->id);
866
qh_errexit (qh_ERRqhull, qh visible_list, NULL);
868
if (facetlist == qh facet_list)
869
vertexlist= qh vertex_list;
870
else if (facetlist == qh newfacet_list)
871
vertexlist= qh newvertex_list;
874
FORALLvertex_(vertexlist) {
878
FORALLfacet_(facetlist) {
881
if (facet->simplicial)
882
numridges += qh hull_dim;
884
numridges += qh_setsize (facet->ridges);
885
FOREACHvertex_(facet->vertices) {
890
if (qh_pointid (vertex->point) == -1) {
891
fprintf (qh ferr, "qhull internal error (qh_checkpolygon): unknown point %p for vertex v%d first_point %p\n",
892
vertex->point, vertex->id, qh first_point);
898
qh vertex_visit += numfacets;
899
if (facetlist == qh facet_list) {
900
if (numfacets != qh num_facets - qh num_visible) {
901
fprintf(qh ferr, "qhull internal error (qh_checkpolygon): actual number of facets is %d, cumulative facet count is %d - %d visible facets\n",
902
numfacets, qh num_facets, qh num_visible);
906
if (qh VERTEXneighbors) {
908
qh_setcheck (vertex->neighbors, "neighbors for v", vertex->id);
911
totvneighbors += qh_setsize (vertex->neighbors);
913
FORALLfacet_(facetlist)
914
totvertices += qh_setsize (facet->vertices);
915
if (totvneighbors != totvertices) {
916
fprintf(qh ferr, "qhull internal error (qh_checkpolygon): vertex neighbors inconsistent. Totvneighbors %d, totvertices %d\n",
917
totvneighbors, totvertices);
921
if (numvertices != qh num_vertices - qh_setsize(qh del_vertices)) {
922
fprintf(qh ferr, "qhull internal error (qh_checkpolygon): actual number of vertices is %d, cumulative vertex count is %d\n",
923
numvertices, qh num_vertices - qh_setsize(qh del_vertices));
926
if (qh hull_dim == 2 && numvertices != numfacets) {
927
fprintf (qh ferr, "qhull internal error (qh_checkpolygon): #vertices %d != #facets %d\n",
928
numvertices, numfacets);
931
if (qh hull_dim == 3 && numvertices + numfacets - numridges/2 != 2) {
932
fprintf (qh ferr, "qhull warning: #vertices %d + #facets %d - #edges %d != 2\n\
933
A vertex appears twice in a edge list. May occur during merging.",
934
numvertices, numfacets, numridges/2);
935
/* occurs if lots of merging and a vertex ends up twice in an edge list. e.g., RBOX 1000 s W1e-13 t995849315 D2 | QHULL d Tc Tv */
939
qh_errexit(qh_ERRqhull, NULL, NULL);
943
/*-<a href="qh-poly.htm#TOC"
944
>-------------------------------</a><a name="checkvertex">-</a>
946
qh_checkvertex( vertex )
947
check vertex for consistency
948
checks vertex->neighbors
951
neighbors checked efficiently in checkpolygon
953
void qh_checkvertex (vertexT *vertex) {
954
boolT waserror= False;
955
facetT *neighbor, **neighborp, *errfacet=NULL;
957
if (qh_pointid (vertex->point) == -1) {
958
fprintf (qh ferr, "qhull internal error (qh_checkvertex): unknown point id %p\n", vertex->point);
961
if (vertex->id >= qh vertex_id) {
962
fprintf (qh ferr, "qhull internal error (qh_checkvertex): unknown vertex id %d\n", vertex->id);
965
if (!waserror && !vertex->deleted) {
966
if (qh_setsize (vertex->neighbors)) {
967
FOREACHneighbor_(vertex) {
968
if (!qh_setin (neighbor->vertices, vertex)) {
969
fprintf (qh ferr, "qhull internal error (qh_checkvertex): neighbor f%d does not contain v%d\n", neighbor->id, vertex->id);
977
qh_errprint ("ERRONEOUS", NULL, NULL, NULL, vertex);
978
qh_errexit (qh_ERRqhull, errfacet, NULL);
982
/*-<a href="qh-poly.htm#TOC"
983
>-------------------------------</a><a name="clearcenters">-</a>
985
qh_clearcenters( type )
986
clear old data from facet->center
990
nop if CENTERtype is the same
992
void qh_clearcenters (qh_CENTER type) {
995
if (qh CENTERtype != type) {
997
if (qh CENTERtype == qh_ASvoronoi){
999
qh_memfree (facet->center, qh center_size);
1000
facet->center= NULL;
1002
}else /* qh CENTERtype == qh_AScentrum */ {
1003
if (facet->center) {
1004
qh_memfree (facet->center, qh normal_size);
1005
facet->center= NULL;
1009
qh CENTERtype= type;
1011
trace2((qh ferr, "qh_clearcenters: switched to center type %d\n", type));
1012
} /* clearcenters */
1014
/*-<a href="qh-poly.htm#TOC"
1015
>-------------------------------</a><a name="createsimplex">-</a>
1017
qh_createsimplex( vertices )
1018
creates a simplex from a set of vertices
1021
initializes qh.facet_list to the simplex
1022
initializes qh.newfacet_list, .facet_tail
1023
initializes qh.vertex_list, .newvertex_list, .vertex_tail
1030
create its neighbor set
1032
void qh_createsimplex(setT *vertices) {
1033
facetT *facet= NULL, *newfacet;
1034
boolT toporient= True;
1035
int vertex_i, vertex_n, nth;
1036
setT *newfacets= qh_settemp (qh hull_dim+1);
1039
qh facet_list= qh newfacet_list= qh facet_tail= qh_newfacet();
1040
qh num_facets= qh num_vertices= qh num_visible= 0;
1041
qh vertex_list= qh newvertex_list= qh vertex_tail= qh_newvertex(NULL);
1042
FOREACHvertex_i_(vertices) {
1043
newfacet= qh_newfacet();
1044
newfacet->vertices= qh_setnew_delnthsorted (vertices, vertex_n,
1046
newfacet->toporient= toporient;
1047
qh_appendfacet(newfacet);
1048
newfacet->newfacet= True;
1049
qh_appendvertex (vertex);
1050
qh_setappend (&newfacets, newfacet);
1055
FORALLfacet_(qh newfacet_list) {
1056
if (facet != newfacet)
1057
SETelem_(newfacet->neighbors, nth++)= facet;
1059
qh_settruncate (newfacet->neighbors, qh hull_dim);
1061
qh_settempfree (&newfacets);
1062
trace1((qh ferr, "qh_createsimplex: created simplex\n"));
1063
} /* createsimplex */
1065
/*-<a href="qh-poly.htm#TOC"
1066
>-------------------------------</a><a name="delridge">-</a>
1068
qh_delridge( ridge )
1069
deletes ridge from data structures it belongs to
1073
in merge.c, caller sets vertex->delridge for each vertex
1074
ridges also freed in qh_freeqhull
1076
void qh_delridge(ridgeT *ridge) {
1077
void **freelistp; /* used !qh_NOmem */
1079
qh_setdel(ridge->top->ridges, ridge);
1080
qh_setdel(ridge->bottom->ridges, ridge);
1081
qh_setfree(&(ridge->vertices));
1082
qh_memfree_(ridge, sizeof(ridgeT), freelistp);
1086
/*-<a href="qh-poly.htm#TOC"
1087
>-------------------------------</a><a name="delvertex">-</a>
1089
qh_delvertex( vertex )
1090
deletes a vertex and frees its memory
1093
assumes vertex->adjacencies have been updated if needed
1094
unlinks from vertex_list
1096
void qh_delvertex (vertexT *vertex) {
1098
if (vertex == qh tracevertex)
1099
qh tracevertex= NULL;
1100
qh_removevertex (vertex);
1101
qh_setfree (&vertex->neighbors);
1102
qh_memfree(vertex, sizeof(vertexT));
1106
/*-<a href="qh-poly.htm#TOC"
1107
>-------------------------------</a><a name="facet3vertex">-</a>
1110
return temporary set of 3-d vertices in qh_ORIENTclock order
1114
build set from facet->vertices with facet->toporient
1116
for each ridge in order
1117
build set from ridge's vertices
1119
setT *qh_facet3vertex (facetT *facet) {
1120
ridgeT *ridge, *firstridge;
1122
int cntvertices, cntprojected=0;
1125
cntvertices= qh_setsize(facet->vertices);
1126
vertices= qh_settemp (cntvertices);
1127
if (facet->simplicial) {
1128
if (cntvertices != 3) {
1129
fprintf (qh ferr, "qhull internal error (qh_facet3vertex): only %d vertices for simplicial facet f%d\n",
1130
cntvertices, facet->id);
1131
qh_errexit(qh_ERRqhull, facet, NULL);
1133
qh_setappend (&vertices, SETfirst_(facet->vertices));
1134
if (facet->toporient ^ qh_ORIENTclock)
1135
qh_setappend (&vertices, SETsecond_(facet->vertices));
1137
qh_setaddnth (&vertices, 0, SETsecond_(facet->vertices));
1138
qh_setappend (&vertices, SETelem_(facet->vertices, 2));
1140
ridge= firstridge= SETfirstt_(facet->ridges, ridgeT); /* no infinite */
1141
while ((ridge= qh_nextridge3d (ridge, facet, &vertex))) {
1142
qh_setappend (&vertices, vertex);
1143
if (++cntprojected > cntvertices || ridge == firstridge)
1146
if (!ridge || cntprojected != cntvertices) {
1147
fprintf (qh ferr, "qhull internal error (qh_facet3vertex): ridges for facet %d don't match up. got at least %d\n",
1148
facet->id, cntprojected);
1149
qh_errexit(qh_ERRqhull, facet, ridge);
1153
} /* facet3vertex */
1155
/*-<a href="qh-poly.htm#TOC"
1156
>-------------------------------</a><a name="findbestfacet">-</a>
1158
qh_findbestfacet( point, bestoutside, bestdist, isoutside )
1159
find facet that is furthest below a point
1161
for Delaunay triangulations,
1162
Use qh_setdelaunay() to lift point to paraboloid and scale by 'Qbb' if needed
1163
Do not use options 'Qbk', 'QBk', or 'QbB' since they scale the coordinates.
1166
if bestoutside is set (e.g., qh_ALL)
1167
returns best facet that is not upperdelaunay
1168
if Delaunay and inside, point is outside circumsphere of bestfacet
1170
returns first facet below point
1171
if point is inside, returns nearest, !upperdelaunay facet
1173
isoutside set if outside of facet
1176
this works for all distributions
1177
if inside, qh_findbestfacet performs an exhaustive search
1178
this may be too conservative. Sometimes it is clearly required.
1179
qh_findbestfacet is not used by qhull.
1180
uses qh.visit_id and qh.coplanarset
1183
<a href="geom.c#findbest">qh_findbest</a>
1185
facetT *qh_findbestfacet (pointT *point, boolT bestoutside,
1186
realT *bestdist, boolT *isoutside) {
1187
facetT *bestfacet= NULL;
1188
int numpart, totpart= 0;
1190
bestfacet= qh_findbest (point, qh facet_list,
1191
bestoutside, !qh_ISnewfacets, bestoutside /* qh_NOupper */,
1192
bestdist, isoutside, &totpart);
1193
if (*bestdist < -qh DISTround) {
1194
bestfacet= qh_findfacet_all (point, bestdist, isoutside, &numpart);
1196
if ((isoutside && bestoutside)
1197
|| (!isoutside && bestfacet->upperdelaunay)) {
1198
bestfacet= qh_findbest (point, bestfacet,
1199
bestoutside, False, bestoutside,
1200
bestdist, isoutside, &totpart);
1204
trace3((qh ferr, "qh_findbestfacet: f%d dist %2.2g isoutside %d totpart %d\n",
1205
bestfacet->id, *bestdist, *isoutside, totpart));
1207
} /* findbestfacet */
1209
/*-<a href="qh-poly.htm#TOC"
1210
>-------------------------------</a><a name="findfacet_all">-</a>
1212
qh_findfacet_all( point, bestdist, isoutside, numpart )
1213
exhaustive search for facet below a point
1215
for Delaunay triangulations,
1216
Use qh_setdelaunay() to lift point to paraboloid and scale by 'Qbb' if needed
1217
Do not use options 'Qbk', 'QBk', or 'QbB' since they scale the coordinates.
1220
returns first facet below point
1222
returns nearest facet
1224
isoutside if point is outside of the hull
1225
number of distance tests
1227
facetT *qh_findfacet_all (pointT *point, realT *bestdist, boolT *isoutside,
1229
facetT *bestfacet= NULL, *facet;
1236
if (facet->flipped || !facet->normal)
1239
qh_distplane (point, facet, &dist);
1240
if (dist > *bestdist) {
1243
if (dist > qh MINoutside) {
1250
trace3((qh ferr, "qh_findfacet_all: f%d dist %2.2g isoutside %d totpart %d\n",
1251
getid_(bestfacet), *bestdist, *isoutside, totpart));
1253
} /* findfacet_all */
1255
/*-<a href="qh-poly.htm#TOC"
1256
>-------------------------------</a><a name="findgood">-</a>
1258
qh_findgood( facetlist, goodhorizon )
1259
identify good facets for qh.PRINTgood
1261
facet includes point as vertex
1262
if !match, returns goodhorizon
1263
inactive if qh.MERGING
1265
facet is visible or coplanar (>0) or not visible (<0)
1267
facet->normal matches threshold
1268
if !goodhorizon and !match,
1269
selects facet with closest angle
1273
number of new, good facets found
1274
determines facet->good
1275
may update qh.GOODclosest
1278
qh_findgood_all further reduces the good region
1282
mark good facets for qh.GOODpoint
1283
mark good facets for qh.GOODthreshold
1285
update qh.GOODclosest
1287
int qh_findgood (facetT *facetlist, int goodhorizon) {
1288
facetT *facet, *bestfacet= NULL;
1289
realT angle, bestangle= REALmax, dist;
1292
FORALLfacet_(facetlist) {
1296
if (qh GOODvertex>0 && !qh MERGING) {
1297
FORALLfacet_(facetlist) {
1298
if (!qh_isvertex (qh GOODvertexp, facet->vertices)) {
1304
if (qh GOODpoint && numgood) {
1305
FORALLfacet_(facetlist) {
1306
if (facet->good && facet->normal) {
1308
qh_distplane (qh GOODpointp, facet, &dist);
1309
if ((qh GOODpoint > 0) ^ (dist > 0.0)) {
1316
if (qh GOODthreshold && (numgood || goodhorizon || qh GOODclosest)) {
1317
FORALLfacet_(facetlist) {
1318
if (facet->good && facet->normal) {
1319
if (!qh_inthresholds (facet->normal, &angle)) {
1322
if (angle < bestangle) {
1329
if (!numgood && (!goodhorizon || qh GOODclosest)) {
1330
if (qh GOODclosest) {
1331
if (qh GOODclosest->visible)
1332
qh GOODclosest= NULL;
1334
qh_inthresholds (qh GOODclosest->normal, &angle);
1335
if (angle < bestangle)
1336
bestfacet= qh GOODclosest;
1339
if (bestfacet && bestfacet != qh GOODclosest) {
1341
qh GOODclosest->good= False;
1342
qh GOODclosest= bestfacet;
1343
bestfacet->good= True;
1345
trace2((qh ferr, "qh_findgood: f%d is closest (%2.2g) to thresholds\n",
1346
bestfacet->id, bestangle));
1349
}else if (qh GOODclosest) { /* numgood > 0 */
1350
qh GOODclosest->good= False;
1351
qh GOODclosest= NULL;
1354
zadd_(Zgoodfacet, numgood);
1355
trace2((qh ferr, "qh_findgood: found %d good facets with %d good horizon\n",
1356
numgood, goodhorizon));
1357
if (!numgood && qh GOODvertex>0 && !qh MERGING)
1362
/*-<a href="qh-poly.htm#TOC"
1363
>-------------------------------</a><a name="findgood_all">-</a>
1365
qh_findgood_all( facetlist )
1366
apply other constraints for good facets (used by qh.PRINTgood)
1368
facet includes (>0) or doesn't include (<0) point as vertex
1369
if last good facet and ONLYgood, prints warning and continues
1370
if qh.SPLITthresholds
1371
facet->normal matches threshold, or if none, the closest one
1373
nop if good not used
1376
clears facet->good if not good
1380
this is like qh_findgood but more restrictive
1383
uses qh_findgood to mark good facets
1384
marks facets for qh.GOODvertex
1385
marks facets for qh.SPLITthreholds
1387
void qh_findgood_all (facetT *facetlist) {
1388
facetT *facet, *bestfacet=NULL;
1389
realT angle, bestangle= REALmax;
1390
int numgood=0, startgood;
1392
if (!qh GOODvertex && !qh GOODthreshold && !qh GOODpoint
1393
&& !qh SPLITthresholds)
1396
qh_findgood (qh facet_list, 0);
1397
FORALLfacet_(facetlist) {
1401
if (qh GOODvertex <0 || (qh GOODvertex > 0 && qh MERGING)) {
1402
FORALLfacet_(facetlist) {
1403
if (facet->good && ((qh GOODvertex > 0) ^ !!qh_isvertex (qh GOODvertexp, facet->vertices))) {
1406
fprintf (qh ferr, "qhull warning: good vertex p%d does not match last good facet f%d. Ignored.\n",
1407
qh_pointid(qh GOODvertexp), facet->id);
1409
}else if (qh GOODvertex > 0)
1410
fprintf (qh ferr, "qhull warning: point p%d is not a vertex ('QV%d').\n",
1411
qh GOODvertex-1, qh GOODvertex-1);
1413
fprintf (qh ferr, "qhull warning: point p%d is a vertex for every facet ('QV-%d').\n",
1414
-qh GOODvertex - 1, -qh GOODvertex - 1);
1421
if (qh SPLITthresholds) {
1422
FORALLfacet_(facetlist) {
1424
if (!qh_inthresholds (facet->normal, &angle)) {
1427
if (angle < bestangle) {
1434
if (!numgood && bestfacet) {
1435
bestfacet->good= True;
1437
trace0((qh ferr, "qh_findgood_all: f%d is closest (%2.2g) to thresholds\n",
1438
bestfacet->id, bestangle));
1442
qh num_good= numgood;
1443
trace0((qh ferr, "qh_findgood_all: %d good facets remain out of %d facets\n",
1444
numgood, startgood));
1445
} /* findgood_all */
1447
/*-<a href="qh-poly.htm#TOC"
1448
>-------------------------------</a><a name="furthestnext">-</a>
1451
set qh.facet_next to facet with furthest of all furthest points
1452
searches all facets on qh.facet_list
1455
this may help avoid precision problems
1457
void qh_furthestnext (void /* qh facet_list */) {
1458
facetT *facet, *bestfacet= NULL;
1459
realT dist, bestdist= -REALmax;
1462
if (facet->outsideset) {
1463
#if qh_COMPUTEfurthest
1465
furthest= (pointT*)qh_setlast (facet->outsideset);
1466
zinc_(Zcomputefurthest);
1467
qh_distplane (furthest, facet, &dist);
1469
dist= facet->furthestdist;
1471
if (dist > bestdist) {
1478
qh_removefacet (bestfacet);
1479
qh_prependfacet (bestfacet, &qh facet_next);
1480
trace1((qh ferr, "qh_furthestnext: made f%d next facet (dist %.2g)\n",
1481
bestfacet->id, bestdist));
1483
} /* furthestnext */
1485
/*-<a href="qh-poly.htm#TOC"
1486
>-------------------------------</a><a name="furthestout">-</a>
1488
qh_furthestout( facet )
1489
make furthest outside point the last point of outsideset
1492
updates facet->outsideset
1493
clears facet->notfurthest
1494
sets facet->furthestdist
1497
determine best point of outsideset
1498
make it the last point of outsideset
1500
void qh_furthestout (facetT *facet) {
1501
pointT *point, **pointp, *bestpoint= NULL;
1502
realT dist, bestdist= -REALmax;
1504
FOREACHpoint_(facet->outsideset) {
1505
qh_distplane (point, facet, &dist);
1506
zinc_(Zcomputefurthest);
1507
if (dist > bestdist) {
1513
qh_setdel (facet->outsideset, point);
1514
qh_setappend (&facet->outsideset, point);
1515
#if !qh_COMPUTEfurthest
1516
facet->furthestdist= bestdist;
1519
facet->notfurthest= False;
1520
trace3((qh ferr, "qh_furthestout: p%d is furthest outside point of f%d\n",
1521
qh_pointid (point), facet->id));
1525
/*-<a href="qh-qhull.htm#TOC"
1526
>-------------------------------</a><a name="infiniteloop">-</a>
1528
qh_infiniteloop( facet )
1529
report infinite loop error due to facet
1531
void qh_infiniteloop (facetT *facet) {
1533
fprintf (qh ferr, "qhull internal error (qh_infiniteloop): potential infinite loop detected\n");
1534
qh_errexit (qh_ERRqhull, facet, NULL);
1535
} /* qh_infiniteloop */
1537
/*-<a href="qh-poly.htm#TOC"
1538
>-------------------------------</a><a name="initbuild">-</a>
1541
initialize hull and outside sets with point array
1542
qh.FIRSTpoint/qh.NUMpoints is point array
1544
adds qh.GOODpoint to initial hull
1547
qh_facetlist with initial hull
1548
points partioned into outside sets, coplanar sets, or inside
1549
initializes qh.GOODpointp, qh.GOODvertexp,
1552
initialize global variables used during qh_buildhull
1553
determine precision constants and points with max/min coordinate values
1554
if qh.SCALElast, scale last coordinate (for 'd')
1555
build initial simplex
1556
partition input points into facets of initial simplex
1560
add qh.GOODvertex if defined
1562
void qh_initbuild( void) {
1563
setT *maxpoints, *vertices;
1571
qh facet_id= qh vertex_id= qh ridge_id= 0;
1572
qh visit_id= qh vertex_visit= 0;
1573
qh maxoutdone= False;
1575
if (qh GOODpoint > 0)
1576
qh GOODpointp= qh_point (qh GOODpoint-1);
1577
else if (qh GOODpoint < 0)
1578
qh GOODpointp= qh_point (-qh GOODpoint-1);
1579
if (qh GOODvertex > 0)
1580
qh GOODvertexp= qh_point (qh GOODvertex-1);
1581
else if (qh GOODvertex < 0)
1582
qh GOODvertexp= qh_point (-qh GOODvertex-1);
1584
&& (qh GOODpointp < qh first_point /* also catches !GOODpointp */
1585
|| qh GOODpointp > qh_point (qh num_points-1)))
1587
&& (qh GOODvertexp < qh first_point /* also catches !GOODvertexp */
1588
|| qh GOODvertexp > qh_point (qh num_points-1)))) {
1589
fprintf (qh ferr, "qhull input error: either QGn or QVn point is > p%d\n",
1591
qh_errexit (qh_ERRinput, NULL, NULL);
1593
maxpoints= qh_maxmin(qh first_point, qh num_points, qh hull_dim);
1595
qh_scalelast (qh first_point, qh num_points, qh hull_dim,
1596
qh MINlastcoord, qh MAXlastcoord, qh MAXwidth);
1598
if (qh DELAUNAY && qh upper_threshold[qh hull_dim-1] > REALmax/2
1599
&& qh lower_threshold[qh hull_dim-1] < -REALmax/2) {
1600
for (i= qh_PRINTEND; i--; ) {
1601
if (qh PRINTout[i] == qh_PRINTgeom && qh DROPdim < 0
1602
&& !qh GOODthreshold && !qh SPLITthresholds)
1603
break; /* in this case, don't set upper_threshold */
1606
if (qh UPPERdelaunay) { /* matches qh.upperdelaunay in qh_setfacetplane */
1607
qh lower_threshold[qh hull_dim-1]= qh ANGLEround * qh_ZEROdelaunay;
1608
qh GOODthreshold= True;
1610
qh upper_threshold[qh hull_dim-1]= -qh ANGLEround * qh_ZEROdelaunay;
1611
if (!qh GOODthreshold)
1612
qh SPLITthresholds= True; /* build upper-convex hull even if Qg */
1613
/* qh_initqhull_globals errors if Qg without Pdk/etc. */
1617
vertices= qh_initialvertices(qh hull_dim, maxpoints, qh first_point, qh num_points);
1618
qh_initialhull (vertices); /* initial qh facet_list */
1619
qh_partitionall (vertices, qh first_point, qh num_points);
1620
if (qh PRINToptions1st || qh TRACElevel || qh IStracing) {
1621
if (qh TRACElevel || qh IStracing)
1622
fprintf (qh ferr, "\nTrace level %d for %s | %s\n",
1623
qh IStracing ? qh IStracing : qh TRACElevel, qh rbox_command, qh qhull_command);
1624
fprintf (qh ferr, "Options selected for Qhull %s:\n%s\n", qh_VERSION, qh qhull_options);
1626
qh_resetlists (False, qh_RESETvisible /*qh visible_list newvertex_list newfacet_list */);
1627
qh facet_next= qh facet_list;
1628
qh_furthestnext (/* qh facet_list */);
1630
qh cos_max= qh premerge_cos;
1631
qh centrum_radius= qh premerge_centrum;
1634
if (qh GOODvertex > 0 && qh MERGING) {
1635
fprintf (qh ferr, "qhull input error: 'Qg QVn' (only good vertex) does not work with merging.\nUse 'QJ' to joggle the input or 'Q0' to turn off merging.\n");
1636
qh_errexit (qh_ERRinput, NULL, NULL);
1638
if (!(qh GOODthreshold || qh GOODpoint
1639
|| (!qh MERGEexact && !qh PREmerge && qh GOODvertexp))) {
1640
fprintf (qh ferr, "qhull input error: 'Qg' (ONLYgood) needs a good threshold ('Pd0D0'), a\n\
1641
good point (QGn or QG-n), or a good vertex with 'QJ' or 'Q0' (QVn).\n");
1642
qh_errexit (qh_ERRinput, NULL, NULL);
1644
if (qh GOODvertex > 0 && !qh MERGING /* matches qh_partitionall */
1645
&& !qh_isvertex (qh GOODvertexp, vertices)) {
1646
facet= qh_findbestnew (qh GOODvertexp, qh facet_list,
1647
&dist, !qh_ALL, &isoutside, &numpart);
1648
zadd_(Zdistgood, numpart);
1650
fprintf (qh ferr, "qhull input error: point for QV%d is inside initial simplex. It can not be made a vertex.\n",
1651
qh_pointid(qh GOODvertexp));
1652
qh_errexit (qh_ERRinput, NULL, NULL);
1654
if (!qh_addpoint (qh GOODvertexp, facet, False)) {
1655
qh_settempfree(&vertices);
1656
qh_settempfree(&maxpoints);
1660
qh_findgood (qh facet_list, 0);
1662
qh_settempfree(&vertices);
1663
qh_settempfree(&maxpoints);
1664
trace1((qh ferr, "qh_initbuild: initial hull created and points partitioned\n"));
1667
/*-<a href="qh-poly.htm#TOC"
1668
>-------------------------------</a><a name="initialhull">-</a>
1670
qh_initialhull( vertices )
1671
constructs the initial hull as a DIM3 simplex of vertices
1674
creates a simplex (initializes lists)
1675
determines orientation of simplex
1676
sets hyperplanes for facets
1677
doubles checks orientation (in case of axis-parallel facets with Gaussian elimination)
1678
checks for flipped facets and qh.NARROWhull
1681
void qh_initialhull(setT *vertices) {
1682
facetT *facet, *firstfacet, *neighbor, **neighborp;
1683
realT dist, angle, minangle= REALmax;
1688
qh_createsimplex(vertices); /* qh facet_list */
1689
qh_resetlists (False, qh_RESETvisible);
1690
qh facet_next= qh facet_list; /* advance facet when processed */
1691
qh interior_point= qh_getcenter(vertices);
1692
firstfacet= qh facet_list;
1693
qh_setfacetplane(firstfacet);
1694
zinc_(Znumvisibility); /* needs to be in printsummary */
1695
qh_distplane(qh interior_point, firstfacet, &dist);
1698
facet->toporient ^= True;
1701
qh_setfacetplane(facet);
1703
if (!qh_checkflipped (facet, NULL, qh_ALL)) {/* due to axis-parallel facet */
1704
trace1((qh ferr, "qh_initialhull: initial orientation incorrect. Correct all facets\n"));
1705
facet->flipped= False;
1707
facet->toporient ^= True;
1708
qh_orientoutside (facet);
1714
if (!qh_checkflipped (facet, NULL, !qh_ALL)) { /* can happen with 'R0.1' */
1715
qh_precision ("initial facet is coplanar with interior point");
1716
fprintf (qh ferr, "qhull precision error: initial facet %d is coplanar with the interior point\n",
1718
qh_errexit (qh_ERRsingular, facet, NULL);
1720
FOREACHneighbor_(facet) {
1721
angle= qh_getangle (facet->normal, neighbor->normal);
1722
minimize_( minangle, angle);
1725
if (minangle < qh_MAXnarrow && !qh NOnarrow) {
1726
realT diff= 1.0 + minangle;
1728
qh NARROWhull= True;
1729
qh_option ("_narrow-hull", NULL, &diff);
1730
if (minangle < qh_WARNnarrow && !qh RERUN && qh PRINTprecision)
1731
fprintf (qh ferr, "qhull precision warning: \n\
1732
The initial hull is narrow (cosine of min. angle is %.16f).\n\
1733
A coplanar point may lead to a wide facet. Options 'QbB' (scale to unit box)\n\
1734
or 'Qbb' (scale last coordinate) may remove this warning. Use 'Pp' to skip\n\
1735
this warning. See 'Limitations' in qh-impre.htm.\n",
1736
-minangle); /* convert from angle between normals to angle between facets */
1738
zzval_(Zprocessed)= qh hull_dim+1;
1739
qh_checkpolygon (qh facet_list);
1740
qh_checkconvex(qh facet_list, qh_DATAfault);
1742
if (qh IStracing >= 1) {
1743
fprintf(qh ferr, "qh_initialhull: simplex constructed, interior point:");
1744
for (k=0; k < qh hull_dim; k++)
1745
fprintf (qh ferr, " %6.4g", qh interior_point[k]);
1746
fprintf (qh ferr, "\n");
1751
/*-<a href="qh-poly.htm#TOC"
1752
>-------------------------------</a><a name="initialvertices">-</a>
1754
qh_initialvertices( dim, maxpoints, points, numpoints )
1755
determines a non-singular set of initial vertices
1756
maxpoints may include duplicate points
1759
temporary set of dim+1 vertices in descending order by vertex id
1760
if qh.RANDOMoutside && !qh.ALLpoints
1762
if dim >= qh_INITIALmax,
1763
uses min/max x and max points with non-zero determinants
1766
unless qh.ALLpoints,
1767
uses maxpoints as long as determinate is non-zero
1769
setT *qh_initialvertices(int dim, setT *maxpoints, pointT *points, int numpoints) {
1770
pointT *point, **pointp;
1771
setT *vertices, *simplex, *tested;
1773
int index, point_i, point_n, k;
1774
boolT nearzero= False;
1776
vertices= qh_settemp (dim + 1);
1777
simplex= qh_settemp (dim+1);
1779
qh_maxsimplex (dim, NULL, points, numpoints, &simplex);
1780
else if (qh RANDOMoutside) {
1781
while (qh_setsize (simplex) != dim+1) {
1782
randr= qh_RANDOMint;
1783
randr= randr/(qh_RANDOMmax+1);
1784
index= (int)floor(qh num_points * randr);
1785
while (qh_setin (simplex, qh_point (index))) {
1786
index++; /* in case qh_RANDOMint always returns the same value */
1787
index= index < qh num_points ? index : 0;
1789
qh_setappend (&simplex, qh_point (index));
1791
}else if (qh hull_dim >= qh_INITIALmax) {
1792
tested= qh_settemp (dim+1);
1793
qh_setappend (&simplex, SETfirst_(maxpoints)); /* max and min X coord */
1794
qh_setappend (&simplex, SETsecond_(maxpoints));
1795
qh_maxsimplex (fmin_(qh_INITIALsearch, dim), maxpoints, points, numpoints, &simplex);
1796
k= qh_setsize (simplex);
1797
FOREACHpoint_i_(maxpoints) {
1798
if (point_i & 0x1) { /* first pick up max. coord. points */
1799
if (!qh_setin (simplex, point) && !qh_setin (tested, point)){
1800
qh_detsimplex(point, simplex, k, &nearzero);
1802
qh_setappend (&tested, point);
1804
qh_setappend (&simplex, point);
1805
if (++k == dim) /* use search for last point */
1811
while (k != dim && (point= (pointT*)qh_setdellast (maxpoints))) {
1812
if (!qh_setin (simplex, point) && !qh_setin (tested, point)){
1813
qh_detsimplex (point, simplex, k, &nearzero);
1815
qh_setappend (&tested, point);
1817
qh_setappend (&simplex, point);
1823
while (k != dim && (point= qh_point (index++))) {
1824
if (!qh_setin (simplex, point) && !qh_setin (tested, point)){
1825
qh_detsimplex (point, simplex, k, &nearzero);
1827
qh_setappend (&simplex, point);
1832
qh_settempfree (&tested);
1833
qh_maxsimplex (dim, maxpoints, points, numpoints, &simplex);
1835
qh_maxsimplex (dim, maxpoints, points, numpoints, &simplex);
1836
FOREACHpoint_(simplex)
1837
qh_setaddnth (&vertices, 0, qh_newvertex(point)); /* descending order */
1838
qh_settempfree (&simplex);
1840
} /* initialvertices */
1843
/*-<a href="qh-poly.htm#TOC"
1844
>-------------------------------</a><a name="isvertex">-</a>
1847
returns vertex if point is in vertex set, else returns NULL
1852
vertexT *qh_isvertex (pointT *point, setT *vertices) {
1853
vertexT *vertex, **vertexp;
1855
FOREACHvertex_(vertices) {
1856
if (vertex->point == point)
1862
/*-<a href="qh-poly.htm#TOC"
1863
>-------------------------------</a><a name="makenewfacets">-</a>
1865
qh_makenewfacets( point )
1866
make new facets from point and qh.visible_list
1869
qh.newfacet_list= list of new facets with hyperplanes and ->newfacet
1870
qh.newvertex_list= list of vertices in new facets with ->newlist set
1873
newfacets reference horizon facets, but not vice versa
1874
ridges reference non-simplicial horizon ridges, but not vice versa
1875
does not change existing facets
1878
new facets attached to horizon facets and ridges
1880
visible->r.replace is corresponding new facet
1883
qh_makenewplanes() -- make hyperplanes for facets
1884
qh_attachnewfacets() -- attachnewfacets if not done here (qh ONLYgood)
1885
qh_matchnewfacets() -- match up neighbors
1886
qh_updatevertices() -- update vertex neighbors and delvertices
1887
qh_deletevisible() -- delete visible facets
1888
qh_checkpolygon() --check the result
1889
qh_triangulate() -- triangulate a non-simplicial facet
1892
for each visible facet
1893
make new facets to its horizon facets
1894
update its f.replace
1895
clear its neighbor set
1897
vertexT *qh_makenewfacets (pointT *point /*visible_list*/) {
1898
facetT *visible, *newfacet= NULL, *newfacet2= NULL, *neighbor, **neighborp;
1902
qh newfacet_list= qh facet_tail;
1903
qh newvertex_list= qh vertex_tail;
1904
apex= qh_newvertex(point);
1905
qh_appendvertex (apex);
1909
FORALLvisible_facets {
1910
FOREACHneighbor_(visible)
1911
neighbor->seen= False;
1912
if (visible->ridges) {
1913
visible->visitid= qh visit_id;
1914
newfacet2= qh_makenew_nonsimplicial (visible, apex, &numnew);
1916
if (visible->simplicial)
1917
newfacet= qh_makenew_simplicial (visible, apex, &numnew);
1919
if (newfacet2) /* newfacet is null if all ridges defined */
1920
newfacet= newfacet2;
1922
visible->f.replace= newfacet;
1924
zinc_(Zinsidevisible);
1925
SETfirst_(visible->neighbors)= NULL;
1928
trace1((qh ferr, "qh_makenewfacets: created %d new facets from point p%d to horizon\n",
1929
numnew, qh_pointid(point)));
1930
if (qh IStracing >= 4)
1931
qh_printfacetlist (qh newfacet_list, NULL, qh_ALL);
1933
} /* makenewfacets */
1935
/*-<a href="qh-poly.htm#TOC"
1936
>-------------------------------</a><a name="matchduplicates">-</a>
1938
qh_matchduplicates( atfacet, atskip, hashsize, hashcount )
1939
match duplicate ridges in qh.hash_table for atfacet/atskip
1940
duplicates marked with ->dupridge and qh_DUPLICATEridge
1943
picks match with worst merge (min distance apart)
1952
compute hash value for atfacet and atskip
1953
repeat twice -- once to make best matches, once to match the rest
1954
for each possible facet in qh.hash_table
1955
if it is a matching facet and pass 2
1957
unless tricoplanar, mark match for merging (qh_MERGEridge)
1958
[e.g., tricoplanar RBOX s 1000 t993602376 | QHULL C-1e-3 d Qbb FA Qt]
1959
if it is a matching facet and pass 1
1960
test if this is a better match
1962
make best match (it will not be merged)
1965
void qh_matchduplicates (facetT *atfacet, int atskip, int hashsize, int *hashcount) {
1966
boolT same, ismatch;
1968
facetT *facet, *newfacet, *maxmatch= NULL, *maxmatch2= NULL, *nextfacet;
1969
int skip, newskip, nextskip= 0, maxskip= 0, maxskip2= 0, makematch;
1970
realT maxdist= -REALmax, mindist, dist2, low, high;
1972
hash= (int)qh_gethash (hashsize, atfacet->vertices, qh hull_dim, 1,
1973
SETelem_(atfacet->vertices, atskip));
1974
trace2((qh ferr, "qh_matchduplicates: find duplicate matches for f%d skip %d hash %d hashcount %d\n",
1975
atfacet->id, atskip, hash, *hashcount));
1976
for (makematch= 0; makematch < 2; makematch++) {
1978
for (newfacet= atfacet, newskip= atskip; newfacet; newfacet= nextfacet, newskip= nextskip) {
1981
newfacet->visitid= qh visit_id;
1982
for (scan= hash; (facet= SETelemt_(qh hash_table, scan, facetT));
1983
scan= (++scan >= hashsize ? 0 : scan)) {
1984
if (!facet->dupridge || facet->visitid == qh visit_id)
1987
if (qh_matchvertices (1, newfacet->vertices, newskip, facet->vertices, &skip, &same)) {
1988
ismatch= (same == (newfacet->toporient ^ facet->toporient));
1989
if (SETelemt_(facet->neighbors, skip, facetT) != qh_DUPLICATEridge) {
1991
fprintf (qh ferr, "qhull internal error (qh_matchduplicates): missing dupridge at f%d skip %d for new f%d skip %d hash %d\n",
1992
facet->id, skip, newfacet->id, newskip, hash);
1993
qh_errexit2 (qh_ERRqhull, facet, newfacet);
1995
}else if (ismatch && makematch) {
1996
if (SETelemt_(newfacet->neighbors, newskip, facetT) == qh_DUPLICATEridge) {
1997
SETelem_(facet->neighbors, skip)= newfacet;
1998
if (newfacet->tricoplanar)
1999
SETelem_(newfacet->neighbors, newskip)= facet;
2001
SETelem_(newfacet->neighbors, newskip)= qh_MERGEridge;
2002
*hashcount -= 2; /* removed two unmatched facets */
2003
trace4((qh ferr, "qh_matchduplicates: duplicate f%d skip %d matched with new f%d skip %d merge\n",
2004
facet->id, skip, newfacet->id, newskip));
2006
}else if (ismatch) {
2007
mindist= qh_getdistance (facet, newfacet, &low, &high);
2008
dist2= qh_getdistance (newfacet, facet, &low, &high);
2009
minimize_(mindist, dist2);
2010
if (mindist > maxdist) {
2014
maxmatch2= newfacet;
2017
trace3((qh ferr, "qh_matchduplicates: duplicate f%d skip %d new f%d skip %d at dist %2.2g, max is now f%d f%d\n",
2018
facet->id, skip, newfacet->id, newskip, mindist,
2019
maxmatch->id, maxmatch2->id));
2020
}else { /* !ismatch */
2025
if (makematch && !facet
2026
&& SETelemt_(facet->neighbors, skip, facetT) == qh_DUPLICATEridge) {
2027
fprintf (qh ferr, "qhull internal error (qh_matchduplicates): no MERGEridge match for duplicate f%d skip %d at hash %d\n",
2028
newfacet->id, newskip, hash);
2029
qh_errexit (qh_ERRqhull, newfacet, NULL);
2032
} /* end of for each new facet at hash */
2035
fprintf (qh ferr, "qhull internal error (qh_matchduplicates): no maximum match at duplicate f%d skip %d at hash %d\n",
2036
atfacet->id, atskip, hash);
2037
qh_errexit (qh_ERRqhull, atfacet, NULL);
2039
SETelem_(maxmatch->neighbors, maxskip)= maxmatch2;
2040
SETelem_(maxmatch2->neighbors, maxskip2)= maxmatch;
2041
*hashcount -= 2; /* removed two unmatched facets */
2042
zzinc_(Zmultiridge);
2043
trace0((qh ferr, "qh_matchduplicates: duplicate f%d skip %d matched with new f%d skip %d keep\n",
2044
maxmatch->id, maxskip, maxmatch2->id, maxskip2));
2045
qh_precision ("ridge with multiple neighbors");
2046
if (qh IStracing >= 4)
2047
qh_errprint ("DUPLICATED/MATCH", maxmatch, maxmatch2, NULL, NULL);
2050
} /* matchduplicates */
2052
/*-<a href="qh-poly.htm#TOC"
2053
>-------------------------------</a><a name="nearcoplanar">-</a>
2056
for all facets, remove near-inside points from facet->coplanarset</li>
2057
coplanar points defined by innerplane from qh_outerinner()
2060
if qh KEEPcoplanar && !qh KEEPinside
2061
facet->coplanarset only contains coplanar points
2063
drops inner plane by another qh.JOGGLEmax diagonal since a
2064
vertex could shift out while a coplanar point shifts in
2067
used for qh.PREmerge and qh.JOGGLEmax
2068
must agree with computation of qh.NEARcoplanar in qh_detroundoff()
2070
if not keeping coplanar or inside points
2071
free all coplanar sets
2072
else if not keeping both coplanar and inside points
2073
remove !coplanar or !inside points from coplanar sets
2075
void qh_nearcoplanar ( void /* qh.facet_list */) {
2077
pointT *point, **pointp;
2079
realT dist, innerplane;
2081
if (!qh KEEPcoplanar && !qh KEEPinside) {
2083
if (facet->coplanarset)
2084
qh_setfree( &facet->coplanarset);
2086
}else if (!qh KEEPcoplanar || !qh KEEPinside) {
2087
qh_outerinner (NULL, NULL, &innerplane);
2088
if (qh JOGGLEmax < REALmax/2)
2089
innerplane -= qh JOGGLEmax * sqrt (qh hull_dim);
2092
if (facet->coplanarset) {
2093
FOREACHpoint_(facet->coplanarset) {
2095
qh_distplane (point, facet, &dist);
2096
if (dist < innerplane) {
2098
SETref_(point)= NULL;
2099
}else if (!qh KEEPcoplanar)
2100
SETref_(point)= NULL;
2102
qh_setcompact (facet->coplanarset);
2105
zzadd_(Zcheckpart, numpart);
2107
} /* nearcoplanar */
2109
/*-<a href="qh-poly.htm#TOC"
2110
>-------------------------------</a><a name="nearvertex">-</a>
2112
qh_nearvertex( facet, point, bestdist )
2113
return nearest vertex in facet to point
2116
vertex and its distance
2120
distance is measured in the input set
2121
searches neighboring tricoplanar facets (requires vertexneighbors)
2122
Slow implementation. Recomputes vertex set for each point.
2123
The vertex set could be stored in the qh.keepcentrum facet.
2125
vertexT *qh_nearvertex (facetT *facet, pointT *point, realT *bestdistp) {
2126
realT bestdist= REALmax, dist;
2127
vertexT *bestvertex= NULL, *vertex, **vertexp, *apex;
2129
facetT *neighbor, **neighborp;
2131
int dim= qh hull_dim;
2135
if (facet->tricoplanar) {
2136
if (!qh VERTEXneighbors || !facet->center) {
2137
fprintf(qh ferr, "qhull internal error (qh_nearvertex): qh.VERTEXneighbors and facet->center required for tricoplanar facets\n");
2138
qh_errexit(qh_ERRqhull, NULL, NULL);
2140
vertices= qh_settemp (qh TEMPsize);
2141
apex= SETfirst_(facet->vertices);
2142
center= facet->center;
2143
FOREACHneighbor_(apex) {
2144
if (neighbor->center == center) {
2145
FOREACHvertex_(neighbor->vertices)
2146
qh_setappend(&vertices, vertex);
2150
vertices= facet->vertices;
2151
FOREACHvertex_(vertices) {
2152
dist= qh_pointdist (vertex->point, point, -dim);
2153
if (dist < bestdist) {
2158
if (facet->tricoplanar)
2159
qh_settempfree (&vertices);
2160
*bestdistp= sqrt (bestdist);
2164
/*-<a href="qh-poly.htm#TOC"
2165
>-------------------------------</a><a name="newhashtable">-</a>
2167
qh_newhashtable( newsize )
2168
returns size of qh.hash_table of at least newsize slots
2171
assumes qh.hash_table is NULL
2172
qh_HASHfactor determines the number of extra slots
2173
size is not divisible by 2, 3, or 5
2175
int qh_newhashtable(int newsize) {
2178
size= ((newsize+1)*qh_HASHfactor) | 0x1; /* odd number */
2180
if ((size%3) && (size%5))
2183
/* loop terminates because there is an infinite number of primes */
2185
qh hash_table= qh_setnew (size);
2186
qh_setzero (qh hash_table, 0, size);
2188
} /* newhashtable */
2190
/*-<a href="qh-poly.htm#TOC"
2191
>-------------------------------</a><a name="newvertex">-</a>
2193
qh_newvertex( point )
2194
returns a new vertex for point
2196
vertexT *qh_newvertex(pointT *point) {
2199
zinc_(Ztotvertices);
2200
vertex= (vertexT *)qh_memalloc(sizeof(vertexT));
2201
memset ((char *) vertex, 0, sizeof (vertexT));
2202
if (qh vertex_id == 0xFFFFFF) {
2203
fprintf(qh ferr, "qhull input error: more than %d vertices. ID field overflows and two vertices\n\
2204
may have the same identifier. Vertices not sorted correctly.\n", 0xFFFFFF);
2205
qh_errexit(qh_ERRinput, NULL, NULL);
2207
if (qh vertex_id == qh tracevertex_id)
2208
qh tracevertex= vertex;
2209
vertex->id= qh vertex_id++;
2210
vertex->point= point;
2211
trace4((qh ferr, "qh_newvertex: vertex p%d (v%d) created\n", qh_pointid(vertex->point),
2216
/*-<a href="qh-poly.htm#TOC"
2217
>-------------------------------</a><a name="nextridge3d">-</a>
2219
qh_nextridge3d( atridge, facet, vertex )
2220
return next ridge and vertex for a 3d facet
2223
in qh_ORIENTclock order
2224
this is a O(n^2) implementation to trace all ridges
2225
be sure to stop on any 2nd visit
2229
exit if it is the ridge after atridge
2231
ridgeT *qh_nextridge3d (ridgeT *atridge, facetT *facet, vertexT **vertexp) {
2232
vertexT *atvertex, *vertex, *othervertex;
2233
ridgeT *ridge, **ridgep;
2235
if ((atridge->top == facet) ^ qh_ORIENTclock)
2236
atvertex= SETsecondt_(atridge->vertices, vertexT);
2238
atvertex= SETfirstt_(atridge->vertices, vertexT);
2239
FOREACHridge_(facet->ridges) {
2240
if (ridge == atridge)
2242
if ((ridge->top == facet) ^ qh_ORIENTclock) {
2243
othervertex= SETsecondt_(ridge->vertices, vertexT);
2244
vertex= SETfirstt_(ridge->vertices, vertexT);
2246
vertex= SETsecondt_(ridge->vertices, vertexT);
2247
othervertex= SETfirstt_(ridge->vertices, vertexT);
2249
if (vertex == atvertex) {
2251
*vertexp= othervertex;
2257
#else /* qh_NOmerge */
2258
void qh_matchduplicates (facetT *atfacet, int atskip, int hashsize, int *hashcount) {
2260
ridgeT *qh_nextridge3d (ridgeT *atridge, facetT *facet, vertexT **vertexp) {
2264
#endif /* qh_NOmerge */
2266
/*-<a href="qh-poly.htm#TOC"
2267
>-------------------------------</a><a name="outcoplanar">-</a>
2270
move points from all facets' outsidesets to their coplanarsets
2273
for post-processing under qh.NARROWhull
2277
for each outside point for facet
2278
partition point into coplanar set
2280
void qh_outcoplanar (void /* facet_list */) {
2281
pointT *point, **pointp;
2285
trace1((qh ferr, "qh_outcoplanar: move outsideset to coplanarset for qh NARROWhull\n"));
2287
FOREACHpoint_(facet->outsideset) {
2289
if (qh KEEPcoplanar || qh KEEPnearinside) {
2290
qh_distplane (point, facet, &dist);
2292
qh_partitioncoplanar (point, facet, &dist);
2295
qh_setfree (&facet->outsideset);
2299
/*-<a href="qh-poly.htm#TOC"
2300
>-------------------------------</a><a name="point">-</a>
2303
return point for a point id, or NULL if unknown
2306
return ((pointT *)((unsigned long)qh.first_point
2307
+ (unsigned long)((id)*qh.normal_size)));
2309
pointT *qh_point (int id) {
2313
if (id < qh num_points)
2314
return qh first_point + id * qh hull_dim;
2315
id -= qh num_points;
2316
if (id < qh_setsize (qh other_points))
2317
return SETelemt_(qh other_points, id, pointT);
2321
/*-<a href="qh-poly.htm#TOC"
2322
>-------------------------------</a><a name="point_add">-</a>
2324
qh_point_add( set, point, elem )
2325
stores elem at set[point.id]
2328
access function for qh_pointfacet and qh_pointvertex
2333
void qh_point_add (setT *set, pointT *point, void *elem) {
2336
SETreturnsize_(set, size);
2337
if ((id= qh_pointid(point)) < 0)
2338
fprintf (qh ferr, "qhull internal warning (point_add): unknown point %p id %d\n",
2340
else if (id >= size) {
2341
fprintf (qh ferr, "qhull internal errror (point_add): point p%d is out of bounds (%d)\n",
2343
qh_errexit (qh_ERRqhull, NULL, NULL);
2345
SETelem_(set, id)= elem;
2349
/*-<a href="qh-poly.htm#TOC"
2350
>-------------------------------</a><a name="pointfacet">-</a>
2353
return temporary set of facet for each point
2354
the set is indexed by point id
2357
vertices assigned to one of the facets
2358
coplanarset assigned to the facet
2359
outside set assigned to the facet
2360
NULL if no facet for point (inside)
2361
includes qh.GOODpointp
2364
FOREACHfacet_i_(facets) { ... }
2370
add each coplanar point
2371
add each outside point
2373
setT *qh_pointfacet (void /*qh facet_list*/) {
2374
int numpoints= qh num_points + qh_setsize (qh other_points);
2377
vertexT *vertex, **vertexp;
2378
pointT *point, **pointp;
2380
facets= qh_settemp (numpoints);
2381
qh_setzero (facets, 0, numpoints);
2384
FOREACHvertex_(facet->vertices) {
2385
if (vertex->visitid != qh vertex_visit) {
2386
vertex->visitid= qh vertex_visit;
2387
qh_point_add (facets, vertex->point, facet);
2390
FOREACHpoint_(facet->coplanarset)
2391
qh_point_add (facets, point, facet);
2392
FOREACHpoint_(facet->outsideset)
2393
qh_point_add (facets, point, facet);
2398
/*-<a href="qh-poly.htm#TOC"
2399
>-------------------------------</a><a name="pointvertex">-</a>
2402
return temporary set of vertices indexed by point id
2403
entry is NULL if no vertex for a point
2404
this will include qh.GOODpointp
2407
FOREACHvertex_i_(vertices) { ... }
2408
SETelem_(vertices, i)
2410
setT *qh_pointvertex (void /*qh facet_list*/) {
2411
int numpoints= qh num_points + qh_setsize (qh other_points);
2415
vertices= qh_settemp (numpoints);
2416
qh_setzero (vertices, 0, numpoints);
2418
qh_point_add (vertices, vertex->point, vertex);
2423
/*-<a href="qh-poly.htm#TOC"
2424
>-------------------------------</a><a name="prependfacet">-</a>
2426
qh_prependfacet( facet, facetlist )
2427
prepend facet to the start of a facetlist
2430
increments qh.numfacets
2431
updates facetlist, qh.facet_list, facet_next
2434
be careful of prepending since it can lose a pointer.
2435
e.g., can lose _next by deleting and then prepending before _next
2437
void qh_prependfacet(facetT *facet, facetT **facetlist) {
2438
facetT *prevfacet, *list;
2441
trace4((qh ferr, "qh_prependfacet: prepend f%d before f%d\n",
2442
facet->id, getid_(*facetlist)));
2444
(*facetlist)= qh facet_tail;
2446
prevfacet= list->previous;
2447
facet->previous= prevfacet;
2449
prevfacet->next= facet;
2450
list->previous= facet;
2451
facet->next= *facetlist;
2452
if (qh facet_list == list) /* this may change *facetlist */
2453
qh facet_list= facet;
2454
if (qh facet_next == list)
2455
qh facet_next= facet;
2458
} /* prependfacet */
2461
/*-<a href="qh-poly.htm#TOC"
2462
>-------------------------------</a><a name="printhashtable">-</a>
2464
qh_printhashtable( fp )
2465
print hash table to fp
2468
not in I/O to avoid bringing io.c in
2473
if unmatched or will merge (NULL, qh_MERGEridge, qh_DUPLICATEridge)
2474
print entry and neighbors
2476
void qh_printhashtable(FILE *fp) {
2477
facetT *facet, *neighbor;
2478
int id, facet_i, facet_n, neighbor_i= 0, neighbor_n= 0;
2479
vertexT *vertex, **vertexp;
2481
FOREACHfacet_i_(qh hash_table) {
2483
FOREACHneighbor_i_(facet) {
2484
if (!neighbor || neighbor == qh_MERGEridge || neighbor == qh_DUPLICATEridge)
2487
if (neighbor_i == neighbor_n)
2489
fprintf (fp, "hash %d f%d ", facet_i, facet->id);
2490
FOREACHvertex_(facet->vertices)
2491
fprintf (fp, "v%d ", vertex->id);
2492
fprintf (fp, "\n neighbors:");
2493
FOREACHneighbor_i_(facet) {
2494
if (neighbor == qh_MERGEridge)
2496
else if (neighbor == qh_DUPLICATEridge)
2499
id= getid_(neighbor);
2500
fprintf (fp, " %d", id);
2505
} /* printhashtable */
2508
/*-<a href="qh-poly.htm#TOC"
2509
>-------------------------------</a><a name="printlists">-</a>
2512
print out facet and vertex list for debugging (without 'f/v' tags)
2514
void qh_printlists (void) {
2519
fprintf (qh ferr, "qh_printlists: facets:");
2521
if (++count % 100 == 0)
2522
fprintf (qh ferr, "\n ");
2523
fprintf (qh ferr, " %d", facet->id);
2525
fprintf (qh ferr, "\n new facets %d visible facets %d next facet for qh_addpoint %d\n vertices (new %d):",
2526
getid_(qh newfacet_list), getid_(qh visible_list), getid_(qh facet_next),
2527
getid_(qh newvertex_list));
2530
if (++count % 100 == 0)
2531
fprintf (qh ferr, "\n ");
2532
fprintf (qh ferr, " %d", vertex->id);
2534
fprintf (qh ferr, "\n");
2537
/*-<a href="qh-poly.htm#TOC"
2538
>-------------------------------</a><a name="resetlists">-</a>
2540
qh_resetlists( stats, qh_RESETvisible )
2541
reset newvertex_list, newfacet_list, visible_list
2543
maintains statistics
2546
visible_list is empty if qh_deletevisible was called
2548
void qh_resetlists (boolT stats, boolT resetVisible /*qh newvertex_list newfacet_list visible_list*/) {
2550
facetT *newfacet, *visible;
2551
int totnew=0, totver=0;
2554
FORALLvertex_(qh newvertex_list)
2558
zadd_(Zvisvertextot, totver);
2559
zmax_(Zvisvertexmax, totver);
2560
zadd_(Znewfacettot, totnew);
2561
zmax_(Znewfacetmax, totnew);
2563
FORALLvertex_(qh newvertex_list)
2564
vertex->newlist= False;
2565
qh newvertex_list= NULL;
2567
newfacet->newfacet= False;
2568
qh newfacet_list= NULL;
2570
FORALLvisible_facets {
2571
visible->f.replace= NULL;
2572
visible->visible= False;
2576
qh visible_list= NULL; /* may still have visible facets via qh_triangulate */
2577
qh NEWfacets= False;
2580
/*-<a href="qh-poly.htm#TOC"
2581
>-------------------------------</a><a name="setvoronoi_all">-</a>
2584
compute Voronoi centers for all facets
2585
includes upperDelaunay facets if qh.UPPERdelaunay ('Qu')
2588
facet->center is the Voronoi center
2591
this is unused/untested code
2592
please email bradb@shore.net if this works ok for you
2595
FORALLvertices {...} to locate the vertex for a point.
2596
FOREACHneighbor_(vertex) {...} to visit the Voronoi centers for a Voronoi cell.
2598
void qh_setvoronoi_all (void) {
2601
qh_clearcenters (qh_ASvoronoi);
2602
qh_vertexneighbors();
2605
if (!facet->normal || !facet->upperdelaunay || qh UPPERdelaunay) {
2607
facet->center= qh_facetcenter (facet->vertices);
2610
} /* setvoronoi_all */
2614
/*-<a href="qh-poly.htm#TOC"
2615
>-------------------------------</a><a name="triangulate">-</a>
2618
triangulate non-simplicial facets on qh.facet_list,
2619
if qh.CENTERtype=qh_ASvoronoi, sets Voronoi centers of non-simplicial facets
2622
all facets simplicial
2623
each tricoplanar facet has ->f.triowner == owner of ->center,normal,etc.
2626
call after qh_check_output since may switch to Voronoi centers
2627
Output may overwrite ->f.triowner with ->f.area
2629
void qh_triangulate (void /*qh facet_list*/) {
2630
facetT *facet, *nextfacet, *owner;
2631
int onlygood= qh ONLYgood;
2632
facetT *neighbor, *visible= NULL, *facet1, *facet2, *new_facet_list= NULL;
2633
facetT *orig_neighbor= NULL, *otherfacet;
2634
vertexT *new_vertex_list= NULL;
2636
mergeType mergetype;
2637
int neighbor_i, neighbor_n;
2639
trace1((qh ferr, "qh_triangulate: triangulate non-simplicial facets\n"));
2640
if (qh hull_dim == 2)
2642
if (qh VORONOI) { /* otherwise lose Voronoi centers [could rebuild vertex set from tricoplanar] */
2643
qh_clearcenters (qh_ASvoronoi);
2644
qh_vertexneighbors();
2646
qh ONLYgood= False; /* for makenew_nonsimplicial */
2649
qh degen_mergeset= qh_settemp (qh TEMPsize);
2650
qh newvertex_list= qh vertex_tail;
2651
for (facet= qh facet_list; facet && facet->next; facet= nextfacet) { /* non-simplicial facets moved to end */
2652
nextfacet= facet->next;
2653
if (facet->visible || facet->simplicial)
2655
/* triangulate all non-simplicial facets, otherwise merging does not work, e.g., RBOX c P-0.1 P+0.1 P+0.1 D3 | QHULL d Qt Tv */
2656
if (!new_facet_list)
2657
new_facet_list= facet; /* will be moved to end */
2658
qh_triangulate_facet (facet, &new_vertex_list);
2660
trace2((qh ferr, "qh_triangulate: delete null facets from f%d -- apex same as second vertex\n", getid_(new_facet_list)));
2661
for (facet= new_facet_list; facet && facet->next; facet= nextfacet) { /* null facets moved to end */
2662
nextfacet= facet->next;
2665
if (facet->ridges) {
2666
if (qh_setsize(facet->ridges) > 0) {
2667
fprintf( qh ferr, "qhull error (qh_triangulate): ridges still defined for f%d\n", facet->id);
2668
qh_errexit (qh_ERRqhull, facet, NULL);
2670
qh_setfree (&facet->ridges);
2672
if (SETfirst_(facet->vertices) == SETsecond_(facet->vertices)) {
2674
qh_triangulate_null (facet);
2677
trace2((qh ferr, "qh_triangulate: delete %d or more mirror facets -- same vertices and neighbors\n", qh_setsize(qh degen_mergeset)));
2678
qh visible_list= qh facet_tail;
2679
while ((merge= (mergeT*)qh_setdellast (qh degen_mergeset))) {
2680
facet1= merge->facet1;
2681
facet2= merge->facet2;
2682
mergetype= merge->type;
2683
qh_memfree (merge, sizeof(mergeT));
2684
if (mergetype == MRGmirror) {
2686
qh_triangulate_mirror (facet1, facet2);
2689
qh_settempfree(&qh degen_mergeset);
2690
trace2((qh ferr, "qh_triangulate: update neighbor lists for vertices from v%d\n", getid_(new_vertex_list)));
2691
qh newvertex_list= new_vertex_list; /* all vertices of new facets */
2692
qh visible_list= NULL;
2693
qh_updatevertices(/*qh newvertex_list, empty newfacet_list and visible_list*/);
2694
qh_resetlists (False, !qh_RESETvisible /*qh newvertex_list, empty newfacet_list and visible_list*/);
2696
trace2((qh ferr, "qh_triangulate: identify degenerate tricoplanar facets from f%d\n", getid_(new_facet_list)));
2697
trace2((qh ferr, "qh_triangulate: and replace facet->f.triowner with tricoplanar facets that own center, normal, etc.\n"));
2698
FORALLfacet_(new_facet_list) {
2699
if (facet->tricoplanar && !facet->visible) {
2700
FOREACHneighbor_i_(facet) {
2701
if (neighbor_i == 0) { /* first iteration */
2702
if (neighbor->tricoplanar)
2703
orig_neighbor= neighbor->f.triowner;
2705
orig_neighbor= neighbor;
2707
if (neighbor->tricoplanar)
2708
otherfacet= neighbor->f.triowner;
2710
otherfacet= neighbor;
2711
if (orig_neighbor == otherfacet) {
2713
facet->degenerate= True;
2721
trace2((qh ferr, "qh_triangulate: delete visible facets -- non-simplicial, null, and mirrored facets\n"));
2724
for (facet= new_facet_list; facet && facet->next; facet= nextfacet) { /* may delete facet */
2725
nextfacet= facet->next;
2726
if (facet->visible) {
2727
if (facet->tricoplanar) { /* a null or mirrored facet */
2730
}else { /* a non-simplicial facet followed by its tricoplanars */
2731
if (visible && !owner) {
2732
/* RBOX 200 s D5 t1001471447 | QHULL Qt C-0.01 Qx Qc Tv Qt -- f4483 had 6 vertices/neighbors and 8 ridges */
2733
trace2((qh ferr, "qh_triangulate: all tricoplanar facets degenerate for non-simplicial facet f%d\n",
2735
qh_delfacet(visible);
2741
}else if (facet->tricoplanar) {
2742
if (facet->f.triowner != visible) {
2743
fprintf( qh ferr, "qhull error (qh_triangulate): tricoplanar facet f%d not owned by its visible, non-simplicial facet f%d\n", facet->id, getid_(visible));
2744
qh_errexit2 (qh_ERRqhull, facet, visible);
2747
facet->f.triowner= owner;
2748
else if (!facet->degenerate) {
2750
nextfacet= visible->next; /* rescan tricoplanar facets with owner */
2751
facet->keepcentrum= True; /* one facet owns ->normal, etc. */
2752
facet->coplanarset= visible->coplanarset;
2753
facet->outsideset= visible->outsideset;
2754
visible->coplanarset= NULL;
2755
visible->outsideset= NULL;
2756
if (!qh TRInormals) { /* center and normal copied to tricoplanar facets */
2757
visible->center= NULL;
2758
visible->normal= NULL;
2760
qh_delfacet(visible);
2765
if (visible && !owner) {
2766
trace2((qh ferr, "qh_triangulate: all tricoplanar facets degenerate for last non-simplicial facet f%d\n",
2768
qh_delfacet(visible);
2771
qh NEWfacets= False;
2772
qh ONLYgood= onlygood; /* restore value */
2773
if (qh CHECKfrequently)
2774
qh_checkpolygon (qh facet_list);
2778
/*-<a href="qh-poly.htm#TOC"
2779
>-------------------------------</a><a name="triangulate_facet">-</a>
2781
qh_triangulate_facet (facetA)
2782
triangulate a non-simplicial facet
2783
if qh.CENTERtype=qh_ASvoronoi, sets its Voronoi center
2785
qh.newfacet_list == simplicial facets
2786
facet->tricoplanar set and ->keepcentrum false
2787
facet->degenerate set if duplicated apex
2788
facet->f.trivisible set to facetA
2789
facet->center copied from facetA (created if qh_ASvoronoi)
2790
qh_eachvoronoi, qh_detvridge, qh_detvridge3 assume centers copied
2791
facet->normal,offset,maxoutside copied from facetA
2794
qh_makenew_nonsimplicial uses neighbor->seen for the same
2797
qh_addpoint() -- add a point
2798
qh_makenewfacets() -- construct a cone of facets for a new vertex
2802
compute Voronoi center (facet->center)
2803
select first vertex (highest ID to preserve ID ordering of ->vertices)
2804
triangulate from vertex to ridges
2805
copy facet->center, normal, offset
2806
update vertex neighbors
2808
void qh_triangulate_facet (facetT *facetA, vertexT **first_vertex) {
2810
facetT *neighbor, **neighborp;
2814
trace3((qh ferr, "qh_triangulate_facet: triangulate facet f%d\n", facetA->id));
2816
if (qh IStracing >= 4)
2817
qh_printfacet (qh ferr, facetA);
2818
FOREACHneighbor_(facetA) {
2819
neighbor->seen= False;
2820
neighbor->coplanar= False;
2822
if (qh CENTERtype == qh_ASvoronoi && !facetA->center /* matches upperdelaunay in qh_setfacetplane() */
2823
&& fabs_(facetA->normal[qh hull_dim -1]) >= qh ANGLEround * qh_ZEROdelaunay) {
2824
facetA->center= qh_facetcenter (facetA->vertices);
2826
qh_willdelete (facetA, NULL);
2827
qh newfacet_list= qh facet_tail;
2828
facetA->visitid= qh visit_id;
2829
apex= SETfirst_(facetA->vertices);
2830
qh_makenew_nonsimplicial (facetA, apex, &numnew);
2831
SETfirst_(facetA->neighbors)= NULL;
2833
newfacet->tricoplanar= True;
2834
newfacet->f.trivisible= facetA;
2835
newfacet->degenerate= False;
2836
newfacet->upperdelaunay= facetA->upperdelaunay;
2837
newfacet->good= facetA->good;
2838
if (qh TRInormals) {
2839
newfacet->keepcentrum= True;
2840
newfacet->normal= qh_copypoints (facetA->normal, 1, qh hull_dim);
2841
if (qh CENTERtype == qh_AScentrum)
2842
newfacet->center= qh_getcentrum (newfacet);
2844
newfacet->center= qh_copypoints (facetA->center, 1, qh hull_dim);
2846
newfacet->keepcentrum= False;
2847
newfacet->normal= facetA->normal;
2848
newfacet->center= facetA->center;
2850
newfacet->offset= facetA->offset;
2852
newfacet->maxoutside= facetA->maxoutside;
2855
qh_matchnewfacets(/*qh newfacet_list*/);
2856
zinc_(Ztricoplanar);
2857
zadd_(Ztricoplanartot, numnew);
2858
zmax_(Ztricoplanarmax, numnew);
2859
qh visible_list= NULL;
2860
if (!(*first_vertex))
2861
(*first_vertex)= qh newvertex_list;
2862
qh newvertex_list= NULL;
2863
qh_updatevertices(/*qh newfacet_list, empty visible_list and newvertex_list*/);
2864
qh_resetlists (False, !qh_RESETvisible /*qh newfacet_list, empty visible_list and newvertex_list*/);
2865
} /* triangulate_facet */
2867
/*-<a href="qh-poly.htm#TOC"
2868
>-------------------------------</a><a name="triangulate_link">-</a>
2870
qh_triangulate_link (oldfacetA, facetA, oldfacetB, facetB)
2871
relink facetA to facetB via oldfacets
2873
adds mirror facets to qh degen_mergeset (4-d and up only)
2875
if they are already neighbors, the opposing neighbors become MRGmirror facets
2877
void qh_triangulate_link (facetT *oldfacetA, facetT *facetA, facetT *oldfacetB, facetT *facetB) {
2878
int errmirror= False;
2880
trace3((qh ferr, "qh_triangulate_link: relink old facets f%d and f%d between neighbors f%d and f%d\n",
2881
oldfacetA->id, oldfacetB->id, facetA->id, facetB->id));
2882
if (qh_setin (facetA->neighbors, facetB)) {
2883
if (!qh_setin (facetB->neighbors, facetA))
2886
qh_appendmergeset (facetA, facetB, MRGmirror, NULL);
2887
}else if (qh_setin (facetB->neighbors, facetA))
2890
fprintf( qh ferr, "qhull error (qh_triangulate_link): mirror facets f%d and f%d do not match for old facets f%d and f%d\n",
2891
facetA->id, facetB->id, oldfacetA->id, oldfacetB->id);
2892
qh_errexit2 (qh_ERRqhull, facetA, facetB);
2894
qh_setreplace (facetB->neighbors, oldfacetB, facetA);
2895
qh_setreplace (facetA->neighbors, oldfacetA, facetB);
2896
} /* triangulate_link */
2898
/*-<a href="qh-poly.htm#TOC"
2899
>-------------------------------</a><a name="triangulate_mirror">-</a>
2901
qh_triangulate_mirror (facetA, facetB)
2902
delete mirrored facets from qh_triangulate_null() and qh_triangulate_mirror
2903
a mirrored facet shares the same vertices of a logical ridge
2905
since a null facet duplicates the first two vertices, the opposing neighbors absorb the null facet
2906
if they are already neighbors, the opposing neighbors become MRGmirror facets
2908
void qh_triangulate_mirror (facetT *facetA, facetT *facetB) {
2909
facetT *neighbor, *neighborB;
2910
int neighbor_i, neighbor_n;
2912
trace3((qh ferr, "qh_triangulate_mirror: delete mirrored facets f%d and f%d\n",
2913
facetA->id, facetB->id));
2914
FOREACHneighbor_i_(facetA) {
2915
neighborB= SETelemt_(facetB->neighbors, neighbor_i, facetT);
2916
if (neighbor == neighborB)
2917
continue; /* occurs twice */
2918
qh_triangulate_link (facetA, neighbor, facetB, neighborB);
2920
qh_willdelete (facetA, NULL);
2921
qh_willdelete (facetB, NULL);
2922
} /* triangulate_mirror */
2924
/*-<a href="qh-poly.htm#TOC"
2925
>-------------------------------</a><a name="triangulate_null">-</a>
2927
qh_triangulate_null (facetA)
2928
remove null facetA from qh_triangulate_facet()
2929
a null facet has vertex #1 (apex) == vertex #2
2931
adds facetA to ->visible for deletion after qh_updatevertices
2932
qh degen_mergeset contains mirror facets (4-d and up only)
2934
since a null facet duplicates the first two vertices, the opposing neighbors absorb the null facet
2935
if they are already neighbors, the opposing neighbors become MRGmirror facets
2937
void qh_triangulate_null (facetT *facetA) {
2938
facetT *neighbor, *otherfacet;
2940
trace3((qh ferr, "qh_triangulate_null: delete null facet f%d\n", facetA->id));
2941
neighbor= SETfirst_(facetA->neighbors);
2942
otherfacet= SETsecond_(facetA->neighbors);
2943
qh_triangulate_link (facetA, neighbor, facetA, otherfacet);
2944
qh_willdelete (facetA, NULL);
2945
} /* triangulate_null */
2947
#else /* qh_NOmerge */
2948
void qh_triangulate (void) {
2950
#endif /* qh_NOmerge */
2952
/*-<a href="qh-poly.htm#TOC"
2953
>-------------------------------</a><a name="vertexintersect">-</a>
2955
qh_vertexintersect( vertexsetA, vertexsetB )
2956
intersects two vertex sets (inverse id ordered)
2957
vertexsetA is a temporary set at the top of qhmem.tempstack
2960
replaces vertexsetA with the intersection
2963
could overwrite vertexsetA if currently too slow
2965
void qh_vertexintersect(setT **vertexsetA,setT *vertexsetB) {
2968
intersection= qh_vertexintersect_new (*vertexsetA, vertexsetB);
2969
qh_settempfree (vertexsetA);
2970
*vertexsetA= intersection;
2971
qh_settemppush (intersection);
2972
} /* vertexintersect */
2974
/*-<a href="qh-poly.htm#TOC"
2975
>-------------------------------</a><a name="vertexintersect_new">-</a>
2977
qh_vertexintersect_new( )
2978
intersects two vertex sets (inverse id ordered)
2983
setT *qh_vertexintersect_new (setT *vertexsetA,setT *vertexsetB) {
2984
setT *intersection= qh_setnew (qh hull_dim - 1);
2985
vertexT **vertexA= SETaddr_(vertexsetA, vertexT);
2986
vertexT **vertexB= SETaddr_(vertexsetB, vertexT);
2988
while (*vertexA && *vertexB) {
2989
if (*vertexA == *vertexB) {
2990
qh_setappend(&intersection, *vertexA);
2991
vertexA++; vertexB++;
2993
if ((*vertexA)->id > (*vertexB)->id)
2999
return intersection;
3000
} /* vertexintersect_new */
3002
/*-<a href="qh-poly.htm#TOC"
3003
>-------------------------------</a><a name="vertexneighbors">-</a>
3005
qh_vertexneighbors()
3006
for each vertex in qh.facet_list,
3007
determine its neighboring facets
3010
sets qh.VERTEXneighbors
3011
nop if qh.VERTEXneighbors already set
3012
qh_addpoint() will maintain them
3015
assumes all vertex->neighbors are NULL
3020
append facet to vertex->neighbors
3022
void qh_vertexneighbors (void /*qh facet_list*/) {
3024
vertexT *vertex, **vertexp;
3026
if (qh VERTEXneighbors)
3028
trace1((qh ferr, "qh_vertexneighbors: determing neighboring facets for each vertex\n"));
3033
FOREACHvertex_(facet->vertices) {
3034
if (vertex->visitid != qh vertex_visit) {
3035
vertex->visitid= qh vertex_visit;
3036
vertex->neighbors= qh_setnew (qh hull_dim);
3038
qh_setappend (&vertex->neighbors, facet);
3041
qh VERTEXneighbors= True;
3042
} /* vertexneighbors */
3044
/*-<a href="qh-poly.htm#TOC"
3045
>-------------------------------</a><a name="vertexsubset">-</a>
3047
qh_vertexsubset( vertexsetA, vertexsetB )
3048
returns True if vertexsetA is a subset of vertexsetB
3049
assumes vertexsets are sorted
3052
empty set is a subset of any other set
3054
boolT qh_vertexsubset(setT *vertexsetA, setT *vertexsetB) {
3055
vertexT **vertexA= (vertexT **) SETaddr_(vertexsetA, vertexT);
3056
vertexT **vertexB= (vertexT **) SETaddr_(vertexsetB, vertexT);
3063
if ((*vertexA)->id > (*vertexB)->id)
3065
if (*vertexA == *vertexB)
3069
return False; /* avoid warnings */
3070
} /* vertexsubset */
added
removed
extern/qhull/src/poly2.c
