qhull.c

关于网格剖分的

      dist= facet->furthestdist;
#endif
      if (dist < qh MINoutside) { /* remainder of outside set is coplanar for qh_outcoplanar */
	qh facet_next= facet->next;
	continue;
      }
    }
    if (!qh RANDOMoutside && !qh VIRTUALmemory) {
      if (qh PICKfurthest) {
	qh_furthestnext (/* qh facet_list */);
	facet= qh facet_next;
      }
      *visible= facet;
      return ((pointT*)qh_setdellast (facet->outsideset));
    }
    if (qh RANDOMoutside) {
      int outcoplanar = 0;
      if (qh NARROWhull) {
        FORALLfacets {
	  if (facet == qh facet_next)
	    break;
	  if (facet->outsideset)
  	    outcoplanar += qh_setsize( facet->outsideset);
	}
      }
      randr= qh_RANDOMint;
      randr= randr/(qh_RANDOMmax+1);
      index= (int)floor((qh num_outside - outcoplanar) * randr);
      FORALLfacet_(qh facet_next) {
        if (facet->outsideset) {
          SETreturnsize_(facet->outsideset, size);
          if (!size)
            qh_setfree (&facet->outsideset);
          else if (size > index) {
            *visible= facet;
            return ((pointT*)qh_setdelnth (facet->outsideset, index));
          }else
            index -= size;
        }
      }
      fprintf (qh ferr, "qhull internal error (qh_nextfurthest): num_outside %d is too low\nby at least %d, or a random real %g >= 1.0\n",
              qh num_outside, index+1, randr);
      qh_errexit (qh_ERRqhull, NULL, NULL);
    }else { /* VIRTUALmemory */
      facet= qh facet_tail->previous;
      if (!(furthest= (pointT*)qh_setdellast(facet->outsideset))) {
        if (facet->outsideset)
          qh_setfree (&facet->outsideset);
        qh_removefacet (facet);
        qh_prependfacet (facet, &qh facet_list);
        continue;
      }
      *visible= facet;
      return furthest;
    }
  }
  return NULL;
} /* nextfurthest */

/*-<a                             href="qh-c.htm#qhull"
  >-------------------------------</a><a name="partitionall">-</a>
  
  qh_partitionall( vertices, points, numpoints )
    partitions all points in points/numpoints to the outsidesets of facets
    vertices= vertices in qh.facet_list (not partitioned)

  returns:
    builds facet->outsideset
    does not partition qh.GOODpoint
    if qh.ONLYgood && !qh.MERGING, 
      does not partition qh.GOODvertex
    sets facet->newfacet for qh_findbestnew() in qh_partitionpoint()

  notes:
    faster if qh.facet_list sorted by anticipated size of outside set

  design:
    initialize pointset with all points
    remove vertices from pointset
    remove qh.GOODpointp from pointset (unless it's qh.STOPcone or qh.STOPpoint)
    for all facets
      for all remaining points in pointset
        compute distance from point to facet
        if point is outside facet
          remove point from pointset (by not reappending)
          update bestpoint
          append point or old bestpoint to facet's outside set
      append bestpoint to facet's outside set (furthest)
    for all points remaining in pointset
      partition point into facets' outside sets and coplanar sets
*/
void qh_partitionall(setT *vertices, pointT *points, int numpoints){
  setT *pointset;
  vertexT *vertex, **vertexp;
  pointT *point, **pointp, *bestpoint;
  int size, point_i, point_n, point_end, remaining, i, id;
  facetT *facet;
  realT bestdist= -REALmax, dist, distoutside;
    
  trace1((qh ferr, "qh_partitionall: partition all points into outside sets\n"));
  pointset= qh_settemp (numpoints);
  qh num_outside= 0;
  pointp= SETaddr_(pointset, pointT);
  for (i=numpoints, point= points; i--; point += qh hull_dim)
    *(pointp++)= point;
  qh_settruncate (pointset, numpoints);
  FOREACHvertex_(vertices) {
    if ((id= qh_pointid(vertex->point)) >= 0)
      SETelem_(pointset, id)= NULL;
  }
  id= qh_pointid (qh GOODpointp);
  if (id >=0 && qh STOPcone-1 != id && -qh STOPpoint-1 != id)
    SETelem_(pointset, id)= NULL;
  if (qh GOODvertexp && qh ONLYgood && !qh MERGING) { /* matches qhull()*/
    if ((id= qh_pointid(qh GOODvertexp)) >= 0)
      SETelem_(pointset, id)= NULL;
  }
  if (!qh BESToutside) {  /* matches conditional for qh_partitionpoint below */
    if (qh MERGING)
      distoutside= qh_DISToutside; /* defined in user.h */
    else
      distoutside= qh MINoutside;
    zval_(Ztotpartition)= qh num_points - qh hull_dim - 1; /*misses GOOD... */
    remaining= qh num_facets;
    point_end= numpoints;
    FORALLfacets {
      size= point_end/(remaining--) + 100;
      facet->outsideset= qh_setnew (size);
      bestpoint= NULL;
      point_end= 0;
      FOREACHpoint_i_(pointset) {
        if (point) {
          zzinc_(Zpartitionall);
          qh_distplane (point, facet, &dist);
          if (dist < distoutside)
            SETelem_(pointset, point_end++)= point;
          else {
	    qh num_outside++;
            if (!bestpoint) {
              bestpoint= point;
              bestdist= dist;
            }else if (dist > bestdist) {
              qh_setappend (&facet->outsideset, bestpoint);
              bestpoint= point;
              bestdist= dist;
            }else 
              qh_setappend (&facet->outsideset, point);
          }
        }
      }
      if (bestpoint) {
        qh_setappend (&facet->outsideset, bestpoint);
#if !qh_COMPUTEfurthest
	facet->furthestdist= bestdist;
#endif
      }else
        qh_setfree (&facet->outsideset);
      qh_settruncate (pointset, point_end);
    }
  }
  if (qh BESToutside || qh MERGING || qh KEEPcoplanar || qh KEEPinside) {
    qh findbestnew= True;
    FOREACHpoint_i_(pointset) { 
      if (point)
        qh_partitionpoint(point, qh facet_list);
    }
    qh findbestnew= False;
  }
  zzadd_(Zpartitionall, zzval_(Zpartition));
  zzval_(Zpartition)= 0;
  qh_settempfree(&pointset);
  if (qh IStracing >= 4)
    qh_printfacetlist (qh facet_list, NULL, True);
} /* partitionall */


/*-<a                             href="qh-c.htm#qhull"
  >-------------------------------</a><a name="partitioncoplanar">-</a>
  
  qh_partitioncoplanar( point, facet, dist )
    partition coplanar point to a facet
    dist is distance from point to facet
    if dist NULL, 
      searches for bestfacet and does nothing if inside
    if qh.findbestnew set, 
      searches new facets instead of using qh_findbest()

  returns:
    qh.max_ouside updated
    if qh.KEEPcoplanar or qh.KEEPinside
      point assigned to best coplanarset
  
  notes:
    facet->maxoutside is updated at end by qh_check_maxout

  design:
    if dist undefined
      find best facet for point
      if point sufficiently below facet (depends on qh.NEARinside and qh.KEEPinside)
        exit
    if keeping coplanar/nearinside/inside points
      if point is above furthest coplanar point
        append point to coplanar set (it is the new furthest)
        update qh.max_outside
      else
        append point one before end of coplanar set
    else
      update qh.max_outside
*/
void qh_partitioncoplanar (pointT *point, facetT *facet, realT *dist) {
  facetT *bestfacet;
  pointT *oldfurthest;
  realT bestdist, dist2;
  int numpart= 0;
  boolT isoutside, istrace= False;

  qh WAScoplanar= True;
  if (!dist) {
    if (qh findbestnew)
      bestfacet= qh_findbestnew (point, facet, 
			  &bestdist, NULL, &numpart);
    else
      bestfacet= qh_findbest (point, facet, qh_ALL, False, !qh_NOupper, 
                          &bestdist, &isoutside, &numpart);
    zinc_(Ztotpartcoplanar);
    zzadd_(Zpartcoplanar, numpart);
    if (!qh KEEPinside) {
      if (qh KEEPnearinside) {
        if (bestdist < -qh NEARinside) { 
          zinc_(Zcoplanarinside);
          return;
        }
      }else if (bestdist < -qh MAXcoplanar) {
        zinc_(Zcoplanarinside);
        return;
      }
    }
  }else {
    bestfacet= facet;
    bestdist= *dist;
  }
  if (qh KEEPcoplanar + qh KEEPinside + qh KEEPnearinside) {
    oldfurthest= (pointT*)qh_setlast (bestfacet->coplanarset);
    if (oldfurthest) {
      zinc_(Zcomputefurthest);
      qh_distplane (oldfurthest, bestfacet, &dist2);
    }
    if (!oldfurthest || dist2 < bestdist) {
      qh_setappend(&bestfacet->coplanarset, point);
      if (bestdist > qh max_outside) {
	qh max_outside= bestdist;
	if (bestdist > qh TRACEdist)
	  istrace= True;
      }
    }else
      qh_setappend2ndlast(&bestfacet->coplanarset, point);
  }else { /* !KEEPcoplanar && !KEEPinside */
    if (bestdist > qh max_outside) {
      qh max_outside= bestdist;
      if (bestdist > qh TRACEdist) 
	istrace= True;
    }
  }
  if (istrace) {
    fprintf (qh ferr, "qh_partitioncoplanar: ====== p%d increases max_outside to %2.2g of f%d last p%d\n",
		   qh_pointid(point), bestdist, bestfacet->id, qh furthest_id);
    qh_errprint ("DISTANT", bestfacet, NULL, NULL, NULL);
  }
  trace4((qh ferr, "qh_partitioncoplanar: point p%d is coplanar with facet f%d (or inside) dist %2.2g\n",
	  qh_pointid(point), bestfacet->id, bestdist));
} /* partitioncoplanar */

/*-<a                             href="qh-c.htm#qhull"
  >-------------------------------</a><a name="partitionpoint">-</a>
  
  qh_partitionpoint( point, facet )
    assigns point to an outside set, coplanar set, or inside set (i.e., dropt)
    if qh.findbestnew
      uses qh_findbestnew() to search all new facets
    else
      uses qh_findbest()
  
  notes:
    after qh_distplane(), this and qh_findbest() are most expensive in 3-d

  design:
    find best facet for point 
      (either exhaustive search of new facets or directed search from facet)
    if qh.NARROWhull
      retain coplanar and nearinside points as outside points
    if point is outside bestfacet
      if point above furthest point for bestfacet
        append point to outside set (it becomes the new furthest)
        if outside set was empty
          move bestfacet to end of qh.facet_list (i.e., after qh.facet_next)
        update bestfacet->furthestdist
      else
        append point one before end of outside set
    else if point is coplanar to bestfacet
      if keeping coplanar points or need to update qh.max_outside
        partition coplanar point into bestfacet
    else if near-inside point        
      partition as coplanar point into bestfacet
    else is an inside point
      if keeping inside points 
        partition as coplanar point into bestfacet
*/
void qh_partitionpoint (pointT *point, facetT *facet) {
  realT bestdist;
  boolT isoutside;
  facetT *bestfacet;
  int numpart;
#if qh_COMPUTEfurthest
  realT dist;
#endif

  if (qh findbestnew)
    bestfacet= qh_findbestnew (point, facet,
			  &bestdist, &isoutside, &numpart);
  else
    bestfacet= qh_findbest (point, facet, qh BESToutside, True, !qh_NOupper,
			  &bestdist, &isoutside, &numpart);
  zinc_(Ztotpartition);
  zzadd_(Zpartition, numpart);
  if (qh NARROWhull) {
    if (qh DELAUNAY && !isoutside && bestdist >= -qh MAXcoplanar)
      qh_precision ("nearly incident point (narrow hull)");
    if (qh KEEPnearinside) {
      if (bestdist >= -qh NEARinside)
	isoutside= True;
    }else if (bestdist >= -qh MAXcoplanar)