00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include "qhull_a.h"
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 void qh_addhash(void* newelem, setT *hashtable, int hashsize, int hash) {
00027 int scan;
00028 void *elem;
00029
00030 for (scan= (int)hash; (elem= SETelem_(hashtable, scan));
00031 scan= (++scan >= hashsize ? 0 : scan)) {
00032 if (elem == newelem)
00033 break;
00034 }
00035
00036 if (!elem)
00037 SETelem_(hashtable, scan)= newelem;
00038 }
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066 void qh_check_bestdist(void) {
00067 boolT waserror= False, unassigned;
00068 facetT *facet, *bestfacet, *errfacet1= NULL, *errfacet2= NULL;
00069 facetT *facetlist;
00070 realT dist, maxoutside, maxdist= -REALmax;
00071 pointT *point;
00072 int numpart= 0, facet_i, facet_n, notgood= 0, notverified= 0;
00073 setT *facets;
00074
00075 trace1((qh ferr, 1020, "qh_check_bestdist: check points below nearest facet. Facet_list f%d\n",
00076 qh facet_list->id));
00077 maxoutside= qh_maxouter();
00078 maxoutside += qh DISTround;
00079
00080 trace1((qh ferr, 1021, "qh_check_bestdist: check that all points are within %2.2g of best facet\n", maxoutside));
00081 facets= qh_pointfacet();
00082 if (!qh_QUICKhelp && qh PRINTprecision)
00083 qh_fprintf(qh ferr, 8091, "\n\
00084 qhull output completed. Verifying that %d points are\n\
00085 below %2.2g of the nearest %sfacet.\n",
00086 qh_setsize(facets), maxoutside, (qh ONLYgood ? "good " : ""));
00087 FOREACHfacet_i_(facets) {
00088 if (facet)
00089 unassigned= False;
00090 else {
00091 unassigned= True;
00092 facet= qh facet_list;
00093 }
00094 point= qh_point(facet_i);
00095 if (point == qh GOODpointp)
00096 continue;
00097 qh_distplane(point, facet, &dist);
00098 numpart++;
00099 bestfacet= qh_findbesthorizon(!qh_IScheckmax, point, facet, qh_NOupper, &dist, &numpart);
00100
00101 maximize_(maxdist, dist);
00102 if (dist > maxoutside) {
00103 if (qh ONLYgood && !bestfacet->good
00104 && !((bestfacet= qh_findgooddist(point, bestfacet, &dist, &facetlist))
00105 && dist > maxoutside))
00106 notgood++;
00107 else {
00108 waserror= True;
00109 qh_fprintf(qh ferr, 6109, "qhull precision error: point p%d is outside facet f%d, distance= %6.8g maxoutside= %6.8g\n",
00110 facet_i, bestfacet->id, dist, maxoutside);
00111 if (errfacet1 != bestfacet) {
00112 errfacet2= errfacet1;
00113 errfacet1= bestfacet;
00114 }
00115 }
00116 }else if (unassigned && dist < -qh MAXcoplanar)
00117 notverified++;
00118 }
00119 qh_settempfree(&facets);
00120 if (notverified && !qh DELAUNAY && !qh_QUICKhelp && qh PRINTprecision)
00121 qh_fprintf(qh ferr, 8092, "\n%d points were well inside the hull. If the hull contains\n\
00122 a lens-shaped component, these points were not verified. Use\n\
00123 options 'Qci Tv' to verify all points.\n", notverified);
00124 if (maxdist > qh outside_err) {
00125 qh_fprintf(qh ferr, 6110, "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",
00126 maxdist, qh outside_err);
00127 qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);
00128 }else if (waserror && qh outside_err > REALmax/2)
00129 qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);
00130
00131 trace0((qh ferr, 20, "qh_check_bestdist: max distance outside %2.2g\n", maxdist));
00132 }
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166 #ifndef qh_NOmerge
00167 void qh_check_maxout(void) {
00168 facetT *facet, *bestfacet, *neighbor, **neighborp, *facetlist;
00169 realT dist, maxoutside, minvertex, old_maxoutside;
00170 pointT *point;
00171 int numpart= 0, facet_i, facet_n, notgood= 0;
00172 setT *facets, *vertices;
00173 vertexT *vertex;
00174
00175 trace1((qh ferr, 1022, "qh_check_maxout: check and update maxoutside for each facet.\n"));
00176 maxoutside= minvertex= 0;
00177 if (qh VERTEXneighbors
00178 && (qh PRINTsummary || qh KEEPinside || qh KEEPcoplanar
00179 || qh TRACElevel || qh PRINTstatistics
00180 || qh PRINTout[0] == qh_PRINTsummary || qh PRINTout[0] == qh_PRINTnone)) {
00181 trace1((qh ferr, 1023, "qh_check_maxout: determine actual maxoutside and minvertex\n"));
00182 vertices= qh_pointvertex();
00183 FORALLvertices {
00184 FOREACHneighbor_(vertex) {
00185 zinc_(Zdistvertex);
00186 qh_distplane(vertex->point, neighbor, &dist);
00187 minimize_(minvertex, dist);
00188 if (-dist > qh TRACEdist || dist > qh TRACEdist
00189 || neighbor == qh tracefacet || vertex == qh tracevertex)
00190 qh_fprintf(qh ferr, 8093, "qh_check_maxout: p%d(v%d) is %.2g from f%d\n",
00191 qh_pointid(vertex->point), vertex->id, dist, neighbor->id);
00192 }
00193 }
00194 if (qh MERGING) {
00195 wmin_(Wminvertex, qh min_vertex);
00196 }
00197 qh min_vertex= minvertex;
00198 qh_settempfree(&vertices);
00199 }
00200 facets= qh_pointfacet();
00201 do {
00202 old_maxoutside= fmax_(qh max_outside, maxoutside);
00203 FOREACHfacet_i_(facets) {
00204 if (facet) {
00205 point= qh_point(facet_i);
00206 if (point == qh GOODpointp)
00207 continue;
00208 zzinc_(Ztotcheck);
00209 qh_distplane(point, facet, &dist);
00210 numpart++;
00211 bestfacet= qh_findbesthorizon(qh_IScheckmax, point, facet, !qh_NOupper, &dist, &numpart);
00212 if (bestfacet && dist > maxoutside) {
00213 if (qh ONLYgood && !bestfacet->good
00214 && !((bestfacet= qh_findgooddist(point, bestfacet, &dist, &facetlist))
00215 && dist > maxoutside))
00216 notgood++;
00217 else
00218 maxoutside= dist;
00219 }
00220 if (dist > qh TRACEdist || (bestfacet && bestfacet == qh tracefacet))
00221 qh_fprintf(qh ferr, 8094, "qh_check_maxout: p%d is %.2g above f%d\n",
00222 qh_pointid(point), dist, bestfacet->id);
00223 }
00224 }
00225 }while
00226 (maxoutside > 2*old_maxoutside);
00227
00228
00229 zzadd_(Zcheckpart, numpart);
00230 qh_settempfree(&facets);
00231 wval_(Wmaxout)= maxoutside - qh max_outside;
00232 wmax_(Wmaxoutside, qh max_outside);
00233 qh max_outside= maxoutside;
00234 qh_nearcoplanar();
00235 qh maxoutdone= True;
00236 trace1((qh ferr, 1024, "qh_check_maxout: maxoutside %2.2g, min_vertex %2.2g, outside of not good %d\n",
00237 maxoutside, qh min_vertex, notgood));
00238 }
00239 #else
00240 void qh_check_maxout(void) {
00241 }
00242 #endif
00243
00244
00245
00246
00247
00248
00249
00250
00251 void qh_check_output(void) {
00252 int i;
00253
00254 if (qh STOPcone)
00255 return;
00256 if (qh VERIFYoutput | qh IStracing | qh CHECKfrequently) {
00257 qh_checkpolygon(qh facet_list);
00258 qh_checkflipped_all(qh facet_list);
00259 qh_checkconvex(qh facet_list, qh_ALGORITHMfault);
00260 }else if (!qh MERGING && qh_newstats(qhstat precision, &i)) {
00261 qh_checkflipped_all(qh facet_list);
00262 qh_checkconvex(qh facet_list, qh_ALGORITHMfault);
00263 }
00264 }
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274 void qh_check_point(pointT *point, facetT *facet, realT *maxoutside, realT *maxdist, facetT **errfacet1, facetT **errfacet2) {
00275 realT dist;
00276
00277
00278 qh_distplane(point, facet, &dist);
00279 if (dist > *maxoutside) {
00280 if (*errfacet1 != facet) {
00281 *errfacet2= *errfacet1;
00282 *errfacet1= facet;
00283 }
00284 qh_fprintf(qh ferr, 6111, "qhull precision error: point p%d is outside facet f%d, distance= %6.8g maxoutside= %6.8g\n",
00285 qh_pointid(point), facet->id, dist, *maxoutside);
00286 }
00287 maximize_(*maxdist, dist);
00288 }
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314 void qh_check_points(void) {
00315 facetT *facet, *errfacet1= NULL, *errfacet2= NULL;
00316 realT total, maxoutside, maxdist= -REALmax;
00317 pointT *point, **pointp, *pointtemp;
00318 boolT testouter;
00319
00320 maxoutside= qh_maxouter();
00321 maxoutside += qh DISTround;
00322
00323 trace1((qh ferr, 1025, "qh_check_points: check all points below %2.2g of all facet planes\n",
00324 maxoutside));
00325 if (qh num_good)
00326 total= (float)qh num_good * (float)qh num_points;
00327 else
00328 total= (float)qh num_facets * (float)qh num_points;
00329 if (total >= qh_VERIFYdirect && !qh maxoutdone) {
00330 if (!qh_QUICKhelp && qh SKIPcheckmax && qh MERGING)
00331 qh_fprintf(qh ferr, 7075, "qhull input warning: merging without checking outer planes('Q5' or 'Po').\n\
00332 Verify may report that a point is outside of a facet.\n");
00333 qh_check_bestdist();
00334 }else {
00335 if (qh_MAXoutside && qh maxoutdone)
00336 testouter= True;
00337 else
00338 testouter= False;
00339 if (!qh_QUICKhelp) {
00340 if (qh MERGEexact)
00341 qh_fprintf(qh ferr, 7076, "qhull input warning: exact merge ('Qx'). Verify may report that a point\n\
00342 is outside of a facet. See qh-optq.htm#Qx\n");
00343 else if (qh SKIPcheckmax || qh NOnearinside)
00344 qh_fprintf(qh ferr, 7077, "qhull input warning: no outer plane check ('Q5') or no processing of\n\
00345 near-inside points ('Q8'). Verify may report that a point is outside\n\
00346 of a facet.\n");
00347 }
00348 if (qh PRINTprecision) {
00349 if (testouter)
00350 qh_fprintf(qh ferr, 8098, "\n\
00351 Output completed. Verifying that all points are below outer planes of\n\
00352 all %sfacets. Will make %2.0f distance computations.\n",
00353 (qh ONLYgood ? "good " : ""), total);
00354 else
00355 qh_fprintf(qh ferr, 8099, "\n\
00356 Output completed. Verifying that all points are below %2.2g of\n\
00357 all %sfacets. Will make %2.0f distance computations.\n",
00358 maxoutside, (qh ONLYgood ? "good " : ""), total);
00359 }
00360 FORALLfacets {
00361 if (!facet->good && qh ONLYgood)
00362 continue;
00363 if (facet->flipped)
00364 continue;
00365 if (!facet->normal) {
00366 qh_fprintf(qh ferr, 7061, "qhull warning (qh_check_points): missing normal for facet f%d\n", facet->id);
00367 continue;
00368 }
00369 if (testouter) {
00370 #if qh_MAXoutside
00371 maxoutside= facet->maxoutside + 2* qh DISTround;
00372
00373 #endif
00374 }
00375 FORALLpoints {
00376 if (point != qh GOODpointp)
00377 qh_check_point(point, facet, &maxoutside, &maxdist, &errfacet1, &errfacet2);
00378 }
00379 FOREACHpoint_(qh other_points) {
00380 if (point != qh GOODpointp)
00381 qh_check_point(point, facet, &maxoutside, &maxdist, &errfacet1, &errfacet2);
00382 }
00383 }
00384 if (maxdist > qh outside_err) {
00385 qh_fprintf(qh ferr, 6112, "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",
00386 maxdist, qh outside_err );
00387 qh_errexit2( qh_ERRprec, errfacet1, errfacet2 );
00388 }else if (errfacet1 && qh outside_err > REALmax/2)
00389 qh_errexit2( qh_ERRprec, errfacet1, errfacet2 );
00390
00391 trace0((qh ferr, 21, "qh_check_points: max distance outside %2.2g\n", maxdist));
00392 }
00393 }
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424 void qh_checkconvex(facetT *facetlist, int fault) {
00425 facetT *facet, *neighbor, **neighborp, *errfacet1=NULL, *errfacet2=NULL;
00426 vertexT *vertex;
00427 realT dist;
00428 pointT *centrum;
00429 boolT waserror= False, centrum_warning= False, tempcentrum= False, allsimplicial;
00430 int neighbor_i;
00431
00432 trace1((qh ferr, 1026, "qh_checkconvex: check all ridges are convex\n"));
00433 if (!qh RERUN) {
00434 zzval_(Zconcaveridges)= 0;
00435 zzval_(Zcoplanarridges)= 0;
00436 }
00437 FORALLfacet_(facetlist) {
00438 if (facet->flipped) {
00439 qh_precision("flipped facet");
00440 qh_fprintf(qh ferr, 6113, "qhull precision error: f%d is flipped(interior point is outside)\n",
00441 facet->id);
00442 errfacet1= facet;
00443 waserror= True;
00444 continue;
00445 }
00446 if (qh MERGING && (!qh ZEROcentrum || !facet->simplicial || facet->tricoplanar))
00447 allsimplicial= False;
00448 else {
00449 allsimplicial= True;
00450 neighbor_i= 0;
00451 FOREACHneighbor_(facet) {
00452 vertex= SETelemt_(facet->vertices, neighbor_i++, vertexT);
00453 if (!neighbor->simplicial || neighbor->tricoplanar) {
00454 allsimplicial= False;
00455 continue;
00456 }
00457 qh_distplane(vertex->point, neighbor, &dist);
00458 if (dist > -qh DISTround) {
00459 if (fault == qh_DATAfault) {
00460 qh_precision("coplanar or concave ridge");
00461 qh_fprintf(qh ferr, 6114, "qhull precision error: initial simplex is not convex. Distance=%.2g\n", dist);
00462 qh_errexit(qh_ERRsingular, NULL, NULL);
00463 }
00464 if (dist > qh DISTround) {
00465 zzinc_(Zconcaveridges);
00466 qh_precision("concave ridge");
00467 qh_fprintf(qh ferr, 6115, "qhull precision error: f%d is concave to f%d, since p%d(v%d) is %6.4g above\n",
00468 facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist);
00469 errfacet1= facet;
00470 errfacet2= neighbor;
00471 waserror= True;
00472 }else if (qh ZEROcentrum) {
00473 if (dist > 0) {
00474 zzinc_(Zcoplanarridges);
00475 qh_precision("coplanar ridge");
00476 qh_fprintf(qh ferr, 6116, "qhull precision error: f%d is clearly not convex to f%d, since p%d(v%d) is %6.4g above\n",
00477 facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist);
00478 errfacet1= facet;
00479 errfacet2= neighbor;
00480 waserror= True;
00481 }
00482 }else {
00483 zzinc_(Zcoplanarridges);
00484 qh_precision("coplanar ridge");
00485 trace0((qh ferr, 22, "qhull precision error: f%d may be coplanar to f%d, since p%d(v%d) is within %6.4g during p%d\n",
00486 facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist, qh furthest_id));
00487 }
00488 }
00489 }
00490 }
00491 if (!allsimplicial) {
00492 if (qh CENTERtype == qh_AScentrum) {
00493 if (!facet->center)
00494 facet->center= qh_getcentrum(facet);
00495 centrum= facet->center;
00496 }else {
00497 if (!centrum_warning && (!facet->simplicial || facet->tricoplanar)) {
00498 centrum_warning= True;
00499 qh_fprintf(qh ferr, 7062, "qhull warning: recomputing centrums for convexity test. This may lead to false, precision errors.\n");
00500 }
00501 centrum= qh_getcentrum(facet);
00502 tempcentrum= True;
00503 }
00504 FOREACHneighbor_(facet) {
00505 if (qh ZEROcentrum && facet->simplicial && neighbor->simplicial)
00506 continue;
00507 if (facet->tricoplanar || neighbor->tricoplanar)
00508 continue;
00509 zzinc_(Zdistconvex);
00510 qh_distplane(centrum, neighbor, &dist);
00511 if (dist > qh DISTround) {
00512 zzinc_(Zconcaveridges);
00513 qh_precision("concave ridge");
00514 qh_fprintf(qh ferr, 6117, "qhull precision error: f%d is concave to f%d. Centrum of f%d is %6.4g above f%d\n",
00515 facet->id, neighbor->id, facet->id, dist, neighbor->id);
00516 errfacet1= facet;
00517 errfacet2= neighbor;
00518 waserror= True;
00519 }else if (dist >= 0.0) {
00520
00521 zzinc_(Zcoplanarridges);
00522 qh_precision("coplanar ridge");
00523 qh_fprintf(qh ferr, 6118, "qhull precision error: f%d is coplanar or concave to f%d. Centrum of f%d is %6.4g above f%d\n",
00524 facet->id, neighbor->id, facet->id, dist, neighbor->id);
00525 errfacet1= facet;
00526 errfacet2= neighbor;
00527 waserror= True;
00528 }
00529 }
00530 if (tempcentrum)
00531 qh_memfree(centrum, qh normal_size);
00532 }
00533 }
00534 if (waserror && !qh FORCEoutput)
00535 qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);
00536 }
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576 void qh_checkfacet(facetT *facet, boolT newmerge, boolT *waserrorp) {
00577 facetT *neighbor, **neighborp, *errother=NULL;
00578 ridgeT *ridge, **ridgep, *errridge= NULL, *ridge2;
00579 vertexT *vertex, **vertexp;
00580 unsigned previousid= INT_MAX;
00581 int numneighbors, numvertices, numridges=0, numRvertices=0;
00582 boolT waserror= False;
00583 int skipA, skipB, ridge_i, ridge_n, i;
00584 setT *intersection;
00585
00586 if (facet->visible) {
00587 qh_fprintf(qh ferr, 6119, "qhull internal error (qh_checkfacet): facet f%d is on the visible_list\n",
00588 facet->id);
00589 qh_errexit(qh_ERRqhull, facet, NULL);
00590 }
00591 if (!facet->normal) {
00592 qh_fprintf(qh ferr, 6120, "qhull internal error (qh_checkfacet): facet f%d does not have a normal\n",
00593 facet->id);
00594 waserror= True;
00595 }
00596 qh_setcheck(facet->vertices, "vertices for f", facet->id);
00597 qh_setcheck(facet->ridges, "ridges for f", facet->id);
00598 qh_setcheck(facet->outsideset, "outsideset for f", facet->id);
00599 qh_setcheck(facet->coplanarset, "coplanarset for f", facet->id);
00600 qh_setcheck(facet->neighbors, "neighbors for f", facet->id);
00601 FOREACHvertex_(facet->vertices) {
00602 if (vertex->deleted) {
00603 qh_fprintf(qh ferr, 6121, "qhull internal error (qh_checkfacet): deleted vertex v%d in f%d\n", vertex->id, facet->id);
00604 qh_errprint("ERRONEOUS", NULL, NULL, NULL, vertex);
00605 waserror= True;
00606 }
00607 if (vertex->id >= previousid) {
00608 qh_fprintf(qh ferr, 6122, "qhull internal error (qh_checkfacet): vertices of f%d are not in descending id order at v%d\n", facet->id, vertex->id);
00609 waserror= True;
00610 break;
00611 }
00612 previousid= vertex->id;
00613 }
00614 numneighbors= qh_setsize(facet->neighbors);
00615 numvertices= qh_setsize(facet->vertices);
00616 numridges= qh_setsize(facet->ridges);
00617 if (facet->simplicial) {
00618 if (numvertices+numneighbors != 2*qh hull_dim
00619 && !facet->degenerate && !facet->redundant) {
00620 qh_fprintf(qh ferr, 6123, "qhull internal error (qh_checkfacet): for simplicial facet f%d, #vertices %d + #neighbors %d != 2*qh hull_dim\n",
00621 facet->id, numvertices, numneighbors);
00622 qh_setprint(qh ferr, "", facet->neighbors);
00623 waserror= True;
00624 }
00625 }else {
00626 if (!newmerge
00627 &&(numvertices < qh hull_dim || numneighbors < qh hull_dim)
00628 && !facet->degenerate && !facet->redundant) {
00629 qh_fprintf(qh ferr, 6124, "qhull internal error (qh_checkfacet): for facet f%d, #vertices %d or #neighbors %d < qh hull_dim\n",
00630 facet->id, numvertices, numneighbors);
00631 waserror= True;
00632 }
00633
00634 if (numridges < numneighbors
00635 ||(qh hull_dim == 3 && numvertices > numridges && !qh NEWfacets)
00636 ||(qh hull_dim == 2 && numridges + numvertices + numneighbors != 6)) {
00637 if (!facet->degenerate && !facet->redundant) {
00638 qh_fprintf(qh ferr, 6125, "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",
00639 facet->id, numridges, numneighbors, numvertices);
00640 waserror= True;
00641 }
00642 }
00643 }
00644 FOREACHneighbor_(facet) {
00645 if (neighbor == qh_MERGEridge || neighbor == qh_DUPLICATEridge) {
00646 qh_fprintf(qh ferr, 6126, "qhull internal error (qh_checkfacet): facet f%d still has a MERGE or DUP neighbor\n", facet->id);
00647 qh_errexit(qh_ERRqhull, facet, NULL);
00648 }
00649 neighbor->seen= True;
00650 }
00651 FOREACHneighbor_(facet) {
00652 if (!qh_setin(neighbor->neighbors, facet)) {
00653 qh_fprintf(qh ferr, 6127, "qhull internal error (qh_checkfacet): facet f%d has neighbor f%d, but f%d does not have neighbor f%d\n",
00654 facet->id, neighbor->id, neighbor->id, facet->id);
00655 errother= neighbor;
00656 waserror= True;
00657 }
00658 if (!neighbor->seen) {
00659 qh_fprintf(qh ferr, 6128, "qhull internal error (qh_checkfacet): facet f%d has a duplicate neighbor f%d\n",
00660 facet->id, neighbor->id);
00661 errother= neighbor;
00662 waserror= True;
00663 }
00664 neighbor->seen= False;
00665 }
00666 FOREACHridge_(facet->ridges) {
00667 qh_setcheck(ridge->vertices, "vertices for r", ridge->id);
00668 ridge->seen= False;
00669 }
00670 FOREACHridge_(facet->ridges) {
00671 if (ridge->seen) {
00672 qh_fprintf(qh ferr, 6129, "qhull internal error (qh_checkfacet): facet f%d has a duplicate ridge r%d\n",
00673 facet->id, ridge->id);
00674 errridge= ridge;
00675 waserror= True;
00676 }
00677 ridge->seen= True;
00678 numRvertices= qh_setsize(ridge->vertices);
00679 if (numRvertices != qh hull_dim - 1) {
00680 qh_fprintf(qh ferr, 6130, "qhull internal error (qh_checkfacet): ridge between f%d and f%d has %d vertices\n",
00681 ridge->top->id, ridge->bottom->id, numRvertices);
00682 errridge= ridge;
00683 waserror= True;
00684 }
00685 neighbor= otherfacet_(ridge, facet);
00686 neighbor->seen= True;
00687 if (!qh_setin(facet->neighbors, neighbor)) {
00688 qh_fprintf(qh ferr, 6131, "qhull internal error (qh_checkfacet): for facet f%d, neighbor f%d of ridge r%d not in facet\n",
00689 facet->id, neighbor->id, ridge->id);
00690 errridge= ridge;
00691 waserror= True;
00692 }
00693 }
00694 if (!facet->simplicial) {
00695 FOREACHneighbor_(facet) {
00696 if (!neighbor->seen) {
00697 qh_fprintf(qh ferr, 6132, "qhull internal error (qh_checkfacet): facet f%d does not have a ridge for neighbor f%d\n",
00698 facet->id, neighbor->id);
00699 errother= neighbor;
00700 waserror= True;
00701 }
00702 intersection= qh_vertexintersect_new(facet->vertices, neighbor->vertices);
00703 qh_settemppush(intersection);
00704 FOREACHvertex_(facet->vertices) {
00705 vertex->seen= False;
00706 vertex->seen2= False;
00707 }
00708 FOREACHvertex_(intersection)
00709 vertex->seen= True;
00710 FOREACHridge_(facet->ridges) {
00711 if (neighbor != otherfacet_(ridge, facet))
00712 continue;
00713 FOREACHvertex_(ridge->vertices) {
00714 if (!vertex->seen) {
00715 qh_fprintf(qh ferr, 6133, "qhull internal error (qh_checkfacet): vertex v%d in r%d not in f%d intersect f%d\n",
00716 vertex->id, ridge->id, facet->id, neighbor->id);
00717 qh_errexit(qh_ERRqhull, facet, ridge);
00718 }
00719 vertex->seen2= True;
00720 }
00721 }
00722 if (!newmerge) {
00723 FOREACHvertex_(intersection) {
00724 if (!vertex->seen2) {
00725 if (qh IStracing >=3 || !qh MERGING) {
00726 qh_fprintf(qh ferr, 6134, "qhull precision error (qh_checkfacet): vertex v%d in f%d intersect f%d but\n\
00727 not in a ridge. This is ok under merging. Last point was p%d\n",
00728 vertex->id, facet->id, neighbor->id, qh furthest_id);
00729 if (!qh FORCEoutput && !qh MERGING) {
00730 qh_errprint("ERRONEOUS", facet, neighbor, NULL, vertex);
00731 if (!qh MERGING)
00732 qh_errexit(qh_ERRqhull, NULL, NULL);
00733 }
00734 }
00735 }
00736 }
00737 }
00738 qh_settempfree(&intersection);
00739 }
00740 }else {
00741 FOREACHneighbor_(facet) {
00742 if (neighbor->simplicial) {
00743 skipA= SETindex_(facet->neighbors, neighbor);
00744 skipB= qh_setindex(neighbor->neighbors, facet);
00745 if (!qh_setequal_skip(facet->vertices, skipA, neighbor->vertices, skipB)) {
00746 qh_fprintf(qh ferr, 6135, "qhull internal error (qh_checkfacet): facet f%d skip %d and neighbor f%d skip %d do not match \n",
00747 facet->id, skipA, neighbor->id, skipB);
00748 errother= neighbor;
00749 waserror= True;
00750 }
00751 }
00752 }
00753 }
00754 if (qh hull_dim < 5 && (qh IStracing > 2 || qh CHECKfrequently)) {
00755 FOREACHridge_i_(facet->ridges) {
00756 for (i=ridge_i+1; i < ridge_n; i++) {
00757 ridge2= SETelemt_(facet->ridges, i, ridgeT);
00758 if (qh_setequal(ridge->vertices, ridge2->vertices)) {
00759 qh_fprintf(qh ferr, 6227, "Qhull internal error (qh_checkfacet): ridges r%d and r%d have the same vertices\n",
00760 ridge->id, ridge2->id);
00761 errridge= ridge;
00762 waserror= True;
00763 }
00764 }
00765 }
00766 }
00767 if (waserror) {
00768 qh_errprint("ERRONEOUS", facet, errother, errridge, NULL);
00769 *waserrorp= True;
00770 }
00771 }
00772
00773
00774
00775
00776
00777
00778
00779
00780 void qh_checkflipped_all(facetT *facetlist) {
00781 facetT *facet;
00782 boolT waserror= False;
00783 realT dist;
00784
00785 if (facetlist == qh facet_list)
00786 zzval_(Zflippedfacets)= 0;
00787 FORALLfacet_(facetlist) {
00788 if (facet->normal && !qh_checkflipped(facet, &dist, !qh_ALL)) {
00789 qh_fprintf(qh ferr, 6136, "qhull precision error: facet f%d is flipped, distance= %6.12g\n",
00790 facet->id, dist);
00791 if (!qh FORCEoutput) {
00792 qh_errprint("ERRONEOUS", facet, NULL, NULL, NULL);
00793 waserror= True;
00794 }
00795 }
00796 }
00797 if (waserror) {
00798 qh_fprintf(qh ferr, 8101, "\n\
00799 A flipped facet occurs when its distance to the interior point is\n\
00800 greater than %2.2g, the maximum roundoff error.\n", -qh DISTround);
00801 qh_errexit(qh_ERRprec, NULL, NULL);
00802 }
00803 }
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827 void qh_checkpolygon(facetT *facetlist) {
00828 facetT *facet;
00829 vertexT *vertex, **vertexp, *vertexlist;
00830 int numfacets= 0, numvertices= 0, numridges= 0;
00831 int totvneighbors= 0, totvertices= 0;
00832 boolT waserror= False, nextseen= False, visibleseen= False;
00833
00834 trace1((qh ferr, 1027, "qh_checkpolygon: check all facets from f%d\n", facetlist->id));
00835 if (facetlist != qh facet_list || qh ONLYgood)
00836 nextseen= True;
00837 FORALLfacet_(facetlist) {
00838 if (facet == qh visible_list)
00839 visibleseen= True;
00840 if (!facet->visible) {
00841 if (!nextseen) {
00842 if (facet == qh facet_next)
00843 nextseen= True;
00844 else if (qh_setsize(facet->outsideset)) {
00845 if (!qh NARROWhull
00846 #if !qh_COMPUTEfurthest
00847 || facet->furthestdist >= qh MINoutside
00848 #endif
00849 ) {
00850 qh_fprintf(qh ferr, 6137, "qhull internal error (qh_checkpolygon): f%d has outside points before qh facet_next\n",
00851 facet->id);
00852 qh_errexit(qh_ERRqhull, facet, NULL);
00853 }
00854 }
00855 }
00856 numfacets++;
00857 qh_checkfacet(facet, False, &waserror);
00858 }
00859 }
00860 if (qh visible_list && !visibleseen && facetlist == qh facet_list) {
00861 qh_fprintf(qh ferr, 6138, "qhull internal error (qh_checkpolygon): visible list f%d no longer on facet list\n", qh visible_list->id);
00862 qh_printlists();
00863 qh_errexit(qh_ERRqhull, qh visible_list, NULL);
00864 }
00865 if (facetlist == qh facet_list)
00866 vertexlist= qh vertex_list;
00867 else if (facetlist == qh newfacet_list)
00868 vertexlist= qh newvertex_list;
00869 else
00870 vertexlist= NULL;
00871 FORALLvertex_(vertexlist) {
00872 vertex->seen= False;
00873 vertex->visitid= 0;
00874 }
00875 FORALLfacet_(facetlist) {
00876 if (facet->visible)
00877 continue;
00878 if (facet->simplicial)
00879 numridges += qh hull_dim;
00880 else
00881 numridges += qh_setsize(facet->ridges);
00882 FOREACHvertex_(facet->vertices) {
00883 vertex->visitid++;
00884 if (!vertex->seen) {
00885 vertex->seen= True;
00886 numvertices++;
00887 if (qh_pointid(vertex->point) == -1) {
00888 qh_fprintf(qh ferr, 6139, "qhull internal error (qh_checkpolygon): unknown point %p for vertex v%d first_point %p\n",
00889 vertex->point, vertex->id, qh first_point);
00890 waserror= True;
00891 }
00892 }
00893 }
00894 }
00895 qh vertex_visit += (unsigned int)numfacets;
00896 if (facetlist == qh facet_list) {
00897 if (numfacets != qh num_facets - qh num_visible) {
00898 qh_fprintf(qh ferr, 6140, "qhull internal error (qh_checkpolygon): actual number of facets is %d, cumulative facet count is %d - %d visible facets\n",
00899 numfacets, qh num_facets, qh num_visible);
00900 waserror= True;
00901 }
00902 qh vertex_visit++;
00903 if (qh VERTEXneighbors) {
00904 FORALLvertices {
00905 qh_setcheck(vertex->neighbors, "neighbors for v", vertex->id);
00906 if (vertex->deleted)
00907 continue;
00908 totvneighbors += qh_setsize(vertex->neighbors);
00909 }
00910 FORALLfacet_(facetlist)
00911 totvertices += qh_setsize(facet->vertices);
00912 if (totvneighbors != totvertices) {
00913 qh_fprintf(qh ferr, 6141, "qhull internal error (qh_checkpolygon): vertex neighbors inconsistent. Totvneighbors %d, totvertices %d\n",
00914 totvneighbors, totvertices);
00915 waserror= True;
00916 }
00917 }
00918 if (numvertices != qh num_vertices - qh_setsize(qh del_vertices)) {
00919 qh_fprintf(qh ferr, 6142, "qhull internal error (qh_checkpolygon): actual number of vertices is %d, cumulative vertex count is %d\n",
00920 numvertices, qh num_vertices - qh_setsize(qh del_vertices));
00921 waserror= True;
00922 }
00923 if (qh hull_dim == 2 && numvertices != numfacets) {
00924 qh_fprintf(qh ferr, 6143, "qhull internal error (qh_checkpolygon): #vertices %d != #facets %d\n",
00925 numvertices, numfacets);
00926 waserror= True;
00927 }
00928 if (qh hull_dim == 3 && numvertices + numfacets - numridges/2 != 2) {
00929 qh_fprintf(qh ferr, 7063, "qhull warning: #vertices %d + #facets %d - #edges %d != 2\n\
00930 A vertex appears twice in a edge list. May occur during merging.",
00931 numvertices, numfacets, numridges/2);
00932
00933 }
00934 }
00935 if (waserror)
00936 qh_errexit(qh_ERRqhull, NULL, NULL);
00937 }
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950 void qh_checkvertex(vertexT *vertex) {
00951 boolT waserror= False;
00952 facetT *neighbor, **neighborp, *errfacet=NULL;
00953
00954 if (qh_pointid(vertex->point) == -1) {
00955 qh_fprintf(qh ferr, 6144, "qhull internal error (qh_checkvertex): unknown point id %p\n", vertex->point);
00956 waserror= True;
00957 }
00958 if (vertex->id >= qh vertex_id) {
00959 qh_fprintf(qh ferr, 6145, "qhull internal error (qh_checkvertex): unknown vertex id %d\n", vertex->id);
00960 waserror= True;
00961 }
00962 if (!waserror && !vertex->deleted) {
00963 if (qh_setsize(vertex->neighbors)) {
00964 FOREACHneighbor_(vertex) {
00965 if (!qh_setin(neighbor->vertices, vertex)) {
00966 qh_fprintf(qh ferr, 6146, "qhull internal error (qh_checkvertex): neighbor f%d does not contain v%d\n", neighbor->id, vertex->id);
00967 errfacet= neighbor;
00968 waserror= True;
00969 }
00970 }
00971 }
00972 }
00973 if (waserror) {
00974 qh_errprint("ERRONEOUS", NULL, NULL, NULL, vertex);
00975 qh_errexit(qh_ERRqhull, errfacet, NULL);
00976 }
00977 }
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989 void qh_clearcenters(qh_CENTER type) {
00990 facetT *facet;
00991
00992 if (qh CENTERtype != type) {
00993 FORALLfacets {
00994 if (facet->tricoplanar && !facet->keepcentrum)
00995 facet->center= NULL;
00996 else if (qh CENTERtype == qh_ASvoronoi){
00997 if (facet->center) {
00998 qh_memfree(facet->center, qh center_size);
00999 facet->center= NULL;
01000 }
01001 }else {
01002 if (facet->center) {
01003 qh_memfree(facet->center, qh normal_size);
01004 facet->center= NULL;
01005 }
01006 }
01007 }
01008 qh CENTERtype= type;
01009 }
01010 trace2((qh ferr, 2043, "qh_clearcenters: switched to center type %d\n", type));
01011 }
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031 void qh_createsimplex(setT *vertices) {
01032 facetT *facet= NULL, *newfacet;
01033 boolT toporient= True;
01034 int vertex_i, vertex_n, nth;
01035 setT *newfacets= qh_settemp(qh hull_dim+1);
01036 vertexT *vertex;
01037
01038 qh facet_list= qh newfacet_list= qh facet_tail= qh_newfacet();
01039 qh num_facets= qh num_vertices= qh num_visible= 0;
01040 qh vertex_list= qh newvertex_list= qh vertex_tail= qh_newvertex(NULL);
01041 FOREACHvertex_i_(vertices) {
01042 newfacet= qh_newfacet();
01043 newfacet->vertices= qh_setnew_delnthsorted(vertices, vertex_n,
01044 vertex_i, 0);
01045 newfacet->toporient= (unsigned char)toporient;
01046 qh_appendfacet(newfacet);
01047 newfacet->newfacet= True;
01048 qh_appendvertex(vertex);
01049 qh_setappend(&newfacets, newfacet);
01050 toporient ^= True;
01051 }
01052 FORALLnew_facets {
01053 nth= 0;
01054 FORALLfacet_(qh newfacet_list) {
01055 if (facet != newfacet)
01056 SETelem_(newfacet->neighbors, nth++)= facet;
01057 }
01058 qh_settruncate(newfacet->neighbors, qh hull_dim);
01059 }
01060 qh_settempfree(&newfacets);
01061 trace1((qh ferr, 1028, "qh_createsimplex: created simplex\n"));
01062 }
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075 void qh_delridge(ridgeT *ridge) {
01076 void **freelistp;
01077
01078 qh_setdel(ridge->top->ridges, ridge);
01079 qh_setdel(ridge->bottom->ridges, ridge);
01080 qh_setfree(&(ridge->vertices));
01081 qh_memfree_(ridge, (int)sizeof(ridgeT), freelistp);
01082 }
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095 void qh_delvertex(vertexT *vertex) {
01096
01097 if (vertex == qh tracevertex)
01098 qh tracevertex= NULL;
01099 qh_removevertex(vertex);
01100 qh_setfree(&vertex->neighbors);
01101 qh_memfree(vertex, (int)sizeof(vertexT));
01102 }
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118 setT *qh_facet3vertex(facetT *facet) {
01119 ridgeT *ridge, *firstridge;
01120 vertexT *vertex;
01121 int cntvertices, cntprojected=0;
01122 setT *vertices;
01123
01124 cntvertices= qh_setsize(facet->vertices);
01125 vertices= qh_settemp(cntvertices);
01126 if (facet->simplicial) {
01127 if (cntvertices != 3) {
01128 qh_fprintf(qh ferr, 6147, "qhull internal error (qh_facet3vertex): only %d vertices for simplicial facet f%d\n",
01129 cntvertices, facet->id);
01130 qh_errexit(qh_ERRqhull, facet, NULL);
01131 }
01132 qh_setappend(&vertices, SETfirst_(facet->vertices));
01133 if (facet->toporient ^ qh_ORIENTclock)
01134 qh_setappend(&vertices, SETsecond_(facet->vertices));
01135 else
01136 qh_setaddnth(&vertices, 0, SETsecond_(facet->vertices));
01137 qh_setappend(&vertices, SETelem_(facet->vertices, 2));
01138 }else {
01139 ridge= firstridge= SETfirstt_(facet->ridges, ridgeT);
01140 while ((ridge= qh_nextridge3d(ridge, facet, &vertex))) {
01141 qh_setappend(&vertices, vertex);
01142 if (++cntprojected > cntvertices || ridge == firstridge)
01143 break;
01144 }
01145 if (!ridge || cntprojected != cntvertices) {
01146 qh_fprintf(qh ferr, 6148, "qhull internal error (qh_facet3vertex): ridges for facet %d don't match up. got at least %d\n",
01147 facet->id, cntprojected);
01148 qh_errexit(qh_ERRqhull, facet, ridge);
01149 }
01150 }
01151 return vertices;
01152 }
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188 facetT *qh_findbestfacet(pointT *point, boolT bestoutside,
01189 realT *bestdist, boolT *isoutside) {
01190 facetT *bestfacet= NULL;
01191 int numpart, totpart= 0;
01192
01193 bestfacet= qh_findbest(point, qh facet_list,
01194 bestoutside, !qh_ISnewfacets, bestoutside ,
01195 bestdist, isoutside, &totpart);
01196 if (*bestdist < -qh DISTround) {
01197 bestfacet= qh_findfacet_all(point, bestdist, isoutside, &numpart);
01198 totpart += numpart;
01199 if ((isoutside && bestoutside)
01200 || (!isoutside && bestfacet->upperdelaunay)) {
01201 bestfacet= qh_findbest(point, bestfacet,
01202 bestoutside, False, bestoutside,
01203 bestdist, isoutside, &totpart);
01204 totpart += numpart;
01205 }
01206 }
01207 trace3((qh ferr, 3014, "qh_findbestfacet: f%d dist %2.2g isoutside %d totpart %d\n",
01208 bestfacet->id, *bestdist, *isoutside, totpart));
01209 return bestfacet;
01210 }
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227 facetT *qh_findbestlower(facetT *upperfacet, pointT *point, realT *bestdistp, int *numpart) {
01228 facetT *neighbor, **neighborp, *bestfacet= NULL;
01229 realT bestdist= -REALmax/2 ;
01230 realT dist;
01231 vertexT *vertex;
01232
01233 zinc_(Zbestlower);
01234 FOREACHneighbor_(upperfacet) {
01235 if (neighbor->upperdelaunay || neighbor->flipped)
01236 continue;
01237 (*numpart)++;
01238 qh_distplane(point, neighbor, &dist);
01239 if (dist > bestdist) {
01240 bestfacet= neighbor;
01241 bestdist= dist;
01242 }
01243 }
01244 if (!bestfacet) {
01245 zinc_(Zbestlowerv);
01246
01247 vertex= qh_nearvertex(upperfacet, point, &dist);
01248 qh_vertexneighbors();
01249 FOREACHneighbor_(vertex) {
01250 if (neighbor->upperdelaunay || neighbor->flipped)
01251 continue;
01252 (*numpart)++;
01253 qh_distplane(point, neighbor, &dist);
01254 if (dist > bestdist) {
01255 bestfacet= neighbor;
01256 bestdist= dist;
01257 }
01258 }
01259 }
01260 if (!bestfacet) {
01261 qh_fprintf(qh ferr, 6228, "\n\
01262 Qhull internal error (qh_findbestlower): all neighbors of facet %d are flipped or upper Delaunay.\n\
01263 Please report this error to qhull_bug@qhull.org with the input and all of the output.\n",
01264 upperfacet->id);
01265 qh_errexit(qh_ERRqhull, upperfacet, NULL);
01266 }
01267 *bestdistp= bestdist;
01268 trace3((qh ferr, 3015, "qh_findbestlower: f%d dist %2.2g for f%d p%d\n",
01269 bestfacet->id, bestdist, upperfacet->id, qh_pointid(point)));
01270 return bestfacet;
01271 }
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294 facetT *qh_findfacet_all(pointT *point, realT *bestdist, boolT *isoutside,
01295 int *numpart) {
01296 facetT *bestfacet= NULL, *facet;
01297 realT dist;
01298 int totpart= 0;
01299
01300 *bestdist= -REALmax;
01301 *isoutside= False;
01302 FORALLfacets {
01303 if (facet->flipped || !facet->normal)
01304 continue;
01305 totpart++;
01306 qh_distplane(point, facet, &dist);
01307 if (dist > *bestdist) {
01308 *bestdist= dist;
01309 bestfacet= facet;
01310 if (dist > qh MINoutside) {
01311 *isoutside= True;
01312 break;
01313 }
01314 }
01315 }
01316 *numpart= totpart;
01317 trace3((qh ferr, 3016, "qh_findfacet_all: f%d dist %2.2g isoutside %d totpart %d\n",
01318 getid_(bestfacet), *bestdist, *isoutside, totpart));
01319 return bestfacet;
01320 }
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354 int qh_findgood(facetT *facetlist, int goodhorizon) {
01355 facetT *facet, *bestfacet= NULL;
01356 realT angle, bestangle= REALmax, dist;
01357 int numgood=0;
01358
01359 FORALLfacet_(facetlist) {
01360 if (facet->good)
01361 numgood++;
01362 }
01363 if (qh GOODvertex>0 && !qh MERGING) {
01364 FORALLfacet_(facetlist) {
01365 if (!qh_isvertex(qh GOODvertexp, facet->vertices)) {
01366 facet->good= False;
01367 numgood--;
01368 }
01369 }
01370 }
01371 if (qh GOODpoint && numgood) {
01372 FORALLfacet_(facetlist) {
01373 if (facet->good && facet->normal) {
01374 zinc_(Zdistgood);
01375 qh_distplane(qh GOODpointp, facet, &dist);
01376 if ((qh GOODpoint > 0) ^ (dist > 0.0)) {
01377 facet->good= False;
01378 numgood--;
01379 }
01380 }
01381 }
01382 }
01383 if (qh GOODthreshold && (numgood || goodhorizon || qh GOODclosest)) {
01384 FORALLfacet_(facetlist) {
01385 if (facet->good && facet->normal) {
01386 if (!qh_inthresholds(facet->normal, &angle)) {
01387 facet->good= False;
01388 numgood--;
01389 if (angle < bestangle) {
01390 bestangle= angle;
01391 bestfacet= facet;
01392 }
01393 }
01394 }
01395 }
01396 if (!numgood && (!goodhorizon || qh GOODclosest)) {
01397 if (qh GOODclosest) {
01398 if (qh GOODclosest->visible)
01399 qh GOODclosest= NULL;
01400 else {
01401 qh_inthresholds(qh GOODclosest->normal, &angle);
01402 if (angle < bestangle)
01403 bestfacet= qh GOODclosest;
01404 }
01405 }
01406 if (bestfacet && bestfacet != qh GOODclosest) {
01407 if (qh GOODclosest)
01408 qh GOODclosest->good= False;
01409 qh GOODclosest= bestfacet;
01410 bestfacet->good= True;
01411 numgood++;
01412 trace2((qh ferr, 2044, "qh_findgood: f%d is closest(%2.2g) to thresholds\n",
01413 bestfacet->id, bestangle));
01414 return numgood;
01415 }
01416 }else if (qh GOODclosest) {
01417 qh GOODclosest->good= False;
01418 qh GOODclosest= NULL;
01419 }
01420 }
01421 zadd_(Zgoodfacet, numgood);
01422 trace2((qh ferr, 2045, "qh_findgood: found %d good facets with %d good horizon\n",
01423 numgood, goodhorizon));
01424 if (!numgood && qh GOODvertex>0 && !qh MERGING)
01425 return goodhorizon;
01426 return numgood;
01427 }
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454 void qh_findgood_all(facetT *facetlist) {
01455 facetT *facet, *bestfacet=NULL;
01456 realT angle, bestangle= REALmax;
01457 int numgood=0, startgood;
01458
01459 if (!qh GOODvertex && !qh GOODthreshold && !qh GOODpoint
01460 && !qh SPLITthresholds)
01461 return;
01462 if (!qh ONLYgood)
01463 qh_findgood(qh facet_list, 0);
01464 FORALLfacet_(facetlist) {
01465 if (facet->good)
01466 numgood++;
01467 }
01468 if (qh GOODvertex <0 || (qh GOODvertex > 0 && qh MERGING)) {
01469 FORALLfacet_(facetlist) {
01470 if (facet->good && ((qh GOODvertex > 0) ^ !!qh_isvertex(qh GOODvertexp, facet->vertices))) {
01471 if (!--numgood) {
01472 if (qh ONLYgood) {
01473 qh_fprintf(qh ferr, 7064, "qhull warning: good vertex p%d does not match last good facet f%d. Ignored.\n",
01474 qh_pointid(qh GOODvertexp), facet->id);
01475 return;
01476 }else if (qh GOODvertex > 0)
01477 qh_fprintf(qh ferr, 7065, "qhull warning: point p%d is not a vertex('QV%d').\n",
01478 qh GOODvertex-1, qh GOODvertex-1);
01479 else
01480 qh_fprintf(qh ferr, 7066, "qhull warning: point p%d is a vertex for every facet('QV-%d').\n",
01481 -qh GOODvertex - 1, -qh GOODvertex - 1);
01482 }
01483 facet->good= False;
01484 }
01485 }
01486 }
01487 startgood= numgood;
01488 if (qh SPLITthresholds) {
01489 FORALLfacet_(facetlist) {
01490 if (facet->good) {
01491 if (!qh_inthresholds(facet->normal, &angle)) {
01492 facet->good= False;
01493 numgood--;
01494 if (angle < bestangle) {
01495 bestangle= angle;
01496 bestfacet= facet;
01497 }
01498 }
01499 }
01500 }
01501 if (!numgood && bestfacet) {
01502 bestfacet->good= True;
01503 numgood++;
01504 trace0((qh ferr, 23, "qh_findgood_all: f%d is closest(%2.2g) to thresholds\n",
01505 bestfacet->id, bestangle));
01506 return;
01507 }
01508 }
01509 qh num_good= numgood;
01510 trace0((qh ferr, 24, "qh_findgood_all: %d good facets remain out of %d facets\n",
01511 numgood, startgood));
01512 }
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524 void qh_furthestnext(void ) {
01525 facetT *facet, *bestfacet= NULL;
01526 realT dist, bestdist= -REALmax;
01527
01528 FORALLfacets {
01529 if (facet->outsideset) {
01530 #if qh_COMPUTEfurthest
01531 pointT *furthest;
01532 furthest= (pointT*)qh_setlast(facet->outsideset);
01533 zinc_(Zcomputefurthest);
01534 qh_distplane(furthest, facet, &dist);
01535 #else
01536 dist= facet->furthestdist;
01537 #endif
01538 if (dist > bestdist) {
01539 bestfacet= facet;
01540 bestdist= dist;
01541 }
01542 }
01543 }
01544 if (bestfacet) {
01545 qh_removefacet(bestfacet);
01546 qh_prependfacet(bestfacet, &qh facet_next);
01547 trace1((qh ferr, 1029, "qh_furthestnext: made f%d next facet(dist %.2g)\n",
01548 bestfacet->id, bestdist));
01549 }
01550 }
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567 void qh_furthestout(facetT *facet) {
01568 pointT *point, **pointp, *bestpoint= NULL;
01569 realT dist, bestdist= -REALmax;
01570
01571 FOREACHpoint_(facet->outsideset) {
01572 qh_distplane(point, facet, &dist);
01573 zinc_(Zcomputefurthest);
01574 if (dist > bestdist) {
01575 bestpoint= point;
01576 bestdist= dist;
01577 }
01578 }
01579 if (bestpoint) {
01580 qh_setdel(facet->outsideset, point);
01581 qh_setappend(&facet->outsideset, point);
01582 #if !qh_COMPUTEfurthest
01583 facet->furthestdist= bestdist;
01584 #endif
01585 }
01586 facet->notfurthest= False;
01587 trace3((qh ferr, 3017, "qh_furthestout: p%d is furthest outside point of f%d\n",
01588 qh_pointid(point), facet->id));
01589 }
01590
01591
01592
01593
01594
01595
01596
01597
01598 void qh_infiniteloop(facetT *facet) {
01599
01600 qh_fprintf(qh ferr, 6149, "qhull internal error (qh_infiniteloop): potential infinite loop detected\n");
01601 qh_errexit(qh_ERRqhull, facet, NULL);
01602 }
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629 void qh_initbuild( void) {
01630 setT *maxpoints, *vertices;
01631 facetT *facet;
01632 int i, numpart;
01633 realT dist;
01634 boolT isoutside;
01635
01636 qh furthest_id= -1;
01637 qh lastreport= 0;
01638 qh facet_id= qh vertex_id= qh ridge_id= 0;
01639 qh visit_id= qh vertex_visit= 0;
01640 qh maxoutdone= False;
01641
01642 if (qh GOODpoint > 0)
01643 qh GOODpointp= qh_point(qh GOODpoint-1);
01644 else if (qh GOODpoint < 0)
01645 qh GOODpointp= qh_point(-qh GOODpoint-1);
01646 if (qh GOODvertex > 0)
01647 qh GOODvertexp= qh_point(qh GOODvertex-1);
01648 else if (qh GOODvertex < 0)
01649 qh GOODvertexp= qh_point(-qh GOODvertex-1);
01650 if ((qh GOODpoint
01651 && (qh GOODpointp < qh first_point
01652 || qh GOODpointp > qh_point(qh num_points-1)))
01653 || (qh GOODvertex
01654 && (qh GOODvertexp < qh first_point
01655 || qh GOODvertexp > qh_point(qh num_points-1)))) {
01656 qh_fprintf(qh ferr, 6150, "qhull input error: either QGn or QVn point is > p%d\n",
01657 qh num_points-1);
01658 qh_errexit(qh_ERRinput, NULL, NULL);
01659 }
01660 maxpoints= qh_maxmin(qh first_point, qh num_points, qh hull_dim);
01661 if (qh SCALElast)
01662 qh_scalelast(qh first_point, qh num_points, qh hull_dim,
01663 qh MINlastcoord, qh MAXlastcoord, qh MAXwidth);
01664 qh_detroundoff();
01665 if (qh DELAUNAY && qh upper_threshold[qh hull_dim-1] > REALmax/2
01666 && qh lower_threshold[qh hull_dim-1] < -REALmax/2) {
01667 for (i=qh_PRINTEND; i--; ) {
01668 if (qh PRINTout[i] == qh_PRINTgeom && qh DROPdim < 0
01669 && !qh GOODthreshold && !qh SPLITthresholds)
01670 break;
01671 }
01672 if (i < 0) {
01673 if (qh UPPERdelaunay) {
01674 qh lower_threshold[qh hull_dim-1]= qh ANGLEround * qh_ZEROdelaunay;
01675 qh GOODthreshold= True;
01676 }else {
01677 qh upper_threshold[qh hull_dim-1]= -qh ANGLEround * qh_ZEROdelaunay;
01678 if (!qh GOODthreshold)
01679 qh SPLITthresholds= True;
01680
01681 }
01682 }
01683 }
01684 vertices= qh_initialvertices(qh hull_dim, maxpoints, qh first_point, qh num_points);
01685 qh_initialhull(vertices);
01686 qh_partitionall(vertices, qh first_point, qh num_points);
01687 if (qh PRINToptions1st || qh TRACElevel || qh IStracing) {
01688 if (qh TRACElevel || qh IStracing)
01689 qh_fprintf(qh ferr, 8103, "\nTrace level %d for %s | %s\n",
01690 qh IStracing ? qh IStracing : qh TRACElevel, qh rbox_command, qh qhull_command);
01691 qh_fprintf(qh ferr, 8104, "Options selected for Qhull %s:\n%s\n", qh_version, qh qhull_options);
01692 }
01693 qh_resetlists(False, qh_RESETvisible );
01694 qh facet_next= qh facet_list;
01695 qh_furthestnext();
01696 if (qh PREmerge) {
01697 qh cos_max= qh premerge_cos;
01698 qh centrum_radius= qh premerge_centrum;
01699 }
01700 if (qh ONLYgood) {
01701 if (qh GOODvertex > 0 && qh MERGING) {
01702 qh_fprintf(qh ferr, 6151, "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");
01703 qh_errexit(qh_ERRinput, NULL, NULL);
01704 }
01705 if (!(qh GOODthreshold || qh GOODpoint
01706 || (!qh MERGEexact && !qh PREmerge && qh GOODvertexp))) {
01707 qh_fprintf(qh ferr, 6152, "qhull input error: 'Qg' (ONLYgood) needs a good threshold('Pd0D0'), a\n\
01708 good point(QGn or QG-n), or a good vertex with 'QJ' or 'Q0' (QVn).\n");
01709 qh_errexit(qh_ERRinput, NULL, NULL);
01710 }
01711 if (qh GOODvertex > 0 && !qh MERGING
01712 && !qh_isvertex(qh GOODvertexp, vertices)) {
01713 facet= qh_findbestnew(qh GOODvertexp, qh facet_list,
01714 &dist, !qh_ALL, &isoutside, &numpart);
01715 zadd_(Zdistgood, numpart);
01716 if (!isoutside) {
01717 qh_fprintf(qh ferr, 6153, "qhull input error: point for QV%d is inside initial simplex. It can not be made a vertex.\n",
01718 qh_pointid(qh GOODvertexp));
01719 qh_errexit(qh_ERRinput, NULL, NULL);
01720 }
01721 if (!qh_addpoint(qh GOODvertexp, facet, False)) {
01722 qh_settempfree(&vertices);
01723 qh_settempfree(&maxpoints);
01724 return;
01725 }
01726 }
01727 qh_findgood(qh facet_list, 0);
01728 }
01729 qh_settempfree(&vertices);
01730 qh_settempfree(&maxpoints);
01731 trace1((qh ferr, 1030, "qh_initbuild: initial hull created and points partitioned\n"));
01732 }
01733
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748 void qh_initialhull(setT *vertices) {
01749 facetT *facet, *firstfacet, *neighbor, **neighborp;
01750 realT dist, angle, minangle= REALmax;
01751 #ifndef qh_NOtrace
01752 int k;
01753 #endif
01754
01755 qh_createsimplex(vertices);
01756 qh_resetlists(False, qh_RESETvisible);
01757 qh facet_next= qh facet_list;
01758 qh interior_point= qh_getcenter(vertices);
01759 firstfacet= qh facet_list;
01760 qh_setfacetplane(firstfacet);
01761 zinc_(Znumvisibility);
01762 qh_distplane(qh interior_point, firstfacet, &dist);
01763 if (dist > 0) {
01764 FORALLfacets
01765 facet->toporient ^= (unsigned char)True;
01766 }
01767 FORALLfacets
01768 qh_setfacetplane(facet);
01769 FORALLfacets {
01770 if (!qh_checkflipped(facet, NULL, qh_ALL)) {
01771 trace1((qh ferr, 1031, "qh_initialhull: initial orientation incorrect. Correct all facets\n"));
01772 facet->flipped= False;
01773 FORALLfacets {
01774 facet->toporient ^= (unsigned char)True;
01775 qh_orientoutside(facet);
01776 }
01777 break;
01778 }
01779 }
01780 FORALLfacets {
01781 if (!qh_checkflipped(facet, NULL, !qh_ALL)) {
01782 qh_precision("initial facet is coplanar with interior point");
01783 qh_fprintf(qh ferr, 6154, "qhull precision error: initial facet %d is coplanar with the interior point\n",
01784 facet->id);
01785 qh_errexit(qh_ERRsingular, facet, NULL);
01786 }
01787 FOREACHneighbor_(facet) {
01788 angle= qh_getangle(facet->normal, neighbor->normal);
01789 minimize_( minangle, angle);
01790 }
01791 }
01792 if (minangle < qh_MAXnarrow && !qh NOnarrow) {
01793 realT diff= 1.0 + minangle;
01794
01795 qh NARROWhull= True;
01796 qh_option("_narrow-hull", NULL, &diff);
01797 if (minangle < qh_WARNnarrow && !qh RERUN && qh PRINTprecision)
01798 qh_printhelp_narrowhull(qh ferr, minangle);
01799 }
01800 zzval_(Zprocessed)= qh hull_dim+1;
01801 qh_checkpolygon(qh facet_list);
01802 qh_checkconvex(qh facet_list, qh_DATAfault);
01803 #ifndef qh_NOtrace
01804 if (qh IStracing >= 1) {
01805 qh_fprintf(qh ferr, 8105, "qh_initialhull: simplex constructed, interior point:");
01806 for (k=0; k < qh hull_dim; k++)
01807 qh_fprintf(qh ferr, 8106, " %6.4g", qh interior_point[k]);
01808 qh_fprintf(qh ferr, 8107, "\n");
01809 }
01810 #endif
01811 }
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831 setT *qh_initialvertices(int dim, setT *maxpoints, pointT *points, int numpoints) {
01832 pointT *point, **pointp;
01833 setT *vertices, *simplex, *tested;
01834 realT randr;
01835 int idx, point_i, point_n, k;
01836 boolT nearzero= False;
01837
01838 vertices= qh_settemp(dim + 1);
01839 simplex= qh_settemp(dim+1);
01840 if (qh ALLpoints)
01841 qh_maxsimplex(dim, NULL, points, numpoints, &simplex);
01842 else if (qh RANDOMoutside) {
01843 while (qh_setsize(simplex) != dim+1) {
01844 randr= qh_RANDOMint;
01845 randr= randr/(qh_RANDOMmax+1);
01846 idx= (int)floor(qh num_points * randr);
01847 while (qh_setin(simplex, qh_point(idx))) {
01848 idx++;
01849 idx= idx < qh num_points ? idx : 0;
01850 }
01851 qh_setappend(&simplex, qh_point(idx));
01852 }
01853 }else if (qh hull_dim >= qh_INITIALmax) {
01854 tested= qh_settemp(dim+1);
01855 qh_setappend(&simplex, SETfirst_(maxpoints));
01856 qh_setappend(&simplex, SETsecond_(maxpoints));
01857 qh_maxsimplex(fmin_(qh_INITIALsearch, dim), maxpoints, points, numpoints, &simplex);
01858 k= qh_setsize(simplex);
01859 FOREACHpoint_i_(maxpoints) {
01860 if (point_i & 0x1) {
01861 if (!qh_setin(simplex, point) && !qh_setin(tested, point)){
01862 qh_detsimplex(point, simplex, k, &nearzero);
01863 if (nearzero)
01864 qh_setappend(&tested, point);
01865 else {
01866 qh_setappend(&simplex, point);
01867 if (++k == dim)
01868 break;
01869 }
01870 }
01871 }
01872 }
01873 while (k != dim && (point= (pointT*)qh_setdellast(maxpoints))) {
01874 if (!qh_setin(simplex, point) && !qh_setin(tested, point)){
01875 qh_detsimplex(point, simplex, k, &nearzero);
01876 if (nearzero)
01877 qh_setappend(&tested, point);
01878 else {
01879 qh_setappend(&simplex, point);
01880 k++;
01881 }
01882 }
01883 }
01884 idx= 0;
01885 while (k != dim && (point= qh_point(idx++))) {
01886 if (!qh_setin(simplex, point) && !qh_setin(tested, point)){
01887 qh_detsimplex(point, simplex, k, &nearzero);
01888 if (!nearzero){
01889 qh_setappend(&simplex, point);
01890 k++;
01891 }
01892 }
01893 }
01894 qh_settempfree(&tested);
01895 qh_maxsimplex(dim, maxpoints, points, numpoints, &simplex);
01896 }else
01897 qh_maxsimplex(dim, maxpoints, points, numpoints, &simplex);
01898 FOREACHpoint_(simplex)
01899 qh_setaddnth(&vertices, 0, qh_newvertex(point));
01900 qh_settempfree(&simplex);
01901 return vertices;
01902 }
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914 vertexT *qh_isvertex(pointT *point, setT *vertices) {
01915 vertexT *vertex, **vertexp;
01916
01917 FOREACHvertex_(vertices) {
01918 if (vertex->point == point)
01919 return vertex;
01920 }
01921 return NULL;
01922 }
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959 vertexT *qh_makenewfacets(pointT *point ) {
01960 facetT *visible, *newfacet= NULL, *newfacet2= NULL, *neighbor, **neighborp;
01961 vertexT *apex;
01962 int numnew=0;
01963
01964 qh newfacet_list= qh facet_tail;
01965 qh newvertex_list= qh vertex_tail;
01966 apex= qh_newvertex(point);
01967 qh_appendvertex(apex);
01968 qh visit_id++;
01969 if (!qh ONLYgood)
01970 qh NEWfacets= True;
01971 FORALLvisible_facets {
01972 FOREACHneighbor_(visible)
01973 neighbor->seen= False;
01974 if (visible->ridges) {
01975 visible->visitid= qh visit_id;
01976 newfacet2= qh_makenew_nonsimplicial(visible, apex, &numnew);
01977 }
01978 if (visible->simplicial)
01979 newfacet= qh_makenew_simplicial(visible, apex, &numnew);
01980 if (!qh ONLYgood) {
01981 if (newfacet2)
01982 newfacet= newfacet2;
01983 if (newfacet)
01984 visible->f.replace= newfacet;
01985 else
01986 zinc_(Zinsidevisible);
01987 SETfirst_(visible->neighbors)= NULL;
01988 }
01989 }
01990 trace1((qh ferr, 1032, "qh_makenewfacets: created %d new facets from point p%d to horizon\n",
01991 numnew, qh_pointid(point)));
01992 if (qh IStracing >= 4)
01993 qh_printfacetlist(qh newfacet_list, NULL, qh_ALL);
01994 return apex;
01995 }
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005
02006
02007
02008
02009
02010
02011
02012
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026 #ifndef qh_NOmerge
02027 void qh_matchduplicates(facetT *atfacet, int atskip, int hashsize, int *hashcount) {
02028 boolT same, ismatch;
02029 int hash, scan;
02030 facetT *facet, *newfacet, *maxmatch= NULL, *maxmatch2= NULL, *nextfacet;
02031 int skip, newskip, nextskip= 0, maxskip= 0, maxskip2= 0, makematch;
02032 realT maxdist= -REALmax, mindist, dist2, low, high;
02033
02034 hash= qh_gethash(hashsize, atfacet->vertices, qh hull_dim, 1,
02035 SETelem_(atfacet->vertices, atskip));
02036 trace2((qh ferr, 2046, "qh_matchduplicates: find duplicate matches for f%d skip %d hash %d hashcount %d\n",
02037 atfacet->id, atskip, hash, *hashcount));
02038 for (makematch= 0; makematch < 2; makematch++) {
02039 qh visit_id++;
02040 for (newfacet= atfacet, newskip= atskip; newfacet; newfacet= nextfacet, newskip= nextskip) {
02041 zinc_(Zhashlookup);
02042 nextfacet= NULL;
02043 newfacet->visitid= qh visit_id;
02044 for (scan= hash; (facet= SETelemt_(qh hash_table, scan, facetT));
02045 scan= (++scan >= hashsize ? 0 : scan)) {
02046 if (!facet->dupridge || facet->visitid == qh visit_id)
02047 continue;
02048 zinc_(Zhashtests);
02049 if (qh_matchvertices(1, newfacet->vertices, newskip, facet->vertices, &skip, &same)) {
02050 ismatch= (same == (boolT)(newfacet->toporient ^ facet->toporient));
02051 if (SETelemt_(facet->neighbors, skip, facetT) != qh_DUPLICATEridge) {
02052 if (!makematch) {
02053 qh_fprintf(qh ferr, 6155, "qhull internal error (qh_matchduplicates): missing dupridge at f%d skip %d for new f%d skip %d hash %d\n",
02054 facet->id, skip, newfacet->id, newskip, hash);
02055 qh_errexit2 (qh_ERRqhull, facet, newfacet);
02056 }
02057 }else if (ismatch && makematch) {
02058 if (SETelemt_(newfacet->neighbors, newskip, facetT) == qh_DUPLICATEridge) {
02059 SETelem_(facet->neighbors, skip)= newfacet;
02060 if (newfacet->tricoplanar)
02061 SETelem_(newfacet->neighbors, newskip)= facet;
02062 else
02063 SETelem_(newfacet->neighbors, newskip)= qh_MERGEridge;
02064 *hashcount -= 2;
02065 trace4((qh ferr, 4059, "qh_matchduplicates: duplicate f%d skip %d matched with new f%d skip %d merge\n",
02066 facet->id, skip, newfacet->id, newskip));
02067 }
02068 }else if (ismatch) {
02069 mindist= qh_getdistance(facet, newfacet, &low, &high);
02070 dist2= qh_getdistance(newfacet, facet, &low, &high);
02071 minimize_(mindist, dist2);
02072 if (mindist > maxdist) {
02073 maxdist= mindist;
02074 maxmatch= facet;
02075 maxskip= skip;
02076 maxmatch2= newfacet;
02077 maxskip2= newskip;
02078 }
02079 trace3((qh ferr, 3018, "qh_matchduplicates: duplicate f%d skip %d new f%d skip %d at dist %2.2g, max is now f%d f%d\n",
02080 facet->id, skip, newfacet->id, newskip, mindist,
02081 maxmatch->id, maxmatch2->id));
02082 }else {
02083 nextfacet= facet;
02084 nextskip= skip;
02085 }
02086 }
02087 if (makematch && !facet
02088 && SETelemt_(facet->neighbors, skip, facetT) == qh_DUPLICATEridge) {
02089 qh_fprintf(qh ferr, 6156, "qhull internal error (qh_matchduplicates): no MERGEridge match for duplicate f%d skip %d at hash %d\n",
02090 newfacet->id, newskip, hash);
02091 qh_errexit(qh_ERRqhull, newfacet, NULL);
02092 }
02093 }
02094 }
02095 if (!makematch) {
02096 if (!maxmatch) {
02097 qh_fprintf(qh ferr, 6157, "qhull internal error (qh_matchduplicates): no maximum match at duplicate f%d skip %d at hash %d\n",
02098 atfacet->id, atskip, hash);
02099 qh_errexit(qh_ERRqhull, atfacet, NULL);
02100 }
02101 SETelem_(maxmatch->neighbors, maxskip)= maxmatch2;
02102 SETelem_(maxmatch2->neighbors, maxskip2)= maxmatch;
02103 *hashcount -= 2;
02104 zzinc_(Zmultiridge);
02105 trace0((qh ferr, 25, "qh_matchduplicates: duplicate f%d skip %d matched with new f%d skip %d keep\n",
02106 maxmatch->id, maxskip, maxmatch2->id, maxskip2));
02107 qh_precision("ridge with multiple neighbors");
02108 if (qh IStracing >= 4)
02109 qh_errprint("DUPLICATED/MATCH", maxmatch, maxmatch2, NULL, NULL);
02110 }
02111 }
02112 }
02113
02114
02115
02116
02117
02118
02119
02120
02121
02122
02123
02124
02125
02126
02127
02128
02129
02130
02131
02132
02133
02134
02135
02136
02137 void qh_nearcoplanar(void ) {
02138 facetT *facet;
02139 pointT *point, **pointp;
02140 int numpart;
02141 realT dist, innerplane;
02142
02143 if (!qh KEEPcoplanar && !qh KEEPinside) {
02144 FORALLfacets {
02145 if (facet->coplanarset)
02146 qh_setfree( &facet->coplanarset);
02147 }
02148 }else if (!qh KEEPcoplanar || !qh KEEPinside) {
02149 qh_outerinner(NULL, NULL, &innerplane);
02150 if (qh JOGGLEmax < REALmax/2)
02151 innerplane -= qh JOGGLEmax * sqrt((realT)qh hull_dim);
02152 numpart= 0;
02153 FORALLfacets {
02154 if (facet->coplanarset) {
02155 FOREACHpoint_(facet->coplanarset) {
02156 numpart++;
02157 qh_distplane(point, facet, &dist);
02158 if (dist < innerplane) {
02159 if (!qh KEEPinside)
02160 SETref_(point)= NULL;
02161 }else if (!qh KEEPcoplanar)
02162 SETref_(point)= NULL;
02163 }
02164 qh_setcompact(facet->coplanarset);
02165 }
02166 }
02167 zzadd_(Zcheckpart, numpart);
02168 }
02169 }
02170
02171
02172
02173
02174
02175
02176
02177
02178
02179
02180
02181
02182
02183
02184
02185
02186
02187 vertexT *qh_nearvertex(facetT *facet, pointT *point, realT *bestdistp) {
02188 realT bestdist= REALmax, dist;
02189 vertexT *bestvertex= NULL, *vertex, **vertexp, *apex;
02190 coordT *center;
02191 facetT *neighbor, **neighborp;
02192 setT *vertices;
02193 int dim= qh hull_dim;
02194
02195 if (qh DELAUNAY)
02196 dim--;
02197 if (facet->tricoplanar) {
02198 if (!qh VERTEXneighbors || !facet->center) {
02199 qh_fprintf(qh ferr, 6158, "qhull internal error (qh_nearvertex): qh.VERTEXneighbors and facet->center required for tricoplanar facets\n");
02200 qh_errexit(qh_ERRqhull, facet, NULL);
02201 }
02202 vertices= qh_settemp(qh TEMPsize);
02203 apex= SETfirstt_(facet->vertices, vertexT);
02204 center= facet->center;
02205 FOREACHneighbor_(apex) {
02206 if (neighbor->center == center) {
02207 FOREACHvertex_(neighbor->vertices)
02208 qh_setappend(&vertices, vertex);
02209 }
02210 }
02211 }else
02212 vertices= facet->vertices;
02213 FOREACHvertex_(vertices) {
02214 dist= qh_pointdist(vertex->point, point, -dim);
02215 if (dist < bestdist) {
02216 bestdist= dist;
02217 bestvertex= vertex;
02218 }
02219 }
02220 if (facet->tricoplanar)
02221 qh_settempfree(&vertices);
02222 *bestdistp= sqrt(bestdist);
02223 trace3((qh ferr, 3019, "qh_nearvertex: v%d dist %2.2g for f%d p%d\n",
02224 bestvertex->id, *bestdistp, facet->id, qh_pointid(point)));
02225 return bestvertex;
02226 }
02227
02228
02229
02230
02231
02232
02233
02234
02235
02236
02237
02238
02239 int qh_newhashtable(int newsize) {
02240 int size;
02241
02242 size= ((newsize+1)*qh_HASHfactor) | 0x1;
02243 while (True) {
02244 if (newsize<0 || size<0) {
02245 qh_fprintf(qhmem.ferr, 6236, "qhull error (qh_newhashtable): negative request (%d) or size (%d). Did int overflow due to high-D?\n", newsize, size);
02246 qh_errexit(qhmem_ERRmem, NULL, NULL);
02247 }
02248 if ((size%3) && (size%5))
02249 break;
02250 size += 2;
02251
02252 }
02253 qh hash_table= qh_setnew(size);
02254 qh_setzero(qh hash_table, 0, size);
02255 return size;
02256 }
02257
02258
02259
02260
02261
02262
02263
02264 vertexT *qh_newvertex(pointT *point) {
02265 vertexT *vertex;
02266
02267 zinc_(Ztotvertices);
02268 vertex= (vertexT *)qh_memalloc((int)sizeof(vertexT));
02269 memset((char *) vertex, (size_t)0, sizeof(vertexT));
02270 if (qh vertex_id == 0xFFFFFF) {
02271 qh_fprintf(qh ferr, 6159, "qhull input error: more than %d vertices. ID field overflows and two vertices\n\
02272 may have the same identifier. Vertices not sorted correctly.\n", 0xFFFFFF);
02273 qh_errexit(qh_ERRinput, NULL, NULL);
02274 }
02275 if (qh vertex_id == qh tracevertex_id)
02276 qh tracevertex= vertex;
02277 vertex->id= qh vertex_id++;
02278 vertex->point= point;
02279 vertex->dim= (unsigned char)(qh hull_dim <= MAX_vdim ? qh hull_dim : 0);
02280 trace4((qh ferr, 4060, "qh_newvertex: vertex p%d(v%d) created\n", qh_pointid(vertex->point),
02281 vertex->id));
02282 return(vertex);
02283 }
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304 ridgeT *qh_nextridge3d(ridgeT *atridge, facetT *facet, vertexT **vertexp) {
02305 vertexT *atvertex, *vertex, *othervertex;
02306 ridgeT *ridge, **ridgep;
02307
02308 if ((atridge->top == facet) ^ qh_ORIENTclock)
02309 atvertex= SETsecondt_(atridge->vertices, vertexT);
02310 else
02311 atvertex= SETfirstt_(atridge->vertices, vertexT);
02312 FOREACHridge_(facet->ridges) {
02313 if (ridge == atridge)
02314 continue;
02315 if ((ridge->top == facet) ^ qh_ORIENTclock) {
02316 othervertex= SETsecondt_(ridge->vertices, vertexT);
02317 vertex= SETfirstt_(ridge->vertices, vertexT);
02318 }else {
02319 vertex= SETsecondt_(ridge->vertices, vertexT);
02320 othervertex= SETfirstt_(ridge->vertices, vertexT);
02321 }
02322 if (vertex == atvertex) {
02323 if (vertexp)
02324 *vertexp= othervertex;
02325 return ridge;
02326 }
02327 }
02328 return NULL;
02329 }
02330 #else
02331 void qh_matchduplicates(facetT *atfacet, int atskip, int hashsize, int *hashcount) {
02332 }
02333 ridgeT *qh_nextridge3d(ridgeT *atridge, facetT *facet, vertexT **vertexp) {
02334
02335 return NULL;
02336 }
02337 #endif
02338
02339
02340
02341
02342
02343
02344
02345
02346
02347
02348
02349
02350
02351
02352
02353 void qh_outcoplanar(void ) {
02354 pointT *point, **pointp;
02355 facetT *facet;
02356 realT dist;
02357
02358 trace1((qh ferr, 1033, "qh_outcoplanar: move outsideset to coplanarset for qh NARROWhull\n"));
02359 FORALLfacets {
02360 FOREACHpoint_(facet->outsideset) {
02361 qh num_outside--;
02362 if (qh KEEPcoplanar || qh KEEPnearinside) {
02363 qh_distplane(point, facet, &dist);
02364 zinc_(Zpartition);
02365 qh_partitioncoplanar(point, facet, &dist);
02366 }
02367 }
02368 qh_setfree(&facet->outsideset);
02369 }
02370 }
02371
02372
02373
02374
02375
02376
02377
02378
02379
02380
02381
02382 pointT *qh_point(int id) {
02383
02384 if (id < 0)
02385 return NULL;
02386 if (id < qh num_points)
02387 return qh first_point + id * qh hull_dim;
02388 id -= qh num_points;
02389 if (id < qh_setsize(qh other_points))
02390 return SETelemt_(qh other_points, id, pointT);
02391 return NULL;
02392 }
02393
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403
02404
02405
02406 void qh_point_add(setT *set, pointT *point, void *elem) {
02407 int id, size;
02408
02409 SETreturnsize_(set, size);
02410 if ((id= qh_pointid(point)) < 0)
02411 qh_fprintf(qh ferr, 7067, "qhull internal warning (point_add): unknown point %p id %d\n",
02412 point, id);
02413 else if (id >= size) {
02414 qh_fprintf(qh ferr, 6160, "qhull internal errror(point_add): point p%d is out of bounds(%d)\n",
02415 id, size);
02416 qh_errexit(qh_ERRqhull, NULL, NULL);
02417 }else
02418 SETelem_(set, id)= elem;
02419 }
02420
02421
02422
02423
02424
02425
02426
02427
02428
02429
02430
02431
02432
02433
02434
02435
02436
02437
02438
02439
02440
02441
02442
02443
02444
02445
02446 setT *qh_pointfacet(void ) {
02447 int numpoints= qh num_points + qh_setsize(qh other_points);
02448 setT *facets;
02449 facetT *facet;
02450 vertexT *vertex, **vertexp;
02451 pointT *point, **pointp;
02452
02453 facets= qh_settemp(numpoints);
02454 qh_setzero(facets, 0, numpoints);
02455 qh vertex_visit++;
02456 FORALLfacets {
02457 FOREACHvertex_(facet->vertices) {
02458 if (vertex->visitid != qh vertex_visit) {
02459 vertex->visitid= qh vertex_visit;
02460 qh_point_add(facets, vertex->point, facet);
02461 }
02462 }
02463 FOREACHpoint_(facet->coplanarset)
02464 qh_point_add(facets, point, facet);
02465 FOREACHpoint_(facet->outsideset)
02466 qh_point_add(facets, point, facet);
02467 }
02468 return facets;
02469 }
02470
02471
02472
02473
02474
02475
02476
02477
02478
02479
02480
02481
02482
02483 setT *qh_pointvertex(void ) {
02484 int numpoints= qh num_points + qh_setsize(qh other_points);
02485 setT *vertices;
02486 vertexT *vertex;
02487
02488 vertices= qh_settemp(numpoints);
02489 qh_setzero(vertices, 0, numpoints);
02490 FORALLvertices
02491 qh_point_add(vertices, vertex->point, vertex);
02492 return vertices;
02493 }
02494
02495
02496
02497
02498
02499
02500
02501
02502
02503
02504
02505
02506
02507
02508
02509
02510 void qh_prependfacet(facetT *facet, facetT **facetlist) {
02511 facetT *prevfacet, *list;
02512
02513
02514 trace4((qh ferr, 4061, "qh_prependfacet: prepend f%d before f%d\n",
02515 facet->id, getid_(*facetlist)));
02516 if (!*facetlist)
02517 (*facetlist)= qh facet_tail;
02518 list= *facetlist;
02519 prevfacet= list->previous;
02520 facet->previous= prevfacet;
02521 if (prevfacet)
02522 prevfacet->next= facet;
02523 list->previous= facet;
02524 facet->next= *facetlist;
02525 if (qh facet_list == list)
02526 qh facet_list= facet;
02527 if (qh facet_next == list)
02528 qh facet_next= facet;
02529 *facetlist= facet;
02530 qh num_facets++;
02531 }
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543
02544
02545
02546
02547
02548
02549 void qh_printhashtable(FILE *fp) {
02550 facetT *facet, *neighbor;
02551 int id, facet_i, facet_n, neighbor_i= 0, neighbor_n= 0;
02552 vertexT *vertex, **vertexp;
02553
02554 FOREACHfacet_i_(qh hash_table) {
02555 if (facet) {
02556 FOREACHneighbor_i_(facet) {
02557 if (!neighbor || neighbor == qh_MERGEridge || neighbor == qh_DUPLICATEridge)
02558 break;
02559 }
02560 if (neighbor_i == neighbor_n)
02561 continue;
02562 qh_fprintf(fp, 9283, "hash %d f%d ", facet_i, facet->id);
02563 FOREACHvertex_(facet->vertices)
02564 qh_fprintf(fp, 9284, "v%d ", vertex->id);
02565 qh_fprintf(fp, 9285, "\n neighbors:");
02566 FOREACHneighbor_i_(facet) {
02567 if (neighbor == qh_MERGEridge)
02568 id= -3;
02569 else if (neighbor == qh_DUPLICATEridge)
02570 id= -2;
02571 else
02572 id= getid_(neighbor);
02573 qh_fprintf(fp, 9286, " %d", id);
02574 }
02575 qh_fprintf(fp, 9287, "\n");
02576 }
02577 }
02578 }
02579
02580
02581
02582
02583
02584
02585
02586
02587 void qh_printlists(void) {
02588 facetT *facet;
02589 vertexT *vertex;
02590 int count= 0;
02591
02592 qh_fprintf(qh ferr, 8108, "qh_printlists: facets:");
02593 FORALLfacets {
02594 if (++count % 100 == 0)
02595 qh_fprintf(qh ferr, 8109, "\n ");
02596 qh_fprintf(qh ferr, 8110, " %d", facet->id);
02597 }
02598 qh_fprintf(qh ferr, 8111, "\n new facets %d visible facets %d next facet for qh_addpoint %d\n vertices(new %d):",
02599 getid_(qh newfacet_list), getid_(qh visible_list), getid_(qh facet_next),
02600 getid_(qh newvertex_list));
02601 count = 0;
02602 FORALLvertices {
02603 if (++count % 100 == 0)
02604 qh_fprintf(qh ferr, 8112, "\n ");
02605 qh_fprintf(qh ferr, 8113, " %d", vertex->id);
02606 }
02607 qh_fprintf(qh ferr, 8114, "\n");
02608 }
02609
02610
02611
02612
02613
02614
02615
02616
02617
02618
02619
02620
02621 void qh_resetlists(boolT stats, boolT resetVisible ) {
02622 vertexT *vertex;
02623 facetT *newfacet, *visible;
02624 int totnew=0, totver=0;
02625
02626 if (stats) {
02627 FORALLvertex_(qh newvertex_list)
02628 totver++;
02629 FORALLnew_facets
02630 totnew++;
02631 zadd_(Zvisvertextot, totver);
02632 zmax_(Zvisvertexmax, totver);
02633 zadd_(Znewfacettot, totnew);
02634 zmax_(Znewfacetmax, totnew);
02635 }
02636 FORALLvertex_(qh newvertex_list)
02637 vertex->newlist= False;
02638 qh newvertex_list= NULL;
02639 FORALLnew_facets
02640 newfacet->newfacet= False;
02641 qh newfacet_list= NULL;
02642 if (resetVisible) {
02643 FORALLvisible_facets {
02644 visible->f.replace= NULL;
02645 visible->visible= False;
02646 }
02647 qh num_visible= 0;
02648 }
02649 qh visible_list= NULL;
02650 qh NEWfacets= False;
02651 }
02652
02653
02654
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664
02665
02666
02667
02668
02669
02670
02671 void qh_setvoronoi_all(void) {
02672 facetT *facet;
02673
02674 qh_clearcenters(qh_ASvoronoi);
02675 qh_vertexneighbors();
02676
02677 FORALLfacets {
02678 if (!facet->normal || !facet->upperdelaunay || qh UPPERdelaunay) {
02679 if (!facet->center)
02680 facet->center= qh_facetcenter(facet->vertices);
02681 }
02682 }
02683 }
02684
02685 #ifndef qh_NOmerge
02686
02687
02688
02689
02690
02691
02692
02693
02694
02695
02696
02697
02698
02699
02700
02701
02702
02703 void qh_triangulate(void ) {
02704 facetT *facet, *nextfacet, *owner;
02705 int onlygood= qh ONLYgood;
02706 facetT *neighbor, *visible= NULL, *facet1, *facet2, *new_facet_list= NULL;
02707 facetT *orig_neighbor= NULL, *otherfacet;
02708 vertexT *new_vertex_list= NULL;
02709 mergeT *merge;
02710 mergeType mergetype;
02711 int neighbor_i, neighbor_n;
02712
02713 if (qh hasTriangulation)
02714 return;
02715 trace1((qh ferr, 1034, "qh_triangulate: triangulate non-simplicial facets\n"));
02716 if (qh hull_dim == 2)
02717 return;
02718 if (qh VORONOI) {
02719 qh_clearcenters(qh_ASvoronoi);
02720 qh_vertexneighbors();
02721 }
02722 qh ONLYgood= False;
02723 qh visit_id++;
02724 qh NEWfacets= True;
02725 qh degen_mergeset= qh_settemp(qh TEMPsize);
02726 qh newvertex_list= qh vertex_tail;
02727 for (facet= qh facet_list; facet && facet->next; facet= nextfacet) {
02728 nextfacet= facet->next;
02729 if (facet->visible || facet->simplicial)
02730 continue;
02731
02732 if (!new_facet_list)
02733 new_facet_list= facet;
02734 qh_triangulate_facet(facet, &new_vertex_list);
02735 }
02736 trace2((qh ferr, 2047, "qh_triangulate: delete null facets from f%d -- apex same as second vertex\n", getid_(new_facet_list)));
02737 for (facet= new_facet_list; facet && facet->next; facet= nextfacet) {
02738 nextfacet= facet->next;
02739 if (facet->visible)
02740 continue;
02741 if (facet->ridges) {
02742 if (qh_setsize(facet->ridges) > 0) {
02743 qh_fprintf(qh ferr, 6161, "qhull error (qh_triangulate): ridges still defined for f%d\n", facet->id);
02744 qh_errexit(qh_ERRqhull, facet, NULL);
02745 }
02746 qh_setfree(&facet->ridges);
02747 }
02748 if (SETfirst_(facet->vertices) == SETsecond_(facet->vertices)) {
02749 zinc_(Ztrinull);
02750 qh_triangulate_null(facet);
02751 }
02752 }
02753 trace2((qh ferr, 2048, "qh_triangulate: delete %d or more mirror facets -- same vertices and neighbors\n", qh_setsize(qh degen_mergeset)));
02754 qh visible_list= qh facet_tail;
02755 while ((merge= (mergeT*)qh_setdellast(qh degen_mergeset))) {
02756 facet1= merge->facet1;
02757 facet2= merge->facet2;
02758 mergetype= merge->type;
02759 qh_memfree(merge, (int)sizeof(mergeT));
02760 if (mergetype == MRGmirror) {
02761 zinc_(Ztrimirror);
02762 qh_triangulate_mirror(facet1, facet2);
02763 }
02764 }
02765 qh_settempfree(&qh degen_mergeset);
02766 trace2((qh ferr, 2049, "qh_triangulate: update neighbor lists for vertices from v%d\n", getid_(new_vertex_list)));
02767 qh newvertex_list= new_vertex_list;
02768 qh visible_list= NULL;
02769 qh_updatevertices();
02770 qh_resetlists(False, !qh_RESETvisible );
02771
02772 trace2((qh ferr, 2050, "qh_triangulate: identify degenerate tricoplanar facets from f%d\n", getid_(new_facet_list)));
02773 trace2((qh ferr, 2051, "qh_triangulate: and replace facet->f.triowner with tricoplanar facets that own center, normal, etc.\n"));
02774 FORALLfacet_(new_facet_list) {
02775 if (facet->tricoplanar && !facet->visible) {
02776 FOREACHneighbor_i_(facet) {
02777 if (neighbor_i == 0) {
02778 if (neighbor->tricoplanar)
02779 orig_neighbor= neighbor->f.triowner;
02780 else
02781 orig_neighbor= neighbor;
02782 }else {
02783 if (neighbor->tricoplanar)
02784 otherfacet= neighbor->f.triowner;
02785 else
02786 otherfacet= neighbor;
02787 if (orig_neighbor == otherfacet) {
02788 zinc_(Ztridegen);
02789 facet->degenerate= True;
02790 break;
02791 }
02792 }
02793 }
02794 }
02795 }
02796
02797 trace2((qh ferr, 2052, "qh_triangulate: delete visible facets -- non-simplicial, null, and mirrored facets\n"));
02798 owner= NULL;
02799 visible= NULL;
02800 for (facet= new_facet_list; facet && facet->next; facet= nextfacet) {
02801 nextfacet= facet->next;
02802 if (facet->visible) {
02803 if (facet->tricoplanar) {
02804 qh_delfacet(facet);
02805 qh num_visible--;
02806 }else {
02807 if (visible && !owner) {
02808
02809 trace2((qh ferr, 2053, "qh_triangulate: all tricoplanar facets degenerate for non-simplicial facet f%d\n",
02810 visible->id));
02811 qh_delfacet(visible);
02812 qh num_visible--;
02813 }
02814 visible= facet;
02815 owner= NULL;
02816 }
02817 }else if (facet->tricoplanar) {
02818 if (facet->f.triowner != visible) {
02819 qh_fprintf(qh ferr, 6162, "qhull error (qh_triangulate): tricoplanar facet f%d not owned by its visible, non-simplicial facet f%d\n", facet->id, getid_(visible));
02820 qh_errexit2 (qh_ERRqhull, facet, visible);
02821 }
02822 if (owner)
02823 facet->f.triowner= owner;
02824 else if (!facet->degenerate) {
02825 owner= facet;
02826 nextfacet= visible->next;
02827 facet->keepcentrum= True;
02828 facet->coplanarset= visible->coplanarset;
02829 facet->outsideset= visible->outsideset;
02830 visible->coplanarset= NULL;
02831 visible->outsideset= NULL;
02832 if (!qh TRInormals) {
02833 visible->center= NULL;
02834 visible->normal= NULL;
02835 }
02836 qh_delfacet(visible);
02837 qh num_visible--;
02838 }
02839 }
02840 }
02841 if (visible && !owner) {
02842 trace2((qh ferr, 2054, "qh_triangulate: all tricoplanar facets degenerate for last non-simplicial facet f%d\n",
02843 visible->id));
02844 qh_delfacet(visible);
02845 qh num_visible--;
02846 }
02847 qh NEWfacets= False;
02848 qh ONLYgood= onlygood;
02849 if (qh CHECKfrequently)
02850 qh_checkpolygon(qh facet_list);
02851 qh hasTriangulation= True;
02852 }
02853
02854
02855
02856
02857
02858
02859
02860
02861
02862
02863
02864
02865
02866
02867
02868
02869
02870
02871
02872
02873
02874
02875
02876
02877
02878
02879
02880
02881
02882
02883
02884
02885 void qh_triangulate_facet(facetT *facetA, vertexT **first_vertex) {
02886 facetT *newfacet;
02887 facetT *neighbor, **neighborp;
02888 vertexT *apex;
02889 int numnew=0;
02890
02891 trace3((qh ferr, 3020, "qh_triangulate_facet: triangulate facet f%d\n", facetA->id));
02892
02893 if (qh IStracing >= 4)
02894 qh_printfacet(qh ferr, facetA);
02895 FOREACHneighbor_(facetA) {
02896 neighbor->seen= False;
02897 neighbor->coplanar= False;
02898 }
02899 if (qh CENTERtype == qh_ASvoronoi && !facetA->center
02900 && fabs_(facetA->normal[qh hull_dim -1]) >= qh ANGLEround * qh_ZEROdelaunay) {
02901 facetA->center= qh_facetcenter(facetA->vertices);
02902 }
02903 qh_willdelete(facetA, NULL);
02904 qh newfacet_list= qh facet_tail;
02905 facetA->visitid= qh visit_id;
02906 apex= SETfirstt_(facetA->vertices, vertexT);
02907 qh_makenew_nonsimplicial(facetA, apex, &numnew);
02908 SETfirst_(facetA->neighbors)= NULL;
02909 FORALLnew_facets {
02910 newfacet->tricoplanar= True;
02911 newfacet->f.trivisible= facetA;
02912 newfacet->degenerate= False;
02913 newfacet->upperdelaunay= facetA->upperdelaunay;
02914 newfacet->good= facetA->good;
02915 if (qh TRInormals) {
02916 newfacet->keepcentrum= True;
02917 newfacet->normal= qh_copypoints(facetA->normal, 1, qh hull_dim);
02918 if (qh CENTERtype == qh_AScentrum)
02919 newfacet->center= qh_getcentrum(newfacet);
02920 else
02921 newfacet->center= qh_copypoints(facetA->center, 1, qh hull_dim);
02922 }else {
02923 newfacet->keepcentrum= False;
02924 newfacet->normal= facetA->normal;
02925 newfacet->center= facetA->center;
02926 }
02927 newfacet->offset= facetA->offset;
02928 #if qh_MAXoutside
02929 newfacet->maxoutside= facetA->maxoutside;
02930 #endif
02931 }
02932 qh_matchnewfacets();
02933 zinc_(Ztricoplanar);
02934 zadd_(Ztricoplanartot, numnew);
02935 zmax_(Ztricoplanarmax, numnew);
02936 qh visible_list= NULL;
02937 if (!(*first_vertex))
02938 (*first_vertex)= qh newvertex_list;
02939 qh newvertex_list= NULL;
02940 qh_updatevertices();
02941 qh_resetlists(False, !qh_RESETvisible );
02942 }
02943
02944
02945
02946
02947
02948
02949
02950
02951
02952
02953
02954 void qh_triangulate_link(facetT *oldfacetA, facetT *facetA, facetT *oldfacetB, facetT *facetB) {
02955 int errmirror= False;
02956
02957 trace3((qh ferr, 3021, "qh_triangulate_link: relink old facets f%d and f%d between neighbors f%d and f%d\n",
02958 oldfacetA->id, oldfacetB->id, facetA->id, facetB->id));
02959 if (qh_setin(facetA->neighbors, facetB)) {
02960 if (!qh_setin(facetB->neighbors, facetA))
02961 errmirror= True;
02962 else
02963 qh_appendmergeset(facetA, facetB, MRGmirror, NULL);
02964 }else if (qh_setin(facetB->neighbors, facetA))
02965 errmirror= True;
02966 if (errmirror) {
02967 qh_fprintf(qh ferr, 6163, "qhull error (qh_triangulate_link): mirror facets f%d and f%d do not match for old facets f%d and f%d\n",
02968 facetA->id, facetB->id, oldfacetA->id, oldfacetB->id);
02969 qh_errexit2 (qh_ERRqhull, facetA, facetB);
02970 }
02971 qh_setreplace(facetB->neighbors, oldfacetB, facetA);
02972 qh_setreplace(facetA->neighbors, oldfacetA, facetB);
02973 }
02974
02975
02976
02977
02978
02979
02980
02981
02982
02983
02984
02985 void qh_triangulate_mirror(facetT *facetA, facetT *facetB) {
02986 facetT *neighbor, *neighborB;
02987 int neighbor_i, neighbor_n;
02988
02989 trace3((qh ferr, 3022, "qh_triangulate_mirror: delete mirrored facets f%d and f%d\n",
02990 facetA->id, facetB->id));
02991 FOREACHneighbor_i_(facetA) {
02992 neighborB= SETelemt_(facetB->neighbors, neighbor_i, facetT);
02993 if (neighbor == neighborB)
02994 continue;
02995 qh_triangulate_link(facetA, neighbor, facetB, neighborB);
02996 }
02997 qh_willdelete(facetA, NULL);
02998 qh_willdelete(facetB, NULL);
02999 }
03000
03001
03002
03003
03004
03005
03006
03007
03008
03009
03010
03011
03012
03013
03014 void qh_triangulate_null(facetT *facetA) {
03015 facetT *neighbor, *otherfacet;
03016
03017 trace3((qh ferr, 3023, "qh_triangulate_null: delete null facet f%d\n", facetA->id));
03018 neighbor= SETfirstt_(facetA->neighbors, facetT);
03019 otherfacet= SETsecondt_(facetA->neighbors, facetT);
03020 qh_triangulate_link(facetA, neighbor, facetA, otherfacet);
03021 qh_willdelete(facetA, NULL);
03022 }
03023
03024 #else
03025 void qh_triangulate(void) {
03026 }
03027 #endif
03028
03029
03030
03031
03032
03033
03034
03035
03036
03037
03038
03039
03040
03041
03042 void qh_vertexintersect(setT **vertexsetA,setT *vertexsetB) {
03043 setT *intersection;
03044
03045 intersection= qh_vertexintersect_new(*vertexsetA, vertexsetB);
03046 qh_settempfree(vertexsetA);
03047 *vertexsetA= intersection;
03048 qh_settemppush(intersection);
03049 }
03050
03051
03052
03053
03054
03055
03056
03057
03058
03059
03060 setT *qh_vertexintersect_new(setT *vertexsetA,setT *vertexsetB) {
03061 setT *intersection= qh_setnew(qh hull_dim - 1);
03062 vertexT **vertexA= SETaddr_(vertexsetA, vertexT);
03063 vertexT **vertexB= SETaddr_(vertexsetB, vertexT);
03064
03065 while (*vertexA && *vertexB) {
03066 if (*vertexA == *vertexB) {
03067 qh_setappend(&intersection, *vertexA);
03068 vertexA++; vertexB++;
03069 }else {
03070 if ((*vertexA)->id > (*vertexB)->id)
03071 vertexA++;
03072 else
03073 vertexB++;
03074 }
03075 }
03076 return intersection;
03077 }
03078
03079
03080
03081
03082
03083
03084
03085
03086
03087
03088
03089
03090
03091
03092
03093
03094
03095
03096
03097
03098
03099 void qh_vertexneighbors(void ) {
03100 facetT *facet;
03101 vertexT *vertex, **vertexp;
03102
03103 if (qh VERTEXneighbors)
03104 return;
03105 trace1((qh ferr, 1035, "qh_vertexneighbors: determing neighboring facets for each vertex\n"));
03106 qh vertex_visit++;
03107 FORALLfacets {
03108 if (facet->visible)
03109 continue;
03110 FOREACHvertex_(facet->vertices) {
03111 if (vertex->visitid != qh vertex_visit) {
03112 vertex->visitid= qh vertex_visit;
03113 vertex->neighbors= qh_setnew(qh hull_dim);
03114 }
03115 qh_setappend(&vertex->neighbors, facet);
03116 }
03117 }
03118 qh VERTEXneighbors= True;
03119 }
03120
03121
03122
03123
03124
03125
03126
03127
03128
03129
03130
03131 boolT qh_vertexsubset(setT *vertexsetA, setT *vertexsetB) {
03132 vertexT **vertexA= (vertexT **) SETaddr_(vertexsetA, vertexT);
03133 vertexT **vertexB= (vertexT **) SETaddr_(vertexsetB, vertexT);
03134
03135 while (True) {
03136 if (!*vertexA)
03137 return True;
03138 if (!*vertexB)
03139 return False;
03140 if ((*vertexA)->id > (*vertexB)->id)
03141 return False;
03142 if (*vertexA == *vertexB)
03143 vertexA++;
03144 vertexB++;
03145 }
03146 return False;
03147 }
