graspit: poly2.c Source File

00001 /*<html><pre>  -<a                             href="qh-poly.htm"
00002   >-------------------------------</a><a name="TOP">-</a>
00003 
00004    poly2.c
00005    implements polygons and simplices
00006 
00007    see qh-poly.htm, poly.h and libqhull.h
00008 
00009    frequently used code is in poly.c
00010 
00011    Copyright (c) 1993-2012 The Geometry Center.
00012    $Id: //main/2011/qhull/src/libqhull/poly2.c#5 $$Change: 1490 $
00013    $DateTime: 2012/02/19 20:27:01 $$Author: bbarber $
00014 */
00015 
00016 #include "qhull_a.h"
00017 
00018 /*======== functions in alphabetical order ==========*/
00019 
00020 /*-<a                             href="qh-poly.htm#TOC"
00021   >-------------------------------</a><a name="addhash">-</a>
00022 
00023   qh_addhash( newelem, hashtable, hashsize, hash )
00024     add newelem to linear hash table at hash if not already there
00025 */
00026 void qh_addhash(void* newelem, setT *hashtable, int hashsize, int hash) {
00027   int scan;
00028   void *elem;
00029 
00030   for (scan= (int)hash; (elem= SETelem_(hashtable, scan));
00031        scan= (++scan >= hashsize ? 0 : scan)) {
00032     if (elem == newelem)
00033       break;
00034   }
00035   /* loop terminates because qh_HASHfactor >= 1.1 by qh_initbuffers */
00036   if (!elem)
00037     SETelem_(hashtable, scan)= newelem;
00038 } /* addhash */
00039 
00040 /*-<a                             href="qh-poly.htm#TOC"
00041   >-------------------------------</a><a name="check_bestdist">-</a>
00042 
00043   qh_check_bestdist()
00044     check that all points are within max_outside of the nearest facet
00045     if qh.ONLYgood,
00046       ignores !good facets
00047 
00048   see:
00049     qh_check_maxout(), qh_outerinner()
00050 
00051   notes:
00052     only called from qh_check_points()
00053       seldom used since qh.MERGING is almost always set
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, unassigned;
00068   facetT *facet, *bestfacet, *errfacet1= NULL, *errfacet2= NULL;
00069   facetT *facetlist;
00070   realT dist, maxoutside, maxdist= -REALmax;
00071   pointT *point;
00072   int numpart= 0, facet_i, facet_n, notgood= 0, notverified= 0;
00073   setT *facets;
00074 
00075   trace1((qh ferr, 1020, "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, 1021, "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     qh_fprintf(qh ferr, 8091, "\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     qh_distplane(point, facet, &dist);
00098     numpart++;
00099     bestfacet= qh_findbesthorizon(!qh_IScheckmax, point, facet, qh_NOupper, &dist, &numpart);
00100     /* occurs after statistics reported */
00101     maximize_(maxdist, dist);
00102     if (dist > maxoutside) {
00103       if (qh ONLYgood && !bestfacet->good
00104           && !((bestfacet= qh_findgooddist(point, bestfacet, &dist, &facetlist))
00105                && dist > maxoutside))
00106         notgood++;
00107       else {
00108         waserror= True;
00109         qh_fprintf(qh ferr, 6109, "qhull precision error: point p%d is outside facet f%d, distance= %6.8g maxoutside= %6.8g\n",
00110                 facet_i, bestfacet->id, dist, maxoutside);
00111         if (errfacet1 != bestfacet) {
00112           errfacet2= errfacet1;
00113           errfacet1= bestfacet;
00114         }
00115       }
00116     }else if (unassigned && dist < -qh MAXcoplanar)
00117       notverified++;
00118   }
00119   qh_settempfree(&facets);
00120   if (notverified && !qh DELAUNAY && !qh_QUICKhelp && qh PRINTprecision)
00121     qh_fprintf(qh ferr, 8092, "\n%d points were well inside the hull.  If the hull contains\n\
00122 a lens-shaped component, these points were not verified.  Use\n\
00123 options 'Qci Tv' to verify all points.\n", notverified);
00124   if (maxdist > qh outside_err) {
00125     qh_fprintf(qh ferr, 6110, "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",
00126               maxdist, qh outside_err);
00127     qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);
00128   }else if (waserror && qh outside_err > REALmax/2)
00129     qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);
00130   /* else if waserror, the error was logged to qh.ferr but does not effect the output */
00131   trace0((qh ferr, 20, "qh_check_bestdist: max distance outside %2.2g\n", maxdist));
00132 } /* check_bestdist */
00133 
00134 /*-<a                             href="qh-poly.htm#TOC"
00135   >-------------------------------</a><a name="check_maxout">-</a>
00136 
00137   qh_check_maxout()
00138     updates qh.max_outside by checking all points against bestfacet
00139     if qh.ONLYgood, ignores !good facets
00140 
00141   returns:
00142     updates facet->maxoutside via qh_findbesthorizon()
00143     sets qh.maxoutdone
00144     if printing qh.min_vertex (qh_outerinner),
00145       it is updated to the current vertices
00146     removes inside/coplanar points from coplanarset as needed
00147 
00148   notes:
00149     defines coplanar as min_vertex instead of MAXcoplanar
00150     may not need to check near-inside points because of qh.MAXcoplanar
00151       and qh.KEEPnearinside (before it was -DISTround)
00152 
00153   see also:
00154     qh_check_bestdist()
00155 
00156   design:
00157     if qh.min_vertex is needed
00158       for all neighbors of all vertices
00159         test distance from vertex to neighbor
00160     determine facet for each point (if any)
00161     for each point with an assigned facet
00162       find the best facet for the point and check all coplanar facets
00163         (updates outer planes)
00164     remove near-inside points from coplanar sets
00165 */
00166 #ifndef qh_NOmerge
00167 void qh_check_maxout(void) {
00168   facetT *facet, *bestfacet, *neighbor, **neighborp, *facetlist;
00169   realT dist, maxoutside, minvertex, old_maxoutside;
00170   pointT *point;
00171   int numpart= 0, facet_i, facet_n, notgood= 0;
00172   setT *facets, *vertices;
00173   vertexT *vertex;
00174 
00175   trace1((qh ferr, 1022, "qh_check_maxout: check and update maxoutside for each facet.\n"));
00176   maxoutside= minvertex= 0;
00177   if (qh VERTEXneighbors
00178   && (qh PRINTsummary || qh KEEPinside || qh KEEPcoplanar
00179         || qh TRACElevel || qh PRINTstatistics
00180         || qh PRINTout[0] == qh_PRINTsummary || qh PRINTout[0] == qh_PRINTnone)) {
00181     trace1((qh ferr, 1023, "qh_check_maxout: determine actual maxoutside and minvertex\n"));
00182     vertices= qh_pointvertex(/*qh facet_list*/);
00183     FORALLvertices {
00184       FOREACHneighbor_(vertex) {
00185         zinc_(Zdistvertex);  /* distance also computed by main loop below */
00186         qh_distplane(vertex->point, neighbor, &dist);
00187         minimize_(minvertex, dist);
00188         if (-dist > qh TRACEdist || dist > qh TRACEdist
00189         || neighbor == qh tracefacet || vertex == qh tracevertex)
00190           qh_fprintf(qh ferr, 8093, "qh_check_maxout: p%d(v%d) is %.2g from f%d\n",
00191                     qh_pointid(vertex->point), vertex->id, dist, neighbor->id);
00192       }
00193     }
00194     if (qh MERGING) {
00195       wmin_(Wminvertex, qh min_vertex);
00196     }
00197     qh min_vertex= minvertex;
00198     qh_settempfree(&vertices);
00199   }
00200   facets= qh_pointfacet(/*qh facet_list*/);
00201   do {
00202     old_maxoutside= fmax_(qh max_outside, maxoutside);
00203     FOREACHfacet_i_(facets) {     /* for each point with facet assignment */
00204       if (facet) {
00205         point= qh_point(facet_i);
00206         if (point == qh GOODpointp)
00207           continue;
00208         zzinc_(Ztotcheck);
00209         qh_distplane(point, facet, &dist);
00210         numpart++;
00211         bestfacet= qh_findbesthorizon(qh_IScheckmax, point, facet, !qh_NOupper, &dist, &numpart);
00212         if (bestfacet && dist > maxoutside) {
00213           if (qh ONLYgood && !bestfacet->good
00214           && !((bestfacet= qh_findgooddist(point, bestfacet, &dist, &facetlist))
00215                && dist > maxoutside))
00216             notgood++;
00217           else
00218             maxoutside= dist;
00219         }
00220         if (dist > qh TRACEdist || (bestfacet && bestfacet == qh tracefacet))
00221           qh_fprintf(qh ferr, 8094, "qh_check_maxout: p%d is %.2g above f%d\n",
00222                      qh_pointid(point), dist, bestfacet->id);
00223       }
00224     }
00225   }while
00226     (maxoutside > 2*old_maxoutside);
00227     /* if qh.maxoutside increases substantially, qh_SEARCHdist is not valid
00228           e.g., RBOX 5000 s Z1 G1e-13 t1001200614 | qhull */
00229   zzadd_(Zcheckpart, numpart);
00230   qh_settempfree(&facets);
00231   wval_(Wmaxout)= maxoutside - qh max_outside;
00232   wmax_(Wmaxoutside, qh max_outside);
00233   qh max_outside= maxoutside;
00234   qh_nearcoplanar(/*qh.facet_list*/);
00235   qh maxoutdone= True;
00236   trace1((qh ferr, 1024, "qh_check_maxout: maxoutside %2.2g, min_vertex %2.2g, outside of not good %d\n",
00237        maxoutside, qh min_vertex, notgood));
00238 } /* check_maxout */
00239 #else /* qh_NOmerge */
00240 void qh_check_maxout(void) {
00241 }
00242 #endif
00243 
00244 /*-<a                             href="qh-poly.htm#TOC"
00245   >-------------------------------</a><a name="check_output">-</a>
00246 
00247   qh_check_output()
00248     performs the checks at the end of qhull algorithm
00249     Maybe called after voronoi output.  Will recompute otherwise centrums are Voronoi centers instead
00250 */
00251 void qh_check_output(void) {
00252   int i;
00253 
00254   if (qh STOPcone)
00255     return;
00256   if (qh VERIFYoutput | qh IStracing | qh CHECKfrequently) {
00257     qh_checkpolygon(qh facet_list);
00258     qh_checkflipped_all(qh facet_list);
00259     qh_checkconvex(qh facet_list, qh_ALGORITHMfault);
00260   }else if (!qh MERGING && qh_newstats(qhstat precision, &i)) {
00261     qh_checkflipped_all(qh facet_list);
00262     qh_checkconvex(qh facet_list, qh_ALGORITHMfault);
00263   }
00264 } /* check_output */
00265 
00266 
00267 
00268 /*-<a                             href="qh-poly.htm#TOC"
00269   >-------------------------------</a><a name="check_point">-</a>
00270 
00271   qh_check_point( point, facet, maxoutside, maxdist, errfacet1, errfacet2 )
00272     check that point is less than maxoutside from facet
00273 */
00274 void qh_check_point(pointT *point, facetT *facet, realT *maxoutside, realT *maxdist, facetT **errfacet1, facetT **errfacet2) {
00275   realT dist;
00276 
00277   /* occurs after statistics reported */
00278   qh_distplane(point, facet, &dist);
00279   if (dist > *maxoutside) {
00280     if (*errfacet1 != facet) {
00281       *errfacet2= *errfacet1;
00282       *errfacet1= facet;
00283     }
00284     qh_fprintf(qh ferr, 6111, "qhull precision error: point p%d is outside facet f%d, distance= %6.8g maxoutside= %6.8g\n",
00285               qh_pointid(point), facet->id, dist, *maxoutside);
00286   }
00287   maximize_(*maxdist, dist);
00288 } /* qh_check_point */
00289 
00290 
00291 /*-<a                             href="qh-poly.htm#TOC"
00292   >-------------------------------</a><a name="check_points">-</a>
00293 
00294   qh_check_points()
00295     checks that all points are inside all facets
00296 
00297   notes:
00298     if many points and qh_check_maxout not called (i.e., !qh.MERGING),
00299        calls qh_findbesthorizon (seldom done).
00300     ignores flipped facets
00301     maxoutside includes 2 qh.DISTrounds
00302       one qh.DISTround for the computed distances in qh_check_points
00303     qh_printafacet and qh_printsummary needs only one qh.DISTround
00304     the computation for qh.VERIFYdirect does not account for qh.other_points
00305 
00306   design:
00307     if many points
00308       use qh_check_bestdist()
00309     else
00310       for all facets
00311         for all points
00312           check that point is inside facet
00313 */
00314 void qh_check_points(void) {
00315   facetT *facet, *errfacet1= NULL, *errfacet2= NULL;
00316   realT total, maxoutside, maxdist= -REALmax;
00317   pointT *point, **pointp, *pointtemp;
00318   boolT testouter;
00319 
00320   maxoutside= qh_maxouter();
00321   maxoutside += qh DISTround;
00322   /* one more qh.DISTround for check computation */
00323   trace1((qh ferr, 1025, "qh_check_points: check all points below %2.2g of all facet planes\n",
00324           maxoutside));
00325   if (qh num_good)   /* miss counts other_points and !good facets */
00326      total= (float)qh num_good * (float)qh num_points;
00327   else
00328      total= (float)qh num_facets * (float)qh num_points;
00329   if (total >= qh_VERIFYdirect && !qh maxoutdone) {
00330     if (!qh_QUICKhelp && qh SKIPcheckmax && qh MERGING)
00331       qh_fprintf(qh ferr, 7075, "qhull input warning: merging without checking outer planes('Q5' or 'Po').\n\
00332 Verify may report that a point is outside of a facet.\n");
00333     qh_check_bestdist();
00334   }else {
00335     if (qh_MAXoutside && qh maxoutdone)
00336       testouter= True;
00337     else
00338       testouter= False;
00339     if (!qh_QUICKhelp) {
00340       if (qh MERGEexact)
00341         qh_fprintf(qh ferr, 7076, "qhull input warning: exact merge ('Qx').  Verify may report that a point\n\
00342 is outside of a facet.  See qh-optq.htm#Qx\n");
00343       else if (qh SKIPcheckmax || qh NOnearinside)
00344         qh_fprintf(qh ferr, 7077, "qhull input warning: no outer plane check ('Q5') or no processing of\n\
00345 near-inside points ('Q8').  Verify may report that a point is outside\n\
00346 of a facet.\n");
00347     }
00348     if (qh PRINTprecision) {
00349       if (testouter)
00350         qh_fprintf(qh ferr, 8098, "\n\
00351 Output completed.  Verifying that all points are below outer planes of\n\
00352 all %sfacets.  Will make %2.0f distance computations.\n",
00353               (qh ONLYgood ?  "good " : ""), total);
00354       else
00355         qh_fprintf(qh ferr, 8099, "\n\
00356 Output completed.  Verifying that all points are below %2.2g of\n\
00357 all %sfacets.  Will make %2.0f distance computations.\n",
00358               maxoutside, (qh ONLYgood ?  "good " : ""), total);
00359     }
00360     FORALLfacets {
00361       if (!facet->good && qh ONLYgood)
00362         continue;
00363       if (facet->flipped)
00364         continue;
00365       if (!facet->normal) {
00366         qh_fprintf(qh ferr, 7061, "qhull warning (qh_check_points): missing normal for facet f%d\n", facet->id);
00367         continue;
00368       }
00369       if (testouter) {
00370 #if qh_MAXoutside
00371         maxoutside= facet->maxoutside + 2* qh DISTround;
00372         /* one DISTround to actual point and another to computed point */
00373 #endif
00374       }
00375       FORALLpoints {
00376         if (point != qh GOODpointp)
00377           qh_check_point(point, facet, &maxoutside, &maxdist, &errfacet1, &errfacet2);
00378       }
00379       FOREACHpoint_(qh other_points) {
00380         if (point != qh GOODpointp)
00381           qh_check_point(point, facet, &maxoutside, &maxdist, &errfacet1, &errfacet2);
00382       }
00383     }
00384     if (maxdist > qh outside_err) {
00385       qh_fprintf(qh ferr, 6112, "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",
00386                 maxdist, qh outside_err );
00387       qh_errexit2( qh_ERRprec, errfacet1, errfacet2 );
00388     }else if (errfacet1 && qh outside_err > REALmax/2)
00389         qh_errexit2( qh_ERRprec, errfacet1, errfacet2 );
00390     /* else if errfacet1, the error was logged to qh.ferr but does not effect the output */
00391     trace0((qh ferr, 21, "qh_check_points: max distance outside %2.2g\n", maxdist));
00392   }
00393 } /* check_points */
00394 
00395 
00396 /*-<a                             href="qh-poly.htm#TOC"
00397   >-------------------------------</a><a name="checkconvex">-</a>
00398 
00399   qh_checkconvex( facetlist, fault )
00400     check that each ridge in facetlist is convex
00401     fault = qh_DATAfault if reporting errors
00402           = qh_ALGORITHMfault otherwise
00403 
00404   returns:
00405     counts Zconcaveridges and Zcoplanarridges
00406     errors if concaveridge or if merging an coplanar ridge
00407 
00408   note:
00409     if not merging,
00410       tests vertices for neighboring simplicial facets
00411     else if ZEROcentrum,
00412       tests vertices for neighboring simplicial   facets
00413     else
00414       tests centrums of neighboring facets
00415 
00416   design:
00417     for all facets
00418       report flipped facets
00419       if ZEROcentrum and simplicial neighbors
00420         test vertices for neighboring simplicial facets
00421       else
00422         test centrum against all neighbors
00423 */
00424 void qh_checkconvex(facetT *facetlist, int fault) {
00425   facetT *facet, *neighbor, **neighborp, *errfacet1=NULL, *errfacet2=NULL;
00426   vertexT *vertex;
00427   realT dist;
00428   pointT *centrum;
00429   boolT waserror= False, centrum_warning= False, tempcentrum= False, allsimplicial;
00430   int neighbor_i;
00431 
00432   trace1((qh ferr, 1026, "qh_checkconvex: check all ridges are convex\n"));
00433   if (!qh RERUN) {
00434     zzval_(Zconcaveridges)= 0;
00435     zzval_(Zcoplanarridges)= 0;
00436   }
00437   FORALLfacet_(facetlist) {
00438     if (facet->flipped) {
00439       qh_precision("flipped facet");
00440       qh_fprintf(qh ferr, 6113, "qhull precision error: f%d is flipped(interior point is outside)\n",
00441                facet->id);
00442       errfacet1= facet;
00443       waserror= True;
00444       continue;
00445     }
00446     if (qh MERGING && (!qh ZEROcentrum || !facet->simplicial || facet->tricoplanar))
00447       allsimplicial= False;
00448     else {
00449       allsimplicial= True;
00450       neighbor_i= 0;
00451       FOREACHneighbor_(facet) {
00452         vertex= SETelemt_(facet->vertices, neighbor_i++, vertexT);
00453         if (!neighbor->simplicial || neighbor->tricoplanar) {
00454           allsimplicial= False;
00455           continue;
00456         }
00457         qh_distplane(vertex->point, neighbor, &dist);
00458         if (dist > -qh DISTround) {
00459           if (fault == qh_DATAfault) {
00460             qh_precision("coplanar or concave ridge");
00461             qh_fprintf(qh ferr, 6114, "qhull precision error: initial simplex is not convex. Distance=%.2g\n", dist);
00462             qh_errexit(qh_ERRsingular, NULL, NULL);
00463           }
00464           if (dist > qh DISTround) {
00465             zzinc_(Zconcaveridges);
00466             qh_precision("concave ridge");
00467             qh_fprintf(qh ferr, 6115, "qhull precision error: f%d is concave to f%d, since p%d(v%d) is %6.4g above\n",
00468               facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist);
00469             errfacet1= facet;
00470             errfacet2= neighbor;
00471             waserror= True;
00472           }else if (qh ZEROcentrum) {
00473             if (dist > 0) {     /* qh_checkzero checks that dist < - qh DISTround */
00474               zzinc_(Zcoplanarridges);
00475               qh_precision("coplanar ridge");
00476               qh_fprintf(qh ferr, 6116, "qhull precision error: f%d is clearly not convex to f%d, since p%d(v%d) is %6.4g above\n",
00477                 facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist);
00478               errfacet1= facet;
00479               errfacet2= neighbor;
00480               waserror= True;
00481             }
00482           }else {
00483             zzinc_(Zcoplanarridges);
00484             qh_precision("coplanar ridge");
00485             trace0((qh ferr, 22, "qhull precision error: f%d may be coplanar to f%d, since p%d(v%d) is within %6.4g during p%d\n",
00486               facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist, qh furthest_id));
00487           }
00488         }
00489       }
00490     }
00491     if (!allsimplicial) {
00492       if (qh CENTERtype == qh_AScentrum) {
00493         if (!facet->center)
00494           facet->center= qh_getcentrum(facet);
00495         centrum= facet->center;
00496       }else {
00497         if (!centrum_warning && (!facet->simplicial || facet->tricoplanar)) {
00498            centrum_warning= True;
00499            qh_fprintf(qh ferr, 7062, "qhull warning: recomputing centrums for convexity test.  This may lead to false, precision errors.\n");
00500         }
00501         centrum= qh_getcentrum(facet);
00502         tempcentrum= True;
00503       }
00504       FOREACHneighbor_(facet) {
00505         if (qh ZEROcentrum && facet->simplicial && neighbor->simplicial)
00506           continue;
00507         if (facet->tricoplanar || neighbor->tricoplanar)
00508           continue;
00509         zzinc_(Zdistconvex);
00510         qh_distplane(centrum, neighbor, &dist);
00511         if (dist > qh DISTround) {
00512           zzinc_(Zconcaveridges);
00513           qh_precision("concave ridge");
00514           qh_fprintf(qh ferr, 6117, "qhull precision error: f%d is concave to f%d.  Centrum of f%d is %6.4g above f%d\n",
00515             facet->id, neighbor->id, facet->id, dist, neighbor->id);
00516           errfacet1= facet;
00517           errfacet2= neighbor;
00518           waserror= True;
00519         }else if (dist >= 0.0) {   /* if arithmetic always rounds the same,
00520                                      can test against centrum radius instead */
00521           zzinc_(Zcoplanarridges);
00522           qh_precision("coplanar ridge");
00523           qh_fprintf(qh ferr, 6118, "qhull precision error: f%d is coplanar or concave to f%d.  Centrum of f%d is %6.4g above f%d\n",
00524             facet->id, neighbor->id, facet->id, dist, neighbor->id);
00525           errfacet1= facet;
00526           errfacet2= neighbor;
00527           waserror= True;
00528         }
00529       }
00530       if (tempcentrum)
00531         qh_memfree(centrum, qh normal_size);
00532     }
00533   }
00534   if (waserror && !qh FORCEoutput)
00535     qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);
00536 } /* checkconvex */
00537 
00538 
00539 /*-<a                             href="qh-poly.htm#TOC"
00540   >-------------------------------</a><a name="checkfacet">-</a>
00541 
00542   qh_checkfacet( facet, newmerge, waserror )
00543     checks for consistency errors in facet
00544     newmerge set if from merge.c
00545 
00546   returns:
00547     sets waserror if any error occurs
00548 
00549   checks:
00550     vertex ids are inverse sorted
00551     unless newmerge, at least hull_dim neighbors and vertices (exactly if simplicial)
00552     if non-simplicial, at least as many ridges as neighbors
00553     neighbors are not duplicated
00554     ridges are not duplicated
00555     in 3-d, ridges=verticies
00556     (qh.hull_dim-1) ridge vertices
00557     neighbors are reciprocated
00558     ridge neighbors are facet neighbors and a ridge for every neighbor
00559     simplicial neighbors match facetintersect
00560     vertex intersection matches vertices of common ridges
00561     vertex neighbors and facet vertices agree
00562     all ridges have distinct vertex sets
00563 
00564   notes:
00565     uses neighbor->seen
00566 
00567   design:
00568     check sets
00569     check vertices
00570     check sizes of neighbors and vertices
00571     check for qh_MERGEridge and qh_DUPLICATEridge flags
00572     check neighbor set
00573     check ridge set
00574     check ridges, neighbors, and vertices
00575 */
00576 void qh_checkfacet(facetT *facet, boolT newmerge, boolT *waserrorp) {
00577   facetT *neighbor, **neighborp, *errother=NULL;
00578   ridgeT *ridge, **ridgep, *errridge= NULL, *ridge2;
00579   vertexT *vertex, **vertexp;
00580   unsigned previousid= INT_MAX;
00581   int numneighbors, numvertices, numridges=0, numRvertices=0;
00582   boolT waserror= False;
00583   int skipA, skipB, ridge_i, ridge_n, i;
00584   setT *intersection;
00585 
00586   if (facet->visible) {
00587     qh_fprintf(qh ferr, 6119, "qhull internal error (qh_checkfacet): facet f%d is on the visible_list\n",
00588       facet->id);
00589     qh_errexit(qh_ERRqhull, facet, NULL);
00590   }
00591   if (!facet->normal) {
00592     qh_fprintf(qh ferr, 6120, "qhull internal error (qh_checkfacet): facet f%d does not have  a normal\n",
00593       facet->id);
00594     waserror= True;
00595   }
00596   qh_setcheck(facet->vertices, "vertices for f", facet->id);
00597   qh_setcheck(facet->ridges, "ridges for f", facet->id);
00598   qh_setcheck(facet->outsideset, "outsideset for f", facet->id);
00599   qh_setcheck(facet->coplanarset, "coplanarset for f", facet->id);
00600   qh_setcheck(facet->neighbors, "neighbors for f", facet->id);
00601   FOREACHvertex_(facet->vertices) {
00602     if (vertex->deleted) {
00603       qh_fprintf(qh ferr, 6121, "qhull internal error (qh_checkfacet): deleted vertex v%d in f%d\n", vertex->id, facet->id);
00604       qh_errprint("ERRONEOUS", NULL, NULL, NULL, vertex);
00605       waserror= True;
00606     }
00607     if (vertex->id >= previousid) {
00608       qh_fprintf(qh ferr, 6122, "qhull internal error (qh_checkfacet): vertices of f%d are not in descending id order at v%d\n", facet->id, vertex->id);
00609       waserror= True;
00610       break;
00611     }
00612     previousid= vertex->id;
00613   }
00614   numneighbors= qh_setsize(facet->neighbors);
00615   numvertices= qh_setsize(facet->vertices);
00616   numridges= qh_setsize(facet->ridges);
00617   if (facet->simplicial) {
00618     if (numvertices+numneighbors != 2*qh hull_dim
00619     && !facet->degenerate && !facet->redundant) {
00620       qh_fprintf(qh ferr, 6123, "qhull internal error (qh_checkfacet): for simplicial facet f%d, #vertices %d + #neighbors %d != 2*qh hull_dim\n",
00621                 facet->id, numvertices, numneighbors);
00622       qh_setprint(qh ferr, "", facet->neighbors);
00623       waserror= True;
00624     }
00625   }else { /* non-simplicial */
00626     if (!newmerge
00627     &&(numvertices < qh hull_dim || numneighbors < qh hull_dim)
00628     && !facet->degenerate && !facet->redundant) {
00629       qh_fprintf(qh ferr, 6124, "qhull internal error (qh_checkfacet): for facet f%d, #vertices %d or #neighbors %d < qh hull_dim\n",
00630          facet->id, numvertices, numneighbors);
00631        waserror= True;
00632     }
00633     /* in 3-d, can get a vertex twice in an edge list, e.g., RBOX 1000 s W1e-13 t995849315 D2 | QHULL d Tc Tv TP624 TW1e-13 T4 */
00634     if (numridges < numneighbors
00635     ||(qh hull_dim == 3 && numvertices > numridges && !qh NEWfacets)
00636     ||(qh hull_dim == 2 && numridges + numvertices + numneighbors != 6)) {
00637       if (!facet->degenerate && !facet->redundant) {
00638         qh_fprintf(qh ferr, 6125, "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",
00639             facet->id, numridges, numneighbors, numvertices);
00640         waserror= True;
00641       }
00642     }
00643   }
00644   FOREACHneighbor_(facet) {
00645     if (neighbor == qh_MERGEridge || neighbor == qh_DUPLICATEridge) {
00646       qh_fprintf(qh ferr, 6126, "qhull internal error (qh_checkfacet): facet f%d still has a MERGE or DUP neighbor\n", facet->id);
00647       qh_errexit(qh_ERRqhull, facet, NULL);
00648     }
00649     neighbor->seen= True;
00650   }
00651   FOREACHneighbor_(facet) {
00652     if (!qh_setin(neighbor->neighbors, facet)) {
00653       qh_fprintf(qh ferr, 6127, "qhull internal error (qh_checkfacet): facet f%d has neighbor f%d, but f%d does not have neighbor f%d\n",
00654               facet->id, neighbor->id, neighbor->id, facet->id);
00655       errother= neighbor;
00656       waserror= True;
00657     }
00658     if (!neighbor->seen) {
00659       qh_fprintf(qh ferr, 6128, "qhull internal error (qh_checkfacet): facet f%d has a duplicate neighbor f%d\n",
00660               facet->id, neighbor->id);
00661       errother= neighbor;
00662       waserror= True;
00663     }
00664     neighbor->seen= False;
00665   }
00666   FOREACHridge_(facet->ridges) {
00667     qh_setcheck(ridge->vertices, "vertices for r", ridge->id);
00668     ridge->seen= False;
00669   }
00670   FOREACHridge_(facet->ridges) {
00671     if (ridge->seen) {
00672       qh_fprintf(qh ferr, 6129, "qhull internal error (qh_checkfacet): facet f%d has a duplicate ridge r%d\n",
00673               facet->id, ridge->id);
00674       errridge= ridge;
00675       waserror= True;
00676     }
00677     ridge->seen= True;
00678     numRvertices= qh_setsize(ridge->vertices);
00679     if (numRvertices != qh hull_dim - 1) {
00680       qh_fprintf(qh ferr, 6130, "qhull internal error (qh_checkfacet): ridge between f%d and f%d has %d vertices\n",
00681                 ridge->top->id, ridge->bottom->id, numRvertices);
00682       errridge= ridge;
00683       waserror= True;
00684     }
00685     neighbor= otherfacet_(ridge, facet);
00686     neighbor->seen= True;
00687     if (!qh_setin(facet->neighbors, neighbor)) {
00688       qh_fprintf(qh ferr, 6131, "qhull internal error (qh_checkfacet): for facet f%d, neighbor f%d of ridge r%d not in facet\n",
00689            facet->id, neighbor->id, ridge->id);
00690       errridge= ridge;
00691       waserror= True;
00692     }
00693   }
00694   if (!facet->simplicial) {
00695     FOREACHneighbor_(facet) {
00696       if (!neighbor->seen) {
00697         qh_fprintf(qh ferr, 6132, "qhull internal error (qh_checkfacet): facet f%d does not have a ridge for neighbor f%d\n",
00698               facet->id, neighbor->id);
00699         errother= neighbor;
00700         waserror= True;
00701       }
00702       intersection= qh_vertexintersect_new(facet->vertices, neighbor->vertices);
00703       qh_settemppush(intersection);
00704       FOREACHvertex_(facet->vertices) {
00705         vertex->seen= False;
00706         vertex->seen2= False;
00707       }
00708       FOREACHvertex_(intersection)
00709         vertex->seen= True;
00710       FOREACHridge_(facet->ridges) {
00711         if (neighbor != otherfacet_(ridge, facet))
00712             continue;
00713         FOREACHvertex_(ridge->vertices) {
00714           if (!vertex->seen) {
00715             qh_fprintf(qh ferr, 6133, "qhull internal error (qh_checkfacet): vertex v%d in r%d not in f%d intersect f%d\n",
00716                   vertex->id, ridge->id, facet->id, neighbor->id);
00717             qh_errexit(qh_ERRqhull, facet, ridge);
00718           }
00719           vertex->seen2= True;
00720         }
00721       }
00722       if (!newmerge) {
00723         FOREACHvertex_(intersection) {
00724           if (!vertex->seen2) {
00725             if (qh IStracing >=3 || !qh MERGING) {
00726               qh_fprintf(qh ferr, 6134, "qhull precision error (qh_checkfacet): vertex v%d in f%d intersect f%d but\n\
00727  not in a ridge.  This is ok under merging.  Last point was p%d\n",
00728                      vertex->id, facet->id, neighbor->id, qh furthest_id);
00729               if (!qh FORCEoutput && !qh MERGING) {
00730                 qh_errprint("ERRONEOUS", facet, neighbor, NULL, vertex);
00731                 if (!qh MERGING)
00732                   qh_errexit(qh_ERRqhull, NULL, NULL);
00733               }
00734             }
00735           }
00736         }
00737       }
00738       qh_settempfree(&intersection);
00739     }
00740   }else { /* simplicial */
00741     FOREACHneighbor_(facet) {
00742       if (neighbor->simplicial) {
00743         skipA= SETindex_(facet->neighbors, neighbor);
00744         skipB= qh_setindex(neighbor->neighbors, facet);
00745         if (!qh_setequal_skip(facet->vertices, skipA, neighbor->vertices, skipB)) {
00746           qh_fprintf(qh ferr, 6135, "qhull internal error (qh_checkfacet): facet f%d skip %d and neighbor f%d skip %d do not match \n",
00747                    facet->id, skipA, neighbor->id, skipB);
00748           errother= neighbor;
00749           waserror= True;
00750         }
00751       }
00752     }
00753   }
00754   if (qh hull_dim < 5 && (qh IStracing > 2 || qh CHECKfrequently)) {
00755     FOREACHridge_i_(facet->ridges) {           /* expensive */
00756       for (i=ridge_i+1; i < ridge_n; i++) {
00757         ridge2= SETelemt_(facet->ridges, i, ridgeT);
00758         if (qh_setequal(ridge->vertices, ridge2->vertices)) {
00759           qh_fprintf(qh ferr, 6227, "Qhull internal error (qh_checkfacet): ridges r%d and r%d have the same vertices\n",
00760                   ridge->id, ridge2->id);
00761           errridge= ridge;
00762           waserror= True;
00763         }
00764       }
00765     }
00766   }
00767   if (waserror) {
00768     qh_errprint("ERRONEOUS", facet, errother, errridge, NULL);
00769     *waserrorp= True;
00770   }
00771 } /* checkfacet */
00772 
00773 
00774 /*-<a                             href="qh-poly.htm#TOC"
00775   >-------------------------------</a><a name="checkflipped_all">-</a>
00776 
00777   qh_checkflipped_all( facetlist )
00778     checks orientation of facets in list against interior point
00779 */
00780 void qh_checkflipped_all(facetT *facetlist) {
00781   facetT *facet;
00782   boolT waserror= False;
00783   realT dist;
00784 
00785   if (facetlist == qh facet_list)
00786     zzval_(Zflippedfacets)= 0;
00787   FORALLfacet_(facetlist) {
00788     if (facet->normal && !qh_checkflipped(facet, &dist, !qh_ALL)) {
00789       qh_fprintf(qh ferr, 6136, "qhull precision error: facet f%d is flipped, distance= %6.12g\n",
00790               facet->id, dist);
00791       if (!qh FORCEoutput) {
00792         qh_errprint("ERRONEOUS", facet, NULL, NULL, NULL);
00793         waserror= True;
00794       }
00795     }
00796   }
00797   if (waserror) {
00798     qh_fprintf(qh ferr, 8101, "\n\
00799 A flipped facet occurs when its distance to the interior point is\n\
00800 greater than %2.2g, the maximum roundoff error.\n", -qh DISTround);
00801     qh_errexit(qh_ERRprec, NULL, NULL);
00802   }
00803 } /* checkflipped_all */
00804 
00805 /*-<a                             href="qh-poly.htm#TOC"
00806   >-------------------------------</a><a name="checkpolygon">-</a>
00807 
00808   qh_checkpolygon( facetlist )
00809     checks the correctness of the structure
00810 
00811   notes:
00812     call with either qh.facet_list or qh.newfacet_list
00813     checks num_facets and num_vertices if qh.facet_list
00814 
00815   design:
00816     for each facet
00817       checks facet and outside set
00818     initializes vertexlist
00819     for each facet
00820       checks vertex set
00821     if checking all facets(qh.facetlist)
00822       check facet count
00823       if qh.VERTEXneighbors
00824         check vertex neighbors and count
00825       check vertex count
00826 */
00827 void qh_checkpolygon(facetT *facetlist) {
00828   facetT *facet;
00829   vertexT *vertex, **vertexp, *vertexlist;
00830   int numfacets= 0, numvertices= 0, numridges= 0;
00831   int totvneighbors= 0, totvertices= 0;
00832   boolT waserror= False, nextseen= False, visibleseen= False;
00833 
00834   trace1((qh ferr, 1027, "qh_checkpolygon: check all facets from f%d\n", facetlist->id));
00835   if (facetlist != qh facet_list || qh ONLYgood)
00836     nextseen= True;
00837   FORALLfacet_(facetlist) {
00838     if (facet == qh visible_list)
00839       visibleseen= True;
00840     if (!facet->visible) {
00841       if (!nextseen) {
00842         if (facet == qh facet_next)
00843           nextseen= True;
00844         else if (qh_setsize(facet->outsideset)) {
00845           if (!qh NARROWhull
00846 #if !qh_COMPUTEfurthest
00847                || facet->furthestdist >= qh MINoutside
00848 #endif
00849                         ) {
00850             qh_fprintf(qh ferr, 6137, "qhull internal error (qh_checkpolygon): f%d has outside points before qh facet_next\n",
00851                      facet->id);
00852             qh_errexit(qh_ERRqhull, facet, NULL);
00853           }
00854         }
00855       }
00856       numfacets++;
00857       qh_checkfacet(facet, False, &waserror);
00858     }
00859   }
00860   if (qh visible_list && !visibleseen && facetlist == qh facet_list) {
00861     qh_fprintf(qh ferr, 6138, "qhull internal error (qh_checkpolygon): visible list f%d no longer on facet list\n", qh visible_list->id);
00862     qh_printlists();
00863     qh_errexit(qh_ERRqhull, qh visible_list, NULL);
00864   }
00865   if (facetlist == qh facet_list)
00866     vertexlist= qh vertex_list;
00867   else if (facetlist == qh newfacet_list)
00868     vertexlist= qh newvertex_list;
00869   else
00870     vertexlist= NULL;
00871   FORALLvertex_(vertexlist) {
00872     vertex->seen= False;
00873     vertex->visitid= 0;
00874   }
00875   FORALLfacet_(facetlist) {
00876     if (facet->visible)
00877       continue;
00878     if (facet->simplicial)
00879       numridges += qh hull_dim;
00880     else
00881       numridges += qh_setsize(facet->ridges);
00882     FOREACHvertex_(facet->vertices) {
00883       vertex->visitid++;
00884       if (!vertex->seen) {
00885         vertex->seen= True;
00886         numvertices++;
00887         if (qh_pointid(vertex->point) == -1) {
00888           qh_fprintf(qh ferr, 6139, "qhull internal error (qh_checkpolygon): unknown point %p for vertex v%d first_point %p\n",
00889                    vertex->point, vertex->id, qh first_point);
00890           waserror= True;
00891         }
00892       }
00893     }
00894   }
00895   qh vertex_visit += (unsigned int)numfacets;
00896   if (facetlist == qh facet_list) {
00897     if (numfacets != qh num_facets - qh num_visible) {
00898       qh_fprintf(qh ferr, 6140, "qhull internal error (qh_checkpolygon): actual number of facets is %d, cumulative facet count is %d - %d visible facets\n",
00899               numfacets, qh num_facets, qh num_visible);
00900       waserror= True;
00901     }
00902     qh vertex_visit++;
00903     if (qh VERTEXneighbors) {
00904       FORALLvertices {
00905         qh_setcheck(vertex->neighbors, "neighbors for v", vertex->id);
00906         if (vertex->deleted)
00907           continue;
00908         totvneighbors += qh_setsize(vertex->neighbors);
00909       }
00910       FORALLfacet_(facetlist)
00911         totvertices += qh_setsize(facet->vertices);
00912       if (totvneighbors != totvertices) {
00913         qh_fprintf(qh ferr, 6141, "qhull internal error (qh_checkpolygon): vertex neighbors inconsistent.  Totvneighbors %d, totvertices %d\n",
00914                 totvneighbors, totvertices);
00915         waserror= True;
00916       }
00917     }
00918     if (numvertices != qh num_vertices - qh_setsize(qh del_vertices)) {
00919       qh_fprintf(qh ferr, 6142, "qhull internal error (qh_checkpolygon): actual number of vertices is %d, cumulative vertex count is %d\n",
00920               numvertices, qh num_vertices - qh_setsize(qh del_vertices));
00921       waserror= True;
00922     }
00923     if (qh hull_dim == 2 && numvertices != numfacets) {
00924       qh_fprintf(qh ferr, 6143, "qhull internal error (qh_checkpolygon): #vertices %d != #facets %d\n",
00925         numvertices, numfacets);
00926       waserror= True;
00927     }
00928     if (qh hull_dim == 3 && numvertices + numfacets - numridges/2 != 2) {
00929       qh_fprintf(qh ferr, 7063, "qhull warning: #vertices %d + #facets %d - #edges %d != 2\n\
00930         A vertex appears twice in a edge list.  May occur during merging.",
00931         numvertices, numfacets, numridges/2);
00932       /* occurs if lots of merging and a vertex ends up twice in an edge list.  e.g., RBOX 1000 s W1e-13 t995849315 D2 | QHULL d Tc Tv */
00933     }
00934   }
00935   if (waserror)
00936     qh_errexit(qh_ERRqhull, NULL, NULL);
00937 } /* checkpolygon */
00938 
00939 
00940 /*-<a                             href="qh-poly.htm#TOC"
00941   >-------------------------------</a><a name="checkvertex">-</a>
00942 
00943   qh_checkvertex( vertex )
00944     check vertex for consistency
00945     checks vertex->neighbors
00946 
00947   notes:
00948     neighbors checked efficiently in checkpolygon
00949 */
00950 void qh_checkvertex(vertexT *vertex) {
00951   boolT waserror= False;
00952   facetT *neighbor, **neighborp, *errfacet=NULL;
00953 
00954   if (qh_pointid(vertex->point) == -1) {
00955     qh_fprintf(qh ferr, 6144, "qhull internal error (qh_checkvertex): unknown point id %p\n", vertex->point);
00956     waserror= True;
00957   }
00958   if (vertex->id >= qh vertex_id) {
00959     qh_fprintf(qh ferr, 6145, "qhull internal error (qh_checkvertex): unknown vertex id %d\n", vertex->id);
00960     waserror= True;
00961   }
00962   if (!waserror && !vertex->deleted) {
00963     if (qh_setsize(vertex->neighbors)) {
00964       FOREACHneighbor_(vertex) {
00965         if (!qh_setin(neighbor->vertices, vertex)) {
00966           qh_fprintf(qh ferr, 6146, "qhull internal error (qh_checkvertex): neighbor f%d does not contain v%d\n", neighbor->id, vertex->id);
00967           errfacet= neighbor;
00968           waserror= True;
00969         }
00970       }
00971     }
00972   }
00973   if (waserror) {
00974     qh_errprint("ERRONEOUS", NULL, NULL, NULL, vertex);
00975     qh_errexit(qh_ERRqhull, errfacet, NULL);
00976   }
00977 } /* checkvertex */
00978 
00979 /*-<a                             href="qh-poly.htm#TOC"
00980   >-------------------------------</a><a name="clearcenters">-</a>
00981 
00982   qh_clearcenters( type )
00983     clear old data from facet->center
00984 
00985   notes:
00986     sets new centertype
00987     nop if CENTERtype is the same
00988 */
00989 void qh_clearcenters(qh_CENTER type) {
00990   facetT *facet;
00991 
00992   if (qh CENTERtype != type) {
00993     FORALLfacets {
00994       if (facet->tricoplanar && !facet->keepcentrum)
00995           facet->center= NULL;
00996       else if (qh CENTERtype == qh_ASvoronoi){
00997         if (facet->center) {
00998           qh_memfree(facet->center, qh center_size);
00999           facet->center= NULL;
01000         }
01001       }else /* qh CENTERtype == qh_AScentrum */ {
01002         if (facet->center) {
01003           qh_memfree(facet->center, qh normal_size);
01004           facet->center= NULL;
01005         }
01006       }
01007     }
01008     qh CENTERtype= type;
01009   }
01010   trace2((qh ferr, 2043, "qh_clearcenters: switched to center type %d\n", type));
01011 } /* clearcenters */
01012 
01013 /*-<a                             href="qh-poly.htm#TOC"
01014   >-------------------------------</a><a name="createsimplex">-</a>
01015 
01016   qh_createsimplex( vertices )
01017     creates a simplex from a set of vertices
01018 
01019   returns:
01020     initializes qh.facet_list to the simplex
01021     initializes qh.newfacet_list, .facet_tail
01022     initializes qh.vertex_list, .newvertex_list, .vertex_tail
01023 
01024   design:
01025     initializes lists
01026     for each vertex
01027       create a new facet
01028     for each new facet
01029       create its neighbor set
01030 */
01031 void qh_createsimplex(setT *vertices) {
01032   facetT *facet= NULL, *newfacet;
01033   boolT toporient= True;
01034   int vertex_i, vertex_n, nth;
01035   setT *newfacets= qh_settemp(qh hull_dim+1);
01036   vertexT *vertex;
01037 
01038   qh facet_list= qh newfacet_list= qh facet_tail= qh_newfacet();
01039   qh num_facets= qh num_vertices= qh num_visible= 0;
01040   qh vertex_list= qh newvertex_list= qh vertex_tail= qh_newvertex(NULL);
01041   FOREACHvertex_i_(vertices) {
01042     newfacet= qh_newfacet();
01043     newfacet->vertices= qh_setnew_delnthsorted(vertices, vertex_n,
01044                                                 vertex_i, 0);
01045     newfacet->toporient= (unsigned char)toporient;
01046     qh_appendfacet(newfacet);
01047     newfacet->newfacet= True;
01048     qh_appendvertex(vertex);
01049     qh_setappend(&newfacets, newfacet);
01050     toporient ^= True;
01051   }
01052   FORALLnew_facets {
01053     nth= 0;
01054     FORALLfacet_(qh newfacet_list) {
01055       if (facet != newfacet)
01056         SETelem_(newfacet->neighbors, nth++)= facet;
01057     }
01058     qh_settruncate(newfacet->neighbors, qh hull_dim);
01059   }
01060   qh_settempfree(&newfacets);
01061   trace1((qh ferr, 1028, "qh_createsimplex: created simplex\n"));
01062 } /* createsimplex */
01063 
01064 /*-<a                             href="qh-poly.htm#TOC"
01065   >-------------------------------</a><a name="delridge">-</a>
01066 
01067   qh_delridge( ridge )
01068     deletes ridge from data structures it belongs to
01069     frees up its memory
01070 
01071   notes:
01072     in merge.c, caller sets vertex->delridge for each vertex
01073     ridges also freed in qh_freeqhull
01074 */
01075 void qh_delridge(ridgeT *ridge) {
01076   void **freelistp; /* used !qh_NOmem */
01077 
01078   qh_setdel(ridge->top->ridges, ridge);
01079   qh_setdel(ridge->bottom->ridges, ridge);
01080   qh_setfree(&(ridge->vertices));
01081   qh_memfree_(ridge, (int)sizeof(ridgeT), freelistp);
01082 } /* delridge */
01083 
01084 
01085 /*-<a                             href="qh-poly.htm#TOC"
01086   >-------------------------------</a><a name="delvertex">-</a>
01087 
01088   qh_delvertex( vertex )
01089     deletes a vertex and frees its memory
01090 
01091   notes:
01092     assumes vertex->adjacencies have been updated if needed
01093     unlinks from vertex_list
01094 */
01095 void qh_delvertex(vertexT *vertex) {
01096 
01097   if (vertex == qh tracevertex)
01098     qh tracevertex= NULL;
01099   qh_removevertex(vertex);
01100   qh_setfree(&vertex->neighbors);
01101   qh_memfree(vertex, (int)sizeof(vertexT));
01102 } /* delvertex */
01103 
01104 
01105 /*-<a                             href="qh-poly.htm#TOC"
01106   >-------------------------------</a><a name="facet3vertex">-</a>
01107 
01108   qh_facet3vertex(  )
01109     return temporary set of 3-d vertices in qh_ORIENTclock order
01110 
01111   design:
01112     if simplicial facet
01113       build set from facet->vertices with facet->toporient
01114     else
01115       for each ridge in order
01116         build set from ridge's vertices
01117 */
01118 setT *qh_facet3vertex(facetT *facet) {
01119   ridgeT *ridge, *firstridge;
01120   vertexT *vertex;
01121   int cntvertices, cntprojected=0;
01122   setT *vertices;
01123 
01124   cntvertices= qh_setsize(facet->vertices);
01125   vertices= qh_settemp(cntvertices);
01126   if (facet->simplicial) {
01127     if (cntvertices != 3) {
01128       qh_fprintf(qh ferr, 6147, "qhull internal error (qh_facet3vertex): only %d vertices for simplicial facet f%d\n",
01129                   cntvertices, facet->id);
01130       qh_errexit(qh_ERRqhull, facet, NULL);
01131     }
01132     qh_setappend(&vertices, SETfirst_(facet->vertices));
01133     if (facet->toporient ^ qh_ORIENTclock)
01134       qh_setappend(&vertices, SETsecond_(facet->vertices));
01135     else
01136       qh_setaddnth(&vertices, 0, SETsecond_(facet->vertices));
01137     qh_setappend(&vertices, SETelem_(facet->vertices, 2));
01138   }else {
01139     ridge= firstridge= SETfirstt_(facet->ridges, ridgeT);   /* no infinite */
01140     while ((ridge= qh_nextridge3d(ridge, facet, &vertex))) {
01141       qh_setappend(&vertices, vertex);
01142       if (++cntprojected > cntvertices || ridge == firstridge)
01143         break;
01144     }
01145     if (!ridge || cntprojected != cntvertices) {
01146       qh_fprintf(qh ferr, 6148, "qhull internal error (qh_facet3vertex): ridges for facet %d don't match up.  got at least %d\n",
01147                   facet->id, cntprojected);
01148       qh_errexit(qh_ERRqhull, facet, ridge);
01149     }
01150   }
01151   return vertices;
01152 } /* facet3vertex */
01153 
01154 /*-<a                             href="qh-poly.htm#TOC"
01155   >-------------------------------</a><a name="findbestfacet">-</a>
01156 
01157   qh_findbestfacet( point, bestoutside, bestdist, isoutside )
01158     find facet that is furthest below a point
01159 
01160     for Delaunay triangulations,
01161       Use qh_setdelaunay() to lift point to paraboloid and scale by 'Qbb' if needed
01162       Do not use options 'Qbk', 'QBk', or 'QbB' since they scale the coordinates.
01163 
01164   returns:
01165     if bestoutside is set (e.g., qh_ALL)
01166       returns best facet that is not upperdelaunay
01167       if Delaunay and inside, point is outside circumsphere of bestfacet
01168     else
01169       returns first facet below point
01170       if point is inside, returns nearest, !upperdelaunay facet
01171     distance to facet
01172     isoutside set if outside of facet
01173 
01174   notes:
01175     For tricoplanar facets, this finds one of the tricoplanar facets closest
01176     to the point.  For Delaunay triangulations, the point may be inside a
01177     different tricoplanar facet. See <a href="../html/qh-code.htm#findfacet">locate a facet with qh_findbestfacet()</a>
01178 
01179     If inside, qh_findbestfacet performs an exhaustive search
01180        this may be too conservative.  Sometimes it is clearly required.
01181 
01182     qh_findbestfacet is not used by qhull.
01183     uses qh.visit_id and qh.coplanarset
01184 
01185   see:
01186     <a href="geom.c#findbest">qh_findbest</a>
01187 */
01188 facetT *qh_findbestfacet(pointT *point, boolT bestoutside,
01189            realT *bestdist, boolT *isoutside) {
01190   facetT *bestfacet= NULL;
01191   int numpart, totpart= 0;
01192 
01193   bestfacet= qh_findbest(point, qh facet_list,
01194                             bestoutside, !qh_ISnewfacets, bestoutside /* qh_NOupper */,
01195                             bestdist, isoutside, &totpart);
01196   if (*bestdist < -qh DISTround) {
01197     bestfacet= qh_findfacet_all(point, bestdist, isoutside, &numpart);
01198     totpart += numpart;
01199     if ((isoutside && bestoutside)
01200     || (!isoutside && bestfacet->upperdelaunay)) {
01201       bestfacet= qh_findbest(point, bestfacet,
01202                             bestoutside, False, bestoutside,
01203                             bestdist, isoutside, &totpart);
01204       totpart += numpart;
01205     }
01206   }
01207   trace3((qh ferr, 3014, "qh_findbestfacet: f%d dist %2.2g isoutside %d totpart %d\n",
01208           bestfacet->id, *bestdist, *isoutside, totpart));
01209   return bestfacet;
01210 } /* findbestfacet */
01211 
01212 /*-<a                             href="qh-poly.htm#TOC"
01213   >-------------------------------</a><a name="findbestlower">-</a>
01214 
01215   qh_findbestlower( facet, point, bestdist, numpart )
01216     returns best non-upper, non-flipped neighbor of facet for point
01217     if needed, searches vertex neighbors
01218 
01219   returns:
01220     returns bestdist and updates numpart
01221 
01222   notes:
01223     if Delaunay and inside, point is outside of circumsphere of bestfacet
01224     called by qh_findbest() for points above an upperdelaunay facet
01225 
01226 */
01227 facetT *qh_findbestlower(facetT *upperfacet, pointT *point, realT *bestdistp, int *numpart) {
01228   facetT *neighbor, **neighborp, *bestfacet= NULL;
01229   realT bestdist= -REALmax/2 /* avoid underflow */;
01230   realT dist;
01231   vertexT *vertex;
01232 
01233   zinc_(Zbestlower);
01234   FOREACHneighbor_(upperfacet) {
01235     if (neighbor->upperdelaunay || neighbor->flipped)
01236       continue;
01237     (*numpart)++;
01238     qh_distplane(point, neighbor, &dist);
01239     if (dist > bestdist) {
01240       bestfacet= neighbor;
01241       bestdist= dist;
01242     }
01243   }
01244   if (!bestfacet) {
01245     zinc_(Zbestlowerv);
01246     /* rarely called, numpart does not count nearvertex computations */
01247     vertex= qh_nearvertex(upperfacet, point, &dist);
01248     qh_vertexneighbors();
01249     FOREACHneighbor_(vertex) {
01250       if (neighbor->upperdelaunay || neighbor->flipped)
01251         continue;
01252       (*numpart)++;
01253       qh_distplane(point, neighbor, &dist);
01254       if (dist > bestdist) {
01255         bestfacet= neighbor;
01256         bestdist= dist;
01257       }
01258     }
01259   }
01260   if (!bestfacet) {
01261     qh_fprintf(qh ferr, 6228, "\n\
01262 Qhull internal error (qh_findbestlower): all neighbors of facet %d are flipped or upper Delaunay.\n\
01263 Please report this error to qhull_bug@qhull.org with the input and all of the output.\n",
01264        upperfacet->id);
01265     qh_errexit(qh_ERRqhull, upperfacet, NULL);
01266   }
01267   *bestdistp= bestdist;
01268   trace3((qh ferr, 3015, "qh_findbestlower: f%d dist %2.2g for f%d p%d\n",
01269           bestfacet->id, bestdist, upperfacet->id, qh_pointid(point)));
01270   return bestfacet;
01271 } /* findbestlower */
01272 
01273 /*-<a                             href="qh-poly.htm#TOC"
01274   >-------------------------------</a><a name="findfacet_all">-</a>
01275 
01276   qh_findfacet_all( point, bestdist, isoutside, numpart )
01277     exhaustive search for facet below a point
01278 
01279     for Delaunay triangulations,
01280       Use qh_setdelaunay() to lift point to paraboloid and scale by 'Qbb' if needed
01281       Do not use options 'Qbk', 'QBk', or 'QbB' since they scale the coordinates.
01282 
01283   returns:
01284     returns first facet below point
01285     if point is inside,
01286       returns nearest facet
01287     distance to facet
01288     isoutside if point is outside of the hull
01289     number of distance tests
01290 
01291   notes:
01292     for library users, not used by Qhull
01293 */
01294 facetT *qh_findfacet_all(pointT *point, realT *bestdist, boolT *isoutside,
01295                           int *numpart) {
01296   facetT *bestfacet= NULL, *facet;
01297   realT dist;
01298   int totpart= 0;
01299 
01300   *bestdist= -REALmax;
01301   *isoutside= False;
01302   FORALLfacets {
01303     if (facet->flipped || !facet->normal)
01304       continue;
01305     totpart++;
01306     qh_distplane(point, facet, &dist);
01307     if (dist > *bestdist) {
01308       *bestdist= dist;
01309       bestfacet= facet;
01310       if (dist > qh MINoutside) {
01311         *isoutside= True;
01312         break;
01313       }
01314     }
01315   }
01316   *numpart= totpart;
01317   trace3((qh ferr, 3016, "qh_findfacet_all: f%d dist %2.2g isoutside %d totpart %d\n",
01318           getid_(bestfacet), *bestdist, *isoutside, totpart));
01319   return bestfacet;
01320 } /* findfacet_all */
01321 
01322 /*-<a                             href="qh-poly.htm#TOC"
01323   >-------------------------------</a><a name="findgood">-</a>
01324 
01325   qh_findgood( facetlist, goodhorizon )
01326     identify good facets for qh.PRINTgood
01327     if qh.GOODvertex>0
01328       facet includes point as vertex
01329       if !match, returns goodhorizon
01330       inactive if qh.MERGING
01331     if qh.GOODpoint
01332       facet is visible or coplanar (>0) or not visible (<0)
01333     if qh.GOODthreshold
01334       facet->normal matches threshold
01335     if !goodhorizon and !match,
01336       selects facet with closest angle
01337       sets GOODclosest
01338 
01339   returns:
01340     number of new, good facets found
01341     determines facet->good
01342     may update qh.GOODclosest
01343 
01344   notes:
01345     qh_findgood_all further reduces the good region
01346 
01347   design:
01348     count good facets
01349     mark good facets for qh.GOODpoint
01350     mark good facets for qh.GOODthreshold
01351     if necessary
01352       update qh.GOODclosest
01353 */
01354 int qh_findgood(facetT *facetlist, int goodhorizon) {
01355   facetT *facet, *bestfacet= NULL;
01356   realT angle, bestangle= REALmax, dist;
01357   int  numgood=0;
01358 
01359   FORALLfacet_(facetlist) {
01360     if (facet->good)
01361       numgood++;
01362   }
01363   if (qh GOODvertex>0 && !qh MERGING) {
01364     FORALLfacet_(facetlist) {
01365       if (!qh_isvertex(qh GOODvertexp, facet->vertices)) {
01366         facet->good= False;
01367         numgood--;
01368       }
01369     }
01370   }
01371   if (qh GOODpoint && numgood) {
01372     FORALLfacet_(facetlist) {
01373       if (facet->good && facet->normal) {
01374         zinc_(Zdistgood);
01375         qh_distplane(qh GOODpointp, facet, &dist);
01376         if ((qh GOODpoint > 0) ^ (dist > 0.0)) {
01377           facet->good= False;
01378           numgood--;
01379         }
01380       }
01381     }
01382   }
01383   if (qh GOODthreshold && (numgood || goodhorizon || qh GOODclosest)) {
01384     FORALLfacet_(facetlist) {
01385       if (facet->good && facet->normal) {
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 && (!goodhorizon || qh GOODclosest)) {
01397       if (qh GOODclosest) {
01398         if (qh GOODclosest->visible)
01399           qh GOODclosest= NULL;
01400         else {
01401           qh_inthresholds(qh GOODclosest->normal, &angle);
01402           if (angle < bestangle)
01403             bestfacet= qh GOODclosest;
01404         }
01405       }
01406       if (bestfacet && bestfacet != qh GOODclosest) {
01407         if (qh GOODclosest)
01408           qh GOODclosest->good= False;
01409         qh GOODclosest= bestfacet;
01410         bestfacet->good= True;
01411         numgood++;
01412         trace2((qh ferr, 2044, "qh_findgood: f%d is closest(%2.2g) to thresholds\n",
01413            bestfacet->id, bestangle));
01414         return numgood;
01415       }
01416     }else if (qh GOODclosest) { /* numgood > 0 */
01417       qh GOODclosest->good= False;
01418       qh GOODclosest= NULL;
01419     }
01420   }
01421   zadd_(Zgoodfacet, numgood);
01422   trace2((qh ferr, 2045, "qh_findgood: found %d good facets with %d good horizon\n",
01423                numgood, goodhorizon));
01424   if (!numgood && qh GOODvertex>0 && !qh MERGING)
01425     return goodhorizon;
01426   return numgood;
01427 } /* findgood */
01428 
01429 /*-<a                             href="qh-poly.htm#TOC"
01430   >-------------------------------</a><a name="findgood_all">-</a>
01431 
01432   qh_findgood_all( facetlist )
01433     apply other constraints for good facets (used by qh.PRINTgood)
01434     if qh.GOODvertex
01435       facet includes (>0) or doesn't include (<0) point as vertex
01436       if last good facet and ONLYgood, prints warning and continues
01437     if qh.SPLITthresholds
01438       facet->normal matches threshold, or if none, the closest one
01439     calls qh_findgood
01440     nop if good not used
01441 
01442   returns:
01443     clears facet->good if not good
01444     sets qh.num_good
01445 
01446   notes:
01447     this is like qh_findgood but more restrictive
01448 
01449   design:
01450     uses qh_findgood to mark good facets
01451     marks facets for qh.GOODvertex
01452     marks facets for qh.SPLITthreholds
01453 */
01454 void qh_findgood_all(facetT *facetlist) {
01455   facetT *facet, *bestfacet=NULL;
01456   realT angle, bestangle= REALmax;
01457   int  numgood=0, startgood;
01458 
01459   if (!qh GOODvertex && !qh GOODthreshold && !qh GOODpoint
01460   && !qh SPLITthresholds)
01461     return;
01462   if (!qh ONLYgood)
01463     qh_findgood(qh facet_list, 0);
01464   FORALLfacet_(facetlist) {
01465     if (facet->good)
01466       numgood++;
01467   }
01468   if (qh GOODvertex <0 || (qh GOODvertex > 0 && qh MERGING)) {
01469     FORALLfacet_(facetlist) {
01470       if (facet->good && ((qh GOODvertex > 0) ^ !!qh_isvertex(qh GOODvertexp, facet->vertices))) {
01471         if (!--numgood) {
01472           if (qh ONLYgood) {
01473             qh_fprintf(qh ferr, 7064, "qhull warning: good vertex p%d does not match last good facet f%d.  Ignored.\n",
01474                qh_pointid(qh GOODvertexp), facet->id);
01475             return;
01476           }else if (qh GOODvertex > 0)
01477             qh_fprintf(qh ferr, 7065, "qhull warning: point p%d is not a vertex('QV%d').\n",
01478                 qh GOODvertex-1, qh GOODvertex-1);
01479           else
01480             qh_fprintf(qh ferr, 7066, "qhull warning: point p%d is a vertex for every facet('QV-%d').\n",
01481                 -qh GOODvertex - 1, -qh GOODvertex - 1);
01482         }
01483         facet->good= False;
01484       }
01485     }
01486   }
01487   startgood= numgood;
01488   if (qh SPLITthresholds) {
01489     FORALLfacet_(facetlist) {
01490       if (facet->good) {
01491         if (!qh_inthresholds(facet->normal, &angle)) {
01492           facet->good= False;
01493           numgood--;
01494           if (angle < bestangle) {
01495             bestangle= angle;
01496             bestfacet= facet;
01497           }
01498         }
01499       }
01500     }
01501     if (!numgood && bestfacet) {
01502       bestfacet->good= True;
01503       numgood++;
01504       trace0((qh ferr, 23, "qh_findgood_all: f%d is closest(%2.2g) to thresholds\n",
01505            bestfacet->id, bestangle));
01506       return;
01507     }
01508   }
01509   qh num_good= numgood;
01510   trace0((qh ferr, 24, "qh_findgood_all: %d good facets remain out of %d facets\n",
01511         numgood, startgood));
01512 } /* findgood_all */
01513 
01514 /*-<a                             href="qh-poly.htm#TOC"
01515   >-------------------------------</a><a name="furthestnext">-</a>
01516 
01517   qh_furthestnext()
01518     set qh.facet_next to facet with furthest of all furthest points
01519     searches all facets on qh.facet_list
01520 
01521   notes:
01522     this may help avoid precision problems
01523 */
01524 void qh_furthestnext(void /* qh facet_list */) {
01525   facetT *facet, *bestfacet= NULL;
01526   realT dist, bestdist= -REALmax;
01527 
01528   FORALLfacets {
01529     if (facet->outsideset) {
01530 #if qh_COMPUTEfurthest
01531       pointT *furthest;
01532       furthest= (pointT*)qh_setlast(facet->outsideset);
01533       zinc_(Zcomputefurthest);
01534       qh_distplane(furthest, facet, &dist);
01535 #else
01536       dist= facet->furthestdist;
01537 #endif
01538       if (dist > bestdist) {
01539         bestfacet= facet;
01540         bestdist= dist;
01541       }
01542     }
01543   }
01544   if (bestfacet) {
01545     qh_removefacet(bestfacet);
01546     qh_prependfacet(bestfacet, &qh facet_next);
01547     trace1((qh ferr, 1029, "qh_furthestnext: made f%d next facet(dist %.2g)\n",
01548             bestfacet->id, bestdist));
01549   }
01550 } /* furthestnext */
01551 
01552 /*-<a                             href="qh-poly.htm#TOC"
01553   >-------------------------------</a><a name="furthestout">-</a>
01554 
01555   qh_furthestout( facet )
01556     make furthest outside point the last point of outsideset
01557 
01558   returns:
01559     updates facet->outsideset
01560     clears facet->notfurthest
01561     sets facet->furthestdist
01562 
01563   design:
01564     determine best point of outsideset
01565     make it the last point of outsideset
01566 */
01567 void qh_furthestout(facetT *facet) {
01568   pointT *point, **pointp, *bestpoint= NULL;
01569   realT dist, bestdist= -REALmax;
01570 
01571   FOREACHpoint_(facet->outsideset) {
01572     qh_distplane(point, facet, &dist);
01573     zinc_(Zcomputefurthest);
01574     if (dist > bestdist) {
01575       bestpoint= point;
01576       bestdist= dist;
01577     }
01578   }
01579   if (bestpoint) {
01580     qh_setdel(facet->outsideset, point);
01581     qh_setappend(&facet->outsideset, point);
01582 #if !qh_COMPUTEfurthest
01583     facet->furthestdist= bestdist;
01584 #endif
01585   }
01586   facet->notfurthest= False;
01587   trace3((qh ferr, 3017, "qh_furthestout: p%d is furthest outside point of f%d\n",
01588           qh_pointid(point), facet->id));
01589 } /* furthestout */
01590 
01591 
01592 /*-<a                             href="qh-qhull.htm#TOC"
01593   >-------------------------------</a><a name="infiniteloop">-</a>
01594 
01595   qh_infiniteloop( facet )
01596     report infinite loop error due to facet
01597 */
01598 void qh_infiniteloop(facetT *facet) {
01599 
01600   qh_fprintf(qh ferr, 6149, "qhull internal error (qh_infiniteloop): potential infinite loop detected\n");
01601   qh_errexit(qh_ERRqhull, facet, NULL);
01602 } /* qh_infiniteloop */
01603 
01604 /*-<a                             href="qh-poly.htm#TOC"
01605   >-------------------------------</a><a name="initbuild">-</a>
01606 
01607   qh_initbuild()
01608     initialize hull and outside sets with point array
01609     qh.FIRSTpoint/qh.NUMpoints is point array
01610     if qh.GOODpoint
01611       adds qh.GOODpoint to initial hull
01612 
01613   returns:
01614     qh_facetlist with initial hull
01615     points partioned into outside sets, coplanar sets, or inside
01616     initializes qh.GOODpointp, qh.GOODvertexp,
01617 
01618   design:
01619     initialize global variables used during qh_buildhull
01620     determine precision constants and points with max/min coordinate values
01621       if qh.SCALElast, scale last coordinate(for 'd')
01622     build initial simplex
01623     partition input points into facets of initial simplex
01624     set up lists
01625     if qh.ONLYgood
01626       check consistency
01627       add qh.GOODvertex if defined
01628 */
01629 void qh_initbuild( void) {
01630   setT *maxpoints, *vertices;
01631   facetT *facet;
01632   int i, numpart;
01633   realT dist;
01634   boolT isoutside;
01635 
01636   qh furthest_id= -1;
01637   qh lastreport= 0;
01638   qh facet_id= qh vertex_id= qh ridge_id= 0;
01639   qh visit_id= qh vertex_visit= 0;
01640   qh maxoutdone= False;
01641 
01642   if (qh GOODpoint > 0)
01643     qh GOODpointp= qh_point(qh GOODpoint-1);
01644   else if (qh GOODpoint < 0)
01645     qh GOODpointp= qh_point(-qh GOODpoint-1);
01646   if (qh GOODvertex > 0)
01647     qh GOODvertexp= qh_point(qh GOODvertex-1);
01648   else if (qh GOODvertex < 0)
01649     qh GOODvertexp= qh_point(-qh GOODvertex-1);
01650   if ((qh GOODpoint
01651        && (qh GOODpointp < qh first_point  /* also catches !GOODpointp */
01652            || qh GOODpointp > qh_point(qh num_points-1)))
01653     || (qh GOODvertex
01654         && (qh GOODvertexp < qh first_point  /* also catches !GOODvertexp */
01655             || qh GOODvertexp > qh_point(qh num_points-1)))) {
01656     qh_fprintf(qh ferr, 6150, "qhull input error: either QGn or QVn point is > p%d\n",
01657              qh num_points-1);
01658     qh_errexit(qh_ERRinput, NULL, NULL);
01659   }
01660   maxpoints= qh_maxmin(qh first_point, qh num_points, qh hull_dim);
01661   if (qh SCALElast)
01662     qh_scalelast(qh first_point, qh num_points, qh hull_dim,
01663                qh MINlastcoord, qh MAXlastcoord, qh MAXwidth);
01664   qh_detroundoff();
01665   if (qh DELAUNAY && qh upper_threshold[qh hull_dim-1] > REALmax/2
01666                   && qh lower_threshold[qh hull_dim-1] < -REALmax/2) {
01667     for (i=qh_PRINTEND; i--; ) {
01668       if (qh PRINTout[i] == qh_PRINTgeom && qh DROPdim < 0
01669           && !qh GOODthreshold && !qh SPLITthresholds)
01670         break;  /* in this case, don't set upper_threshold */
01671     }
01672     if (i < 0) {
01673       if (qh UPPERdelaunay) { /* matches qh.upperdelaunay in qh_setfacetplane */
01674         qh lower_threshold[qh hull_dim-1]= qh ANGLEround * qh_ZEROdelaunay;
01675         qh GOODthreshold= True;
01676       }else {
01677         qh upper_threshold[qh hull_dim-1]= -qh ANGLEround * qh_ZEROdelaunay;
01678         if (!qh GOODthreshold)
01679           qh SPLITthresholds= True; /* build upper-convex hull even if Qg */
01680           /* qh_initqhull_globals errors if Qg without Pdk/etc. */
01681       }
01682     }
01683   }
01684   vertices= qh_initialvertices(qh hull_dim, maxpoints, qh first_point, qh num_points);
01685   qh_initialhull(vertices);  /* initial qh facet_list */
01686   qh_partitionall(vertices, qh first_point, qh num_points);
01687   if (qh PRINToptions1st || qh TRACElevel || qh IStracing) {
01688     if (qh TRACElevel || qh IStracing)
01689       qh_fprintf(qh ferr, 8103, "\nTrace level %d for %s | %s\n",
01690          qh IStracing ? qh IStracing : qh TRACElevel, qh rbox_command, qh qhull_command);
01691     qh_fprintf(qh ferr, 8104, "Options selected for Qhull %s:\n%s\n", qh_version, qh qhull_options);
01692   }
01693   qh_resetlists(False, qh_RESETvisible /*qh visible_list newvertex_list newfacet_list */);
01694   qh facet_next= qh facet_list;
01695   qh_furthestnext(/* qh facet_list */);
01696   if (qh PREmerge) {
01697     qh cos_max= qh premerge_cos;
01698     qh centrum_radius= qh premerge_centrum;
01699   }
01700   if (qh ONLYgood) {
01701     if (qh GOODvertex > 0 && qh MERGING) {
01702       qh_fprintf(qh ferr, 6151, "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");
01703       qh_errexit(qh_ERRinput, NULL, NULL);
01704     }
01705     if (!(qh GOODthreshold || qh GOODpoint
01706          || (!qh MERGEexact && !qh PREmerge && qh GOODvertexp))) {
01707       qh_fprintf(qh ferr, 6152, "qhull input error: 'Qg' (ONLYgood) needs a good threshold('Pd0D0'), a\n\
01708 good point(QGn or QG-n), or a good vertex with 'QJ' or 'Q0' (QVn).\n");
01709       qh_errexit(qh_ERRinput, NULL, NULL);
01710     }
01711     if (qh GOODvertex > 0  && !qh MERGING  /* matches qh_partitionall */
01712         && !qh_isvertex(qh GOODvertexp, vertices)) {
01713       facet= qh_findbestnew(qh GOODvertexp, qh facet_list,
01714                           &dist, !qh_ALL, &isoutside, &numpart);
01715       zadd_(Zdistgood, numpart);
01716       if (!isoutside) {
01717         qh_fprintf(qh ferr, 6153, "qhull input error: point for QV%d is inside initial simplex.  It can not be made a vertex.\n",
01718                qh_pointid(qh GOODvertexp));
01719         qh_errexit(qh_ERRinput, NULL, NULL);
01720       }
01721       if (!qh_addpoint(qh GOODvertexp, facet, False)) {
01722         qh_settempfree(&vertices);
01723         qh_settempfree(&maxpoints);
01724         return;
01725       }
01726     }
01727     qh_findgood(qh facet_list, 0);
01728   }
01729   qh_settempfree(&vertices);
01730   qh_settempfree(&maxpoints);
01731   trace1((qh ferr, 1030, "qh_initbuild: initial hull created and points partitioned\n"));
01732 } /* initbuild */
01733 
01734 /*-<a                             href="qh-poly.htm#TOC"
01735   >-------------------------------</a><a name="initialhull">-</a>
01736 
01737   qh_initialhull( vertices )
01738     constructs the initial hull as a DIM3 simplex of vertices
01739 
01740   design:
01741     creates a simplex (initializes lists)
01742     determines orientation of simplex
01743     sets hyperplanes for facets
01744     doubles checks orientation (in case of axis-parallel facets with Gaussian elimination)
01745     checks for flipped facets and qh.NARROWhull
01746     checks the result
01747 */
01748 void qh_initialhull(setT *vertices) {
01749   facetT *facet, *firstfacet, *neighbor, **neighborp;
01750   realT dist, angle, minangle= REALmax;
01751 #ifndef qh_NOtrace
01752   int k;
01753 #endif
01754 
01755   qh_createsimplex(vertices);  /* qh facet_list */
01756   qh_resetlists(False, qh_RESETvisible);
01757   qh facet_next= qh facet_list;      /* advance facet when processed */
01758   qh interior_point= qh_getcenter(vertices);
01759   firstfacet= qh facet_list;
01760   qh_setfacetplane(firstfacet);
01761   zinc_(Znumvisibility); /* needs to be in printsummary */
01762   qh_distplane(qh interior_point, firstfacet, &dist);
01763   if (dist > 0) {
01764     FORALLfacets
01765       facet->toporient ^= (unsigned char)True;
01766   }
01767   FORALLfacets
01768     qh_setfacetplane(facet);
01769   FORALLfacets {
01770     if (!qh_checkflipped(facet, NULL, qh_ALL)) {/* due to axis-parallel facet */
01771       trace1((qh ferr, 1031, "qh_initialhull: initial orientation incorrect.  Correct all facets\n"));
01772       facet->flipped= False;
01773       FORALLfacets {
01774         facet->toporient ^= (unsigned char)True;
01775         qh_orientoutside(facet);
01776       }
01777       break;
01778     }
01779   }
01780   FORALLfacets {
01781     if (!qh_checkflipped(facet, NULL, !qh_ALL)) {  /* can happen with 'R0.1' */
01782       if (qh DELAUNAY && ! qh ATinfinity) {
01783         if (qh UPPERdelaunay)
01784           qh_fprintf(qh ferr, 6240, "Qhull input error: Can not compute the upper Delaunay triangulation or upper Voronoi diagram of cocircular/cospherical points.\n");
01785         else
01786           qh_fprintf(qh ferr, 6239, "Qhull input error: Use option 'Qz' for the Delaunay triangulation or Voronoi diagram of cocircular/cospherical points.  Option 'Qz' adds a point \"at infinity\" (above the corresponding paraboloid).\n");
01787         qh_errexit(qh_ERRinput, NULL, NULL);
01788       }
01789       qh_precision("initial facet is coplanar with interior point");
01790       qh_fprintf(qh ferr, 6154, "qhull precision error: initial facet %d is coplanar with the interior point\n",
01791                    facet->id);
01792       qh_errexit(qh_ERRsingular, facet, NULL);
01793     }
01794     FOREACHneighbor_(facet) {
01795       angle= qh_getangle(facet->normal, neighbor->normal);
01796       minimize_( minangle, angle);
01797     }
01798   }
01799   if (minangle < qh_MAXnarrow && !qh NOnarrow) {
01800     realT diff= 1.0 + minangle;
01801 
01802     qh NARROWhull= True;
01803     qh_option("_narrow-hull", NULL, &diff);
01804     if (minangle < qh_WARNnarrow && !qh RERUN && qh PRINTprecision)
01805       qh_printhelp_narrowhull(qh ferr, minangle);
01806   }
01807   zzval_(Zprocessed)= qh hull_dim+1;
01808   qh_checkpolygon(qh facet_list);
01809   qh_checkconvex(qh facet_list,   qh_DATAfault);
01810 #ifndef qh_NOtrace
01811   if (qh IStracing >= 1) {
01812     qh_fprintf(qh ferr, 8105, "qh_initialhull: simplex constructed, interior point:");
01813     for (k=0; k < qh hull_dim; k++)
01814       qh_fprintf(qh ferr, 8106, " %6.4g", qh interior_point[k]);
01815     qh_fprintf(qh ferr, 8107, "\n");
01816   }
01817 #endif
01818 } /* initialhull */
01819 
01820 /*-<a                             href="qh-poly.htm#TOC"
01821   >-------------------------------</a><a name="initialvertices">-</a>
01822 
01823   qh_initialvertices( dim, maxpoints, points, numpoints )
01824     determines a non-singular set of initial vertices
01825     maxpoints may include duplicate points
01826 
01827   returns:
01828     temporary set of dim+1 vertices in descending order by vertex id
01829     if qh.RANDOMoutside && !qh.ALLpoints
01830       picks random points
01831     if dim >= qh_INITIALmax,
01832       uses min/max x and max points with non-zero determinants
01833 
01834   notes:
01835     unless qh.ALLpoints,
01836       uses maxpoints as long as determinate is non-zero
01837 */
01838 setT *qh_initialvertices(int dim, setT *maxpoints, pointT *points, int numpoints) {
01839   pointT *point, **pointp;
01840   setT *vertices, *simplex, *tested;
01841   realT randr;
01842   int idx, point_i, point_n, k;
01843   boolT nearzero= False;
01844 
01845   vertices= qh_settemp(dim + 1);
01846   simplex= qh_settemp(dim+1);
01847   if (qh ALLpoints)
01848     qh_maxsimplex(dim, NULL, points, numpoints, &simplex);
01849   else if (qh RANDOMoutside) {
01850     while (qh_setsize(simplex) != dim+1) {
01851       randr= qh_RANDOMint;
01852       randr= randr/(qh_RANDOMmax+1);
01853       idx= (int)floor(qh num_points * randr);
01854       while (qh_setin(simplex, qh_point(idx))) {
01855             idx++; /* in case qh_RANDOMint always returns the same value */
01856         idx= idx < qh num_points ? idx : 0;
01857       }
01858       qh_setappend(&simplex, qh_point(idx));
01859     }
01860   }else if (qh hull_dim >= qh_INITIALmax) {
01861     tested= qh_settemp(dim+1);
01862     qh_setappend(&simplex, SETfirst_(maxpoints));   /* max and min X coord */
01863     qh_setappend(&simplex, SETsecond_(maxpoints));
01864     qh_maxsimplex(fmin_(qh_INITIALsearch, dim), maxpoints, points, numpoints, &simplex);
01865     k= qh_setsize(simplex);
01866     FOREACHpoint_i_(maxpoints) {
01867       if (point_i & 0x1) {     /* first pick up max. coord. points */
01868         if (!qh_setin(simplex, point) && !qh_setin(tested, point)){
01869           qh_detsimplex(point, simplex, k, &nearzero);
01870           if (nearzero)
01871             qh_setappend(&tested, point);
01872           else {
01873             qh_setappend(&simplex, point);
01874             if (++k == dim)  /* use search for last point */
01875               break;
01876           }
01877         }
01878       }
01879     }
01880     while (k != dim && (point= (pointT*)qh_setdellast(maxpoints))) {
01881       if (!qh_setin(simplex, point) && !qh_setin(tested, point)){
01882         qh_detsimplex(point, simplex, k, &nearzero);
01883         if (nearzero)
01884           qh_setappend(&tested, point);
01885         else {
01886           qh_setappend(&simplex, point);
01887           k++;
01888         }
01889       }
01890     }
01891     idx= 0;
01892     while (k != dim && (point= qh_point(idx++))) {
01893       if (!qh_setin(simplex, point) && !qh_setin(tested, point)){
01894         qh_detsimplex(point, simplex, k, &nearzero);
01895         if (!nearzero){
01896           qh_setappend(&simplex, point);
01897           k++;
01898         }
01899       }
01900     }
01901     qh_settempfree(&tested);
01902     qh_maxsimplex(dim, maxpoints, points, numpoints, &simplex);
01903   }else
01904     qh_maxsimplex(dim, maxpoints, points, numpoints, &simplex);
01905   FOREACHpoint_(simplex)
01906     qh_setaddnth(&vertices, 0, qh_newvertex(point)); /* descending order */
01907   qh_settempfree(&simplex);
01908   return vertices;
01909 } /* initialvertices */
01910 
01911 
01912 /*-<a                             href="qh-poly.htm#TOC"
01913   >-------------------------------</a><a name="isvertex">-</a>
01914 
01915   qh_isvertex(  )
01916     returns vertex if point is in vertex set, else returns NULL
01917 
01918   notes:
01919     for qh.GOODvertex
01920 */
01921 vertexT *qh_isvertex(pointT *point, setT *vertices) {
01922   vertexT *vertex, **vertexp;
01923 
01924   FOREACHvertex_(vertices) {
01925     if (vertex->point == point)
01926       return vertex;
01927   }
01928   return NULL;
01929 } /* isvertex */
01930 
01931 /*-<a                             href="qh-poly.htm#TOC"
01932   >-------------------------------</a><a name="makenewfacets">-</a>
01933 
01934   qh_makenewfacets( point )
01935     make new facets from point and qh.visible_list
01936 
01937   returns:
01938     qh.newfacet_list= list of new facets with hyperplanes and ->newfacet
01939     qh.newvertex_list= list of vertices in new facets with ->newlist set
01940 
01941     if (qh.ONLYgood)
01942       newfacets reference horizon facets, but not vice versa
01943       ridges reference non-simplicial horizon ridges, but not vice versa
01944       does not change existing facets
01945     else
01946       sets qh.NEWfacets
01947       new facets attached to horizon facets and ridges
01948       for visible facets,
01949         visible->r.replace is corresponding new facet
01950 
01951   see also:
01952     qh_makenewplanes() -- make hyperplanes for facets
01953     qh_attachnewfacets() -- attachnewfacets if not done here(qh ONLYgood)
01954     qh_matchnewfacets() -- match up neighbors
01955     qh_updatevertices() -- update vertex neighbors and delvertices
01956     qh_deletevisible() -- delete visible facets
01957     qh_checkpolygon() --check the result
01958     qh_triangulate() -- triangulate a non-simplicial facet
01959 
01960   design:
01961     for each visible facet
01962       make new facets to its horizon facets
01963       update its f.replace
01964       clear its neighbor set
01965 */
01966 vertexT *qh_makenewfacets(pointT *point /*visible_list*/) {
01967   facetT *visible, *newfacet= NULL, *newfacet2= NULL, *neighbor, **neighborp;
01968   vertexT *apex;
01969   int numnew=0;
01970 
01971   qh newfacet_list= qh facet_tail;
01972   qh newvertex_list= qh vertex_tail;
01973   apex= qh_newvertex(point);
01974   qh_appendvertex(apex);
01975   qh visit_id++;
01976   if (!qh ONLYgood)
01977     qh NEWfacets= True;
01978   FORALLvisible_facets {
01979     FOREACHneighbor_(visible)
01980       neighbor->seen= False;
01981     if (visible->ridges) {
01982       visible->visitid= qh visit_id;
01983       newfacet2= qh_makenew_nonsimplicial(visible, apex, &numnew);
01984     }
01985     if (visible->simplicial)
01986       newfacet= qh_makenew_simplicial(visible, apex, &numnew);
01987     if (!qh ONLYgood) {
01988       if (newfacet2)  /* newfacet is null if all ridges defined */
01989         newfacet= newfacet2;
01990       if (newfacet)
01991         visible->f.replace= newfacet;
01992       else
01993         zinc_(Zinsidevisible);
01994       SETfirst_(visible->neighbors)= NULL;
01995     }
01996   }
01997   trace1((qh ferr, 1032, "qh_makenewfacets: created %d new facets from point p%d to horizon\n",
01998           numnew, qh_pointid(point)));
01999   if (qh IStracing >= 4)
02000     qh_printfacetlist(qh newfacet_list, NULL, qh_ALL);
02001   return apex;
02002 } /* makenewfacets */
02003 
02004 /*-<a                             href="qh-poly.htm#TOC"
02005   >-------------------------------</a><a name="matchduplicates">-</a>
02006 
02007   qh_matchduplicates( atfacet, atskip, hashsize, hashcount )
02008     match duplicate ridges in qh.hash_table for atfacet/atskip
02009     duplicates marked with ->dupridge and qh_DUPLICATEridge
02010 
02011   returns:
02012     picks match with worst merge (min distance apart)
02013     updates hashcount
02014 
02015   see also:
02016     qh_matchneighbor
02017 
02018   notes:
02019 
02020   design:
02021     compute hash value for atfacet and atskip
02022     repeat twice -- once to make best matches, once to match the rest
02023       for each possible facet in qh.hash_table
02024         if it is a matching facet and pass 2
02025           make match
02026           unless tricoplanar, mark match for merging (qh_MERGEridge)
02027           [e.g., tricoplanar RBOX s 1000 t993602376 | QHULL C-1e-3 d Qbb FA Qt]
02028         if it is a matching facet and pass 1
02029           test if this is a better match
02030       if pass 1,
02031         make best match (it will not be merged)
02032 */
02033 #ifndef qh_NOmerge
02034 void qh_matchduplicates(facetT *atfacet, int atskip, int hashsize, int *hashcount) {
02035   boolT same, ismatch;
02036   int hash, scan;
02037   facetT *facet, *newfacet, *maxmatch= NULL, *maxmatch2= NULL, *nextfacet;
02038   int skip, newskip, nextskip= 0, maxskip= 0, maxskip2= 0, makematch;
02039   realT maxdist= -REALmax, mindist, dist2, low, high;
02040 
02041   hash= qh_gethash(hashsize, atfacet->vertices, qh hull_dim, 1,
02042                      SETelem_(atfacet->vertices, atskip));
02043   trace2((qh ferr, 2046, "qh_matchduplicates: find duplicate matches for f%d skip %d hash %d hashcount %d\n",
02044           atfacet->id, atskip, hash, *hashcount));
02045   for (makematch= 0; makematch < 2; makematch++) {
02046     qh visit_id++;
02047     for (newfacet= atfacet, newskip= atskip; newfacet; newfacet= nextfacet, newskip= nextskip) {
02048       zinc_(Zhashlookup);
02049       nextfacet= NULL;
02050       newfacet->visitid= qh visit_id;
02051       for (scan= hash; (facet= SETelemt_(qh hash_table, scan, facetT));
02052            scan= (++scan >= hashsize ? 0 : scan)) {
02053         if (!facet->dupridge || facet->visitid == qh visit_id)
02054           continue;
02055         zinc_(Zhashtests);
02056         if (qh_matchvertices(1, newfacet->vertices, newskip, facet->vertices, &skip, &same)) {
02057           ismatch= (same == (boolT)(newfacet->toporient ^ facet->toporient));
02058           if (SETelemt_(facet->neighbors, skip, facetT) != qh_DUPLICATEridge) {
02059             if (!makematch) {
02060               qh_fprintf(qh ferr, 6155, "qhull internal error (qh_matchduplicates): missing dupridge at f%d skip %d for new f%d skip %d hash %d\n",
02061                      facet->id, skip, newfacet->id, newskip, hash);
02062               qh_errexit2 (qh_ERRqhull, facet, newfacet);
02063             }
02064           }else if (ismatch && makematch) {
02065             if (SETelemt_(newfacet->neighbors, newskip, facetT) == qh_DUPLICATEridge) {
02066               SETelem_(facet->neighbors, skip)= newfacet;
02067               if (newfacet->tricoplanar)
02068                 SETelem_(newfacet->neighbors, newskip)= facet;
02069               else
02070                 SETelem_(newfacet->neighbors, newskip)= qh_MERGEridge;
02071               *hashcount -= 2; /* removed two unmatched facets */
02072               trace4((qh ferr, 4059, "qh_matchduplicates: duplicate f%d skip %d matched with new f%d skip %d merge\n",
02073                     facet->id, skip, newfacet->id, newskip));
02074             }
02075           }else if (ismatch) {
02076             mindist= qh_getdistance(facet, newfacet, &low, &high);
02077             dist2= qh_getdistance(newfacet, facet, &low, &high);
02078             minimize_(mindist, dist2);
02079             if (mindist > maxdist) {
02080               maxdist= mindist;
02081               maxmatch= facet;
02082               maxskip= skip;
02083               maxmatch2= newfacet;
02084               maxskip2= newskip;
02085             }
02086             trace3((qh ferr, 3018, "qh_matchduplicates: duplicate f%d skip %d new f%d skip %d at dist %2.2g, max is now f%d f%d\n",
02087                     facet->id, skip, newfacet->id, newskip, mindist,
02088                     maxmatch->id, maxmatch2->id));
02089           }else { /* !ismatch */
02090             nextfacet= facet;
02091             nextskip= skip;
02092           }
02093         }
02094         if (makematch && !facet
02095         && SETelemt_(facet->neighbors, skip, facetT) == qh_DUPLICATEridge) {
02096           qh_fprintf(qh ferr, 6156, "qhull internal error (qh_matchduplicates): no MERGEridge match for duplicate f%d skip %d at hash %d\n",
02097                      newfacet->id, newskip, hash);
02098           qh_errexit(qh_ERRqhull, newfacet, NULL);
02099         }
02100       }
02101     } /* end of for each new facet at hash */
02102     if (!makematch) {
02103       if (!maxmatch) {
02104         qh_fprintf(qh ferr, 6157, "qhull internal error (qh_matchduplicates): no maximum match at duplicate f%d skip %d at hash %d\n",
02105                      atfacet->id, atskip, hash);
02106         qh_errexit(qh_ERRqhull, atfacet, NULL);
02107       }
02108       SETelem_(maxmatch->neighbors, maxskip)= maxmatch2;
02109       SETelem_(maxmatch2->neighbors, maxskip2)= maxmatch;
02110       *hashcount -= 2; /* removed two unmatched facets */
02111       zzinc_(Zmultiridge);
02112       trace0((qh ferr, 25, "qh_matchduplicates: duplicate f%d skip %d matched with new f%d skip %d keep\n",
02113               maxmatch->id, maxskip, maxmatch2->id, maxskip2));
02114       qh_precision("ridge with multiple neighbors");
02115       if (qh IStracing >= 4)
02116         qh_errprint("DUPLICATED/MATCH", maxmatch, maxmatch2, NULL, NULL);
02117     }
02118   }
02119 } /* matchduplicates */
02120 
02121 /*-<a                             href="qh-poly.htm#TOC"
02122   >-------------------------------</a><a name="nearcoplanar">-</a>
02123 
02124   qh_nearcoplanar()
02125     for all facets, remove near-inside points from facet->coplanarset</li>
02126     coplanar points defined by innerplane from qh_outerinner()
02127 
02128   returns:
02129     if qh KEEPcoplanar && !qh KEEPinside
02130       facet->coplanarset only contains coplanar points
02131     if qh.JOGGLEmax
02132       drops inner plane by another qh.JOGGLEmax diagonal since a
02133         vertex could shift out while a coplanar point shifts in
02134 
02135   notes:
02136     used for qh.PREmerge and qh.JOGGLEmax
02137     must agree with computation of qh.NEARcoplanar in qh_detroundoff()
02138   design:
02139     if not keeping coplanar or inside points
02140       free all coplanar sets
02141     else if not keeping both coplanar and inside points
02142       remove !coplanar or !inside points from coplanar sets
02143 */
02144 void qh_nearcoplanar(void /* qh.facet_list */) {
02145   facetT *facet;
02146   pointT *point, **pointp;
02147   int numpart;
02148   realT dist, innerplane;
02149 
02150   if (!qh KEEPcoplanar && !qh KEEPinside) {
02151     FORALLfacets {
02152       if (facet->coplanarset)
02153         qh_setfree( &facet->coplanarset);
02154     }
02155   }else if (!qh KEEPcoplanar || !qh KEEPinside) {
02156     qh_outerinner(NULL, NULL, &innerplane);
02157     if (qh JOGGLEmax < REALmax/2)
02158       innerplane -= qh JOGGLEmax * sqrt((realT)qh hull_dim);
02159     numpart= 0;
02160     FORALLfacets {
02161       if (facet->coplanarset) {
02162         FOREACHpoint_(facet->coplanarset) {
02163           numpart++;
02164           qh_distplane(point, facet, &dist);
02165           if (dist < innerplane) {
02166             if (!qh KEEPinside)
02167               SETref_(point)= NULL;
02168           }else if (!qh KEEPcoplanar)
02169             SETref_(point)= NULL;
02170         }
02171         qh_setcompact(facet->coplanarset);
02172       }
02173     }
02174     zzadd_(Zcheckpart, numpart);
02175   }
02176 } /* nearcoplanar */
02177 
02178 /*-<a                             href="qh-poly.htm#TOC"
02179   >-------------------------------</a><a name="nearvertex">-</a>
02180 
02181   qh_nearvertex( facet, point, bestdist )
02182     return nearest vertex in facet to point
02183 
02184   returns:
02185     vertex and its distance
02186 
02187   notes:
02188     if qh.DELAUNAY
02189       distance is measured in the input set
02190     searches neighboring tricoplanar facets (requires vertexneighbors)
02191       Slow implementation.  Recomputes vertex set for each point.
02192     The vertex set could be stored in the qh.keepcentrum facet.
02193 */
02194 vertexT *qh_nearvertex(facetT *facet, pointT *point, realT *bestdistp) {
02195   realT bestdist= REALmax, dist;
02196   vertexT *bestvertex= NULL, *vertex, **vertexp, *apex;
02197   coordT *center;
02198   facetT *neighbor, **neighborp;
02199   setT *vertices;
02200   int dim= qh hull_dim;
02201 
02202   if (qh DELAUNAY)
02203     dim--;
02204   if (facet->tricoplanar) {
02205     if (!qh VERTEXneighbors || !facet->center) {
02206       qh_fprintf(qh ferr, 6158, "qhull internal error (qh_nearvertex): qh.VERTEXneighbors and facet->center required for tricoplanar facets\n");
02207       qh_errexit(qh_ERRqhull, facet, NULL);
02208     }
02209     vertices= qh_settemp(qh TEMPsize);
02210     apex= SETfirstt_(facet->vertices, vertexT);
02211     center= facet->center;
02212     FOREACHneighbor_(apex) {
02213       if (neighbor->center == center) {
02214         FOREACHvertex_(neighbor->vertices)
02215           qh_setappend(&vertices, vertex);
02216       }
02217     }
02218   }else
02219     vertices= facet->vertices;
02220   FOREACHvertex_(vertices) {
02221     dist= qh_pointdist(vertex->point, point, -dim);
02222     if (dist < bestdist) {
02223       bestdist= dist;
02224       bestvertex= vertex;
02225     }
02226   }
02227   if (facet->tricoplanar)
02228     qh_settempfree(&vertices);
02229   *bestdistp= sqrt(bestdist);
02230   trace3((qh ferr, 3019, "qh_nearvertex: v%d dist %2.2g for f%d p%d\n",
02231         bestvertex->id, *bestdistp, facet->id, qh_pointid(point)));
02232   return bestvertex;
02233 } /* nearvertex */
02234 
02235 /*-<a                             href="qh-poly.htm#TOC"
02236   >-------------------------------</a><a name="newhashtable">-</a>
02237 
02238   qh_newhashtable( newsize )
02239     returns size of qh.hash_table of at least newsize slots
02240 
02241   notes:
02242     assumes qh.hash_table is NULL
02243     qh_HASHfactor determines the number of extra slots
02244     size is not divisible by 2, 3, or 5
02245 */
02246 int qh_newhashtable(int newsize) {
02247   int size;
02248 
02249   size= ((newsize+1)*qh_HASHfactor) | 0x1;  /* odd number */
02250   while (True) {
02251     if (newsize<0 || size<0) {
02252         qh_fprintf(qhmem.ferr, 6236, "qhull error (qh_newhashtable): negative request (%d) or size (%d).  Did int overflow due to high-D?\n", newsize, size); /* WARN64 */
02253         qh_errexit(qhmem_ERRmem, NULL, NULL);
02254     }
02255     if ((size%3) && (size%5))
02256       break;
02257     size += 2;
02258     /* loop terminates because there is an infinite number of primes */
02259   }
02260   qh hash_table= qh_setnew(size);
02261   qh_setzero(qh hash_table, 0, size);
02262   return size;
02263 } /* newhashtable */
02264 
02265 /*-<a                             href="qh-poly.htm#TOC"
02266   >-------------------------------</a><a name="newvertex">-</a>
02267 
02268   qh_newvertex( point )
02269     returns a new vertex for point
02270 */
02271 vertexT *qh_newvertex(pointT *point) {
02272   vertexT *vertex;
02273 
02274   zinc_(Ztotvertices);
02275   vertex= (vertexT *)qh_memalloc((int)sizeof(vertexT));
02276   memset((char *) vertex, (size_t)0, sizeof(vertexT));
02277   if (qh vertex_id == 0xFFFFFF) {
02278     qh_fprintf(qh ferr, 6159, "qhull error: more than %d vertices.  ID field overflows and two vertices\n\
02279 may have the same identifier.  Vertices will not be sorted correctly.\n", 0xFFFFFF);
02280     qh_errexit(qh_ERRqhull, NULL, NULL);
02281   }
02282   if (qh vertex_id == qh tracevertex_id)
02283     qh tracevertex= vertex;
02284   vertex->id= qh vertex_id++;
02285   vertex->point= point;
02286   vertex->dim= (unsigned char)(qh hull_dim <= MAX_vdim ? qh hull_dim : 0);
02287   trace4((qh ferr, 4060, "qh_newvertex: vertex p%d(v%d) created\n", qh_pointid(vertex->point),
02288           vertex->id));
02289   return(vertex);
02290 } /* newvertex */
02291 
02292 /*-<a                             href="qh-poly.htm#TOC"
02293   >-------------------------------</a><a name="nextridge3d">-</a>
02294 
02295   qh_nextridge3d( atridge, facet, vertex )
02296     return next ridge and vertex for a 3d facet
02297     returns NULL on error
02298     [for QhullFacet::nextRidge3d] Does not call qh_errexit nor access qh_qh.
02299 
02300   notes:
02301     in qh_ORIENTclock order
02302     this is a O(n^2) implementation to trace all ridges
02303     be sure to stop on any 2nd visit
02304     same as QhullRidge::nextRidge3d
02305     does not use qh_qh or qh_errexit [QhullFacet.cpp]
02306 
02307   design:
02308     for each ridge
02309       exit if it is the ridge after atridge
02310 */
02311 ridgeT *qh_nextridge3d(ridgeT *atridge, facetT *facet, vertexT **vertexp) {
02312   vertexT *atvertex, *vertex, *othervertex;
02313   ridgeT *ridge, **ridgep;
02314 
02315   if ((atridge->top == facet) ^ qh_ORIENTclock)
02316     atvertex= SETsecondt_(atridge->vertices, vertexT);
02317   else
02318     atvertex= SETfirstt_(atridge->vertices, vertexT);
02319   FOREACHridge_(facet->ridges) {
02320     if (ridge == atridge)
02321       continue;
02322     if ((ridge->top == facet) ^ qh_ORIENTclock) {
02323       othervertex= SETsecondt_(ridge->vertices, vertexT);
02324       vertex= SETfirstt_(ridge->vertices, vertexT);
02325     }else {
02326       vertex= SETsecondt_(ridge->vertices, vertexT);
02327       othervertex= SETfirstt_(ridge->vertices, vertexT);
02328     }
02329     if (vertex == atvertex) {
02330       if (vertexp)
02331         *vertexp= othervertex;
02332       return ridge;
02333     }
02334   }
02335   return NULL;
02336 } /* nextridge3d */
02337 #else /* qh_NOmerge */
02338 void qh_matchduplicates(facetT *atfacet, int atskip, int hashsize, int *hashcount) {
02339 }
02340 ridgeT *qh_nextridge3d(ridgeT *atridge, facetT *facet, vertexT **vertexp) {
02341 
02342   return NULL;
02343 }
02344 #endif /* qh_NOmerge */
02345 
02346 /*-<a                             href="qh-poly.htm#TOC"
02347   >-------------------------------</a><a name="outcoplanar">-</a>
02348 
02349   qh_outcoplanar()
02350     move points from all facets' outsidesets to their coplanarsets
02351 
02352   notes:
02353     for post-processing under qh.NARROWhull
02354 
02355   design:
02356     for each facet
02357       for each outside point for facet
02358         partition point into coplanar set
02359 */
02360 void qh_outcoplanar(void /* facet_list */) {
02361   pointT *point, **pointp;
02362   facetT *facet;
02363   realT dist;
02364 
02365   trace1((qh ferr, 1033, "qh_outcoplanar: move outsideset to coplanarset for qh NARROWhull\n"));
02366   FORALLfacets {
02367     FOREACHpoint_(facet->outsideset) {
02368       qh num_outside--;
02369       if (qh KEEPcoplanar || qh KEEPnearinside) {
02370         qh_distplane(point, facet, &dist);
02371         zinc_(Zpartition);
02372         qh_partitioncoplanar(point, facet, &dist);
02373       }
02374     }
02375     qh_setfree(&facet->outsideset);
02376   }
02377 } /* outcoplanar */
02378 
02379 /*-<a                             href="qh-poly.htm#TOC"
02380   >-------------------------------</a><a name="point">-</a>
02381 
02382   qh_point( id )
02383     return point for a point id, or NULL if unknown
02384 
02385   alternative code:
02386     return((pointT *)((unsigned   long)qh.first_point
02387            + (unsigned long)((id)*qh.normal_size)));
02388 */
02389 pointT *qh_point(int id) {
02390 
02391   if (id < 0)
02392     return NULL;
02393   if (id < qh num_points)
02394     return qh first_point + id * qh hull_dim;
02395   id -= qh num_points;
02396   if (id < qh_setsize(qh other_points))
02397     return SETelemt_(qh other_points, id, pointT);
02398   return NULL;
02399 } /* point */
02400 
02401 /*-<a                             href="qh-poly.htm#TOC"
02402   >-------------------------------</a><a name="point_add">-</a>
02403 
02404   qh_point_add( set, point, elem )
02405     stores elem at set[point.id]
02406 
02407   returns:
02408     access function for qh_pointfacet and qh_pointvertex
02409 
02410   notes:
02411     checks point.id
02412 */
02413 void qh_point_add(setT *set, pointT *point, void *elem) {
02414   int id, size;
02415 
02416   SETreturnsize_(set, size);
02417   if ((id= qh_pointid(point)) < 0)
02418     qh_fprintf(qh ferr, 7067, "qhull internal warning (point_add): unknown point %p id %d\n",
02419       point, id);
02420   else if (id >= size) {
02421     qh_fprintf(qh ferr, 6160, "qhull internal errror(point_add): point p%d is out of bounds(%d)\n",
02422              id, size);
02423     qh_errexit(qh_ERRqhull, NULL, NULL);
02424   }else
02425     SETelem_(set, id)= elem;
02426 } /* point_add */
02427 
02428 
02429 /*-<a                             href="qh-poly.htm#TOC"
02430   >-------------------------------</a><a name="pointfacet">-</a>
02431 
02432   qh_pointfacet()
02433     return temporary set of facet for each point
02434     the set is indexed by point id
02435 
02436   notes:
02437     vertices assigned to one of the facets
02438     coplanarset assigned to the facet
02439     outside set assigned to the facet
02440     NULL if no facet for point (inside)
02441       includes qh.GOODpointp
02442 
02443   access:
02444     FOREACHfacet_i_(facets) { ... }
02445     SETelem_(facets, i)
02446 
02447   design:
02448     for each facet
02449       add each vertex
02450       add each coplanar point
02451       add each outside point
02452 */
02453 setT *qh_pointfacet(void /*qh facet_list*/) {
02454   int numpoints= qh num_points + qh_setsize(qh other_points);
02455   setT *facets;
02456   facetT *facet;
02457   vertexT *vertex, **vertexp;
02458   pointT *point, **pointp;
02459 
02460   facets= qh_settemp(numpoints);
02461   qh_setzero(facets, 0, numpoints);
02462   qh vertex_visit++;
02463   FORALLfacets {
02464     FOREACHvertex_(facet->vertices) {
02465       if (vertex->visitid != qh vertex_visit) {
02466         vertex->visitid= qh vertex_visit;
02467         qh_point_add(facets, vertex->point, facet);
02468       }
02469     }
02470     FOREACHpoint_(facet->coplanarset)
02471       qh_point_add(facets, point, facet);
02472     FOREACHpoint_(facet->outsideset)
02473       qh_point_add(facets, point, facet);
02474   }
02475   return facets;
02476 } /* pointfacet */
02477 
02478 /*-<a                             href="qh-poly.htm#TOC"
02479   >-------------------------------</a><a name="pointvertex">-</a>
02480 
02481   qh_pointvertex(  )
02482     return temporary set of vertices indexed by point id
02483     entry is NULL if no vertex for a point
02484       this will include qh.GOODpointp
02485 
02486   access:
02487     FOREACHvertex_i_(vertices) { ... }
02488     SETelem_(vertices, i)
02489 */
02490 setT *qh_pointvertex(void /*qh facet_list*/) {
02491   int numpoints= qh num_points + qh_setsize(qh other_points);
02492   setT *vertices;
02493   vertexT *vertex;
02494 
02495   vertices= qh_settemp(numpoints);
02496   qh_setzero(vertices, 0, numpoints);
02497   FORALLvertices
02498     qh_point_add(vertices, vertex->point, vertex);
02499   return vertices;
02500 } /* pointvertex */
02501 
02502 
02503 /*-<a                             href="qh-poly.htm#TOC"
02504   >-------------------------------</a><a name="prependfacet">-</a>
02505 
02506   qh_prependfacet( facet, facetlist )
02507     prepend facet to the start of a facetlist
02508 
02509   returns:
02510     increments qh.numfacets
02511     updates facetlist, qh.facet_list, facet_next
02512 
02513   notes:
02514     be careful of prepending since it can lose a pointer.
02515       e.g., can lose _next by deleting and then prepending before _next
02516 */
02517 void qh_prependfacet(facetT *facet, facetT **facetlist) {
02518   facetT *prevfacet, *list;
02519 
02520 
02521   trace4((qh ferr, 4061, "qh_prependfacet: prepend f%d before f%d\n",
02522           facet->id, getid_(*facetlist)));
02523   if (!*facetlist)
02524     (*facetlist)= qh facet_tail;
02525   list= *facetlist;
02526   prevfacet= list->previous;
02527   facet->previous= prevfacet;
02528   if (prevfacet)
02529     prevfacet->next= facet;
02530   list->previous= facet;
02531   facet->next= *facetlist;
02532   if (qh facet_list == list)  /* this may change *facetlist */
02533     qh facet_list= facet;
02534   if (qh facet_next == list)
02535     qh facet_next= facet;
02536   *facetlist= facet;
02537   qh num_facets++;
02538 } /* prependfacet */
02539 
02540 
02541 /*-<a                             href="qh-poly.htm#TOC"
02542   >-------------------------------</a><a name="printhashtable">-</a>
02543 
02544   qh_printhashtable( fp )
02545     print hash table to fp
02546 
02547   notes:
02548     not in I/O to avoid bringing io.c in
02549 
02550   design:
02551     for each hash entry
02552       if defined
02553         if unmatched or will merge (NULL, qh_MERGEridge, qh_DUPLICATEridge)
02554           print entry and neighbors
02555 */
02556 void qh_printhashtable(FILE *fp) {
02557   facetT *facet, *neighbor;
02558   int id, facet_i, facet_n, neighbor_i= 0, neighbor_n= 0;
02559   vertexT *vertex, **vertexp;
02560 
02561   FOREACHfacet_i_(qh hash_table) {
02562     if (facet) {
02563       FOREACHneighbor_i_(facet) {
02564         if (!neighbor || neighbor == qh_MERGEridge || neighbor == qh_DUPLICATEridge)
02565           break;
02566       }
02567       if (neighbor_i == neighbor_n)
02568         continue;
02569       qh_fprintf(fp, 9283, "hash %d f%d ", facet_i, facet->id);
02570       FOREACHvertex_(facet->vertices)
02571         qh_fprintf(fp, 9284, "v%d ", vertex->id);
02572       qh_fprintf(fp, 9285, "\n neighbors:");
02573       FOREACHneighbor_i_(facet) {
02574         if (neighbor == qh_MERGEridge)
02575           id= -3;
02576         else if (neighbor == qh_DUPLICATEridge)
02577           id= -2;
02578         else
02579           id= getid_(neighbor);
02580         qh_fprintf(fp, 9286, " %d", id);
02581       }
02582       qh_fprintf(fp, 9287, "\n");
02583     }
02584   }
02585 } /* printhashtable */
02586 
02587 
02588 /*-<a                             href="qh-poly.htm#TOC"
02589   >-------------------------------</a><a name="printlists">-</a>
02590 
02591   qh_printlists( fp )
02592     print out facet and vertex list for debugging (without 'f/v' tags)
02593 */
02594 void qh_printlists(void) {
02595   facetT *facet;
02596   vertexT *vertex;
02597   int count= 0;
02598 
02599   qh_fprintf(qh ferr, 8108, "qh_printlists: facets:");
02600   FORALLfacets {
02601     if (++count % 100 == 0)
02602       qh_fprintf(qh ferr, 8109, "\n     ");
02603     qh_fprintf(qh ferr, 8110, " %d", facet->id);
02604   }
02605   qh_fprintf(qh ferr, 8111, "\n  new facets %d visible facets %d next facet for qh_addpoint %d\n  vertices(new %d):",
02606      getid_(qh newfacet_list), getid_(qh visible_list), getid_(qh facet_next),
02607      getid_(qh newvertex_list));
02608   count = 0;
02609   FORALLvertices {
02610     if (++count % 100 == 0)
02611       qh_fprintf(qh ferr, 8112, "\n     ");
02612     qh_fprintf(qh ferr, 8113, " %d", vertex->id);
02613   }
02614   qh_fprintf(qh ferr, 8114, "\n");
02615 } /* printlists */
02616 
02617 /*-<a                             href="qh-poly.htm#TOC"
02618   >-------------------------------</a><a name="resetlists">-</a>
02619 
02620   qh_resetlists( stats, qh_RESETvisible )
02621     reset newvertex_list, newfacet_list, visible_list
02622     if stats,
02623       maintains statistics
02624 
02625   returns:
02626     visible_list is empty if qh_deletevisible was called
02627 */
02628 void qh_resetlists(boolT stats, boolT resetVisible /*qh newvertex_list newfacet_list visible_list*/) {
02629   vertexT *vertex;
02630   facetT *newfacet, *visible;
02631   int totnew=0, totver=0;
02632 
02633   if (stats) {
02634     FORALLvertex_(qh newvertex_list)
02635       totver++;
02636     FORALLnew_facets
02637       totnew++;
02638     zadd_(Zvisvertextot, totver);
02639     zmax_(Zvisvertexmax, totver);
02640     zadd_(Znewfacettot, totnew);
02641     zmax_(Znewfacetmax, totnew);
02642   }
02643   FORALLvertex_(qh newvertex_list)
02644     vertex->newlist= False;
02645   qh newvertex_list= NULL;
02646   FORALLnew_facets
02647     newfacet->newfacet= False;
02648   qh newfacet_list= NULL;
02649   if (resetVisible) {
02650     FORALLvisible_facets {
02651       visible->f.replace= NULL;
02652       visible->visible= False;
02653     }
02654     qh num_visible= 0;
02655   }
02656   qh visible_list= NULL; /* may still have visible facets via qh_triangulate */
02657   qh NEWfacets= False;
02658 } /* resetlists */
02659 
02660 /*-<a                             href="qh-poly.htm#TOC"
02661   >-------------------------------</a><a name="setvoronoi_all">-</a>
02662 
02663   qh_setvoronoi_all()
02664     compute Voronoi centers for all facets
02665     includes upperDelaunay facets if qh.UPPERdelaunay ('Qu')
02666 
02667   returns:
02668     facet->center is the Voronoi center
02669 
02670   notes:
02671     this is unused/untested code
02672       please email bradb@shore.net if this works ok for you
02673 
02674   use:
02675     FORALLvertices {...} to locate the vertex for a point.
02676     FOREACHneighbor_(vertex) {...} to visit the Voronoi centers for a Voronoi cell.
02677 */
02678 void qh_setvoronoi_all(void) {
02679   facetT *facet;
02680 
02681   qh_clearcenters(qh_ASvoronoi);
02682   qh_vertexneighbors();
02683 
02684   FORALLfacets {
02685     if (!facet->normal || !facet->upperdelaunay || qh UPPERdelaunay) {
02686       if (!facet->center)
02687         facet->center= qh_facetcenter(facet->vertices);
02688     }
02689   }
02690 } /* setvoronoi_all */
02691 
02692 #ifndef qh_NOmerge
02693 
02694 /*-<a                             href="qh-poly.htm#TOC"
02695   >-------------------------------</a><a name="triangulate">-</a>
02696 
02697   qh_triangulate()
02698     triangulate non-simplicial facets on qh.facet_list,
02699     if qh VORONOI, sets Voronoi centers of non-simplicial facets
02700     nop if hasTriangulation
02701 
02702   returns:
02703     all facets simplicial
02704     each tricoplanar facet has ->f.triowner == owner of ->center,normal,etc.
02705 
02706   notes:
02707     call after qh_check_output since may switch to Voronoi centers
02708     Output may overwrite ->f.triowner with ->f.area
02709 */
02710 void qh_triangulate(void /*qh facet_list*/) {
02711   facetT *facet, *nextfacet, *owner;
02712   int onlygood= qh ONLYgood;
02713   facetT *neighbor, *visible= NULL, *facet1, *facet2, *new_facet_list= NULL;
02714   facetT *orig_neighbor= NULL, *otherfacet;
02715   vertexT *new_vertex_list= NULL;
02716   mergeT *merge;
02717   mergeType mergetype;
02718   int neighbor_i, neighbor_n;
02719 
02720   if (qh hasTriangulation)
02721       return;
02722   trace1((qh ferr, 1034, "qh_triangulate: triangulate non-simplicial facets\n"));
02723   if (qh hull_dim == 2)
02724     return;
02725   if (qh VORONOI) {  /* otherwise lose Voronoi centers [could rebuild vertex set from tricoplanar] */
02726     qh_clearcenters(qh_ASvoronoi);
02727     qh_vertexneighbors();
02728   }
02729   qh ONLYgood= False; /* for makenew_nonsimplicial */
02730   qh visit_id++;
02731   qh NEWfacets= True;
02732   qh degen_mergeset= qh_settemp(qh TEMPsize);
02733   qh newvertex_list= qh vertex_tail;
02734   for (facet= qh facet_list; facet && facet->next; facet= nextfacet) { /* non-simplicial facets moved to end */
02735     nextfacet= facet->next;
02736     if (facet->visible || facet->simplicial)
02737       continue;
02738     /* triangulate all non-simplicial facets, otherwise merging does not work, e.g., RBOX c P-0.1 P+0.1 P+0.1 D3 | QHULL d Qt Tv */
02739     if (!new_facet_list)
02740       new_facet_list= facet;  /* will be moved to end */
02741     qh_triangulate_facet(facet, &new_vertex_list);
02742   }
02743   trace2((qh ferr, 2047, "qh_triangulate: delete null facets from f%d -- apex same as second vertex\n", getid_(new_facet_list)));
02744   for (facet= new_facet_list; facet && facet->next; facet= nextfacet) { /* null facets moved to end */
02745     nextfacet= facet->next;
02746     if (facet->visible)
02747       continue;
02748     if (facet->ridges) {
02749       if (qh_setsize(facet->ridges) > 0) {
02750         qh_fprintf(qh ferr, 6161, "qhull error (qh_triangulate): ridges still defined for f%d\n", facet->id);
02751         qh_errexit(qh_ERRqhull, facet, NULL);
02752       }
02753       qh_setfree(&facet->ridges);
02754     }
02755     if (SETfirst_(facet->vertices) == SETsecond_(facet->vertices)) {
02756       zinc_(Ztrinull);
02757       qh_triangulate_null(facet);
02758     }
02759   }
02760   trace2((qh ferr, 2048, "qh_triangulate: delete %d or more mirror facets -- same vertices and neighbors\n", qh_setsize(qh degen_mergeset)));
02761   qh visible_list= qh facet_tail;
02762   while ((merge= (mergeT*)qh_setdellast(qh degen_mergeset))) {
02763     facet1= merge->facet1;
02764     facet2= merge->facet2;
02765     mergetype= merge->type;
02766     qh_memfree(merge, (int)sizeof(mergeT));
02767     if (mergetype == MRGmirror) {
02768       zinc_(Ztrimirror);
02769       qh_triangulate_mirror(facet1, facet2);
02770     }
02771   }
02772   qh_settempfree(&qh degen_mergeset);
02773   trace2((qh ferr, 2049, "qh_triangulate: update neighbor lists for vertices from v%d\n", getid_(new_vertex_list)));
02774   qh newvertex_list= new_vertex_list;  /* all vertices of new facets */
02775   qh visible_list= NULL;
02776   qh_updatevertices(/*qh newvertex_list, empty newfacet_list and visible_list*/);
02777   qh_resetlists(False, !qh_RESETvisible /*qh newvertex_list, empty newfacet_list and visible_list*/);
02778 
02779   trace2((qh ferr, 2050, "qh_triangulate: identify degenerate tricoplanar facets from f%d\n", getid_(new_facet_list)));
02780   trace2((qh ferr, 2051, "qh_triangulate: and replace facet->f.triowner with tricoplanar facets that own center, normal, etc.\n"));
02781   FORALLfacet_(new_facet_list) {
02782     if (facet->tricoplanar && !facet->visible) {
02783       FOREACHneighbor_i_(facet) {
02784         if (neighbor_i == 0) {  /* first iteration */
02785           if (neighbor->tricoplanar)
02786             orig_neighbor= neighbor->f.triowner;
02787           else
02788             orig_neighbor= neighbor;
02789         }else {
02790           if (neighbor->tricoplanar)
02791             otherfacet= neighbor->f.triowner;
02792           else
02793             otherfacet= neighbor;
02794           if (orig_neighbor == otherfacet) {
02795             zinc_(Ztridegen);
02796             facet->degenerate= True;
02797             break;
02798           }
02799         }
02800       }
02801     }
02802   }
02803 
02804   trace2((qh ferr, 2052, "qh_triangulate: delete visible facets -- non-simplicial, null, and mirrored facets\n"));
02805   owner= NULL;
02806   visible= NULL;
02807   for (facet= new_facet_list; facet && facet->next; facet= nextfacet) { /* may delete facet */
02808     nextfacet= facet->next;
02809     if (facet->visible) {
02810       if (facet->tricoplanar) { /* a null or mirrored facet */
02811         qh_delfacet(facet);
02812         qh num_visible--;
02813       }else {  /* a non-simplicial facet followed by its tricoplanars */
02814         if (visible && !owner) {
02815           /*  RBOX 200 s D5 t1001471447 | QHULL Qt C-0.01 Qx Qc Tv Qt -- f4483 had 6 vertices/neighbors and 8 ridges */
02816           trace2((qh ferr, 2053, "qh_triangulate: all tricoplanar facets degenerate for non-simplicial facet f%d\n",
02817                        visible->id));
02818           qh_delfacet(visible);
02819           qh num_visible--;
02820         }
02821         visible= facet;
02822         owner= NULL;
02823       }
02824     }else if (facet->tricoplanar) {
02825       if (facet->f.triowner != visible) {
02826         qh_fprintf(qh ferr, 6162, "qhull error (qh_triangulate): tricoplanar facet f%d not owned by its visible, non-simplicial facet f%d\n", facet->id, getid_(visible));
02827         qh_errexit2 (qh_ERRqhull, facet, visible);
02828       }
02829       if (owner)
02830         facet->f.triowner= owner;
02831       else if (!facet->degenerate) {
02832         owner= facet;
02833         nextfacet= visible->next; /* rescan tricoplanar facets with owner */
02834         facet->keepcentrum= True;  /* one facet owns ->normal, etc. */
02835         facet->coplanarset= visible->coplanarset;
02836         facet->outsideset= visible->outsideset;
02837         visible->coplanarset= NULL;
02838         visible->outsideset= NULL;
02839         if (!qh TRInormals) { /* center and normal copied to tricoplanar facets */
02840           visible->center= NULL;
02841           visible->normal= NULL;
02842         }
02843         qh_delfacet(visible);
02844         qh num_visible--;
02845       }
02846     }
02847   }
02848   if (visible && !owner) {
02849     trace2((qh ferr, 2054, "qh_triangulate: all tricoplanar facets degenerate for last non-simplicial facet f%d\n",
02850                  visible->id));
02851     qh_delfacet(visible);
02852     qh num_visible--;
02853   }
02854   qh NEWfacets= False;
02855   qh ONLYgood= onlygood; /* restore value */
02856   if (qh CHECKfrequently)
02857     qh_checkpolygon(qh facet_list);
02858   qh hasTriangulation= True;
02859 } /* triangulate */
02860 
02861 
02862 /*-<a                             href="qh-poly.htm#TOC"
02863   >-------------------------------</a><a name="triangulate_facet">-</a>
02864 
02865   qh_triangulate_facet(facetA)
02866     triangulate a non-simplicial facet
02867       if qh.CENTERtype=qh_ASvoronoi, sets its Voronoi center
02868   returns:
02869     qh.newfacet_list == simplicial facets
02870       facet->tricoplanar set and ->keepcentrum false
02871       facet->degenerate set if duplicated apex
02872       facet->f.trivisible set to facetA
02873       facet->center copied from facetA (created if qh_ASvoronoi)
02874         qh_eachvoronoi, qh_detvridge, qh_detvridge3 assume centers copied
02875       facet->normal,offset,maxoutside copied from facetA
02876 
02877   notes:
02878       qh_makenew_nonsimplicial uses neighbor->seen for the same
02879 
02880   see also:
02881       qh_addpoint() -- add a point
02882       qh_makenewfacets() -- construct a cone of facets for a new vertex
02883 
02884   design:
02885       if qh_ASvoronoi,
02886          compute Voronoi center (facet->center)
02887       select first vertex (highest ID to preserve ID ordering of ->vertices)
02888       triangulate from vertex to ridges
02889       copy facet->center, normal, offset
02890       update vertex neighbors
02891 */
02892 void qh_triangulate_facet(facetT *facetA, vertexT **first_vertex) {
02893   facetT *newfacet;
02894   facetT *neighbor, **neighborp;
02895   vertexT *apex;
02896   int numnew=0;
02897 
02898   trace3((qh ferr, 3020, "qh_triangulate_facet: triangulate facet f%d\n", facetA->id));
02899 
02900   if (qh IStracing >= 4)
02901     qh_printfacet(qh ferr, facetA);
02902   FOREACHneighbor_(facetA) {
02903     neighbor->seen= False;
02904     neighbor->coplanar= False;
02905   }
02906   if (qh CENTERtype == qh_ASvoronoi && !facetA->center  /* matches upperdelaunay in qh_setfacetplane() */
02907         && fabs_(facetA->normal[qh hull_dim -1]) >= qh ANGLEround * qh_ZEROdelaunay) {
02908     facetA->center= qh_facetcenter(facetA->vertices);
02909   }
02910   qh_willdelete(facetA, NULL);
02911   qh newfacet_list= qh facet_tail;
02912   facetA->visitid= qh visit_id;
02913   apex= SETfirstt_(facetA->vertices, vertexT);
02914   qh_makenew_nonsimplicial(facetA, apex, &numnew);
02915   SETfirst_(facetA->neighbors)= NULL;
02916   FORALLnew_facets {
02917     newfacet->tricoplanar= True;
02918     newfacet->f.trivisible= facetA;
02919     newfacet->degenerate= False;
02920     newfacet->upperdelaunay= facetA->upperdelaunay;
02921     newfacet->good= facetA->good;
02922     if (qh TRInormals) {
02923       newfacet->keepcentrum= True;
02924       newfacet->normal= qh_copypoints(facetA->normal, 1, qh hull_dim);
02925       if (qh CENTERtype == qh_AScentrum)
02926         newfacet->center= qh_getcentrum(newfacet);
02927       else
02928         newfacet->center= qh_copypoints(facetA->center, 1, qh hull_dim);
02929     }else {
02930       newfacet->keepcentrum= False;
02931       newfacet->normal= facetA->normal;
02932       newfacet->center= facetA->center;
02933     }
02934     newfacet->offset= facetA->offset;
02935 #if qh_MAXoutside
02936     newfacet->maxoutside= facetA->maxoutside;
02937 #endif
02938   }
02939   qh_matchnewfacets(/*qh newfacet_list*/);
02940   zinc_(Ztricoplanar);
02941   zadd_(Ztricoplanartot, numnew);
02942   zmax_(Ztricoplanarmax, numnew);
02943   qh visible_list= NULL;
02944   if (!(*first_vertex))
02945     (*first_vertex)= qh newvertex_list;
02946   qh newvertex_list= NULL;
02947   qh_updatevertices(/*qh newfacet_list, empty visible_list and newvertex_list*/);
02948   qh_resetlists(False, !qh_RESETvisible /*qh newfacet_list, empty visible_list and newvertex_list*/);
02949 } /* triangulate_facet */
02950 
02951 /*-<a                             href="qh-poly.htm#TOC"
02952   >-------------------------------</a><a name="triangulate_link">-</a>
02953 
02954   qh_triangulate_link(oldfacetA, facetA, oldfacetB, facetB)
02955     relink facetA to facetB via oldfacets
02956   returns:
02957     adds mirror facets to qh degen_mergeset (4-d and up only)
02958   design:
02959     if they are already neighbors, the opposing neighbors become MRGmirror facets
02960 */
02961 void qh_triangulate_link(facetT *oldfacetA, facetT *facetA, facetT *oldfacetB, facetT *facetB) {
02962   int errmirror= False;
02963 
02964   trace3((qh ferr, 3021, "qh_triangulate_link: relink old facets f%d and f%d between neighbors f%d and f%d\n",
02965          oldfacetA->id, oldfacetB->id, facetA->id, facetB->id));
02966   if (qh_setin(facetA->neighbors, facetB)) {
02967     if (!qh_setin(facetB->neighbors, facetA))
02968       errmirror= True;
02969     else
02970       qh_appendmergeset(facetA, facetB, MRGmirror, NULL);
02971   }else if (qh_setin(facetB->neighbors, facetA))
02972     errmirror= True;
02973   if (errmirror) {
02974     qh_fprintf(qh ferr, 6163, "qhull error (qh_triangulate_link): mirror facets f%d and f%d do not match for old facets f%d and f%d\n",
02975        facetA->id, facetB->id, oldfacetA->id, oldfacetB->id);
02976     qh_errexit2 (qh_ERRqhull, facetA, facetB);
02977   }
02978   qh_setreplace(facetB->neighbors, oldfacetB, facetA);
02979   qh_setreplace(facetA->neighbors, oldfacetA, facetB);
02980 } /* triangulate_link */
02981 
02982 /*-<a                             href="qh-poly.htm#TOC"
02983   >-------------------------------</a><a name="triangulate_mirror">-</a>
02984 
02985   qh_triangulate_mirror(facetA, facetB)
02986     delete mirrored facets from qh_triangulate_null() and qh_triangulate_mirror
02987       a mirrored facet shares the same vertices of a logical ridge
02988   design:
02989     since a null facet duplicates the first two vertices, the opposing neighbors absorb the null facet
02990     if they are already neighbors, the opposing neighbors become MRGmirror facets
02991 */
02992 void qh_triangulate_mirror(facetT *facetA, facetT *facetB) {
02993   facetT *neighbor, *neighborB;
02994   int neighbor_i, neighbor_n;
02995 
02996   trace3((qh ferr, 3022, "qh_triangulate_mirror: delete mirrored facets f%d and f%d\n",
02997          facetA->id, facetB->id));
02998   FOREACHneighbor_i_(facetA) {
02999     neighborB= SETelemt_(facetB->neighbors, neighbor_i, facetT);
03000     if (neighbor == neighborB)
03001       continue; /* occurs twice */
03002     qh_triangulate_link(facetA, neighbor, facetB, neighborB);
03003   }
03004   qh_willdelete(facetA, NULL);
03005   qh_willdelete(facetB, NULL);
03006 } /* triangulate_mirror */
03007 
03008 /*-<a                             href="qh-poly.htm#TOC"
03009   >-------------------------------</a><a name="triangulate_null">-</a>
03010 
03011   qh_triangulate_null(facetA)
03012     remove null facetA from qh_triangulate_facet()
03013       a null facet has vertex #1 (apex) == vertex #2
03014   returns:
03015     adds facetA to ->visible for deletion after qh_updatevertices
03016     qh degen_mergeset contains mirror facets (4-d and up only)
03017   design:
03018     since a null facet duplicates the first two vertices, the opposing neighbors absorb the null facet
03019     if they are already neighbors, the opposing neighbors become MRGmirror facets
03020 */
03021 void qh_triangulate_null(facetT *facetA) {
03022   facetT *neighbor, *otherfacet;
03023 
03024   trace3((qh ferr, 3023, "qh_triangulate_null: delete null facet f%d\n", facetA->id));
03025   neighbor= SETfirstt_(facetA->neighbors, facetT);
03026   otherfacet= SETsecondt_(facetA->neighbors, facetT);
03027   qh_triangulate_link(facetA, neighbor, facetA, otherfacet);
03028   qh_willdelete(facetA, NULL);
03029 } /* triangulate_null */
03030 
03031 #else /* qh_NOmerge */
03032 void qh_triangulate(void) {
03033 }
03034 #endif /* qh_NOmerge */
03035 
03036    /*-<a                             href="qh-poly.htm#TOC"
03037   >-------------------------------</a><a name="vertexintersect">-</a>
03038 
03039   qh_vertexintersect( vertexsetA, vertexsetB )
03040     intersects two vertex sets (inverse id ordered)
03041     vertexsetA is a temporary set at the top of qhmem.tempstack
03042 
03043   returns:
03044     replaces vertexsetA with the intersection
03045 
03046   notes:
03047     could overwrite vertexsetA if currently too slow
03048 */
03049 void qh_vertexintersect(setT **vertexsetA,setT *vertexsetB) {
03050   setT *intersection;
03051 
03052   intersection= qh_vertexintersect_new(*vertexsetA, vertexsetB);
03053   qh_settempfree(vertexsetA);
03054   *vertexsetA= intersection;
03055   qh_settemppush(intersection);
03056 } /* vertexintersect */
03057 
03058 /*-<a                             href="qh-poly.htm#TOC"
03059   >-------------------------------</a><a name="vertexintersect_new">-</a>
03060 
03061   qh_vertexintersect_new(  )
03062     intersects two vertex sets (inverse id ordered)
03063 
03064   returns:
03065     a new set
03066 */
03067 setT *qh_vertexintersect_new(setT *vertexsetA,setT *vertexsetB) {
03068   setT *intersection= qh_setnew(qh hull_dim - 1);
03069   vertexT **vertexA= SETaddr_(vertexsetA, vertexT);
03070   vertexT **vertexB= SETaddr_(vertexsetB, vertexT);
03071 
03072   while (*vertexA && *vertexB) {
03073     if (*vertexA  == *vertexB) {
03074       qh_setappend(&intersection, *vertexA);
03075       vertexA++; vertexB++;
03076     }else {
03077       if ((*vertexA)->id > (*vertexB)->id)
03078         vertexA++;
03079       else
03080         vertexB++;
03081     }
03082   }
03083   return intersection;
03084 } /* vertexintersect_new */
03085 
03086 /*-<a                             href="qh-poly.htm#TOC"
03087   >-------------------------------</a><a name="vertexneighbors">-</a>
03088 
03089   qh_vertexneighbors()
03090     for each vertex in qh.facet_list,
03091       determine its neighboring facets
03092 
03093   returns:
03094     sets qh.VERTEXneighbors
03095       nop if qh.VERTEXneighbors already set
03096       qh_addpoint() will maintain them
03097 
03098   notes:
03099     assumes all vertex->neighbors are NULL
03100 
03101   design:
03102     for each facet
03103       for each vertex
03104         append facet to vertex->neighbors
03105 */
03106 void qh_vertexneighbors(void /*qh facet_list*/) {
03107   facetT *facet;
03108   vertexT *vertex, **vertexp;
03109 
03110   if (qh VERTEXneighbors)
03111     return;
03112   trace1((qh ferr, 1035, "qh_vertexneighbors: determing neighboring facets for each vertex\n"));
03113   qh vertex_visit++;
03114   FORALLfacets {
03115     if (facet->visible)
03116       continue;
03117     FOREACHvertex_(facet->vertices) {
03118       if (vertex->visitid != qh vertex_visit) {
03119         vertex->visitid= qh vertex_visit;
03120         vertex->neighbors= qh_setnew(qh hull_dim);
03121       }
03122       qh_setappend(&vertex->neighbors, facet);
03123     }
03124   }
03125   qh VERTEXneighbors= True;
03126 } /* vertexneighbors */
03127 
03128 /*-<a                             href="qh-poly.htm#TOC"
03129   >-------------------------------</a><a name="vertexsubset">-</a>
03130 
03131   qh_vertexsubset( vertexsetA, vertexsetB )
03132     returns True if vertexsetA is a subset of vertexsetB
03133     assumes vertexsets are sorted
03134 
03135   note:
03136     empty set is a subset of any other set
03137 */
03138 boolT qh_vertexsubset(setT *vertexsetA, setT *vertexsetB) {
03139   vertexT **vertexA= (vertexT **) SETaddr_(vertexsetA, vertexT);
03140   vertexT **vertexB= (vertexT **) SETaddr_(vertexsetB, vertexT);
03141 
03142   while (True) {
03143     if (!*vertexA)
03144       return True;
03145     if (!*vertexB)
03146       return False;
03147     if ((*vertexA)->id > (*vertexB)->id)
03148       return False;
03149     if (*vertexA  == *vertexB)
03150       vertexA++;
03151     vertexB++;
03152   }
03153   return False; /* avoid warnings */
03154 } /* vertexsubset */