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 errorsreturns:  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;  if (REALmin < REALepsilon && REALmin < REALmax && REALmin > -REALmax  && REALmax > 0.0 && -REALmax < 0.0)    ; /* all ok */  else {    fprintf (qh ferr, "qhull error: floating point constants in user.h are wrong\n\REALepsilon %g REALmin %g REALmax %g -REALmax %g\n",	     REALepsilon, REALmin, REALmax, -REALmax);    qh_errexit (qh_ERRinput, NULL, NULL);  }  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_option ("Error-roundoff", NULL, &qh DISTround);  }  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;                                              /* for inner product */  qh ANGLEround= 1.01 * qh hull_dim * REALepsilon;  if (qh RANDOMdist)    qh ANGLEround += qh RANDOMfactor;  if (qh premerge_cos < REALmax/2) {    qh premerge_cos -= qh ANGLEround;    if (qh RANDOMdist)       qh_option ("Angle-premerge-with-random", NULL, &qh premerge_cos);  }  if (qh postmerge_cos < REALmax/2) {    qh postmerge_cos -= qh ANGLEround;    if (qh RANDOMdist)      qh_option ("Angle-postmerge-with-random", NULL, &qh postmerge_cos);  }  qh premerge_centrum += 2 * qh DISTround;    /*2 for centrum and distplane()*/  qh postmerge_centrum += 2 * qh DISTround;  if (qh RANDOMdist && (qh MERGEexact || qh PREmerge))    qh_option ("Centrum-premerge-with-random", NULL, &qh premerge_centrum);  if (qh RANDOMdist && qh POSTmerge)    qh_option ("Centrum-postmerge-with-random", NULL, &qh postmerge_centrum);  qh_option ("_max-coord", NULL, &maxpos);  qh_option ("_min-coord", NULL, &maxneg);  { /* compute ONEmerge, max vertex offset for merging simplicial facets */    realT maxangle= 1.0, maxrho;        minimize_(maxangle, qh premerge_cos);    minimize_(maxangle, qh postmerge_cos);    /* max diameter * sin theta + DISTround for vertex to its hyperplane */    qh ONEmerge= sqrt (qh hull_dim) * (maxpos - maxneg) *      sqrt (1.0 - maxangle * maxangle) + qh DISTround;      maxrho= qh hull_dim * qh premerge_centrum + qh DISTround;    maximize_(qh ONEmerge, maxrho);    maxrho= qh hull_dim * qh postmerge_centrum + qh DISTround;    maximize_(qh ONEmerge, maxrho);    if (qh MERGING)      qh_option ("_one-merge", NULL, &qh ONEmerge);  }  qh NEARinside= qh ONEmerge * qh_RATIOnearinside; /* only used if qh KEEPnearinside */  if (qh KEEPnearinside)    qh_option ("_near-inside", NULL, &qh NEARinside);  qh MINnorm= sqrt( REALepsilon) * qh maxmaxcoord; /* FIXUP: what is correct?*/  if (qh hull_dim <= 4) /* used in qh_sethyperplane_det */    qh_option ("_min-norm", NULL, &qh MINnorm);  if (qh MINvisible > REALmax/2) {    if (!qh MERGING)      qh MINvisible= qh DISTround;    else if (qh hull_dim <= 3)      qh MINvisible= qh premerge_centrum;    else      qh MINvisible= qh_COPLANARratio * qh premerge_centrum;    if (qh APPROXhull && qh MINvisible > qh MINoutside)      qh MINvisible= qh MINoutside;    qh_option ("Visible-distance", NULL, &qh MINvisible);  }  if (qh MAXcoplanar > REALmax/2) {    qh MAXcoplanar= qh MINvisible;    qh_option ("U-coplanar-distance", NULL, &qh MAXcoplanar);  }  if (!qh APPROXhull) {             /* user may specify qh MINoutside */    qh MINoutside= 2 * qh MINvisible;    if (qh premerge_cos < REALmax/2)       maximize_(qh MINoutside, (1- qh premerge_cos) * qh maxmaxcoord);    qh_option ("Width-outside", NULL, &qh MINoutside);  }  qh WIDEfacet= qh MINoutside;  maximize_(qh WIDEfacet, qh_WIDEcoplanar * qh MAXcoplanar);   maximize_(qh WIDEfacet, qh_WIDEcoplanar * qh MINvisible);   qh_option ("_wide-facet", NULL, &qh WIDEfacet);  if (qh MINvisible > qh MINoutside + 3 * REALepsilon   && !qh BESToutside && !qh FORCEoutput)    fprintf (qh ferr, "qhull input warning: minimum visibility V%.2g is greater than \nminimum outside W%.2g.  Flipped facets are likely.\n",	     qh MINvisible, qh MINoutside);  qh max_vertex= qh DISTround;  qh min_vertex= -qh DISTround;  if (qh IStracing >=1)    qh_printpoints (qh ferr, "qh_maxmin: found the max and min points (by dim):", set);  /* numeric constants reported in printsummary */  return(set);} /* maxmin *//*--------------------------------------------------maxsimplex- determines maximum simplex for a set of points   assumes at least pointsneeded points in points  skips qh GOODpointp (assumes that it isn't in maxpoints)  starts from points already in simplexreturns:  temporary set of dim+1 pointsnotes:  maximizes determinate for x,y,z,w, etc.  uses maxpoints as long as determinate is clearly non-zero*/void qh_maxsimplex (int dim, setT *maxpoints, pointT *points, int numpoints, setT **simplex) {  pointT *point, **pointp, *pointtemp, *maxpoint, *minx=NULL, *maxx=NULL;  boolT nearzero, maxnearzero= False;  int k, sizinit;  realT maxdet= -REALmax, det, mincoord= REALmax, maxcoord= -REALmax;  sizinit= qh_setsize (*simplex);  if (sizinit < 2) {    if (qh_setsize (maxpoints) >= 2) {      FOREACHpoint_(maxpoints) {	        if (maxcoord < point[0]) {          maxcoord= point[0];          maxx= point;        }	if (mincoord > point[0]) {          mincoord= point[0];          minx= point;        }      }    }else {      FORALLpoint_(points, numpoints) {	if (point == qh GOODpointp)	  continue;        if (maxcoord < point[0]) {	  maxcoord= point[0];          maxx= point;        }	if (mincoord > point[0]) {          mincoord= point[0];          minx= point;	}      }    }    qh_setunique (simplex, minx);    if (qh_setsize (*simplex) < 2)      qh_setunique (simplex, maxx);    sizinit= qh_setsize (*simplex);    if (sizinit < 2) {      if (zzval_(Zsetplane) > qh hull_dim+1) {	fprintf (qh ferr, "qhull precision error (qh_maxsimplex for voronoi_center):\n%d points with the same x coordinate.\n",		 qh_setsize(maxpoints)+numpoints);	qh_errexit (qh_ERRprec, NULL, NULL);      }else {	fprintf (qh ferr, "qhull input error: input is less than %d-dimensional since it has the same x coordinate\n", qh hull_dim);	qh_errexit (qh_ERRinput, NULL, NULL);      }    }  }  for(k= sizinit; k < dim+1; k++) {    maxpoint= NULL;    maxdet= -REALmax;    FOREACHpoint_(maxpoints) {      if (!qh_setin (*simplex, point)) {        det= qh_detsimplex(point, *simplex, k, &nearzero);        if ((det= fabs_(det)) > maxdet) {	  maxdet= det;          maxpoint= point;	  maxnearzero= nearzero;        }      }    }    if (!maxpoint || maxnearzero) {      zinc_(Zsearchpoints);      if (!maxpoint) {        trace0((qh ferr, "qh_maxsimplex: searching all points for %d-th initial vertex.\n", k));      }else {        trace0((qh ferr, "qh_maxsimplex: searching all points for %d-th initial vertex, better than p%d det %2.2g\n",		k+1, qh_pointid(maxpoint), maxdet));      }      FORALLpoint_(points, numpoints) {	if (point == qh GOODpointp)	  continue;        if (!qh_setin (*simplex, point)) {          det= qh_detsimplex(point, *simplex, k, &nearzero);          if ((det= fabs_(det)) > maxdet) {	    maxdet= det;            maxpoint= point;	    maxnearzero= nearzero;	  }        }      }    } /* !maxpoint */    if (!maxpoint) {      fprintf (qh ferr, "qhull internal error (qh_maxsimplex): not enough points available\n");      qh_errexit (qh_ERRqhull, NULL, NULL);    }    qh_setappend(simplex, maxpoint);    trace1((qh ferr, "qh_maxsimplex: selected point p%d for %d`th initial vertex, det=%2.2g\n",	    qh_pointid(maxpoint), k, maxdet));  } /* k */ } /* maxsimplex *//*---------------------------------------------------minabsval -- return min absolute value of a dim vector*/realT qh_minabsval (realT *normal, int dim) {  realT minval= 0;  realT maxval= 0;  realT *colp;  int k;  for (k= dim, colp= normal; k--; colp++) {    maximize_(maxval, *colp);    minimize_(minval, *colp);  }  return fmax_(maxval, -minval);} /* minabsval *//*---------------------------------------------------mindiff -- return index of min abs. difference of two vectors*/int qh_mindiff (realT *vecA, realT *vecB, int dim) {  realT mindiff= REALmax, diff;  realT *vecAp= vecA, *vecBp= vecB;  int k, mink;  for (k= 0; k<dim; k++) {    diff= *vecAp++ - *vecBp++;    diff= fabs_(diff);    if (diff < mindiff) {      mindiff= diff;      mink= k;    }  }  return k;} /* mindiff *//*--------------------------------------------orientoutside- make facet outside oriented via qh interior_point  returns True if reversed orientation.*/boolT qh_orientoutside (facetT *facet) {  int k;  realT dist;  qh_distplane (qh interior_point, facet, &dist);  if (dist > 0) {    for (k= qh hull_dim; k--; )      facet->normal[k]= -facet->normal[k];    facet->offset= -facet->offset;    return True;  }  return False;} /* orientoutside *//*--------------------------------------------pointdist- distance between points  returns distance squared if 'dim' is negative*/coordT qh_pointdist(pointT *point1, pointT *point2, int dim) {  coordT dist, diff;  int k;    dist= 0.0;  for (k= (dim > 0 ? dim : -dim); k--; ) {    diff= *point1++ - *point2++;    dist += diff * diff;  }  if (dim > 0)    return(sqrt(dist));  return dist;} /* pointdist *//*--------------------------------------------------printmatrix- print matrix given by row vectors  print a vector by (fp, "", &vect, 1, len)*/void qh_printmatrix (FILE *fp, char *string, realT **rows, int numrow, int numcol) {  realT *rowp;  realT r; /*bug fix*/  int i,k;  fprintf (fp, "%s\n", string);  for (i= 0; i<numrow; i++) {    rowp= rows[i];    for (k= 0; k<numcol; k++)      fprintf (fp, "%6.3g ", r=*rowp++);    fprintf (fp, "\n");  }} /* printmatrix */  /*--------------------------------------------------printpoints- print pointids for a set of points starting at index   prints string and 'p' if defined*/void qh_printpoints (FILE *fp, char *string, setT *points) {  pointT *point, **pointp;  if (string) {    fprintf (fp, "%s", string);    FOREACHpoint_(points)       fprintf (fp, " p%d", qh_pointid(point));    fprintf (fp, "\n");  }else {    FOREACHpoint_(points)       fprintf (fp, " %d", qh_pointid(point));    fprintf (fp, "\n");  }} /* printpoints */  /*--------------------------------------------------projectinput- project input points using qh DELAUNAY and qh lower_bound/upper_bound  input points in qh first_point, num_points, input_dim     if POINTSmalloc, will free old point array  if low[k]=high[k]= 0, removes dimension k      checks that hull_dim agrees with input_dim, PROJECTinput, and DELAUNAY  if DELAUNAY     projects points to paraboloidreturns:  new point array in first_point of qh hull_dim coordinates  sets POINTSmalloc  lowbound/highbound is also projected*/void qh_projectinput (void) {  int k,i;  int newdim= qh input_dim, newnum= qh num_points;  signed char *project;  int size= (qh input_dim+1)*sizeof(*project);  pointT *newpoints, *coord, *infinity;  realT paraboloid, maxboloid= 0;    project= (signed char*)qh_memalloc (size);  memset ((char*)project, 0, size);  for (k= 0; k<qh input_dim; k++) {   /* skip Delaunay bound */    if (qh lower_bound[k] == 0 && qh upper_bound[k] == 0) {      project[k]= -1;      newdim--;    }  }  if (qh DELAUNAY) {    project[k]= 1;    newdim++;    if (qh ATinfinity)      newnum++;  }  if (newdim != qh hull_dim) {    fprintf(qh ferr, "qhull internal error (qh_projectinput): dimension after projection %d != hull_dim %d\n", newdim, qh hull_dim);    qh_errexit(qh_ERRqhull, NULL, NULL);  }  if (!(newpoints=(coordT*)malloc(newnum*newdim*sizeof(coordT)))){    fprintf(qh ferr, "qhull error: insufficient memory to project %d points\n",           qh num_points);    qh_errexit(qh_ERRmem, NULL, NULL);  }  qh_projectpoints (project, qh input_dim+1, qh first_point,                    qh num_points, qh input_dim, newpoints, newdim);  trace1((qh ferr, "qh_projectinput: updating lower and upper_bound\n"));  qh_projectpoints (project, qh input_dim+1, qh lower_bound,                    1, qh input_dim+1, qh lower_bound, newdim+1);  qh_projectpoints (project, qh input_dim+1, qh upper_bound,                    1, qh input_dim+1, qh upper_bound, newdim+1);  qh_memfree(project, ((qh input_dim+1)*sizeof(*project)));  if (qh POINTSmalloc)    free (qh first_point);  qh first_point= newpoints;  qh POINTSmalloc= True;  if (qh DELAUNAY && qh ATinfinity) {