geom.c

关于网格剖分的

    dist is distance to facet
    isoutside is true if point is outside of facet
    numpart is number of distance tests

  notes:
    if qh.BESToutside or !isoutside
       stops at furthest facet
    if qh.MERGING 
       stops when distance > qh_DISToutside (max(4*MINoutside, 2*max_outside))
    else 
       stops when distance > MINoutside (DISTround in precise case)
    searches newfacets then searchs neighbors of best facet.
    avoids upperdelaunay facet unless none other or (isoutside and outside)

    uses visit_id and seen flags
    caller traces the results
  
  see also:
    qh_partitionall() and qh_findbest()

  design:
    for each new facet
      test distance from point to facet
      select best facet
    test each horizon facet of best facet
    return best facet 
*/
facetT *qh_findbestnew (pointT *point, facetT *startfacet,
	   realT *dist, boolT *isoutside, int *numpart) {
  realT bestdist= -REALmax, bestdist2= -REALmax;
  facetT *neighbor, **neighborp, *bestfacet= NULL, *newfacet, *facet;
  facetT *bestfacet2= NULL;
  int oldtrace= qh IStracing, i;
  realT distoutside;

  if (!startfacet) {
    if (qh MERGING)
      fprintf(qh ferr, "qhull precision error (qh_findbestnew): merging has formed and deleted an independent cycle of facets.  Can not continue.\n");
    else
      fprintf(qh ferr, "qhull internal error (qh_findbestnew): no new facets for point p%d\n",
      	      qh furthest_id);      
    qh_errexit (qh_ERRqhull, NULL, NULL);
  }
  if (qh BESToutside || !isoutside)
    distoutside= REALmax;
  else if (qh MERGING)
    distoutside= qh_DISToutside; /* defined in user.h */
  else
    distoutside= qh MINoutside;
  if (qh TRACElevel && qh TRACEpoint >= 0 && qh TRACEpoint == qh_pointid (point)) {
    qh IStracing= qh TRACElevel;
    fprintf(qh ferr, "qh_findbestnew: point p%d facet f%d. Stop if dist > %2.2g\n",
	     qh TRACEpoint, startfacet->id, distoutside);
    fprintf(qh ferr, "  Last point added to hull was p%d.", qh furthest_id);
    fprintf(qh ferr, "  Last merge was #%d.\n", zzval_(Ztotmerge));
  }
  if (isoutside)
    *isoutside= True;
  *numpart= 0;

  /* visit all new facets starting with startfacet */
  for (i= 0, facet= startfacet; i < 2; i++, facet= qh newfacet_list) {
    FORALLfacet_(facet) {
      if (facet == startfacet && i)
	break;
      qh_distplane (point, facet, dist);
      (*numpart)++;
      if (facet->upperdelaunay) {
	if (*dist > bestdist2) {
	  bestdist2= *dist;
	  bestfacet2= facet;
	  if (*dist >= distoutside) {
	    bestfacet= facet;
	    goto LABELreturn_bestnew;
	  }
	}
      }else if (*dist > bestdist) {
	bestdist= *dist;
	bestfacet= facet;
	if (*dist >= distoutside) 
	  goto LABELreturn_bestnew;
      }
    }
  }
  newfacet= bestfacet ? bestfacet : bestfacet2;
  	/* !bestfacet only occurs if 'd' creates incorrect upper-delaunay facets */
  FOREACHneighbor_(newfacet) {
    if (!neighbor->newfacet) {
      qh_distplane (point, neighbor, dist);
      (*numpart)++;
      if (neighbor->upperdelaunay) {
	if (*dist > bestdist2) {
	  bestdist2= *dist;
	  bestfacet2= neighbor;
	}
      }else if (*dist > bestdist) {
	bestdist= *dist;
	bestfacet= neighbor;
      }
    }
  }
  if (!bestfacet  
  || (isoutside && bestdist2 >= qh MINoutside && bestdist2 > bestdist)) {
    *dist= bestdist2;
    bestfacet= bestfacet2;
  }else 
    *dist= bestdist;
  if (isoutside && *dist < qh MINoutside)
    *isoutside= False;
LABELreturn_bestnew:
  qh IStracing= oldtrace;
  return bestfacet;
}  /* findbestnew */


/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="gausselim">-</a>
  
  qh_gausselim( rows, numrow, numcol, sign )
    Gaussian elimination with partial pivoting

  returns:
    rows is upper triangular (includes row exchanges)
    flips sign for each row exchange
    sets nearzero if pivot[k] < qh.NEARzero[k], else clears it

  notes:
    if nearzero, the determinant's sign may be incorrect.
    assumes numrow <= numcol

  design:
    for each row
      determine pivot and exchange rows if necessary
      test for near zero
      perform gaussian elimination step
*/
void qh_gausselim(realT **rows, int numrow, int numcol, boolT *sign, boolT *nearzero) {
  realT *ai, *ak, *rowp, *pivotrow;
  realT n, pivot, pivot_abs= 0.0, temp;
  int i, j, k, pivoti, flip=0;
  
  *nearzero= False;
  for(k= 0; k < numrow; k++) {
    pivot_abs= fabs_((rows[k])[k]);
    pivoti= k;
    for(i= k+1; i < numrow; i++) {
      if ((temp= fabs_((rows[i])[k])) > pivot_abs) {
	pivot_abs= temp;
	pivoti= i;
      }
    }
    if (pivoti != k) {
      rowp= rows[pivoti]; 
      rows[pivoti]= rows[k]; 
      rows[k]= rowp; 
      *sign ^= 1;
      flip ^= 1;
    }
    if (pivot_abs <= qh NEARzero[k]) {
      *nearzero= True;
      if (pivot_abs == 0.0) {   /* remainder of column == 0 */
	if (qh IStracing >= 4) {
	  fprintf (qh ferr, "qh_gausselim: 0 pivot at column %d. (%2.2g < %2.2g)\n", k, pivot_abs, qh DISTround);
	  qh_printmatrix (qh ferr, "Matrix:", rows, numrow, numcol);
	}
	zzinc_(Zgauss0);
        qh_precision ("zero pivot for Gaussian elimination");
	goto LABELnextcol;
      }
    }
    pivotrow= rows[k] + k;
    pivot= *pivotrow++;  /* signed value of pivot, and remainder of row */
    for(i= k+1; i < numrow; i++) {
      ai= rows[i] + k;
      ak= pivotrow;
      n= (*ai++)/pivot;   /* divzero() not needed since |pivot| >= |*ai| */
      for(j= numcol - (k+1); j--; )
	*ai++ -= n * *ak++;
    }
  LABELnextcol:
    ;
  }
  wmin_(Wmindenom, pivot_abs);  /* last pivot element */
  if (qh IStracing >= 5)
    qh_printmatrix (qh ferr, "qh_gausselem: result", rows, numrow, numcol);
} /* gausselim */


/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="getangle">-</a>
  
  qh_getangle( vect1, vect2 )
    returns the dot product of two, DIM3 vectors
    if qh.RANDOMdist, joggles result

  notes:
    the angle may be > 1.0 or < -1.0 because of roundoff errors

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


/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="getcenter">-</a>
  
  qh_getcenter( vertices )
    returns arithmetic center of a set of vertices as a new point

  notes:
    allocates point array for center
*/
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 */


/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="getcentrum">-</a>
  
  qh_getcentrum( facet )
    returns the centrum for a facet as a new point

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


/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="getdistance">-</a>
  
  qh_getdistance( facet, neighbor, mindist, maxdist )
    returns the maxdist and mindist distance of any vertex from neighbor

  returns:
    the max absolute value

  design:
    for each vertex of facet that is not in neighbor
      test the distance from vertex to neighbor
*/
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 */


/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="normalize">-</a>

  qh_normalize( normal, dim, toporient )
    normalize a vector and report if too small
    does not use min norm
  
  see:
    qh_normalize2
*/
void qh_normalize (coordT *normal, int dim, boolT toporient) {
  qh_normalize2( normal, dim, toporient, NULL, NULL);
} /* normalize */

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="normalize2">-</a>
  
  qh_normalize2( normal, dim, toporient, minnorm, ismin )
    normalize a vector and report if too small
    qh.MINdenom/MINdenom1 are the upper limits for divide overflow

  returns:
    normalized vector
    flips sign if !toporient
    if minnorm non-NULL, 
      sets ismin if normal < minnorm

  notes:
    if zero norm
       sets all elements to sqrt(1.0/dim)
    if divide by zero (divzero ())
       sets largest element to   +/-1
       bumps Znearlysingular
      
  design:
    computes norm
    test for minnorm
    if not near zero
      normalizes normal
    else if zero norm
      sets normal to standard value
    else
      uses qh_divzero to normalize
      if nearzero
        sets norm to direction of maximum value
*/
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)