qhull_geom2.cxx

hl2 source code. Do not use it illegal.

  
  qh_minabsval( normal, dim )
    return minimum 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 */


/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="mindiff">-</a>
  
  qh_mindif( vecA, vecB, dim )
    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= 0;

  for (k= 0; k < dim; k++) {
    diff= *vecAp++ - *vecBp++;
    diff= fabs_(diff);
    if (diff < mindiff) {
      mindiff= diff;
      mink= k;
    }
  }
  return mink;
} /* mindiff */



/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="orientoutside">-</a>
  
  qh_orientoutside( facet  )
    make facet outside oriented via qh.interior_point

  returns:
    True if facet 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 */

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="outerinner">-</a>
  
  qh_outerinner( facet, outerplane, innerplane  )
    if facet
      returns outer and inner plane for facet
      requires qh.maxoutdone, i.e., qh_check_maxout()
    else
      returns maximum outer and inner plane
    accounts for qh.JOGGLEmax

  see:
    qh_maxouter(), qh_check_bestdist(), qh_check_points()

  notes:
    outerplaner or innerplane may be NULL
    
    includes qh.DISTround for actual points
    add another qh.DISTround if testing with floating point arithmetic
*/
void qh_outerinner (facetT *facet, realT *outerplane, realT *innerplane) {
  realT dist, mindist;
  vertexT *vertex, **vertexp;

  if (outerplane) {
    if (!qh_MAXoutside || !facet || !qh maxoutdone) {
      *outerplane= qh_maxouter();       /* includes qh.DISTround */
    }else { /* qh_MAXoutside ... */
#if qh_MAXoutside 
      *outerplane= facet->maxoutside + qh DISTround;
#endif
      
    }
    if (qh JOGGLEmax < REALmax/2)
      *outerplane += qh JOGGLEmax * sqrt (qh hull_dim);
  }
  if (innerplane) {
    if (facet) {
      mindist= REALmax;
      FOREACHvertex_(facet->vertices) {
        zinc_(Zdistio);
        qh_distplane (vertex->point, facet, &dist);
        minimize_(mindist, dist);
      }
      *innerplane= mindist - qh DISTround;
    }else 
      *innerplane= qh min_vertex - qh DISTround;
    if (qh JOGGLEmax < REALmax/2)
      *innerplane -= qh JOGGLEmax * sqrt (qh hull_dim);
  }
} /* outerinner */

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="pointdist">-</a>
  
  qh_pointdist( point1, point2, dim )
    return distance between two points

  notes:
    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 */


/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="printmatrix">-</a>
  
  qh_printmatrix( fp, string, rows, numrow, numcol )
    print matrix to fp given by row vectors
    print string as header

  notes:
    print a vector by qh_printmatrix(fp, "", &vect, 1, len)
*/
void qh_printmatrix (FILE *fp, const 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++) {
      r= *rowp++;
      fprintf (fp, "%6.3g ", r);
    }
    fprintf (fp, "\n");
  }
} /* printmatrix */

  
/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="printpoints">-</a>
  
  qh_printpoints( fp, string, points )
    print pointids to fp for a set of points
    if string, prints string and 'p' point ids
*/
void qh_printpoints (FILE *fp, const 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 */

  
/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="projectinput">-</a>
  
  qh_projectinput()
    project input points using qh.lower_bound/upper_bound and qh DELAUNAY
    if qh.lower_bound[k]=qh.upper_bound[k]= 0, 
      removes dimension k 
    input points in qh first_point, num_points, input_dim

  returns:
    new point array in qh first_point of qh hull_dim coordinates
    sets qh POINTSmalloc
    if qh DELAUNAY 
      projects points to paraboloid
      lowbound/highbound is also projected
    if qh ATinfinity
      adds point "at-infinity"
    if qh POINTSmalloc 
      frees old point array

  notes:
    checks that qh.hull_dim agrees with qh.input_dim, PROJECTinput, and DELAUNAY


  design:
    sets project[k] to -1 (delete), 0 (keep), 1 (add for Delaunay)
    determines newdim and newnum for qh hull_dim and qh num_points
    projects points to newpoints
    projects qh.lower_bound to itself
    projects qh.upper_bound to itself
    if qh DELAUNAY
      if qh ATINFINITY
        projects points to paraboloid
        computes "infinity" point as vertex average and 10% above all points 
      else
        uses qh_setdelaunay to project points to paraboloid
*/
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*)p_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)
    P_FREE (qh first_point);
  qh first_point= newpoints;
  qh POINTSmalloc= True;
  if (qh DELAUNAY && qh ATinfinity) {
    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);
    }
    /* coord == infinity */
    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"));
  }else if (qh DELAUNAY)  /* !qh ATinfinity */
    qh_setdelaunay( qh hull_dim, qh num_points, qh first_point);
} /* projectinput */

  
/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="projectpoints">-</a>
  
  qh_projectpoints( project, n, points, numpoints, dim, newpoints, newdim )
    project points/numpoints/dim to newpoints/newdim
    if project[k] == -1
      delete dimension k 
    if project[k] == 1 
      add dimension k by duplicating previous column
    n is size of project

  notes:
    newpoints may be points if only adding dimension at end

  design:
    check that 'project' and 'newdim' agree
    for each dimension
      if project == -1
        skip dimension
      else
        determine start of column in newpoints
        determine start of column in points 
          if project == +1, duplicate previous column
        copy dimension (column) from points to newpoints
*/
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) {
    ivp_message( "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 */
        

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="rand">-</a>
  
  qh_rand() 
  qh_srand( seed )
    generate pseudo-random number between 1 and 2^31 -2

  notes:
    from Park & Miller's minimimal standard random number generator
       Communications of the ACM, 31:1192-1201, 1988.
    does not use 0 or 2^31 -1
       this is silently enforced by qh_srand()
    can make 'Rn' much faster by moving qh_rand to qh_distplane
*/
int qh_rand_seed= 1;  /* define as global variable instead of using qh */

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;
  /* seed = seed < qh_RANDOMmax/2 ? 0 : qh_RANDOMmax;  for testing */
  /* seed = qh_RANDOMmax;  for testing */
  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 */

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="randomfactor">-</a>
  
  qh_randomfactor()
    return a random factor within qh.RANDOMmax of 1.0

  notes:
    qh.RANDOMa/b are defined in global.c