io.c

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

  if (string)
    fprintf (fp, string, facet->id);
  if (qh CENTERtype == qh_ASvoronoi) {
    num= qh hull_dim-1;
    if (!facet->normal || !facet->upperdelaunay || !qh ATinfinity) {
      if (!facet->center)
        facet->center= qh_facetcenter (facet->vertices);
      for (k=0; k<num; k++)
        fprintf (fp, qh_REAL_1, facet->center[k]);
    }else {
      for (k=0; k<num; k++)
        fprintf (fp, qh_REAL_1, qh_INFINITE);
    }
  }else /* qh CENTERtype == qh_AScentrum */ {
    num= qh hull_dim;
    if (!facet->center)
      facet->center= qh_getcentrum (facet);
    for (k=0; k<num; k++)
      fprintf (fp, qh_REAL_1, facet->center[k]);
  }
  if (format == qh_PRINTgeom && num == 2)
    fprintf (fp, " 0\n");
  else
    fprintf (fp, "\n");
} /* printcenter */

/*-----------------------------------------
-printcentrum- print centrum for a facet in OOGL format 2-d or 3-d only
  defines ->center if needed
*/
void qh_printcentrum (FILE *fp, facetT *facet, realT radius) {
  pointT *centrum, *projpt;
  boolT tempcentrum= False;
  realT xaxis[4], yaxis[4], normal[4], dist;
  realT green[3]={0, 1, 0};
  vertexT *apex;
  int k;
  
  if (qh CENTERtype == qh_AScentrum) {
    if (!facet->center)
      facet->center= qh_getcentrum (facet);
    centrum= facet->center;
  }else {
    centrum= qh_getcentrum (facet);
    tempcentrum= True;
  }
  fprintf (fp, "{appearance {-normal -edge normscale 0} ");
  if (qh firstcentrum) {
    qh firstcentrum= False;
    fprintf (fp, "{INST geom { define centrum CQUAD  # f%d\n\
-0.3 -0.3 0.0001     0 0 1 1\n\
 0.3 -0.3 0.0001     0 0 1 1\n\
 0.3  0.3 0.0001     0 0 1 1\n\
-0.3  0.3 0.0001     0 0 1 1 } transform { \n", facet->id);
  }else
    fprintf (fp, "{INST geom { : centrum } transform { # f%d\n", facet->id);
  apex= SETfirst_(facet->vertices);
  qh_distplane(apex->point, facet, &dist);
  projpt= qh_projectpoint(apex->point, facet, dist);
  for (k= qh hull_dim; k--; ) {
    xaxis[k]= projpt[k] - centrum[k];
    normal[k]= facet->normal[k];
  }
  if (qh hull_dim == 2) {
    xaxis[2]= 0;
    normal[2]= 0;
  }else if (qh hull_dim == 4) {
    qh_projectdim3 (xaxis, xaxis);
    qh_projectdim3 (normal, normal);
    qh_normalize2 (normal, qh PRINTdim, True, NULL, NULL);
  }
  qh_crossproduct (3, xaxis, normal, yaxis);
  fprintf (fp, "%8.4g %8.4g %8.4g 0\n", xaxis[0], xaxis[1], xaxis[2]);
  fprintf (fp, "%8.4g %8.4g %8.4g 0\n", yaxis[0], yaxis[1], yaxis[2]);
  fprintf (fp, "%8.4g %8.4g %8.4g 0\n", normal[0], normal[1], normal[2]);
  qh_printpoint3 (fp, centrum);
  fprintf (fp, "1 }}}\n"); 
  qh_memfree (projpt, qh normal_size);
  qh_printpointvect (fp, centrum, facet->normal, NULL, radius, green);
  if (tempcentrum)
    qh_memfree (centrum, qh normal_size);
} /* printcentrum */
  
/*-----------------------------------------
-printend- prints trailer for all output formats
*/
void qh_printend (FILE *fp, int format, facetT *facetlist, setT *facets, boolT printall) {
  int num;
  facetT *facet, **facetp;

  if (!qh printoutnum)
    fprintf (qh ferr, "qhull warning: no facets printed\n");
  switch (format) {
  case qh_PRINTgeom:
    if (qh hull_dim == 4 && qh DROPdim < 0  && !qh PRINTnoplanes) {
      qh visit_id++;
      num= 0;
      FORALLfacet_(facetlist)
        qh_printend4geom (fp, facet,&num, printall);
      FOREACHfacet_(facets) 
        qh_printend4geom (fp, facet, &num, printall);
      if (num != qh ridgeoutnum || qh printoutvar != qh ridgeoutnum) {
	fprintf (qh ferr, "qhull internal error (qh_printend): number of ridges %d != number printed %d and at end %d\n", qh ridgeoutnum, qh printoutvar, num);
	qh_errexit (qh_ERRqhull, NULL, NULL);
      }
    }else
      fprintf(fp, "}\n");
    break;
  case qh_PRINTinner:
  case qh_PRINTnormals:
  case qh_PRINTouter:
    if (qh CDDoutput) 
      fprintf (fp, "end\n");
    break;
  case qh_PRINTmathematica:
    fprintf(fp, "}\n");
    break;
  case qh_PRINTpoints:
    if (qh CDDoutput)
      fprintf (fp, "end\n");
    break;
  }
} /* printend */

/*------------------------------------------
-printend4geom- helper function for printbegin/printend
  just counts printed ridges if fp=NULL
  uses visitid
  must agree with printfacet4geom...
*/
void qh_printend4geom (FILE *fp, facetT *facet, int *nump, boolT printall) {
  realT color[3];
  int i, num= *nump;
  facetT *neighbor, **neighborp;
  ridgeT *ridge, **ridgep;
  
  if (!printall && qh_skipfacet(facet))
    return;
  if (qh PRINTnoplanes || (facet->visible && qh NEWfacets))
    return;
  if (!facet->normal)
    return;
  if (fp) {
    for (i=0; i<3; i++) {
      color[i]= (facet->normal[i]+1.0)/2.0;
      maximize_(color[i], -1.0);
      minimize_(color[i], +1.0);
    }
  }
  facet->visitid= qh visit_id;
  if (facet->simplicial) {
    FOREACHneighbor_(facet) {
      if (neighbor->visitid != qh visit_id) {
	if (fp)
          fprintf (fp, "3 %d %d %d %8.4g %8.4g %8.4g 1 # f%d f%d\n",
		 3*num, 3*num+1, 3*num+2, color[0], color[1], color[2],
		 facet->id, neighbor->id);
	num++;
      }
    }
  }else {
    FOREACHridge_(facet->ridges) {
      neighbor= otherfacet_(ridge, facet);
      if (neighbor->visitid != qh visit_id) {
	if (fp)
          fprintf (fp, "3 %d %d %d %8.4g %8.4g %8.4g 1 #r%d f%d f%d\n",
		 3*num, 3*num+1, 3*num+2, color[0], color[1], color[2],
		 ridge->id, facet->id, neighbor->id);
	num++;
      }
    }
  }
  *nump= num;
} /* printend4geom */

/*-----------------------------------------
-printfacet- prints all fields of a facet to fp
  ridges printed in neighbor order
*/
void qh_printfacet(FILE *fp, facetT *facet) {

  qh_printfacetheader (fp, facet);
  if (facet->ridges)
    qh_printfacetridges (fp, facet);
} /* printfacet */


/*----------------------------------------
-printfacet2geom- print facet as part of a 2-d VECT for Geomview
*/
void qh_printfacet2geom(FILE *fp, facetT *facet, realT color[3]) {
  pointT *point0, *point1;
  realT mindist, maxdist;
  int k;

  qh_facet2point (facet, &point0, &point1, &mindist);
  /*
    assume precise calculations in io.c with roundoff covered by qh_GEOMepsilon
    mindist is calculated within io.c.  maxoutside is calculated elsewhere
    so a DISTround error may have occured.
  */
  if (qh MERGING) {
    mindist -= qh PRINTradius;
#if qh_MAXoutside
    maxdist= facet->maxoutside+qh PRINTradius+ qh DISTround;
#else
    maxdist= qh max_outside+qh PRINTradius+ qh DISTround;
#endif
    if (qh PRINTcoplanar || qh PRINTspheres)
      maxdist += qh maxmaxcoord * qh_GEOMepsilon;
  }else 
    mindist= maxdist= 0;
  if (qh PRINTouter || (!qh PRINTnoplanes && !qh PRINTinner))
    qh_printfacet2geom_points(fp, point0, point1, facet, maxdist, color);
  if (qh PRINTinner || (!qh PRINTnoplanes && !qh PRINTouter &&
                maxdist - mindist > 2 * qh maxmaxcoord * qh_GEOMepsilon)) {
    for(k= 3; k--; )
      color[k]= 1.0 - color[k];
    qh_printfacet2geom_points(fp, point0, point1, facet, mindist, color);
  }
  qh_memfree (point1, qh normal_size);
  qh_memfree (point0, qh normal_size); 
} /* printfacet2geom */

/*-----------------------------------------
-printfacet2geom_points- prints a 2-d facet as a VECT with 2 points at
   some offset.   The points are on the facet's plane.
*/
void qh_printfacet2geom_points(FILE *fp, pointT *point1, pointT *point2,
			       facetT *facet, realT offset, realT color[3]) {
  pointT *p1= point1, *p2= point2;

  fprintf(fp, "VECT 1 2 1 2 1 # f%d\n", facet->id);
  if (offset != 0.0) {
    p1= qh_projectpoint (p1, facet, -offset);
    p2= qh_projectpoint (p2, facet, -offset);
  }
  fprintf(fp, "%8.4g %8.4g %8.4g\n%8.4g %8.4g %8.4g\n",
           p1[0], p1[1], 0.0, p2[0], p2[1], 0.0);
  if (offset != 0.0) {
    qh_memfree (p1, qh normal_size);
    qh_memfree (p2, qh normal_size);
  }
  fprintf(fp, "%8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]);
} /* printfacet2geom_points */



/*----------------------------------------
-printfacet2math- print 2-d Mathematica output for a facet
  may be non-simplicial
*/
void qh_printfacet2math(FILE *fp, facetT *facet, int notfirst) {
  pointT *point0, *point1;
  realT mindist;
  
  qh_facet2point (facet, &point0, &point1, &mindist);
  if (notfirst)
    fprintf(fp, ",");
  fprintf(fp, "Line[{{%10.8g, %10.8g}, {%10.8g, %10.8g}}]\n",
	  point0[0], point0[1], point1[0], point1[1]);
  qh_memfree (point1, qh normal_size);
  qh_memfree (point0, qh normal_size);
} /* printfacet2math */


/*-----------------------------------------
-printfacet3geom_nonsimplicial- print Geomview OFF for a 3-d nonsimplicial facet.
  if DOintersections, prints ridges to unvisited neighbors (qh visit_id) 
  uses visitid for intersections and ridges
*/
void qh_printfacet3geom_nonsimplicial(FILE *fp, facetT *facet, realT color[3]) {
  ridgeT *ridge, **ridgep;
  setT *projectedpoints, *vertices;
  vertexT *vertex, **vertexp, *vertexA, *vertexB;
  pointT *projpt, *point, **pointp;
  facetT *neighbor;
  realT dist, mindist= REALmax, maxdist;
  int cntvertices, k;
  realT black[3]={0, 0, 0}, green[3]={0, 1, 0};
 
  vertices= qh_facet3vertex (facet); /* oriented */
  cntvertices= qh_setsize(vertices);
  projectedpoints= qh_settemp(cntvertices);
  FOREACHvertex_(vertices) {
    zinc_(Zdistio);
    qh_distplane(vertex->point, facet, &dist);
    projpt= qh_projectpoint(vertex->point, facet, dist);
    qh_setappend (&projectedpoints, projpt);
    minimize_(mindist, dist);
  }
  /*
    assume precise calculations in io.c with roundoff covered by qh_GEOMepsilon
    mindist is calculated within io.c.  maxoutside is calculated elsewhere
    so a DISTround error may have occured.
  */
  if (qh MERGING) {
    mindist -= qh PRINTradius;
#if qh_MAXoutside
    maxdist= facet->maxoutside+qh PRINTradius+ qh DISTround;
#else
    maxdist= qh max_outside+qh PRINTradius+ qh DISTround;
#endif
    if (qh PRINTcoplanar || qh PRINTspheres)
      maxdist += qh maxmaxcoord * qh_GEOMepsilon;
  }else
    mindist= maxdist= 0;
  if (qh PRINTouter || (!qh PRINTnoplanes && !qh PRINTinner))
    qh_printfacet3geom_points(fp, projectedpoints, facet, maxdist, color);
  if (qh PRINTinner || (!qh PRINTnoplanes && !qh PRINTouter &&
                maxdist - mindist > 2 * qh maxmaxcoord * qh_GEOMepsilon)) {
    for (k=3; k--; )
      color[k]= 1.0 - color[k];
    qh_printfacet3geom_points(fp, projectedpoints, facet, mindist, color);
  }
  FOREACHpoint_(projectedpoints)
    qh_memfree (point, qh normal_size);
  qh_settempfree(&projectedpoints);
  qh_settempfree(&vertices);
  if ((qh DOintersections || qh PRINTridges)
  && (!facet->visible || !qh NEWfacets)) {
    facet->visitid= qh visit_id;
    FOREACHridge_(facet->ridges) {
      neighbor= otherfacet_(ridge, facet);
      if (neighbor->visitid != qh visit_id) {
        if (qh DOintersections)
          qh_printhyperplaneintersection(fp, facet, neighbor, ridge->vertices, black);
        if (qh PRINTridges) {
          vertexA= SETfirst_(ridge->vertices);
          vertexB= SETsecond_(ridge->vertices);
          qh_printline3geom (fp, vertexA->point, vertexB->point, green);
        }
      }
    }
  }
} /* printfacet3geom_nonsimplicial */

/*-----------------------------------------
-printfacet3geom_points- prints a 3-d facet as OFF Geomview object. 
  Facet is determined as a list of points
*/
void qh_printfacet3geom_points(FILE *fp, setT *points, facetT *facet, realT offset, realT color[3]) {
  int k, n= qh_setsize(points), i;
  pointT *point, **pointp;
  setT *printpoints;

  fprintf(fp, "{ OFF %d 1 1 # f%d\n", n, facet->id);
  if (offset != 0.0) {
    printpoints= qh_settemp (n);
    FOREACHpoint_(points) 
      qh_setappend (&printpoints, qh_projectpoint(point, facet, -offset));
  }else
    printpoints= points;
  FOREACHpoint_(printpoints) {
    for (k=0; k<qh hull_dim; k++) {
      if (k == qh DROPdim)
        fprintf(fp, "0 ");
      else
        fprintf(fp, "%8.4g ", point[k]);
    }
    if (printpoints != points)
      qh_memfree (point, qh normal_size);
    fprintf (fp, "\n");
  }
  if (printpoints != points)
    qh_settempfree (&printpoints);
  fprintf(fp, "%d ", n);
  for(i= 0; i < n; i++)
    fprintf(fp, "%d ", i);
  fprintf(fp, "%8.4g %8.4g %8.4g 1.0 }\n", color[0], color[1], color[2]);
} /* printfacet3geom_points */


/*-----------------------------------------
-printfacet3geom_simplicial- print Geomview OFF for a 3-d simplicial facet.
  may flip color
  uses visitid for intersections and ridges
*/
void qh_printfacet3geom_simplicial(FILE *fp, facetT *facet, realT color[3]) {
  setT *points, *vertices;
  vertexT *vertex, **vertexp, *vertexA, *vertexB;
  facetT *neighbor, **neighborp;
  realT maxdist, mindist;
  realT black[3]={0, 0, 0}, green[3]={0, 1, 0};
  int k;

  vertices= qh_facet3vertex (facet);
  points= qh_settemp (qh TEMPsize);
  FOREACHvertex_(vertices)
    qh_setappend(&points, vertex->point);