io.c

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

    qh_printfacetheader (fp, facet);
    break;
  case qh_PRINTgeom:  /* either 2 , 3, or 4-d by qh_printbegin */
    if (!facet->normal)
      break;
    for (k= qh hull_dim; k--; ) {
      color[k]= (facet->normal[k]+1.0)/2.0;
      maximize_(color[k], -1.0);
      minimize_(color[k], +1.0);
    }
    qh_projectdim3 (color, color);
    if (qh PRINTdim != qh hull_dim)
      qh_normalize2 (color, 3, True, NULL, NULL);
    if (qh hull_dim <= 2)
      qh_printfacet2geom (fp, facet, color);
    else if (qh hull_dim == 3) {
      if (facet->simplicial)
        qh_printfacet3geom_simplicial (fp, facet, color);
      else
        qh_printfacet3geom_nonsimplicial (fp, facet, color);
    }else {
      if (facet->simplicial)
        qh_printfacet4geom_simplicial (fp, facet, color);
      else
        qh_printfacet4geom_nonsimplicial (fp, facet, color);
    }
    break;
  case qh_PRINTids:
    fprintf (fp, "%d\n", facet->id);
    break;
  case qh_PRINToff:
  case qh_PRINTincidences:
    if (qh hull_dim == 3) 
      qh_printfacet3vertex (fp, facet, format);
    else if (facet->simplicial || qh hull_dim == 2 || format == qh_PRINToff)
      qh_printfacetNvertex_simplicial (fp, facet, format);
    else
      qh_printfacetNvertex_nonsimplicial (fp, facet, qh printoutvar++);
    break;
  case qh_PRINTinner:
    mindist= REALmax;
    FOREACHvertex_(facet->vertices) {
      qh_distplane (vertex->point, facet, &dist);
      minimize_(mindist, dist);
    }
    offset= facet->offset - mindist - qh DISTround;
    goto LABELprintnorm;
  case qh_PRINTmerges:
    fprintf (fp, "%d\n", facet->nummerge);
    break;
  case qh_PRINTnormals:
    offset= facet->offset;
    goto LABELprintnorm;
  case qh_PRINTouter:
#if qh_MAXoutside
    offset= -facet->maxoutside;
#endif
    if (!qh MERGING || !qh_MAXoutside || qh SKIPcheckmax)
      offset= -fmax_(qh max_outside, qh DISTround);
    offset -= qh DISTround;      /* agrees with qh_check_points */
    /* 1 DISTround to actual point */
    offset += facet->offset;
  LABELprintnorm:
    if (!facet->normal) {
      fprintf (fp, "no normal for facet f%d\n", facet->id);
      break;
    }
    if (qh CDDoutput) 
      fprintf (fp, qh_REAL_1, -offset);
    for (k=0; k<qh hull_dim; k++) 
      fprintf (fp, qh_REAL_1, facet->normal[k]);
    if (!qh CDDoutput) 
      fprintf (fp, qh_REAL_1, offset);
    fprintf (fp, "\n");
    break;
  case qh_PRINTmathematica:  /* either 2 or 3-d by qh_printbegin */
    if (qh hull_dim == 2)
      qh_printfacet2math (fp, facet, qh printoutvar++);
    else 
      qh_printfacet3math (fp, facet, qh printoutvar++);
    break;
  case qh_PRINTneighbors:
    fprintf (fp, "%d", qh_setsize (facet->neighbors));
    FOREACHneighbor_(facet)
      fprintf (fp, " %d", 
	       neighbor->visitid ? neighbor->visitid - 1: - neighbor->id);
    fprintf (fp, "\n");
    break;
  case qh_PRINTpointintersect:
    if (!qh feasible_point) {
      fprintf (fp, "qhull input error (qh_printafacet): option 'Fp' needs qh feasible_point\n");
      qh_errexit( qh_ERRinput, NULL, NULL);
    }
    if (facet->offset > 0)
      goto LABELprintinfinite;
    point= coordp= (coordT*)qh_memalloc (qh normal_size);
    normp= facet->normal;
    feasiblep= qh feasible_point;
    if (facet->offset < -qh MINdenom) {
      for (k= qh hull_dim; k--; )
        *(coordp++)= (*(normp++) / - facet->offset) + *(feasiblep++);
    }else {
      for (k= qh hull_dim; k--; ) {
        *(coordp++)= qh_divzero (*(normp++), facet->offset, qh MINdenom_1,
				 &zerodiv) + *(feasiblep++);
        if (zerodiv) {
          qh_memfree (point, qh normal_size);
          goto LABELprintinfinite;
        }
      }
    }
    qh_printpoint (fp, NULL, point);
    qh_memfree (point, qh normal_size);
    break;
  LABELprintinfinite:
    for (k= qh hull_dim; k--; )
      fprintf (fp, qh_REAL_1, qh_INFINITE);
    fprintf (fp, "\n");   
    break;
  case qh_PRINTpointnearest:
    FOREACHpoint_(facet->coplanarset) {
      int id, id2;
      vertex= qh_nearvertex (facet, point, &dist);
      id= qh_pointid (vertex->point);
      id2= qh_pointid (point);
      fprintf (fp, "%d %d %d " qh_REAL_1 "\n", id, id2, facet->id, dist);
    }
    break;
  case qh_PRINTpoints:  /* VORONOI only by qh_printbegin */
    if (qh CDDoutput)
      fprintf (fp, "1 ");
    qh_printcenter (fp, format, NULL, facet);
    break;
  case qh_PRINTvertices:
    fprintf (fp, "%d", qh_setsize (facet->vertices));
    FOREACHvertex_(facet->vertices)
      fprintf (fp, " %d", qh_pointid (vertex->point));
    fprintf (fp, "\n");
    break;
  }
} /* printafacet */

/*-----------------------------------------
-printbegin- prints header for all output formats
  checks for valid format
  uses qh visit_id for 3/4off
  changes qh interior_point if printing centrums
*/
void qh_printbegin (FILE *fp, int format, facetT *facetlist, setT *facets, boolT printall) {
  int numfacets, numsimplicial, numridges, totneighbors, numcoplanars;
  int i, num;
  facetT *facet, **facetp;
  vertexT *vertex, **vertexp;
  setT *vertices;
  pointT *point, **pointp, *pointtemp;

  qh printoutnum= 0;
  qh_countfacets (facetlist, facets, printall, &numfacets, &numsimplicial, 
      &totneighbors, &numridges, &numcoplanars);
  switch (format) {
  case qh_PRINTnone:
    break;
  case qh_PRINTarea:
    fprintf (fp, "%d\n", numfacets);
    break;
  case qh_PRINTcoplanars:
    fprintf (fp, "%d\n", numfacets);
    break;
  case qh_PRINTcentrums:
    if (qh CENTERtype == qh_ASnone)
      qh_clearcenters (qh_AScentrum);
    fprintf (fp, "%d\n%d\n", qh hull_dim, numfacets);
    break;
  case qh_PRINTfacets:
  case qh_PRINTfacets_xridge:
    if (facetlist)
      qh_printvertexlist (fp, "Vertices and facets:\n", facetlist, facets, printall);
    break;
  case qh_PRINTgeom: 
    if (qh hull_dim > 4)  /* qh_initqhull_globals also checks */
      goto LABELnoformat;
    if (qh VORONOI && qh hull_dim > 3)  /* PRINTdim == DROPdim == hull_dim-1 */
      goto LABELnoformat;
    if (qh hull_dim == 2 && (qh PRINTridges || qh DOintersections))
      fprintf (qh ferr, "qhull warning: output for ridges and intersections not implemented in 2-d\n");
    if (qh hull_dim == 4 && (qh PRINTinner || qh PRINTouter ||
			     (qh PRINTdim == 4 && qh PRINTcentrums)))
      fprintf (qh ferr, "qhull warning: output for outer/inner planes and centrums not implemented in 4-d\n");
    if (qh PRINTdim == 4 && (qh PRINTspheres))
      fprintf (qh ferr, "qhull warning: output for vertices not implemented in 4-d\n");
    if (qh PRINTdim == 4 && qh DOintersections && qh PRINTnoplanes)
      fprintf (qh ferr, "qhull warning: 'Gnh' generates no output in 4-d\n");
    if (qh PRINTdim == 2) {
      fprintf(fp, "{appearance {linewidth 3} LIST # %s | %s\n",
	      qh rbox_command, qh qhull_command);
    }else if (qh PRINTdim == 3) {
      fprintf(fp, "{appearance {+edge -evert linewidth 2} LIST # %s | %s\n",
	      qh rbox_command, qh qhull_command);
    }else if (qh PRINTdim == 4) {
      qh visit_id++;
      num= 0;
      FORALLfacet_(facetlist)    /* get number of ridges to be printed */
        qh_printend4geom (NULL, facet, &num, printall);
      FOREACHfacet_(facets)
        qh_printend4geom (NULL, facet, &num, printall);
      qh ridgeoutnum= num;
      qh printoutvar= 0;  /* counts number of ridges in output */
      fprintf (fp, "LIST # %s | %s\n", qh rbox_command, qh qhull_command);
    }
    if (qh PRINTdots) {
      qh printoutnum++;
      num= qh num_points + qh_setsize (qh other_points);
      if (qh DELAUNAY && qh ATinfinity)
	num--;
      if (qh PRINTdim == 4)
        fprintf (fp, "4VECT %d %d 1\n", num, num);
      else
	fprintf (fp, "VECT %d %d 1\n", num, num);
      for (i= num; i--; ) {
        if (i % 20 == 0)
          fprintf (fp, "\n");
	fprintf (fp, "1 ");
      }
      fprintf (fp, "# 1 point per line\n1 ");
      for (i= num-1; i--; ) {
        if (i % 20 == 0)
          fprintf (fp, "\n");
	fprintf (fp, "0 ");
      }
      fprintf (fp, "# 1 color for all\n");
      FORALLpoints {
        if (!qh DELAUNAY || !qh ATinfinity || qh_pointid(point) != qh num_points-1) {
	  if (qh PRINTdim == 4)
	    qh_printpoint (fp, NULL, point);
	  else
	    qh_printpoint3 (fp, point);
	}
      }
      FOREACHpoint_(qh other_points) {
	if (qh PRINTdim == 4)
	  qh_printpoint (fp, NULL, point);
	else
	  qh_printpoint3 (fp, point);
      }
      fprintf (fp, "0 1 1 1  # color of points\n");
    }
    if (qh PRINTdim == 4  && !qh PRINTnoplanes)
      /* 4dview loads up multiple 4OFF objects slowly */
      fprintf(fp, "4OFF %d %d 1\n", 3*qh ridgeoutnum, qh ridgeoutnum);
    qh PRINTcradius= 2 * qh DISTround;  /* include test DISTround */
    if (qh PREmerge) {
      maximize_(qh PRINTcradius, qh premerge_centrum + qh DISTround);
    }else if (qh POSTmerge)
      maximize_(qh PRINTcradius, qh postmerge_centrum + qh DISTround);
    qh PRINTradius= qh PRINTcradius;
    if (qh PRINTspheres + qh PRINTcoplanar)
      maximize_(qh PRINTradius, qh maxmaxcoord * qh_MINradius);
    if (qh premerge_cos < REALmax/2) {
      maximize_(qh PRINTradius, (1- qh premerge_cos) * qh maxmaxcoord);
    }else if (!qh PREmerge && qh POSTmerge && qh postmerge_cos < REALmax/2) {
      maximize_(qh PRINTradius, (1- qh postmerge_cos) * qh maxmaxcoord);
    }
    maximize_(qh PRINTradius, qh MINvisible); 
    if (qh PRINTdim != 4 &&
	(qh PRINTcoplanar || qh PRINTspheres || qh PRINTcentrums)) {
      vertices= qh_facetvertices (facetlist, facets, printall);
      if (qh PRINTspheres && qh PRINTdim <= 3)
         qh_printspheres (fp, vertices, qh PRINTradius);
      if (qh PRINTcoplanar || qh PRINTcentrums) {
        qh firstcentrum= True;
        if (qh PRINTcoplanar&& !qh PRINTspheres) {
          FOREACHvertex_(vertices) 
            qh_printpointvect2 (fp, vertex->point, NULL,
				qh interior_point, qh PRINTradius);
	}
        FORALLfacet_(facetlist) {
	  if (!printall && qh_skipfacet(facet))
	    continue;
	  if (!facet->normal)
	    continue;
          if (qh PRINTcentrums && qh PRINTdim <= 3)
            qh_printcentrum (fp, facet, qh PRINTcradius);
          FOREACHpoint_(facet->coplanarset)
            qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius);
          FOREACHpoint_(facet->outsideset)
            qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius);
        }
        FOREACHfacet_(facets) {
	  if (!printall && qh_skipfacet(facet))
	    continue;
	  if (!facet->normal)
	    continue;
          if (qh PRINTcentrums && qh PRINTdim <= 3)
            qh_printcentrum (fp, facet, qh PRINTcradius);
          FOREACHpoint_(facet->coplanarset)
            qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius);
          FOREACHpoint_(facet->outsideset)
            qh_printpointvect2 (fp, point, facet->normal, NULL, qh PRINTradius);
        }
      }
      qh_settempfree (&vertices);
    }
    qh visit_id++; /* for printing hyperplane intersections */
    break;
  case qh_PRINTids:
    fprintf (fp, "%d\n", numfacets);
    break;
  case qh_PRINTincidences:
    if (qh VORONOI)
      fprintf (qh ferr, "qhull warning: writing Delaunay.  Use 'p' or 'o' for Voronoi centers\n");
    qh printoutvar= qh vertex_id;  /* centrum id for non-simplicial facets */
    if (qh hull_dim <= 3)
      fprintf(fp, "%d\n", numfacets);
    else
      fprintf(fp, "%d\n", numsimplicial+numridges);
    break;
  case qh_PRINTinner:
  case qh_PRINTnormals:
  case qh_PRINTouter:
    if (qh CDDoutput)
      fprintf (fp, "%s | %s\nbegin\n    %d %d real\n", qh rbox_command, 
              qh qhull_command, numfacets, qh hull_dim+1);
    else
      fprintf (fp, "%d\n%d\n", qh hull_dim+1, numfacets);
    break;
  case qh_PRINTmathematica:  
    if (qh hull_dim > 3)  /* qh_initbuffers also checks */
      goto LABELnoformat;
    if (qh VORONOI)
      fprintf (qh ferr, "qhull warning: output is the Delaunay triangulation\n");
    fprintf(fp, "{\n");
    qh printoutvar= 0;   /* counts number of facets for notfirst */
    break;
  case qh_PRINTmerges:
    fprintf (fp, "%d\n", numfacets);
    break;
  case qh_PRINTpointintersect:
    fprintf (fp, "%d\n%d\n", qh hull_dim, numfacets);
    break;
  case qh_PRINTneighbors:
    fprintf (fp, "%d\n", numfacets);
    break;
  case qh_PRINToff:
    if (qh VORONOI)
      goto LABELnoformat;
    fprintf (fp, "%d\n%d %d %d\n", qh hull_dim,
      qh num_points+qh_setsize (qh other_points), numfacets, totneighbors/2);
    FORALLpoints
      qh_printpoint (qh fout, NULL, point);
    FOREACHpoint_(qh other_points)
      qh_printpoint (qh fout, NULL, point);
    break;
  case qh_PRINTpointnearest:
    fprintf (fp, "%d\n", numcoplanars);
    break;
  case qh_PRINTpoints:
    if (!qh VORONOI)
      goto LABELnoformat;
    if (qh CDDoutput)
      fprintf (fp, "%s | %s\nbegin\n%d %d real\n", qh rbox_command,
             qh qhull_command, numfacets, qh hull_dim);
    else
      fprintf (fp, "%d\n%d\n", qh hull_dim-1, numfacets);
    break;
  case qh_PRINTvertices:
    fprintf (fp, "%d\n", numfacets);
    break;
  case qh_PRINTsummary:
  default:
  LABELnoformat:
    fprintf (qh ferr, "qhull internal error (qh_printbegin): can not use this format for dimension %d\n",
         qh hull_dim);
    qh_errexit (qh_ERRqhull, NULL, NULL);
  }
} /* printbegin */

/*-----------------------------------------
-printcenter- print facet center as either centrum or Voronoi center
  nop if qh CENTERtype neither CENTERvoronoi nor CENTERcentrum
  if upper envelope of Delaunay triangulation and point at-infinity
    prints qh_INFINITE instead;
  defines facet->center if needed
  string may be NULL or include %d for id.  Don't include other '%' codes.
  if format=PRINTgeom, adds a 0 if otherwise 2-d
*/
void qh_printcenter (FILE *fp, int format, char *string, facetT *facet) {
  int k, num;
  
  if (qh CENTERtype != qh_ASvoronoi && qh CENTERtype != qh_AScentrum)
    return;