qhull_geom2.cxx

hl2 source code. Do not use it illegal.

*/
realT qh_randomfactor (void) {
  realT randr;

  randr= qh_RANDOMint;
  return randr * qh RANDOMa + qh RANDOMb;
} /* randomfactor */

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="randommatrix">-</a>
  
  qh_randommatrix( buffer, dim, rows )
    generate a random dim X dim matrix in range [-1,1]
    assumes buffer is [dim+1, dim]

  returns:
    sets buffer to random numbers
    sets rows to rows of buffer
      sets row[dim] as scratch row
*/
void qh_randommatrix (realT *buffer, int dim, realT **rows) {
  int i, k;
  realT **rowi, *coord, realr;

  coord= buffer;
  rowi= rows;
  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 */

        
/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="rotateinput">-</a>
  
  qh_rotateinput( rows )
    rotate input using row matrix
    input points given by qh first_point, num_points, hull_dim
    assumes rows[dim] is a scratch buffer
    if qh POINTSmalloc, overwrites input points, else mallocs a new array

  returns:
    rotated input
    sets qh POINTSmalloc

  design:
    see qh_rotatepoints
*/
void qh_rotateinput (realT **rows) {

  if (!qh POINTSmalloc) {
    qh first_point= qh_copypoints (qh first_point, qh num_points, qh hull_dim);
    qh POINTSmalloc= True;
  }
  qh_rotatepoints (qh first_point, qh num_points, qh hull_dim, rows);
}  /* rotateinput */

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="rotatepoints">-</a>
  
  qh_rotatepoints( points, numpoints, dim, row )
    rotate numpoints points by a d-dim row matrix
    assumes rows[dim] is a scratch buffer

  returns:
    rotated points in place

  design:
    for each point
      for each coordinate
        use row[dim] to compute partial inner product
      for each coordinate
        rotate by partial inner product
*/
void qh_rotatepoints (realT *points, int numpoints, int dim, realT **row) {
  realT *point, *rowi, *coord= NULL, sum, *newval;
  int i,j,k;

  if (qh IStracing >= 1)
    qh_printmatrix (qh ferr, "qh_rotatepoints: rotate points by", row, dim, dim);
  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 */  
  

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="scaleinput">-</a>
  
  qh_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 array

  returns:
    scales coordinates of points to low_bound[k], high_bound[k]
    sets qh POINTSmalloc

  design:
    see qh_scalepoints
*/
void qh_scaleinput (void) {

  if (!qh POINTSmalloc) {
    qh first_point= qh_copypoints (qh first_point, qh num_points, qh hull_dim);
    qh POINTSmalloc= True;
  }
  qh_scalepoints (qh first_point, qh num_points, qh hull_dim,
       qh lower_bound, qh upper_bound);
}  /* scaleinput */
  
/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="scalelast">-</a>
  
  qh_scalelast( points, numpoints, dim, low, high, newhigh )
    scale last coordinate to [0,m] for Delaunay triangulations
    input points given by points, numpoints, dim

  returns:
    changes scale of last coordinate from [low, high] to [0, newhigh]
    overwrites last coordinate of each point
    saves low/high/newhigh in qh.last_low, etc. for qh_setdelaunay()

  notes:
    when called by qh_setdelaunay, low/high may not match actual data
    
  design:
    compute scale and shift factors
    apply to last coordinate of each point
*/
void qh_scalelast (coordT *points, int numpoints, int dim, coordT low,
		   coordT high, coordT newhigh) {
  realT scale, shift;
  coordT *coord;
  int i;
  boolT nearzero= False;

  trace4((qh ferr, "qh_scalelast: scale last coordinate from [%2.2g, %2.2g] to [0,%2.2g]\n",
    low, high, newhigh));
  qh last_low= low;
  qh last_high= high;
  qh last_newhigh= newhigh;
  scale= qh_divzero (newhigh, high - low,
                  qh MINdenom_1, &nearzero);
  if (nearzero) {
    ivp_message( "qhull input error: last coordinate's new bounds [0, %2.2g] too wide for\nexisting bounds [%2.2g, %2.2g] with width %2.2g\n",
              newhigh, low, high, high-low);
    qh_errexit (qh_ERRinput, NULL, NULL);
  }
  shift= - low * newhigh / (high-low);
  coord= points + dim - 1;
  for (i= numpoints; i--; coord += dim)
    *coord= *coord * scale + shift;
} /* scalelast */

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="scalepoints">-</a>
  
  qh_scalepoints( points, numpoints, dim, newlows, newhighs )
    scale points to new lowbound and highbound
    retains old bound when newlow= -REALmax or newhigh= +REALmax

  returns:
    scaled points
    overwrites old points

  design:
    for each coordinate
      compute current low and high bound
      compute scale and shift factors
      scale all points
      enforce new low and high bound for all 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;
    if (qh DELAUNAY && k == dim-1 && newhigh < newlow) {
      ivp_message( "qhull input error: 'Qb%d' or 'QB%d' inverts paraboloid since high bound %.2g < low bound %.2g\n",
	       k, k, newhigh, newlow);
      qh_errexit (qh_ERRinput, NULL, NULL);
    }
    scale= qh_divzero (newhigh - newlow, high - low,
                  qh MINdenom_1, &nearzero);
    if (nearzero) {
      ivp_message( "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 */    

       
/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="setdelaunay">-</a>
  
  qh_setdelaunay( dim, count, points )
    project count points to dim-d paraboloid for Delaunay triangulation
    
    dim is one more than the dimension of the input set
    assumes dim is at least 3 (i.e., at least a 2-d Delaunay triangulation)

    points is a dim*count realT array.  The first dim-1 coordinates
    are the coordinates of the first input point.  array[dim] is
    the first coordinate of the second input point.  array[2*dim] is
    the first coordinate of the third input point.

    if qh.last_low defined (i.e., 'Qbb' called qh_scalelast)
      calls qh_scalelast to scale the last coordinate the same as the other points

  returns:
    for each point
      sets point[dim-1] to sum of squares of coordinates
    scale points to 'Qbb' if needed
      
  notes:
    to project one point, use
      qh_setdelaunay (qh hull_dim, 1, point)
      
    Do not use options 'Qbk', 'QBk', or 'QbB' since they scale 
    the coordinates after the original projection.

*/
void qh_setdelaunay (int dim, int count, pointT *points) {
  int i, k;
  coordT *coordp, coord;
  realT paraboloid;

  trace0((qh ferr, "qh_setdelaunay: project %d points to paraboloid for Delaunay triangulation\n", count));
  coordp= points;
  for (i= 0; i < count; i++) {
    coord= *coordp++;
    paraboloid= coord*coord;
    for (k= dim-2; k--; ) {
      coord= *coordp++;
      paraboloid += coord*coord;
    }
    *coordp++ = paraboloid;
  }
  if (qh last_low < REALmax/2) 
    qh_scalelast (points, count, dim, qh last_low, qh last_high, qh last_newhigh);
} /* setdelaunay */

  
/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="sethalfspace">-</a>
  
  qh_sethalfspace( dim, coords, nextp, normal, offset, feasible )
    set point to dual of halfspace relative to feasible point
    halfspace is normal coefficients and offset.

  returns:
    false if feasible point is outside of hull (error message already reported)
    overwrites coordinates for point at dim coords
    nextp= next point (coords)

  design:
    compute distance from feasible point to halfspace
    divide each normal coefficient by -dist
*/
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) {
    ivp_message( "qh_sethalfspace: halfspace at offset %6.2g to point: ", *offset);
    for (k= dim, coordp= coords; k--; ) {
      r= *coordp++;
      ivp_message( " %6.2g", r);
    }
    ivp_message( "\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--; )
    ivp_message( qh_REAL_1, r=*(feasiblep++));
  ivp_message( "\n     halfspace: "); 
  for (k= dim; k--; )
    ivp_message( qh_REAL_1, r=*(normp++));
  ivp_message( "\n     at offset: ");
  ivp_message( qh_REAL_1, *offset);
  ivp_message( " and distance: ");
  ivp_message( qh_REAL_1, dist);
  ivp_message( "\n");
  return False;
} /* sethalfspace */

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="sethalfspace_all">-</a>
  
  qh_sethalfspace_all( dim, count, halfspaces, feasible )
    generate dual for halfspace intersection with feasible point
    array of count halfspaces
      each halfspace is normal coefficients followed by offset 
      the origin is inside the halfspace if the offset is negative

  returns:
    malloc'd array of count X dim-1 points

  notes:
    call before qh_init_B or qh_initqhull_globals 
    unused/untested code: please email [e-mail] if this works ok for you
    If using option 'Fp', also set qh feasible_point. It is a malloc'd array 
      that is freed by qh_freebuffers.

  design:
    see qh_sethalfspace
*/
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)) {
      ivp_message( "The halfspace was at index %d\n", i);
      qh_errexit (qh_ERRinput, NULL, NULL);
    }
    normalp= offsetp + 1;
  }
  return newpoints;
} /* sethalfspace_all */

  
/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="voronoi_center">-</a>
  
  qh_v