libqhull.c
Go to the documentation of this file.
00001 /*<html><pre>  -<a                             href="qh-qhull.htm"
00002   >-------------------------------</a><a name="TOP">-</a>
00003 
00004    libqhull.c
00005    Quickhull algorithm for convex hulls
00006 
00007    qhull() and top-level routines
00008 
00009    see qh-qhull.htm, libqhull.h, unix.c
00010 
00011    see qhull_a.h for internal functions
00012 
00013    Copyright (c) 1993-2011 The Geometry Center.
00014    $Id: //main/2011/qhull/src/libqhull/libqhull.c#2 $$Change: 1342 $
00015    $DateTime: 2011/03/07 21:55:47 $$Author: bbarber $
00016 */
00017 
00018 #include "qhull_a.h"
00019 
00020 /*============= functions in alphabetic order after qhull() =======*/
00021 
00022 /*-<a                             href="qh-qhull.htm#TOC"
00023   >-------------------------------</a><a name="qhull">-</a>
00024 
00025   qh_qhull()
00026     compute DIM3 convex hull of qh.num_points starting at qh.first_point
00027     qh contains all global options and variables
00028 
00029   returns:
00030     returns polyhedron
00031       qh.facet_list, qh.num_facets, qh.vertex_list, qh.num_vertices,
00032 
00033     returns global variables
00034       qh.hulltime, qh.max_outside, qh.interior_point, qh.max_vertex, qh.min_vertex
00035 
00036     returns precision constants
00037       qh.ANGLEround, centrum_radius, cos_max, DISTround, MAXabs_coord, ONEmerge
00038 
00039   notes:
00040     unless needed for output
00041       qh.max_vertex and qh.min_vertex are max/min due to merges
00042 
00043   see:
00044     to add individual points to either qh.num_points
00045       use qh_addpoint()
00046 
00047     if qh.GETarea
00048       qh_produceoutput() returns qh.totarea and qh.totvol via qh_getarea()
00049 
00050   design:
00051     record starting time
00052     initialize hull and partition points
00053     build convex hull
00054     unless early termination
00055       update facet->maxoutside for vertices, coplanar, and near-inside points
00056     error if temporary sets exist
00057     record end time
00058 */
00059 
00060 void qh_qhull(void) {
00061   int numoutside;
00062 
00063   qh hulltime= qh_CPUclock;
00064   if (qh RERUN || qh JOGGLEmax < REALmax/2)
00065     qh_build_withrestart();
00066   else {
00067     qh_initbuild();
00068     qh_buildhull();
00069   }
00070   if (!qh STOPpoint && !qh STOPcone) {
00071     if (qh ZEROall_ok && !qh TESTvneighbors && qh MERGEexact)
00072       qh_checkzero( qh_ALL);
00073     if (qh ZEROall_ok && !qh TESTvneighbors && !qh WAScoplanar) {
00074       trace2((qh ferr, 2055, "qh_qhull: all facets are clearly convex and no coplanar points.  Post-merging and check of maxout not needed.\n"));
00075       qh DOcheckmax= False;
00076     }else {
00077       if (qh MERGEexact || (qh hull_dim > qh_DIMreduceBuild && qh PREmerge))
00078         qh_postmerge("First post-merge", qh premerge_centrum, qh premerge_cos,
00079              (qh POSTmerge ? False : qh TESTvneighbors));
00080       else if (!qh POSTmerge && qh TESTvneighbors)
00081         qh_postmerge("For testing vertex neighbors", qh premerge_centrum,
00082              qh premerge_cos, True);
00083       if (qh POSTmerge)
00084         qh_postmerge("For post-merging", qh postmerge_centrum,
00085              qh postmerge_cos, qh TESTvneighbors);
00086       if (qh visible_list == qh facet_list) { /* i.e., merging done */
00087         qh findbestnew= True;
00088         qh_partitionvisible(/*visible_list, newfacet_list*/ !qh_ALL, &numoutside);
00089         qh findbestnew= False;
00090         qh_deletevisible(/*qh visible_list*/);
00091         qh_resetlists(False, qh_RESETvisible /*qh visible_list newvertex_list newfacet_list */);
00092       }
00093     }
00094     if (qh DOcheckmax){
00095       if (qh REPORTfreq) {
00096         qh_buildtracing(NULL, NULL);
00097         qh_fprintf(qh ferr, 8115, "\nTesting all coplanar points.\n");
00098       }
00099       qh_check_maxout();
00100     }
00101     if (qh KEEPnearinside && !qh maxoutdone)
00102       qh_nearcoplanar();
00103   }
00104   if (qh_setsize(qhmem.tempstack) != 0) {
00105     qh_fprintf(qh ferr, 6164, "qhull internal error (qh_qhull): temporary sets not empty(%d)\n",
00106              qh_setsize(qhmem.tempstack));
00107     qh_errexit(qh_ERRqhull, NULL, NULL);
00108   }
00109   qh hulltime= qh_CPUclock - qh hulltime;
00110   qh QHULLfinished= True;
00111   trace1((qh ferr, 1036, "Qhull: algorithm completed\n"));
00112 } /* qhull */
00113 
00114 /*-<a                             href="qh-qhull.htm#TOC"
00115   >-------------------------------</a><a name="addpoint">-</a>
00116 
00117   qh_addpoint( furthest, facet, checkdist )
00118     add point (usually furthest point) above facet to hull
00119     if checkdist,
00120       check that point is above facet.
00121       if point is not outside of the hull, uses qh_partitioncoplanar()
00122       assumes that facet is defined by qh_findbestfacet()
00123     else if facet specified,
00124       assumes that point is above facet (major damage if below)
00125     for Delaunay triangulations,
00126       Use qh_setdelaunay() to lift point to paraboloid and scale by 'Qbb' if needed
00127       Do not use options 'Qbk', 'QBk', or 'QbB' since they scale the coordinates.
00128 
00129   returns:
00130     returns False if user requested an early termination
00131      qh.visible_list, newfacet_list, delvertex_list, NEWfacets may be defined
00132     updates qh.facet_list, qh.num_facets, qh.vertex_list, qh.num_vertices
00133     clear qh.maxoutdone (will need to call qh_check_maxout() for facet->maxoutside)
00134     if unknown point, adds a pointer to qh.other_points
00135       do not deallocate the point's coordinates
00136 
00137   notes:
00138     assumes point is near its best facet and not at a local minimum of a lens
00139       distributions.  Use qh_findbestfacet to avoid this case.
00140     uses qh.visible_list, qh.newfacet_list, qh.delvertex_list, qh.NEWfacets
00141 
00142   see also:
00143     qh_triangulate() -- triangulate non-simplicial facets
00144 
00145   design:
00146     add point to other_points if needed
00147     if checkdist
00148       if point not above facet
00149         partition coplanar point
00150         exit
00151     exit if pre STOPpoint requested
00152     find horizon and visible facets for point
00153     make new facets for point to horizon
00154     make hyperplanes for point
00155     compute balance statistics
00156     match neighboring new facets
00157     update vertex neighbors and delete interior vertices
00158     exit if STOPcone requested
00159     merge non-convex new facets
00160     if merge found, many merges, or 'Qf'
00161        use qh_findbestnew() instead of qh_findbest()
00162     partition outside points from visible facets
00163     delete visible facets
00164     check polyhedron if requested
00165     exit if post STOPpoint requested
00166     reset working lists of facets and vertices
00167 */
00168 boolT qh_addpoint(pointT *furthest, facetT *facet, boolT checkdist) {
00169   int goodvisible, goodhorizon;
00170   vertexT *vertex;
00171   facetT *newfacet;
00172   realT dist, newbalance, pbalance;
00173   boolT isoutside= False;
00174   int numpart, numpoints, numnew, firstnew;
00175 
00176   qh maxoutdone= False;
00177   if (qh_pointid(furthest) == -1)
00178     qh_setappend(&qh other_points, furthest);
00179   if (!facet) {
00180     qh_fprintf(qh ferr, 6213, "qhull internal error (qh_addpoint): NULL facet.  Need to call qh_findbestfacet first\n");
00181     qh_errexit(qh_ERRqhull, NULL, NULL);
00182   }
00183   if (checkdist) {
00184     facet= qh_findbest(furthest, facet, !qh_ALL, !qh_ISnewfacets, !qh_NOupper,
00185                         &dist, &isoutside, &numpart);
00186     zzadd_(Zpartition, numpart);
00187     if (!isoutside) {
00188       zinc_(Znotmax);  /* last point of outsideset is no longer furthest. */
00189       facet->notfurthest= True;
00190       qh_partitioncoplanar(furthest, facet, &dist);
00191       return True;
00192     }
00193   }
00194   qh_buildtracing(furthest, facet);
00195   if (qh STOPpoint < 0 && qh furthest_id == -qh STOPpoint-1) {
00196     facet->notfurthest= True;
00197     return False;
00198   }
00199   qh_findhorizon(furthest, facet, &goodvisible, &goodhorizon);
00200   if (qh ONLYgood && !(goodvisible+goodhorizon) && !qh GOODclosest) {
00201     zinc_(Znotgood);
00202     facet->notfurthest= True;
00203     /* last point of outsideset is no longer furthest.  This is ok
00204        since all points of the outside are likely to be bad */
00205     qh_resetlists(False, qh_RESETvisible /*qh visible_list newvertex_list newfacet_list */);
00206     return True;
00207   }
00208   zzinc_(Zprocessed);
00209   firstnew= qh facet_id;
00210   vertex= qh_makenewfacets(furthest /*visible_list, attaches if !ONLYgood */);
00211   qh_makenewplanes(/* newfacet_list */);
00212   numnew= qh facet_id - firstnew;
00213   newbalance= numnew - (realT) (qh num_facets-qh num_visible)
00214                          * qh hull_dim/qh num_vertices;
00215   wadd_(Wnewbalance, newbalance);
00216   wadd_(Wnewbalance2, newbalance * newbalance);
00217   if (qh ONLYgood
00218   && !qh_findgood(qh newfacet_list, goodhorizon) && !qh GOODclosest) {
00219     FORALLnew_facets
00220       qh_delfacet(newfacet);
00221     qh_delvertex(vertex);
00222     qh_resetlists(True, qh_RESETvisible /*qh visible_list newvertex_list newfacet_list */);
00223     zinc_(Znotgoodnew);
00224     facet->notfurthest= True;
00225     return True;
00226   }
00227   if (qh ONLYgood)
00228     qh_attachnewfacets(/*visible_list*/);
00229   qh_matchnewfacets();
00230   qh_updatevertices();
00231   if (qh STOPcone && qh furthest_id == qh STOPcone-1) {
00232     facet->notfurthest= True;
00233     return False;  /* visible_list etc. still defined */
00234   }
00235   qh findbestnew= False;
00236   if (qh PREmerge || qh MERGEexact) {
00237     qh_premerge(vertex, qh premerge_centrum, qh premerge_cos);
00238     if (qh_USEfindbestnew)
00239       qh findbestnew= True;
00240     else {
00241       FORALLnew_facets {
00242         if (!newfacet->simplicial) {
00243           qh findbestnew= True;  /* use qh_findbestnew instead of qh_findbest*/
00244           break;
00245         }
00246       }
00247     }
00248   }else if (qh BESToutside)
00249     qh findbestnew= True;
00250   qh_partitionvisible(/*visible_list, newfacet_list*/ !qh_ALL, &numpoints);
00251   qh findbestnew= False;
00252   qh findbest_notsharp= False;
00253   zinc_(Zpbalance);
00254   pbalance= numpoints - (realT) qh hull_dim /* assumes all points extreme */
00255                 * (qh num_points - qh num_vertices)/qh num_vertices;
00256   wadd_(Wpbalance, pbalance);
00257   wadd_(Wpbalance2, pbalance * pbalance);
00258   qh_deletevisible(/*qh visible_list*/);
00259   zmax_(Zmaxvertex, qh num_vertices);
00260   qh NEWfacets= False;
00261   if (qh IStracing >= 4) {
00262     if (qh num_facets < 2000)
00263       qh_printlists();
00264     qh_printfacetlist(qh newfacet_list, NULL, True);
00265     qh_checkpolygon(qh facet_list);
00266   }else if (qh CHECKfrequently) {
00267     if (qh num_facets < 50)
00268       qh_checkpolygon(qh facet_list);
00269     else
00270       qh_checkpolygon(qh newfacet_list);
00271   }
00272   if (qh STOPpoint > 0 && qh furthest_id == qh STOPpoint-1)
00273     return False;
00274   qh_resetlists(True, qh_RESETvisible /*qh visible_list newvertex_list newfacet_list */);
00275   /* qh_triangulate(); to test qh.TRInormals */
00276   trace2((qh ferr, 2056, "qh_addpoint: added p%d new facets %d new balance %2.2g point balance %2.2g\n",
00277     qh_pointid(furthest), numnew, newbalance, pbalance));
00278   return True;
00279 } /* addpoint */
00280 
00281 /*-<a                             href="qh-qhull.htm#TOC"
00282   >-------------------------------</a><a name="build_withrestart">-</a>
00283 
00284   qh_build_withrestart()
00285     allow restarts due to qh.JOGGLEmax while calling qh_buildhull()
00286     qh.FIRSTpoint/qh.NUMpoints is point array
00287         it may be moved by qh_joggleinput()
00288 */
00289 void qh_build_withrestart(void) {
00290   int restart;
00291 
00292   qh ALLOWrestart= True;
00293   while (True) {
00294     restart= setjmp(qh restartexit); /* simple statement for CRAY J916 */
00295     if (restart) {       /* only from qh_precision() */
00296       zzinc_(Zretry);
00297       wmax_(Wretrymax, qh JOGGLEmax);
00298       /* QH7078 warns about using 'TCn' with 'QJn' */
00299       qh STOPcone= -1; /* if break from joggle, prevents normal output */
00300     }
00301     if (!qh RERUN && qh JOGGLEmax < REALmax/2) {
00302       if (qh build_cnt > qh_JOGGLEmaxretry) {
00303         qh_fprintf(qh ferr, 6229, "qhull precision error: %d attempts to construct a convex hull\n\
00304         with joggled input.  Increase joggle above 'QJ%2.2g'\n\
00305         or modify qh_JOGGLE... parameters in user.h\n",
00306            qh build_cnt, qh JOGGLEmax);
00307         qh_errexit(qh_ERRqhull, NULL, NULL);
00308       }
00309       if (qh build_cnt && !restart)
00310         break;
00311     }else if (qh build_cnt && qh build_cnt >= qh RERUN)
00312       break;
00313     qh STOPcone= 0;
00314     qh_freebuild(True);  /* first call is a nop */
00315     qh build_cnt++;
00316     if (!qh qhull_optionsiz)
00317       qh qhull_optionsiz= (int)strlen(qh qhull_options);   /* WARN64 */
00318     else {
00319       qh qhull_options [qh qhull_optionsiz]= '\0';
00320       qh qhull_optionlen= qh_OPTIONline;  /* starts a new line */
00321     }
00322     qh_option("_run", &qh build_cnt, NULL);
00323     if (qh build_cnt == qh RERUN) {
00324       qh IStracing= qh TRACElastrun;  /* duplicated from qh_initqhull_globals */
00325       if (qh TRACEpoint != -1 || qh TRACEdist < REALmax/2 || qh TRACEmerge) {
00326         qh TRACElevel= (qh IStracing? qh IStracing : 3);
00327         qh IStracing= 0;
00328       }
00329       qhmem.IStracing= qh IStracing;
00330     }
00331     if (qh JOGGLEmax < REALmax/2)
00332       qh_joggleinput();
00333     qh_initbuild();
00334     qh_buildhull();
00335     if (qh JOGGLEmax < REALmax/2 && !qh MERGING)
00336       qh_checkconvex(qh facet_list, qh_ALGORITHMfault);
00337   }
00338   qh ALLOWrestart= False;
00339 } /* qh_build_withrestart */
00340 
00341 /*-<a                             href="qh-qhull.htm#TOC"
00342   >-------------------------------</a><a name="buildhull">-</a>
00343 
00344   qh_buildhull()
00345     construct a convex hull by adding outside points one at a time
00346 
00347   returns:
00348 
00349   notes:
00350     may be called multiple times
00351     checks facet and vertex lists for incorrect flags
00352     to recover from STOPcone, call qh_deletevisible and qh_resetlists
00353 
00354   design:
00355     check visible facet and newfacet flags
00356     check newlist vertex flags and qh.STOPcone/STOPpoint
00357     for each facet with a furthest outside point
00358       add point to facet
00359       exit if qh.STOPcone or qh.STOPpoint requested
00360     if qh.NARROWhull for initial simplex
00361       partition remaining outside points to coplanar sets
00362 */
00363 void qh_buildhull(void) {
00364   facetT *facet;
00365   pointT *furthest;
00366   vertexT *vertex;
00367   int id;
00368 
00369   trace1((qh ferr, 1037, "qh_buildhull: start build hull\n"));
00370   FORALLfacets {
00371     if (facet->visible || facet->newfacet) {
00372       qh_fprintf(qh ferr, 6165, "qhull internal error (qh_buildhull): visible or new facet f%d in facet list\n",
00373                    facet->id);
00374       qh_errexit(qh_ERRqhull, facet, NULL);
00375     }
00376   }
00377   FORALLvertices {
00378     if (vertex->newlist) {
00379       qh_fprintf(qh ferr, 6166, "qhull internal error (qh_buildhull): new vertex f%d in vertex list\n",
00380                    vertex->id);
00381       qh_errprint("ERRONEOUS", NULL, NULL, NULL, vertex);
00382       qh_errexit(qh_ERRqhull, NULL, NULL);
00383     }
00384     id= qh_pointid(vertex->point);
00385     if ((qh STOPpoint>0 && id == qh STOPpoint-1) ||
00386         (qh STOPpoint<0 && id == -qh STOPpoint-1) ||
00387         (qh STOPcone>0 && id == qh STOPcone-1)) {
00388       trace1((qh ferr, 1038,"qh_buildhull: stop point or cone P%d in initial hull\n", id));
00389       return;
00390     }
00391   }
00392   qh facet_next= qh facet_list;      /* advance facet when processed */
00393   while ((furthest= qh_nextfurthest(&facet))) {
00394     qh num_outside--;  /* if ONLYmax, furthest may not be outside */
00395     if (!qh_addpoint(furthest, facet, qh ONLYmax))
00396       break;
00397   }
00398   if (qh NARROWhull) /* move points from outsideset to coplanarset */
00399     qh_outcoplanar( /* facet_list */ );
00400   if (qh num_outside && !furthest) {
00401     qh_fprintf(qh ferr, 6167, "qhull internal error (qh_buildhull): %d outside points were never processed.\n", qh num_outside);
00402     qh_errexit(qh_ERRqhull, NULL, NULL);
00403   }
00404   trace1((qh ferr, 1039, "qh_buildhull: completed the hull construction\n"));
00405 } /* buildhull */
00406 
00407 
00408 /*-<a                             href="qh-qhull.htm#TOC"
00409   >-------------------------------</a><a name="buildtracing">-</a>
00410 
00411   qh_buildtracing( furthest, facet )
00412     trace an iteration of qh_buildhull() for furthest point and facet
00413     if !furthest, prints progress message
00414 
00415   returns:
00416     tracks progress with qh.lastreport
00417     updates qh.furthest_id (-3 if furthest is NULL)
00418     also resets visit_id, vertext_visit on wrap around
00419 
00420   see:
00421     qh_tracemerging()
00422 
00423   design:
00424     if !furthest
00425       print progress message
00426       exit
00427     if 'TFn' iteration
00428       print progress message
00429     else if tracing
00430       trace furthest point and facet
00431     reset qh.visit_id and qh.vertex_visit if overflow may occur
00432     set qh.furthest_id for tracing
00433 */
00434 void qh_buildtracing(pointT *furthest, facetT *facet) {
00435   realT dist= 0;
00436   float cpu;
00437   int total, furthestid;
00438   time_t timedata;
00439   struct tm *tp;
00440   vertexT *vertex;
00441 
00442   qh old_randomdist= qh RANDOMdist;
00443   qh RANDOMdist= False;
00444   if (!furthest) {
00445     time(&timedata);
00446     tp= localtime(&timedata);
00447     cpu= (float)qh_CPUclock - (float)qh hulltime;
00448     cpu /= (float)qh_SECticks;
00449     total= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
00450     qh_fprintf(qh ferr, 8118, "\n\
00451 At %02d:%02d:%02d & %2.5g CPU secs, qhull has created %d facets and merged %d.\n\
00452  The current hull contains %d facets and %d vertices.  Last point was p%d\n",
00453       tp->tm_hour, tp->tm_min, tp->tm_sec, cpu, qh facet_id -1,
00454       total, qh num_facets, qh num_vertices, qh furthest_id);
00455     return;
00456   }
00457   furthestid= qh_pointid(furthest);
00458   if (qh TRACEpoint == furthestid) {
00459     qh IStracing= qh TRACElevel;
00460     qhmem.IStracing= qh TRACElevel;
00461   }else if (qh TRACEpoint != -1 && qh TRACEdist < REALmax/2) {
00462     qh IStracing= 0;
00463     qhmem.IStracing= 0;
00464   }
00465   if (qh REPORTfreq && (qh facet_id-1 > qh lastreport+qh REPORTfreq)) {
00466     qh lastreport= qh facet_id-1;
00467     time(&timedata);
00468     tp= localtime(&timedata);
00469     cpu= (float)qh_CPUclock - (float)qh hulltime;
00470     cpu /= (float)qh_SECticks;
00471     total= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
00472     zinc_(Zdistio);
00473     qh_distplane(furthest, facet, &dist);
00474     qh_fprintf(qh ferr, 8119, "\n\
00475 At %02d:%02d:%02d & %2.5g CPU secs, qhull has created %d facets and merged %d.\n\
00476  The current hull contains %d facets and %d vertices.  There are %d\n\
00477  outside points.  Next is point p%d(v%d), %2.2g above f%d.\n",
00478       tp->tm_hour, tp->tm_min, tp->tm_sec, cpu, qh facet_id -1,
00479       total, qh num_facets, qh num_vertices, qh num_outside+1,
00480       furthestid, qh vertex_id, dist, getid_(facet));
00481   }else if (qh IStracing >=1) {
00482     cpu= (float)qh_CPUclock - (float)qh hulltime;
00483     cpu /= (float)qh_SECticks;
00484     qh_distplane(furthest, facet, &dist);
00485     qh_fprintf(qh ferr, 8120, "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",
00486       furthestid, qh vertex_id, qh num_facets, dist,
00487       getid_(facet), qh num_outside+1, cpu, qh furthest_id);
00488   }
00489   zmax_(Zvisit2max, (int)qh visit_id/2);
00490   if (qh visit_id > (unsigned) INT_MAX) {
00491     zinc_(Zvisit);
00492     qh visit_id= 0;
00493     FORALLfacets
00494       facet->visitid= 0;
00495   }
00496   zmax_(Zvvisit2max, (int)qh vertex_visit/2);
00497   if (qh vertex_visit > (unsigned) INT_MAX/2) { /* 31 bits */
00498     zinc_(Zvvisit);
00499     qh vertex_visit= 0;
00500     FORALLvertices
00501       vertex->visitid= 0;
00502   }
00503   qh furthest_id= furthestid;
00504   qh RANDOMdist= qh old_randomdist;
00505 } /* buildtracing */
00506 
00507 /*-<a                             href="qh-qhull.htm#TOC"
00508   >-------------------------------</a><a name="errexit2">-</a>
00509 
00510   qh_errexit2( exitcode, facet, otherfacet )
00511     return exitcode to system after an error
00512     report two facets
00513 
00514   returns:
00515     assumes exitcode non-zero
00516 
00517   see:
00518     normally use qh_errexit() in user.c(reports a facet and a ridge)
00519 */
00520 void qh_errexit2(int exitcode, facetT *facet, facetT *otherfacet) {
00521 
00522   qh_errprint("ERRONEOUS", facet, otherfacet, NULL, NULL);
00523   qh_errexit(exitcode, NULL, NULL);
00524 } /* errexit2 */
00525 
00526 
00527 /*-<a                             href="qh-qhull.htm#TOC"
00528   >-------------------------------</a><a name="findhorizon">-</a>
00529 
00530   qh_findhorizon( point, facet, goodvisible, goodhorizon )
00531     given a visible facet, find the point's horizon and visible facets
00532     for all facets, !facet-visible
00533 
00534   returns:
00535     returns qh.visible_list/num_visible with all visible facets
00536       marks visible facets with ->visible
00537     updates count of good visible and good horizon facets
00538     updates qh.max_outside, qh.max_vertex, facet->maxoutside
00539 
00540   see:
00541     similar to qh_delpoint()
00542 
00543   design:
00544     move facet to qh.visible_list at end of qh.facet_list
00545     for all visible facets
00546      for each unvisited neighbor of a visible facet
00547        compute distance of point to neighbor
00548        if point above neighbor
00549          move neighbor to end of qh.visible_list
00550        else if point is coplanar with neighbor
00551          update qh.max_outside, qh.max_vertex, neighbor->maxoutside
00552          mark neighbor coplanar (will create a samecycle later)
00553          update horizon statistics
00554 */
00555 void qh_findhorizon(pointT *point, facetT *facet, int *goodvisible, int *goodhorizon) {
00556   facetT *neighbor, **neighborp, *visible;
00557   int numhorizon= 0, coplanar= 0;
00558   realT dist;
00559 
00560   trace1((qh ferr, 1040,"qh_findhorizon: find horizon for point p%d facet f%d\n",qh_pointid(point),facet->id));
00561   *goodvisible= *goodhorizon= 0;
00562   zinc_(Ztotvisible);
00563   qh_removefacet(facet);  /* visible_list at end of qh facet_list */
00564   qh_appendfacet(facet);
00565   qh num_visible= 1;
00566   if (facet->good)
00567     (*goodvisible)++;
00568   qh visible_list= facet;
00569   facet->visible= True;
00570   facet->f.replace= NULL;
00571   if (qh IStracing >=4)
00572     qh_errprint("visible", facet, NULL, NULL, NULL);
00573   qh visit_id++;
00574   FORALLvisible_facets {
00575     if (visible->tricoplanar && !qh TRInormals) {
00576       qh_fprintf(qh ferr, 6230, "Qhull internal error (qh_findhorizon): does not work for tricoplanar facets.  Use option 'Q11'\n");
00577       qh_errexit(qh_ERRqhull, visible, NULL);
00578     }
00579     visible->visitid= qh visit_id;
00580     FOREACHneighbor_(visible) {
00581       if (neighbor->visitid == qh visit_id)
00582         continue;
00583       neighbor->visitid= qh visit_id;
00584       zzinc_(Znumvisibility);
00585       qh_distplane(point, neighbor, &dist);
00586       if (dist > qh MINvisible) {
00587         zinc_(Ztotvisible);
00588         qh_removefacet(neighbor);  /* append to end of qh visible_list */
00589         qh_appendfacet(neighbor);
00590         neighbor->visible= True;
00591         neighbor->f.replace= NULL;
00592         qh num_visible++;
00593         if (neighbor->good)
00594           (*goodvisible)++;
00595         if (qh IStracing >=4)
00596           qh_errprint("visible", neighbor, NULL, NULL, NULL);
00597       }else {
00598         if (dist > - qh MAXcoplanar) {
00599           neighbor->coplanar= True;
00600           zzinc_(Zcoplanarhorizon);
00601           qh_precision("coplanar horizon");
00602           coplanar++;
00603           if (qh MERGING) {
00604             if (dist > 0) {
00605               maximize_(qh max_outside, dist);
00606               maximize_(qh max_vertex, dist);
00607 #if qh_MAXoutside
00608               maximize_(neighbor->maxoutside, dist);
00609 #endif
00610             }else
00611               minimize_(qh min_vertex, dist);  /* due to merge later */
00612           }
00613           trace2((qh ferr, 2057, "qh_findhorizon: point p%d is coplanar to horizon f%d, dist=%2.7g < qh MINvisible(%2.7g)\n",
00614               qh_pointid(point), neighbor->id, dist, qh MINvisible));
00615         }else
00616           neighbor->coplanar= False;
00617         zinc_(Ztothorizon);
00618         numhorizon++;
00619         if (neighbor->good)
00620           (*goodhorizon)++;
00621         if (qh IStracing >=4)
00622           qh_errprint("horizon", neighbor, NULL, NULL, NULL);
00623       }
00624     }
00625   }
00626   if (!numhorizon) {
00627     qh_precision("empty horizon");
00628     qh_fprintf(qh ferr, 6168, "qhull precision error (qh_findhorizon): empty horizon\n\
00629 QhullPoint p%d was above all facets.\n", qh_pointid(point));
00630     qh_printfacetlist(qh facet_list, NULL, True);
00631     qh_errexit(qh_ERRprec, NULL, NULL);
00632   }
00633   trace1((qh ferr, 1041, "qh_findhorizon: %d horizon facets(good %d), %d visible(good %d), %d coplanar\n",
00634        numhorizon, *goodhorizon, qh num_visible, *goodvisible, coplanar));
00635   if (qh IStracing >= 4 && qh num_facets < 50)
00636     qh_printlists();
00637 } /* findhorizon */
00638 
00639 /*-<a                             href="qh-qhull.htm#TOC"
00640   >-------------------------------</a><a name="nextfurthest">-</a>
00641 
00642   qh_nextfurthest( visible )
00643     returns next furthest point and visible facet for qh_addpoint()
00644     starts search at qh.facet_next
00645 
00646   returns:
00647     removes furthest point from outside set
00648     NULL if none available
00649     advances qh.facet_next over facets with empty outside sets
00650 
00651   design:
00652     for each facet from qh.facet_next
00653       if empty outside set
00654         advance qh.facet_next
00655       else if qh.NARROWhull
00656         determine furthest outside point
00657         if furthest point is not outside
00658           advance qh.facet_next(point will be coplanar)
00659     remove furthest point from outside set
00660 */
00661 pointT *qh_nextfurthest(facetT **visible) {
00662   facetT *facet;
00663   int size, idx;
00664   realT randr, dist;
00665   pointT *furthest;
00666 
00667   while ((facet= qh facet_next) != qh facet_tail) {
00668     if (!facet->outsideset) {
00669       qh facet_next= facet->next;
00670       continue;
00671     }
00672     SETreturnsize_(facet->outsideset, size);
00673     if (!size) {
00674       qh_setfree(&facet->outsideset);
00675       qh facet_next= facet->next;
00676       continue;
00677     }
00678     if (qh NARROWhull) {
00679       if (facet->notfurthest)
00680         qh_furthestout(facet);
00681       furthest= (pointT*)qh_setlast(facet->outsideset);
00682 #if qh_COMPUTEfurthest
00683       qh_distplane(furthest, facet, &dist);
00684       zinc_(Zcomputefurthest);
00685 #else
00686       dist= facet->furthestdist;
00687 #endif
00688       if (dist < qh MINoutside) { /* remainder of outside set is coplanar for qh_outcoplanar */
00689         qh facet_next= facet->next;
00690         continue;
00691       }
00692     }
00693     if (!qh RANDOMoutside && !qh VIRTUALmemory) {
00694       if (qh PICKfurthest) {
00695         qh_furthestnext(/* qh facet_list */);
00696         facet= qh facet_next;
00697       }
00698       *visible= facet;
00699       return((pointT*)qh_setdellast(facet->outsideset));
00700     }
00701     if (qh RANDOMoutside) {
00702       int outcoplanar = 0;
00703       if (qh NARROWhull) {
00704         FORALLfacets {
00705           if (facet == qh facet_next)
00706             break;
00707           if (facet->outsideset)
00708             outcoplanar += qh_setsize( facet->outsideset);
00709         }
00710       }
00711       randr= qh_RANDOMint;
00712       randr= randr/(qh_RANDOMmax+1);
00713       idx= (int)floor((qh num_outside - outcoplanar) * randr);
00714       FORALLfacet_(qh facet_next) {
00715         if (facet->outsideset) {
00716           SETreturnsize_(facet->outsideset, size);
00717           if (!size)
00718             qh_setfree(&facet->outsideset);
00719           else if (size > idx) {
00720             *visible= facet;
00721             return((pointT*)qh_setdelnth(facet->outsideset, idx));
00722           }else
00723             idx -= size;
00724         }
00725       }
00726       qh_fprintf(qh ferr, 6169, "qhull internal error (qh_nextfurthest): num_outside %d is too low\nby at least %d, or a random real %g >= 1.0\n",
00727               qh num_outside, idx+1, randr);
00728       qh_errexit(qh_ERRqhull, NULL, NULL);
00729     }else { /* VIRTUALmemory */
00730       facet= qh facet_tail->previous;
00731       if (!(furthest= (pointT*)qh_setdellast(facet->outsideset))) {
00732         if (facet->outsideset)
00733           qh_setfree(&facet->outsideset);
00734         qh_removefacet(facet);
00735         qh_prependfacet(facet, &qh facet_list);
00736         continue;
00737       }
00738       *visible= facet;
00739       return furthest;
00740     }
00741   }
00742   return NULL;
00743 } /* nextfurthest */
00744 
00745 /*-<a                             href="qh-qhull.htm#TOC"
00746   >-------------------------------</a><a name="partitionall">-</a>
00747 
00748   qh_partitionall( vertices, points, numpoints )
00749     partitions all points in points/numpoints to the outsidesets of facets
00750     vertices= vertices in qh.facet_list(!partitioned)
00751 
00752   returns:
00753     builds facet->outsideset
00754     does not partition qh.GOODpoint
00755     if qh.ONLYgood && !qh.MERGING,
00756       does not partition qh.GOODvertex
00757 
00758   notes:
00759     faster if qh.facet_list sorted by anticipated size of outside set
00760 
00761   design:
00762     initialize pointset with all points
00763     remove vertices from pointset
00764     remove qh.GOODpointp from pointset (unless it's qh.STOPcone or qh.STOPpoint)
00765     for all facets
00766       for all remaining points in pointset
00767         compute distance from point to facet
00768         if point is outside facet
00769           remove point from pointset (by not reappending)
00770           update bestpoint
00771           append point or old bestpoint to facet's outside set
00772       append bestpoint to facet's outside set (furthest)
00773     for all points remaining in pointset
00774       partition point into facets' outside sets and coplanar sets
00775 */
00776 void qh_partitionall(setT *vertices, pointT *points, int numpoints){
00777   setT *pointset;
00778   vertexT *vertex, **vertexp;
00779   pointT *point, **pointp, *bestpoint;
00780   int size, point_i, point_n, point_end, remaining, i, id;
00781   facetT *facet;
00782   realT bestdist= -REALmax, dist, distoutside;
00783 
00784   trace1((qh ferr, 1042, "qh_partitionall: partition all points into outside sets\n"));
00785   pointset= qh_settemp(numpoints);
00786   qh num_outside= 0;
00787   pointp= SETaddr_(pointset, pointT);
00788   for (i=numpoints, point= points; i--; point += qh hull_dim)
00789     *(pointp++)= point;
00790   qh_settruncate(pointset, numpoints);
00791   FOREACHvertex_(vertices) {
00792     if ((id= qh_pointid(vertex->point)) >= 0)
00793       SETelem_(pointset, id)= NULL;
00794   }
00795   id= qh_pointid(qh GOODpointp);
00796   if (id >=0 && qh STOPcone-1 != id && -qh STOPpoint-1 != id)
00797     SETelem_(pointset, id)= NULL;
00798   if (qh GOODvertexp && qh ONLYgood && !qh MERGING) { /* matches qhull()*/
00799     if ((id= qh_pointid(qh GOODvertexp)) >= 0)
00800       SETelem_(pointset, id)= NULL;
00801   }
00802   if (!qh BESToutside) {  /* matches conditional for qh_partitionpoint below */
00803     distoutside= qh_DISToutside; /* multiple of qh.MINoutside & qh.max_outside, see user.h */
00804     zval_(Ztotpartition)= qh num_points - qh hull_dim - 1; /*misses GOOD... */
00805     remaining= qh num_facets;
00806     point_end= numpoints;
00807     FORALLfacets {
00808       size= point_end/(remaining--) + 100;
00809       facet->outsideset= qh_setnew(size);
00810       bestpoint= NULL;
00811       point_end= 0;
00812       FOREACHpoint_i_(pointset) {
00813         if (point) {
00814           zzinc_(Zpartitionall);
00815           qh_distplane(point, facet, &dist);
00816           if (dist < distoutside)
00817             SETelem_(pointset, point_end++)= point;
00818           else {
00819             qh num_outside++;
00820             if (!bestpoint) {
00821               bestpoint= point;
00822               bestdist= dist;
00823             }else if (dist > bestdist) {
00824               qh_setappend(&facet->outsideset, bestpoint);
00825               bestpoint= point;
00826               bestdist= dist;
00827             }else
00828               qh_setappend(&facet->outsideset, point);
00829           }
00830         }
00831       }
00832       if (bestpoint) {
00833         qh_setappend(&facet->outsideset, bestpoint);
00834 #if !qh_COMPUTEfurthest
00835         facet->furthestdist= bestdist;
00836 #endif
00837       }else
00838         qh_setfree(&facet->outsideset);
00839       qh_settruncate(pointset, point_end);
00840     }
00841   }
00842   /* if !qh BESToutside, pointset contains points not assigned to outsideset */
00843   if (qh BESToutside || qh MERGING || qh KEEPcoplanar || qh KEEPinside) {
00844     qh findbestnew= True;
00845     FOREACHpoint_i_(pointset) {
00846       if (point)
00847         qh_partitionpoint(point, qh facet_list);
00848     }
00849     qh findbestnew= False;
00850   }
00851   zzadd_(Zpartitionall, zzval_(Zpartition));
00852   zzval_(Zpartition)= 0;
00853   qh_settempfree(&pointset);
00854   if (qh IStracing >= 4)
00855     qh_printfacetlist(qh facet_list, NULL, True);
00856 } /* partitionall */
00857 
00858 
00859 /*-<a                             href="qh-qhull.htm#TOC"
00860   >-------------------------------</a><a name="partitioncoplanar">-</a>
00861 
00862   qh_partitioncoplanar( point, facet, dist )
00863     partition coplanar point to a facet
00864     dist is distance from point to facet
00865     if dist NULL,
00866       searches for bestfacet and does nothing if inside
00867     if qh.findbestnew set,
00868       searches new facets instead of using qh_findbest()
00869 
00870   returns:
00871     qh.max_ouside updated
00872     if qh.KEEPcoplanar or qh.KEEPinside
00873       point assigned to best coplanarset
00874 
00875   notes:
00876     facet->maxoutside is updated at end by qh_check_maxout
00877 
00878   design:
00879     if dist undefined
00880       find best facet for point
00881       if point sufficiently below facet (depends on qh.NEARinside and qh.KEEPinside)
00882         exit
00883     if keeping coplanar/nearinside/inside points
00884       if point is above furthest coplanar point
00885         append point to coplanar set (it is the new furthest)
00886         update qh.max_outside
00887       else
00888         append point one before end of coplanar set
00889     else if point is clearly outside of qh.max_outside and bestfacet->coplanarset
00890     and bestfacet is more than perpendicular to facet
00891       repartition the point using qh_findbest() -- it may be put on an outsideset
00892     else
00893       update qh.max_outside
00894 */
00895 void qh_partitioncoplanar(pointT *point, facetT *facet, realT *dist) {
00896   facetT *bestfacet;
00897   pointT *oldfurthest;
00898   realT bestdist, dist2= 0, angle;
00899   int numpart= 0, oldfindbest;
00900   boolT isoutside;
00901 
00902   qh WAScoplanar= True;
00903   if (!dist) {
00904     if (qh findbestnew)
00905       bestfacet= qh_findbestnew(point, facet, &bestdist, qh_ALL, &isoutside, &numpart);
00906     else
00907       bestfacet= qh_findbest(point, facet, qh_ALL, !qh_ISnewfacets, qh DELAUNAY,
00908                           &bestdist, &isoutside, &numpart);
00909     zinc_(Ztotpartcoplanar);
00910     zzadd_(Zpartcoplanar, numpart);
00911     if (!qh DELAUNAY && !qh KEEPinside) { /*  for 'd', bestdist skips upperDelaunay facets */
00912       if (qh KEEPnearinside) {
00913         if (bestdist < -qh NEARinside) {
00914           zinc_(Zcoplanarinside);
00915           trace4((qh ferr, 4062, "qh_partitioncoplanar: point p%d is more than near-inside facet f%d dist %2.2g findbestnew %d\n",
00916                   qh_pointid(point), bestfacet->id, bestdist, qh findbestnew));
00917           return;
00918         }
00919       }else if (bestdist < -qh MAXcoplanar) {
00920           trace4((qh ferr, 4063, "qh_partitioncoplanar: point p%d is inside facet f%d dist %2.2g findbestnew %d\n",
00921                   qh_pointid(point), bestfacet->id, bestdist, qh findbestnew));
00922         zinc_(Zcoplanarinside);
00923         return;
00924       }
00925     }
00926   }else {
00927     bestfacet= facet;
00928     bestdist= *dist;
00929   }
00930   if (bestdist > qh max_outside) {
00931     if (!dist && facet != bestfacet) {
00932       zinc_(Zpartangle);
00933       angle= qh_getangle(facet->normal, bestfacet->normal);
00934       if (angle < 0) {
00935         /* typically due to deleted vertex and coplanar facets, e.g.,
00936              RBOX 1000 s Z1 G1e-13 t1001185205 | QHULL Tv */
00937         zinc_(Zpartflip);
00938         trace2((qh ferr, 2058, "qh_partitioncoplanar: repartition point p%d from f%d.  It is above flipped facet f%d dist %2.2g\n",
00939                 qh_pointid(point), facet->id, bestfacet->id, bestdist));
00940         oldfindbest= qh findbestnew;
00941         qh findbestnew= False;
00942         qh_partitionpoint(point, bestfacet);
00943         qh findbestnew= oldfindbest;
00944         return;
00945       }
00946     }
00947     qh max_outside= bestdist;
00948     if (bestdist > qh TRACEdist) {
00949       qh_fprintf(qh ferr, 8122, "qh_partitioncoplanar: ====== p%d from f%d increases max_outside to %2.2g of f%d last p%d\n",
00950                      qh_pointid(point), facet->id, bestdist, bestfacet->id, qh furthest_id);
00951       qh_errprint("DISTANT", facet, bestfacet, NULL, NULL);
00952     }
00953   }
00954   if (qh KEEPcoplanar + qh KEEPinside + qh KEEPnearinside) {
00955     oldfurthest= (pointT*)qh_setlast(bestfacet->coplanarset);
00956     if (oldfurthest) {
00957       zinc_(Zcomputefurthest);
00958       qh_distplane(oldfurthest, bestfacet, &dist2);
00959     }
00960     if (!oldfurthest || dist2 < bestdist)
00961       qh_setappend(&bestfacet->coplanarset, point);
00962     else
00963       qh_setappend2ndlast(&bestfacet->coplanarset, point);
00964   }
00965   trace4((qh ferr, 4064, "qh_partitioncoplanar: point p%d is coplanar with facet f%d(or inside) dist %2.2g\n",
00966           qh_pointid(point), bestfacet->id, bestdist));
00967 } /* partitioncoplanar */
00968 
00969 /*-<a                             href="qh-qhull.htm#TOC"
00970   >-------------------------------</a><a name="partitionpoint">-</a>
00971 
00972   qh_partitionpoint( point, facet )
00973     assigns point to an outside set, coplanar set, or inside set (i.e., dropt)
00974     if qh.findbestnew
00975       uses qh_findbestnew() to search all new facets
00976     else
00977       uses qh_findbest()
00978 
00979   notes:
00980     after qh_distplane(), this and qh_findbest() are most expensive in 3-d
00981 
00982   design:
00983     find best facet for point
00984       (either exhaustive search of new facets or directed search from facet)
00985     if qh.NARROWhull
00986       retain coplanar and nearinside points as outside points
00987     if point is outside bestfacet
00988       if point above furthest point for bestfacet
00989         append point to outside set (it becomes the new furthest)
00990         if outside set was empty
00991           move bestfacet to end of qh.facet_list (i.e., after qh.facet_next)
00992         update bestfacet->furthestdist
00993       else
00994         append point one before end of outside set
00995     else if point is coplanar to bestfacet
00996       if keeping coplanar points or need to update qh.max_outside
00997         partition coplanar point into bestfacet
00998     else if near-inside point
00999       partition as coplanar point into bestfacet
01000     else is an inside point
01001       if keeping inside points
01002         partition as coplanar point into bestfacet
01003 */
01004 void qh_partitionpoint(pointT *point, facetT *facet) {
01005   realT bestdist;
01006   boolT isoutside;
01007   facetT *bestfacet;
01008   int numpart;
01009 #if qh_COMPUTEfurthest
01010   realT dist;
01011 #endif
01012 
01013   if (qh findbestnew)
01014     bestfacet= qh_findbestnew(point, facet, &bestdist, qh BESToutside, &isoutside, &numpart);
01015   else
01016     bestfacet= qh_findbest(point, facet, qh BESToutside, qh_ISnewfacets, !qh_NOupper,
01017                           &bestdist, &isoutside, &numpart);
01018   zinc_(Ztotpartition);
01019   zzadd_(Zpartition, numpart);
01020   if (qh NARROWhull) {
01021     if (qh DELAUNAY && !isoutside && bestdist >= -qh MAXcoplanar)
01022       qh_precision("nearly incident point(narrow hull)");
01023     if (qh KEEPnearinside) {
01024       if (bestdist >= -qh NEARinside)
01025         isoutside= True;
01026     }else if (bestdist >= -qh MAXcoplanar)
01027       isoutside= True;
01028   }
01029 
01030   if (isoutside) {
01031     if (!bestfacet->outsideset
01032     || !qh_setlast(bestfacet->outsideset)) {
01033       qh_setappend(&(bestfacet->outsideset), point);
01034       if (!bestfacet->newfacet) {
01035         qh_removefacet(bestfacet);  /* make sure it's after qh facet_next */
01036         qh_appendfacet(bestfacet);
01037       }
01038 #if !qh_COMPUTEfurthest
01039       bestfacet->furthestdist= bestdist;
01040 #endif
01041     }else {
01042 #if qh_COMPUTEfurthest
01043       zinc_(Zcomputefurthest);
01044       qh_distplane(oldfurthest, bestfacet, &dist);
01045       if (dist < bestdist)
01046         qh_setappend(&(bestfacet->outsideset), point);
01047       else
01048         qh_setappend2ndlast(&(bestfacet->outsideset), point);
01049 #else
01050       if (bestfacet->furthestdist < bestdist) {
01051         qh_setappend(&(bestfacet->outsideset), point);
01052         bestfacet->furthestdist= bestdist;
01053       }else
01054         qh_setappend2ndlast(&(bestfacet->outsideset), point);
01055 #endif
01056     }
01057     qh num_outside++;
01058     trace4((qh ferr, 4065, "qh_partitionpoint: point p%d is outside facet f%d new? %d(or narrowhull)\n",
01059           qh_pointid(point), bestfacet->id, bestfacet->newfacet));
01060   }else if (qh DELAUNAY || bestdist >= -qh MAXcoplanar) { /* for 'd', bestdist skips upperDelaunay facets */
01061     zzinc_(Zcoplanarpart);
01062     if (qh DELAUNAY)
01063       qh_precision("nearly incident point");
01064     if ((qh KEEPcoplanar + qh KEEPnearinside) || bestdist > qh max_outside)
01065       qh_partitioncoplanar(point, bestfacet, &bestdist);
01066     else {
01067       trace4((qh ferr, 4066, "qh_partitionpoint: point p%d is coplanar to facet f%d(dropped)\n",
01068           qh_pointid(point), bestfacet->id));
01069     }
01070   }else if (qh KEEPnearinside && bestdist > -qh NEARinside) {
01071     zinc_(Zpartnear);
01072     qh_partitioncoplanar(point, bestfacet, &bestdist);
01073   }else {
01074     zinc_(Zpartinside);
01075     trace4((qh ferr, 4067, "qh_partitionpoint: point p%d is inside all facets, closest to f%d dist %2.2g\n",
01076           qh_pointid(point), bestfacet->id, bestdist));
01077     if (qh KEEPinside)
01078       qh_partitioncoplanar(point, bestfacet, &bestdist);
01079   }
01080 } /* partitionpoint */
01081 
01082 /*-<a                             href="qh-qhull.htm#TOC"
01083   >-------------------------------</a><a name="partitionvisible">-</a>
01084 
01085   qh_partitionvisible( allpoints, numoutside )
01086     partitions points in visible facets to qh.newfacet_list
01087     qh.visible_list= visible facets
01088     for visible facets
01089       1st neighbor (if any) points to a horizon facet or a new facet
01090     if allpoints(!used),
01091       repartitions coplanar points
01092 
01093   returns:
01094     updates outside sets and coplanar sets of qh.newfacet_list
01095     updates qh.num_outside (count of outside points)
01096 
01097   notes:
01098     qh.findbest_notsharp should be clear (extra work if set)
01099 
01100   design:
01101     for all visible facets with outside set or coplanar set
01102       select a newfacet for visible facet
01103       if outside set
01104         partition outside set into new facets
01105       if coplanar set and keeping coplanar/near-inside/inside points
01106         if allpoints
01107           partition coplanar set into new facets, may be assigned outside
01108         else
01109           partition coplanar set into coplanar sets of new facets
01110     for each deleted vertex
01111       if allpoints
01112         partition vertex into new facets, may be assigned outside
01113       else
01114         partition vertex into coplanar sets of new facets
01115 */
01116 void qh_partitionvisible(/*visible_list*/ boolT allpoints, int *numoutside) {
01117   facetT *visible, *newfacet;
01118   pointT *point, **pointp;
01119   int coplanar=0, size;
01120   unsigned count;
01121   vertexT *vertex, **vertexp;
01122 
01123   if (qh ONLYmax)
01124     maximize_(qh MINoutside, qh max_vertex);
01125   *numoutside= 0;
01126   FORALLvisible_facets {
01127     if (!visible->outsideset && !visible->coplanarset)
01128       continue;
01129     newfacet= visible->f.replace;
01130     count= 0;
01131     while (newfacet && newfacet->visible) {
01132       newfacet= newfacet->f.replace;
01133       if (count++ > qh facet_id)
01134         qh_infiniteloop(visible);
01135     }
01136     if (!newfacet)
01137       newfacet= qh newfacet_list;
01138     if (newfacet == qh facet_tail) {
01139       qh_fprintf(qh ferr, 6170, "qhull precision error (qh_partitionvisible): all new facets deleted as\n        degenerate facets. Can not continue.\n");
01140       qh_errexit(qh_ERRprec, NULL, NULL);
01141     }
01142     if (visible->outsideset) {
01143       size= qh_setsize(visible->outsideset);
01144       *numoutside += size;
01145       qh num_outside -= size;
01146       FOREACHpoint_(visible->outsideset)
01147         qh_partitionpoint(point, newfacet);
01148     }
01149     if (visible->coplanarset && (qh KEEPcoplanar + qh KEEPinside + qh KEEPnearinside)) {
01150       size= qh_setsize(visible->coplanarset);
01151       coplanar += size;
01152       FOREACHpoint_(visible->coplanarset) {
01153         if (allpoints) /* not used */
01154           qh_partitionpoint(point, newfacet);
01155         else
01156           qh_partitioncoplanar(point, newfacet, NULL);
01157       }
01158     }
01159   }
01160   FOREACHvertex_(qh del_vertices) {
01161     if (vertex->point) {
01162       if (allpoints) /* not used */
01163         qh_partitionpoint(vertex->point, qh newfacet_list);
01164       else
01165         qh_partitioncoplanar(vertex->point, qh newfacet_list, NULL);
01166     }
01167   }
01168   trace1((qh ferr, 1043,"qh_partitionvisible: partitioned %d points from outsidesets and %d points from coplanarsets\n", *numoutside, coplanar));
01169 } /* partitionvisible */
01170 
01171 
01172 
01173 /*-<a                             href="qh-qhull.htm#TOC"
01174   >-------------------------------</a><a name="precision">-</a>
01175 
01176   qh_precision( reason )
01177     restart on precision errors if not merging and if 'QJn'
01178 */
01179 void qh_precision(const char *reason) {
01180 
01181   if (qh ALLOWrestart && !qh PREmerge && !qh MERGEexact) {
01182     if (qh JOGGLEmax < REALmax/2) {
01183       trace0((qh ferr, 26, "qh_precision: qhull restart because of %s\n", reason));
01184       longjmp(qh restartexit, qh_ERRprec);
01185     }
01186   }
01187 } /* qh_precision */
01188 
01189 /*-<a                             href="qh-qhull.htm#TOC"
01190   >-------------------------------</a><a name="printsummary">-</a>
01191 
01192   qh_printsummary( fp )
01193     prints summary to fp
01194 
01195   notes:
01196     not in io.c so that user_eg.c can prevent io.c from loading
01197     qh_printsummary and qh_countfacets must match counts
01198 
01199   design:
01200     determine number of points, vertices, and coplanar points
01201     print summary
01202 */
01203 void qh_printsummary(FILE *fp) {
01204   realT ratio, outerplane, innerplane;
01205   float cpu;
01206   int size, id, nummerged, numvertices, numcoplanars= 0, nonsimplicial=0;
01207   int goodused;
01208   facetT *facet;
01209   const char *s;
01210   int numdel= zzval_(Zdelvertextot);
01211   int numtricoplanars= 0;
01212 
01213   size= qh num_points + qh_setsize(qh other_points);
01214   numvertices= qh num_vertices - qh_setsize(qh del_vertices);
01215   id= qh_pointid(qh GOODpointp);
01216   FORALLfacets {
01217     if (facet->coplanarset)
01218       numcoplanars += qh_setsize( facet->coplanarset);
01219     if (facet->good) {
01220       if (facet->simplicial) {
01221         if (facet->keepcentrum && facet->tricoplanar)
01222           numtricoplanars++;
01223       }else if (qh_setsize(facet->vertices) != qh hull_dim)
01224         nonsimplicial++;
01225     }
01226   }
01227   if (id >=0 && qh STOPcone-1 != id && -qh STOPpoint-1 != id)
01228     size--;
01229   if (qh STOPcone || qh STOPpoint)
01230       qh_fprintf(fp, 9288, "\nAt a premature exit due to 'TVn', 'TCn', 'TRn', or precision error with 'QJn'.");
01231   if (qh UPPERdelaunay)
01232     goodused= qh GOODvertex + qh GOODpoint + qh SPLITthresholds;
01233   else if (qh DELAUNAY)
01234     goodused= qh GOODvertex + qh GOODpoint + qh GOODthreshold;
01235   else
01236     goodused= qh num_good;
01237   nummerged= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
01238   if (qh VORONOI) {
01239     if (qh UPPERdelaunay)
01240       qh_fprintf(fp, 9289, "\n\
01241 Furthest-site Voronoi vertices by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
01242     else
01243       qh_fprintf(fp, 9290, "\n\
01244 Voronoi diagram by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
01245     qh_fprintf(fp, 9291, "  Number of Voronoi regions%s: %d\n",
01246               qh ATinfinity ? " and at-infinity" : "", numvertices);
01247     if (numdel)
01248       qh_fprintf(fp, 9292, "  Total number of deleted points due to merging: %d\n", numdel);
01249     if (numcoplanars - numdel > 0)
01250       qh_fprintf(fp, 9293, "  Number of nearly incident points: %d\n", numcoplanars - numdel);
01251     else if (size - numvertices - numdel > 0)
01252       qh_fprintf(fp, 9294, "  Total number of nearly incident points: %d\n", size - numvertices - numdel);
01253     qh_fprintf(fp, 9295, "  Number of%s Voronoi vertices: %d\n",
01254               goodused ? " 'good'" : "", qh num_good);
01255     if (nonsimplicial)
01256       qh_fprintf(fp, 9296, "  Number of%s non-simplicial Voronoi vertices: %d\n",
01257               goodused ? " 'good'" : "", nonsimplicial);
01258   }else if (qh DELAUNAY) {
01259     if (qh UPPERdelaunay)
01260       qh_fprintf(fp, 9297, "\n\
01261 Furthest-site Delaunay triangulation by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
01262     else
01263       qh_fprintf(fp, 9298, "\n\
01264 Delaunay triangulation by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
01265     qh_fprintf(fp, 9299, "  Number of input sites%s: %d\n",
01266               qh ATinfinity ? " and at-infinity" : "", numvertices);
01267     if (numdel)
01268       qh_fprintf(fp, 9300, "  Total number of deleted points due to merging: %d\n", numdel);
01269     if (numcoplanars - numdel > 0)
01270       qh_fprintf(fp, 9301, "  Number of nearly incident points: %d\n", numcoplanars - numdel);
01271     else if (size - numvertices - numdel > 0)
01272       qh_fprintf(fp, 9302, "  Total number of nearly incident points: %d\n", size - numvertices - numdel);
01273     qh_fprintf(fp, 9303, "  Number of%s Delaunay regions: %d\n",
01274               goodused ? " 'good'" : "", qh num_good);
01275     if (nonsimplicial)
01276       qh_fprintf(fp, 9304, "  Number of%s non-simplicial Delaunay regions: %d\n",
01277               goodused ? " 'good'" : "", nonsimplicial);
01278   }else if (qh HALFspace) {
01279     qh_fprintf(fp, 9305, "\n\
01280 Halfspace intersection by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
01281     qh_fprintf(fp, 9306, "  Number of halfspaces: %d\n", size);
01282     qh_fprintf(fp, 9307, "  Number of non-redundant halfspaces: %d\n", numvertices);
01283     if (numcoplanars) {
01284       if (qh KEEPinside && qh KEEPcoplanar)
01285         s= "similar and redundant";
01286       else if (qh KEEPinside)
01287         s= "redundant";
01288       else
01289         s= "similar";
01290       qh_fprintf(fp, 9308, "  Number of %s halfspaces: %d\n", s, numcoplanars);
01291     }
01292     qh_fprintf(fp, 9309, "  Number of intersection points: %d\n", qh num_facets - qh num_visible);
01293     if (goodused)
01294       qh_fprintf(fp, 9310, "  Number of 'good' intersection points: %d\n", qh num_good);
01295     if (nonsimplicial)
01296       qh_fprintf(fp, 9311, "  Number of%s non-simplicial intersection points: %d\n",
01297               goodused ? " 'good'" : "", nonsimplicial);
01298   }else {
01299     qh_fprintf(fp, 9312, "\n\
01300 Convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
01301     qh_fprintf(fp, 9313, "  Number of vertices: %d\n", numvertices);
01302     if (numcoplanars) {
01303       if (qh KEEPinside && qh KEEPcoplanar)
01304         s= "coplanar and interior";
01305       else if (qh KEEPinside)
01306         s= "interior";
01307       else
01308         s= "coplanar";
01309       qh_fprintf(fp, 9314, "  Number of %s points: %d\n", s, numcoplanars);
01310     }
01311     qh_fprintf(fp, 9315, "  Number of facets: %d\n", qh num_facets - qh num_visible);
01312     if (goodused)
01313       qh_fprintf(fp, 9316, "  Number of 'good' facets: %d\n", qh num_good);
01314     if (nonsimplicial)
01315       qh_fprintf(fp, 9317, "  Number of%s non-simplicial facets: %d\n",
01316               goodused ? " 'good'" : "", nonsimplicial);
01317   }
01318   if (numtricoplanars)
01319       qh_fprintf(fp, 9318, "  Number of triangulated facets: %d\n", numtricoplanars);
01320   qh_fprintf(fp, 9319, "\nStatistics for: %s | %s",
01321                       qh rbox_command, qh qhull_command);
01322   if (qh ROTATErandom != INT_MIN)
01323     qh_fprintf(fp, 9320, " QR%d\n\n", qh ROTATErandom);
01324   else
01325     qh_fprintf(fp, 9321, "\n\n");
01326   qh_fprintf(fp, 9322, "  Number of points processed: %d\n", zzval_(Zprocessed));
01327   qh_fprintf(fp, 9323, "  Number of hyperplanes created: %d\n", zzval_(Zsetplane));
01328   if (qh DELAUNAY)
01329     qh_fprintf(fp, 9324, "  Number of facets in hull: %d\n", qh num_facets - qh num_visible);
01330   qh_fprintf(fp, 9325, "  Number of distance tests for qhull: %d\n", zzval_(Zpartition)+
01331       zzval_(Zpartitionall)+zzval_(Znumvisibility)+zzval_(Zpartcoplanar));
01332 #if 0  /* NOTE: must print before printstatistics() */
01333   {realT stddev, ave;
01334   qh_fprintf(fp, 9326, "  average new facet balance: %2.2g\n",
01335           wval_(Wnewbalance)/zval_(Zprocessed));
01336   stddev= qh_stddev(zval_(Zprocessed), wval_(Wnewbalance),
01337                                  wval_(Wnewbalance2), &ave);
01338   qh_fprintf(fp, 9327, "  new facet standard deviation: %2.2g\n", stddev);
01339   qh_fprintf(fp, 9328, "  average partition balance: %2.2g\n",
01340           wval_(Wpbalance)/zval_(Zpbalance));
01341   stddev= qh_stddev(zval_(Zpbalance), wval_(Wpbalance),
01342                                  wval_(Wpbalance2), &ave);
01343   qh_fprintf(fp, 9329, "  partition standard deviation: %2.2g\n", stddev);
01344   }
01345 #endif
01346   if (nummerged) {
01347     qh_fprintf(fp, 9330,"  Number of distance tests for merging: %d\n",zzval_(Zbestdist)+
01348           zzval_(Zcentrumtests)+zzval_(Zdistconvex)+zzval_(Zdistcheck)+
01349           zzval_(Zdistzero));
01350     qh_fprintf(fp, 9331,"  Number of distance tests for checking: %d\n",zzval_(Zcheckpart));
01351     qh_fprintf(fp, 9332,"  Number of merged facets: %d\n", nummerged);
01352   }
01353   if (!qh RANDOMoutside && qh QHULLfinished) {
01354     cpu= (float)qh hulltime;
01355     cpu /= (float)qh_SECticks;
01356     wval_(Wcpu)= cpu;
01357     qh_fprintf(fp, 9333, "  CPU seconds to compute hull (after input): %2.4g\n", cpu);
01358   }
01359   if (qh RERUN) {
01360     if (!qh PREmerge && !qh MERGEexact)
01361       qh_fprintf(fp, 9334, "  Percentage of runs with precision errors: %4.1f\n",
01362            zzval_(Zretry)*100.0/qh build_cnt);  /* careful of order */
01363   }else if (qh JOGGLEmax < REALmax/2) {
01364     if (zzval_(Zretry))
01365       qh_fprintf(fp, 9335, "  After %d retries, input joggled by: %2.2g\n",
01366          zzval_(Zretry), qh JOGGLEmax);
01367     else
01368       qh_fprintf(fp, 9336, "  Input joggled by: %2.2g\n", qh JOGGLEmax);
01369   }
01370   if (qh totarea != 0.0)
01371     qh_fprintf(fp, 9337, "  %s facet area:   %2.8g\n",
01372             zzval_(Ztotmerge) ? "Approximate" : "Total", qh totarea);
01373   if (qh totvol != 0.0)
01374     qh_fprintf(fp, 9338, "  %s volume:       %2.8g\n",
01375             zzval_(Ztotmerge) ? "Approximate" : "Total", qh totvol);
01376   if (qh MERGING) {
01377     qh_outerinner(NULL, &outerplane, &innerplane);
01378     if (outerplane > 2 * qh DISTround) {
01379       qh_fprintf(fp, 9339, "  Maximum distance of %spoint above facet: %2.2g",
01380             (qh QHULLfinished ? "" : "merged "), outerplane);
01381       ratio= outerplane/(qh ONEmerge + qh DISTround);
01382       /* don't report ratio if MINoutside is large */
01383       if (ratio > 0.05 && 2* qh ONEmerge > qh MINoutside && qh JOGGLEmax > REALmax/2)
01384         qh_fprintf(fp, 9340, " (%.1fx)\n", ratio);
01385       else
01386         qh_fprintf(fp, 9341, "\n");
01387     }
01388     if (innerplane < -2 * qh DISTround) {
01389       qh_fprintf(fp, 9342, "  Maximum distance of %svertex below facet: %2.2g",
01390             (qh QHULLfinished ? "" : "merged "), innerplane);
01391       ratio= -innerplane/(qh ONEmerge+qh DISTround);
01392       if (ratio > 0.05 && qh JOGGLEmax > REALmax/2)
01393         qh_fprintf(fp, 9343, " (%.1fx)\n", ratio);
01394       else
01395         qh_fprintf(fp, 9344, "\n");
01396     }
01397   }
01398   qh_fprintf(fp, 9345, "\n");
01399 } /* printsummary */
01400 
01401
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

libqhull
Author(s): Robert Krug
autogenerated on Tue Jun 18 2013 12:38:49