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

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