qhull.c

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

/* qhull - Quickhull algorithm for convex hulls

   qhull() and top-level routines

   see README, qhull.h, unix.c and mac.c

   see qhull_a.h for internal functions
   
   copyright (c) 1993-1995 The Geometry Center        
*/

#include "qhull_a.h" 


/*-------------------------------------------------
-qhull- hull_dim convex hull of num_points starting at first_point
returns:
  returns facet_list, numfacets, etc. 
*/
void qh_qhull (void) {
  setT *maxpoints, *vertices;
  facetT *facet;
  int numpart, i, numoutside;
  realT dist;
  boolT isoutside;

  qh hulltime= (unsigned)clock();
  if (qh DELAUNAY && qh upper_threshold[qh hull_dim-1] > REALmax/2
                  && qh lower_threshold[qh hull_dim-1] < -REALmax/2) {
    for (i= qh_PRINTEND; i--; ) {
      if (qh PRINTout[i] == qh_PRINTgeom && qh DROPdim < 0 
 	  && !qh GOODthreshold && !qh SPLITthresholds)
	break;  /* in this case, don't set upper_threshold */
    }
    if (i < 0) {
      qh upper_threshold[qh hull_dim-1]= 0.0;
      if (!qh GOODthreshold)
	qh SPLITthresholds= True;
    }
  }
  maxpoints= qh_maxmin(qh first_point, qh num_points, qh hull_dim);
  /* qh_maxmin sets DISTround and other precision constants */
  if (qh PRINToptions1st || qh TRACElevel || qh IStracing) {
    if (qh TRACElevel || qh IStracing)
      fprintf (qh ferr, "\nTrace level %d for %s | %s\n", 
         qh IStracing ? qh IStracing : qh TRACElevel, qh rbox_command, qh qhull_command);
    fprintf (qh ferr, "Options selected for qhull %s:\n%s\n", qh_version, qh qhull_options);
  }
  vertices= qh_initialvertices(qh hull_dim, maxpoints, qh first_point, qh num_points); 
  qh_initialhull (vertices);  /* initial qh facet_list */
  qh_partitionall (vertices, qh first_point, qh num_points);
  if (qh PREmerge) {
    qh cos_max= qh premerge_cos;
    qh centrum_radius= qh premerge_centrum;
  }
  if (qh ONLYgood) {
    if (!(qh GOODthreshold || qh GOODpoint
	  || (qh GOODvertex > 0 && !qh MERGING))) {
      fprintf (qh ferr, "qhull input error: 'Qg' (ONLYgood) needs a good threshold ('Pd0D0'), a\n\
good point (QGn or QG-n), or a good vertex without merging (QVn).\n");
      qh_errexit (qh_ERRinput, NULL, NULL);
    }
    if (qh GOODvertex > 0  && !qh MERGING  /* matches qh_partitionall */
	&& !qh_isvertex (qh GOODvertexp, vertices)) {
      facet= qh_findbestnew (qh GOODvertexp, qh facet_list, 
			  &dist, &isoutside, &numpart);
      zadd_(Zdistgood, numpart);
      if (!isoutside) {
        fprintf (qh ferr, "qhull input error: point for QV%d is inside initial simplex\n",
	       qh_pointid(qh GOODvertexp));
        qh_errexit (qh_ERRinput, NULL, NULL);
      }
      if (!qh_addpoint (qh GOODvertexp, facet, False)) {
	qh_settempfree(&vertices);
	qh_settempfree(&maxpoints);
	return;
      }
    }
    qh_findgood (qh facet_list, 0);
  }
  qh_settempfree(&vertices);
  qh_settempfree(&maxpoints);
  qh_resetlists (False /*qh visible_list newvertex_list newfacet_list */);
  qh_buildhull();
  if (!qh STOPpoint && !qh STOPcone) {
    if (qh ZEROall_ok && !qh TESTvneighbors && qh MERGEexact)
      qh_checkzero( qh_ALL);
    if (qh ZEROall_ok && !qh TESTvneighbors) {
      trace2((qh ferr, "qh_qhull: all facets are clearly convex.  Post-merging not needed.\n"));
    }else {
      if (qh MERGEexact || (qh hull_dim > qh_DIMreduceBuild && qh PREmerge))
        qh_postmerge ("First post-merge", qh premerge_centrum, qh premerge_cos, 
             (qh POSTmerge ? False : qh TESTvneighbors));
      else if (!qh POSTmerge && qh TESTvneighbors) 
        qh_postmerge ("For testing vertex neighbors", qh premerge_centrum,
             qh premerge_cos, True); 
      if (qh POSTmerge)
        qh_postmerge ("For post-merging", qh postmerge_centrum, 
             qh postmerge_cos, qh TESTvneighbors);
      if (qh visible_list == qh facet_list) { /* i.e., merging done */
        qh findbestnew= True;
        qh_partitionvisible (/*visible_list, newfacet_list*/ !qh_ALL, &numoutside);
        qh findbestnew= False;
        qh_deletevisible (/*qh visible_list*/);
        qh_resetlists (False /*qh visible_list newvertex_list newfacet_list */);
      }
      if (qh DOcheckmax){
        if (qh REPORTfreq) {
	  qh_buildtracing (NULL, NULL); 
  	  fprintf (qh ferr, "\nTesting all coplanar points.\n");
        }
        qh_check_maxout();
      }
    }
  }
  if (qh_setsize (qhmem.tempstack) != 0) {
    fprintf (qh ferr, "qhull internal error (qh_qhull): temporary sets not empty (%d)\n",
	     qh_setsize (qhmem.tempstack));
    qh_errexit (qh_ERRqhull, NULL, NULL);
  }
  qh hulltime= (unsigned)clock() - qh hulltime;
  qh QHULLfinished= True;
  trace1((qh ferr, "qh_qhull: algorithm completed\n"));
} /* qhull */

/*-------------------------------------------------
-addpoint-  add point to hull above a facet
  if !facet, locates a facet for the point
    !facet works for Delaunay triangulations
    !facet does not work for lens-shaped hulls
    if point is not outside of the hull, uses qh_partitioncoplanar()
  if checkdist, checks that point is outside of the facet.
    if not outside, uses qh_partitioncoplanar()
    if facet->upperdelaunay, assumes that the point is above facet
  if !checkdist and facet, assumes point is above facet (major damage if below)
returns:
  if unknown point, adds it to qh other_points
  returns False if user requested break
     visible_list, newfacet_list, delvertex_list, NEWfacets may be defined
*/
boolT qh_addpoint (pointT *furthest, facetT *facet, boolT checkdist) {
  int goodvisible, goodhorizon;
  vertexT *vertex;
  facetT *newfacet;
  realT dist, newbalance, pbalance;
  boolT isoutside= False;
  int numpart, numpoints, numnew, firstnew;

  if (qh_pointid (furthest) == -1)
    qh_setappend (&qh other_points, furthest);
  if (!facet 
  || (checkdist && !facet->upperdelaunay)) { /* else missed by qh_findbest */
    if (!facet)
      facet= qh_nonupper( qh facet_list);
    facet= qh_findbest (furthest, facet, !qh_ALL, False, &dist, &isoutside, &numpart);
    zzadd_(Zpartition, numpart);
    if (!isoutside) {
      zinc_(Znotmax);  /* last point of outsideset is no longer furthest. */
      qh_partitioncoplanar (furthest, facet, &dist);
      return True;
    }
  }
  qh_buildtracing (furthest, facet);
  if (qh STOPpoint < 0 && qh furthest_id == -qh STOPpoint-1)
    return False;
  qh_findhorizon (furthest, facet, &goodvisible, &goodhorizon); 
  if (qh ONLYgood && !(goodvisible+goodhorizon)) {
    zinc_(Znotgood);  
    /* last point of outsideset is no longer furthest.  This is ok
       since all points of the outside are likely to be bad */
    qh_resetlists (False /*qh visible_list newvertex_list newfacet_list */);
    return True;
  }
  zzinc_(Zprocessed);
  firstnew= qh facet_id;
  vertex= qh_makenewfacets (furthest /*visible_list, attaches if !ONLYgood */);
  qh_makenewplanes (/* newfacet_list */);
  numnew= qh facet_id - firstnew;
  newbalance= numnew - (realT) (qh num_facets-qh num_visible)
                         * qh hull_dim/qh num_vertices;
  wadd_(Wnewbalance, newbalance);
  wadd_(Wnewbalance2, newbalance * newbalance);
  if (qh ONLYgood && !qh_findgood (qh newfacet_list, goodhorizon)) {
    FORALLnew_facets 
      qh_delfacet (newfacet);
    qh_delvertex (vertex);
    qh_resetlists (True /*qh visible_list newvertex_list newfacet_list */);
    zinc_(Znotgoodnew);
    return True;
  }
  if (qh ONLYgood)
    qh_attachnewfacets(/*visible_list*/);
  qh_matchnewfacets();
  qh_updatevertices();
  if (qh STOPcone && qh furthest_id == qh STOPcone-1)
    return False;  /* visible_list etc. still defined */
  if (qh PREmerge || qh MERGEexact) {
    qh_premerge (vertex, qh premerge_centrum, qh premerge_cos);
    if (zzval_(Ztotmerge) > qh_USEfindbestnew)
      qh findbestnew= True;
    else {
      FORALLnew_facets {
	if (!newfacet->simplicial) {
	  qh findbestnew= True;  /* qh_findbest can not be used */
	  break;
	}
      }
    }
  }else if (qh BESToutside)
    qh findbestnew= True;
  qh_partitionvisible (/*visible_list, newfacet_list*/ !qh_ALL, &numpoints);
  qh findbestnew= False;
  qh findbest_notsharp= False;
  zinc_(Zpbalance);
  pbalance= numpoints - (realT) qh hull_dim /* assumes all points extreme */
                * (qh num_points - qh num_vertices)/qh num_vertices;
  wadd_(Wpbalance, pbalance);
  wadd_(Wpbalance2, pbalance * pbalance);
  qh_deletevisible (/*qh visible_list*/);
  zmax_(Zmaxvertex, qh num_vertices);
  qh NEWfacets= False;
  if (qh IStracing >= 4)
    qh_printfacetlist (qh newfacet_list, NULL, True);
  if (qh CHECKfrequently) {
    if (qh num_facets < 50)
      qh_checkpolygon (qh facet_list);
    else
      qh_checkpolygon (qh newfacet_list);
  }
  if (qh STOPpoint > 0 && qh furthest_id == qh STOPpoint-1)
    return False; 
  qh_resetlists (True /*qh visible_list newvertex_list newfacet_list */);
  trace2((qh ferr, "qh_addpoint: added p%d new facets %d new balance %2.2g point balance %2.2g\n",
    qh_pointid (furthest), numnew, newbalance, pbalance));
  return True;
} /* addpoint */

/*-------------------------------------------------
-buildhull- constructs a hull by adding outside points one at a time
  may be called multiple times
  checks facet and vertex lists for 'visible', 'newfacet', and 'newlist'
notes:
  to recover from STOPcone, call qh_deletevisible and qh_resetlists
*/
void qh_buildhull(void) {
  facetT *facet;
  pointT *furthest;
  vertexT *vertex;
  int id;
  
  trace1((qh ferr, "qh_buildhull: start build hull\n"));
  FORALLfacets {
    if (facet->visible || facet->newfacet) {
      fprintf (qh ferr, "qhull internal error (qh_buildhull): visible or new facet f%d in facet list\n",
                   facet->id);    
      qh_errexit (qh_ERRqhull, facet, NULL);
    }
  }
  FORALLvertices {
    if (vertex->newlist) {
      fprintf (qh ferr, "qhull internal error (qh_buildhull): new vertex f%d in vertex list\n",
                   vertex->id);
      qh_errprint ("ERRONEOUS", NULL, NULL, NULL, vertex);
      qh_errexit (qh_ERRqhull, NULL, NULL);
    }
    id= qh_pointid (vertex->point);
    if ((qh STOPpoint>0 && id == qh STOPpoint-1) ||
	(qh STOPpoint<0 && id == -qh STOPpoint-1) ||
	(qh STOPcone>0 && id == qh STOPcone-1)) {
      trace1((qh ferr,"qh_buildhull: stop point or cone P%d in initial hull\n", id));
      return;
    }
  }
  qh facet_next= qh facet_list;      /* advance facet when processed */
  while ((furthest= qh_nextfurthest (&facet))) {
    qh num_outside--;  /* if ONLYmax, furthest may not be outside */
    if (!qh_addpoint (furthest, facet, qh ONLYmax))
      break;
  }
  if (qh num_outside && !furthest) {
    fprintf (qh ferr, "qhull internal error (qh_buildhull): %d outside points were never processed.\n", qh num_outside);
    qh_errexit (qh_ERRqhull, NULL, NULL);
  }
  trace1((qh ferr, "qh_buildhull: completed the hull construction\n"));
} /* buildhull */
  

/*-------------------------------------------
-buildtracing- for tracing execution of buildhull
  tracks progress with qh lastreport
  updates qh furthest_id (-3 if furthest is NULL)
  also resets visit_id, vertext_visit on wrap around
  if !furthest, prints basic message
  see also qh_tracemerging()
*/
void qh_buildtracing (pointT *furthest, facetT *facet) {
  realT dist= 0;
  float cpu;
  int total, furthestid;
  time_t timedata;
  struct tm *tp;
  vertexT *vertex;

  if (!furthest) {
    time (&timedata);
    tp= localtime (&timedata);
    cpu= (unsigned)clock() - qh hulltime;
    cpu /= qh_SECticks;
    total= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
    fprintf (qh ferr, "\n\