geom2.c

Quick hull implementation

      threshold= qh upper_threshold[k];
      if (normal[k] > threshold)
        within= False;
      if (angle) {
	threshold -= normal[k];
	*angle += fabs_(threshold);
      }
    }
  }
  return within;
} /* inthresholds */
    

/*-<a                             href="qh-geom.htm#TOC"
  >-------------------------------</a><a name="joggleinput">-</a>
  
  qh_joggleinput()
    randomly joggle input to Qhull by qh.JOGGLEmax
    initial input is qh.first_point/qh.num_points of qh.hull_dim
      repeated calls use qh.input_points/qh.num_points
 
  returns:
    joggles points at qh.first_point/qh.num_points
    copies data to qh.input_points/qh.input_malloc if first time
    determines qh.JOGGLEmax if it was zero
    if qh.DELAUNAY
      computes the Delaunay projection of the joggled points

  notes:
    if qh.DELAUNAY, unnecessarily joggles the last coordinate
    the initial 'QJn' may be set larger than qh_JOGGLEmaxincrease

  design:
    if qh.DELAUNAY
      set qh.SCALElast for reduced precision errors
    if first call
      initialize qh.input_points to the original input points
      if qh.JOGGLEmax == 0
        determine default qh.JOGGLEmax
    else
      increase qh.JOGGLEmax according to qh.build_cnt
    joggle the input by adding a random number in [-qh.JOGGLEmax,qh.JOGGLEmax]
    if qh.DELAUNAY
      sets the Delaunay projection
*/
void qh_joggleinput (void) {
  int size, i, seed;
  coordT *coordp, *inputp;
  realT randr, randa, randb;

  if (!qh input_points) { /* first call */
    qh input_points= qh first_point;
    qh input_malloc= qh POINTSmalloc;
    size= qh num_points * qh hull_dim * sizeof(coordT);
    if (!(qh first_point=(coordT*)malloc(size))) {
      fprintf(qh ferr, "qhull error: insufficient memory to joggle %d points\n",
          qh num_points);
      qh_errexit(qh_ERRmem, NULL, NULL);
    }
    qh POINTSmalloc= True;
    if (qh JOGGLEmax == 0.0) {
      qh JOGGLEmax= qh_detjoggle (qh input_points, qh num_points, qh hull_dim);
      qh_option ("QJoggle", NULL, &qh JOGGLEmax);
    }
  }else {                 /* repeated call */
    if (!qh RERUN && qh build_cnt > qh_JOGGLEretry) {
      if (((qh build_cnt-qh_JOGGLEretry-1) % qh_JOGGLEagain) == 0) {
	realT maxjoggle= qh MAXwidth * qh_JOGGLEmaxincrease;
	if (qh JOGGLEmax < maxjoggle) {
	  qh JOGGLEmax *= qh_JOGGLEincrease;
	  minimize_(qh JOGGLEmax, maxjoggle); 
	}
      }
    }
    qh_option ("QJoggle", NULL, &qh JOGGLEmax);
  }
  if (qh build_cnt > 1 && qh JOGGLEmax > fmax_(qh MAXwidth/4, 0.1)) {
      fprintf (qh ferr, "qhull error: the current joggle for 'QJn', %.2g, is too large for the width\nof the input.  If possible, recompile Qhull with higher-precision reals.\n",
	        qh JOGGLEmax);
      qh_errexit (qh_ERRqhull, NULL, NULL);
  }
  /* for some reason, using qh ROTATErandom and qh_RANDOMseed does not repeat the run. Use 'TRn' instead */
  seed= qh_RANDOMint;
  qh_option ("_joggle-seed", &seed, NULL);
  trace0((qh ferr, "qh_joggleinput: joggle input by %2.2g with seed %d\n", 
    qh JOGGLEmax, seed));
  inputp= qh input_points;
  coordp= qh first_point;
  randa= 2.0 * qh JOGGLEmax/qh_RANDOMmax;
  randb= -qh JOGGLEmax;
  size= qh num_points * qh hull_dim;
  for (i= size; i--; ) {
    randr= qh_RANDOMint;
    *(coordp++)= *(inputp++) + (randr * randa + randb);
  }
  if (qh DELAUNAY) {
    qh last_low= qh last_high= qh last_newhigh= REALmax;
    qh_setdelaunay (qh hull_dim, qh num_points, qh first_point);
  }
} /* joggleinput */

/*-<a                             href="qh-geom.htm#TOC"
  >-------------------------------</a><a name="maxabsval">-</a>
  
  qh_maxabsval( normal, dim )
    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 */


/*-<a                             href="qh-geom.htm#TOC"
  >-------------------------------</a><a name="maxmin">-</a>
  
  qh_maxmin( points, numpoints, dimension )
    return max/min points for each dimension      
    determine max and min coordinates

  returns:
    returns a temporary set of max and min points
      may include duplicate points. Does not include qh.GOODpoint
    sets qh.NEARzero, qh.MAXabs_coord, qh.MAXsumcoord, qh.MAXwidth
         qh.MAXlastcoord, qh.MINlastcoord
    initializes qh.max_outside, qh.min_vertex, qh.WAScoplanar, qh.ZEROall_ok

  notes:
    loop duplicated in qh_detjoggle()

  design:
    initialize global precision variables
    checks definition of REAL...
    for each dimension
      for each point
        collect maximum and minimum point
      collect maximum of maximums and minimum of minimums
      determine qh.NEARzero for Gaussian Elimination
*/
setT *qh_maxmin(pointT *points, int numpoints, int dimension) {
  int k;
  realT maxcoord, temp;
  pointT *minimum, *maximum, *point, *pointtemp;
  setT *set;

  qh max_outside= 0.0;
  qh MAXabs_coord= 0.0;
  qh MAXwidth= -REALmax;
  qh MAXsumcoord= 0.0;
  qh min_vertex= 0.0;
  qh WAScoplanar= False;
  if (qh ZEROcentrum)
    qh ZEROall_ok= True;
  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 + dimension;
    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;
    }
    if (k == dimension-1) {
      qh MINlastcoord= minimum[k];
      qh MAXlastcoord= maximum[k];
    }
    if (qh SCALElast && k == dimension-1)
      maxcoord= qh MAXwidth;
    else {
      maxcoord= fmax_(maximum[k], -minimum[k]);
      if (qh GOODpointp) {
        temp= fmax_(qh GOODpointp[k], -qh GOODpointp[k]);
        maximize_(maxcoord, temp);
      }
      temp= maximum[k] - minimum[k];
      maximize_(qh MAXwidth, temp);
    }
    maximize_(qh MAXabs_coord, maxcoord);
    qh MAXsumcoord += maxcoord;
    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 * qh MAXsumcoord * REALepsilon;
  }
  if (qh IStracing >=1)
    qh_printpoints (qh ferr, "qh_maxmin: found the max and min points (by dim):", set);
  return(set);
} /* maxmin */

/*-<a                             href="qh-geom.htm#TOC"
  >-------------------------------</a><a name="maxouter">-</a>

  qh_maxouter()
    return maximum distance from facet to outer plane
    normally this is qh.max_outside+qh.DISTround
    does not include qh.JOGGLEmax

  see:
    qh_outerinner()
    
  notes:
    need to add another qh.DISTround if testing actual point with computation

  for joggle:
    qh_setfacetplane() updated qh.max_outer for Wnewvertexmax (max distance to vertex)
    need to use Wnewvertexmax since could have a coplanar point for a high 
      facet that is replaced by a low facet
    need to add qh.JOGGLEmax if testing input points
*/
realT qh_maxouter (void) {
  realT dist;

  dist= fmax_(qh max_outside, qh DISTround);
  dist += qh DISTround;
  trace4((qh ferr, "qh_maxouter: max distance from facet to outer plane is %2.2g max_outside is %2.2g\n", dist, qh max_outside));
  return dist;
} /* maxouter */

/*-<a                             href="qh-geom.htm#TOC"
  >-------------------------------</a><a name="maxsimplex">-</a>
  
  qh_maxsimplex( dim, maxpoints, points, numpoints, simplex )
    determines maximum simplex for a set of points 
    starts from points already in simplex
    skips qh.GOODpointp (assumes that it isn't in maxpoints)
  
  returns:
    simplex with dim+1 points

  notes:
    assumes at least pointsneeded points in points
    maximizes determinate for x,y,z,w, etc.
    uses maxpoints as long as determinate is clearly non-zero

  design:
    initialize simplex with at least two points
      (find points with max or min x coordinate)
    for each remaining dimension
      add point that maximizes the determinate
        (use points from maxpoints first)    
*/
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) {
      qh_precision ("input has same x coordinate");
      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+1));
      }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+1, maxdet));
  } /* k */ 
} /* maxsimplex */

/*-<a                             href="qh-geom.htm#TOC"
  >-------------------------------</a><a name="minabsval">-</a>
  
  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-geom.htm#TOC"
  >-------------------------------</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-geom.htm#TOC"
  >-------------------------------</a><a name="orientoutside">-</a>
  
  qh_orientoutside( facet  )