geom.c

关于网格剖分的

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


/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="projectpoint">-</a>
  
  qh_projectpoint( point, facet, dist )
    project point onto a facet by dist

  returns:
    returns a new point
    
  notes:
    if dist= distplane(point,facet)
      this projects point to hyperplane
    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; /* used !qh_NOmem */
  
  qh_memalloc_(normsize, freelistp, newpoint, pointT);
  np= newpoint;
  normal= facet->normal;
  for(k= qh hull_dim; k--; )
    *(np++)= *point++ - dist * *normal++;
  return(newpoint);
} /* projectpoint */

  
/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="setfacetplane">-</a>
  
  qh_setfacetplane( facet )
    sets the hyperplane for a facet
    if qh.RANDOMdist, joggles hyperplane

  notes:
    uses global buffers qh.gm_matrix and qh.gm_row
    overwrites facet->normal if already defined
    updates Wnewvertex if PRINTstatistics
    sets facet->upperdelaunay if upper envelope of Delaunay triangulation

  design:
    copy vertex coordinates to qh.gm_matrix/gm_row
    compute determinate
    if nearzero
      recompute determinate with gaussian elimination
      if nearzero
        force outside orientation by testing interior point
*/
void qh_setfacetplane(facetT *facet) {
  pointT *point;
  vertexT *vertex, **vertexp;
  int k,i, normsize= qh normal_size, oldtrace= 0;
  realT dist;
  void **freelistp; /* used !qh_NOmem */
  coordT *coord, *gmcoord;
  pointT *point0= SETfirstt_(facet->vertices, vertexT)->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) {
      gmcoord= qh gm_matrix;
      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;
    gmcoord= qh gm_matrix;
    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. */
      }
    }
  }
  facet->upperdelaunay= False;
  if (qh DELAUNAY) {
    if (qh UPPERdelaunay) {     /* matches qh.lower_threshold in qh_initbuild */
      if (facet->normal[qh hull_dim -1] >= qh ANGLEround * qh_ZEROdelaunay)
        facet->upperdelaunay= True;
    }else {
      if (facet->normal[qh hull_dim -1] > -qh ANGLEround * qh_ZEROdelaunay)
        facet->upperdelaunay= True;
    }
  }
  if (qh PRINTstatistics || qh IStracing || qh TRACElevel || qh JOGGLEmax < REALmax) {
    qh old_randomdist= qh RANDOMdist;
    qh RANDOMdist= False;
    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 > wwval_(Wnewvertexmax)) {
          wwval_(Wnewvertexmax)= dist;
	  if (dist > qh max_outside) {
	    qh max_outside= dist;  /* used by qh_maxouter() */
	    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);
	}
      }
    }
    qh RANDOMdist= qh old_randomdist;
  }
  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 */


/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="sethyperplane_det">-</a>
  
  qh_sethyperplane_det( dim, rows, point0, toporient, normal, offset, nearzero )
    given dim X dim array indexed by rows[], one row per point, 
        toporient (flips all signs),
        and point0 (any row)
    set normalized hyperplane equation from oriented simplex

  returns:
    normal (normalized)
    offset (places point0 on the hyperplane)
    sets nearzero if hyperplane not through points

  notes:
    only defined for dim == 2..4
    rows[] is not modified
    solves det(P-V_0, V_n-V_0, ..., V_1-V_0)=0, i.e. every point is on hyperplane
    see Bower & Woodworth, A programmer's geometry, Butterworths 1983.

  derivation of 3-d minnorm
    Goal: all vertices V_i within qh.one_merge of hyperplane
    Plan: exactly translate the facet so that V_0 is the origin
          exactly rotate the facet so that V_1 is on the x-axis and y_2=0.
          exactly rotate the effective perturbation to only effect n_0
	     this introduces a factor of sqrt(3)
    n_0 = ((y_2-y_0)*(z_1-z_0) - (z_2-z_0)*(y_1-y_0)) / norm
    Let M_d be the max coordinate difference
    Let M_a be the greater of M_d and the max abs. coordinate
    Let u be machine roundoff and distround be max error for distance computation
    The max error for n_0 is sqrt(3) u M_a M_d / norm.  n_1 is approx. 1 and n_2 is approx. 0
    The max error for distance of V_1 is sqrt(3) u M_a M_d M_d / norm.  Offset=0 at origin
    Then minnorm = 1.8 u M_a M_d M_d / qh.ONEmerge
    Note that qh.one_merge is approx. 45.5 u M_a and norm is usually about M_d M_d

  derivation of 4-d minnorm
    same as above except rotate the facet so that V_1 on x-axis and w_2, y_3, w_3=0
     [if two vertices fixed on x-axis, can rotate the other two in yzw.]
    n_0 = det3_(...) = y_2 det2_(z_1, w_1, z_3, w_3) = - y_2 w_1 z_3
     [all other terms contain at least two factors nearly zero.]
    The max error for n_0 is sqrt(4) u M_a M_d M_d / norm
    Then minnorm = 2 u M_a M_d M_d M_d / qh.ONEmerge
    Note that qh.one_merge is approx. 82 u M_a and norm is usually about M_d M_d M_d
*/
void qh_sethyperplane_det (int dim, coordT **rows, coordT *point0, 
          boolT toporient, coordT *normal, realT *offset, boolT *nearzero) {
  realT maxround, dist;
  int i;
  pointT *point;


  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, NULL, NULL);
    *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
	       + point0[2]*normal[2]);
    maxround= qh DISTround;
    for (i=dim; i--; ) {
      point= rows[i];
      if (point != point0) {
        dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
	       + point[2]*normal[2]);
        if (dist > maxround || dist < -maxround) {
  	  *nearzero= True;
	  break;
	}
      }
    }
  }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, NULL, NULL);
    *offset= -(point0[0]*normal[0] + point0[1]*normal[1]
	       + point0[2]*normal[2] + point0[3]*normal[3]);
    maxround= qh DISTround;
    for (i=dim; i--; ) {
      point= rows[i];
      if (point != point0) {
        dist= *offset + (point[0]*normal[0] + point[1]*normal[1]
	       + point[2]*normal[2] + point[3]*normal[3]);
        if (dist > maxround || dist < -maxround) {
  	  *nearzero= True;
	  break;
	}
      }
    }
  }
  if (*nearzero) {
    zzinc_(Zminnorm);
    trace0((qh ferr, "qh_sethyperplane_det: degenerate norm during p%d.\n", qh furthest_id));
    zzinc_(Znearlysingular);
  }
} /* sethyperplane_det */


/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="sethyperplane_gauss">-</a>
  
  qh_sethyperplane_gauss( dim, rows, point0, toporient, normal, offset, nearzero )
    given (dim-1) X dim array of rows[i]= V_{i+1} - V_0 (point0)
    set normalized hyperplane equation from oriented simplex

  returns:
    normal (normalized)
    offset (places point0 on the hyperplane)

  notes:
    if nearzero
      orientation may be incorrect because of incorrect sign flips in gausselim
    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

  design:
    perform gaussian elimination
    flip sign for negative values
    perform back substitution 
    normalize result
    compute offset
*/
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 */