00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include "qhull_a.h"
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 void qh_addhash (void* newelem, setT *hashtable, int hashsize, unsigned hash) {
00025 int scan;
00026 void *elem;
00027
00028 for (scan= (int)hash; (elem= SETelem_(hashtable, scan));
00029 scan= (++scan >= hashsize ? 0 : scan)) {
00030 if (elem == newelem)
00031 break;
00032 }
00033
00034 if (!elem)
00035 SETelem_(hashtable, scan)= newelem;
00036 }
00037
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 void qh_check_bestdist (void) {
00063 boolT waserror= False, isoutside, unassigned;
00064 facetT *facet, *bestfacet, *errfacet1= NULL, *errfacet2= NULL;
00065 facetT *facetlist;
00066 realT dist, maxoutside, maxdist= -REALmax;
00067 pointT *point;
00068 int numpart, facet_i, facet_n, notgood= 0, notverified= 0;
00069 setT *facets;
00070
00071 trace1((qh ferr, "qh_check_bestdist: check points below nearest facet. Facet_list f%d\n",
00072 qh facet_list->id));
00073 maxoutside= qh_maxouter();
00074 maxoutside += qh DISTround;
00075
00076 trace1((qh ferr, "qh_check_bestdist: check that all points are within %2.2g of best facet\n", maxoutside));
00077 facets= qh_pointfacet ();
00078 if (!qh_QUICKhelp && qh PRINTprecision)
00079 fprintf (qh ferr, "\n\
00080 qhull output completed. Verifying that %d points are\n\
00081 below %2.2g of the nearest %sfacet.\n",
00082 qh_setsize(facets), maxoutside, (qh ONLYgood ? "good " : ""));
00083 FOREACHfacet_i_(facets) {
00084 if (facet)
00085 unassigned= False;
00086 else {
00087 unassigned= True;
00088 facet= qh facet_list;
00089 }
00090 point= qh_point(facet_i);
00091 if (point == qh GOODpointp)
00092 continue;
00093 bestfacet= qh_findbest (point, facet, qh_ALL, False, !qh_NOupper,
00094 &dist, &isoutside, &numpart);
00095
00096 maximize_(maxdist, dist);
00097 if (dist > maxoutside) {
00098 if (qh ONLYgood && !bestfacet->good
00099 && !((bestfacet= qh_findgooddist (point, bestfacet, &dist, &facetlist))
00100 && dist > maxoutside))
00101 notgood++;
00102 else {
00103 waserror= True;
00104 fprintf(qh ferr, "qhull precision error: point p%d is outside facet f%d, distance= %6.8g maxoutside= %6.8g\n",
00105 facet_i, bestfacet->id, dist, maxoutside);
00106 errfacet2= errfacet1;
00107 errfacet1= bestfacet;
00108 }
00109 }else if (unassigned && dist < -qh MAXcoplanar)
00110 notverified++;
00111 }
00112 qh_settempfree (&facets);
00113 if (notverified && !qh DELAUNAY && !qh_QUICKhelp && qh PRINTprecision)
00114 fprintf(qh ferr, "\n%d points were well inside the hull. If the hull contains\n\
00115 a lens-shaped component, these points were not verified. Use\n\
00116 options 'Qci Tv' to verify all points.\n", notverified);
00117 if (maxdist > qh outside_err) {
00118 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",
00119 maxdist, qh outside_err);
00120 qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);
00121 }else if (waserror && qh outside_err > REALmax/2)
00122 qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);
00123 else if (waserror)
00124 ;
00125 trace0((qh ferr, "qh_check_bestdist: max distance outside %2.2g\n", maxdist));
00126 }
00127
00128
00129
00130
00131
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 #ifndef qh_NOmerge
00161 void qh_check_maxout (void) {
00162 facetT *facet, *bestfacet, *neighbor, **neighborp, *facetlist;
00163 realT dist, maxoutside, minvertex;
00164 pointT *point;
00165 int numpart, facet_i, facet_n, notgood= 0;
00166 setT *facets, *vertices;
00167 vertexT *vertex;
00168
00169 trace1((qh ferr, "qh_check_maxout: check and update maxoutside for each facet.\n"));
00170 maxoutside= minvertex= 0;
00171 if (qh VERTEXneighbors
00172 && (qh PRINTsummary || qh KEEPinside || qh KEEPcoplanar
00173 || qh TRACElevel || qh PRINTstatistics
00174 || qh PRINTout[0] == qh_PRINTsummary || qh PRINTout[0] == qh_PRINTnone)) {
00175 trace1((qh ferr, "qh_check_maxout: determine actual maxoutside and minvertex\n"));
00176 vertices= qh_pointvertex ();
00177 FORALLvertices {
00178 FOREACHneighbor_(vertex) {
00179 zinc_(Zdistvertex);
00180 qh_distplane (vertex->point, neighbor, &dist);
00181 minimize_(minvertex, dist);
00182 if (-dist > qh TRACEdist || dist > qh TRACEdist
00183 || neighbor == qh tracefacet || vertex == qh tracevertex)
00184 fprintf (qh ferr, "qh_check_maxout: p%d (v%d) is %.2g from f%d\n",
00185 qh_pointid (vertex->point), vertex->id, dist, neighbor->id);
00186 }
00187 }
00188 if (qh MERGING) {
00189 wmin_(Wminvertex, qh min_vertex);
00190 }
00191 qh min_vertex= minvertex;
00192 qh_settempfree (&vertices);
00193 }
00194 facets= qh_pointfacet ();
00195 FOREACHfacet_i_(facets) {
00196 if (facet) {
00197 point= qh_point(facet_i);
00198 if (point == qh GOODpointp)
00199 continue;
00200 zinc_(Ztotcheck);
00201 bestfacet= qh_findbest (point, facet, qh_ALL, False, !qh_NOupper,
00202 &dist, NULL, &numpart);
00203 zadd_(Zcheckpart, numpart);
00204 if (bestfacet && dist > maxoutside) {
00205 if (qh ONLYgood && !bestfacet->good
00206 && !((bestfacet= qh_findgooddist (point, bestfacet, &dist, &facetlist))
00207 && dist > maxoutside))
00208 notgood++;
00209 else
00210 maxoutside= dist;
00211 }
00212 if (dist > qh TRACEdist || (bestfacet && bestfacet == qh tracefacet))
00213 fprintf (qh ferr, "qh_check_maxout: p%d is %.2g above f%d\n",
00214 qh_pointid (point), dist, bestfacet->id);
00215 }
00216 }
00217 qh_settempfree (&facets);
00218 wval_(Wmaxout)= maxoutside - qh max_outside;
00219 wmax_(Wmaxoutside, qh max_outside);
00220 qh max_outside= maxoutside;
00221 qh_nearcoplanar ();
00222 qh maxoutdone= True;
00223 trace1((qh ferr, "qh_check_maxout: maxoutside %2.2g, min_vertex %2.2g, outside of not good %d\n",
00224 maxoutside, qh min_vertex, notgood));
00225 }
00226 #else
00227 void qh_check_maxout (void) {
00228 }
00229 #endif
00230
00231
00232
00233
00234
00235
00236
00237 void qh_check_output (void) {
00238 int i;
00239
00240 if (qh STOPcone)
00241 return;
00242 if (qh VERIFYoutput | qh IStracing | qh CHECKfrequently) {
00243 qh_checkpolygon (qh facet_list);
00244 qh_checkflipped_all (qh facet_list);
00245 qh_checkconvex (qh facet_list, qh_ALGORITHMfault);
00246 }else if (!qh MERGING && qh_newstats (qhstat precision, &i)) {
00247 qh_checkflipped_all (qh facet_list);
00248 qh_checkconvex (qh facet_list, qh_ALGORITHMfault);
00249 }
00250 }
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260 void qh_check_point (pointT *point, facetT *facet, realT *maxoutside, realT *maxdist, facetT **errfacet1, facetT **errfacet2) {
00261 realT dist;
00262
00263
00264 qh_distplane(point, facet, &dist);
00265 if (dist > *maxoutside) {
00266 *errfacet2= *errfacet1;
00267 *errfacet1= facet;
00268 fprintf(qh ferr, "qhull precision error: point p%d is outside facet f%d, distance= %6.8g maxoutside= %6.8g\n",
00269 qh_pointid(point), facet->id, dist, *maxoutside);
00270 }
00271 maximize_(*maxdist, dist);
00272 }
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 void qh_check_points (void) {
00298 facetT *facet, *errfacet1= NULL, *errfacet2= NULL;
00299 realT total, maxoutside, maxdist= -REALmax;
00300 pointT *point, **pointp, *pointtemp;
00301 boolT testouter;
00302
00303 maxoutside= qh_maxouter();
00304 maxoutside += qh DISTround;
00305
00306 trace1((qh ferr, "qh_check_points: check all points below %2.2g of all facet planes\n",
00307 maxoutside));
00308 if (qh num_good)
00309 total= (float) qh num_good * qh num_points;
00310 else
00311 total= (float) qh num_facets * qh num_points;
00312 if (total >= qh_VERIFYdirect && !qh maxoutdone) {
00313 if (!qh_QUICKhelp && qh SKIPcheckmax && qh MERGING)
00314 fprintf (qh ferr, "\n\
00315 qhull input warning: merging without checking outer planes ('Q5').\n\
00316 Verify may report that a point is outside of a facet.\n");
00317 qh_check_bestdist();
00318 }else {
00319 if (qh_MAXoutside && qh maxoutdone)
00320 testouter= True;
00321 else
00322 testouter= False;
00323 if (!qh_QUICKhelp) {
00324 if (qh MERGEexact)
00325 fprintf (qh ferr, "\n\
00326 qhull input warning: exact merge ('Qx'). Verify may report that a point\n\
00327 is outside of a facet. See qh-optq.htm#Qx\n");
00328 else if (qh SKIPcheckmax || qh NOnearinside)
00329 fprintf (qh ferr, "\n\
00330 qhull input warning: no outer plane check ('Q5') or no processing of\n\
00331 near-inside points ('Q8'). Verify may report that a point is outside\n\
00332 of a facet.\n");
00333 }
00334 if (qh PRINTprecision) {
00335 if (testouter)
00336 fprintf (qh ferr, "\n\
00337 Output completed. Verifying that all points are below outer planes of\n\
00338 all %sfacets. Will make %2.0f distance computations.\n",
00339 (qh ONLYgood ? "good " : ""), total);
00340 else
00341 fprintf (qh ferr, "\n\
00342 Output completed. Verifying that all points are below %2.2g of\n\
00343 all %sfacets. Will make %2.0f distance computations.\n",
00344 maxoutside, (qh ONLYgood ? "good " : ""), total);
00345 }
00346 FORALLfacets {
00347 if (!facet->good && qh ONLYgood)
00348 continue;
00349 if (facet->flipped)
00350 continue;
00351 if (testouter) {
00352 #if qh_MAXoutside
00353 maxoutside= facet->maxoutside + 2* qh DISTround;
00354
00355 #endif
00356 }
00357 FORALLpoints {
00358 if (point != qh GOODpointp)
00359 qh_check_point (point, facet, &maxoutside, &maxdist, &errfacet1, &errfacet2);
00360 }
00361 FOREACHpoint_(qh other_points) {
00362 if (point != qh GOODpointp)
00363 qh_check_point (point, facet, &maxoutside, &maxdist, &errfacet1, &errfacet2);
00364 }
00365 }
00366 if (maxdist > qh outside_err) {
00367 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",
00368 maxdist, qh outside_err );
00369 qh_errexit2( qh_ERRprec, errfacet1, errfacet2 );
00370 }else if (errfacet1 && qh outside_err > REALmax/2)
00371 qh_errexit2( qh_ERRprec, errfacet1, errfacet2 );
00372 else if (errfacet1)
00373 ;
00374 trace0((qh ferr, "qh_check_points: max distance outside %2.2g\n", maxdist));
00375 }
00376 }
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407 void qh_checkconvex(facetT *facetlist, int fault) {
00408 facetT *facet, *neighbor, **neighborp, *errfacet1=NULL, *errfacet2=NULL;
00409 vertexT *vertex;
00410 realT dist;
00411 pointT *centrum;
00412 boolT waserror= False, tempcentrum= False, allsimplicial;
00413 int neighbor_i;
00414
00415 trace1((qh ferr, "qh_checkconvex: check all ridges are convex\n"));
00416 if (!qh RERUN) {
00417 zzval_(Zconcaveridges)= 0;
00418 zzval_(Zcoplanarridges)= 0;
00419 }
00420 FORALLfacet_(facetlist) {
00421 if (facet->flipped) {
00422 qh_precision ("flipped facet");
00423 fprintf (qh ferr, "qhull precision error: f%d is flipped (interior point is outside)\n",
00424 facet->id);
00425 errfacet1= facet;
00426 waserror= True;
00427 continue;
00428 }
00429 if (qh MERGING && (!qh ZEROcentrum || !facet->simplicial))
00430 allsimplicial= False;
00431 else {
00432 allsimplicial= True;
00433 neighbor_i= 0;
00434 FOREACHneighbor_(facet) {
00435 vertex= SETelemt_(facet->vertices, neighbor_i++, vertexT);
00436 if (!neighbor->simplicial) {
00437 allsimplicial= False;
00438 continue;
00439 }
00440 qh_distplane (vertex->point, neighbor, &dist);
00441 if (dist > -qh DISTround) {
00442 if (fault == qh_DATAfault) {
00443 qh_precision ("coplanar or concave ridge");
00444 fprintf (qh ferr, "qhull precision error: initial simplex is not convex. Distance=%.2g\n", dist);
00445 qh_errexit(qh_ERRsingular, NULL, NULL);
00446 }
00447 if (dist > qh DISTround) {
00448 zzinc_(Zconcaveridges);
00449 qh_precision ("concave ridge");
00450 fprintf (qh ferr, "qhull precision error: f%d is concave to f%d, since p%d (v%d) is %6.4g above\n",
00451 facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist);
00452 errfacet1= facet;
00453 errfacet2= neighbor;
00454 waserror= True;
00455 }else if (qh ZEROcentrum) {
00456 if (dist > 0) {
00457 zzinc_(Zcoplanarridges);
00458 qh_precision ("coplanar ridge");
00459 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",
00460 facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist);
00461 errfacet1= facet;
00462 errfacet2= neighbor;
00463 waserror= True;
00464 }
00465 }else {
00466 zzinc_(Zcoplanarridges);
00467 qh_precision ("coplanar ridge");
00468 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",
00469 facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist, qh furthest_id));
00470 }
00471 }
00472 }
00473 }
00474 if (!allsimplicial) {
00475 if (qh CENTERtype == qh_AScentrum) {
00476 if (!facet->center)
00477 facet->center= qh_getcentrum (facet);
00478 centrum= facet->center;
00479 }else {
00480 centrum= qh_getcentrum(facet);
00481 tempcentrum= True;
00482 }
00483 FOREACHneighbor_(facet) {
00484 if (qh ZEROcentrum && facet->simplicial && neighbor->simplicial)
00485 continue;
00486 zzinc_(Zdistconvex);
00487 qh_distplane (centrum, neighbor, &dist);
00488 if (dist > qh DISTround) {
00489 zzinc_(Zconcaveridges);
00490 qh_precision ("concave ridge");
00491 fprintf (qh ferr, "qhull precision error: f%d is concave to f%d. Centrum of f%d is %6.4g above f%d\n",
00492 facet->id, neighbor->id, facet->id, dist, neighbor->id);
00493 errfacet1= facet;
00494 errfacet2= neighbor;
00495 waserror= True;
00496 }else if (dist >= 0.0) {
00497
00498 zzinc_(Zcoplanarridges);
00499 qh_precision ("coplanar ridge");
00500 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",
00501 facet->id, neighbor->id, facet->id, dist, neighbor->id);
00502 errfacet1= facet;
00503 errfacet2= neighbor;
00504 waserror= True;
00505 }
00506 }
00507 if (tempcentrum)
00508 qh_memfree(centrum, qh normal_size);
00509 }
00510 }
00511 if (waserror && !qh FORCEoutput)
00512 qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);
00513 }
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553 void qh_checkfacet(facetT *facet, boolT newmerge, boolT *waserrorp) {
00554 facetT *neighbor, **neighborp, *errother=NULL;
00555 ridgeT *ridge, **ridgep, *errridge= NULL, *ridge2;
00556 vertexT *vertex, **vertexp;
00557 unsigned previousid= INT_MAX;
00558 int numneighbors, numvertices, numridges=0, numRvertices=0;
00559 boolT waserror= False;
00560 int skipA, skipB, ridge_i, ridge_n, i;
00561 setT *intersection;
00562
00563 if (facet->visible) {
00564 fprintf (qh ferr, "qhull internal error (qh_checkfacet): facet f%d is on the visible_list\n",
00565 facet->id);
00566 qh_errexit (qh_ERRqhull, facet, NULL);
00567 }
00568 if (!facet->normal) {
00569 fprintf (qh ferr, "qhull internal error (qh_checkfacet): facet f%d does not have a normal\n",
00570 facet->id);
00571 waserror= True;
00572 }
00573 qh_setcheck (facet->vertices, "vertices for f", facet->id);
00574 qh_setcheck (facet->ridges, "ridges for f", facet->id);
00575 qh_setcheck (facet->outsideset, "outsideset for f", facet->id);
00576 qh_setcheck (facet->coplanarset, "coplanarset for f", facet->id);
00577 qh_setcheck (facet->neighbors, "neighbors for f", facet->id);
00578 FOREACHvertex_(facet->vertices) {
00579 if (vertex->deleted) {
00580 fprintf(qh ferr, "qhull internal error (qh_checkfacet): deleted vertex v%d in f%d\n", vertex->id, facet->id);
00581 qh_errprint ("ERRONEOUS", NULL, NULL, NULL, vertex);
00582 waserror= True;
00583 }
00584 if (vertex->id >= previousid) {
00585 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);
00586 waserror= True;
00587 break;
00588 }
00589 previousid= vertex->id;
00590 }
00591 numneighbors= qh_setsize(facet->neighbors);
00592 numvertices= qh_setsize(facet->vertices);
00593 numridges= qh_setsize(facet->ridges);
00594 if (facet->simplicial) {
00595 if (numvertices+numneighbors != 2*qh hull_dim
00596 && !facet->degenerate && !facet->redundant) {
00597 fprintf(qh ferr, "qhull internal error (qh_checkfacet): for simplicial facet f%d, #vertices %d + #neighbors %d != 2*qh hull_dim\n",
00598 facet->id, numvertices, numneighbors);
00599 qh_setprint (qh ferr, "", facet->neighbors);
00600 waserror= True;
00601 }
00602 }else {
00603 if (!newmerge
00604 &&(numvertices < qh hull_dim || numneighbors < qh hull_dim)
00605 && !facet->degenerate && !facet->redundant) {
00606 fprintf(qh ferr, "qhull internal error (qh_checkfacet): for facet f%d, #vertices %d or #neighbors %d < qh hull_dim\n",
00607 facet->id, numvertices, numneighbors);
00608 waserror= True;
00609 }
00610 if (numridges < numneighbors
00611 ||(qh hull_dim == 3 && numvertices != numridges && !qh NEWfacets)
00612 ||(qh hull_dim == 2 && numridges + numvertices + numneighbors != 6)) {
00613 if (!facet->degenerate && !facet->redundant) {
00614 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",
00615 facet->id, numridges, numneighbors, numvertices);
00616 waserror= True;
00617 }
00618 }
00619 }
00620 FOREACHneighbor_(facet) {
00621 if (neighbor == qh_MERGEridge || neighbor == qh_DUPLICATEridge) {
00622 fprintf(qh ferr, "qhull internal error (qh_checkfacet): facet f%d still has a MERGE or DUP neighbor\n", facet->id);
00623 qh_errexit (qh_ERRqhull, facet, NULL);
00624 }
00625 neighbor->seen= True;
00626 }
00627 FOREACHneighbor_(facet) {
00628 if (!qh_setin(neighbor->neighbors, facet)) {
00629 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",
00630 facet->id, neighbor->id, neighbor->id, facet->id);
00631 errother= neighbor;
00632 waserror= True;
00633 }
00634 if (!neighbor->seen) {
00635 fprintf(qh ferr, "qhull internal error (qh_checkfacet): facet f%d has a duplicate neighbor f%d\n",
00636 facet->id, neighbor->id);
00637 errother= neighbor;
00638 waserror= True;
00639 }
00640 neighbor->seen= False;
00641 }
00642 FOREACHridge_(facet->ridges) {
00643 qh_setcheck (ridge->vertices, "vertices for r", ridge->id);
00644 ridge->seen= False;
00645 }
00646 FOREACHridge_(facet->ridges) {
00647 if (ridge->seen) {
00648 fprintf(qh ferr, "qhull internal error (qh_checkfacet): facet f%d has a duplicate ridge r%d\n",
00649 facet->id, ridge->id);
00650 errridge= ridge;
00651 waserror= True;
00652 }
00653 ridge->seen= True;
00654 numRvertices= qh_setsize(ridge->vertices);
00655 if (numRvertices != qh hull_dim - 1) {
00656 fprintf(qh ferr, "qhull internal error (qh_checkfacet): ridge between f%d and f%d has %d vertices\n",
00657 ridge->top->id, ridge->bottom->id, numRvertices);
00658 errridge= ridge;
00659 waserror= True;
00660 }
00661 neighbor= otherfacet_(ridge, facet);
00662 neighbor->seen= True;
00663 if (!qh_setin(facet->neighbors, neighbor)) {
00664 fprintf(qh ferr, "qhull internal error (qh_checkfacet): for facet f%d, neighbor f%d of ridge r%d not in facet\n",
00665 facet->id, neighbor->id, ridge->id);
00666 errridge= ridge;
00667 waserror= True;
00668 }
00669 }
00670 if (!facet->simplicial) {
00671 FOREACHneighbor_(facet) {
00672 if (!neighbor->seen) {
00673 fprintf(qh ferr, "qhull internal error (qh_checkfacet): facet f%d does not have a ridge for neighbor f%d\n",
00674 facet->id, neighbor->id);
00675 errother= neighbor;
00676 waserror= True;
00677 }
00678 intersection= qh_vertexintersect_new(facet->vertices, neighbor->vertices);
00679 qh_settemppush (intersection);
00680 FOREACHvertex_(facet->vertices) {
00681 vertex->seen= False;
00682 vertex->seen2= False;
00683 }
00684 FOREACHvertex_(intersection)
00685 vertex->seen= True;
00686 FOREACHridge_(facet->ridges) {
00687 if (neighbor != otherfacet_(ridge, facet))
00688 continue;
00689 FOREACHvertex_(ridge->vertices) {
00690 if (!vertex->seen) {
00691 fprintf (qh ferr, "qhull internal error (qh_checkfacet): vertex v%d in r%d not in f%d intersect f%d\n",
00692 vertex->id, ridge->id, facet->id, neighbor->id);
00693 qh_errexit (qh_ERRqhull, facet, ridge);
00694 }
00695 vertex->seen2= True;
00696 }
00697 }
00698 if (!newmerge) {
00699 FOREACHvertex_(intersection) {
00700 if (!vertex->seen2) {
00701 if (qh IStracing >=3 || !qh MERGING) {
00702 fprintf (qh ferr, "qhull precision error (qh_checkfacet): vertex v%d in f%d intersect f%d but\n\
00703 not in a ridge. This is ok under merging. Last point was p%d\n",
00704 vertex->id, facet->id, neighbor->id, qh furthest_id);
00705 if (!qh FORCEoutput && !qh MERGING) {
00706 qh_errprint ("ERRONEOUS", facet, neighbor, NULL, vertex);
00707 if (!qh MERGING)
00708 qh_errexit (qh_ERRqhull, NULL, NULL);
00709 }
00710 }
00711 }
00712 }
00713 }
00714 qh_settempfree (&intersection);
00715 }
00716 }else {
00717 FOREACHneighbor_(facet) {
00718 if (neighbor->simplicial) {
00719 skipA= SETindex_(facet->neighbors, neighbor);
00720 skipB= qh_setindex (neighbor->neighbors, facet);
00721 if (!qh_setequal_skip (facet->vertices, skipA, neighbor->vertices, skipB)) {
00722 fprintf (qh ferr, "qhull internal error (qh_checkfacet): facet f%d skip %d and neighbor f%d skip %d do not match \n",
00723 facet->id, skipA, neighbor->id, skipB);
00724 errother= neighbor;
00725 waserror= True;
00726 }
00727 }
00728 }
00729 }
00730 if (qh hull_dim < 5 && (qh IStracing > 2 || qh CHECKfrequently)) {
00731 FOREACHridge_i_(facet->ridges) {
00732 for (i= ridge_i+1; i < ridge_n; i++) {
00733 ridge2= SETelemt_(facet->ridges, i, ridgeT);
00734 if (qh_setequal (ridge->vertices, ridge2->vertices)) {
00735 fprintf (qh ferr, "qh_checkfacet: ridges r%d and r%d have the same vertices\n",
00736 ridge->id, ridge2->id);
00737 errridge= ridge;
00738 waserror= True;
00739 }
00740 }
00741 }
00742 }
00743 if (waserror) {
00744 qh_errprint("ERRONEOUS", facet, errother, errridge, NULL);
00745 *waserrorp= True;
00746 }
00747 }
00748
00749
00750
00751
00752
00753
00754
00755
00756 void qh_checkflipped_all (facetT *facetlist) {
00757 facetT *facet;
00758 boolT waserror= False;
00759 realT dist;
00760
00761 if (facetlist == qh facet_list)
00762 zzval_(Zflippedfacets)= 0;
00763 FORALLfacet_(facetlist) {
00764 if (facet->normal && !qh_checkflipped (facet, &dist, !qh_ALL)) {
00765 fprintf(qh ferr, "qhull precision error: facet f%d is flipped, distance= %6.12g\n",
00766 facet->id, dist);
00767 if (!qh FORCEoutput) {
00768 qh_errprint("ERRONEOUS", facet, NULL, NULL, NULL);
00769 waserror= True;
00770 }
00771 }
00772 }
00773 if (waserror) {
00774 fprintf (qh ferr, "\n\
00775 A flipped facet occurs when its distance to the interior point is\n\
00776 greater than %2.2g, the maximum roundoff error.\n", -qh DISTround);
00777 qh_errexit(qh_ERRprec, NULL, NULL);
00778 }
00779 }
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803 void qh_checkpolygon(facetT *facetlist) {
00804 facetT *facet;
00805 vertexT *vertex, **vertexp, *vertexlist;
00806 int numfacets= 0, numvertices= 0, numridges= 0;
00807 int totvneighbors= 0, totvertices= 0;
00808 boolT waserror= False, nextseen= False, visibleseen= False;
00809
00810 trace1((qh ferr, "qh_checkpolygon: check all facets from f%d\n", facetlist->id));
00811 if (facetlist != qh facet_list || qh ONLYgood)
00812 nextseen= True;
00813 FORALLfacet_(facetlist) {
00814 if (facet == qh visible_list)
00815 visibleseen= True;
00816 if (!facet->visible) {
00817 if (!nextseen) {
00818 if (facet == qh facet_next)
00819 nextseen= True;
00820 else if (qh_setsize (facet->outsideset)
00821 && !(qh NARROWhull && qh CHECKfrequently)) {
00822 fprintf (qh ferr, "qhull internal error (qh_checkpolygon): f%d has outside set before qh facet_next\n",
00823 facet->id);
00824 qh_errexit (qh_ERRqhull, facet, NULL);
00825 }
00826 }
00827 numfacets++;
00828 qh_checkfacet(facet, False, &waserror);
00829 }
00830 }
00831 if (qh visible_list && !visibleseen && facetlist == qh facet_list) {
00832 fprintf (qh ferr, "qhull internal error (qh_checkpolygon): visible list f%d no longer on facet list\n", qh visible_list->id);
00833 qh_printlists();
00834 qh_errexit (qh_ERRqhull, qh visible_list, NULL);
00835 }
00836 if (facetlist == qh facet_list)
00837 vertexlist= qh vertex_list;
00838 else if (facetlist == qh newfacet_list)
00839 vertexlist= qh newvertex_list;
00840 else
00841 vertexlist= NULL;
00842 FORALLvertex_(vertexlist) {
00843 vertex->seen= False;
00844 vertex->visitid= 0;
00845 }
00846 FORALLfacet_(facetlist) {
00847 if (facet->visible)
00848 continue;
00849 if (facet->simplicial)
00850 numridges += qh hull_dim;
00851 else
00852 numridges += qh_setsize (facet->ridges);
00853 FOREACHvertex_(facet->vertices) {
00854 vertex->visitid++;
00855 if (!vertex->seen) {
00856 vertex->seen= True;
00857 numvertices++;
00858 if (qh_pointid (vertex->point) == -1) {
00859 fprintf (qh ferr, "qhull internal error (qh_checkpolygon): unknown point %p for vertex v%d first_point %p\n",
00860 vertex->point, vertex->id, qh first_point);
00861 waserror= True;
00862 }
00863 }
00864 }
00865 }
00866 qh vertex_visit += numfacets;
00867 if (facetlist == qh facet_list) {
00868 if (numfacets != qh num_facets - qh num_visible) {
00869 fprintf(qh ferr, "qhull internal error (qh_checkpolygon): actual number of facets is %d, cumulative facet count is %d\n",
00870 numfacets, qh num_facets- qh num_visible);
00871 waserror= True;
00872 }
00873 qh vertex_visit++;
00874 if (qh VERTEXneighbors) {
00875 FORALLvertices {
00876 qh_setcheck (vertex->neighbors, "neighbors for v", vertex->id);
00877 if (vertex->deleted)
00878 continue;
00879 totvneighbors += qh_setsize (vertex->neighbors);
00880 }
00881 FORALLfacet_(facetlist)
00882 totvertices += qh_setsize (facet->vertices);
00883 if (totvneighbors != totvertices) {
00884 fprintf(qh ferr, "qhull internal error (qh_checkpolygon): vertex neighbors inconsistent. Totvneighbors %d, totvertices %d\n",
00885 totvneighbors, totvertices);
00886 waserror= True;
00887 }
00888 }
00889 if (numvertices != qh num_vertices - qh_setsize(qh del_vertices)) {
00890 fprintf(qh ferr, "qhull internal error (qh_checkpolygon): actual number of vertices is %d, cumulative vertex count is %d\n",
00891 numvertices, qh num_vertices - qh_setsize(qh del_vertices));
00892 waserror= True;
00893 }
00894 if (qh hull_dim == 2 && numvertices != numfacets) {
00895 fprintf (qh ferr, "qhull internal error (qh_checkpolygon): #vertices %d != #facets %d\n",
00896 numvertices, numfacets);
00897 waserror= True;
00898 }
00899 if (qh hull_dim == 3 && numvertices + numfacets - numridges/2 != 2) {
00900 fprintf (qh ferr, "qhull internal error (qh_checkpolygon): #vertices %d + #facets %d - #edges %d != 2\n",
00901 numvertices, numfacets, numridges/2);
00902 waserror= True;
00903 }
00904 }
00905 if (waserror)
00906 qh_errexit(qh_ERRqhull, NULL, NULL);
00907 }
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920 void qh_checkvertex (vertexT *vertex) {
00921 boolT waserror= False;
00922 facetT *neighbor, **neighborp, *errfacet=NULL;
00923
00924 if (qh_pointid (vertex->point) == -1) {
00925 fprintf (qh ferr, "qhull internal error (qh_checkvertex): unknown point id %p\n", vertex->point);
00926 waserror= True;
00927 }
00928 if (vertex->id >= qh vertex_id) {
00929 fprintf (qh ferr, "qhull internal error (qh_checkvertex): unknown vertex id %d\n", vertex->id);
00930 waserror= True;
00931 }
00932 if (!waserror && !vertex->deleted) {
00933 if (qh_setsize (vertex->neighbors)) {
00934 FOREACHneighbor_(vertex) {
00935 if (!qh_setin (neighbor->vertices, vertex)) {
00936 fprintf (qh ferr, "qhull internal error (qh_checkvertex): neighbor f%d does not contain v%d\n", neighbor->id, vertex->id);
00937 errfacet= neighbor;
00938 waserror= True;
00939 }
00940 }
00941 }
00942 }
00943 if (waserror) {
00944 qh_errprint ("ERRONEOUS", NULL, NULL, NULL, vertex);
00945 qh_errexit (qh_ERRqhull, errfacet, NULL);
00946 }
00947 }
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959 void qh_clearcenters (qh_CENTER type) {
00960 facetT *facet;
00961
00962 if (qh CENTERtype != type) {
00963 FORALLfacets {
00964 if (qh CENTERtype == qh_ASvoronoi){
00965 if (facet->center) {
00966 qh_memfree (facet->center, qh center_size);
00967 facet->center= NULL;
00968 }
00969 }else {
00970 if (facet->center) {
00971 qh_memfree (facet->center, qh normal_size);
00972 facet->center= NULL;
00973 }
00974 }
00975 }
00976 qh CENTERtype= type;
00977 }
00978 trace2((qh ferr, "qh_clearcenters: switched to center type %d\n", type));
00979 }
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999 void qh_createsimplex(setT *vertices) {
01000 facetT *facet= NULL, *newfacet;
01001 boolT toporient= True;
01002 int vertex_i, vertex_n, nth;
01003 setT *newfacets= qh_settemp (qh hull_dim+1);
01004 vertexT *vertex;
01005
01006 qh facet_list= qh newfacet_list= qh facet_tail= qh_newfacet();
01007 qh num_facets= qh num_vertices= 0;
01008 qh vertex_list= qh newvertex_list= qh vertex_tail= qh_newvertex(NULL);
01009 FOREACHvertex_i_(vertices) {
01010 newfacet= qh_newfacet();
01011 newfacet->vertices= qh_setnew_delnthsorted (vertices, vertex_n,
01012 vertex_i, 0);
01013 newfacet->toporient= toporient;
01014 qh_appendfacet(newfacet);
01015 newfacet->newfacet= True;
01016 qh_appendvertex (vertex);
01017 qh_setappend (&newfacets, newfacet);
01018 toporient ^= True;
01019 }
01020 FORALLnew_facets {
01021 nth= 0;
01022 FORALLfacet_(qh newfacet_list) {
01023 if (facet != newfacet)
01024 SETelem_(newfacet->neighbors, nth++)= facet;
01025 }
01026 qh_settruncate (newfacet->neighbors, qh hull_dim);
01027 }
01028 qh_settempfree (&newfacets);
01029 trace1((qh ferr, "qh_createsimplex: created simplex\n"));
01030 }
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043 void qh_delridge(ridgeT *ridge) {
01044 void **freelistp;
01045
01046 qh_setdel(ridge->top->ridges, ridge);
01047 qh_setdel(ridge->bottom->ridges, ridge);
01048 qh_setfree(&(ridge->vertices));
01049 qh_memfree_(ridge, sizeof(ridgeT), freelistp);
01050 }
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063 void qh_delvertex (vertexT *vertex) {
01064
01065 if (vertex == qh tracevertex)
01066 qh tracevertex= NULL;
01067 qh_removevertex (vertex);
01068 qh_setfree (&vertex->neighbors);
01069 qh_memfree(vertex, sizeof(vertexT));
01070 }
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086 setT *qh_facet3vertex (facetT *facet) {
01087 ridgeT *ridge, *firstridge;
01088 vertexT *vertex;
01089 int cntvertices, cntprojected=0;
01090 setT *vertices;
01091
01092 cntvertices= qh_setsize(facet->vertices);
01093 vertices= qh_settemp (cntvertices);
01094 if (facet->simplicial) {
01095 if (cntvertices != 3) {
01096 fprintf (qh ferr, "qhull internal error (qh_facet3vertex): only %d vertices for simplicial facet f%d\n",
01097 cntvertices, facet->id);
01098 qh_errexit(qh_ERRqhull, facet, NULL);
01099 }
01100 qh_setappend (&vertices, SETfirst_(facet->vertices));
01101 if (facet->toporient ^ qh_ORIENTclock)
01102 qh_setappend (&vertices, SETsecond_(facet->vertices));
01103 else
01104 qh_setaddnth (&vertices, 0, SETsecond_(facet->vertices));
01105 qh_setappend (&vertices, SETelem_(facet->vertices, 2));
01106 }else {
01107 ridge= firstridge= SETfirstt_(facet->ridges, ridgeT);
01108 while ((ridge= qh_nextridge3d (ridge, facet, &vertex))) {
01109 qh_setappend (&vertices, vertex);
01110 if (++cntprojected > cntvertices || ridge == firstridge)
01111 break;
01112 }
01113 if (!ridge || cntprojected != cntvertices) {
01114 fprintf (qh ferr, "qhull internal error (qh_facet3vertex): ridges for facet %d don't match up. got at least %d\n",
01115 facet->id, cntprojected);
01116 qh_errexit(qh_ERRqhull, facet, ridge);
01117 }
01118 }
01119 return vertices;
01120 }
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151 facetT *qh_findbestfacet (pointT *point, boolT bestoutside,
01152 realT *bestdist, boolT *isoutside) {
01153 facetT *bestfacet= NULL;
01154 int numpart, totpart= 0;
01155
01156 bestfacet= qh_findbest (point, qh facet_list,
01157 bestoutside, False, bestoutside,
01158 bestdist, isoutside, &totpart);
01159 if (!bestfacet) {
01160 fprintf (qh ferr, "qh_findbestfacet: all facets are flipped or upper Delaunay\n");
01161 qh_errexit (qh_ERRqhull, NULL, NULL);
01162 }
01163 if (*bestdist < -qh DISTround) {
01164 bestfacet= qh_findfacet_all (point, bestdist, isoutside, &numpart);
01165 totpart += numpart;
01166 if ((isoutside && bestoutside)
01167 || (!isoutside && bestfacet->upperdelaunay)) {
01168 bestfacet= qh_findbest (point, bestfacet,
01169 bestoutside, False, bestoutside,
01170 bestdist, isoutside, &totpart);
01171 totpart += numpart;
01172 }
01173 }
01174 trace3((qh ferr, "qh_findbestfacet: f%d dist %2.2g isoutside %d totpart %d\n",
01175 bestfacet->id, *bestdist, *isoutside, totpart));
01176 return bestfacet;
01177 }
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197 facetT *qh_findfacet_all (pointT *point, realT *bestdist, boolT *isoutside,
01198 int *numpart) {
01199 facetT *bestfacet= NULL, *facet;
01200 realT dist;
01201 int totpart= 0;
01202
01203 *bestdist= REALmin;
01204 *isoutside= False;
01205 FORALLfacets {
01206 if (facet->flipped || !facet->normal)
01207 continue;
01208 totpart++;
01209 qh_distplane (point, facet, &dist);
01210 if (dist > *bestdist) {
01211 *bestdist= dist;
01212 bestfacet= facet;
01213 if (dist > qh MINoutside) {
01214 *isoutside= True;
01215 break;
01216 }
01217 }
01218 }
01219 *numpart= totpart;
01220 trace3((qh ferr, "qh_findfacet_all: f%d dist %2.2g isoutside %d totpart %d\n",
01221 getid_(bestfacet), *bestdist, *isoutside, totpart));
01222 return bestfacet;
01223 }
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257 int qh_findgood (facetT *facetlist, int goodhorizon) {
01258 facetT *facet, *bestfacet= NULL;
01259 realT angle, bestangle= REALmax, dist;
01260 int numgood=0;
01261
01262 FORALLfacet_(facetlist) {
01263 if (facet->good)
01264 numgood++;
01265 }
01266 if (qh GOODvertex>0 && !qh MERGING) {
01267 FORALLfacet_(facetlist) {
01268 if (!qh_isvertex (qh GOODvertexp, facet->vertices)) {
01269 facet->good= False;
01270 numgood--;
01271 }
01272 }
01273 }
01274 if (qh GOODpoint && numgood) {
01275 FORALLfacet_(facetlist) {
01276 if (facet->good && facet->normal) {
01277 zinc_(Zdistgood);
01278 qh_distplane (qh GOODpointp, facet, &dist);
01279 if ((qh GOODpoint > 0) ^ (dist > 0.0)) {
01280 facet->good= False;
01281 numgood--;
01282 }
01283 }
01284 }
01285 }
01286 if (qh GOODthreshold && (numgood || goodhorizon || qh GOODclosest)) {
01287 FORALLfacet_(facetlist) {
01288 if (facet->good && facet->normal) {
01289 if (!qh_inthresholds (facet->normal, &angle)) {
01290 facet->good= False;
01291 numgood--;
01292 if (angle < bestangle) {
01293 bestangle= angle;
01294 bestfacet= facet;
01295 }
01296 }
01297 }
01298 }
01299 if (!numgood && (!goodhorizon || qh GOODclosest)) {
01300 if (qh GOODclosest) {
01301 if (qh GOODclosest->visible)
01302 qh GOODclosest= NULL;
01303 else {
01304 qh_inthresholds (qh GOODclosest->normal, &angle);
01305 if (angle < bestangle)
01306 bestfacet= qh GOODclosest;
01307 }
01308 }
01309 if (bestfacet && bestfacet != qh GOODclosest) {
01310 if (qh GOODclosest)
01311 qh GOODclosest->good= False;
01312 qh GOODclosest= bestfacet;
01313 bestfacet->good= True;
01314 numgood++;
01315 trace2((qh ferr, "qh_findgood: f%d is closest (%2.2g) to thresholds\n",
01316 bestfacet->id, bestangle));
01317 return numgood;
01318 }
01319 }else if (qh GOODclosest) {
01320 qh GOODclosest->good= False;
01321 qh GOODclosest= NULL;
01322 }
01323 }
01324 zadd_(Zgoodfacet, numgood);
01325 trace2((qh ferr, "qh_findgood: found %d good facets with %d good horizon\n",
01326 numgood, goodhorizon));
01327 if (!numgood && qh GOODvertex>0 && !qh MERGING)
01328 return goodhorizon;
01329 return numgood;
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
01355
01356
01357 void qh_findgood_all (facetT *facetlist) {
01358 facetT *facet, *bestfacet=NULL;
01359 realT angle, bestangle= REALmax;
01360 int numgood=0, startgood;
01361
01362 if (!qh GOODvertex && !qh GOODthreshold && !qh GOODpoint
01363 && !qh SPLITthresholds)
01364 return;
01365 if (!qh ONLYgood)
01366 qh_findgood (qh facet_list, 0);
01367 FORALLfacet_(facetlist) {
01368 if (facet->good)
01369 numgood++;
01370 }
01371 if (qh GOODvertex <0 || (qh GOODvertex > 0 && qh MERGING)) {
01372 FORALLfacet_(facetlist) {
01373 if (facet->good && ((qh GOODvertex > 0) ^ !!qh_isvertex (qh GOODvertexp, facet->vertices))) {
01374 if (!--numgood) {
01375 if (qh ONLYgood) {
01376 fprintf (qh ferr, "qhull warning: good vertex p%d does not match last good facet f%d. Ignored.\n",
01377 qh_pointid(qh GOODvertexp), facet->id);
01378 return;
01379 }else if (qh GOODvertex > 0)
01380 fprintf (qh ferr, "qhull warning: point p%d is not a vertex ('QV%d').\n",
01381 qh GOODvertex-1, qh GOODvertex-1);
01382 else
01383 fprintf (qh ferr, "qhull warning: point p%d is a vertex for every facet ('QV-%d').\n",
01384 -qh GOODvertex - 1, -qh GOODvertex - 1);
01385 }
01386 facet->good= False;
01387 }
01388 }
01389 }
01390 startgood= numgood;
01391 if (qh SPLITthresholds) {
01392 FORALLfacet_(facetlist) {
01393 if (facet->good) {
01394 if (!qh_inthresholds (facet->normal, &angle)) {
01395 facet->good= False;
01396 numgood--;
01397 if (angle < bestangle) {
01398 bestangle= angle;
01399 bestfacet= facet;
01400 }
01401 }
01402 }
01403 }
01404 if (!numgood && bestfacet) {
01405 bestfacet->good= True;
01406 numgood++;
01407 trace0((qh ferr, "qh_findgood_all: f%d is closest (%2.2g) to thresholds\n",
01408 bestfacet->id, bestangle));
01409 return;
01410 }
01411 }
01412 qh num_good= numgood;
01413 trace0((qh ferr, "qh_findgood_all: %d good facets remain out of %d facets\n",
01414 numgood, startgood));
01415 }
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427 void qh_furthestnext (void ) {
01428 facetT *facet, *bestfacet= NULL;
01429 realT dist, bestdist= -REALmax;
01430
01431 FORALLfacets {
01432 if (facet->outsideset) {
01433 #if qh_COMPUTEfurthest
01434 pointT *furthest;
01435 furthest= (pointT*)qh_setlast (facet->outsideset);
01436 zinc_(Zcomputefurthest);
01437 qh_distplane (furthest, facet, &dist);
01438 #else
01439 dist= facet->furthestdist;
01440 #endif
01441 if (dist > bestdist) {
01442 bestfacet= facet;
01443 bestdist= dist;
01444 }
01445 }
01446 }
01447 if (bestfacet) {
01448 qh_removefacet (bestfacet);
01449 qh_prependfacet (bestfacet, &qh facet_next);
01450 trace1((qh ferr, "qh_furthestnext: made f%d next facet (dist %.2g)\n",
01451 bestfacet->id, bestdist));
01452 }
01453 }
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470 void qh_furthestout (facetT *facet) {
01471 pointT *point, **pointp, *bestpoint= NULL;
01472 realT dist, bestdist= -REALmax;
01473
01474 FOREACHpoint_(facet->outsideset) {
01475 qh_distplane (point, facet, &dist);
01476 zinc_(Zcomputefurthest);
01477 if (dist > bestdist) {
01478 bestpoint= point;
01479 bestdist= dist;
01480 }
01481 }
01482 if (bestpoint) {
01483 qh_setdel (facet->outsideset, point);
01484 qh_setappend (&facet->outsideset, point);
01485 #if !qh_COMPUTEfurthest
01486 facet->furthestdist= bestdist;
01487 #endif
01488 }
01489 facet->notfurthest= False;
01490 trace3((qh ferr, "qh_furthestout: p%d is furthest outside point of f%d\n",
01491 qh_pointid (point), facet->id));
01492 }
01493
01494
01495
01496
01497
01498
01499
01500
01501 void qh_infiniteloop (facetT *facet) {
01502
01503 fprintf (qh ferr, "qhull internal error (qh_infiniteloop): potential infinite loop detected\n");
01504 qh_errexit (qh_ERRqhull, facet, NULL);
01505 }
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532 void qh_initbuild( void) {
01533 setT *maxpoints, *vertices;
01534 facetT *facet;
01535 int i, numpart;
01536 realT dist;
01537 boolT isoutside;
01538
01539 qh furthest_id= -1;
01540 qh lastreport= 0;
01541 qh facet_id= qh vertex_id= qh ridge_id= 0;
01542 qh visit_id= qh vertex_visit= 0;
01543 qh maxoutdone= False;
01544
01545 if (qh GOODpoint > 0)
01546 qh GOODpointp= qh_point (qh GOODpoint-1);
01547 else if (qh GOODpoint < 0)
01548 qh GOODpointp= qh_point (-qh GOODpoint-1);
01549 if (qh GOODvertex > 0)
01550 qh GOODvertexp= qh_point (qh GOODvertex-1);
01551 else if (qh GOODvertex < 0)
01552 qh GOODvertexp= qh_point (-qh GOODvertex-1);
01553 if ((qh GOODpoint
01554 && (qh GOODpointp < qh first_point
01555 || qh GOODpointp > qh_point (qh num_points-1)))
01556 || (qh GOODvertex
01557 && (qh GOODvertexp < qh first_point
01558 || qh GOODvertexp > qh_point (qh num_points-1)))) {
01559 fprintf (qh ferr, "qhull input error: either QGn or QVn point is > p%d\n",
01560 qh num_points-1);
01561 qh_errexit (qh_ERRinput, NULL, NULL);
01562 }
01563 maxpoints= qh_maxmin(qh first_point, qh num_points, qh hull_dim);
01564 if (qh SCALElast)
01565 qh_scalelast (qh first_point, qh num_points, qh hull_dim,
01566 qh MINlastcoord, qh MAXlastcoord, qh MAXwidth);
01567 qh_detroundoff();
01568 if (qh DELAUNAY && qh upper_threshold[qh hull_dim-1] > REALmax/2
01569 && qh lower_threshold[qh hull_dim-1] < -REALmax/2) {
01570 for (i= qh_PRINTEND; i--; ) {
01571 if (qh PRINTout[i] == qh_PRINTgeom && qh DROPdim < 0
01572 && !qh GOODthreshold && !qh SPLITthresholds)
01573 break;
01574 }
01575 if (i < 0) {
01576 if (qh UPPERdelaunay) {
01577 qh lower_threshold[qh hull_dim-1]= qh ANGLEround * qh_ZEROdelaunay;
01578 qh GOODthreshold= True;
01579 }else {
01580 qh upper_threshold[qh hull_dim-1]= -qh ANGLEround * qh_ZEROdelaunay;
01581 if (!qh GOODthreshold)
01582 qh SPLITthresholds= True;
01583
01584 }
01585 }
01586 }
01587 vertices= qh_initialvertices(qh hull_dim, maxpoints, qh first_point, qh num_points);
01588 qh_initialhull (vertices);
01589 qh_partitionall (vertices, qh first_point, qh num_points);
01590 if (qh PRINToptions1st || qh TRACElevel || qh IStracing) {
01591 if (qh TRACElevel || qh IStracing)
01592 fprintf (qh ferr, "\nTrace level %d for %s | %s\n",
01593 qh IStracing ? qh IStracing : qh TRACElevel, qh rbox_command, qh qhull_command);
01594 fprintf (qh ferr, "Options selected for Qhull %s:\n%s\n", qh_version, qh qhull_options);
01595 }
01596 qh_resetlists (False );
01597 qh facet_next= qh facet_list;
01598 qh_furthestnext ();
01599 if (qh PREmerge) {
01600 qh cos_max= qh premerge_cos;
01601 qh centrum_radius= qh premerge_centrum;
01602 }
01603 if (qh ONLYgood) {
01604 if (qh GOODvertex > 0 && qh MERGING) {
01605 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");
01606 qh_errexit (qh_ERRinput, NULL, NULL);
01607 }
01608 if (!(qh GOODthreshold || qh GOODpoint
01609 || (!qh MERGEexact && !qh PREmerge && qh GOODvertexp))) {
01610 fprintf (qh ferr, "qhull input error: 'Qg' (ONLYgood) needs a good threshold ('Pd0D0'), a\n\
01611 good point (QGn or QG-n), or a good vertex with 'QJ' or 'Q0' (QVn).\n");
01612 qh_errexit (qh_ERRinput, NULL, NULL);
01613 }
01614 if (qh GOODvertex > 0 && !qh MERGING
01615 && !qh_isvertex (qh GOODvertexp, vertices)) {
01616 facet= qh_findbestnew (qh GOODvertexp, qh facet_list,
01617 &dist, &isoutside, &numpart);
01618 zadd_(Zdistgood, numpart);
01619 if (!isoutside) {
01620 fprintf (qh ferr, "qhull input error: point for QV%d is inside initial simplex. It can not be made a vertex.\n",
01621 qh_pointid(qh GOODvertexp));
01622 qh_errexit (qh_ERRinput, NULL, NULL);
01623 }
01624 if (!qh_addpoint (qh GOODvertexp, facet, False)) {
01625 qh_settempfree(&vertices);
01626 qh_settempfree(&maxpoints);
01627 return;
01628 }
01629 }
01630 qh_findgood (qh facet_list, 0);
01631 }
01632 qh_settempfree(&vertices);
01633 qh_settempfree(&maxpoints);
01634 trace1((qh ferr, "qh_initbuild: initial hull created and points partitioned\n"));
01635 }
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651 void qh_initialhull(setT *vertices) {
01652 facetT *facet, *firstfacet, *neighbor, **neighborp;
01653 realT dist, angle, minangle= REALmax;
01654 #ifndef qh_NOtrace
01655 int k;
01656 #endif
01657
01658 qh_createsimplex(vertices);
01659 qh_resetlists (False);
01660 qh facet_next= qh facet_list;
01661 qh interior_point= qh_getcenter(vertices);
01662 firstfacet= qh facet_list;
01663 qh_setfacetplane(firstfacet);
01664 zinc_(Znumvisibility);
01665 qh_distplane(qh interior_point, firstfacet, &dist);
01666 if (dist > 0) {
01667 FORALLfacets
01668 facet->toporient ^= True;
01669 }
01670 FORALLfacets
01671 qh_setfacetplane(facet);
01672 FORALLfacets {
01673 if (!qh_checkflipped (facet, NULL, qh_ALL)) {
01674 trace1((qh ferr, "qh_initialhull: initial orientation incorrect. Correct all facets\n"));
01675 facet->flipped= False;
01676 FORALLfacets {
01677 facet->toporient ^= True;
01678 qh_orientoutside (facet);
01679 }
01680 break;
01681 }
01682 }
01683 FORALLfacets {
01684 if (!qh_checkflipped (facet, NULL, !qh_ALL)) {
01685 qh_precision ("initial facet is coplanar with interior point");
01686 fprintf (qh ferr, "qhull precision error: initial facet %d is coplanar with the interior point\n",
01687 facet->id);
01688 qh_errexit (qh_ERRsingular, facet, NULL);
01689 }
01690 FOREACHneighbor_(facet) {
01691 angle= qh_getangle (facet->normal, neighbor->normal);
01692 minimize_( minangle, angle);
01693 }
01694 }
01695 if (minangle < qh_MAXnarrow) {
01696 realT diff= 1.0 + minangle;
01697
01698 qh NARROWhull= True;
01699 qh_option ("_narrow-hull", NULL, &diff);
01700 if (minangle < qh_WARNnarrow && !qh RERUN && qh PRINTprecision)
01701 fprintf (qh ferr, "qhull precision warning: \n\
01702 The initial hull is narrow (the cosine of the minimum angle is %.9g).\n\
01703 A coplanar point may lead to a wide facet. Options 'Qs' (search for best\n\
01704 initial hull), 'QbB' (scale to unit box), or 'Qbb' (scale last coordinate)\n\
01705 may remove this warning. Use 'Pp' to ignore this warning.\n\
01706 See 'Limitations' in qh-impre.htm.\n",
01707 -minangle);
01708 }
01709 zzval_(Zprocessed)= qh hull_dim+1;
01710 qh_checkpolygon (qh facet_list);
01711 qh_checkconvex(qh facet_list, qh_DATAfault);
01712 #ifndef qh_NOtrace
01713 if (qh IStracing >= 1) {
01714 fprintf(qh ferr, "qh_initialhull: simplex constructed, interior point:");
01715 for (k=0; k < qh hull_dim; k++)
01716 fprintf (qh ferr, " %6.4g", qh interior_point[k]);
01717 fprintf (qh ferr, "\n");
01718 }
01719 #endif
01720 }
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740 setT *qh_initialvertices(int dim, setT *maxpoints, pointT *points, int numpoints) {
01741 pointT *point, **pointp;
01742 setT *vertices, *simplex, *tested;
01743 realT randr;
01744 int index, point_i, point_n, k;
01745 boolT nearzero= False;
01746
01747 vertices= qh_settemp (dim + 1);
01748 simplex= qh_settemp (dim+1);
01749 if (qh ALLpoints)
01750 qh_maxsimplex (dim, NULL, points, numpoints, &simplex);
01751 else if (qh RANDOMoutside) {
01752 while (qh_setsize (simplex) != dim+1) {
01753 randr= qh_RANDOMint;
01754 randr= randr/(qh_RANDOMmax+1);
01755 index= (int)floor(qh num_points * randr);
01756 while (qh_setin (simplex, qh_point (index))) {
01757 index++;
01758 index= index < qh num_points ? index : 0;
01759 }
01760 qh_setappend (&simplex, qh_point (index));
01761 }
01762 }else if (qh hull_dim >= qh_INITIALmax) {
01763 tested= qh_settemp (dim+1);
01764 qh_setappend (&simplex, SETfirst_(maxpoints));
01765 qh_setappend (&simplex, SETsecond_(maxpoints));
01766 qh_maxsimplex (fmin_(qh_INITIALsearch, dim), maxpoints, points, numpoints, &simplex);
01767 k= qh_setsize (simplex);
01768 FOREACHpoint_i_(maxpoints) {
01769 if (point_i & 0x1) {
01770 if (!qh_setin (simplex, point) && !qh_setin (tested, point)){
01771 qh_detsimplex(point, simplex, k, &nearzero);
01772 if (nearzero)
01773 qh_setappend (&tested, point);
01774 else {
01775 qh_setappend (&simplex, point);
01776 if (++k == dim)
01777 break;
01778 }
01779 }
01780 }
01781 }
01782 while (k != dim && (point= (pointT*)qh_setdellast (maxpoints))) {
01783 if (!qh_setin (simplex, point) && !qh_setin (tested, point)){
01784 qh_detsimplex (point, simplex, k, &nearzero);
01785 if (nearzero)
01786 qh_setappend (&tested, point);
01787 else {
01788 qh_setappend (&simplex, point);
01789 k++;
01790 }
01791 }
01792 }
01793 index= 0;
01794 while (k != dim && (point= qh_point (index++))) {
01795 if (!qh_setin (simplex, point) && !qh_setin (tested, point)){
01796 qh_detsimplex (point, simplex, k, &nearzero);
01797 if (!nearzero){
01798 qh_setappend (&simplex, point);
01799 k++;
01800 }
01801 }
01802 }
01803 qh_settempfree (&tested);
01804 qh_maxsimplex (dim, maxpoints, points, numpoints, &simplex);
01805 }else
01806 qh_maxsimplex (dim, maxpoints, points, numpoints, &simplex);
01807 FOREACHpoint_(simplex)
01808 qh_setaddnth (&vertices, 0, qh_newvertex(point));
01809 qh_settempfree (&simplex);
01810 return vertices;
01811 }
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823 vertexT *qh_isvertex (pointT *point, setT *vertices) {
01824 vertexT *vertex, **vertexp;
01825
01826 FOREACHvertex_(vertices) {
01827 if (vertex->point == point)
01828 return vertex;
01829 }
01830 return NULL;
01831 }
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859 vertexT *qh_makenewfacets (pointT *point ) {
01860 facetT *visible, *newfacet= NULL, *newfacet2= NULL, *neighbor, **neighborp;
01861 vertexT *apex;
01862 int numnew=0;
01863
01864 qh newfacet_list= qh facet_tail;
01865 qh newvertex_list= qh vertex_tail;
01866 apex= qh_newvertex(point);
01867 qh_appendvertex (apex);
01868 qh visit_id++;
01869 if (!qh ONLYgood)
01870 qh NEWfacets= True;
01871 FORALLvisible_facets {
01872 FOREACHneighbor_(visible)
01873 neighbor->seen= False;
01874 if (visible->ridges) {
01875 visible->visitid= qh visit_id;
01876 newfacet2= qh_makenew_nonsimplicial (visible, apex, &numnew);
01877 }
01878 if (visible->simplicial)
01879 newfacet= qh_makenew_simplicial (visible, apex, &numnew);
01880 if (!qh ONLYgood) {
01881 if (newfacet2)
01882 newfacet= newfacet2;
01883 if (newfacet)
01884 visible->f.replace= newfacet;
01885 else
01886 zinc_(Zinsidevisible);
01887 SETfirst_(visible->neighbors)= NULL;
01888 }
01889 }
01890 trace1((qh ferr, "qh_makenewfacets: created %d new facets from point p%d to horizon\n",
01891 numnew, qh_pointid(point)));
01892 if (qh IStracing >= 4)
01893 qh_printfacetlist (qh newfacet_list, NULL, qh_ALL);
01894 return apex;
01895 }
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922 #ifndef qh_NOmerge
01923 void qh_matchduplicates (facetT *atfacet, int atskip, int hashsize, int *hashcount) {
01924 boolT same, ismatch;
01925 int hash, scan;
01926 facetT *facet, *newfacet, *maxmatch= NULL, *maxmatch2= NULL, *nextfacet;
01927 int skip, newskip, nextskip= 0, maxskip= 0, maxskip2= 0, makematch;
01928 realT maxdist= -REALmax, mindist, dist2, low, high;
01929
01930 hash= (int)qh_gethash (hashsize, atfacet->vertices, qh hull_dim, 1,
01931 SETelem_(atfacet->vertices, atskip));
01932 trace2((qh ferr, "qh_matchduplicates: find duplicate matches for f%d skip %d hash %d hashcount %d\n",
01933 atfacet->id, atskip, hash, *hashcount));
01934 for (makematch= 0; makematch < 2; makematch++) {
01935 qh visit_id++;
01936 for (newfacet= atfacet, newskip= atskip; newfacet; newfacet= nextfacet, newskip= nextskip) {
01937 zinc_(Zhashlookup);
01938 nextfacet= NULL;
01939 newfacet->visitid= qh visit_id;
01940 for (scan= hash; (facet= SETelemt_(qh hash_table, scan, facetT));
01941 scan= (++scan >= hashsize ? 0 : scan)) {
01942 if (!facet->dupridge || facet->visitid == qh visit_id)
01943 continue;
01944 zinc_(Zhashtests);
01945 if (qh_matchvertices (1, newfacet->vertices, newskip, facet->vertices, &skip, &same)) {
01946 ismatch= (same == (newfacet->toporient ^ facet->toporient));
01947 if (SETelemt_(facet->neighbors, skip, facetT) != qh_DUPLICATEridge) {
01948 if (!makematch) {
01949 fprintf (qh ferr, "qhull internal error (qh_matchduplicates): missing dupridge at f%d skip %d for new f%d skip %d hash %d\n",
01950 facet->id, skip, newfacet->id, newskip, hash);
01951 qh_errexit2 (qh_ERRqhull, facet, newfacet);
01952 }
01953 }else if (ismatch && makematch) {
01954 if (SETelemt_(newfacet->neighbors, newskip, facetT) == qh_DUPLICATEridge) {
01955 SETelem_(facet->neighbors, skip)= newfacet;
01956 SETelem_(newfacet->neighbors, newskip)= qh_MERGEridge;
01957 *hashcount -= 2;
01958 trace4((qh ferr, "qh_matchduplicates: duplicate f%d skip %d matched with new f%d skip %d merge\n",
01959 facet->id, skip, newfacet->id, newskip));
01960 }
01961 }else if (ismatch) {
01962 mindist= qh_getdistance (facet, newfacet, &low, &high);
01963 dist2= qh_getdistance (newfacet, facet, &low, &high);
01964 minimize_(mindist, dist2);
01965 if (mindist > maxdist) {
01966 maxdist= mindist;
01967 maxmatch= facet;
01968 maxskip= skip;
01969 maxmatch2= newfacet;
01970 maxskip2= newskip;
01971 }
01972 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",
01973 facet->id, skip, newfacet->id, newskip, mindist,
01974 maxmatch->id, maxmatch2->id));
01975 }else {
01976 nextfacet= facet;
01977 nextskip= skip;
01978 }
01979 }
01980 if (makematch && !facet
01981 && SETelemt_(facet->neighbors, skip, facetT) == qh_DUPLICATEridge) {
01982 fprintf (qh ferr, "qhull internal error (qh_matchduplicates): no MERGEridge match for duplicate f%d skip %d at hash %d\n",
01983 newfacet->id, newskip, hash);
01984 qh_errexit (qh_ERRqhull, newfacet, NULL);
01985 }
01986 }
01987 }
01988 if (!makematch) {
01989 if (!maxmatch) {
01990 fprintf (qh ferr, "qhull internal error (qh_matchduplicates): no maximum match at duplicate f%d skip %d at hash %d\n",
01991 atfacet->id, atskip, hash);
01992 qh_errexit (qh_ERRqhull, atfacet, NULL);
01993 }
01994 SETelem_(maxmatch->neighbors, maxskip)= maxmatch2;
01995 SETelem_(maxmatch2->neighbors, maxskip2)= maxmatch;
01996 *hashcount -= 2;
01997 zzinc_(Zmultiridge);
01998 trace0((qh ferr, "qh_matchduplicates: duplicate f%d skip %d matched with new f%d skip %d keep\n",
01999 maxmatch->id, maxskip, maxmatch2->id, maxskip2));
02000 qh_precision ("ridge with multiple neighbors");
02001 if (qh IStracing >= 4)
02002 qh_errprint ("DUPLICATED/MATCH", maxmatch, maxmatch2, NULL, NULL);
02003 }
02004 }
02005 }
02006
02007
02008
02009
02010
02011
02012
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030 void qh_nearcoplanar ( void ) {
02031 facetT *facet;
02032 pointT *point, **pointp;
02033 int numpart;
02034 realT dist, innerplane;
02035
02036 if (!qh KEEPcoplanar && !qh KEEPinside) {
02037 FORALLfacets {
02038 if (facet->coplanarset)
02039 qh_setfree( &facet->coplanarset);
02040 }
02041 }else if (!qh KEEPcoplanar || !qh KEEPinside) {
02042 qh_outerinner (NULL, NULL, &innerplane);
02043 if (qh JOGGLEmax < REALmax/2)
02044 innerplane -= qh JOGGLEmax * sqrt (qh hull_dim);
02045 numpart= 0;
02046 FORALLfacets {
02047 if (facet->coplanarset) {
02048 FOREACHpoint_(facet->coplanarset) {
02049 numpart++;
02050 qh_distplane (point, facet, &dist);
02051 if (dist < innerplane) {
02052 if (!qh KEEPinside)
02053 SETref_(point)= NULL;
02054 }else if (!qh KEEPcoplanar)
02055 SETref_(point)= NULL;
02056 }
02057 qh_setcompact (facet->coplanarset);
02058 }
02059 }
02060 zadd_(Zcheckpart, numpart);
02061 }
02062 }
02063
02064
02065
02066
02067
02068
02069
02070
02071
02072
02073
02074
02075
02076
02077 vertexT *qh_nearvertex (facetT *facet, pointT *point, realT *bestdistp) {
02078 realT bestdist= REALmax, dist;
02079 vertexT *bestvertex= NULL, *vertex, **vertexp;
02080 int dim= qh hull_dim;
02081
02082 if (qh DELAUNAY)
02083 dim--;
02084 FOREACHvertex_(facet->vertices) {
02085 dist= qh_pointdist (vertex->point, point, -dim);
02086 if (dist < bestdist) {
02087 bestdist= dist;
02088 bestvertex= vertex;
02089 }
02090 }
02091 *bestdistp= sqrt (bestdist);
02092 return bestvertex;
02093 }
02094
02095
02096
02097
02098
02099
02100
02101
02102
02103
02104
02105
02106 int qh_newhashtable(int newsize) {
02107 int size;
02108
02109 size= ((newsize+1)*qh_HASHfactor) | 0x1;
02110 while (True) {
02111 if ((size%3) && (size%5))
02112 break;
02113 size += 2;
02114
02115 }
02116 qh hash_table= qh_setnew (size);
02117 qh_setzero (qh hash_table, 0, size);
02118 return size;
02119 }
02120
02121
02122
02123
02124
02125
02126
02127 vertexT *qh_newvertex(pointT *point) {
02128 vertexT *vertex;
02129
02130 zinc_(Ztotvertices);
02131 vertex= (vertexT *)qh_memalloc(sizeof(vertexT));
02132 memset ((char *) vertex, 0, sizeof (vertexT));
02133 if (qh vertex_id == 0xFFFFFF) {
02134 fprintf(qh ferr, "qhull input error: more than %d vertices. ID field overflows and two vertices\n\
02135 may have the same identifier. Vertices not sorted correctly.\n", 0xFFFFFF);
02136 qh_errexit(qh_ERRinput, NULL, NULL);
02137 }
02138 if (qh vertex_id == qh tracevertex_id)
02139 qh tracevertex= vertex;
02140 vertex->id= qh vertex_id++;
02141 vertex->point= point;
02142 trace4((qh ferr, "qh_newvertex: vertex p%d (v%d) created\n", qh_pointid(vertex->point),
02143 vertex->id));
02144 return (vertex);
02145 }
02146
02147
02148
02149
02150
02151
02152
02153
02154
02155
02156
02157
02158
02159
02160
02161
02162 ridgeT *qh_nextridge3d (ridgeT *atridge, facetT *facet, vertexT **vertexp) {
02163 vertexT *atvertex, *vertex, *othervertex;
02164 ridgeT *ridge, **ridgep;
02165
02166 if ((atridge->top == facet) ^ qh_ORIENTclock)
02167 atvertex= SETsecondt_(atridge->vertices, vertexT);
02168 else
02169 atvertex= SETfirstt_(atridge->vertices, vertexT);
02170 FOREACHridge_(facet->ridges) {
02171 if (ridge == atridge)
02172 continue;
02173 if ((ridge->top == facet) ^ qh_ORIENTclock) {
02174 othervertex= SETsecondt_(ridge->vertices, vertexT);
02175 vertex= SETfirstt_(ridge->vertices, vertexT);
02176 }else {
02177 vertex= SETsecondt_(ridge->vertices, vertexT);
02178 othervertex= SETfirstt_(ridge->vertices, vertexT);
02179 }
02180 if (vertex == atvertex) {
02181 if (vertexp)
02182 *vertexp= othervertex;
02183 return ridge;
02184 }
02185 }
02186 return NULL;
02187 }
02188 #else
02189 void qh_matchduplicates (facetT *atfacet, int atskip, int hashsize, int *hashcount) {
02190 }
02191 ridgeT *qh_nextridge3d (ridgeT *atridge, facetT *facet, vertexT **vertexp) {
02192
02193 return NULL;
02194 }
02195 #endif
02196
02197
02198
02199
02200
02201
02202
02203
02204
02205
02206
02207
02208
02209
02210
02211 void qh_outcoplanar (void ) {
02212 pointT *point, **pointp;
02213 facetT *facet;
02214 realT dist;
02215
02216 trace1((qh ferr, "qh_outcoplanar: move outsideset to coplanarset for qh NARROWhull\n"));
02217 FORALLfacets {
02218 FOREACHpoint_(facet->outsideset) {
02219 qh num_outside--;
02220 if (qh KEEPcoplanar || qh KEEPnearinside) {
02221 qh_distplane (point, facet, &dist);
02222 zinc_(Zpartition);
02223 qh_partitioncoplanar (point, facet, &dist);
02224 }
02225 }
02226 qh_setfree (&facet->outsideset);
02227 }
02228 }
02229
02230
02231
02232
02233
02234
02235
02236
02237
02238
02239
02240 pointT *qh_point (int id) {
02241
02242 if (id < 0)
02243 return NULL;
02244 if (id < qh num_points)
02245 return qh first_point + id * qh hull_dim;
02246 id -= qh num_points;
02247 if (id < qh_setsize (qh other_points))
02248 return SETelemt_(qh other_points, id, pointT);
02249 return NULL;
02250 }
02251
02252
02253
02254
02255
02256
02257
02258
02259
02260
02261
02262
02263
02264 void qh_point_add (setT *set, pointT *point, void *elem) {
02265 int id, size;
02266
02267 SETreturnsize_(set, size);
02268 if ((id= qh_pointid(point)) < 0)
02269 fprintf (qh ferr, "qhull internal warning (point_add): unknown point %p id %d\n",
02270 point, id);
02271 else if (id >= size) {
02272 fprintf (qh ferr, "qhull internal errror (point_add): point p%d is out of bounds (%d)\n",
02273 id, size);
02274 qh_errexit (qh_ERRqhull, NULL, NULL);
02275 }else
02276 SETelem_(set, id)= elem;
02277 }
02278
02279
02280
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304 setT *qh_pointfacet (void ) {
02305 int numpoints= qh num_points + qh_setsize (qh other_points);
02306 setT *facets;
02307 facetT *facet;
02308 vertexT *vertex, **vertexp;
02309 pointT *point, **pointp;
02310
02311 facets= qh_settemp (numpoints);
02312 qh_setzero (facets, 0, numpoints);
02313 qh vertex_visit++;
02314 FORALLfacets {
02315 FOREACHvertex_(facet->vertices) {
02316 if (vertex->visitid != qh vertex_visit) {
02317 vertex->visitid= qh vertex_visit;
02318 qh_point_add (facets, vertex->point, facet);
02319 }
02320 }
02321 FOREACHpoint_(facet->coplanarset)
02322 qh_point_add (facets, point, facet);
02323 FOREACHpoint_(facet->outsideset)
02324 qh_point_add (facets, point, facet);
02325 }
02326 return facets;
02327 }
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338
02339
02340
02341 setT *qh_pointvertex (void ) {
02342 int numpoints= qh num_points + qh_setsize (qh other_points);
02343 setT *vertices;
02344 vertexT *vertex;
02345
02346 vertices= qh_settemp (numpoints);
02347 qh_setzero (vertices, 0, numpoints);
02348 FORALLvertices
02349 qh_point_add (vertices, vertex->point, vertex);
02350 return vertices;
02351 }
02352
02353
02354
02355
02356
02357
02358
02359
02360
02361
02362
02363
02364
02365
02366
02367
02368 void qh_prependfacet(facetT *facet, facetT **facetlist) {
02369 facetT *prevfacet, *list= *facetlist;
02370
02371
02372 trace4((qh ferr, "qh_prependfacet: prepend f%d before f%d\n",
02373 facet->id, list->id));
02374 prevfacet= list->previous;
02375 facet->previous= prevfacet;
02376 if (prevfacet)
02377 prevfacet->next= facet;
02378 list->previous= facet;
02379 facet->next= *facetlist;
02380 if (qh facet_list == list)
02381 qh facet_list= facet;
02382 if (qh facet_next == list)
02383 qh facet_next= facet;
02384 *facetlist= facet;
02385 qh num_facets++;
02386 }
02387
02388
02389
02390
02391
02392
02393
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403
02404 void qh_printhashtable(FILE *fp) {
02405 facetT *facet, *neighbor;
02406 int id, facet_i, facet_n, neighbor_i= 0, neighbor_n= 0;
02407 vertexT *vertex, **vertexp;
02408
02409 FOREACHfacet_i_(qh hash_table) {
02410 if (facet) {
02411 FOREACHneighbor_i_(facet) {
02412 if (!neighbor || neighbor == qh_MERGEridge || neighbor == qh_DUPLICATEridge)
02413 break;
02414 }
02415 if (neighbor_i == neighbor_n)
02416 continue;
02417 fprintf (fp, "hash %d f%d ", facet_i, facet->id);
02418 FOREACHvertex_(facet->vertices)
02419 fprintf (fp, "v%d ", vertex->id);
02420 fprintf (fp, "\n neighbors:");
02421 FOREACHneighbor_i_(facet) {
02422 if (neighbor == qh_MERGEridge)
02423 id= -3;
02424 else if (neighbor == qh_DUPLICATEridge)
02425 id= -2;
02426 else
02427 id= getid_(neighbor);
02428 fprintf (fp, " %d", id);
02429 }
02430 fprintf (fp, "\n");
02431 }
02432 }
02433 }
02434
02435
02436
02437
02438
02439
02440
02441
02442 void qh_printlists (void) {
02443 facetT *facet;
02444 vertexT *vertex;
02445
02446 fprintf (qh ferr, "qh_printlists: facets:");
02447 FORALLfacets
02448 fprintf (qh ferr, " %d", facet->id);
02449 fprintf (qh ferr, "\n new facets %d visible facets %d next facet for addpoint %d\n vertices (new %d):",
02450 getid_(qh newfacet_list), getid_(qh visible_list), getid_(qh facet_next),
02451 getid_(qh newvertex_list));
02452 FORALLvertices
02453 fprintf (qh ferr, " %d", vertex->id);
02454 fprintf (qh ferr, "\n");
02455 }
02456
02457
02458
02459
02460
02461
02462
02463
02464
02465
02466
02467
02468
02469 void qh_resetlists (boolT stats ) {
02470 vertexT *vertex;
02471 facetT *newfacet, *visible;
02472 int totnew=0, totver=0;
02473
02474 if (stats) {
02475 FORALLvertex_(qh newvertex_list)
02476 totver++;
02477 FORALLnew_facets
02478 totnew++;
02479 zadd_(Zvisvertextot, totver);
02480 zmax_(Zvisvertexmax, totver);
02481 zadd_(Znewfacettot, totnew);
02482 zmax_(Znewfacetmax, totnew);
02483 }
02484 FORALLvertex_(qh newvertex_list)
02485 vertex->newlist= False;
02486 qh newvertex_list= NULL;
02487 FORALLnew_facets
02488 newfacet->newfacet= False;
02489 qh newfacet_list= NULL;
02490 FORALLvisible_facets {
02491 visible->f.replace= NULL;
02492 visible->visible= False;
02493 }
02494 qh visible_list= NULL;
02495 qh num_visible= 0;
02496 qh NEWfacets= False;
02497 }
02498
02499
02500
02501
02502
02503
02504
02505
02506
02507
02508
02509
02510
02511
02512
02513
02514
02515
02516
02517 void qh_setvoronoi_all (void) {
02518 facetT *facet;
02519
02520 qh_clearcenters (qh_ASvoronoi);
02521 qh_vertexneighbors();
02522
02523 FORALLfacets {
02524 if (!facet->normal || !facet->upperdelaunay || qh UPPERdelaunay) {
02525 if (!facet->center)
02526 facet->center= qh_facetcenter (facet->vertices);
02527 }
02528 }
02529 }
02530
02531
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543
02544 void qh_vertexintersect(setT **vertexsetA,setT *vertexsetB) {
02545 setT *intersection;
02546
02547 intersection= qh_vertexintersect_new (*vertexsetA, vertexsetB);
02548 qh_settempfree (vertexsetA);
02549 *vertexsetA= intersection;
02550 qh_settemppush (intersection);
02551 }
02552
02553
02554
02555
02556
02557
02558
02559
02560
02561
02562 setT *qh_vertexintersect_new (setT *vertexsetA,setT *vertexsetB) {
02563 setT *intersection= qh_setnew (qh hull_dim - 1);
02564 vertexT **vertexA= SETaddr_(vertexsetA, vertexT);
02565 vertexT **vertexB= SETaddr_(vertexsetB, vertexT);
02566
02567 while (*vertexA && *vertexB) {
02568 if (*vertexA == *vertexB) {
02569 qh_setappend(&intersection, *vertexA);
02570 vertexA++; vertexB++;
02571 }else {
02572 if ((*vertexA)->id > (*vertexB)->id)
02573 vertexA++;
02574 else
02575 vertexB++;
02576 }
02577 }
02578 return intersection;
02579 }
02580
02581
02582
02583
02584
02585
02586
02587
02588
02589
02590
02591
02592
02593
02594
02595
02596
02597
02598
02599
02600
02601 void qh_vertexneighbors (void ) {
02602 facetT *facet;
02603 vertexT *vertex, **vertexp;
02604
02605 if (qh VERTEXneighbors)
02606 return;
02607 trace1((qh ferr, "qh_vertexneighbors: determing neighboring facets for each vertex\n"));
02608 qh vertex_visit++;
02609 FORALLfacets {
02610 if (facet->visible)
02611 continue;
02612 FOREACHvertex_(facet->vertices) {
02613 if (vertex->visitid != qh vertex_visit) {
02614 vertex->visitid= qh vertex_visit;
02615 vertex->neighbors= qh_setnew (qh hull_dim);
02616 }
02617 qh_setappend (&vertex->neighbors, facet);
02618 }
02619 }
02620 qh VERTEXneighbors= True;
02621 }
02622
02623
02624
02625
02626
02627
02628
02629
02630
02631
02632
02633 boolT qh_vertexsubset(setT *vertexsetA, setT *vertexsetB) {
02634 vertexT **vertexA= (vertexT **) SETaddr_(vertexsetA, vertexT);
02635 vertexT **vertexB= (vertexT **) SETaddr_(vertexsetB, vertexT);
02636
02637 while (True) {
02638 if (!*vertexA)
02639 return True;
02640 if (!*vertexB)
02641 return False;
02642 if ((*vertexA)->id > (*vertexB)->id)
02643 return False;
02644 if (*vertexA == *vertexB)
02645 vertexA++;
02646 vertexB++;
02647 }
02648 return False;
02649 }
