Doxygen Source Code Documentation

Main Page Alphabetical List Data Structures File List Data Fields Globals Search
qhull.c

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

Sections

Personal tools

Navigation

Doxygen Source Code Documentation

qhull.c
