geom2.c

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

    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+1Xdim
returns:
  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 buffer
returns:
  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 array
returns:
  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 negative
returns:
  unused/untested code: please email [e-mail] 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 points
notes:
  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 */