geom2.c

这是一个新的用于图像处理的工具箱

      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 errors
returns:
  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 simplex
returns:
  temporary set of dim+1 points
notes:
  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 paraboloid
returns:
  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) {