poly2.c

Quick hull implementation

        facet->good= False;
      }
    }
  }
  startgood= numgood;
  if (qh SPLITthresholds) {
    FORALLfacet_(facetlist) {
      if (facet->good) {
        if (!qh_inthresholds (facet->normal, &angle)) {
          facet->good= False;
          numgood--;
          if (angle < bestangle) {
            bestangle= angle;
            bestfacet= facet;
          }
        }
      }
    }
    if (!numgood && bestfacet) {
      bestfacet->good= True;
      numgood++;
      trace0((qh ferr, "qh_findgood_all: f%d is closest (%2.2g) to thresholds\n", 
           bestfacet->id, bestangle));
      return;
    }
  }
  qh num_good= numgood;
  trace0((qh ferr, "qh_findgood_all: %d good facets remain out of %d facets\n",
        numgood, startgood));
} /* findgood_all */

/*-<a                             href="qh-poly.htm#TOC"
  >-------------------------------</a><a name="furthestnext">-</a>
  
  qh_furthestnext()
    set qh.facet_next to facet with furthest of all furthest points
    searches all facets on qh.facet_list

  notes:
    this may help avoid precision problems
*/
void qh_furthestnext (void /* qh facet_list */) {
  facetT *facet, *bestfacet= NULL;
  realT dist, bestdist= -REALmax;

  FORALLfacets {
    if (facet->outsideset) {
#if qh_COMPUTEfurthest
      pointT *furthest;
      furthest= (pointT*)qh_setlast (facet->outsideset);
      zinc_(Zcomputefurthest);
      qh_distplane (furthest, facet, &dist);
#else
      dist= facet->furthestdist;
#endif
      if (dist > bestdist) {
	bestfacet= facet;
	bestdist= dist;
      }
    }
  }
  if (bestfacet) {
    qh_removefacet (bestfacet);
    qh_prependfacet (bestfacet, &qh facet_next);
    trace1((qh ferr, "qh_furthestnext: made f%d next facet (dist %.2g)\n",
	    bestfacet->id, bestdist));
  }
} /* furthestnext */

/*-<a                             href="qh-poly.htm#TOC"
  >-------------------------------</a><a name="furthestout">-</a>
  
  qh_furthestout( facet )
    make furthest outside point the last point of outsideset

  returns:
    updates facet->outsideset
    clears facet->notfurthest
    sets facet->furthestdist

  design:
    determine best point of outsideset
    make it the last point of outsideset
*/
void qh_furthestout (facetT *facet) {
  pointT *point, **pointp, *bestpoint= NULL;
  realT dist, bestdist= -REALmax;

  FOREACHpoint_(facet->outsideset) {
    qh_distplane (point, facet, &dist);
    zinc_(Zcomputefurthest);
    if (dist > bestdist) {
      bestpoint= point;
      bestdist= dist;
    }
  }
  if (bestpoint) {
    qh_setdel (facet->outsideset, point);
    qh_setappend (&facet->outsideset, point);
#if !qh_COMPUTEfurthest
    facet->furthestdist= bestdist;
#endif
  }
  facet->notfurthest= False;
  trace3((qh ferr, "qh_furthestout: p%d is furthest outside point of f%d\n",
	  qh_pointid (point), facet->id));
} /* furthestout */


/*-<a                             href="qh-qhull.htm#TOC"
  >-------------------------------</a><a name="infiniteloop">-</a>
  
  qh_infiniteloop( facet )
    report infinite loop error due to facet
*/
void qh_infiniteloop (facetT *facet) {

  fprintf (qh ferr, "qhull internal error (qh_infiniteloop): potential infinite loop detected\n");
  qh_errexit (qh_ERRqhull, facet, NULL);
} /* qh_infiniteloop */

/*-<a                             href="qh-poly.htm#TOC"
  >-------------------------------</a><a name="initbuild">-</a>
  
  qh_initbuild()
    initialize hull and outside sets with point array
    qh.FIRSTpoint/qh.NUMpoints is point array
    if qh.GOODpoint
      adds qh.GOODpoint to initial hull

  returns:
    qh_facetlist with initial hull
    points partioned into outside sets, coplanar sets, or inside
    initializes qh.GOODpointp, qh.GOODvertexp,

  design:
    initialize global variables used during qh_buildhull
    determine precision constants and points with max/min coordinate values
      if qh.SCALElast, scale last coordinate (for 'd')
    build initial simplex
    partition input points into facets of initial simplex
    set up lists
    if qh.ONLYgood
      check consistency  
      add qh.GOODvertex if defined
*/
void qh_initbuild( void) {
  setT *maxpoints, *vertices;
  facetT *facet;
  int i, numpart;
  realT dist;
  boolT isoutside;

  qh furthest_id= -1;
  qh lastreport= 0;
  qh facet_id= qh vertex_id= qh ridge_id= 0;
  qh visit_id= qh vertex_visit= 0;
  qh maxoutdone= False;

  if (qh GOODpoint > 0) 
    qh GOODpointp= qh_point (qh GOODpoint-1);
  else if (qh GOODpoint < 0) 
    qh GOODpointp= qh_point (-qh GOODpoint-1);
  if (qh GOODvertex > 0)
    qh GOODvertexp= qh_point (qh GOODvertex-1);
  else if (qh GOODvertex < 0) 
    qh GOODvertexp= qh_point (-qh GOODvertex-1);
  if ((qh GOODpoint  
       && (qh GOODpointp < qh first_point  /* also catches !GOODpointp */
	   || qh GOODpointp > qh_point (qh num_points-1)))
    || (qh GOODvertex
	&& (qh GOODvertexp < qh first_point  /* also catches !GOODvertexp */
	    || qh GOODvertexp > qh_point (qh num_points-1)))) {
    fprintf (qh ferr, "qhull input error: either QGn or QVn point is > p%d\n",
	     qh num_points-1);
    qh_errexit (qh_ERRinput, NULL, NULL);
  }
  maxpoints= qh_maxmin(qh first_point, qh num_points, qh hull_dim);
  if (qh SCALElast)
    qh_scalelast (qh first_point, qh num_points, qh hull_dim,
               qh MINlastcoord, qh MAXlastcoord, qh MAXwidth);
  qh_detroundoff();
  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) {
      if (qh UPPERdelaunay) { /* matches qh.upperdelaunay in qh_setfacetplane */
	qh lower_threshold[qh hull_dim-1]= qh ANGLEround * qh_ZEROdelaunay;
	qh GOODthreshold= True;
      }else { 
	qh upper_threshold[qh hull_dim-1]= -qh ANGLEround * qh_ZEROdelaunay;
        if (!qh GOODthreshold) 
	  qh SPLITthresholds= True; /* build upper-convex hull even if Qg */
          /* qh_initqhull_globals errors if Qg without Pdk/etc. */
      }
    }
  }
  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 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);
  }
  qh_resetlists (False, qh_RESETvisible /*qh visible_list newvertex_list newfacet_list */);
  qh facet_next= qh facet_list;
  qh_furthestnext (/* qh facet_list */);
  if (qh PREmerge) {
    qh cos_max= qh premerge_cos;
    qh centrum_radius= qh premerge_centrum;
  }
  if (qh ONLYgood) {
    if (qh GOODvertex > 0 && qh MERGING) {
      fprintf (qh ferr, "qhull input error: 'Qg QVn' (only good vertex) does not work with merging.\nUse 'QJ' to joggle the input or 'Q0' to turn off merging.\n");
      qh_errexit (qh_ERRinput, NULL, NULL);
    }
    if (!(qh GOODthreshold || qh GOODpoint
         || (!qh MERGEexact && !qh PREmerge && qh GOODvertexp))) {
      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 with 'QJ' or 'Q0' (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, !qh_ALL, &isoutside, &numpart);
      zadd_(Zdistgood, numpart);
      if (!isoutside) {
        fprintf (qh ferr, "qhull input error: point for QV%d is inside initial simplex.  It can not be made a vertex.\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);
  trace1((qh ferr, "qh_initbuild: initial hull created and points partitioned\n"));
} /* initbuild */

/*-<a                             href="qh-poly.htm#TOC"
  >-------------------------------</a><a name="initialhull">-</a>
  
  qh_initialhull( vertices )
    constructs the initial hull as a DIM3 simplex of vertices

  design:
    creates a simplex (initializes lists)
    determines orientation of simplex
    sets hyperplanes for facets
    doubles checks orientation (in case of axis-parallel facets with Gaussian elimination)
    checks for flipped facets and qh.NARROWhull
    checks the result   
*/
void qh_initialhull(setT *vertices) {
  facetT *facet, *firstfacet, *neighbor, **neighborp;
  realT dist, angle, minangle= REALmax;
#ifndef qh_NOtrace
  int k;
#endif

  qh_createsimplex(vertices);  /* qh facet_list */
  qh_resetlists (False, qh_RESETvisible);
  qh facet_next= qh facet_list;      /* advance facet when processed */
  qh interior_point= qh_getcenter(vertices);
  firstfacet= qh facet_list;
  qh_setfacetplane(firstfacet);
  zinc_(Znumvisibility); /* needs to be in printsummary */
  qh_distplane(qh interior_point, firstfacet, &dist);
  if (dist > 0) {  
    FORALLfacets
      facet->toporient ^= True;
  }
  FORALLfacets
    qh_setfacetplane(facet);
  FORALLfacets {
    if (!qh_checkflipped (facet, NULL, qh_ALL)) {/* due to axis-parallel facet */
      trace1((qh ferr, "qh_initialhull: initial orientation incorrect.  Correct all facets\n"));
      facet->flipped= False;
      FORALLfacets {
	facet->toporient ^= True;
	qh_orientoutside (facet);
      }
      break;
    }
  }
  FORALLfacets {
    if (!qh_checkflipped (facet, NULL, !qh_ALL)) {  /* can happen with 'R0.1' */
      qh_precision ("initial facet is coplanar with interior point");
      fprintf (qh ferr, "qhull precision error: initial facet %d is coplanar with the interior point\n",
                   facet->id);
      qh_errexit (qh_ERRsingular, facet, NULL);
    }
    FOREACHneighbor_(facet) {
      angle= qh_getangle (facet->normal, neighbor->normal);
      minimize_( minangle, angle);
    }
  }
  if (minangle < qh_MAXnarrow && !qh NOnarrow) { 
    realT diff= 1.0 + minangle;

    qh NARROWhull= True;
    qh_option ("_narrow-hull", NULL, &diff);
    if (minangle < qh_WARNnarrow && !qh RERUN && qh PRINTprecision)
      fprintf (qh ferr, "qhull precision warning: \n\
The initial hull is narrow (cosine of min. angle is %.16f).\n\
A coplanar point may lead to a wide facet.  Options 'QbB' (scale to unit box)\n\
or 'Qbb' (scale last coordinate) may remove this warning.  Use 'Pp' to skip\n\
this warning.  See 'Limitations' in qh-impre.htm.\n",
          -minangle);   /* convert from angle between normals to angle between facets */
  }
  zzval_(Zprocessed)= qh hull_dim+1;
  qh_checkpolygon (qh facet_list);
  qh_checkconvex(qh facet_list,   qh_DATAfault);
#ifndef qh_NOtrace
  if (qh IStracing >= 1) {
    fprintf(qh ferr, "qh_initialhull: simplex constructed, interior point:");
    for (k=0; k < qh hull_dim; k++) 
      fprintf (qh ferr, " %6.4g", qh interior_point[k]);
    fprintf (qh ferr, "\n");
  }
#endif
} /* initialhull */

/*-<a                             href="qh-poly.htm#TOC"
  >-------------------------------</a><a name="initialvertices">-</a>
  
  qh_initialvertices( dim, maxpoints, points, numpoints )
    determines a non-singular set of initial vertices
    maxpoints may include duplicate points

  returns:
    temporary set of dim+1 vertices in descending order by vertex id
    if qh.RANDOMoutside && !qh.ALLpoints
      picks random points
    if dim >= qh_INITIALmax, 
      uses min/max x and max points with non-zero determinants

  notes:
    unless qh.ALLpoints, 
      uses maxpoints as long as determinate is non-zero
*/
setT *qh_initialvertices(int dim, setT *maxpoints, pointT *points, int numpoints) {
  pointT *point, **pointp;
  setT *vertices, *simplex, *tested;
  realT randr;
  int index, point_i, point_n, k;
  boolT nearzero= False;
  
  vertices= qh_settemp (dim + 1);
  simplex= qh_settemp (dim+1);
  if (qh ALLpoints) 
    qh_maxsimplex (dim, NULL, points, numpoints, &simplex);
  else if (qh RANDOMoutside) {
    while (qh_setsize (simplex) != dim+1) {
      randr= qh_RANDOMint;
      randr= randr/(qh_RANDOMmax+1);
      index= (int