geom.c

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

  wmin_(Wmindenom, pivot_abs);  /* last pivot element */
  if (qh IStracing >= 5)
    qh_printmatrix (qh ferr, "qh_gausselem: result", rows, numrow, numcol);
} /* gausselim */


/*----------------------------------------------
-getangle- returns the dot product of two, qh hull_dim vectors
  may be > 1.0 or < -1.0
*/
realT qh_getangle(pointT *vect1, pointT *vect2) {
  realT angle= 0, randr;
  int k;

  for(k= qh hull_dim; k--; )
    angle += *vect1++ * *vect2++;
  if (qh RANDOMdist) {
    randr= qh_RANDOMint;
    angle += (2.0 * randr / qh_RANDOMmax - 1.0) *
      qh RANDOMfactor;
  }
  trace4((qh ferr, "qh_getangle: %2.2g\n", angle));
  return(angle);
} /* getangle */


/*----------------------------------------------
-getcenter-  gets arithmetic center of a set of vertices as a new point
*/
pointT *qh_getcenter(setT *vertices) {
  int k;
  pointT *center, *coord;
  vertexT *vertex, **vertexp;
  int count= qh_setsize(vertices);

  if (count < 2) {
    fprintf (qh ferr, "qhull internal error (qh_getcenter): not defined for %d points\n", count);
    qh_errexit (qh_ERRqhull, NULL, NULL);
  }
  center= (pointT *)qh_memalloc(qh normal_size);
  for (k=0; k < qh hull_dim; k++) {
    coord= center+k;
    *coord= 0.0;
    FOREACHvertex_(vertices)
      *coord += vertex->point[k];
    *coord /= count;
  }
  return(center);
} /* getcenter */


/*----------------------------------------------
-getcentrum- returns the centrum for a facet as a new point
  assumes qh_memfree_() is valid for normal_size
*/
pointT *qh_getcentrum(facetT *facet) {
  realT dist;
  pointT *centrum, *point;

  point= qh_getcenter(facet->vertices);
  zzinc_(Zcentrumtests);
  qh_distplane (point, facet, &dist);
  centrum= qh_projectpoint(point, facet, dist);
  qh_memfree(point, qh normal_size);
  trace4((qh ferr, "qh_getcentrum: for f%d, %d vertices dist= %2.2g\n",
	  facet->id, qh_setsize(facet->vertices), dist));
  return centrum;
} /* getcentrum */


/*-------------------------------------------------
-getdistance- returns the max and min distance of any vertex from neighbor
  returns the max absolute value
*/
realT qh_getdistance(facetT *facet, facetT *neighbor, realT *mindist, realT *maxdist) {
  vertexT *vertex, **vertexp;
  realT dist, maxd, mind;
  
  FOREACHvertex_(facet->vertices)
    vertex->seen= False;
  FOREACHvertex_(neighbor->vertices)
    vertex->seen= True;
  mind= 0.0;
  maxd= 0.0;
  FOREACHvertex_(facet->vertices) {
    if (!vertex->seen) {
      zzinc_(Zbestdist);
      qh_distplane(vertex->point, neighbor, &dist);
      if (dist < mind)
	mind= dist;
      else if (dist > maxd)
	maxd= dist;
    }
  }
  *mindist= mind;
  *maxdist= maxd;
  mind= -mind;
  if (maxd > mind)
    return maxd;
  else
    return mind;
} /* getdistance */


/*--------------------------------------------------
-normalize -- normalize a vector w/o min norm
*/
void qh_normalize (coordT *normal, int dim, boolT toporient) {
  qh_normalize2( normal, dim, toporient, NULL, NULL);
} /* normalize */

/*--------------------------------------------------
-normalize2 -- normalize a vector and report if too small
   qh MINdenom/MINdenom1 upper limits for divide overflow
   if minnorm non-NULL, sets ismin if normal is smaller
returns:
    normalized vector
    flips sign if !toporient
    if zero norm
       sets all elements to sqrt(1.0/dim)
    if divide by zero (divzero ())
       sets largest element to +/-1
       bumps Znearlysingular
*/
void qh_normalize2 (coordT *normal, int dim, boolT toporient, 
            realT *minnorm, boolT *ismin) {
  int k;
  realT *colp, *maxp, norm= 0, temp, *norm1, *norm2, *norm3;
  boolT zerodiv;

  norm1= normal+1;
  norm2= normal+2;
  norm3= normal+3;
  if (dim == 2)
    norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1));
  else if (dim == 3)
    norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2));
  else if (dim == 4) {
    norm= sqrt((*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2) 
               + (*norm3)*(*norm3));
  }else if (dim > 4) {
    norm= (*normal)*(*normal) + (*norm1)*(*norm1) + (*norm2)*(*norm2) 
               + (*norm3)*(*norm3);
    for (k= dim-4, colp= normal+4; k--; colp++)
      norm += (*colp) * (*colp);
    norm= sqrt(norm);
  }
  if (minnorm) {
    if (norm < *minnorm) 
      *ismin= True;
    else
      *ismin= False;
  }
  wmin_(Wmindenom, norm);
  if (norm > qh MINdenom) {
    if (!toporient)
      norm= -norm;
    *normal /= norm;
    *norm1 /= norm;
    if (dim == 2)
      ; /* all done */
    else if (dim == 3)
      *norm2 /= norm;
    else if (dim == 4) {
      *norm2 /= norm;
      *norm3 /= norm;
    }else if (dim >4) {
      *norm2 /= norm;
      *norm3 /= norm;
      for (k= dim-4, colp= normal+4; k--; )
        *colp++ /= norm;
    }
  }else if (norm == 0.0) {
    temp= sqrt (1.0/dim);
    for (k= dim, colp= normal; k--; )
      *colp++ = temp;
  }else {
    if (!toporient)
      norm= -norm;
    for (k= dim, colp= normal; k--; colp++) { /* k used below */
      temp= qh_divzero (*colp, norm, qh MINdenom_1, &zerodiv);
      if (!zerodiv)
	*colp= temp;
      else {
	maxp= qh_maxabsval(normal, dim);
	temp= ((*maxp * norm >= 0.0) ? 1.0 : -1.0);
	for (k= dim, colp= normal; k--; colp++)
	  *colp= 0.0;
	*maxp= temp;
	zzinc_(Znearlysingular);
	trace0((qh ferr, "qh_normalize: norm=%2.2g too small during p%d\n", 
	       norm, qh furthest_id));
	return;
      }
    }
  }
} /* normalize */


/*-------------------------------------------------
-projectpoint- project point onto a facet by dist
  projects point to hyperplane if dist= distplane(point,facet)
returns:
  returns a new point
  assumes qh_memfree_() is valid for normal_size
*/

pointT *qh_projectpoint(pointT *point, facetT *facet, realT dist) {
  pointT *newpoint, *np, *normal;
  int normsize= qh normal_size,k;
  void **freelistp;
  
  qh_memalloc_(normsize, freelistp, newpoint, pointT);
  np= newpoint;
  normal= facet->normal;
  for(k= qh hull_dim; k--; )
    *(np++)= *point++ - dist * *normal++;
  return(newpoint);
} /* projectpoint */

  
/*-------------------------------------------------
-setfacetplane- sets the hyperplane for a facet
   uses global buffers qh gm_matrix and qh gm_row
   overwrites facet->normal if already defined
   sets facet->upperdelaunay if upper envelope of Delaunay triangulation
   updates Wnewvertex if PRINTstatistics
*/
void qh_setfacetplane(facetT *facet) {
  pointT *point;
  vertexT *vertex, **vertexp;
  int k,i, normsize= qh normal_size, oldtrace= 0;
  realT dist;
  void **freelistp;
  coordT *coord, *gmcoord= qh gm_matrix;
  pointT *point0= ((vertexT*)SETfirst_(facet->vertices))->point;
  boolT nearzero= False;

  zzinc_(Zsetplane);
  if (!facet->normal)
    qh_memalloc_(normsize, freelistp, facet->normal, coordT);
  if (facet == qh tracefacet) {
    oldtrace= qh IStracing;
    qh IStracing= 5;
    fprintf (qh ferr, "qh_setfacetplane: facet f%d created.\n", facet->id);
    fprintf (qh ferr, "  Last point added to hull was p%d.", qh furthest_id);
    if (zzval_(Ztotmerge))
      fprintf(qh ferr, "  Last merge was #%d.", zzval_(Ztotmerge));
    fprintf (qh ferr, "\n\nCurrent summary is:\n");
      qh_printsummary (qh ferr);
  }
  if (qh hull_dim <= 4) {
    i= 0;
    if (qh RANDOMdist) {
      FOREACHvertex_(facet->vertices) {
        qh gm_row[i++]= gmcoord;
	coord= vertex->point;
	for (k= qh hull_dim; k--; )
	  *(gmcoord++)= *coord++ * qh_randomfactor();
      }	  
    }else {
      FOREACHvertex_(facet->vertices)
       qh gm_row[i++]= vertex->point;
    }
    qh_sethyperplane_det(qh hull_dim, qh gm_row, point0, facet->toporient,
                facet->normal, &facet->offset, &nearzero);
  }
  if (qh hull_dim > 4 || nearzero) {
    i= 0;
    FOREACHvertex_(facet->vertices) {
      if (vertex->point != point0) {
	qh gm_row[i++]= gmcoord;
	coord= vertex->point;
	point= point0;
	for(k= qh hull_dim; k--; )
	  *(gmcoord++)= *coord++ - *point++;
      }
    }
    qh gm_row[i]= gmcoord;  /* for areasimplex */
    if (qh RANDOMdist) {
      gmcoord= qh gm_matrix;
      for (i= qh hull_dim-1; i--; ) {
	for (k= qh hull_dim; k--; )
	  *(gmcoord++) *= qh_randomfactor();
      }
    }
    qh_sethyperplane_gauss(qh hull_dim, qh gm_row, point0, facet->toporient,
           	facet->normal, &facet->offset, &nearzero);
    if (nearzero) { 
      if (qh_orientoutside (facet)) {
	trace0((qh ferr, "qh_setfacetplane: flipped orientation after testing interior_point during p%d\n", qh furthest_id));
      /* this is part of using Gaussian Elimination.  For example in 5-d
	   1 1 1 1 0
	   1 1 1 1 1
	   0 0 0 1 0
	   0 1 0 0 0
	   1 0 0 0 0
	   norm= 0.38 0.38 -0.76 0.38 0
	 has a determinate of 1, but g.e. after subtracting pt. 0 has
	 0's in the diagonal, even with full pivoting.  It does work
	 if you subtract pt. 4 instead. */
      }
    }
  }
  if (qh DELAUNAY && facet->normal[qh hull_dim -1] > qh ANGLEround)
    facet->upperdelaunay= True;
  else
    facet->upperdelaunay= False;

  if (qh PRINTstatistics || qh IStracing || qh TRACElevel) {
    FOREACHvertex_(facet->vertices) {
      if (vertex->point != point0) {
	boolT istrace= False;
	zinc_(Zdiststat);
        qh_distplane(vertex->point, facet, &dist);
        dist= fabs_(dist);
        zinc_(Znewvertex);
        wadd_(Wnewvertex, dist);
        if (dist > wval_(Wnewvertexmax)) {
          wval_(Wnewvertexmax)= dist;
	  if (dist > qh max_outside) {
	    qh max_outside= dist;
	    if (dist > qh TRACEdist) 
	      istrace= True;
	  }
	}else if (-dist > qh TRACEdist)
	  istrace= True;
	if (istrace) {
	  fprintf (qh ferr, "qh_setfacetplane: ====== vertex p%d (v%d) increases max_outside to %2.2g for new facet f%d last p%d\n",
	        qh_pointid(vertex->point), vertex->id, dist, facet->id, qh furthest_id);
	  qh_errprint ("DISTANT", facet, NULL, NULL, NULL);
	}
      }
    }
  }
  if (qh IStracing >= 3) {
    fprintf (qh ferr, "qh_setfacetplane: f%d offset %2.2g normal: ",
	     facet->id, facet->offset);
    for (k=0; k<qh hull_dim; k++)
      fprintf (qh ferr, "%2.2g ", facet->normal[k]);
    fprintf (qh ferr, "\n");
  }
  if (facet == qh tracefacet)
    qh IStracing= oldtrace;
} /* setfacetplane */


/*-------------------------------------------------
-sethyperplane_det- set normalized hyperplane equation from oriented simplex
  dim X dim array indexed by rows[], one row per point, point0 is any row
  only defined for dim == 2..4
returns:
  offset, normal
  bumps Znearlysingular if normalization fails
  rows[] is not modified
notes:
  solves det(P-V_0, V_n-V_0, ..., V_1-V_0)=0, i.e. every point is on hyperplane
  offset places point0 on the hyperplane
  toporient just flips all signs, so orientation is correct
  see Bower & Woodworth, A programmer's geometry, Butterworths 1983.
*/
void qh_sethyperplane_det (int dim, coordT **rows, coordT *point0, 
          boolT toporient, coordT *normal, realT *offset, boolT *nearzero) {

  if (dim == 2) {
    normal[0]= dY(1,0);
    normal[1]= dX(0,1);
    qh_normalize2 (normal, dim, toporient, NULL, NULL);
    *offset= -(point0[0]*normal[0]+point0[1]*normal[1]);
    *nearzero= False;  /* since nearzero norm => incident points */
  }else if (dim == 3) {
    normal[0]= det2_(dY(2,0), dZ(2,0),
		     dY(1,0), dZ(1,0));
    normal[1]= det2_(dX(1,0), dZ(1,0),
		     dX(2,0), dZ(2,0));
    normal[2]= det2_(dX(2,0), dY(2,0),
		     dX(1,0), dY(1,0));
    qh_normalize2 (normal, dim, toporient, &qh MINnorm, nearzero);
    *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
	       + point0[2]*normal[2]);
  }else if (dim == 4) {
    normal[0]= - det3_(dY(2,0), dZ(2,0), dW(2,0),
			dY(1,0), dZ(1,0), dW(1,0),
			dY(3,0), dZ(3,0), dW(3,0));
    normal[1]=   det3_(dX(2,0), dZ(2,0), dW(2,0),
		        dX(1,0), dZ(1,0), dW(1,0),
		        dX(3,0), dZ(3,0), dW(3,0));
    normal[2]= - det3_(dX(2,0), dY(2,0), dW(2,0),
			dX(1,0), dY(1,0), dW(1,0),
			dX(3,0), dY(3,0), dW(3,0));
    normal[3]=   det3_(dX(2,0), dY(2,0), dZ(2,0),
		        dX(1,0), dY(1,0), dZ(1,0),
		        dX(3,0), dY(3,0), dZ(3,0));
    qh_normalize2 (normal, dim, toporient, &qh MINnorm, nearzero);
    *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
	       + point0[2]*normal[2] + point0[3]*normal[3]);
  }
  if (*nearzero) {
    zzinc_(Zminnorm);
    trace0((qh ferr, "qh_sethyperplane_det: degenerate norm during p%d.\n", qh furthest_id));
  }
} /* sethyperplane_det */


/*-------------------------------------------------
-sethyperplane_gauss- set normalized hyperplane equation from oriented simplex
    (dim-1) X dim array of rows[i]= V_{i+1} - V_0 (point0)
returns:
    offset, normal
    if nearzero, bumps Znearlysingular
      orientation may be incorrect because of incorrect sign flips in gausselim
notes:
    solves [V_n-V_0,...,V_1-V_0, 0 .. 0 1] * N == [0 .. 0 1] 
        or [V_n-V_0,...,V_1-V_0, 0 .. 0 1] * N == [0] 
    i.e., N is normal to the hyperplane, and the unnormalized
        distance to [0 .. 1] is either 1 or 0
    offset places point0 on the hyperplane
*/
void qh_sethyperplane_gauss (int dim, coordT **rows, pointT *point0, 
		boolT toporient, coordT *normal, coordT *offset, boolT *nearzero) {
  coordT *pointcoord, *normalcoef;
  int k;
  boolT sign= toporient, nearzero2= False;
  
  qh_gausselim(rows, dim-1, dim, &sign, nearzero);
  for(k= dim-1; k--; ) {
    if ((rows[k])[k] < 0)
      sign ^= 1;
  }
  if (*nearzero) {
    zzinc_(Znearlysingular);
    trace0((qh ferr, "qh_sethyperplane_gauss: nearly singular or axis parallel hyperplane during p%d.\n", qh furthest_id));
    qh_backnormal(rows, dim-1, dim, sign, normal, &nearzero2);
  }else {
    qh_backnormal(rows, dim-1, dim, sign, normal, &nearzero2);
    if (nearzero2) {
      zzinc_(Znearlysingular);
      trace0((qh ferr, "qh_sethyperplane_gauss: singular or axis parallel hyperplane at normalization during p%d.\n", qh furthest_id));
    }
  }
  if (nearzero2)
    *nearzero= True;
  qh_normalize2(normal, dim, True, NULL, NULL);
  pointcoord= point0;
  normalcoef= normal;
  *offset= -(*pointcoord++ * *normalcoef++);
  for(k= dim-1; k--; )
    *offset -= *pointcoord++ * *normalcoef++;
} /* sethyperplane_gauss */