qhull.c

关于网格剖分的

  
  qh_buildhull()
    construct a convex hull by adding outside points one at a time

  returns:
  
  notes:
    may be called multiple times
    checks facet and vertex lists for incorrect flags
    to recover from STOPcone, call qh_deletevisible and qh_resetlists

  design:
    check visible facet and newfacet flags
    check newlist vertex flags and qh.STOPcone/STOPpoint
    for each facet with a furthest outside point
      add point to facet
      exit if qh.STOPcone or qh.STOPpoint requested
    if qh.NARROWhull for initial simplex
      partition remaining outside points to coplanar sets
*/
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 NARROWhull) /* move points from outsideset to coplanarset */
    qh_outcoplanar( /* facet_list */ );
  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 */
  

/*-<a                             href="qh-c.htm#qhull"
  >-------------------------------</a><a name="buildtracing">-</a>
  
  qh_buildtracing( furthest, facet )
    trace an iteration of qh_buildhull() for furthest point and facet
    if !furthest, prints progress message

  returns:
    tracks progress with qh.lastreport
    updates qh.furthest_id (-3 if furthest is NULL)
    also resets visit_id, vertext_visit on wrap around

  see:
    qh_tracemerging()

  design:
    if !furthest
      print progress message
      exit
    if 'TFn' iteration
      print progress message
    else if tracing
      trace furthest point and facet
    reset qh.visit_id and qh.vertex_visit if overflow may occur
    set qh.furthest_id for tracing
*/
void qh_buildtracing (pointT *furthest, facetT *facet) {
  realT dist= 0;
  float cpu;
  int total, furthestid;
  time_t timedata;
  struct tm *tp;
  vertexT *vertex;

  qh old_randomdist= qh RANDOMdist;
  qh RANDOMdist= False;
  if (!furthest) {
    time (&timedata);
    tp= localtime (&timedata);
    cpu= qh_CPUclock - qh hulltime;
    cpu /= qh_SECticks;
    total= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
    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.  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= qh_CPUclock - 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= qh_CPUclock - qh hulltime;
    cpu /= qh_SECticks;
    qh_distplane (furthest, facet, &dist);
    fprintf (qh ferr, "qh_addpoint: 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;
  qh RANDOMdist= qh old_randomdist;
} /* buildtracing */

/*-<a                             href="qh-c.htm#qhull"
  >-------------------------------</a><a name="errexit2">-</a>
  
  qh_errexit2( exitcode, facet, otherfacet )
    return exitcode to system after an error
    report two facets

  returns:
    assumes exitcode non-zero

  see:
    normally use qh_errexit() in user.c (reports a facet and a ridge)
*/
void qh_errexit2(int exitcode, facetT *facet, facetT *otherfacet) {
  
  qh_errprint("ERRONEOUS", facet, otherfacet, NULL, NULL);
  qh_errexit (exitcode, NULL, NULL);
} /* errexit2 */


/*-<a                             href="qh-c.htm#qhull"
  >-------------------------------</a><a name="findhorizon">-</a>
  
  qh_findhorizon( point, facet, goodvisible, goodhorizon )
    given a visible facet, find the point's horizon and visible facets

  returns:
    returns qh.visible_list/num_visible with all visible facets 
      marks visible facets with ->visible 
    updates count of good visible and good horizon facets
    updates qh.max_outside, qh.max_vertex, facet->maxoutside

  see:
    similar to qh_delpoint()

  design:
    move facet to qh.visible_list at end of qh.facet_list
    for all visible facets
     for each unvisited neighbor of a visible facet
       compute distance of point to neighbor
       if point above neighbor
         move neighbor to end of qh.visible_list
       else if point is coplanar with neighbor
         update qh.max_outside, qh.max_vertex, neighbor->maxoutside
         mark neighbor coplanar (will create a samecycle later)
         update horizon statistics         
*/
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);
          qh_precision ("coplanar horizon");
	  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) {
    qh_precision ("empty horizon");
    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 */

/*-<a                             href="qh-c.htm#qhull"
  >-------------------------------</a><a name="nextfurthest">-</a>
  
  qh_nextfurthest( visible )
    returns next furthest point and visible facet for qh_addpoint()
    starts search at qh.facet_next

  returns:
    removes furthest point from outside set
    NULL if none available
    advances qh.facet_next over facets with empty outside sets  

  design:
    for each facet from qh.facet_next
      if empty outside set
        advance qh.facet_next
      else if qh.NARROWhull
        determine furthest outside point
        if furthest point is not outside
          advance qh.facet_next (point will be coplanar)
    remove furthest point from outside set
*/
pointT *qh_nextfurthest (facetT **visible) {
  facetT *facet;
  int size, index;
  realT randr, dist;
  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 NARROWhull) {
      if (facet->notfurthest) 
	qh_furthestout (facet);
      furthest= (pointT*)qh_setlast (facet->outsideset);
#if qh_COMPUTEfurthest
      qh_distplane (furthest, facet, &dist);
      zinc_(Zcomputefurthest);
#else