poly2.c

这是一个新的用于图像处理的工具箱

/* poly2.c -- implements polygons and simplices

   see README, poly.h and qhull.h
   
   copyright (c) 1993-1995, The Geometry Center

   frequently used code is in poly.c
*/

#include "qhull_a.h"

/*======== functions in alphabetical order ==========*/



/*-----------------------------------------------
-addhash- add hash element to linear hash table if not already there
*/
void qh_addhash (void* newelem, setT *hashtable, int hashsize, unsigned hash) {
  unsigned scan;
  void *elem;

  for (scan= hash; (elem= SETelem_(hashtable, scan)); 
       scan= (++scan >= hashsize ? 0 : scan)) {
    if (elem == newelem)
      break;
  }
  /* loop terminates because qh_HASHfactor >= 1.1 by qh_initbuffers */
  if (!elem)
    SETelem_(hashtable, scan)= newelem;
} /* addhash */

/*-----------------------------------------------
-check_bestdist- check that points are within max_outside of the nearest facet
  if ONLYgood, ignores !good facets
  see: check_maxout
*/
void qh_check_bestdist (void) {
  boolT waserror= False, isoutside, unassigned;
  facetT *facet, *bestfacet, *errfacet1= NULL, *errfacet2= NULL;
  facetT *facetlist; 
  realT dist, maxoutside;
  pointT *point;
  int numpart, facet_i, facet_n, notgood= 0, notverified= 0;
  setT *facets;

  maxoutside= fmax_(qh max_outside, qh DISTround);
  maxoutside += 2 * qh DISTround;
  /* 1 DISTround to actual point and another DISTround to computed point */
  if (qh RANDOMdist) /* repeated computations can differ by 2*distround */
    maxoutside += qh DISTround;
  trace1((qh ferr, "qh_check_bestdist: check that all points are within %2.2g of best facet\n", maxoutside));
  facets= qh_pointfacet (/*qh facet_list*/);
  if (!qh_QUICKhelp && qh PRINTprecision)
    fprintf (qh ferr, "\n\
qhull output completed.  Verifying that %d points are\n\
below %2.2g of the nearest %sfacet.\n",
	     qh_setsize(facets), maxoutside, (qh ONLYgood ?  "good " : ""));
  FOREACHfacet_i_(facets) {  /* for each point with facet assignment */
    if (facet)
      unassigned= False;
    else {
      unassigned= True;
      facet= qh facet_list;
    }
    point= qh_point(facet_i);
    if (point == qh GOODpointp)
      continue;
    bestfacet= qh_findbest (point, facet, qh_ALL, False,
			    &dist, &isoutside, &numpart);
    /* occurs after statistics reported */
    if (dist > maxoutside) {
      if (qh ONLYgood && !bestfacet->good 
	  && !((bestfacet= qh_findgooddist (point, bestfacet, &dist, &facetlist))
	       && dist > maxoutside))
	notgood++;
      else {
	waserror= True;
	fprintf(qh ferr, "qhull precision error: point p%d is outside facet f%d, distance= %6.8g maxoutside= %6.8g\n", 
		facet_i, bestfacet->id, dist, maxoutside);
	errfacet2= errfacet1;
	errfacet1= bestfacet;		    
      }
    }else if (unassigned && dist < -qh MAXcoplanar)
      notverified++;
  }
  qh_settempfree (&facets);
  if (notverified && !qh DELAUNAY && !qh_QUICKhelp && qh PRINTprecision) 
    fprintf(qh ferr, "%d points were well inside the hull.  If the hull contains\n\
a lens-shaped component, these points were not verified.  Use\n\
options 'Qi Tv' to verify all points.\n", notverified); 
  if (waserror)
    qh_errexit2 (qh_ERRprec, errfacet1, errfacet2);
} /* check_bestdist */

/*-----------------------------------------------
-check_maxout- updates max_outside by checking all points against bestfacet
  updates facet->maxoutside via findbest
  if printing min_vertex, it is updated to the current vertices
  if ONLYgood, ignores !good facets
  see check_bestdist
notes:
  may not need to check near-inside points if KEEPcoplanar 
     (since coplanar is now MAXcoplanar instead of -DISTround)
*/
#ifndef qh_NOmerge
void qh_check_maxout (void) {
  facetT *facet, *bestfacet, *neighbor, **neighborp, *facetlist;
  realT dist, maxoutside, minvertex;
  pointT *point, **pointp;
  int numpart, facet_i, facet_n, notgood= 0;
  setT *facets, *vertices;
  vertexT *vertex;

  maxoutside= minvertex= 0;
  trace1((qh ferr, "qh_check_maxout: determine actual maxoutside and minoutside\n"));
  if (qh VERTEXneighbors 
  && (qh PRINTsummary || qh PRINTout[0] == qh_PRINTnone || qh PRINTstatistics
      || qh PRINTout[0] == qh_PRINTsummary || qh TRACElevel)) {
    vertices= qh_pointvertex (/*qh facet_list*/);
    FORALLvertices {
      FOREACHneighbor_(vertex) {
	zinc_(Zdistvertex);  /* distance also computed by main loop below */
	qh_distplane (vertex->point, neighbor, &dist);
	minimize_(minvertex, dist);
	if (-dist > qh TRACEdist || dist > qh TRACEdist 
	    || neighbor == qh tracefacet || vertex == qh tracevertex)
	  fprintf (qh ferr, "qh_check_maxout: p%d (v%d) is %.2g from f%d\n",
		   qh_pointid (vertex->point), vertex->id, dist, neighbor->id);
      }
    }
    if (qh MERGING) {
      wmin_(Wminvertex, qh min_vertex);
    }
    qh min_vertex= minvertex;
    qh_settempfree (&vertices);
  }
  facets= qh_pointfacet (/*qh facet_list*/);
  FOREACHfacet_i_(facets) {     /* for each point with facet assignment */
    if (facet) { 
      zinc_(Ztotcheck);
      point= qh_point(facet_i);
      if (point == qh GOODpointp)
	continue;
      bestfacet= qh_findbest (point, facet, qh_ALL,
			    False, &dist, NULL, &numpart);
      zadd_(Zcheckpart, numpart);
      if (bestfacet && dist > maxoutside) {
        if (qh ONLYgood && !bestfacet->good 
        && !((bestfacet= qh_findgooddist (point, bestfacet, &dist, &facetlist))
             && dist > maxoutside))
          notgood++;
        else
	  maxoutside= dist;
      }
      if (dist > qh TRACEdist || (bestfacet && bestfacet == qh tracefacet))
	fprintf (qh ferr, "qh_check_maxout: p%d is %.2g above f%d\n",
		   qh_pointid (point), dist, bestfacet->id);
    }
  }
  qh_settempfree (&facets);
  wval_(Wmaxout)= maxoutside - qh max_outside;
  wmax_(Wmaxoutside, qh max_outside);
  qh max_outside= maxoutside;
  if (qh KEEPnearinside && !qh KEEPcoplanar && !qh KEEPinside) {
    FORALLfacets {
      if (facet->coplanarset) 
        qh_setfree( &facet->coplanarset);
    }
  }else if (qh KEEPnearinside) {
    numpart= 0;
    FORALLfacets {   /* could be combined with qh_findbest */
      if (facet->coplanarset) {
        FOREACHpoint_(facet->coplanarset) {
  	  numpart++;
	  qh_distplane (point, facet, &dist); 
	  if (dist < -qh MAXcoplanar) {
	    if (!qh KEEPinside)
              SETref_(point)= NULL;
          }else if (!qh KEEPcoplanar)
            SETref_(point)= NULL;
        }
	qh_setcompact (facet->coplanarset);
      }
    }
    zadd_(Zcheckpart, numpart);
  }
  trace1((qh ferr, "qh_check_maxout: maxoutside %2.2g, minvertex %2.2g, outside of not good %d\n",
       maxoutside, minvertex, notgood));
} /* check_maxout */
#else /* qh_NOmerge */
void qh_check_maxout (void) {
}
#endif

/*----------------------------------------
-check_output- performs the checks at the end of qhull algorithm
  does not check points (may take a long time)
*/
void qh_check_output (void) {
  int i;

  if (qh STOPcone)
    return;
  if (qh VERIFYoutput | qh IStracing | qh CHECKfrequently) {
    qh_checkpolygon (qh facet_list);
    qh_checkflipped_all (qh facet_list);
    qh_checkconvex (qh facet_list, qh_ALGORITHMfault);
  }else if (!qh MERGING && qh_newstats (qhstat precision, &i)) {
    qh_checkflipped_all (qh facet_list);
    qh_checkconvex (qh facet_list, qh_ALGORITHMfault);
  }
} /* check_output */



/*-------------------------------------------------------------
-check_point- check that point is not outside facet
  if maxerror, doesn't report an error
*/
void qh_check_point (pointT *point, facetT *facet, realT *maxoutside, facetT **errfacet1, facetT **errfacet2) {
  realT dist;

  /* occurs after statistics reported */
  qh_distplane(point, facet, &dist);
  if (dist > *maxoutside) {
    *errfacet2= *errfacet1;
    *errfacet1= facet;
    fprintf(qh ferr, "qhull precision error: point p%d is outside facet f%d, distance= %6.8g maxoutside= %6.8g\n", 
	      qh_pointid(point), facet->id, dist, *maxoutside);
  }
} /* qh_check_point */


/*-------------------------------------------------
-check_points- checks that all points are inside all facets
     uses findbest if lots of points
     ignores flipped facets
notes:
  maxoutside includes 2 DISTrounds, one for the computed
  distances in qh_check_points
  qh_printafacet and qh_printsummary needs only one DISTround
*/
void qh_check_points (void) {
  facetT *facet, *errfacet1= NULL, *errfacet2= NULL;
  realT total, maxoutside;
  pointT *point, **pointp, *pointtemp;
  boolT testouter;

  maxoutside= fmax_(qh max_outside, qh DISTround);  /* agrees with qh_printafacet, qh_printsummary */
  maxoutside += 2* qh DISTround;
  /* 1 DISTround to actual point and another DISTround to computed point */
  if (qh RANDOMdist) /* repeated computations can differ by 2*distround */
    maxoutside += qh DISTround;
  trace1((qh ferr, "qh_check_points: check all points below %2.2g of all facet planes\n",
	  maxoutside));
  if (qh num_good)   /* miss counts other_points and !good facets */
     total= (float) qh num_good * qh num_points;
  else
     total= (float) qh num_facets * qh num_points;
  if (total >= qh_VERIFYdirect 
  && (!qh MERGING || qh SKIPcheckmax || qh ZEROall_ok)) {
    if (!qh_QUICKhelp && qh SKIPcheckmax && qh MERGING)
      fprintf (qh ferr, "\n\
qhull input warning: merging without checking outer planes ('Q5').\n\
Verify may report that a point is outside of a facet.\n");
    qh_check_bestdist();
  }else {
    if (qh MERGING && qh_MAXoutside && !qh SKIPcheckmax)
      testouter= True;  /* agrees with qh_printafacet() */
    else
      testouter= False;
    if (!qh_QUICKhelp) {
      if (qh MERGEexact || qh SKIPcheckmax || qh NOnearinside)
	fprintf (qh ferr, "\n\
qhull input warning: exact merge ('Qx'), no outer plane check ('Q5'), or\n\
no processing of near-inside points ('Q8').  Verify may report that a point\n\
is outside of a facet.\n");
    }
    if (qh PRINTprecision) {
      if (testouter)
	fprintf (qh ferr, "\n\
Output completed.  Verifying that all points are below outer planes of\n\
all %sfacets.  Will make %2.0f distance computations.\n", 
	      (qh ONLYgood ?  "good " : ""), total);
      else
	fprintf (qh ferr, "\n\
Output completed.  Verifying that all points are below %2.2g of\n\
all %sfacets.  Will make %2.0f distance computations.\n", 
	      maxoutside, (qh ONLYgood ?  "good " : ""), total);
    }
    FORALLfacets {
      if (!facet->good && qh ONLYgood)
        continue;
      if (facet->flipped)
        continue;
      if (testouter) {
#if qh_MAXoutside
	maxoutside= facet->maxoutside + 2* qh DISTround;
	/* 1 DISTround to actual point and another to computed point */
#endif
      }
      FORALLpoints {
	if (point != qh GOODpointp)
	  qh_check_point (point, facet, &maxoutside, &errfacet1, &errfacet2);
      }
      FOREACHpoint_(qh other_points) {
	if (point != qh GOODpointp)
	  qh_check_point (point, facet, &maxoutside, &errfacet1, &errfacet2);
      }
    }
    if (errfacet1)
      qh_errexit2(qh_ERRprec, errfacet1, errfacet2);
  }
} /* check_points */


/*-------------------------------------------------
-checkconvex- check that each ridge in facetlist is convex
returns:
    counts Zconcaveridges and Zcoplanarridges
    errors if concaveridge or if merging an coplanar ridge
note:
    if not merging, tests vertices for neighboring simplicial facets
    else if ZEROcentrum, tests vertices for neighboring simplicial facets
    else tests centrums of neighboring facets
*/
void qh_checkconvex(facetT *facetlist, int fault) {
  facetT *facet, *neighbor, **neighborp, *errfacet1=NULL, *errfacet2=NULL;
  vertexT *vertex;
  realT dist;
  pointT *centrum;
  boolT waserror= False, tempcentrum= False, allsimplicial;
  int neighbor_i;

  trace1((qh ferr, "qh_checkconvex: check all ridges are convex\n"));
  zzval_(Zconcaveridges)= 0;
  zzval_(Zcoplanarridges)= 0;
  FORALLfacet_(facetlist) {
    if (facet->flipped) {
      fprintf (qh ferr, "qhull precision error: f%d is flipped (interior point is outside)\n",
	       facet->id);
      errfacet1= facet;
      waserror= True;
      continue;
    }
    if (qh MERGING && (!qh ZEROcentrum || !facet->simplicial))
      allsimplicial= False;
    else {
      allsimplicial= True;
      neighbor_i= 0;
      FOREACHneighbor_(facet) {
        vertex= SETelem_(facet->vertices, neighbor_i++);
	if (!neighbor->simplicial) {
	  allsimplicial= False;
	  continue;
	}
        qh_distplane (vertex->point, neighbor, &dist);
        if (dist > -qh DISTround) {
	  if (fault == qh_DATAfault) {
	    fprintf (qh ferr, "qhull precision error: initial simplex is not convex. Distance=%.2g\n", dist);
	    qh_errexit(qh_ERRsingular, NULL, NULL);
	  }
          if (dist > qh DISTround) {
            zzinc_(Zconcaveridges);
            fprintf (qh ferr, "qhull precision error: f%d is concave to f%d, since p%d (v%d) is %6.4g above\n",
              facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist);
            errfacet1= facet;
            errfacet2= neighbor;
            waserror= True;
          }else if (qh ZEROcentrum) {
            if (dist > 0) {     /* qh_checkzero checks that dist < - qh DISTround */
              zzinc_(Zcoplanarridges); 
              fprintf (qh ferr, "qhull precision error: f%d is clearly not convex to f%d, since p%d (v%d) is %6.4g above\n",
                facet->id, neighbor->id, qh_pointid(vertex->point), vertex->id, dist);