pelqhull.cc

Gambit 是一个游戏库理论软件

  }
  *newp= NULL;
  return(newset);
} /* setnew_delnthsorted */


/*----------------------------------------
-setprint- print set elements to fp
notes:
  never errors
*/
void qh_setprint(FILE *fp, char* string, setT *set) {
  int size, k;

  if (!set)
    fprintf (fp, "%s set is null\n", string);
  else {
    SETreturnsize_(set, size);
    fprintf (fp, "%s set=%x maxsize=%d size=%d elems=",
	     string, (unsigned int)set, set->maxsize, size);
    if (size > (int)set->maxsize)
      size= set->maxsize+1;
    for (k=0; k<size; k++)
      fprintf(fp, " %x", (unsigned int)(set->e[k]));
    fprintf(fp, "\n");
  }
} /* setprint */

/*----------------------------------------
-setreplace- replaces oldelem in set with newelem
   errors if oldelem not in the set
   if newelem is NULL then FOREACH no longer works
*/
void qh_setreplace(setT *set, void *oldelem, void *newelem) {
  void **elemp;
  
  elemp= SETaddr_(set, void);
  while(*elemp != oldelem && *elemp)
    elemp++;
  if (*elemp)
    *elemp= newelem;

  else qhull_fatal(15);

} /* setreplace */


/*----------------------------------------
-setsize- returns the size of a set
  same as SETreturnsize_(set)
*/
int qh_setsize(setT *set) {
  int size, *sizep;
  
  if (!set)
    return (0);
  sizep= SETsizeaddr_(set);
  if ((size= *sizep)) {
    size--;

    if (size > (int)set->maxsize) qhull_fatal(16);

  }else
    size= set->maxsize;
  return size;
} /* setsize */

/*----------------------------------------
-settemp- return a stacked, temporary set
  use settempfree or settempfree_all to release from qhmem.tempstack
  see also qh_setnew
*/
setT *qh_settemp(int setsize) {
  setT *newset;
  
  newset= qh_setnew (setsize);
  qh_setappend ((setT **)&qhmem.tempstack, newset);
  if (qhmem.IStracing >= 5)
    fprintf (qhmem.ferr, "qh_settemp: temp set %x of %d elements, depth %d\n",
       (int)newset, newset->maxsize, qh_setsize ((setT *)qhmem.tempstack));
  return newset;
} /* settemp */

/*----------------------------------------
-settempfree- free temporary set at top of qhmem.tempstack
  nop if NULL
  errors if set not from previous qh_settemp
    locate source by T2 and find mis-matching qh_settemp
*/
void qh_settempfree(setT **set) {
  setT *stackedset;

  if (!*set)
    return;
  stackedset= qh_settemppop ();

  if (stackedset != *set) qhull_fatal(17);

  qh_setfree (set);
} /* settempfree */

/*----------------------------------------
-settempfree_all- free all temporary sets in qhmem.tempstack
*/
void qh_settempfree_all(void) {
  setT *set, **setp;

  FOREACHset_((setT *)qhmem.tempstack) 
    qh_setfree(&set);
  qh_setfree((setT **)&qhmem.tempstack);
} /* settempfree_all */

/*----------------------------------------
-settemppop- pop and return temporary set from qhmem.tempstack (makes it permanent)
*/
setT *qh_settemppop(void) {
  setT *stackedset;
  
  stackedset= (setT *)qh_setdellast((setT *)qhmem.tempstack);

  if (!stackedset) qhull_fatal(18);

  if (qhmem.IStracing >= 5)
    fprintf (qhmem.ferr, "qh_settemppop: depth %d temp set %x of %d elements\n",
       qh_setsize((setT *)qhmem.tempstack)+1, (int)stackedset, qh_setsize(stackedset));
  return stackedset;
} /* settemppop */

/*----------------------------------------
-settemppush- push temporary set unto qhmem.tempstack (makes it temporary)
  duplicates settemp() for tracing
*/
void qh_settemppush(setT *set) {
  
  qh_setappend ((setT**)&qhmem.tempstack, set);
  if (qhmem.IStracing >= 5)
    fprintf (qhmem.ferr, "qh_settemppush: depth %d temp set %x of %d elements\n",
    qh_setsize((setT *)qhmem.tempstack), (int)set, qh_setsize(set));
} /* settemppush */

 
/*----------------------------------------
-settruncate- truncate set to size elements
  set must be defined
*/
void qh_settruncate (setT *set, int size) {

  if (size < 0 || size > (int)set->maxsize) qhull_fatal(19);

  set->e[set->maxsize]= (void *) (size+1);   /* maybe overwritten */
  set->e[size]= NULL;
} /* setruncate */
    
/*----------------------------------------
-setunique- add element if it isn't already
  returns 1 if it's appended
*/
int qh_setunique (setT **set, void *elem) {

  if (!qh_setin (*set, elem)) {
    qh_setappend (set, elem);
    return 1;
  }
  return 0;
} /* setunique */
    
/*----------------------------------------
-setzero- zero remainder of set and set its size
  set must be defined
*/
void qh_setzero (setT *set, int index, int size) {
  int count;

  if (index < 0 || index >= size || size > (int)set->maxsize) qhull_fatal(20);

  (set->e[set->maxsize])=  (void *)(size+1);  /* may be overwritten */
  count= size - index + 1;   /* +1 for NULL terminator */
  memset ((char *)SETelemaddr_(set, index, void), 0, count * sizeof(void *));
} /* setzero */

    
/*************************************************************************/
/****************** implementation code from geom.c **********************/
/*************************************************************************/


/* geom.c -- geometric routines of qhull
   
   see README and geom.h

   copyright (c) 1993-1994 The Geometry Center        
*/
   
/*-------------------------------------------------
-backnormal- solve for normal x using back substitution over rows U
     solves Ux=b where Ax=b and PA=LU
     b= [0,...,0,sign or 0]  (-1 if sign, else +1)
     last row of A= [0,...,0,1]
     assumes numrow == numcol-1
returns:
     normal= x 
     if can't divzero() for later normalization (qh MINdenom_2 and qh MINdenom_1_2),
         sets tail of normal to [...,sign,0,...], i.e., solves for b= [0...]
	 sets nearzero, unless last row (i.e., hyperplane intersects [0,..,1])
notes:
     1) Ly=Pb == y=b since P only permutes the 0's of b
     see Golub & van Loan 4.4-9 for back substitution
*/
void qh_backnormal (realT **rows, int numrow, int numcol, boolT sign,
  	coordT *normal, boolT *nearzero) {
  int i, j;
  coordT *normalp, *normal_tail, *ai, *ak;
  realT diagonal;
  boolT waszero;
  int zerocol=-1;
  
  normalp= normal + numcol - 1;
  *normalp--= (sign ? -1.0 : 1.0);
  for(i= numrow; i--; ) {
    *normalp= 0.0;
    ai= rows[i] + i + 1;
    ak= normalp+1;
    for(j= i+1; j < numcol; j++)
      *normalp -= *ai++ * *ak++;
    diagonal= (rows[i])[i];
    if (fabs_(diagonal) > qh MINdenom_2)
      *(normalp--) /= diagonal;
    else {
      waszero= False;
      *normalp= qh_divzero (*normalp, diagonal, qh MINdenom_1_2, &waszero);
      if (waszero) {
        zerocol= i;
	*(normalp--)= (sign ? -1.0 : 1.0);
	for (normal_tail= normalp+2; normal_tail < normal + numcol; normal_tail++)
	  *normal_tail= 0.0;
      }else
	normalp--;
    }
  }
  if (zerocol != -1) {
    zzinc_(Zback0);
    *nearzero= True;
    trace4((qh ferr, "qh_backnormal: zero diagonal at column %d.\n", i));
  }
} /* backnormal */

/*-------------------------------------------------
-crossproduct- of 2 dim vectors, C= A x B
    from Glasner, Graphics Gems I, p. 639
    NOTE: only defined for dim==3
*/
void qh_crossproduct (int dim, realT vecA[3], realT vecB[3], realT vecC[3]){

  if (dim == 3) {
    vecC[0]=   det2_(vecA[1], vecA[2],
		     vecB[1], vecB[2]);
    vecC[1]= - det2_(vecA[0], vecA[2],
		     vecB[0], vecB[2]);
    vecC[2]=   det2_(vecA[0], vecA[1],
		     vecB[0], vecB[1]);
  }
} /* vcross */

/*-------------------------------------------------
-determinant- compute the determinant of a square matrix
  rows= row vectors
  uses qh NEARzero to test for degenerate matrices
    this does look right, probably no easy way of doing it
returns:
  determinant
  overwrites rows and the matrix
  nearzero set if degenerate
*/
realT qh_determinant (realT **rows, int dim, boolT *nearzero) {
  realT det=0;
  int i;
  boolT sign= False;

  *nearzero= False;
  if (dim < 2) { qhull_fatal(21); }
  else if (dim == 2) {
    det= det2_(rows[0][0], rows[0][1],
		 rows[1][0], rows[1][1]);
    if (fabs_(det) < qh NEARzero[1])  /* not really correct, what should this be? */
      *nearzero= True;
  }else if (dim == 3) {
    det= det3_(rows[0][0], rows[0][1], rows[0][2],
		 rows[1][0], rows[1][1], rows[1][2],
		 rows[2][0], rows[2][1], rows[2][2]);
    if (fabs_(det) < qh NEARzero[2])  /* not really correct, what should this be? */
      *nearzero= True;
  }else {	
    qh_gausselim(rows, dim, dim, &sign, nearzero);  /* if nearzero, diagonal still ok*/
    det= 1.0;
    for (i= dim; i--; )
      det *= (rows[i])[i];
    if (sign)
      det= -det;
  }
  return det;
} /* determinant */

/*-------------------------------------------------
-detsimplex- compute determinant of a simplex with point apex and base points
   uses qh gm_matrix/qh gm_row (assumes they're big enough)
   uses dim coordinates of point and vertex->point
returns:
   if dim == 2 or 3
     nearzero iff determinant < qh NEARzero[dim-1]  (not quite correct)
   if dim >= 4
     nearzero iff diagonal[k] < qh NEARzero[k]
*/
realT qh_detsimplex(pointT *apex, setT *points, int dim, boolT *nearzero) {
  pointT *coorda, *coordp, *gmcoord, *point, **pointp;
  coordT **rows;
  int k,  i=0;
  realT det= 0.0;

  zinc_(Zdetsimplex);
  gmcoord= qh gm_matrix;
  rows= qh gm_row;
  FOREACHpoint_(points) {
    if (i == dim)
      break;
    rows[i++]= gmcoord;
    coordp= point;
    coorda= apex;
    for(k= dim; k--; )
      *(gmcoord++)= *coordp++ - *coorda++;
  }
  if (i < dim) { qhull_fatal(22); }
  det= qh_determinant (rows, dim, nearzero);
  trace2((qh ferr, "qh_detsimplex: det=%2.2g for point p%d, dimension %d, nearzero? %d\n",
	  det, qh_pointid(apex), dim, *nearzero)); 
  return det;
} /* detsimplex */

/*-------------------------------------------
-distplane- get distance from point to facet
returns:
    positive if point is above facet (i.e., outside)
    can not qhull_fatal (for sortfacets)
*/
void qh_distplane (pointT *point, facetT *facet, realT *dist) {
  coordT *normal= facet->normal, *coordp, randr;
  int k;
  
  switch(qh hull_dim){
  case 2:
    *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1];
    break;
  case 3:
    *dist= facet->offset + point[0] * normal[0] + point[1] * normal[1] + point[2] * normal[2];
    break;
  case 4:
    *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3];
    break;
  case 5:
    *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4];
    break;
  case 6:
    *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5];
    break;
  case 7:  
    *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5]+point[6]*normal[6];
    break;
  case 8:
    *dist= facet->offset+point[0]*normal[0]+point[1]*normal[1]+point[2]*normal[2]+point[3]*normal[3]+point[4]*normal[4]+point[5]*normal[5]+point[6]*normal[6]+point[7]*normal[7];
    break;
  default:
    *dist= facet->offset;
    coordp= point;
    for (k= qh hull_dim; k--; )
      *dist += *coordp++ * *normal++;
    break;
  }
  zinc_(Zdistplane);
  if (!qh RANDOMdist && qh IStracing < 4)
    return;
  if (qh RANDOMdist) {
    randr= qh_RANDOMint;
    *dist += (2.0 * randr / qh_RANDOMmax - 1.0) *
      qh RANDOMfactor * qh maxmaxcoord;
  }
  /*
  if (qh IStracing >= 4) {
    fprintf (qh ferr, "qh_distplane: ");
    fprintf (qh ferr, qh_REAL_1, *dist);
    fprintf (qh ferr, "from p%d to f%d\n", qh_pointid(point), facet->id);
  }
  */
  return;
} /* distplane */


/*--------------------------------------------------
-divzero -- divide by a number that's nearly zero
  mindenom1= minimum denominator for dividing into 1.0
returns:
  zerodiv and 0.0 if it would overflow
*/
realT qh_divzero (realT numer, realT denom, realT mindenom1, boolT *zerodiv) {
  realT temp, numerx, denomx;
  

  if (numer < mindenom1 && numer > -mindenom1) {
    numerx= fabs_(numer);
    denomx= fabs_(denom);
    if (numerx < denomx) {
      *zerodiv= False;
      return numer/denom;
    }else {
      *zerodiv= True;
      return 0.0;
    }
  }
  temp= denom/numer;
  if (temp > mindenom1 || temp < -mindenom1) {
    *zerodiv= False;
    return numer/denom;
  }else {
    *zerodiv= True;
    return 0.0;
  }
} /* divzero */
  

/*-------------------------------------------------
-facetcenter- return Voronoi center for a facet's vertices
*/
pointT *qh_facetcenter (setT *vertices) {
  setT *points= qh_settemp (qh_setsize (vertices));
  vertexT *vertex, **vertexp;
  pointT *center;
  
  FOREACHvertex_(vertices) 
    qh_setappend (&points, vertex->point);
  center= qh_voronoi_center (qh hull_dim-1, points);
  qh_settempfree (&points);
  return center;
} /* facetcenter */

/*-------------------------------------------------
-findbest- find best facet for point starting at a facet (not flipped!)
  if bestoutside, searches all facets else stops at first outside
    MINoutside is DISTround in precise case
  if firstid, searches facets with ids >= firstid
    searches old facets if bestoutside || (not outside and imprecise)
  searches all neighbors of coplanar and flipped facets
    searchdist is arbitrarily set to min_vertex+max_outside+DISTround
       max_outside is needed for setting facet->maxoutside
returns:
  if !firstid, updates facet->maxoutside for good, visited facets
  distance to facet
  isoutside true if point is outside of facet
  bumps visit_id and seen flags
notes:
  uses visitid and seen
  statistics collected here for partitions, caller does outside/coplanar
  caller traces the results
  #2 after setfacetplane in D3, optimized for outside points and !bestoutside