poly2.c

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

      return;
    }
  }
  qh num_good= numgood;
  trace0((qh ferr, "qh_findgood_all: %d good facets remain out of %d facets\n",
        numgood, startgood));
} /* findgood_all */

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

/*--------------------------------------------------
-initialhull- constructs the initial hull as a qh hull_dim simplex of vertices
*/
void qh_initialhull(setT *vertices) {
  facetT *facet, *firstfacet;
  realT dist;
#ifndef qh_NOtrace
  int k;
#endif

  qh_createsimplex(vertices);  /* qh facet_list */
  qh interior_point= qh_getcenter(vertices);
  firstfacet= qh facet_list;
  qh_setfacetplane(firstfacet);
  zinc_(Znumvisibility); /* needs to be in printsummary */
  qh_distplane(qh interior_point, firstfacet, &dist);
  if (dist > 0) {  
    FORALLfacets
      facet->toporient ^= True;
  }
  FORALLfacets
    qh_setfacetplane(facet);
  FORALLfacets {
    if (!qh_checkflipped (facet, NULL, qh_ALL)) {/* due to axis-parallel facet */
      trace1((qh ferr, "qh_initialhull: initial orientation incorrect.  Correct all facets\n"));
      facet->flipped= False;
      FORALLfacets {
	facet->toporient ^= True;
	qh_orientoutside (facet);
      }
      break;
    }
  }
  FORALLfacets {
    if (!qh_checkflipped (facet, NULL, !qh_ALL)) {  /* can happen with 'R0.1' */
      fprintf (qh ferr, "qhull precision error: facet %d is coplanar with the interior point\n",
                   facet->id);
      qh_errexit (qh_ERRsingular, facet, NULL);
    }
  }
  zzval_(Zprocessed)= qh hull_dim+1;
  qh_checkpolygon (qh facet_list);
  qh_checkconvex(qh facet_list,   qh_DATAfault);
#ifndef qh_NOtrace
  if (qh IStracing >= 1) {
    fprintf(qh ferr, "qh_initialhull: simplex constructed, interior point:");
    for (k=0; k<qh hull_dim; k++) 
      fprintf (qh ferr, " %6.4g", qh interior_point[k]);
    fprintf (qh ferr, "\n");
  }
#endif
} /* initialhull */

/*-------------------------------------------------
-initialvertices- determines a non-singular set of initial vertices
  maxpoints are not unique
returns:
  temporary set of dim+1 vertices in descending order by vertex id
notes:
  unless qh ALLpoints, uses maxpoints as long as determinate is non-zero
  picks random points if qh RANDOMoutside && !ALLpoints
  if dim >= qh_INITIALmax, uses min/max x and max points with non-zero determinants
*/
setT *qh_initialvertices(int dim, setT *maxpoints, pointT *points, int numpoints) {
  pointT *point, **pointp;
  setT *vertices, *simplex, *tested;
  realT randr, det;
  int index, point_i, point_n, k;
  boolT nearzero= False;
  
  vertices= qh_settemp (dim + 1);
  simplex= qh_settemp (dim+1);
  if (qh ALLpoints) 
    qh_maxsimplex (dim, NULL, points, numpoints, &simplex);
  else if (qh RANDOMoutside) {
    while (qh_setsize (simplex) != dim+1) {
      randr= qh_RANDOMint;
      randr= randr/(qh_RANDOMmax+1);
      index= floor(qh num_points * randr);
      point= qh_point (index);
      qh_setunique (&simplex, point);
    }
  }else if (qh hull_dim >= qh_INITIALmax) {
    tested= qh_settemp (dim+1);
    qh_setappend (&simplex, SETfirst_(maxpoints));   /* max and min X coord */
    qh_setappend (&simplex, SETsecond_(maxpoints));
    qh_maxsimplex (fmin_(qh_INITIALsearch, dim), maxpoints, points, numpoints, &simplex);
    k= qh_setsize (simplex);
    FOREACHpoint_i_(maxpoints) { 
      if (point_i & 0x1) {     /* first pick up max. coord. points */
      	if (!qh_setin (simplex, point) && !qh_setin (tested, point)){
	  det= qh_detsimplex(point, simplex, k, &nearzero);
          if (nearzero)
            qh_setappend (&tested, point);
          else {
            qh_setappend (&simplex, point);
            if (++k == dim)  /* use search for last point */
	      break;
	  }
	}
      }
    }
    while (k != dim && (point= (pointT*)qh_setdellast (maxpoints))) {
      if (!qh_setin (simplex, point) && !qh_setin (tested, point)){
        det= qh_detsimplex (point, simplex, k, &nearzero);
        if (nearzero)
          qh_setappend (&tested, point);
        else {
          qh_setappend (&simplex, point);
          k++;
	}
      }
    }
    index= 0;
    while (k != dim && (point= qh_point (index++))) {
      if (!qh_setin (simplex, point) && !qh_setin (tested, point)){
        det= qh_detsimplex (point, simplex, k, &nearzero);
        if (!nearzero){
          qh_setappend (&simplex, point);
          k++;
	}
      }
    }
    qh_settempfree (&tested);
    qh_maxsimplex (dim, maxpoints, points, numpoints, &simplex);
  }else
    qh_maxsimplex (dim, maxpoints, points, numpoints, &simplex);
  FOREACHpoint_(simplex) 
    qh_setaddnth (&vertices, 0, qh_newvertex(point)); /* descending order */
  qh_settempfree (&simplex);
  return vertices;
} /* initialvertices */


/*-------------------------------------------------
-isvertex- returns vertex if point is in vertex set, else returns NULL
*/
vertexT *qh_isvertex (pointT *point, setT *vertices) {
  vertexT *vertex, **vertexp;

  FOREACHvertex_(vertices) {
    if (vertex->point == point)
      return vertex;
  }
  return NULL;
} /* isvertex */

/*-------------------------------------------------
-makenewfacets- make new facets from point and qh visible_list
returns:
  qh newfacet_list= list of new facets with hyperplanes and ->newfacet
  qh newvertex_list= list of vertices in new facets with ->newlist set
  if (qh ONLYgood)
    newfacets reference horizon facets, but not vice versa
    ridges reference non-simplicial horizon ridges, but not vice versa
    does not change existing facets
  otherwise
    NEWfacets set
    newfacets attached to horizon facets and ridges
    visible->r.replace is corresponding new facet
*/
vertexT *qh_makenewfacets (pointT *point /*visible_list*/) {
  facetT *visible, *newfacet= NULL, *newfacet2= NULL, *neighbor, **neighborp;
  vertexT *apex;
  int numnew=0;

  qh newfacet_list= qh facet_tail;
  qh newvertex_list= qh vertex_tail;
  apex= qh_newvertex(point);
  qh_appendvertex (apex);  
  qh visit_id++;
  if (!qh ONLYgood)
    qh NEWfacets= True;
  FORALLvisible_facets {
    FOREACHneighbor_(visible) 
      neighbor->seen= False;
    if (visible->ridges) {
      visible->visitid= qh visit_id;
      newfacet2= qh_makenew_nonsimplicial (visible, apex, &numnew);
    }
    if (visible->simplicial)
      newfacet= qh_makenew_simplicial (visible, apex, &numnew);
    if (!qh ONLYgood) {
      if (newfacet2)  /* newfacet is null if all ridges defined */
        newfacet= newfacet2;
      if (newfacet)
      	visible->f.replace= newfacet;
      else
        zinc_(Zinsidevisible);
      SETfirst_(visible->neighbors)= NULL;
    }
  }
  trace1((qh ferr, "qh_makenewfacets: created %d new facets from point p%d to horizon\n",
	  numnew, qh_pointid(point)));
  if (qh IStracing >= 4)
    qh_printfacetlist (qh newfacet_list, NULL, qh_ALL);
  return apex;
} /* makenewfacets */

/*-------------------------------------------------
-matchduplicates- match duplicate ridges in hashtable
  marked with ->dupridge and qh_DUPLICATEridge
  picks match with worst merge (min distance apart)
returns:
  updates hashcount
  similar to qh_matchneighbor
*/
#ifndef qh_NOmerge
void qh_matchduplicates (facetT *atfacet, int atskip, int hashsize, int *hashcount) {
  boolT same, ismatch;
  unsigned hash, scan;
  facetT *facet, *newfacet, *maxmatch= NULL, *maxmatch2= NULL, *nextfacet;
  int skip, newskip, nextskip= 0, maxskip= 0, maxskip2= 0, makematch;
  realT maxdist= -REALmax, mindist, dist2, low, high;

  hash= qh_gethash (hashsize, atfacet->vertices, qh hull_dim, 1, 
                     SETelem_(atfacet->vertices, atskip));
  trace2((qh ferr, "qh_matchduplicates: find duplicate matches for f%d skip %d hash %d hashcount %d\n",
	  atfacet->id, atskip, hash, *hashcount));
  for (makematch= 0; makematch < 2; makematch++) {
    qh visit_id++;
    for (newfacet= atfacet, newskip= atskip; newfacet; newfacet= nextfacet, newskip= nextskip) {
      zinc_(Zhashlookup);
      nextfacet= NULL;
      newfacet->visitid= qh visit_id;
      for (scan= hash; (facet= SETelem_(qh hash_table, scan)); 
	   scan= (++scan >= hashsize ? 0 : scan)) {
	if (!facet->dupridge || facet->visitid == qh visit_id)
	  continue;
	zinc_(Zhashtests);
	if (qh_matchvertices (1, newfacet->vertices, newskip, facet->vertices, &skip, &same)) {
	  ismatch= (same == (newfacet->toporient ^ facet->toporient));
	  if (SETelem_(facet->neighbors, skip) != qh_DUPLICATEridge) {
	    if (!makematch) {
	      fprintf (qh ferr, "qhull internal error (qh_matchduplicates): missing dupridge at f%d skip %d for new f%d skip %d hash %d\n",
		     facet->id, skip, newfacet->id, newskip, hash);
	      qh_errexit2 (qh_ERRqhull, facet, newfacet);
	    }
	  }else if (ismatch && makematch) {
	    if (SETelem_(newfacet->neighbors, newskip) == qh_DUPLICATEridge) {
	      SETelem_(facet->neighbors, skip)= newfacet;
	      SETelem_(newfacet->neighbors, newskip)= qh_MERGEridge;
	      *hashcount -= 2; /* removed two unmatched facets */
	      trace4((qh ferr, "qh_matchduplicates: duplicate f%d skip %d matched with new f%d skip %d merge\n",
		    facet->id, skip, newfacet->id, newskip));
	    }
	  }else if (ismatch) {
	    mindist= qh_getdistance (facet, newfacet, &low, &high);
	    dist2= qh_getdistance (newfacet, facet, &low, &high);
	    minimize_(mindist, dist2);
	    if (mindist > maxdist) {
	      maxdist= mindist;
	      maxmatch= facet;
	      maxskip= skip;
	      maxmatch2= newfacet;
	      maxskip2= newskip;
	    }
	    trace3((qh ferr, "qh_matchduplicates: duplicate f%d skip %d new f%d skip %d at dist %2.2g, max is now f%d f%d\n",
		    facet->id, skip, newfacet->id, newskip, mindist, 
		    maxmatch->id, maxmatch2->id));
	  }else { /* !ismatch */
	    nextfacet= facet;
	    nextskip= skip;
	  }
	}
	if (makematch && !facet 
        && SETelem_(facet->neighbors, skip) == qh_DUPLICATEridge) {
	  fprintf (qh ferr, "qhull internal error (qh_matchduplicates): no MERGEridge match for duplicate f%d skip %d at hash %d\n",
		     newfacet->id, newskip, hash);
	  qh_errexit (qh_ERRqhull, newfacet, NULL);
	}
      }
    } /* end of for each new facet at hash */
    if (!makematch) {
      if (!maxmatch) {
	fprintf (qh ferr, "qhull internal error (qh_matchduplicates): no maximum match at duplicate f%d skip %d at hash %d\n",
		     atfacet->id, atskip, hash);
	qh_errexit (qh_ERRqhull, atfacet, NULL);
      }
      SETelem_(maxmatch->neighbors, maxskip)= maxmatch2;
      SETelem_(maxmatch2->neighbors, maxskip2)= maxmatch;
      *hashcount -= 2; /* removed two unmatched facets */
      zzinc_(Zmultiridge);
      trace1((qh ferr, "qh_matchduplicates: duplicate f%d skip %d matched with new f%d skip %d keep\n",
	      maxmatch->id, maxskip, maxmatch2->id, maxskip2));
      if (qh IStracing >= 4)
	qh_errprint ("DUPLICATED/MATCH", maxmatch, maxmatch2, NULL, NULL);
    }
  }
} /* matchduplicates */

/*----------------------------------------
-nearvertex- return nearest vertex to point
returns:
  vertex and distance
*/
vertexT *qh_nearvertex (facetT *facet, pointT *point, realT *bestdistp) {
  realT bestdist= REALmax, dist;
  vertexT *bestvertex= NULL, *vertex, **vertexp;

  FOREACHvertex_(facet->vertices) {
    dist= qh_pointdist (vertex->point, point, -qh hull_dim);
    if (dist < bestdist) {
      bestdist= dist;
      bestvertex= vertex;
    }
  }
  *bestdistp= sqrt (bestdist);
  return bestvertex;
} /* nearvertex */

/*-------------------------------------------------
-newhashtable- returns size of qh hash_table of at least newsize slots
  assumes qh hash_table is NULL
  qh_HASHfactor determines the number of extra slots
*/
int qh_newhashtable(int newsize) {
  int size;

  size= ((newsize+1)*qh_HASHfactor) | 0x1;  /* odd number */
  while (True) { 
    if ((size%3) && (size%5))
      break;
    size += 2;
    /* loop terminates because there is an infinite number of primes */
  }
  qh hash_table= qh_setnew (size);
  qh_setzero (qh hash_table, 0, size);
  return size;
} /* newhashtable */

/*----------------------------------------
-newvertex- creates and allocates space for a vertex
*/
vertexT *qh_newvertex(pointT *point) {
  vertexT *vertex;

  zinc_(Ztotvertices);
  vertex= (vertexT *)qh_memalloc(sizeof(vertexT));
  memset ((char *) vertex, 0, sizeof (vertexT));
  if (qh vertex_id == 0xFFFFFF) {
    fprintf(qh ferr, "qhull input error: more than %d vertices.  Id field overflows and two vertices\n\
may have the same identifier.  Vertices not sorted correctly.\n", 0xFFFFFF);
    qh_errexit(qh_ERRinput, NULL, NULL);
  }
  if (qh vertex_id == qh tracevertex_id)
    qh tracevertex= vertex;
  vertex->id= qh vertex_id++;
  vertex->point= point;
  trace4((qh ferr, "qh_newvertex: vertex p%d (v%d) created\n", qh_pointid(vertex->point), 
	  vertex->id));
  return (vertex);
} /* newvertex */

/*----------------------------------------
-nextridge3d- return next ridge and vertex for a 3d facet
  in qh_ORIENTclock order
  n^2 implementation to trace all ridges