graspit: /opt/ros/diamondback/stacks/graspit_simulator/graspit/graspit

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