/* qhull - Quickhull algorithm for convex hulls   qhull() and top-level routines   see README, qhull.h, unix.c and mac.c   see qhull_a.h for internal functions      copyright (c) 1993-1995 The Geometry Center        */#include "qhull_a.h" /*--------------------------------------------------qhull- hull_dim convex hull of num_points starting at first_pointreturns:  returns facet_list, numfacets, etc. */void qh_qhull (void) {  setT *maxpoints, *vertices;  facetT *facet;  int numpart, i, numoutside;  realT dist;  boolT isoutside;  qh hulltime= (unsigned)clock();  if (qh DELAUNAY && qh upper_threshold[qh hull_dim-1] > REALmax/2                  && qh lower_threshold[qh hull_dim-1] < -REALmax/2) {    for (i= qh_PRINTEND; i--; ) {      if (qh PRINTout[i] == qh_PRINTgeom && qh DROPdim < 0  	  && !qh GOODthreshold && !qh SPLITthresholds)	break;  /* in this case, don't set upper_threshold */    }    if (i < 0) {      qh upper_threshold[qh hull_dim-1]= 0.0;      if (!qh GOODthreshold)	qh SPLITthresholds= True;    }  }  maxpoints= qh_maxmin(qh first_point, qh num_points, qh hull_dim);  /* qh_maxmin sets DISTround and other precision constants */  if (qh PRINToptions1st || qh TRACElevel || qh IStracing) {    if (qh TRACElevel || qh IStracing)      fprintf (qh ferr, "\nTrace level %d for %s | %s\n",          qh IStracing ? qh IStracing : qh TRACElevel, qh rbox_command, qh qhull_command);    fprintf (qh ferr, "Options selected for qhull %s:\n%s\n", qh_version, qh qhull_options);  }  vertices= qh_initialvertices(qh hull_dim, maxpoints, qh first_point, qh num_points);   qh_initialhull (vertices);  /* initial qh facet_list */  qh_partitionall (vertices, qh first_point, qh num_points);  if (qh PREmerge) {    qh cos_max= qh premerge_cos;    qh centrum_radius= qh premerge_centrum;  }  if (qh ONLYgood) {    if (!(qh GOODthreshold || qh GOODpoint	  || (qh GOODvertex > 0 && !qh MERGING))) {      fprintf (qh ferr, "qhull input error: 'Qg' (ONLYgood) needs a good threshold ('Pd0D0'), a\n\good point (QGn or QG-n), or a good vertex without merging (QVn).\n");      qh_errexit (qh_ERRinput, NULL, NULL);    }    if (qh GOODvertex > 0  && !qh MERGING  /* matches qh_partitionall */	&& !qh_isvertex (qh GOODvertexp, vertices)) {      facet= qh_findbestnew (qh GOODvertexp, qh facet_list, 			  &dist, &isoutside, &numpart);      zadd_(Zdistgood, numpart);      if (!isoutside) {        fprintf (qh ferr, "qhull input error: point for QV%d is inside initial simplex\n",	       qh_pointid(qh GOODvertexp));        qh_errexit (qh_ERRinput, NULL, NULL);      }      if (!qh_addpoint (qh GOODvertexp, facet, False)) {	qh_settempfree(&vertices);	qh_settempfree(&maxpoints);	return;      }    }    qh_findgood (qh facet_list, 0);  }  qh_settempfree(&vertices);  qh_settempfree(&maxpoints);  qh_resetlists (False /*qh visible_list newvertex_list newfacet_list */);  qh_buildhull();  if (!qh STOPpoint && !qh STOPcone) {    if (qh ZEROall_ok && !qh TESTvneighbors && qh MERGEexact)      qh_checkzero( qh_ALL);    if (qh ZEROall_ok && !qh TESTvneighbors) {      trace2((qh ferr, "qh_qhull: all facets are clearly convex.  Post-merging not needed.\n"));    }else {      if (qh MERGEexact || (qh hull_dim > qh_DIMreduceBuild && qh PREmerge))        qh_postmerge ("First post-merge", qh premerge_centrum, qh premerge_cos,              (qh POSTmerge ? False : qh TESTvneighbors));      else if (!qh POSTmerge && qh TESTvneighbors)         qh_postmerge ("For testing vertex neighbors", qh premerge_centrum,             qh premerge_cos, True);       if (qh POSTmerge)        qh_postmerge ("For post-merging", qh postmerge_centrum,              qh postmerge_cos, qh TESTvneighbors);      if (qh visible_list == qh facet_list) { /* i.e., merging done */        qh findbestnew= True;        qh_partitionvisible (/*visible_list, newfacet_list*/ !qh_ALL, &numoutside);        qh findbestnew= False;        qh_deletevisible (/*qh visible_list*/);        qh_resetlists (False /*qh visible_list newvertex_list newfacet_list */);      }      if (qh DOcheckmax){        if (qh REPORTfreq) {	  qh_buildtracing (NULL, NULL);   	  fprintf (qh ferr, "\nTesting all coplanar points.\n");        }        qh_check_maxout();      }    }  }  if (qh_setsize (qhmem.tempstack) != 0) {    fprintf (qh ferr, "qhull internal error (qh_qhull): temporary sets not empty (%d)\n",	     qh_setsize (qhmem.tempstack));    qh_errexit (qh_ERRqhull, NULL, NULL);  }  qh hulltime= (unsigned)clock() - qh hulltime;  qh QHULLfinished= True;  trace1((qh ferr, "qh_qhull: algorithm completed\n"));} /* qhull *//*--------------------------------------------------addpoint-  add point to hull above a facet  if !facet, locates a facet for the point    !facet works for Delaunay triangulations    !facet does not work for lens-shaped hulls    if point is not outside of the hull, uses qh_partitioncoplanar()  if checkdist, checks that point is outside of the facet.    if not outside, uses qh_partitioncoplanar()    if facet->upperdelaunay, assumes that the point is above facet  if !checkdist and facet, assumes point is above facet (major damage if below)returns:  if unknown point, adds it to qh other_points  returns False if user requested break     visible_list, newfacet_list, delvertex_list, NEWfacets may be defined*/boolT qh_addpoint (pointT *furthest, facetT *facet, boolT checkdist) {  int goodvisible, goodhorizon;  vertexT *vertex;  facetT *newfacet;  realT dist, newbalance, pbalance;  boolT isoutside= False;  int numpart, numpoints, numnew, firstnew;  if (qh_pointid (furthest) == -1)    qh_setappend (&qh other_points, furthest);  if (!facet   || (checkdist && !facet->upperdelaunay)) { /* else missed by qh_findbest */    if (!facet)      facet= qh_nonupper( qh facet_list);    facet= qh_findbest (furthest, facet, !qh_ALL, False, &dist, &isoutside, &numpart);    zzadd_(Zpartition, numpart);    if (!isoutside) {      zinc_(Znotmax);  /* last point of outsideset is no longer furthest. */      qh_partitioncoplanar (furthest, facet, &dist);      return True;    }  }  qh_buildtracing (furthest, facet);  if (qh STOPpoint < 0 && qh furthest_id == -qh STOPpoint-1)    return False;  qh_findhorizon (furthest, facet, &goodvisible, &goodhorizon);   if (qh ONLYgood && !(goodvisible+goodhorizon)) {    zinc_(Znotgood);      /* last point of outsideset is no longer furthest.  This is ok       since all points of the outside are likely to be bad */    qh_resetlists (False /*qh visible_list newvertex_list newfacet_list */);    return True;  }  zzinc_(Zprocessed);  firstnew= qh facet_id;  vertex= qh_makenewfacets (furthest /*visible_list, attaches if !ONLYgood */);  qh_makenewplanes (/* newfacet_list */);  numnew= qh facet_id - firstnew;  newbalance= numnew - (realT) (qh num_facets-qh num_visible)                         * qh hull_dim/qh num_vertices;  wadd_(Wnewbalance, newbalance);  wadd_(Wnewbalance2, newbalance * newbalance);  if (qh ONLYgood && !qh_findgood (qh newfacet_list, goodhorizon)) {    FORALLnew_facets       qh_delfacet (newfacet);    qh_delvertex (vertex);    qh_resetlists (True /*qh visible_list newvertex_list newfacet_list */);    zinc_(Znotgoodnew);    return True;  }  if (qh ONLYgood)    qh_attachnewfacets(/*visible_list*/);  qh_matchnewfacets();  qh_updatevertices();  if (qh STOPcone && qh furthest_id == qh STOPcone-1)    return False;  /* visible_list etc. still defined */  if (qh PREmerge || qh MERGEexact) {    qh_premerge (vertex, qh premerge_centrum, qh premerge_cos);    if (zzval_(Ztotmerge) > qh_USEfindbestnew)      qh findbestnew= True;    else {      FORALLnew_facets {	if (!newfacet->simplicial) {	  qh findbestnew= True;  /* qh_findbest can not be used */	  break;	}      }    }  }else if (qh BESToutside)    qh findbestnew= True;  qh_partitionvisible (/*visible_list, newfacet_list*/ !qh_ALL, &numpoints);  qh findbestnew= False;  qh findbest_notsharp= False;  zinc_(Zpbalance);  pbalance= numpoints - (realT) qh hull_dim /* assumes all points extreme */                * (qh num_points - qh num_vertices)/qh num_vertices;  wadd_(Wpbalance, pbalance);  wadd_(Wpbalance2, pbalance * pbalance);  qh_deletevisible (/*qh visible_list*/);  zmax_(Zmaxvertex, qh num_vertices);  qh NEWfacets= False;  if (qh IStracing >= 4)    qh_printfacetlist (qh newfacet_list, NULL, True);  if (qh CHECKfrequently) {    if (qh num_facets < 50)      qh_checkpolygon (qh facet_list);    else      qh_checkpolygon (qh newfacet_list);  }  if (qh STOPpoint > 0 && qh furthest_id == qh STOPpoint-1)    return False;   qh_resetlists (True /*qh visible_list newvertex_list newfacet_list */);  trace2((qh ferr, "qh_addpoint: added p%d new facets %d new balance %2.2g point balance %2.2g\n",    qh_pointid (furthest), numnew, newbalance, pbalance));  return True;} /* addpoint *//*--------------------------------------------------buildhull- constructs a hull by adding outside points one at a time  may be called multiple times  checks facet and vertex lists for 'visible', 'newfacet', and 'newlist'notes:  to recover from STOPcone, call qh_deletevisible and qh_resetlists*/void qh_buildhull(void) {  facetT *facet;  pointT *furthest;  vertexT *vertex;  int id;    trace1((qh ferr, "qh_buildhull: start build hull\n"));  FORALLfacets {    if (facet->visible || facet->newfacet) {      fprintf (qh ferr, "qhull internal error (qh_buildhull): visible or new facet f%d in facet list\n",                   facet->id);          qh_errexit (qh_ERRqhull, facet, NULL);    }  }  FORALLvertices {    if (vertex->newlist) {      fprintf (qh ferr, "qhull internal error (qh_buildhull): new vertex f%d in vertex list\n",                   vertex->id);      qh_errprint ("ERRONEOUS", NULL, NULL, NULL, vertex);      qh_errexit (qh_ERRqhull, NULL, NULL);    }    id= qh_pointid (vertex->point);    if ((qh STOPpoint>0 && id == qh STOPpoint-1) ||	(qh STOPpoint<0 && id == -qh STOPpoint-1) ||	(qh STOPcone>0 && id == qh STOPcone-1)) {      trace1((qh ferr,"qh_buildhull: stop point or cone P%d in initial hull\n", id));      return;    }  }  qh facet_next= qh facet_list;      /* advance facet when processed */  while ((furthest= qh_nextfurthest (&facet))) {    qh num_outside--;  /* if ONLYmax, furthest may not be outside */    if (!qh_addpoint (furthest, facet, qh ONLYmax))      break;  }  if (qh num_outside && !furthest) {    fprintf (qh ferr, "qhull internal error (qh_buildhull): %d outside points were never processed.\n", qh num_outside);    qh_errexit (qh_ERRqhull, NULL, NULL);  }  trace1((qh ferr, "qh_buildhull: completed the hull construction\n"));} /* buildhull */  /*--------------------------------------------buildtracing- for tracing execution of buildhull  tracks progress with qh lastreport  updates qh furthest_id (-3 if furthest is NULL)  also resets visit_id, vertext_visit on wrap around  if !furthest, prints basic message  see also qh_tracemerging()*/void qh_buildtracing (pointT *furthest, facetT *facet) {  realT dist= 0;  float cpu;  int total, furthestid;  time_t timedata;  struct tm *tp;  vertexT *vertex;  if (!furthest) {    time (&timedata);    tp= localtime (&timedata);    cpu= (unsigned)clock() - qh hulltime;    cpu /= qh_SECticks;    total= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);    fprintf (qh ferr, "\n\