poly2.c

Quick hull implementation

  design:
    if simplicial facet
      build set from facet->vertices with facet->toporient
    else
      for each ridge in order
        build set from ridge's vertices
*/
setT *qh_facet3vertex (facetT *facet) {
  ridgeT *ridge, *firstridge;
  vertexT *vertex;
  int cntvertices, cntprojected=0;
  setT *vertices;

  cntvertices= qh_setsize(facet->vertices);
  vertices= qh_settemp (cntvertices);
  if (facet->simplicial) {
    if (cntvertices != 3) {
      fprintf (qh ferr, "qhull internal error (qh_facet3vertex): only %d vertices for simplicial facet f%d\n", 
                  cntvertices, facet->id);
      qh_errexit(qh_ERRqhull, facet, NULL);
    }
    qh_setappend (&vertices, SETfirst_(facet->vertices));
    if (facet->toporient ^ qh_ORIENTclock)
      qh_setappend (&vertices, SETsecond_(facet->vertices));
    else
      qh_setaddnth (&vertices, 0, SETsecond_(facet->vertices));
    qh_setappend (&vertices, SETelem_(facet->vertices, 2));
  }else {
    ridge= firstridge= SETfirstt_(facet->ridges, ridgeT);   /* no infinite */
    while ((ridge= qh_nextridge3d (ridge, facet, &vertex))) {
      qh_setappend (&vertices, vertex);
      if (++cntprojected > cntvertices || ridge == firstridge)
        break;
    }
    if (!ridge || cntprojected != cntvertices) {
      fprintf (qh ferr, "qhull internal error (qh_facet3vertex): ridges for facet %d don't match up.  got at least %d\n", 
                  facet->id, cntprojected);
      qh_errexit(qh_ERRqhull, facet, ridge);
    }
  }
  return vertices;
} /* facet3vertex */

/*-<a                             href="qh-poly.htm#TOC"
  >-------------------------------</a><a name="findbestfacet">-</a>
  
  qh_findbestfacet( point, bestoutside, bestdist, isoutside )
    find facet that is furthest 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:
    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:
    For tricoplanar facets, this finds one of the tricoplanar facets closest 
    to the point.  For Delaunay triangulations, the point may be inside a 
    different tricoplanar facet. See <a href="../html/qh-in.htm#findfacet">locate a facet with qh_findbestfacet()</a>
    
    If inside, qh_findbestfacet performs an exhaustive search
       this may be too conservative.  Sometimes it is clearly required.

    qh_findbestfacet is not used by qhull.
    uses qh.visit_id and qh.coplanarset
    
  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, !qh_ISnewfacets, bestoutside /* qh_NOupper */,
			    bestdist, isoutside, &totpart);
  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-poly.htm#TOC"
  >-------------------------------</a><a name="findbestlower">-</a>
  
  qh_findbestlower( facet, point, bestdist, numpart )
    returns best non-upper, non-flipped neighbor of facet for point
    if needed, searches vertex neighbors 

  returns:
    returns bestdist and updates numpart

  notes:
    if Delaunay and inside, point is outside of circumsphere of bestfacet
    called by qh_findbest() for points above an upperdelaunay facet

*/
facetT *qh_findbestlower (facetT *upperfacet, pointT *point, realT *bestdistp, int *numpart) {
  facetT *neighbor, **neighborp, *bestfacet= NULL;
  realT bestdist= -REALmax/2 /* avoid underflow */;
  realT dist;
  vertexT *vertex;

  zinc_(Zbestlower);
  FOREACHneighbor_(upperfacet) {
    if (neighbor->upperdelaunay || neighbor->flipped)
      continue;
    (*numpart)++;
    qh_distplane (point, neighbor, &dist);
    if (dist > bestdist) {
      bestfacet= neighbor;
      bestdist= dist;
    }
  }
  if (!bestfacet) {
    zinc_(Zbestlowerv);
    /* rarely called, numpart does not count nearvertex computations */
    vertex= qh_nearvertex (upperfacet, point, &dist);
    qh_vertexneighbors();
    FOREACHneighbor_(vertex) {
      if (neighbor->upperdelaunay || neighbor->flipped)
	continue;
      (*numpart)++;
      qh_distplane (point, neighbor, &dist);
      if (dist > bestdist) {
	bestfacet= neighbor;
	bestdist= dist;
      }
    }
  }
  if (!bestfacet) {
    fprintf(qh ferr, "\n\
qh_findbestlower: all neighbors of facet %d are flipped or upper Delaunay.\n\
Please report this error to [e-mail] with the input and all of the output.\n",
       upperfacet->id);
    qh_errexit (qh_ERRqhull, upperfacet, NULL);
  }
  *bestdistp= bestdist;
  trace3((qh ferr, "qh_findbestlower: f%d dist %2.2g for f%d p%d\n",
	  bestfacet->id, bestdist, upperfacet->id, qh_pointid(point)));
  return bestfacet;
} /* findbestlower */

/*-<a                             href="qh-poly.htm#TOC"
  >-------------------------------</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-poly.htm#TOC"
  >-------------------------------</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-poly.htm#TOC"
  >-------------------------------</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 and ONLYgood, 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) {
	  if (qh ONLYgood) {
            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;
	  }else if (qh GOODvertex > 0)
            fprintf (qh ferr, "qhull warning: point p%d is not a vertex ('QV%d').\n",
		qh GOODvertex-1, qh GOODvertex-1);
	  else
            fprintf (qh ferr, "qhull warning: point p%d is a vertex for every facet ('QV-%d').\n",
	        -qh GOODvertex - 1, -qh GOODvertex - 1);
        }