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
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 void qh_qhull (void) {
00058 int numoutside;
00059
00060 qh hulltime= qh_CPUclock;
00061 if (qh RERUN || qh JOGGLEmax < REALmax/2)
00062 qh_build_withrestart();
00063 else {
00064 qh_initbuild();
00065 qh_buildhull();
00066 }
00067 if (!qh STOPpoint && !qh STOPcone) {
00068 if (qh ZEROall_ok && !qh TESTvneighbors && qh MERGEexact)
00069 qh_checkzero( qh_ALL);
00070 if (qh ZEROall_ok && !qh TESTvneighbors && !qh WAScoplanar) {
00071 trace2((qh ferr, "qh_qhull: all facets are clearly convex and no coplanar points. Post-merging and check of maxout not needed.\n"));
00072 }else {
00073 if (qh MERGEexact || (qh hull_dim > qh_DIMreduceBuild && qh PREmerge))
00074 qh_postmerge ("First post-merge", qh premerge_centrum, qh premerge_cos,
00075 (qh POSTmerge ? False : qh TESTvneighbors));
00076 else if (!qh POSTmerge && qh TESTvneighbors)
00077 qh_postmerge ("For testing vertex neighbors", qh premerge_centrum,
00078 qh premerge_cos, True);
00079 if (qh POSTmerge)
00080 qh_postmerge ("For post-merging", qh postmerge_centrum,
00081 qh postmerge_cos, qh TESTvneighbors);
00082 if (qh visible_list == qh facet_list) {
00083 qh findbestnew= True;
00084 qh_partitionvisible ( !qh_ALL, &numoutside);
00085 qh findbestnew= False;
00086 qh_deletevisible ();
00087 qh_resetlists (False );
00088 }
00089 if (qh DOcheckmax){
00090 if (qh REPORTfreq) {
00091 qh_buildtracing (NULL, NULL);
00092 fprintf (qh ferr, "\nTesting all coplanar points.\n");
00093 }
00094 qh_check_maxout();
00095 }
00096 }
00097 }
00098 if (qh KEEPnearinside && !qh maxoutdone)
00099 qh_nearcoplanar();
00100 if (qh_setsize ((setT*)qhmem.tempstack) != 0) {
00101 fprintf (qh ferr, "qhull internal error (qh_qhull): temporary sets not empty (%d)\n",
00102 qh_setsize ((setT*)qhmem.tempstack));
00103 qh_errexit (qh_ERRqhull, NULL, NULL);
00104 }
00105 qh hulltime= qh_CPUclock - qh hulltime;
00106 qh QHULLfinished= True;
00107 trace1((qh ferr, "qh_qhull: algorithm completed\n"));
00108 }
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
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 boolT qh_addpoint (pointT *furthest, facetT *facet, boolT checkdist) {
00158 int goodvisible, goodhorizon;
00159 vertexT *vertex;
00160 facetT *newfacet;
00161 realT dist, newbalance, pbalance;
00162 boolT isoutside= False;
00163 int numpart, numpoints, numnew, firstnew;
00164
00165 qh maxoutdone= False;
00166 if (qh_pointid (furthest) == -1)
00167 qh_setappend (&qh other_points, furthest);
00168 if (!facet) {
00169 fprintf (qh ferr, "qh_addpoint: NULL facet. Use qh_findbestfacet\n");
00170 qh_errexit (qh_ERRqhull, NULL, NULL);
00171 }
00172 if (checkdist) {
00173 facet= qh_findbest (furthest, facet, !qh_ALL, False, !qh_NOupper,
00174 &dist, &isoutside, &numpart);
00175 zzadd_(Zpartition, numpart);
00176 if (!isoutside) {
00177 zinc_(Znotmax);
00178 facet->notfurthest= True;
00179 qh_partitioncoplanar (furthest, facet, &dist);
00180 return True;
00181 }
00182 }
00183 qh_buildtracing (furthest, facet);
00184 if (qh STOPpoint < 0 && qh furthest_id == -qh STOPpoint-1) {
00185 facet->notfurthest= True;
00186 return False;
00187 }
00188 qh_findhorizon (furthest, facet, &goodvisible, &goodhorizon);
00189 if (qh ONLYgood && !(goodvisible+goodhorizon) && !qh GOODclosest) {
00190 zinc_(Znotgood);
00191 facet->notfurthest= True;
00192
00193
00194 qh_resetlists (False );
00195 return True;
00196 }
00197 zzinc_(Zprocessed);
00198 firstnew= qh facet_id;
00199 vertex= qh_makenewfacets (furthest );
00200 qh_makenewplanes ();
00201 numnew= qh facet_id - firstnew;
00202 newbalance= numnew - (realT) (qh num_facets-qh num_visible)
00203 * qh hull_dim/qh num_vertices;
00204 wadd_(Wnewbalance, newbalance);
00205 wadd_(Wnewbalance2, newbalance * newbalance);
00206 if (qh ONLYgood
00207 && !qh_findgood (qh newfacet_list, goodhorizon) && !qh GOODclosest) {
00208 FORALLnew_facets
00209 qh_delfacet (newfacet);
00210 qh_delvertex (vertex);
00211 qh_resetlists (True );
00212 zinc_(Znotgoodnew);
00213 facet->notfurthest= True;
00214 return True;
00215 }
00216 if (qh ONLYgood)
00217 qh_attachnewfacets();
00218 qh_matchnewfacets();
00219 qh_updatevertices();
00220 if (qh STOPcone && qh furthest_id == qh STOPcone-1) {
00221 facet->notfurthest= True;
00222 return False;
00223 }
00224 if (qh PREmerge || qh MERGEexact) {
00225 qh_premerge (vertex, qh premerge_centrum, qh premerge_cos);
00226 if (zzval_(Ztotmerge) > qh_USEfindbestnew)
00227 qh findbestnew= True;
00228 else {
00229 FORALLnew_facets {
00230 if (!newfacet->simplicial) {
00231 qh findbestnew= True;
00232 break;
00233 }
00234 }
00235 }
00236 }else if (qh BESToutside)
00237 qh findbestnew= True;
00238 qh_partitionvisible ( !qh_ALL, &numpoints);
00239 qh findbestnew= False;
00240 qh findbest_notsharp= False;
00241 zinc_(Zpbalance);
00242 pbalance= numpoints - (realT) qh hull_dim
00243 * (qh num_points - qh num_vertices)/qh num_vertices;
00244 wadd_(Wpbalance, pbalance);
00245 wadd_(Wpbalance2, pbalance * pbalance);
00246 qh_deletevisible ();
00247 zmax_(Zmaxvertex, qh num_vertices);
00248 qh NEWfacets= False;
00249 if (qh IStracing >= 4)
00250 qh_printfacetlist (qh newfacet_list, NULL, True);
00251 if (qh CHECKfrequently) {
00252 if (qh num_facets < 50)
00253 qh_checkpolygon (qh facet_list);
00254 else
00255 qh_checkpolygon (qh newfacet_list);
00256 }
00257 if (qh STOPpoint > 0 && qh furthest_id == qh STOPpoint-1)
00258 return False;
00259 qh_resetlists (True );
00260 trace2((qh ferr, "qh_addpoint: added p%d new facets %d new balance %2.2g point balance %2.2g\n",
00261 qh_pointid (furthest), numnew, newbalance, pbalance));
00262 if (qh hull_dim > 3 && qh TRACEpoint == qh_pointid (furthest)) {
00263 qh IStracing= 0;
00264 qhmem.IStracing= 0;
00265 }
00266 return True;
00267 }
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277 void qh_build_withrestart (void) {
00278 int restart;
00279
00280 qh ALLOWrestart= True;
00281 while (True) {
00282 restart= setjmp (qh restartexit);
00283 if (restart) {
00284 zzinc_(Zretry);
00285 wmax_(Wretrymax, qh JOGGLEmax);
00286 qh ERREXITcalled= False;
00287 qh STOPcone= True;
00288 }
00289 if (!qh RERUN && qh JOGGLEmax < REALmax/2) {
00290 if (qh build_cnt > qh_JOGGLEmaxretry) {
00291 fprintf(qh ferr, "\n\
00292 qhull precision error: %d attempts to construct a convex hull\n\
00293 with joggled input. Increase joggle above 'QJ%2.2g'\n\
00294 or modify qh_JOGGLE... parameters in user.h\n",
00295 qh build_cnt, qh JOGGLEmax);
00296 qh_errexit (qh_ERRqhull, NULL, NULL);
00297 }
00298 if (qh build_cnt && !restart)
00299 break;
00300 }else if (qh build_cnt && qh build_cnt >= qh RERUN)
00301 break;
00302 qh STOPcone= False;
00303 qh_freebuild (True);
00304 qh build_cnt++;
00305 if (!qh qhull_optionsiz)
00306 qh qhull_optionsiz= strlen (qh qhull_options);
00307 else {
00308 qh qhull_options [qh qhull_optionsiz]= '\0';
00309 qh qhull_optionlen= 80;
00310 }
00311 qh_option("_run", &qh build_cnt, NULL);
00312 if (qh build_cnt == qh RERUN) {
00313 qh IStracing= qh TRACElastrun;
00314 if (qh TRACEpoint != -1 || qh TRACEdist < REALmax/2 || qh TRACEmerge) {
00315 qh TRACElevel= (qh IStracing? qh IStracing : 3);
00316 qh IStracing= 0;
00317 }
00318 qhmem.IStracing= qh IStracing;
00319 }
00320 if (qh JOGGLEmax < REALmax/2)
00321 qh_joggleinput();
00322 qh_initbuild();
00323 qh_buildhull();
00324 if (qh JOGGLEmax < REALmax/2 && !qh MERGING)
00325 qh_checkconvex (qh facet_list, qh_ALGORITHMfault);
00326 }
00327 qh ALLOWrestart= False;
00328 }
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352 void qh_buildhull(void) {
00353 facetT *facet;
00354 pointT *furthest;
00355 vertexT *vertex;
00356 int id;
00357
00358 trace1((qh ferr, "qh_buildhull: start build hull\n"));
00359 FORALLfacets {
00360 if (facet->visible || facet->newfacet) {
00361 fprintf (qh ferr, "qhull internal error (qh_buildhull): visible or new facet f%d in facet list\n",
00362 facet->id);
00363 qh_errexit (qh_ERRqhull, facet, NULL);
00364 }
00365 }
00366 FORALLvertices {
00367 if (vertex->newlist) {
00368 fprintf (qh ferr, "qhull internal error (qh_buildhull): new vertex f%d in vertex list\n",
00369 vertex->id);
00370 qh_errprint ("ERRONEOUS", NULL, NULL, NULL, vertex);
00371 qh_errexit (qh_ERRqhull, NULL, NULL);
00372 }
00373 id= qh_pointid (vertex->point);
00374 if ((qh STOPpoint>0 && id == qh STOPpoint-1) ||
00375 (qh STOPpoint<0 && id == -qh STOPpoint-1) ||
00376 (qh STOPcone>0 && id == qh STOPcone-1)) {
00377 trace1((qh ferr,"qh_buildhull: stop point or cone P%d in initial hull\n", id));
00378 return;
00379 }
00380 }
00381 qh facet_next= qh facet_list;
00382 while ((furthest= qh_nextfurthest (&facet))) {
00383 qh num_outside--;
00384 if (!qh_addpoint (furthest, facet, qh ONLYmax))
00385 break;
00386 }
00387 if (qh NARROWhull)
00388 qh_outcoplanar( );
00389 if (qh num_outside && !furthest) {
00390 fprintf (qh ferr, "qhull internal error (qh_buildhull): %d outside points were never processed.\n", qh num_outside);
00391 qh_errexit (qh_ERRqhull, NULL, NULL);
00392 }
00393 trace1((qh ferr, "qh_buildhull: completed the hull construction\n"));
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 void qh_buildtracing (pointT *furthest, facetT *facet) {
00424 realT dist= 0;
00425 float cpu;
00426 int total, furthestid;
00427 time_t timedata;
00428 struct tm *tp;
00429 vertexT *vertex;
00430
00431 qh old_randomdist= qh RANDOMdist;
00432 qh RANDOMdist= False;
00433 if (!furthest) {
00434 time (&timedata);
00435 tp= localtime (&timedata);
00436 cpu= qh_CPUclock - qh hulltime;
00437 cpu /= qh_SECticks;
00438 total= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
00439 fprintf (qh ferr, "\n\
00440 At %02d:%02d:%02d & %2.5g CPU secs, qhull has created %d facets and merged %d.\n\
00441 The current hull contains %d facets and %d vertices. Last point was p%d\n",
00442 tp->tm_hour, tp->tm_min, tp->tm_sec, cpu, qh facet_id -1,
00443 total, qh num_facets, qh num_vertices, qh furthest_id);
00444 return;
00445 }
00446 furthestid= qh_pointid (furthest);
00447 if (qh TRACEpoint == furthestid) {
00448 qh IStracing= qh TRACElevel;
00449 qhmem.IStracing= qh TRACElevel;
00450 }
00451 if (qh REPORTfreq && (qh facet_id-1 > qh lastreport+qh REPORTfreq)) {
00452 qh lastreport= qh facet_id-1;
00453 time (&timedata);
00454 tp= localtime (&timedata);
00455 cpu= qh_CPUclock - qh hulltime;
00456 cpu /= qh_SECticks;
00457 total= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
00458 zinc_(Zdistio);
00459 qh_distplane (furthest, facet, &dist);
00460 fprintf (qh ferr, "\n\
00461 At %02d:%02d:%02d & %2.5g CPU secs, qhull has created %d facets and merged %d.\n\
00462 The current hull contains %d facets and %d vertices. There are %d\n\
00463 outside points. Next is point p%d (v%d), %2.2g above f%d.\n",
00464 tp->tm_hour, tp->tm_min, tp->tm_sec, cpu, qh facet_id -1,
00465 total, qh num_facets, qh num_vertices, qh num_outside+1,
00466 furthestid, qh vertex_id, dist, getid_(facet));
00467 }else if (qh IStracing >=1) {
00468 cpu= qh_CPUclock - qh hulltime;
00469 cpu /= qh_SECticks;
00470 qh_distplane (furthest, facet, &dist);
00471 fprintf (qh ferr, "qh_addpoint: add p%d (v%d) to hull of %d facets (%2.2g above f%d) and %d outside at %4.4g CPU secs. Previous was p%d.\n",
00472 furthestid, qh vertex_id, qh num_facets, dist,
00473 getid_(facet), qh num_outside+1, cpu, qh furthest_id);
00474 }
00475 if (qh visit_id > (unsigned) INT_MAX) {
00476 qh visit_id= 0;
00477 FORALLfacets
00478 facet->visitid= qh visit_id;
00479 }
00480 if (qh vertex_visit > (unsigned) INT_MAX) {
00481 qh vertex_visit= 0;
00482 FORALLvertices
00483 vertex->visitid= qh vertex_visit;
00484 }
00485 qh furthest_id= furthestid;
00486 qh RANDOMdist= qh old_randomdist;
00487 }
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502 void qh_errexit2(int exitcode, facetT *facet, facetT *otherfacet) {
00503
00504 qh_errprint("ERRONEOUS", facet, otherfacet, NULL, NULL);
00505 qh_errexit (exitcode, NULL, NULL);
00506 }
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536 void qh_findhorizon(pointT *point, facetT *facet, int *goodvisible, int *goodhorizon) {
00537 facetT *neighbor, **neighborp, *visible;
00538 int numhorizon= 0, coplanar= 0;
00539 realT dist;
00540
00541 trace1((qh ferr,"qh_findhorizon: find horizon for point p%d facet f%d\n",qh_pointid(point),facet->id));
00542 *goodvisible= *goodhorizon= 0;
00543 zinc_(Ztotvisible);
00544 qh_removefacet(facet);
00545 qh_appendfacet(facet);
00546 qh num_visible= 1;
00547 if (facet->good)
00548 (*goodvisible)++;
00549 qh visible_list= facet;
00550 facet->visible= True;
00551 facet->f.replace= NULL;
00552 if (qh IStracing >=4)
00553 qh_errprint ("visible", facet, NULL, NULL, NULL);
00554 qh visit_id++;
00555 FORALLvisible_facets {
00556 visible->visitid= qh visit_id;
00557 FOREACHneighbor_(visible) {
00558 if (neighbor->visitid == qh visit_id)
00559 continue;
00560 neighbor->visitid= qh visit_id;
00561 zzinc_(Znumvisibility);
00562 qh_distplane(point, neighbor, &dist);
00563 if (dist > qh MINvisible) {
00564 zinc_(Ztotvisible);
00565 qh_removefacet(neighbor);
00566 qh_appendfacet(neighbor);
00567 neighbor->visible= True;
00568 neighbor->f.replace= NULL;
00569 qh num_visible++;
00570 if (neighbor->good)
00571 (*goodvisible)++;
00572 if (qh IStracing >=4)
00573 qh_errprint ("visible", neighbor, NULL, NULL, NULL);
00574 }else {
00575 if (dist > - qh MAXcoplanar) {
00576 neighbor->coplanar= True;
00577 zzinc_(Zcoplanarhorizon);
00578 qh_precision ("coplanar horizon");
00579 coplanar++;
00580 if (qh MERGING) {
00581 if (dist > 0) {
00582 maximize_(qh max_outside, dist);
00583 maximize_(qh max_vertex, dist);
00584 #if qh_MAXoutside
00585 maximize_(neighbor->maxoutside, dist);
00586 #endif
00587 }else
00588 minimize_(qh min_vertex, dist);
00589 }
00590 trace2((qh ferr, "qh_findhorizon: point p%d is coplanar to horizon f%d, dist=%2.7g < qh MINvisible (%2.7g)\n",
00591 qh_pointid(point), neighbor->id, dist, qh MINvisible));
00592 }else
00593 neighbor->coplanar= False;
00594 zinc_(Ztothorizon);
00595 numhorizon++;
00596 if (neighbor->good)
00597 (*goodhorizon)++;
00598 if (qh IStracing >=4)
00599 qh_errprint ("horizon", neighbor, NULL, NULL, NULL);
00600 }
00601 }
00602 }
00603 if (!numhorizon) {
00604 qh_precision ("empty horizon");
00605 fprintf(qh ferr, "qhull precision error (qh_findhorizon): empty horizon\n\
00606 Point p%d was above all facets.\n", qh_pointid(point));
00607 qh_printfacetlist (qh facet_list, NULL, True);
00608 qh_errexit(qh_ERRprec, NULL, NULL);
00609 }
00610 trace1((qh ferr, "qh_findhorizon: %d horizon facets (good %d), %d visible (good %d), %d coplanar\n",
00611 numhorizon, *goodhorizon, qh num_visible, *goodvisible, coplanar));
00612 if (qh IStracing >= 4 && qh num_facets < 50)
00613 qh_printlists ();
00614 }
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638 pointT *qh_nextfurthest (facetT **visible) {
00639 facetT *facet;
00640 int size, index;
00641 realT randr, dist;
00642 pointT *furthest;
00643
00644 while ((facet= qh facet_next) != qh facet_tail) {
00645 if (!facet->outsideset) {
00646 qh facet_next= facet->next;
00647 continue;
00648 }
00649 SETreturnsize_(facet->outsideset, size);
00650 if (!size) {
00651 qh_setfree (&facet->outsideset);
00652 qh facet_next= facet->next;
00653 continue;
00654 }
00655 if (qh NARROWhull) {
00656 if (facet->notfurthest)
00657 qh_furthestout (facet);
00658 furthest= (pointT*)qh_setlast (facet->outsideset);
00659 #if qh_COMPUTEfurthest
00660 qh_distplane (furthest, facet, &dist);
00661 zinc_(Zcomputefurthest);
00662 #else
00663 dist= facet->furthestdist;
00664 #endif
00665 if (dist < qh MINoutside) {
00666 qh facet_next= facet->next;
00667 continue;
00668 }
00669 }
00670 if (!qh RANDOMoutside && !qh VIRTUALmemory) {
00671 if (qh PICKfurthest) {
00672 qh_furthestnext ();
00673 facet= qh facet_next;
00674 }
00675 *visible= facet;
00676 return ((pointT*)qh_setdellast (facet->outsideset));
00677 }
00678 if (qh RANDOMoutside) {
00679 int outcoplanar = 0;
00680 if (qh NARROWhull) {
00681 FORALLfacets {
00682 if (facet == qh facet_next)
00683 break;
00684 if (facet->outsideset)
00685 outcoplanar += qh_setsize( facet->outsideset);
00686 }
00687 }
00688 randr= qh_RANDOMint;
00689 randr= randr/(qh_RANDOMmax+1);
00690 index= (int)floor((qh num_outside - outcoplanar) * randr);
00691 FORALLfacet_(qh facet_next) {
00692 if (facet->outsideset) {
00693 SETreturnsize_(facet->outsideset, size);
00694 if (!size)
00695 qh_setfree (&facet->outsideset);
00696 else if (size > index) {
00697 *visible= facet;
00698 return ((pointT*)qh_setdelnth (facet->outsideset, index));
00699 }else
00700 index -= size;
00701 }
00702 }
00703 fprintf (qh ferr, "qhull internal error (qh_nextfurthest): num_outside %d is too low\nby at least %d, or a random real %g >= 1.0\n",
00704 qh num_outside, index+1, randr);
00705 qh_errexit (qh_ERRqhull, NULL, NULL);
00706 }else {
00707 facet= qh facet_tail->previous;
00708 if (!(furthest= (pointT*)qh_setdellast(facet->outsideset))) {
00709 if (facet->outsideset)
00710 qh_setfree (&facet->outsideset);
00711 qh_removefacet (facet);
00712 qh_prependfacet (facet, &qh facet_list);
00713 continue;
00714 }
00715 *visible= facet;
00716 return furthest;
00717 }
00718 }
00719 return NULL;
00720 }
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754 void qh_partitionall(setT *vertices, pointT *points, int numpoints){
00755 setT *pointset;
00756 vertexT *vertex, **vertexp;
00757 pointT *point, **pointp, *bestpoint;
00758 int size, point_i, point_n, point_end, remaining, i, id;
00759 facetT *facet;
00760 realT bestdist= -REALmax, dist, distoutside;
00761
00762 trace1((qh ferr, "qh_partitionall: partition all points into outside sets\n"));
00763 pointset= qh_settemp (numpoints);
00764 qh num_outside= 0;
00765 pointp= SETaddr_(pointset, pointT);
00766 for (i=numpoints, point= points; i--; point += qh hull_dim)
00767 *(pointp++)= point;
00768 qh_settruncate (pointset, numpoints);
00769 FOREACHvertex_(vertices) {
00770 if ((id= qh_pointid(vertex->point)) >= 0)
00771 SETelem_(pointset, id)= NULL;
00772 }
00773 id= qh_pointid (qh GOODpointp);
00774 if (id >=0 && qh STOPcone-1 != id && -qh STOPpoint-1 != id)
00775 SETelem_(pointset, id)= NULL;
00776 if (qh GOODvertexp && qh ONLYgood && !qh MERGING) {
00777 if ((id= qh_pointid(qh GOODvertexp)) >= 0)
00778 SETelem_(pointset, id)= NULL;
00779 }
00780 if (!qh BESToutside) {
00781 if (qh MERGING)
00782 distoutside= qh_DISToutside;
00783 else
00784 distoutside= qh MINoutside;
00785 zval_(Ztotpartition)= qh num_points - qh hull_dim - 1;
00786 remaining= qh num_facets;
00787 point_end= numpoints;
00788 FORALLfacets {
00789 size= point_end/(remaining--) + 100;
00790 facet->outsideset= qh_setnew (size);
00791 bestpoint= NULL;
00792 point_end= 0;
00793 FOREACHpoint_i_(pointset) {
00794 if (point) {
00795 zzinc_(Zpartitionall);
00796 qh_distplane (point, facet, &dist);
00797 if (dist < distoutside)
00798 SETelem_(pointset, point_end++)= point;
00799 else {
00800 qh num_outside++;
00801 if (!bestpoint) {
00802 bestpoint= point;
00803 bestdist= dist;
00804 }else if (dist > bestdist) {
00805 qh_setappend (&facet->outsideset, bestpoint);
00806 bestpoint= point;
00807 bestdist= dist;
00808 }else
00809 qh_setappend (&facet->outsideset, point);
00810 }
00811 }
00812 }
00813 if (bestpoint) {
00814 qh_setappend (&facet->outsideset, bestpoint);
00815 #if !qh_COMPUTEfurthest
00816 facet->furthestdist= bestdist;
00817 #endif
00818 }else
00819 qh_setfree (&facet->outsideset);
00820 qh_settruncate (pointset, point_end);
00821 }
00822 }
00823 if (qh BESToutside || qh MERGING || qh KEEPcoplanar || qh KEEPinside) {
00824 qh findbestnew= True;
00825 FOREACHpoint_i_(pointset) {
00826 if (point)
00827 qh_partitionpoint(point, qh facet_list);
00828 }
00829 qh findbestnew= False;
00830 }
00831 zzadd_(Zpartitionall, zzval_(Zpartition));
00832 zzval_(Zpartition)= 0;
00833 qh_settempfree(&pointset);
00834 if (qh IStracing >= 4)
00835 qh_printfacetlist (qh facet_list, NULL, True);
00836 }
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872 void qh_partitioncoplanar (pointT *point, facetT *facet, realT *dist) {
00873 facetT *bestfacet;
00874 pointT *oldfurthest;
00875 realT bestdist, dist2;
00876 int numpart= 0;
00877 boolT isoutside, istrace= False;
00878
00879 qh WAScoplanar= True;
00880 if (!dist) {
00881 if (qh findbestnew)
00882 bestfacet= qh_findbestnew (point, facet,
00883 &bestdist, NULL, &numpart);
00884 else
00885 bestfacet= qh_findbest (point, facet, qh_ALL, False, !qh_NOupper,
00886 &bestdist, &isoutside, &numpart);
00887 zinc_(Ztotpartcoplanar);
00888 zzadd_(Zpartcoplanar, numpart);
00889 if (!qh KEEPinside) {
00890 if (qh KEEPnearinside) {
00891 if (bestdist < -qh NEARinside) {
00892 zinc_(Zcoplanarinside);
00893 return;
00894 }
00895 }else if (bestdist < -qh MAXcoplanar) {
00896 zinc_(Zcoplanarinside);
00897 return;
00898 }
00899 }
00900 }else {
00901 bestfacet= facet;
00902 bestdist= *dist;
00903 }
00904 if (qh KEEPcoplanar + qh KEEPinside + qh KEEPnearinside) {
00905 oldfurthest= (pointT*)qh_setlast (bestfacet->coplanarset);
00906 if (oldfurthest) {
00907 zinc_(Zcomputefurthest);
00908 qh_distplane (oldfurthest, bestfacet, &dist2);
00909 }
00910 if (!oldfurthest || dist2 < bestdist) {
00911 qh_setappend(&bestfacet->coplanarset, point);
00912 if (bestdist > qh max_outside) {
00913 qh max_outside= bestdist;
00914 if (bestdist > qh TRACEdist)
00915 istrace= True;
00916 }
00917 }else
00918 qh_setappend2ndlast(&bestfacet->coplanarset, point);
00919 }else {
00920 if (bestdist > qh max_outside) {
00921 qh max_outside= bestdist;
00922 if (bestdist > qh TRACEdist)
00923 istrace= True;
00924 }
00925 }
00926 if (istrace) {
00927 fprintf (qh ferr, "qh_partitioncoplanar: ====== p%d increases max_outside to %2.2g of f%d last p%d\n",
00928 qh_pointid(point), bestdist, bestfacet->id, qh furthest_id);
00929 qh_errprint ("DISTANT", bestfacet, NULL, NULL, NULL);
00930 }
00931 trace4((qh ferr, "qh_partitioncoplanar: point p%d is coplanar with facet f%d (or inside) dist %2.2g\n",
00932 qh_pointid(point), bestfacet->id, bestdist));
00933 }
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970 void qh_partitionpoint (pointT *point, facetT *facet) {
00971 realT bestdist;
00972 boolT isoutside;
00973 facetT *bestfacet;
00974 int numpart;
00975 #if qh_COMPUTEfurthest
00976 realT dist;
00977 #endif
00978
00979 if (qh findbestnew)
00980 bestfacet= qh_findbestnew (point, facet,
00981 &bestdist, &isoutside, &numpart);
00982 else
00983 bestfacet= qh_findbest (point, facet, qh BESToutside, True, !qh_NOupper,
00984 &bestdist, &isoutside, &numpart);
00985 zinc_(Ztotpartition);
00986 zzadd_(Zpartition, numpart);
00987 if (qh NARROWhull) {
00988 if (qh DELAUNAY && !isoutside && bestdist >= -qh MAXcoplanar)
00989 qh_precision ("nearly incident point (narrow hull)");
00990 if (qh KEEPnearinside) {
00991 if (bestdist >= -qh NEARinside)
00992 isoutside= True;
00993 }else if (bestdist >= -qh MAXcoplanar)
00994 isoutside= True;
00995 }
00996
00997 if (isoutside) {
00998 if (!bestfacet->outsideset
00999 || !qh_setlast (bestfacet->outsideset)) {
01000 qh_setappend(&(bestfacet->outsideset), point);
01001 if (!bestfacet->newfacet) {
01002 qh_removefacet (bestfacet);
01003 qh_appendfacet (bestfacet);
01004 }
01005 #if !qh_COMPUTEfurthest
01006 bestfacet->furthestdist= bestdist;
01007 #endif
01008 }else {
01009 #if qh_COMPUTEfurthest
01010 zinc_(Zcomputefurthest);
01011 qh_distplane (oldfurthest, bestfacet, &dist);
01012 if (dist < bestdist)
01013 qh_setappend(&(bestfacet->outsideset), point);
01014 else
01015 qh_setappend2ndlast(&(bestfacet->outsideset), point);
01016 #else
01017 if (bestfacet->furthestdist < bestdist) {
01018 qh_setappend(&(bestfacet->outsideset), point);
01019 bestfacet->furthestdist= bestdist;
01020 }else
01021 qh_setappend2ndlast(&(bestfacet->outsideset), point);
01022 #endif
01023 }
01024 qh num_outside++;
01025 trace4((qh ferr, "qh_partitionpoint: point p%d is outside facet f%d\n",
01026 qh_pointid(point), bestfacet->id));
01027 }else if (bestdist >= -qh MAXcoplanar) {
01028 zzinc_(Zcoplanarpart);
01029 if (qh DELAUNAY)
01030 qh_precision ("nearly incident point");
01031 if (qh KEEPcoplanar + qh KEEPnearinside || bestdist > qh max_outside)
01032 qh_partitioncoplanar (point, bestfacet, &bestdist);
01033 }else if (qh KEEPnearinside && bestdist > -qh NEARinside) {
01034 zinc_(Zpartnear);
01035 qh_partitioncoplanar (point, bestfacet, &bestdist);
01036 }else {
01037 zinc_(Zpartinside);
01038 trace4((qh ferr, "qh_partitionpoint: point p%d is inside all facets, closest to f%d dist %2.2g\n",
01039 qh_pointid(point), bestfacet->id, bestdist));
01040 if (qh KEEPinside)
01041 qh_partitioncoplanar (point, bestfacet, &bestdist);
01042 }
01043 }
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079 void qh_partitionvisible( boolT allpoints, int *numoutside) {
01080 facetT *visible, *newfacet;
01081 pointT *point, **pointp;
01082 int coplanar=0, size;
01083 unsigned count;
01084 vertexT *vertex, **vertexp;
01085
01086 if (qh ONLYmax)
01087 maximize_(qh MINoutside, qh max_vertex);
01088 *numoutside= 0;
01089 FORALLvisible_facets {
01090 if (!visible->outsideset && !visible->coplanarset)
01091 continue;
01092 newfacet= visible->f.replace;
01093 count= 0;
01094 while (newfacet && newfacet->visible) {
01095 newfacet= newfacet->f.replace;
01096 if (count++ > qh facet_id)
01097 qh_infiniteloop (visible);
01098 }
01099 if (!newfacet)
01100 newfacet= qh newfacet_list;
01101 if (visible->outsideset) {
01102 size= qh_setsize (visible->outsideset);
01103 *numoutside += size;
01104 qh num_outside -= size;
01105 FOREACHpoint_(visible->outsideset)
01106 qh_partitionpoint (point, newfacet);
01107 }
01108 if (visible->coplanarset && (qh KEEPcoplanar + qh KEEPinside + qh KEEPnearinside)) {
01109 size= qh_setsize (visible->coplanarset);
01110 coplanar += size;
01111 FOREACHpoint_(visible->coplanarset) {
01112 if (allpoints)
01113 qh_partitionpoint (point, newfacet);
01114 else
01115 qh_partitioncoplanar (point, newfacet, NULL);
01116 }
01117 }
01118 }
01119 FOREACHvertex_(qh del_vertices) {
01120 if (vertex->point) {
01121 if (allpoints)
01122 qh_partitionpoint (vertex->point, qh newfacet_list);
01123 else
01124 qh_partitioncoplanar (vertex->point, qh newfacet_list, NULL);
01125 }
01126 }
01127 trace1((qh ferr,"qh_partitionvisible: partitioned %d points from outsidesets and %d points from coplanarsets\n", *numoutside, coplanar));
01128 }
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138 void qh_precision (char *reason) {
01139
01140 if (qh ALLOWrestart && !qh PREmerge && !qh MERGEexact) {
01141 if (qh JOGGLEmax < REALmax/2) {
01142 trace0((qh ferr, "qh_precision: qhull restart because of %s\n", reason));
01143 longjmp(qh restartexit, qh_ERRprec);
01144 }
01145 }
01146 }
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161 void qh_printsummary(FILE *fp) {
01162 realT ratio, outerplane, innerplane;
01163 float cpu;
01164 int size, id, nummerged, numvertices, numcoplanars= 0, nonsimplicial=0;
01165 int goodused;
01166 facetT *facet;
01167 char *s;
01168
01169 size= qh num_points + qh_setsize (qh other_points);
01170 numvertices= qh num_vertices - qh_setsize (qh del_vertices);
01171 id= qh_pointid (qh GOODpointp);
01172 FORALLfacets {
01173 if (facet->coplanarset)
01174 numcoplanars += qh_setsize( facet->coplanarset);
01175 if (facet->good) {
01176 if (!facet->simplicial && qh_setsize(facet->vertices) != qh hull_dim)
01177 nonsimplicial++;
01178 }
01179 }
01180 if (id >=0 && qh STOPcone-1 != id && -qh STOPpoint-1 != id)
01181 size--;
01182 if (qh STOPcone || qh STOPpoint)
01183 fprintf (fp, "\nAt a premature exit due to 'TVn', 'TCn', or precision error.");
01184 if (qh UPPERdelaunay)
01185 goodused= qh GOODvertex + qh GOODpoint + qh SPLITthresholds;
01186 else if (qh DELAUNAY)
01187 goodused= qh GOODvertex + qh GOODpoint + qh GOODthreshold;
01188 else
01189 goodused= qh num_good;
01190 nummerged= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
01191 if (qh VORONOI) {
01192 if (qh UPPERdelaunay)
01193 fprintf (fp, "\n\
01194 Furthest-site Voronoi vertices by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
01195 else
01196 fprintf (fp, "\n\
01197 Voronoi diagram by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
01198 fprintf(fp, " Number of Voronoi regions%s: %d\n",
01199 qh ATinfinity ? " and at-infinity" : "", numvertices);
01200 if (numcoplanars)
01201 fprintf(fp, " Number of nearly incident points: %d\n", numcoplanars);
01202 else if (size > numvertices)
01203 fprintf(fp, " Total number of nearly incident points: %d\n", size - numvertices);
01204 fprintf(fp, " Number of%s Voronoi vertices: %d\n",
01205 goodused ? " 'good'" : "", qh num_good);
01206 if (nonsimplicial)
01207 fprintf(fp, " Number of%s non-simplicial Voronoi vertices: %d\n",
01208 goodused ? " 'good'" : "", nonsimplicial);
01209 }else if (qh DELAUNAY) {
01210 if (qh UPPERdelaunay)
01211 fprintf (fp, "\n\
01212 Furthest-site Delaunay triangulation by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
01213 else
01214 fprintf (fp, "\n\
01215 Delaunay triangulation by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
01216 fprintf(fp, " Number of input sites%s: %d\n",
01217 qh ATinfinity ? " and at-infinity" : "", numvertices);
01218 if (numcoplanars)
01219 fprintf(fp, " Number of nearly incident points: %d\n", numcoplanars);
01220 else if (size > numvertices)
01221 fprintf(fp, " Total number of nearly incident points: %d\n", size - numvertices);
01222 fprintf(fp, " Number of%s Delaunay regions: %d\n",
01223 goodused ? " 'good'" : "", qh num_good);
01224 if (nonsimplicial)
01225 fprintf(fp, " Number of%s non-simplicial Delaunay regions: %d\n",
01226 goodused ? " 'good'" : "", nonsimplicial);
01227 }else if (qh HALFspace) {
01228 fprintf (fp, "\n\
01229 Halfspace intersection by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
01230 fprintf(fp, " Number of halfspaces: %d\n", size);
01231 fprintf(fp, " Number of non-redundant halfspaces: %d\n", numvertices);
01232 if (numcoplanars) {
01233 if (qh KEEPinside && qh KEEPcoplanar)
01234 s= "similar and redundant";
01235 else if (qh KEEPinside)
01236 s= "redundant";
01237 else
01238 s= "similar";
01239 fprintf(fp, " Number of %s halfspaces: %d\n", s, numcoplanars);
01240 }
01241 fprintf(fp, " Number of intersection points: %d\n", qh num_facets - qh num_visible);
01242 if (goodused)
01243 fprintf(fp, " Number of 'good' intersection points: %d\n", qh num_good);
01244 if (nonsimplicial)
01245 fprintf(fp, " Number of%s non-simplicial intersection points: %d\n",
01246 goodused ? " 'good'" : "", nonsimplicial);
01247 }else {
01248 fprintf (fp, "\n\
01249 Convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
01250 fprintf(fp, " Number of vertices: %d\n", numvertices);
01251 if (numcoplanars) {
01252 if (qh KEEPinside && qh KEEPcoplanar)
01253 s= "coplanar and interior";
01254 else if (qh KEEPinside)
01255 s= "interior";
01256 else
01257 s= "coplanar";
01258 fprintf(fp, " Number of %s points: %d\n", s, numcoplanars);
01259 }
01260 fprintf(fp, " Number of facets: %d\n", qh num_facets - qh num_visible);
01261 if (goodused)
01262 fprintf(fp, " Number of 'good' facets: %d\n", qh num_good);
01263 if (nonsimplicial)
01264 fprintf(fp, " Number of%s non-simplicial facets: %d\n",
01265 goodused ? " 'good'" : "", nonsimplicial);
01266 }
01267 fprintf(fp, "\nStatistics for: %s | %s",
01268 qh rbox_command, qh qhull_command);
01269 if (qh ROTATErandom != INT_MIN)
01270 fprintf(fp, " QR%d\n\n", qh ROTATErandom);
01271 else
01272 fprintf(fp, "\n\n");
01273 fprintf(fp, " Number of points processed: %d\n", zzval_(Zprocessed));
01274 fprintf(fp, " Number of hyperplanes created: %d\n", zzval_(Zsetplane));
01275 if (qh DELAUNAY)
01276 fprintf(fp, " Number of facets in hull: %d\n", qh num_facets - qh num_visible);
01277 fprintf(fp, " Number of distance tests for qhull: %d\n", zzval_(Zpartition)+
01278 zzval_(Zpartitionall)+zzval_(Znumvisibility)+zzval_(Zpartcoplanar));
01279 #if 0
01280 {realT stddev, ave;
01281 fprintf(fp, " average new facet balance: %2.2g\n",
01282 wval_(Wnewbalance)/zval_(Zprocessed));
01283 stddev= qh_stddev (zval_(Zprocessed), wval_(Wnewbalance),
01284 wval_(Wnewbalance2), &ave);
01285 fprintf(fp, " new facet standard deviation: %2.2g\n", stddev);
01286 fprintf(fp, " average partition balance: %2.2g\n",
01287 wval_(Wpbalance)/zval_(Zpbalance));
01288 stddev= qh_stddev (zval_(Zpbalance), wval_(Wpbalance),
01289 wval_(Wpbalance2), &ave);
01290 fprintf(fp, " partition standard deviation: %2.2g\n", stddev);
01291 }
01292 #endif
01293 if (nummerged) {
01294 fprintf(fp," Number of merged facets: %d\n", nummerged);
01295 fprintf(fp," Number of distance tests for merging: %d\n",zzval_(Zbestdist)+
01296 zzval_(Zcentrumtests)+zzval_(Zdistconvex)+zzval_(Zdistcheck)+
01297 zzval_(Zdistzero));
01298 }
01299 if (!qh RANDOMoutside && qh QHULLfinished) {
01300 cpu= qh hulltime;
01301 cpu /= qh_SECticks;
01302 wval_(Wcpu)= cpu;
01303 fprintf (fp, " CPU seconds to compute hull (after input): %2.4g\n", cpu);
01304 }
01305 if (qh RERUN) {
01306 if (!qh PREmerge && !qh MERGEexact)
01307 fprintf(fp, " Percentage of runs with precision errors: %4.1f\n",
01308 zzval_(Zretry)*100.0/qh build_cnt);
01309 }else if (qh JOGGLEmax < REALmax/2) {
01310 if (zzval_(Zretry))
01311 fprintf(fp, " After %d retries, input joggled by: %2.2g\n",
01312 zzval_(Zretry), qh JOGGLEmax);
01313 else
01314 fprintf(fp, " Input joggled by: %2.2g\n", qh JOGGLEmax);
01315 }
01316 if (qh totarea != 0.0)
01317 fprintf(fp, " %s facet area: %2.8g\n",
01318 zzval_(Ztotmerge) ? "Approximate" : "Total", qh totarea);
01319 if (qh totvol != 0.0)
01320 fprintf(fp, " %s volume: %2.8g\n",
01321 zzval_(Ztotmerge) ? "Approximate" : "Total", qh totvol);
01322 if (qh MERGING) {
01323 qh_outerinner (NULL, &outerplane, &innerplane);
01324 if (outerplane > 2 * qh DISTround) {
01325 fprintf(fp, " Maximum distance of %spoint above facet: %2.2g",
01326 (qh QHULLfinished ? "" : "merged "), outerplane);
01327 ratio= outerplane/(qh ONEmerge + qh DISTround);
01328 if (ratio > 0.05 && qh ONEmerge > qh MINoutside && qh JOGGLEmax > REALmax/2)
01329 fprintf (fp, " (%.1fx)\n", ratio);
01330 else
01331 fprintf (fp, "\n");
01332 }
01333 if (innerplane < -2 * qh DISTround) {
01334 fprintf(fp, " Maximum distance of %svertex below facet: %2.2g",
01335 (qh QHULLfinished ? "" : "merged "), innerplane);
01336 ratio= -innerplane/(qh ONEmerge+qh DISTround);
01337 if (ratio > 0.05 && qh JOGGLEmax > REALmax/2)
01338 fprintf (fp, " (%.1fx)\n", ratio);
01339 else
01340 fprintf (fp, "\n");
01341 }
01342 }
01343 fprintf(fp, "\n");
01344 }
01345
01346