poly2.c

关于网格剖分的

    if bestoutside is set (e.g., qh_ALL)
      returns best facet that is not upperdelaunay
      if Delaunay and inside, point is outside circumsphere of bestfacet
    else
      returns first facet below point
      if point is inside, returns nearest, !upperdelaunay facet
    distance to facet
    isoutside set if outside of facet
    
  notes:
    this works for all distributions
    if inside, qh_findbestfacet performed an exhaustive search
    qh_findbestfacet is not used by qhull.
    uses qh.visit_id, qh.searchset
    
  see:
    <a href="geom.c#findbest">qh_findbest</a>
*/
facetT *qh_findbestfacet (pointT *point, boolT bestoutside,
           realT *bestdist, boolT *isoutside) {
  facetT *bestfacet= NULL;
  int numpart, totpart= 0;
  
  bestfacet= qh_findbest (point, qh facet_list, 
			    bestoutside, False, bestoutside,
			    bestdist, isoutside, &totpart);
  if (!bestfacet) {
    fprintf (qh ferr, "qh_findbestfacet: all facets are flipped or upper Delaunay\n");
    qh_errexit (qh_ERRqhull, NULL, NULL);
  }
  if (*bestdist < -qh DISTround) {
    bestfacet= qh_findfacet_all (point, bestdist, isoutside, &numpart);
    totpart += numpart;
    if ((isoutside && bestoutside)
    || (!isoutside && bestfacet->upperdelaunay)) {
      bestfacet= qh_findbest (point, bestfacet, 
			    bestoutside, False, bestoutside,
			    bestdist, isoutside, &totpart);
      totpart += numpart;
    }
  }
  trace3((qh ferr, "qh_findbestfacet: f%d dist %2.2g isoutside %d totpart %d\n",
	  bestfacet->id, *bestdist, *isoutside, totpart));
  return bestfacet;
} /* findbestfacet */ 
 
/*-<a                             href="qh-c.htm"
  >-------------------------------</a><a name="findfacet_all">-</a>
  
  qh_findfacet_all( point, bestdist, isoutside, numpart )
    exhaustive search for facet below a point 

    for Delaunay triangulations, 
      Use qh_setdelaunay() to lift point to paraboloid and scale by 'Qbb' if needed
      Do not use options 'Qbk', 'QBk', or 'QbB' since they scale the coordinates. 

  returns:
    returns first facet below point
    if point is inside, 
      returns nearest facet
    distance to facet
    isoutside if point is outside of the hull
    number of distance tests
*/
facetT *qh_findfacet_all (pointT *point, realT *bestdist, boolT *isoutside,
			  int *numpart) {
  facetT *bestfacet= NULL, *facet;
  realT dist;
  int totpart= 0;
  
  *bestdist= REALmin;
  *isoutside= False;
  FORALLfacets {
    if (facet->flipped || !facet->normal)
      continue;
    totpart++;
    qh_distplane (point, facet, &dist);
    if (dist > *bestdist) {
      *bestdist= dist;
      bestfacet= facet;
      if (dist > qh MINoutside) {
        *isoutside= True;
        break;
      }
    }
  }
  *numpart= totpart;
  trace3((qh ferr, "qh_findfacet_all: f%d dist %2.2g isoutside %d totpart %d\n",
	  getid_(bestfacet), *bestdist, *isoutside, totpart));
  return bestfacet;
} /* findfacet_all */ 
 
/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="findgood">-</a>
  
  qh_findgood( facetlist, goodhorizon )
    identify good facets for qh.PRINTgood
    if qh.GOODvertex>0
      facet includes point as vertex
      if !match, returns goodhorizon
      inactive if qh.MERGING
    if qh.GOODpoint
      facet is visible or coplanar (>0) or not visible (<0) 
    if qh.GOODthreshold
      facet->normal matches threshold
    if !goodhorizon and !match, 
      selects facet with closest angle
      sets GOODclosest
      
  returns:
    number of new, good facets found
    determines facet->good
    may update qh.GOODclosest
    
  notes:
    qh_findgood_all further reduces the good region

  design:
    count good facets
    mark good facets for qh.GOODpoint  
    mark good facets for qh.GOODthreshold
    if necessary
      update qh.GOODclosest  
*/
int qh_findgood (facetT *facetlist, int goodhorizon) {
  facetT *facet, *bestfacet= NULL;
  realT angle, bestangle= REALmax, dist;
  int  numgood=0;

  FORALLfacet_(facetlist) {
    if (facet->good)
      numgood++;
  }
  if (qh GOODvertex>0 && !qh MERGING) {
    FORALLfacet_(facetlist) {
      if (!qh_isvertex (qh GOODvertexp, facet->vertices)) {
        facet->good= False;
        numgood--;
      }
    }
  }
  if (qh GOODpoint && numgood) {
    FORALLfacet_(facetlist) {
      if (facet->good && facet->normal) {
        zinc_(Zdistgood);
        qh_distplane (qh GOODpointp, facet, &dist);
        if ((qh GOODpoint > 0) ^ (dist > 0.0)) {
          facet->good= False;
          numgood--;
        }
      }
    }
  }
  if (qh GOODthreshold && (numgood || goodhorizon || qh GOODclosest)) {
    FORALLfacet_(facetlist) {
      if (facet->good && facet->normal) {
        if (!qh_inthresholds (facet->normal, &angle)) {
          facet->good= False;
          numgood--;
          if (angle < bestangle) {
            bestangle= angle;
            bestfacet= facet;
          }
        }
      }
    }
    if (!numgood && (!goodhorizon || qh GOODclosest)) {
      if (qh GOODclosest) {
	if (qh GOODclosest->visible)
	  qh GOODclosest= NULL;
	else {
	  qh_inthresholds (qh GOODclosest->normal, &angle);
	  if (angle < bestangle)
	    bestfacet= qh GOODclosest;
	}
      }
      if (bestfacet && bestfacet != qh GOODclosest) {
	if (qh GOODclosest)
	  qh GOODclosest->good= False;
	qh GOODclosest= bestfacet;
	bestfacet->good= True;
	numgood++;
	trace2((qh ferr, "qh_findgood: f%d is closest (%2.2g) to thresholds\n", 
           bestfacet->id, bestangle));
	return numgood;
      }
    }else if (qh GOODclosest) { /* numgood > 0 */
      qh GOODclosest->good= False;
      qh GOODclosest= NULL;
    }
  }
  zadd_(Zgoodfacet, numgood);
  trace2((qh ferr, "qh_findgood: found %d good facets with %d good horizon\n",
               numgood, goodhorizon));
  if (!numgood && qh GOODvertex>0 && !qh MERGING) 
    return goodhorizon;
  return numgood;
} /* findgood */

/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="findgood_all">-</a>
  
  qh_findgood_all( facetlist )
    apply other constraints for good facets (used by qh.PRINTgood)
    if qh.GOODvertex 
      facet includes (>0) or doesn't include (<0) point as vertex
      if last good facet, prints warning and continues
    if qh.SPLITthresholds
      facet->normal matches threshold, or if none, the closest one
    calls qh_findgood
    nop if good not used

  returns:
    clears facet->good if not good
    sets qh.num_good

  notes:
    this is like qh_findgood but more restrictive

  design:
    uses qh_findgood to mark good facets
    marks facets for qh.GOODvertex
    marks facets for qh.SPLITthreholds  
*/
void qh_findgood_all (facetT *facetlist) {
  facetT *facet, *bestfacet=NULL;
  realT angle, bestangle= REALmax;
  int  numgood=0, startgood;

  if (!qh GOODvertex && !qh GOODthreshold && !qh GOODpoint 
  && !qh SPLITthresholds)
    return;
  if (!qh ONLYgood)
    qh_findgood (qh facet_list, 0);
  FORALLfacet_(facetlist) {
    if (facet->good)
      numgood++;
  }
  if (qh GOODvertex <0 || (qh GOODvertex > 0 && qh MERGING)) {
    FORALLfacet_(facetlist) {
      if (facet->good && ((qh GOODvertex > 0) ^ !!qh_isvertex (qh GOODvertexp, facet->vertices))) {
        if (!--numgood) {
          fprintf (qh ferr, "qhull warning: good vertex p%d does not match last good facet f%d.  Ignored.\n",
             qh_pointid(qh GOODvertexp), facet->id);
          return;
        }
        facet->good= False;
      }
    }
  }
  startgood= numgood;
  if (qh SPLITthresholds) {
    FORALLfacet_(facetlist) {
      if (facet->good) {
        if (!qh_inthresholds (facet->normal, &angle)) {
          facet->good= False;
          numgood--;
          if (angle < bestangle) {
            bestangle= angle;
            bestfacet= facet;
          }
        }
      }
    }
    if (!numgood && bestfacet) {
      bestfacet->good= True;
      numgood++;
      trace0((qh ferr, "qh_findgood_all: f%d is closest (%2.2g) to thresholds\n", 
           bestfacet->id, bestangle));
      return;
    }
  }
  qh num_good= numgood;
  trace0((qh ferr, "qh_findgood_all: %d good facets remain out of %d facets\n",
        numgood, startgood));
} /* findgood_all */

/*-<a                             href="qh-c.htm"
  >-------------------------------</a><a name="furthestnext">-</a>
  
  qh_furthestnext()
    set qh.facet_next to facet with furthest of all furthest points
    searches all facets on qh.facet_list

  notes:
    this may help avoid precision problems
*/
void qh_furthestnext (void /* qh facet_list */) {
  facetT *facet, *bestfacet= NULL;
  realT dist, bestdist= -REALmax;

  FORALLfacets {
    if (facet->outsideset) {
#if qh_COMPUTEfurthest
      pointT *furthest;
      furthest= (pointT*)qh_setlast (facet->outsideset);
      zinc_(Zcomputefurthest);
      qh_distplane (furthest, facet, &dist);
#else
      dist= facet->furthestdist;
#endif
      if (dist > bestdist) {
	bestfacet= facet;
	bestdist= dist;
      }
    }
  }
  if (bestfacet) {
    qh_removefacet (bestfacet);
    qh_prependfacet (bestfacet, &qh facet_next);
    trace1((qh ferr, "qh_furthestnext: made f%d next facet (dist %.2g)\n",
	    bestfacet->id, bestdist));
  }
} /* furthestnext */

/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="furthestout">-</a>
  
  qh_furthestout( facet )
    make furthest outside point the last point of outsideset

  returns:
    updates facet->outsideset
    clears facet->notfurthest
    sets facet->furthestdist

  design:
    determine best point of outsideset
    make it the last point of outsideset
*/
void qh_furthestout (facetT *facet) {
  pointT *point, **pointp, *bestpoint= NULL;
  realT dist, bestdist= -REALmax;

  FOREACHpoint_(facet->outsideset) {
    qh_distplane (point, facet, &dist);
    zinc_(Zcomputefurthest);
    if (dist > bestdist) {
      bestpoint= point;
      bestdist= dist;
    }
  }
  if (bestpoint) {
    qh_setdel (facet->outsideset, point);
    qh_setappend (&facet->outsideset, point);
#if !qh_COMPUTEfurthest
    facet->furthestdist= bestdist;
#endif
  }
  facet->notfurthest= False;
  trace3((qh ferr, "qh_furthestout: p%d is furthest outside point of f%d\n",
	  qh_pointid (point), facet->id));
} /* furthestout */


/*-<a                             href="qh-c.htm#qhull"
  >-------------------------------</a><a name="infiniteloop">-</a>
  
  qh_infiniteloop( facet )
    report infinite loop error due to facet
*/
void qh_infiniteloop (facetT *facet) {

  fprintf (qh ferr, "qhull internal error (qh_infiniteloop): potential infinite loop detected\n");
  qh_errexit (qh_ERRqhull, facet, NULL);
} /* qh_infiniteloop */

/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="initbuild">-</a>
  
  qh_initbuild()
    initialize hull and outside sets with point array
    qh.FIRSTpoint/qh.NUMpoints is point array
    if qh.GOODpoint
      adds qh.GOODpoint to initial hull