qhull.c

DT三角化实现

    }
    qh num_outside++;
    trace4((qh ferr, "qh_partitionpoint: point p%d is outside facet f%d new? %d(or narrowhull)\n",
	  qh_pointid(point), bestfacet->id, bestfacet->newfacet));
  }else if (qh DELAUNAY || bestdist >= -qh MAXcoplanar) { /* for 'd', bestdist skips upperDelaunay facets */
    zzinc_(Zcoplanarpart);
    if (qh DELAUNAY)
      qh_precision ("nearly incident point");
    if ((qh KEEPcoplanar + qh KEEPnearinside) || bestdist > qh max_outside) 
      qh_partitioncoplanar (point, bestfacet, &bestdist);
    else {
      trace4((qh ferr, "qh_partitionpoint: point p%d is coplanar to facet f%d (dropped)\n",
	  qh_pointid(point), bestfacet->id));
    }
  }else if (qh KEEPnearinside && bestdist > -qh NEARinside) {
    zinc_(Zpartnear);
    qh_partitioncoplanar (point, bestfacet, &bestdist);
  }else {
    zinc_(Zpartinside);
    trace4((qh ferr, "qh_partitionpoint: point p%d is inside all facets, closest to f%d dist %2.2g\n",
	  qh_pointid(point), bestfacet->id, bestdist));
    if (qh KEEPinside)
      qh_partitioncoplanar (point, bestfacet, &bestdist);
  }
} /* partitionpoint */

/*-<a                             href="qh-qhull.htm#TOC"
  >-------------------------------</a><a name="partitionvisible">-</a>
  
  qh_partitionvisible( allpoints, numoutside )
    partitions points in visible facets to qh.newfacet_list
    qh.visible_list= visible facets
    for visible facets
      1st neighbor (if any) points to a horizon facet or a new facet
    if allpoints (not used),
      repartitions coplanar points

  returns:
    updates outside sets and coplanar sets of qh.newfacet_list
    updates qh.num_outside (count of outside points)
  
  notes:
    qh.findbest_notsharp should be clear (extra work if set)

  design:
    for all visible facets with outside set or coplanar set
      select a newfacet for visible facet
      if outside set
        partition outside set into new facets
      if coplanar set and keeping coplanar/near-inside/inside points
        if allpoints
          partition coplanar set into new facets, may be assigned outside
        else
          partition coplanar set into coplanar sets of new facets
    for each deleted vertex
      if allpoints
        partition vertex into new facets, may be assigned outside
      else
        partition vertex into coplanar sets of new facets
*/
void qh_partitionvisible(/*visible_list*/ boolT allpoints, int *numoutside) {
  facetT *visible, *newfacet;
  pointT *point, **pointp;
  int coplanar=0, size;
  unsigned count;
  vertexT *vertex, **vertexp;
  
  if (qh ONLYmax)
    maximize_(qh MINoutside, qh max_vertex);
  *numoutside= 0;
  FORALLvisible_facets {
    if (!visible->outsideset && !visible->coplanarset)
      continue;
    newfacet= visible->f.replace;
    count= 0;
    while (newfacet && newfacet->visible) {
      newfacet= newfacet->f.replace;
      if (count++ > qh facet_id)
	qh_infiniteloop (visible);
    }
    if (!newfacet)
      newfacet= qh newfacet_list;
    if (newfacet == qh facet_tail) {
      fprintf (qh ferr, "qhull precision error (qh_partitionvisible): all new facets deleted as\n        degenerate facets. Can not continue.\n");
      qh_errexit (qh_ERRprec, NULL, NULL);
    }
    if (visible->outsideset) {
      size= qh_setsize (visible->outsideset);
      *numoutside += size;
      qh num_outside -= size;
      FOREACHpoint_(visible->outsideset) 
        qh_partitionpoint (point, newfacet);
    }
    if (visible->coplanarset && (qh KEEPcoplanar + qh KEEPinside + qh KEEPnearinside)) {
      size= qh_setsize (visible->coplanarset);
      coplanar += size;
      FOREACHpoint_(visible->coplanarset) {
        if (allpoints) /* not used */
          qh_partitionpoint (point, newfacet);
        else
          qh_partitioncoplanar (point, newfacet, NULL);
      }
    }
  }
  FOREACHvertex_(qh del_vertices) {
    if (vertex->point) {
      if (allpoints) /* not used */
        qh_partitionpoint (vertex->point, qh newfacet_list);
      else
        qh_partitioncoplanar (vertex->point, qh newfacet_list, NULL);
    }
  }
  trace1((qh ferr,"qh_partitionvisible: partitioned %d points from outsidesets and %d points from coplanarsets\n", *numoutside, coplanar));
} /* partitionvisible */



/*-<a                             href="qh-qhull.htm#TOC"
  >-------------------------------</a><a name="precision">-</a>
  
  qh_precision( reason )
    restart on precision errors if not merging and if 'QJn'
*/
void qh_precision (char *reason) {

  if (qh ALLOWrestart && !qh PREmerge && !qh MERGEexact) {
    if (qh JOGGLEmax < REALmax/2) {
      trace0((qh ferr, "qh_precision: qhull restart because of %s\n", reason));
      longjmp(qh restartexit, qh_ERRprec);
    }
  }
} /* qh_precision */

/*-<a                             href="qh-qhull.htm#TOC"
  >-------------------------------</a><a name="printsummary">-</a>
  
  qh_printsummary( fp )
    prints summary to fp

  notes:
    not in io.c so that user_eg.c can prevent io.c from loading
    qh_printsummary and qh_countfacets must match counts

  design:
    determine number of points, vertices, and coplanar points
    print summary
*/
void qh_printsummary(FILE *fp) {
  realT ratio, outerplane, innerplane;
  float cpu;
  int size, id, nummerged, numvertices, numcoplanars= 0, nonsimplicial=0;
  int goodused;
  facetT *facet;
  char *s;
  int numdel= zzval_(Zdelvertextot);
  int numtricoplanars= 0;

  size= qh num_points + qh_setsize (qh other_points);
  numvertices= qh num_vertices - qh_setsize (qh del_vertices);
  id= qh_pointid (qh GOODpointp);
  FORALLfacets {
    if (facet->coplanarset)
      numcoplanars += qh_setsize( facet->coplanarset);
    if (facet->good) {
      if (facet->simplicial) {
	if (facet->keepcentrum && facet->tricoplanar)
	  numtricoplanars++;
      }else if (qh_setsize(facet->vertices) != qh hull_dim)
	nonsimplicial++;
    }
  }
  if (id >=0 && qh STOPcone-1 != id && -qh STOPpoint-1 != id)
    size--;
  if (qh STOPcone || qh STOPpoint)
      fprintf (fp, "\nAt a premature exit due to 'TVn', 'TCn', 'TRn', or precision error.");
  if (qh UPPERdelaunay)
    goodused= qh GOODvertex + qh GOODpoint + qh SPLITthresholds;
  else if (qh DELAUNAY)
    goodused= qh GOODvertex + qh GOODpoint + qh GOODthreshold;
  else
    goodused= qh num_good;
  nummerged= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
  if (qh VORONOI) {
    if (qh UPPERdelaunay)
      fprintf (fp, "\n\
Furthest-site Voronoi vertices by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
    else
      fprintf (fp, "\n\
Voronoi diagram by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
    fprintf(fp, "  Number of Voronoi regions%s: %d\n",
              qh ATinfinity ? " and at-infinity" : "", numvertices);
    if (numdel)
      fprintf(fp, "  Total number of deleted points due to merging: %d\n", numdel); 
    if (numcoplanars - numdel > 0) 
      fprintf(fp, "  Number of nearly incident points: %d\n", numcoplanars - numdel); 
    else if (size - numvertices - numdel > 0) 
      fprintf(fp, "  Total number of nearly incident points: %d\n", size - numvertices - numdel); 
    fprintf(fp, "  Number of%s Voronoi vertices: %d\n", 
              goodused ? " 'good'" : "", qh num_good);
    if (nonsimplicial) 
      fprintf(fp, "  Number of%s non-simplicial Voronoi vertices: %d\n", 
              goodused ? " 'good'" : "", nonsimplicial);
  }else if (qh DELAUNAY) {
    if (qh UPPERdelaunay)
      fprintf (fp, "\n\
Furthest-site Delaunay triangulation by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
    else
      fprintf (fp, "\n\
Delaunay triangulation by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
    fprintf(fp, "  Number of input sites%s: %d\n", 
              qh ATinfinity ? " and at-infinity" : "", numvertices);
    if (numdel)
      fprintf(fp, "  Total number of deleted points due to merging: %d\n", numdel); 
    if (numcoplanars - numdel > 0) 
      fprintf(fp, "  Number of nearly incident points: %d\n", numcoplanars - numdel); 
    else if (size - numvertices - numdel > 0)
      fprintf(fp, "  Total number of nearly incident points: %d\n", size - numvertices - numdel); 
    fprintf(fp, "  Number of%s Delaunay regions: %d\n", 
              goodused ? " 'good'" : "", qh num_good);
    if (nonsimplicial) 
      fprintf(fp, "  Number of%s non-simplicial Delaunay regions: %d\n", 
              goodused ? " 'good'" : "", nonsimplicial);
  }else if (qh HALFspace) {
    fprintf (fp, "\n\
Halfspace intersection by the convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
    fprintf(fp, "  Number of halfspaces: %d\n", size);
    fprintf(fp, "  Number of non-redundant halfspaces: %d\n", numvertices);
    if (numcoplanars) {
      if (qh KEEPinside && qh KEEPcoplanar)
      	s= "similar and redundant";
      else if (qh KEEPinside)
        s= "redundant";
      else
        s= "similar"; 
      fprintf(fp, "  Number of %s halfspaces: %d\n", s, numcoplanars);
    } 
    fprintf(fp, "  Number of intersection points: %d\n", qh num_facets - qh num_visible);
    if (goodused)
      fprintf(fp, "  Number of 'good' intersection points: %d\n", qh num_good);
    if (nonsimplicial) 
      fprintf(fp, "  Number of%s non-simplicial intersection points: %d\n", 
              goodused ? " 'good'" : "", nonsimplicial);
  }else {
    fprintf (fp, "\n\
Convex hull of %d points in %d-d:\n\n", size, qh hull_dim);
    fprintf(fp, "  Number of vertices: %d\n", numvertices);
    if (numcoplanars) {
      if (qh KEEPinside && qh KEEPcoplanar)
      	s= "coplanar and interior";
      else if (qh KEEPinside)
        s= "interior";
      else
        s= "coplanar"; 
      fprintf(fp, "  Number of %s points: %d\n", s, numcoplanars);
    } 
    fprintf(fp, "  Number of facets: %d\n", qh num_facets - qh num_visible);
    if (goodused)
      fprintf(fp, "  Number of 'good' facets: %d\n", qh num_good);
    if (nonsimplicial) 
      fprintf(fp, "  Number of%s non-simplicial facets: %d\n", 
              goodused ? " 'good'" : "", nonsimplicial);
  }
  if (numtricoplanars)
      fprintf(fp, "  Number of triangulated facets: %d\n", numtricoplanars);
  fprintf(fp, "\nStatistics for: %s | %s", 
                      qh rbox_command, qh qhull_command);
  if (qh ROTATErandom != INT_MIN)
    fprintf(fp, " QR%d\n\n", qh ROTATErandom);
  else
    fprintf(fp, "\n\n");
  fprintf(fp, "  Number of points processed: %d\n", zzval_(Zprocessed));
  fprintf(fp, "  Number of hyperplanes created: %d\n", zzval_(Zsetplane));
  if (qh DELAUNAY)
    fprintf(fp, "  Number of facets in hull: %d\n", qh num_facets - qh num_visible);
  fprintf(fp, "  Number of distance tests for qhull: %d\n", zzval_(Zpartition)+
      zzval_(Zpartitionall)+zzval_(Znumvisibility)+zzval_(Zpartcoplanar));
#if 0  /* NOTE: must print before printstatistics() */
  {realT stddev, ave;
  fprintf(fp, "  average new facet balance: %2.2g\n",
	  wval_(Wnewbalance)/zval_(Zprocessed));
  stddev= qh_stddev (zval_(Zprocessed), wval_(Wnewbalance), 
                                 wval_(Wnewbalance2), &ave);
  fprintf(fp, "  new facet standard deviation: %2.2g\n", stddev);
  fprintf(fp, "  average partition balance: %2.2g\n",
	  wval_(Wpbalance)/zval_(Zpbalance));
  stddev= qh_stddev (zval_(Zpbalance), wval_(Wpbalance), 
                                 wval_(Wpbalance2), &ave);
  fprintf(fp, "  partition standard deviation: %2.2g\n", stddev);
  }
#endif
  if (nummerged) {
    fprintf(fp,"  Number of distance tests for merging: %d\n",zzval_(Zbestdist)+
          zzval_(Zcentrumtests)+zzval_(Zdistconvex)+zzval_(Zdistcheck)+
          zzval_(Zdistzero));
    fprintf(fp,"  Number of distance tests for checking: %d\n",zzval_(Zcheckpart));
    fprintf(fp,"  Number of merged facets: %d\n", nummerged);
  }
  if (!qh RANDOMoutside && qh QHULLfinished) {
    cpu= qh hulltime;
    cpu /= qh_SECticks;
    wval_(Wcpu)= cpu;
    fprintf (fp, "  CPU seconds to compute hull (after input): %2.4g\n", cpu);
  }
  if (qh RERUN) {
    if (!qh PREmerge && !qh MERGEexact)
      fprintf(fp, "  Percentage of runs with precision errors: %4.1f\n",
	   zzval_(Zretry)*100.0/qh build_cnt);  /* careful of order */
  }else if (qh JOGGLEmax < REALmax/2) {
    if (zzval_(Zretry))
      fprintf(fp, "  After %d retries, input joggled by: %2.2g\n",
         zzval_(Zretry), qh JOGGLEmax);
    else
      fprintf(fp, "  Input joggled by: %2.2g\n", qh JOGGLEmax);
  }
  if (qh totarea != 0.0) 
    fprintf(fp, "  %s facet area:   %2.8g\n", 
	    zzval_(Ztotmerge) ? "Approximate" : "Total", qh totarea);
  if (qh totvol != 0.0) 
    fprintf(fp, "  %s volume:       %2.8g\n", 
	    zzval_(Ztotmerge) ? "Approximate" : "Total", qh totvol);
  if (qh MERGING) {
    qh_outerinner (NULL, &outerplane, &innerplane);
    if (outerplane > 2 * qh DISTround) {
      fprintf(fp, "  Maximum distance of %spoint above facet: %2.2g", 
	    (qh QHULLfinished ? "" : "merged "), outerplane);
      ratio= outerplane/(qh ONEmerge + qh DISTround);
      /* don't report ratio if MINoutside is large */
      if (ratio > 0.05 && 2* qh ONEmerge > qh MINoutside && qh JOGGLEmax > REALmax/2)
        fprintf (fp, " (%.1fx)\n", ratio);
      else
        fprintf (fp, "\n");
    }
    if (innerplane < -2 * qh DISTround) {
      fprintf(fp, "  Maximum distance of %svertex below facet: %2.2g", 
	    (qh QHULLfinished ? "" : "merged "), innerplane);
      ratio= -innerplane/(qh ONEmerge+qh DISTround);
      if (ratio > 0.05 && qh JOGGLEmax > REALmax/2)
        fprintf (fp, " (%.1fx)\n", ratio);
      else
        fprintf (fp, "\n");
    }
  }
  fprintf(fp, "\n");
} /* printsummary */