coord= qh first_point;    infinity= qh first_point + qh hull_dim * qh num_points;    for (k=qh hull_dim-1; k--; )      infinity[k]= 0.0;    for (i=qh num_points; i--; ) {      paraboloid= 0.0;      for (k=qh hull_dim-1; k--; ) {        paraboloid += *coord * *coord;	infinity[k] += *coord;        coord++;      }      *(coord++)= paraboloid;      maximize_(maxboloid, paraboloid);    }    for (k=qh hull_dim-1; k--; )      *(coord++) /= qh num_points;    *(coord++)= maxboloid * 1.1;    qh num_points++;    trace0((qh ferr, "qh_projectinput: projected points to paraboloid for Delaunay\n"));  }} /* projectinput */  /*--------------------------------------------------projectpoints- project along one or more dimensions  delete dimension k if project[k] == -1  add dimension k if project[k] == 1   n is size of project  points, numpoints, dim is old points  newpoints, newdim is buffer for new points (already allocated)    newpoints may be points if only adding dimension at end*/void qh_projectpoints (signed char *project, int n, realT *points,         int numpoints, int dim, realT *newpoints, int newdim) {  int testdim= dim, oldk=0, newk=0, i,j=0,k;  realT *newp, *oldp;    for (k= 0; k<n; k++)    testdim += project[k];  if (testdim != newdim) {    fprintf (qh ferr, "qhull internal error (qh_projectpoints): newdim %d should be %d after projection\n",      newdim, testdim);    qh_errexit (qh_ERRqhull, NULL, NULL);  }  for (j= 0; j<n; j++) {    if (project[j] == -1)      oldk++;    else {      newp= newpoints+newk++;      if (project[j] == +1) {	if (oldk >= dim)	  continue;	oldp= points+oldk;      }else 	oldp= points+oldk++;      for (i=numpoints; i--; ) {        *newp= *oldp;        newp += newdim;        oldp += dim;      }    }    if (oldk >= dim)      break;  }  trace1((qh ferr, "qh_projectpoints: projected %d points from dim %d to dim %d\n",     numpoints, dim, newdim));} /* projectpoints */        /*--------------------------------------------------rand & srand- generate pseudo-random number between 1 and 2^31 -2  from Park & Miller's minimimal standard random number generator  Communications of the ACM, 31:1192-1201, 1988.notes:  does not use 0 or 2^31 -1  this is silently enforced by qh_srand()  copied to rbox.c  can make 'Rn' much faster by moving qh_rand to qh_distplane*/int qh_rand( void) {#define qh_rand_a 16807#define qh_rand_m 2147483647#define qh_rand_q 127773  /* m div a */#define qh_rand_r 2836    /* m mod a */  int lo, hi, test;  int seed = qh rand_seed;  hi = seed / qh_rand_q;  /* seed div q */  lo = seed % qh_rand_q;  /* seed mod q */  test = qh_rand_a * lo - qh_rand_r * hi;  if (test > 0)    seed= test;  else    seed= test + qh_rand_m;  qh rand_seed= seed;  return seed;} /* rand */void qh_srand( int seed) {  if (seed < 1)    qh rand_seed= 1;  else if (seed >= qh_rand_m)    qh rand_seed= qh_rand_m - 1;  else    qh rand_seed= seed;} /* qh_srand *//*--------------------------------------------------randomfactor- return a random factor within qh RANDOMmax of 1.0  RANDOMa/b definedin global.c*/realT qh_randomfactor (void) {  realT randr;  randr= qh_RANDOMint;  return randr * qh RANDOMa + qh RANDOMb;} /* randomfactor *//*--------------------------------------------------randommatrix- generate a random dimXdim matrix in range (-1,1)  assumes buffer is dim+1Xdimreturns:  returns row vector for buffer  plus row[dim] for scratch*/void qh_randommatrix (realT *buffer, int dim, realT **row) {  int i, k;  realT **rowi, *coord, realr;  coord= buffer;  rowi= row;  for (i=0; i<dim; i++) {    *(rowi++)= coord;    for (k=0; k<dim; k++) {      realr= qh_RANDOMint;      *(coord++)= 2.0 * realr/(qh_RANDOMmax+1) - 1.0;    }  }  *rowi= coord;} /* randommatrix */        /*--------------------------------------------------rotateinput- rotate input using row matrix  input points given by qh first_point, num_points, hull_dim  if qh POINTSmalloc, overwrites input points, else mallocs a new array  assumes rows[dim] is a scratch bufferreturns:  sets qh POINTSmalloc*/void qh_rotateinput (realT **rows) {  int size;  pointT *newpoints;  if (!qh POINTSmalloc) {    size= qh num_points*qh hull_dim*sizeof(pointT);    if (!(newpoints=(coordT*)malloc(size))) {      fprintf(qh ferr, "qhull error: insufficient memory to rotate %d points\n",          qh num_points);      qh_errexit(qh_ERRmem, NULL, NULL);    }    memcpy ((char *)newpoints, (char *)qh first_point, size);    qh first_point= newpoints;    qh POINTSmalloc= True;  }  qh_rotatepoints (qh first_point, qh num_points, qh hull_dim, rows);}  /* rotateinput *//*--------------------------------------------------rotatepoints- rotate numpoints points by a row matrix  assumes rows[dim] is a scratch buffer*/void qh_rotatepoints (realT *points, int numpoints, int dim, realT **row) {  realT *point, *rowi, *coord= NULL, sum, *newval;  int i,j,k;  for (point= points, j= numpoints; j--; point += dim) {    newval= row[dim];    for (i= 0; i<dim; i++) {      rowi= row[i];      coord= point;      for (sum= 0.0, k= dim; k--; )        sum += *rowi++ * *coord++;      *(newval++)= sum;    }    for (k= dim; k--; )      *(--coord)= *(--newval);  }} /* rotatepoints */    /*--------------------------------------------------scaleinput- scale input points using qh low_bound/high_bound  input points given by qh first_point, num_points, hull_dim  if qh POINTSmalloc, overwrites input points, else mallocs a new arrayreturns:  scales points to low[k], high[k]  sets qh POINTSmalloc*/void qh_scaleinput (void) {  int size;  pointT *newpoints;  if (!qh POINTSmalloc) {    size= qh num_points*qh hull_dim*sizeof(pointT);    if (!(newpoints=(coordT*)malloc(size))) {      fprintf(qh ferr, "qhull error: insufficient memory to scale %d points\n",          qh num_points);      qh_errexit(qh_ERRmem, NULL, NULL);    }    memcpy ((char *)newpoints, (char *)qh first_point, size);    qh first_point= newpoints;    qh POINTSmalloc= True;  }  qh_scalepoints (qh first_point, qh num_points, qh hull_dim,       qh lower_bound, qh upper_bound);}  /* scaleinput */  /*--------------------------------------------------scalepoints- scale points to new lowbound and highbound  retains old bound when newlow= -REALmax or newhigh= +REALmax  overwrites old points*/void qh_scalepoints (pointT *points, int numpoints, int dim,	realT *newlows, realT *newhighs) {  int i,k;  realT shift, scale, *coord, low, high, newlow, newhigh, mincoord, maxcoord;  boolT nearzero= False;       for (k= 0; k<dim; k++) {    newhigh= newhighs[k];    newlow= newlows[k];    if (newhigh > REALmax/2 && newlow < -REALmax/2)      continue;    low= REALmax;    high= -REALmax;    for (i= numpoints, coord= points+k; i--; coord += dim) {      minimize_(low, *coord);      maximize_(high, *coord);    }    if (newhigh > REALmax/2)      newhigh= high;    if (newlow < -REALmax/2)      newlow= low;    scale= qh_divzero (newhigh - newlow, high - low,                  qh MINdenom_1, &nearzero);    if (nearzero) {      fprintf (qh ferr, "qhull input error: %d'th dimension's new bounds [%2.2g, %2.2g] too wide for\nexisting bounds [%2.2g, %2.2g]\n",              k, newlow, newhigh, low, high);      qh_errexit (qh_ERRinput, NULL, NULL);    }    shift= (newlow * high - low * newhigh)/(high-low);    coord= points+k;    for (i= numpoints; i--; coord += dim)      *coord= *coord * scale + shift;    coord= points+k;    if (newlow < newhigh) {      mincoord= newlow;      maxcoord= newhigh;    }else {      mincoord= newhigh;      maxcoord= newlow;    }    for (i= numpoints; i--; coord += dim) {      minimize_(*coord, maxcoord);  /* because of roundoff error */      maximize_(*coord, mincoord);    }    trace0((qh ferr, "qh_scalepoints: scaled %d'th coordinate [%2.2g, %2.2g] to [%.2g, %.2g] for %d points by %2.2g and shifted %2.2g\n",      k, low, high, newlow, newhigh, numpoints, scale, shift));  }} /* scalepoints */           /*--------------------------------------------sethalfspace- set coords to dual of halfspace relative to feasible point  the halfspace is its normal coefficients and offset.returns:  false if feasible point is outside of hull (error message reported)  next value for coords*/boolT qh_sethalfspace (int dim, coordT *coords, coordT **nextp,          coordT *normal, coordT *offset, coordT *feasible) {  coordT *normp= normal, *feasiblep= feasible, *coordp= coords;  realT dist;  realT r; /*bug fix*/  int k;  boolT zerodiv;  dist= *offset;  for (k= dim; k--; )    dist += *(normp++) * *(feasiblep++);  if (dist > 0)    goto LABELerroroutside;  normp= normal;  if (dist < -qh MINdenom) {    for (k= dim; k--; )      *(coordp++)= *(normp++) / -dist;  }else {    for (k= dim; k--; ) {      *(coordp++)= qh_divzero (*(normp++), -dist, qh MINdenom_1, &zerodiv);      if (zerodiv)         goto LABELerroroutside;    }  }  *nextp= coordp;  if (qh IStracing >= 4) {    fprintf (qh ferr, "qh_sethalfspace: halfspace at offset %6.2g to point: ", *offset);    for (k= dim, coordp= coords; k--; )      fprintf (qh ferr, " %6.2g", r=*coordp++);    fprintf (qh ferr, "\n");  }  return True;LABELerroroutside:  feasiblep= feasible;  normp= normal;  fprintf(qh ferr, "qhull input error: feasible point is not clearly inside halfspace\nfeasible point: ");  for (k= dim; k--; )    fprintf (qh ferr, qh_REAL_1, r=*(feasiblep++));  fprintf (qh ferr, "\n     halfspace: ");   for (k= dim; k--; )    fprintf (qh ferr, qh_REAL_1, r=*(normp++));  fprintf (qh ferr, "\n     at offset: ");  fprintf (qh ferr, qh_REAL_1, *offset);  fprintf (qh ferr, " and distance: ");  fprintf (qh ferr, qh_REAL_1, dist);  fprintf (qh ferr, "\n");  return False;} /* sethalfspace *//*--------------------------------------------------sethalfspace_all- generate dual for halfspace intersection with feasible point     each halfspace is normal coefficients followed by offset      the origin is inside the halfspace if the offset is negativereturns:  unused/untested code: please email barber@tiac.net if this works ok for you  malloc'd array of count X dim-1 points  call before qh_init_B or qh_initqhull_globals notes:  If using option 'Fp', also set qh feasible_point. It is a malloc'd array that  is freed by qh_freebuffers.*/coordT *qh_sethalfspace_all (int dim, int count, coordT *halfspaces, pointT *feasible) {  int i, newdim;  pointT *newpoints;  coordT *coordp, *normalp, *offsetp;  trace0((qh ferr, "qh_sethalfspace_all: compute dual for halfspace intersection\n"));  newdim= dim - 1;  if (!(newpoints=(coordT*)malloc(count*newdim*sizeof(coordT)))){    fprintf(qh ferr, "qhull error: insufficient memory to compute dual of %d halfspaces\n",          count);    qh_errexit(qh_ERRmem, NULL, NULL);  }  coordp= newpoints;  normalp= halfspaces;  for (i= 0; i < count; i++) {    offsetp= normalp + newdim;    if (!qh_sethalfspace (newdim, coordp, &coordp, normalp, offsetp, feasible)) {      fprintf (qh ferr, "The halfspace was at index %d\n", i);      qh_errexit (qh_ERRinput, NULL, NULL);    }    normalp= offsetp + 1;  }  return newpoints;} /* sethalfspace_all */  /*--------------------------------------------voronoi_center- return Voronoi center for a set of points  dim is the orginal dimension of the pointsnotes:  if non-simplicial, returns center for max simplex of points  from Bowyer & Woodwark, A Programmer's Geometry, 1983, p. 65*/pointT *qh_voronoi_center (int dim, setT *points) {  pointT *point, **pointp, *point0;  pointT *center= (pointT*)qh_memalloc (qh center_size);  setT *simplex;  int i, j, k, num, size= qh_setsize(points);  coordT *gmcoord;  realT *diffp, sum2, *sum2row, *sum2p, det, factor;  boolT nearzero, infinite;  if (size == dim+1)    simplex= points;  else if (size < dim+1) {    fprintf (qh ferr, "qhull internal error (qh_voronoi_center):\n  need at least %d points to construct a Voronoi center\n",	     dim+1);    qh_errexit (qh_ERRqhull, NULL, NULL);  }else {    simplex= qh_settemp (dim+1);    qh_maxsimplex (dim, points, NULL, 0, &simplex);  }  num= qh_setsize (simplex);  point0= SETfirst_(simplex);  gmcoord= qh gm_matrix;  for (k=0; k<dim; k++) {    qh gm_row[k]= gmcoord;    FOREACHpoint_(simplex) {      if (point != point0)        *(gmcoord++)= point[k] - point0[k];    }  }  sum2row= gmcoord;  for (i=0; i<dim; i++) {    sum2= 0.0;    for (k= 0; k<dim; k++) {      diffp= qh gm_row[k] + i;      sum2 += *diffp * *diffp;    }    *(gmcoord++)= sum2;  }  det= qh_determinant (qh gm_row, dim, &nearzero);  factor= qh_divzero (0.5, det, qh MINdenom, &infinite);  if (infinite) {    for (k=dim; k--; )      center[k]= qh_INFINITE;    if (qh IStracing)      qh_printpoints (qh ferr, "qh_voronoi_center: at infinity for ", simplex);  }else {    for (i=0; i<dim; i++) {      gmcoord= qh gm_matrix;      sum2p= sum2row;      for (k=0; k<dim; k++) {	qh gm_row[k]= gmcoord;	if (k == i) {	  for (j= dim; j--; )	    *(gmcoord++)= *sum2p++;	}else {	  FOREACHpoint_(simplex) {	    if (point != point0)	      *(gmcoord++)= point[k] - point0[k];	  }	}      }      center[i]= qh_determinant (qh gm_row, dim, &nearzero)*factor + point0[i];    }#ifndef qh_NOtrace    if (qh IStracing >= 3) {      fprintf (qh ferr, "qh_voronoi_center: det %2.2g factor %2.2g ", det, factor);      qh_printmatrix (qh ferr, "center:", &center, 1, dim);      if (qh IStracing >= 5) {	qh_printpoints (qh ferr, "points", simplex);	FOREACHpoint_(simplex)	  fprintf (qh ferr, "p%d dist %.2g, ", qh_pointid (point),		   qh_pointdist (point, center, dim));	fprintf (qh ferr, "\n");      }    }#endif  }  if (simplex != points)    qh_settempfree (&simplex);  return center;} /* voronoi_center */