/* poly2.c -- implements polygons and simplices   see README, poly.h and qhull.h      copyright (c) 1993-1995, The Geometry Center   frequently used code is in poly.c*/#include "qhull_a.h"/*======== functions in alphabetical order ==========*//*------------------------------------------------addhash- add hash element to linear hash table if not already there*/void qh_addhash (void* newelem, setT *hashtable, int hashsize, unsigned hash) {  unsigned scan;  void *elem;  for (scan= hash; (elem= SETelem_(hashtable, scan));        scan= (++scan >= hashsize ? 0 : scan)) {    if (elem == newelem)      break;  }  /* loop terminates because qh_HASHfactor >= 1.1 by qh_initbuffers */  if (!elem)    SETelem_(hashtable, scan)= newelem;} /* addhash *//*------------------------------------------------check_bestdist- check that points are within max_outside of the nearest facet  if ONLYgood, ignores !good facets  see: check_maxout*/void qh_check_bestdist (void) {  boolT waserror= False, isoutside, unassigned;  facetT *facet, *bestfacet, *errfacet1= NULL, *errfacet2= NULL;  facetT *facetlist;   realT dist, maxoutside;  pointT *point;  int numpart, facet_i, facet_n, notgood= 0, notverified= 0;  setT *facets;  maxoutside= fmax_(qh max_outside, qh DISTround);  maxoutside += 2 * qh DISTround;  /* 1 DISTround to actual point and another DISTround to computed point */  if (qh RANDOMdist) /* repeated computations can differ by 2*distround */    maxoutside += qh DISTround;  trace1((qh ferr, "qh_check_bestdist: check that all points are within %2.2g of best facet\n", maxoutside));  facets= qh_pointfacet (/*qh facet_list*/);  if (!qh_QUICKhelp && qh PRINTprecision)    fprintf (qh ferr, "\n\qhull output completed.  Verifying that %d points are\n\below %2.2g of the nearest %sfacet.\n",	     qh_setsize(facets), maxoutside, (qh ONLYgood ?  "good " : ""));  FOREACHfacet_i_(facets) {  /* for each point with facet assignment */    if (facet)      unassigned= False;    else {      unassigned= True;      facet= qh facet_list;    }    point= qh_point(facet_i);    if (point == qh GOODpointp)      continue;    bestfacet= qh_findbest (point, facet, qh_ALL, False,			    &dist, &isoutside, &numpart);    /* occurs after statistics reported */    if (dist > maxoutside) {      if (qh ONLYgood && !bestfacet->good 	  && !((bestfacet= qh_findgooddist (point, bestfacet, &dist, &facetlist))	       && dist > maxoutside))	notgood++;      else {	waserror= True;	fprintf(qh ferr, "qhull precision error: point p%d is outside facet f%d, distance= %6.8g maxoutside= %6.8g\n", 		facet_i, bestfacet->id, dist, maxoutside);	errfacet2= errfacet1;	errfacet1= bestfacet;		          }    }else if (unassigned && dist < -qh MAXcoplanar)      notverified++;  }  qh_settempfree (&facets);  if (notverified && !qh DELAUNAY && !qh_QUICKhelp && qh PRINTprecision)     fprintf(qh ferr, "%d points were well inside the hull.  If the hull contains\n\a lens-shaped component, these points were not verified.  Use\n\options 'Qi Tv' to verify all points.\n", notverified);   if (waserror)    qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);} /* check_bestdist *//*------------------------------------------------check_maxout- updates max_outside by checking all points against bestfacet  updates facet->maxoutside via findbest  if printing min_vertex, it is updated to the current vertices  if ONLYgood, ignores !good facets  see check_bestdistnotes:  may not need to check near-inside points if KEEPcoplanar      (since coplanar is now MAXcoplanar instead of -DISTround)*/#ifndef qh_NOmergevoid qh_check_maxout (void) {  facetT *facet, *bestfacet, *neighbor, **neighborp, *facetlist;  realT dist, maxoutside, minvertex;  pointT *point, **pointp;  int numpart, facet_i, facet_n, notgood= 0;  setT *facets, *vertices;  vertexT *vertex;  maxoutside= minvertex= 0;  trace1((qh ferr, "qh_check_maxout: determine actual maxoutside and minoutside\n"));  if (qh VERTEXneighbors   && (qh PRINTsummary || qh PRINTout[0] == qh_PRINTnone || qh PRINTstatistics      || qh PRINTout[0] == qh_PRINTsummary || qh TRACElevel)) {    vertices= qh_pointvertex (/*qh facet_list*/);    FORALLvertices {      FOREACHneighbor_(vertex) {	zinc_(Zdistvertex);  /* distance also computed by main loop below */	qh_distplane (vertex->point, neighbor, &dist);	minimize_(minvertex, dist);	if (-dist > qh TRACEdist || dist > qh TRACEdist 	    || neighbor == qh tracefacet || vertex == qh tracevertex)	  fprintf (qh ferr, "qh_check_maxout: p%d (v%d) is %.2g from f%d\n",		   qh_pointid (vertex->point), vertex->id, dist, neighbor->id);      }    }    if (qh MERGING) {      wmin_(Wminvertex, qh min_vertex);    }    qh min_vertex= minvertex;    qh_settempfree (&vertices);  }  facets= qh_pointfacet (/*qh facet_list*/);  FOREACHfacet_i_(facets) {     /* for each point with facet assignment */    if (facet) {       zinc_(Ztotcheck);      point= qh_point(facet_i);      if (point == qh GOODpointp)	continue;      bestfacet= qh_findbest (point, facet, qh_ALL,			    False, &dist, NULL, &numpart);      zadd_(Zcheckpart, numpart);      if (bestfacet && dist > maxoutside) {        if (qh ONLYgood && !bestfacet->good         && !((bestfacet= qh_findgooddist (point, bestfacet, &dist, &facetlist))             && dist > maxoutside))          notgood++;        else	  maxoutside= dist;      }      if (dist > qh TRACEdist || (bestfacet && bestfacet == qh tracefacet))	fprintf (qh ferr, "qh_check_maxout: p%d is %.2g above f%d\n",		   qh_pointid (point), dist, bestfacet->id);    }  }  qh_settempfree (&facets);  wval_(Wmaxout)= maxoutside - qh max_outside;  wmax_(Wmaxoutside, qh max_outside);  qh max_outside= maxoutside;  if (qh KEEPnearinside && !qh KEEPcoplanar && !qh KEEPinside) {    FORALLfacets {      if (facet->coplanarset)         qh_setfree( &facet->coplanarset);    }  }else if (qh KEEPnearinside) {    numpart= 0;    FORALLfacets {   /* could be combined with qh_findbest */      if (facet->coplanarset) {        FOREACHpoint_(facet->coplanarset) {  	  numpart++;	  qh_distplane (point, facet, &dist); 	  if (dist < -qh MAXcoplanar) {	    if (!qh KEEPinside)              SETref_(point)= NULL;          }else if (!qh KEEPcoplanar)            SETref_(point)= NULL;        }	qh_setcompact (facet->coplanarset);      }    }    zadd_(Zcheckpart, numpart);  }  trace1((qh ferr, "qh_check_maxout: maxoutside %2.2g, minvertex %2.2g, outside of not good %d\n",       maxoutside, minvertex, notgood));} /* check_maxout */#else /* qh_NOmerge */void qh_check_maxout (void) {}#endif/*-----------------------------------------check_output- performs the checks at the end of qhull algorithm  does not check points (may take a long time)*/void qh_check_output (void) {  int i;  if (qh STOPcone)    return;  if (qh VERIFYoutput | qh IStracing | qh CHECKfrequently) {    qh_checkpolygon (qh facet_list);    qh_checkflipped_all (qh facet_list);    qh_checkconvex (qh facet_list, qh_ALGORITHMfault);  }else if (!qh MERGING && qh_newstats (qhstat precision, &i)) {    qh_checkflipped_all (qh facet_list);    qh_checkconvex (qh facet_list, qh_ALGORITHMfault);  }} /* check_output *//*--------------------------------------------------------------check_point- check that point is not outside facet  if maxerror, doesn't report an error*/void qh_check_point (pointT *point, facetT *facet, realT *maxoutside, facetT **errfacet1, facetT **errfacet2) {  realT dist;  /* occurs after statistics reported */  qh_distplane(point, facet, &dist);  if (dist > *maxoutside) {    *errfacet2= *errfacet1;    *errfacet1= facet;    fprintf(qh ferr, "qhull precision error: point p%d is outside facet f%d, distance= %6.8g maxoutside= %6.8g\n", 	      qh_pointid(point), facet->id, dist, *maxoutside);  }} /* qh_check_point *//*--------------------------------------------------check_points- checks that all points are inside all facets     uses findbest if lots of points     ignores flipped facetsnotes:  maxoutside includes 2 DISTrounds, one for the computed  distances in qh_check_points  qh_printafacet and qh_printsummary needs only one DISTround*/void qh_check_points (void) {  facetT *facet, *errfacet1= NULL, *errfacet2= NULL;  realT total, maxoutside;  pointT *point, **pointp, *pointtemp;  boolT testouter;  maxoutside= fmax_(qh max_outside, qh DISTround);  /* agrees with qh_printafacet, qh_printsummary */  maxoutside += 2* qh DISTround;  /* 1 DISTround to actual point and another DISTround to computed point */  if (qh RANDOMdist) /* repeated computations can differ by 2*distround */    maxoutside += qh DISTround;  trace1((qh ferr, "qh_check_points: check all points below %2.2g of all facet planes\n",	  maxoutside));  if (qh num_good)   /* miss counts other_points and !good facets */     total= (float) qh num_good * qh num_points;  else     total= (float) qh num_facets * qh num_points;  if (total >= qh_VERIFYdirect   && (!qh MERGING || qh SKIPcheckmax || qh ZEROall_ok)) {    if (!qh_QUICKhelp && qh SKIPcheckmax && qh MERGING)      fprintf (qh ferr, "\n\qhull input warning: merging without checking outer planes ('Q5').\n\Verify may report that a point is outside of a facet.\n");    qh_check_bestdist();  }else {    if (qh MERGING && qh_MAXoutside && !qh SKIPcheckmax)      testouter= True;  /* agrees with qh_printafacet() */    else      testouter= False;    if (!qh_QUICKhelp) {      if (qh MERGEexact || qh SKIPcheckmax || qh NOnearinside)	fprintf (qh ferr, "\n\qhull input warning: exact merge ('Qx'), no outer plane check ('Q5'), or\n\no processing of near-inside points ('Q8').  Verify may report that a point\n\is outside of a facet.\n");    }    if (qh PRINTprecision) {      if (testouter)	fprintf (qh ferr, "\n\Output completed.  Verifying that all points are below outer planes of\n\all %sfacets.  Will make %2.0f distance computations.\n", 	      (qh ONLYgood ?  "good " : ""), total);      else	fprintf (qh ferr, "\n\Output completed.  Verifying that all points are below %2.2g of\n\all %sfacets.  Will make %2.0f distance computations.\n", 	      maxoutside, (qh ONLYgood ?  "good " : ""), total);    }    FORALLfacets {      if (!facet->good && qh ONLYgood)        continue;      if (facet->flipped)        continue;      if (testouter) {#if qh_MAXoutside	maxoutside= facet->maxoutside + 2* qh DISTround;	/* 1 DISTround to actual point and another to computed point */#endif      }      FORALLpoints {	if (point != qh GOODpointp)	  qh_check_point (point, facet, &maxoutside, &errfacet1, &errfacet2);      }      FOREACHpoint_(qh other_points) {	if (point != qh GOODpointp)	  qh_check_point (point, facet, &maxoutside, &errfacet1, &errfacet2);      }    }    if (errfacet1)      qh_errexit2(qh_ERRprec, errfacet1, errfacet2);  }} /* check_points *//*--------------------------------------------------checkconvex- check that each ridge in facetlist is convexreturns:    counts Zconcaveridges and Zcoplanarridges    errors if concaveridge or if merging an coplanar ridgenote:    if not merging, tests vertices for neighboring simplicial facets    else if ZEROcentrum, tests vertices for neighboring simplicial facets    else tests centrums of neighboring facets*/void qh_checkconvex(facetT *facetlist, int fault) {  facetT *facet, *neighbor, **neighborp, *errfacet1=NULL, *errfacet2=NULL;  vertexT *vertex;  realT dist;  pointT *centrum;  boolT waserror= False, tempcentrum= False, allsimplicial;  int neighbor_i;  trace1((qh ferr, "qh_checkconvex: check all ridges are convex\n"));  zzval_(Zconcaveridges)= 0;  zzval_(Zcoplanarridges)= 0;  FORALLfacet_(facetlist) {    if (facet->flipped) {      fprintf (qh ferr, "qhull precision error: f%d is flipped (interior point is outside)\n",	       facet->id);      errfacet1= facet;      waserror= True;      continue;    }    if (qh MERGING && (!qh ZEROcentrum || !facet->simplicial))      allsimplicial= False;    else {      allsimplicial= True;      neighbor_i= 0;      FOREACHneighbor_(facet) {        vertex= SETelem_(facet->vertices, neighbor_i++);	if (!neighbor->simplicial) {	  allsimplicial= False;	  continue;	}        qh_distplane (vertex->point, neighbor, &dist);        if (dist > -qh DISTround) {	  if (fault == qh_DATAfault) {	    fprintf (qh ferr, "qhull precision error: initial simplex is not convex. Distance=%.2g\n", dist);	    qh_errexit(qh_ERRsingular, NULL, NULL);	  }          if (dist > qh DISTround) {            zzinc_(Zconcaveridges);            fprintf (qh ferr, "qhull precision error: f%d is concave to f%d, since p%d (v%d) is %6.4g above\n",              facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist);            errfacet1= facet;            errfacet2= neighbor;            waserror= True;          }else if (qh ZEROcentrum) {            if (dist > 0) {     /* qh_checkzero checks that dist < - qh DISTround */              zzinc_(Zcoplanarridges);               fprintf (qh ferr, "qhull precision error: f%d is clearly not convex to f%d, since p%d (v%d) is %6.4g above\n",                facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist);