pelqhull.cc

Gambit 是一个游戏库理论软件

  #1 when merging in D3
  see also partitionall()
notes on searchdist:
  searchdist needed since vertex neighbors can be geometric neighbors of facet 
  if searchdist=DISTround, gets stuck for rbox 50 W1e-3 D7 | qhull A-0.99 W0.2
  if !BESToutside and merging, gets stuck for rbox 1000 W8e-6 | qhull C-0
    because nearly coplanar widens when the point is outside of the facets
  searching all new facets does not prevent !BESToutside getting stuck
  check_maxoutside can also get stuck, should keep coplanars
*/
facetT *qh_findbest (pointT *point, facetT *facet, boolT bestoutside,
	   unsigned firstid, realT *dist, boolT *isoutside, int *numpart) {
  realT bestdist, searchdist;
  facetT *neighbor, **neighborp, *bestfacet;
  setT *search= NULL;
  int oldtrace= qh IStracing;
  boolT checkmax= (boolT)(bestoutside && !firstid && qh_MAXoutside
    && (qh MERGING || qh APPROXhull));

  if (qh TRACEpoint >= 0 && qh TRACEpoint == qh_pointid (point)) {
    qh IStracing= qh TRACElevel;
    /*
    fprintf (qh ferr, "qh_findbest: point p%d starting at f%d bestoutside? %d firstid %d\n",
	     qh TRACEpoint, facet->id, bestoutside, firstid);
    fprintf (qh ferr, "  Last point added to hull was p%d.", qh furthest_id);
    fprintf(qh ferr, "  Last merge was #%d.\n", zzval_(Ztotmerge));
    */
  }
  searchdist= - qh min_vertex + qh max_outside + 2* qh DISTround;
  *isoutside= True;
  *numpart= 1;
  qh_distplane (point, facet, dist);
  bestdist= *dist;
  bestfacet= facet;
  if (!bestoutside &&  *dist >= qh MINoutside) 
    goto LABELreturn_best;
#if qh_MAXoutside
  if (checkmax && (!qh ONLYgood || facet->good) && *dist > facet->maxoutside)
    facet->maxoutside= *dist;
#endif
  facet->visitid= ++qh visit_id;
  facet->seen= False;
  if (True) {                   /* directed search for bestfacet */
LABELrepeat:			   /* facet->seen if clearly worse */
    trace4((qh ferr, "qh_findbest: neighbors of f%d\n", facet->id));
    FOREACHneighbor_(facet) {
      if ((int)neighbor->visitid == qh visit_id)
        continue;
      if (neighbor->id < firstid) {
        neighbor->seen= True;
        continue;
      }
      neighbor->visitid= qh visit_id;
      neighbor->seen= False;
      if (neighbor->flipped)
        continue;
      (*numpart)++;
      qh_distplane (point, neighbor, dist);
      if (!bestoutside && *dist >= qh MINoutside) {
	bestfacet= neighbor;
	goto LABELreturn_best;
      }
#if qh_MAXoutside
      if (checkmax && (!qh ONLYgood || neighbor->good) 
          && *dist > neighbor->maxoutside)
        neighbor->maxoutside= *dist;
#endif
      if (*dist >= bestdist) {  /* >= for exact coplanar */
        bestdist= *dist;
        bestfacet= neighbor;
        if (*dist > bestdist + searchdist)
          facet->seen= True;
        facet= neighbor;
        goto LABELrepeat;
      }else if (*dist < bestdist - searchdist)
        neighbor->seen= True;
    }
  }
  do {                   /* search horizon of facet */
    FOREACHneighbor_(facet) {
      if ((int)neighbor->visitid == qh visit_id) {
        if (!neighbor->seen) {
          neighbor->seen= True;
          if (!search)
            search= qh_settemp (qh TEMPsize);
          qh_setappend (&search, neighbor);
        }
        continue;
      }
      neighbor->visitid= qh visit_id;
      neighbor->seen= True;
      if (neighbor->flipped) {
        if (!search)
          search= qh_settemp (qh TEMPsize);
        qh_setappend (&search, neighbor);
        continue;
      }
      if (neighbor->id < firstid) {
        if (!(bestoutside+qh APPROXhull+qh PREmerge))
          continue;
      }else
        zinc_(Zpartneighbor);
      (*numpart)++;
      qh_distplane (point, neighbor, dist);
      if (!bestoutside && *dist >= qh MINoutside) {
	bestfacet= neighbor;
	goto LABELreturn_best;
      }
#if qh_MAXoutside
      if (checkmax && *dist > neighbor->maxoutside)
        neighbor->maxoutside= *dist;
#endif
      if (*dist >= bestdist - searchdist) {
        if (!search)
          search= qh_settemp (qh TEMPsize);
        qh_setappend (&search, neighbor);
        if (*dist > bestdist) {
          bestdist= *dist;
          bestfacet= neighbor;
        }
      }
    }
  }while ((facet= (facetT *)qh_setdellast (search)));
  *dist= bestdist;
  if (!bestoutside || bestdist < qh MINoutside)
    *isoutside= False;
LABELreturn_best:
  if (search)
    qh_settempfree (&search);
  qh IStracing= oldtrace;
  return bestfacet;
}  /* findbest */


/*-------------------------------------------------
-findgooddist- find best good facet visible for point from facetA
  assumes facetA is visible from point
  uses qh visit_id and qh visible_list (but doesn't set visible)
returns:
  furthest distance to good facet, if any
  bumps visit_id and seen flags
*/
facetT *qh_findgooddist (pointT *point, facetT *facetA, realT *distp) {
  realT bestdist= -REALmax, dist;
  facetT *neighbor, **neighborp, *bestfacet=NULL, *facet;
  boolT goodseen= False;  

  if (facetA->good) {
    zinc_(Zverifypart);
    qh_distplane (point, facetA, &bestdist);
    bestfacet= facetA;
    goodseen= True;
  }
  qh_removefacet (facetA);
  qh_appendfacet (facetA);
  qh visible_list= facetA;
  facetA->visitid= ++qh visit_id;
  FORALLfacet_(qh visible_list) {
    FOREACHneighbor_(facet) {
      if ((int)neighbor->visitid == qh visit_id)
        continue;
      neighbor->visitid= qh visit_id;
      if (goodseen && !neighbor->good)
        continue;
      zinc_(Zverifypart);
      qh_distplane (point, neighbor, &dist);
      if (dist > 0) {
        qh_removefacet (neighbor);
        qh_appendfacet (neighbor);
        if (neighbor->good) {
          goodseen= True;
          if (dist > bestdist) {
            bestdist= dist;
            bestfacet= neighbor;
          }
        }
      }
    }
  }
  if (bestfacet) {
    *distp= bestdist;
    trace2((qh ferr, "qh_findgooddist: p%d is %2.2g above good facet f%d\n",
      qh_pointid(point), bestdist, bestfacet->id));
    return bestfacet;
  }
  trace4((qh ferr, "qh_findgooddist: no good facet for p%d above f%d\n", 
      qh_pointid(point), facetA->id));
  return NULL;
}  /* findgooddist */
    
/*-------------------------------------------------
-gausselim- Gaussian elimination with partial pivoting
  coordT data in rows
  assumes numrow <= numcol
returns:
  rows is upper triangular (includes row exchanges)
  flips sign for each row exchange
  sets nearzero if pivot[k] < qh NEARzero[k], else False.
    if nearzero, the determinant's sign may be incorrect.
*/
void qh_gausselim(realT **rows, int numrow, int numcol, boolT *sign, boolT *nearzero) {
  realT *ai, *ak, *rowp, *pivotrow;
  realT n, pivot, pivot_abs= 0.0, temp;
  int i, j, k, pivoti, flip=0, tempint;
  
  *nearzero= False;
  for(k= 0; k < numrow; k++) {
    pivot_abs= fabs_((rows[k])[k]);
    pivoti= k;
    for(i= k+1; i < numrow; i++) {
      if ((temp= fabs_((rows[i])[k])) > pivot_abs) {
	pivot_abs= temp;
	pivoti= i;
      }
    }
    if (pivoti != k) {
      rowp= rows[pivoti]; 
      rows[pivoti]= rows[k]; 
      rows[k]= rowp; 
      tempint = (int)*sign;
      tempint ^= 1;
      *sign = (boolT)tempint;
      flip ^= 1;
    }
    if (pivot_abs <= qh NEARzero[k]) {
      *nearzero= True;
      if (pivot_abs == 0.0) {   /* remainder of column == 0 */

	/*
	if (qh IStracing >= 4) {
	  fprintf (qh ferr, "qh_gausselim: 0 pivot at column %d. (%2.2g < %2.2g)\n", k, pivot_abs, qh DISTround);
	  qh_printmatrix (qh ferr, "Matrix:", rows, numrow, numcol);
	}
	*/

	zzinc_(Zgauss0);
	goto LABELnextcol;
      }
    }
    pivotrow= rows[k] + k;
    pivot= *pivotrow++;  /* signed value of pivot, and remainder of row */
    for(i= k+1; i < numrow; i++) {
      ai= rows[i] + k;
      ak= pivotrow;
      n= (*ai++)/pivot;   /* divzero() not needed since |pivot| >= |*ai| */
      for(j= numcol - (k+1); j--; )
	*ai++ -= n * *ak++;
    }
  LABELnextcol:
    ;
  }
  wmin_(Wmindenom, pivot_abs);  /* last pivot element */
  if (qh IStracing >= 5)
    qh_printmatrix (qh ferr, "qh_gausselem: result", rows, numrow, numcol);
} /* gausselim */


/*----------------------------------------------
-getangle- returns the dot product of two, qh hull_dim vectors
  may be > 1.0 or < -1.0
*/
realT qh_getangle(pointT *vect1, pointT *vect2) {
  realT angle= 0;
  int k;

  for(k= qh hull_dim; k--; )
    angle += *vect1++ * *vect2++;
  trace4((qh ferr, "qh_getangle: %2.2g\n", angle));
  return(angle);
} /* getangle */


/*----------------------------------------------
-getcenter-  gets arithmetic center of a set of vertices as a new point
  assumes normal_size is in short memory
*/
pointT *qh_getcenter(setT *vertices) {
  int k;
  pointT *center, *coord;
  vertexT *vertex, **vertexp;
  int count= qh_setsize(vertices);

  if (count < 2) qhull_fatal(23);

  center= (pointT *)qh_memalloc(qh normal_size);
  for (k=0; k < qh hull_dim; k++) {
    coord= center+k;
    *coord= 0.0;
    FOREACHvertex_(vertices)
      *coord += vertex->point[k];
    *coord /= count;
  }
  return(center);
} /* getcenter */


/*----------------------------------------------
-getcentrum- returns the centrum for a facet as a new point
  assumes normal_size is in short memory
*/
pointT *qh_getcentrum(facetT *facet) {
  realT dist;
  pointT *centrum, *point;

  point= qh_getcenter(facet->vertices);
  zinc_(Zcentrumtests);
  qh_distplane (point, facet, &dist);
  centrum= qh_projectpoint(point, facet, dist);
  qh_memfree(point, qh normal_size);
  trace4((qh ferr, "qh_getcentrum: for f%d, %d vertices dist= %2.2g\n",
	  facet->id, qh_setsize(facet->vertices), dist));
  return centrum;
} /* getcentrum */


/*-------------------------------------------------
-gram_schmidt- implements Gram-Schmidt orthogonalization by rows
   overwrites rows[dim][dim]
returns:
   false if gets a zero norm
notes:
   see Golub & van Loan Algorithm 6.2-2
   overflow due to small divisors not handled
*/
boolT qh_gram_schmidt(int dim, realT **row) {
  realT *rowi, *rowj, norm;
  int i, j, k;
  
  for(i=0; i < dim; i++) {
    rowi= row[i];
    for (norm= 0.0, k= dim; k--; rowi++)
      norm += *rowi * *rowi;
    norm= sqrt(norm);
    wmin_(Wmindenom, norm);
    if (norm == 0.0)  /* either 0 or overflow due to sqrt */
      return False;
    for(k= dim; k--; )
      *(--rowi) /= norm;  
    for(j= i+1; j < dim; j++) {
      rowj= row[j];
      for(norm= 0.0, k=dim; k--; )
	norm += *rowi++ * *rowj++;
      for(k=dim; k--; )
	*(--rowj) -= *(--rowi) * norm;
    }
  }
  return True;
} /* gram_schmidt */



	
/*--------------------------------------------------
-inthresholds- return True if normal within qh lower_/upper_threshold
returns:
  angle cos to a threshold border (may be NULL, invalid if qh SPLITthresholds)
*/
boolT qh_inthresholds (coordT *normal, realT *angle) {
  boolT within= True;
  int k;

  if (angle)
    *angle= 0.0;
  for(k= 0; k < qh hull_dim; k++) {
    if (qh lower_threshold[k] > -REALmax/2) {
      if (normal[k] < qh lower_threshold[k])
        within= False;
      if (angle)
        *angle += normal[k] * qh lower_threshold[k];
    }
    if (qh upper_threshold[k] < REALmax/2) {
      if (normal[k] > qh upper_threshold[k])
        within= False;
      if (angle)
        *angle += normal[k] * qh upper_threshold[k];
    }
  }
  return within;
} /* inthresholds */
    

/*--------------------------------------------------
-maxabsval -- return pointer to maximum absolute value of a dim vector
   returns NULL if dim==0
*/
realT *qh_maxabsval (realT *normal, int dim) {
  realT maxval= -REALmax;
  realT *maxp= NULL, *colp, absval;
  int k;

  for (k= dim, colp= normal; k--; colp++) {
    absval= fabs_(*colp);
    if (absval > maxval) {
      maxval= absval;
      maxp= colp;
    }
  }
  return maxp;
} /* maxabsval */


/*-------------------------------------------------
-maxmin- collects the maximum and minimum points of input into a set
  determines maximum roundoff errors
returns:
  returns a temporary set, without qh GOODpoint
  points are not unique
*/
setT *qh_maxmin(pointT *points, int numpoints, int dimension) {
  int k;
  realT maxsum= 0.0, maxcoord, temp, maxdistsum;
  realT maxneg= REALmax, maxpos= -REALmax;
  pointT *minimum, *maximum, *point, *pointtemp;
  setT *set;

  set= qh_settemp(2*dimension);
  for(k= 0; k < dimension; k++) {
    if (points == qh GOODpointp)
      minimum= maximum= points + qh hull_dim;
    else
      minimum= maximum= points;
    FORALLpoint_(points, numpoints) {
      if (point == qh GOODpointp)
	continue;
      if (maximum[k] < point[k])
	maximum= point;
      else if (minimum[k] > point[k])
	minimum= point;
    }
    maxcoord= fmax_(maximum[k], -minimum[k]);
    if (qh GOODpointp) {
      temp= fmax_(qh GOODpointp[k], -qh GOODpointp[k]);
      maximize_(maxcoord, temp);
    }
    maximize_(qh maxmaxcoord, maxcoord);
    maxsum += maxcoord;
    maximize_(maxpos, maximum[k]);
    minimize_(maxneg, minimum[k]);
    qh_setappend (&set, maximum);
    qh_setappend (&set, minimum);
    /* calculation of qh NEARzero is based on error formula 4.4-13 of
       Golub & van Loan, authors say n^3 can be ignored and 10 be used in
       place of rho */
    qh NEARzero[k]= 80 * maxsum * REALepsilon;
  }
  /* calculate roundoff error according to
     Lemma 3.2-1 of Golub and van Loan "Matrix Computation"
     use sqrt(dim) since one vector is normalized */
  maxdistsum= sqrt (qh hull_dim) * qh maxmaxcoord;
  if (!qh SETroundoff) {
    qh DISTround= REALepsilon * (qh hull_dim * maxdistsum * 1.01
			   	       + qh maxmaxcoord);  /* for offset */
    if (qh RANDOMdist)
      qh DISTround += qh RANDOMfactor * qh maxmaxcoord;
  }
  qh MINdenom= qh MINdenom_1 * qh maxmaxcoord;
  qh MINdenom_1_2= sqrt (qh MINdenom_1 * qh hull_dim) ;  /* if will be normalized */
  qh MINdenom_2= qh MINdenom_1_2 * qh maxmaxcoord;
  if (qh premerge_cos < REALmax/2)      /* for inner product */
    qh premerge_cos -= 1.01 * qh hull_dim * REALepsilon; 
  if (qh postmerge_cos < REALmax/2)