qhull.c

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

At %02d:%02d:%02d & %2.5g CPU secs, qhull has created %d facets and merged %d.\n\
 The current hull contains %d facets and %d vertices.  Last point was p%d\n",
      tp->tm_hour, tp->tm_min, tp->tm_sec, cpu, qh facet_id -1,
      total, qh num_facets, qh num_vertices, qh furthest_id);
    return;
  }
  furthestid= qh_pointid (furthest);
  if (qh TRACEpoint == furthestid) {
    qh IStracing= qh TRACElevel;
    qhmem.IStracing= qh TRACElevel;
  }
  if (qh REPORTfreq && (qh facet_id-1 > qh lastreport+qh REPORTfreq)) {
    qh lastreport= qh facet_id-1;
    time (&timedata);
    tp= localtime (&timedata);
    cpu= (unsigned)clock() - qh hulltime;
    cpu /= qh_SECticks;
    total= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
    zinc_(Zdistio);
    qh_distplane (furthest, facet, &dist);
    fprintf (qh ferr, "\n\
At %02d:%02d:%02d & %2.5g CPU secs, qhull has created %d facets and merged %d.\n\
 The current hull contains %d facets and %d vertices.  There are %d\n\
 outside points.  Next is point p%d (v%d), %2.2g above f%d.\n",
      tp->tm_hour, tp->tm_min, tp->tm_sec, cpu, qh facet_id -1,
      total, qh num_facets, qh num_vertices, qh num_outside+1,
      furthestid, qh vertex_id, dist, getid_(facet));
  }else if (qh IStracing >=1) {
    cpu= (unsigned)clock() - qh hulltime;
    cpu /= qh_SECticks;
    qh_distplane (furthest, facet, &dist);
    fprintf (qh ferr, "qh_buildhull: add p%d (v%d) to hull of %d facets (%2.2g above f%d) and %d outside at %4.4g CPU secs.  Previous was p%d.\n",
      furthestid, qh vertex_id, qh num_facets, dist,
      getid_(facet), qh num_outside+1, cpu, qh furthest_id);
  }
  if (qh visit_id > (unsigned) INT_MAX) {
    qh visit_id= 0;
    FORALLfacets
      facet->visitid= qh visit_id;
  }
  if (qh vertex_visit > (unsigned) INT_MAX) {
    qh vertex_visit= 0;
    FORALLvertices
      vertex->visitid= qh vertex_visit;
  }
  qh furthest_id= furthestid;
} /* buildtracing */

/*-------------------------------------------
-errexit2- return exitcode to system after an error
  assumes exitcode non-zero
  for two facets, see qh_errexit() in user.c
*/
void qh_errexit2(int exitcode, facetT *facet, facetT *otherfacet) {
  
  qh_errprint("ERRONEOUS", facet, otherfacet, NULL, NULL);
  qh_errexit (exitcode, NULL, NULL);
} /* errexit2 */


/*-------------------------------------------------
-findhorizon- given a visible facet, find the point's horizon and visible facets
returns:
  qh visible_list to all visible facets 
    marks visible facets with ->visible 
    goodvisible counts visible->good
    initializes num_visible
notes:
  similar to delpoint()
*/
void qh_findhorizon(pointT *point, facetT *facet, int *goodvisible, int *goodhorizon) {
  facetT *neighbor, **neighborp, *visible;
  int numhorizon= 0, coplanar= 0;
  realT dist;
  
  trace1((qh ferr,"qh_findhorizon: find horizon for point p%d facet f%d\n",qh_pointid(point),facet->id));
  *goodvisible= *goodhorizon= 0;
  zinc_(Ztotvisible);
  qh_removefacet(facet);  /* visible_list at end of qh facet_list */
  qh_appendfacet(facet);
  qh num_visible= 1;
  if (facet->good)
    (*goodvisible)++;
  qh visible_list= facet;
  facet->visible= True;
  facet->f.replace= NULL;
  if (qh IStracing >=4)
    qh_errprint ("visible", facet, NULL, NULL, NULL);
  qh visit_id++;
  FORALLvisible_facets {
    visible->visitid= qh visit_id;
    FOREACHneighbor_(visible) {
      if (neighbor->visitid == qh visit_id) 
        continue;
      neighbor->visitid= qh visit_id;
      zzinc_(Znumvisibility);
      qh_distplane(point, neighbor, &dist);
      if (dist > qh MINvisible) {
        zinc_(Ztotvisible);
	qh_removefacet(neighbor);  /* append to end of qh visible_list */
	qh_appendfacet(neighbor);
	neighbor->visible= True;
        neighbor->f.replace= NULL;
	qh num_visible++;
	if (neighbor->good)
	  (*goodvisible)++;
        if (qh IStracing >=4)
          qh_errprint ("visible", neighbor, NULL, NULL, NULL);
      }else {
 	if (dist > - qh MAXcoplanar) {
    	  neighbor->coplanar= True;
          zzinc_(Zcoplanarhorizon);
	  coplanar++;
	  if (qh MERGING) {
	    if (dist > 0) {
	      maximize_(qh max_outside, dist);
	      maximize_(qh max_vertex, dist);
#if qh_MAXoutside
	      maximize_(neighbor->maxoutside, dist);
#endif
	    }else
	      minimize_(qh min_vertex, dist);  /* due to merge later */
	  }
      	  trace2((qh ferr, "qh_findhorizon: point p%d is coplanar to horizon f%d, dist=%2.7g < qh MINvisible (%2.7g)\n",
	      qh_pointid(point), neighbor->id, dist, qh MINvisible));
	}else
    	  neighbor->coplanar= False;
    	zinc_(Ztothorizon);
        numhorizon++;
	if (neighbor->good)
	  (*goodhorizon)++;
        if (qh IStracing >=4)
          qh_errprint ("horizon", neighbor, NULL, NULL, NULL);
      }
    }
  }
  if (!numhorizon) {
    fprintf(qh ferr, "qhull precision error (qh_findhorizon): empty horizon\n\
Point p%d was above all facets.\n", qh_pointid(point));
    qh_printfacetlist (qh facet_list, NULL, True);
    qh_errexit(qh_ERRprec, NULL, NULL);
  }
  trace1((qh ferr, "qh_findhorizon: %d horizon facets (good %d), %d visible (good %d), %d coplanar\n", 
       numhorizon, *goodhorizon, qh num_visible, *goodvisible, coplanar));
  if (qh IStracing >= 4 && qh num_facets < 50) 
    qh_printlists ();
} /* findhorizon */


/*------------------------------------------------
-nextfurthest- returns next furthest point for processing
returns:
  NULL if none available
  visible facet for furthest
  removes empty outside sets  
*/
pointT *qh_nextfurthest (facetT **visible) {
  facetT *facet;
  int size, index;
  realT randr;
  pointT *furthest;

  while ((facet= qh facet_next) != qh facet_tail) {
    if (!facet->outsideset) {
      qh facet_next= facet->next;
      continue;
    }
    SETreturnsize_(facet->outsideset, size);
    if (!size) {
      qh_setfree (&facet->outsideset);
      qh facet_next= facet->next;
      continue;
    }
    if (!qh RANDOMoutside && !qh VIRTUALmemory) {
      *visible= facet;
      return ((pointT*)qh_setdellast (facet->outsideset));
    }
    if (qh RANDOMoutside) {
      randr= qh_RANDOMint;
      randr= randr/(qh_RANDOMmax+1);
      index= floor(qh num_outside * 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 incorrect or random %2.2g >= 1.0\n",
              qh num_outside, 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 */

/*-------------------------------------------------
-partitionall- partitions all points into the outsidesets of facets
   vertices= set of vertices used by qh facet_list
     does not partition qh GOODpoint
     if ONLYgood && !MERGING, does not partition GOODvertex
   all facets have ->newfacet for qh_findbestnew in qh_partitionpoint
notes:
   faster if qh facet_list sorted by anticipated size of outside set
*/
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);
  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 */


/*-------------------------------------------------
-partitioncoplanar- partition coplanar point to a facet
  if dist NULL, searches for bestfacet, and does nothing if inside
      if qh findbestnew, searches new facets instead of findbest
returns: