qhull.c

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

  max_ouside updated
  facet->maxoutside's updated at end by qh_check_maxout
  if KEEPcoplanar or KEEPinside
    point assigned to best coplanarset
*/
void qh_partitioncoplanar (pointT *point, facetT *facet, realT *dist) {
  facetT *bestfacet;
  pointT *oldfurthest;
  realT bestdist, dist2;
  int numpart= 0;
  boolT isoutside, istrace= False;

  if (!dist) {
    if (qh findbestnew)
      bestfacet= qh_findbestnew (point, facet, 
			  &bestdist, NULL, &numpart);
    else
      bestfacet= qh_findbest (point, facet, qh_ALL, False, 
                          &bestdist, &isoutside, &numpart);
    zinc_(Ztotpartcoplanar);
    zzadd_(Zpartcoplanar, numpart);
    if (qh KEEPnearinside) {
      if (bestdist < -qh NEARinside) { 
        zinc_(Zcoplanarinside);
        return;
      }
    }else if (!qh KEEPinside && bestdist < -qh MAXcoplanar) {
      zinc_(Zcoplanarinside);
      return;
    }
  }else {
    bestfacet= facet;
    bestdist= *dist;
  }
  if (qh KEEPcoplanar + qh KEEPinside + qh KEEPnearinside) {
    oldfurthest= (pointT*)qh_setlast (bestfacet->coplanarset);
    if (oldfurthest) {
      zinc_(Zcomputefurthest);
      qh_distplane (oldfurthest, bestfacet, &dist2);
    }
    if (!oldfurthest || dist2 < bestdist) {
      qh_setappend(&bestfacet->coplanarset, point);
      if (bestdist > qh max_outside) {
	qh max_outside= bestdist;
	if (bestdist > qh TRACEdist)
	  istrace= True;
      }
    }else
      qh_setappend2ndlast(&bestfacet->coplanarset, point);
  }else { /* !KEEPcoplanar && !KEEPinside */
    if (bestdist > qh max_outside) {
      qh max_outside= bestdist;
      if (bestdist > qh TRACEdist) 
	istrace= True;
    }
  }
  if (istrace) {
    fprintf (qh ferr, "qh_partitioncoplanar: ====== p%d increases max_outside to %2.2g of f%d last p%d\n",
		   qh_pointid(point), bestdist, bestfacet->id, qh furthest_id);
    qh_errprint ("DISTANT", bestfacet, NULL, NULL, NULL);
  }
  trace4((qh ferr, "qh_partitioncoplanar: point p%d is coplanar with facet f%d (or inside) dist %2.2g\n",
	  qh_pointid(point), bestfacet->id, bestdist));
} /* partitioncoplanar */


/*-------------------------------------------------
-partitionpoint- assigns point to a visible facet
  uses qh_findbest with qh BESToutside and newfacets
  qh findbestnew if search new facets instead of findbest()
notes:
  after qh_distplane, this and qh_findbest are most expensive in 3-d
*/
void qh_partitionpoint (pointT *point, facetT *facet) {
  realT bestdist;
  pointT *oldfurthest;
  boolT isoutside;
  facetT *bestfacet;
  int numpart;
#if qh_COMPUTEfurthest
  realT dist;
#endif

  if (qh findbestnew)
    bestfacet= qh_findbestnew (point, facet,
			  &bestdist, &isoutside, &numpart);
  else
    bestfacet= qh_findbest (point, facet, qh BESToutside, True,
			  &bestdist, &isoutside, &numpart);
  zinc_(Ztotpartition);
  zzadd_(Zpartition, numpart);
  if (isoutside) {
    if (!bestfacet->outsideset 
    || !(oldfurthest= (pointT*)qh_setlast (bestfacet->outsideset))) {
      qh_setappend(&(bestfacet->outsideset), point);
      if (!bestfacet->newfacet) {
        qh_removefacet (bestfacet);  /* make sure it's after qh facet_next */
        qh_appendfacet (bestfacet);
      }
#if !qh_COMPUTEfurthest
      bestfacet->furthestdist= bestdist;
#endif
    }else {
#if qh_COMPUTEfurthest
      zinc_(Zcomputefurthest);
      qh_distplane (oldfurthest, bestfacet, &dist);
      if (dist < bestdist) 
	qh_setappend(&(bestfacet->outsideset), point);
      else
	qh_setappend2ndlast(&(bestfacet->outsideset), point);
#else
      if (bestfacet->furthestdist < bestdist) {
	qh_setappend(&(bestfacet->outsideset), point);
	bestfacet->furthestdist= bestdist;
      }else
	qh_setappend2ndlast(&(bestfacet->outsideset), point);
#endif
    }
    qh num_outside++;
    trace4((qh ferr, "qh_partitionpoint: point p%d is outside facet f%d\n",
	  qh_pointid(point), bestfacet->id));
  }else if (bestdist >= -qh MAXcoplanar) {
    zzinc_(Zcoplanarpart);
    if (qh KEEPcoplanar + qh KEEPnearinside || bestdist > qh max_outside) 
      qh_partitioncoplanar (point, bestfacet, &bestdist);
  }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 */

/*-------------------------------------------------
-partitionvisible- partitions points in visible_list to newfacet_list
  1st neighbor (if any) points to a horizon facet or a new facet
  repartitions coplanar points if allpoints (not used)
  qh findbestnew if search new facets instead of findbest()
  qh findbest_notsharp should be clear
*/
void qh_partitionvisible(/*visible_list*/ boolT allpoints, int *numoutside) {
  facetT *visible, *newfacet;
  pointT *point, **pointp;
  int coplanar=0, size, 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 (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)
          qh_partitionpoint (point, newfacet);
        else
          qh_partitioncoplanar (point, newfacet, NULL);
      }
    }
  }
  FOREACHvertex_(qh del_vertices) {
    if (vertex->point) {
      if (allpoints)
        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 */



/*----------------------------------------
-printsummary- prints the summary about the computation
  not in io.c so that user_eg.c can prevent io.c from loading
*/
void qh_printsummary(FILE *fp) {
  realT ratio, dist;
  float cpu;
  int size, id, total, numvertices, numcoplanars= 0;
  facetT *facet;

  size= qh num_points + qh_setsize (qh other_points);
  numvertices= qh num_vertices - qh_setsize (qh del_vertices);
  id= qh_pointid (qh GOODpointp);
  if (qh DELAUNAY) 
    numcoplanars= size - numvertices;
  else {
    FORALLfacets {
      if (facet->coplanarset)
	numcoplanars += qh_setsize( facet->coplanarset);
    }
  }
  if (id >=0 && qh STOPcone-1 != id && -qh STOPpoint-1 != id)
    size--;
  if (qh VORONOI)
    fprintf (fp, "\nVoronoi diagram by the convex");
  else if (qh DELAUNAY)
    fprintf (fp, "\nDelaunay triangulation by the convex");
  else if (qh HALFspace)
    fprintf (fp, "\nHalfspace intersection by the convex");
  else
    fprintf (fp, "\nConvex");
  fprintf(fp, " hull of %d points in %d-d:\n\n", size, qh hull_dim);
  fprintf(fp, "  Number of vertices%s: %d\n",
	  qh HALFspace ? " (halfspaces)" : "", numvertices);
  if (numcoplanars) 
    fprintf(fp, "  Number of %s points: %d\n", 
	    qh DELAUNAY ? "similar" : "coplanar", numcoplanars); 
  fprintf(fp, "  Number of facets%s: %d\n", 
	  qh HALFspace ? " (intersections)" : "",
	  qh num_facets - qh num_visible);
  if (qh num_good)
    fprintf(fp, "  Number of good facets: %d\n", qh num_good);
  fprintf(fp, "\nStatistics for: %s | %s", 
                      qh rbox_command, qh qhull_command);
  if (qh ROTATErandom > 0)
    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));
  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 (qh MERGING) {
    total= zzval_(Ztotmerge) - zzval_(Zcyclehorizon) + zzval_(Zcyclefacettot);
    fprintf(fp,"  Number of merged facets: %d\n", total);
    fprintf(fp,"  Number of distance tests for merging: %d\n",zzval_(Zbestdist)+
          zzval_(Zcentrumtests)+zzval_(Zdistconvex)+zzval_(Zdistcheck)+
          zzval_(Zdistzero));
  }
  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 totarea != 0.0) {
    fprintf(fp, "  %s facet area:   %2.8g\n", 
	    zzval_(Ztotmerge) ? "Approximate" : "Total", qh totarea);
    fprintf(fp, "  %s volume:       %2.8g\n", 
	    zzval_(Ztotmerge) ? "Approximate" : "Total", qh totvol);
  }
  if (!qh FORCEoutput && qh max_outside > qh DISTround) {
    dist= qh max_outside + qh DISTround;   /* agrees with qh_check_points */
    /* 1 DISTround to actual point */
    fprintf(fp, "  Maximum distance of %spoint above facet: %2.2g", 
	    (qh QHULLfinished ? "" : "merged "), dist);
    ratio= dist/(qh ONEmerge+ qh DISTround);
    if (qh MERGING && ratio > 0.05 && (qh ONEmerge > qh MINoutside))
      fprintf (fp, " (%.1fx)\n", ratio);
    else
      fprintf (fp, "\n");
  }
  if (!qh FORCEoutput && qh min_vertex < -qh DISTround) {
    dist= qh min_vertex - qh DISTround;   /* agrees with qh_check_points */
    /* 1 DISTround to actual point */
    fprintf(fp, "  Maximum distance of %svertex below facet: %2.2g", 
	    (qh QHULLfinished ? "" : "merged "), dist);
    ratio= -dist/(qh ONEmerge+qh DISTround);
    if (qh MERGING && ratio > 0.05) 
      fprintf (fp, " (%.1fx)\n", ratio);
    else
      fprintf (fp, "\n");
  }
  fprintf(fp, "\n");
} /* printsummary */