qhull_poly2.cxx

hl2 source code. Do not use it illegal.

  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)))) {
    ivp_message( "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)
      ivp_message( "\nTrace level %d for %s | %s\n", 
         qh IStracing ? qh IStracing : qh TRACElevel, qh rbox_command, qh qhull_command);
    ivp_message( "Options selected for qhull %s:\n%s\n", qh_version, qh qhull_options);
  }
  qh_resetlists (False /*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) {
      ivp_message( "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))) {
      ivp_message( "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, &isoutside, &numpart);
      zadd_(Zdistgood, numpart);
      if (!isoutside) {
        ivp_message( "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-c.htm#poly"
  >-------------------------------</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 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");
      if (qh ferr) {ivp_message( "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) {
    realT diff= 1.0 + minangle;

    qh NARROWhull= True;
    qh_option ("_narrow-hull", NULL, &diff);
    if (minangle < qh_WARNnarrow && !qh RERUN && qh PRINTprecision)
      ivp_message( "qhull precision warning: \n\
The initial hull is narrow (the cosine of the minimum angle is %.9g).\n\
A coplanar point may lead to a wide facet.  Options 'Qs' (search for best\n\
initial hull), 'QbB' (scale to unit box), or 'Qbb' (scale last coordinate)\n\
may remove this warning.  Use 'Pp' to ignore this warning.\n\
See 'Limitations' in qh-impre.htm.\n",
          minangle);
  }
  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) {
    ivp_message( "qh_initialhull: simplex constructed, interior point:");
    for (k=0; k < qh hull_dim; k++) 
      ivp_message( " %6.4g", qh interior_point[k]);
    ivp_message( "\n");
  }
#endif
} /* initialhull */

/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</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)floor(qh num_points * randr);
      while (qh_setin (simplex, qh_point (index))) {
	index++; /* in case qh_RANDOMint always returns the same value */
        index= index < qh num_points ? index : 0;
      }
      qh_setappend (&simplex, qh_point (index));
    }
  }else if (qh hull_dim >= qh_INITIALmax) {
    tested= qh_settemp (dim+1);
    qh_setappend (&simplex, SETfirst_(maxpoints));   /* max and min X coord */
    qh_setappend (&simplex, SETsecond_(maxpoints));
    qh_maxsimplex (fmin_(qh_INITIALsearch, dim), maxpoints, points, numpoints, &simplex);
    k= qh_setsize (simplex);
    FOREACHpoint_i_(maxpoints) { 
      if (point_i & 0x1) {     /* first pick up max. coord. points */
      	if (!qh_setin (simplex, point) && !qh_setin (tested, point)){
	  qh_detsimplex(point, simplex, k, &nearzero);
          if (nearzero)
            qh_setappend (&tested, point);
          else {
            qh_setappend (&simplex, point);
            if (++k == dim)  /* use search for last point */
	      break;
	  }
	}
      }
    }
    while (k != dim && (point= (pointT*)qh_setdellast (maxpoints))) {
      if (!qh_setin (simplex, point) && !qh_setin (tested, point)){
        qh_detsimplex (point, simplex, k, &nearzero);
        if (nearzero)
          qh_setappend (&tested, point);
        else {
          qh_setappend (&simplex, point);
          k++;
	}
      }
    }
    index= 0;
    while (k != dim && (point= qh_point (index++))) {
      if (!qh_setin (simplex, point) && !qh_setin (tested, point)){
        qh_detsimplex (point, simplex, k, &nearzero);
        if (!nearzero){
          qh_setappend (&simplex, point);
          k++;
	}
      }
    }
    qh_settempfree (&tested);
    qh_maxsimplex (dim, maxpoints, points, numpoints, &simplex);
  }else
    qh_maxsimplex (dim, maxpoints, points, numpoints, &simplex);
  FOREACHpoint_(simplex) 
    qh_setaddnth (&vertices, 0, qh_newvertex(point)); /* descending order */
  qh_settempfree (&simplex);
  return vertices;
} /* initialvertices */


/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="isvertex">-</a>
  
  qh_isvertex(  )
    returns vertex if point is in vertex set, else returns NULL

  notes:
    for qh.GOODvertex
*/
vertexT *qh_isvertex (pointT *point, setT *vertices) {
  vertexT *vertex, **vertexp;

  FOREACHvertex_(vertices) {
    if (vertex->point == point)
      return vertex;
  }
  return NULL;
} /* isvertex */

/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="makenewfacets">-</a>
  
  qh_makenewfacets( point )
    make new facets from point and qh.visible_list

  returns:
    qh.newfacet_list= list of new facets with hyperplanes and ->newfacet
    qh.newvertex_list= list of vertices in new facets with ->newlist set
    
    if (qh.ONLYgood)
      newfacets reference horizon facets, but not vice versa
      ridges reference non-simplicial horizon ridges, but not vice versa
      does not change existing facets
    else
      sets qh.NEWfacets
      new facets attached to horizon facets and ridges
      for visible facets, 
        visible->r.replace is corresponding new facet

  design:
    for each visible facet
      make new facets to its horizon facets
      update its f.replace 
      clear its neighbor set
*/
vertexT *qh_makenewfacets (pointT *point /*visible_list*/) {
  facetT *visible, *newfacet= NULL, *newfacet2= NULL, *neighbor, **neighborp;
  vertexT *apex;
  int numnew=0;

  qh newfacet_list= qh facet_tail;
  qh newvertex_list= qh vertex_tail;
  apex= qh_newvertex(point);
  qh_appendvertex (apex);  
  qh visit_id++;
  if (!qh ONLYgood)
    qh NEWfacets= True;
  FORALLvisible_facets {
    FOREACHneighbor_(visible) 
      neighbor->seen= False;
    if (visible->ridges) {
      visible->visitid= qh visit_id;
      newfacet2= qh_makenew_nonsimplicial (visible, apex, &numnew);
    }
    if (visible->simplicial)
      newfacet= qh_makenew_simplicial (visible, apex, &numnew);
    if (!qh ONLYgood) {
      if (newfacet2)  /* newfacet is null if all ridges defined */
        newfacet= newfacet2;
      if (newfacet)
      	visible->f.replace= newfacet;
      else
        zinc_(Zinsidevisible);
      SETfirst_(visible->neighbors)= NULL;
    }
  }
  trace1((qh ferr, "qh_makenewfacets: created %d new facets from point p%d to horizon\n",
	  numnew, qh_pointid(p