io.c

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

	ridge= qh_nextridge3d (ridge, facet, NULL);
	}
    }else {
      FOREACHneighbor_(facet) {
	FOREACHridge_(facet->ridges) {
	  if (otherfacet_(ridge,facet) == neighbor) {
	    ridge->seen= True;
	    qh_printridge(fp, ridge);
	    numridges++;
	  }
	}
      }
    }
    if (numridges != qh_setsize (facet->ridges)) {
      fprintf (fp, "     - all ridges:");
      FOREACHridge_(facet->ridges) 
	fprintf (fp, " r%d", ridge->id);
        fprintf (fp, "\n");
    }
    FOREACHridge_(facet->ridges) {
      if (!ridge->seen) 
	qh_printridge(fp, ridge);
    }
  }
} /* printfacetridges */

/*-----------------------------------------
-printfacets- prints facetlist and/or facet set in output format
*/
void qh_printfacets(FILE *fp, int format, facetT *facetlist, setT *facets, boolT printall) {
  int numfacets, numsimplicial, numridges, totneighbors, numcoplanars;
  facetT *facet, **facetp;
  setT *vertices;
  coordT *center;

  qh old_randomdist= qh RANDOMdist;
  qh RANDOMdist= False;
  if (format == qh_PRINTnone)
    ; /* print nothing */
  else if (format == qh_PRINTaverage) {
    vertices= qh_facetvertices (facetlist, facets, printall);
    center= qh_getcenter (vertices);
    fprintf (fp, "%d 1\n", qh hull_dim);
    qh_printpoint (fp, NULL, center);
    qh_memfree (center, qh normal_size);
    qh_settempfree (&vertices);
  }else if (format == qh_PRINToptions)
    fprintf(fp, "Options selected for qhull %s:\n%s\n", qh_version, qh qhull_options);
  else if (format == qh_PRINTpoints && !qh VORONOI)
    qh_printpoints_out (fp, facetlist, facets, printall);
  else if (format == qh_PRINTqhull)
    fprintf (fp, "%s | %s\n", qh rbox_command, qh qhull_command);
  else if (format == qh_PRINTsize) {
    fprintf (fp, "0\n2 ");
    fprintf (fp, qh_REAL_1, qh totarea);
    fprintf (fp, qh_REAL_1, qh totvol);
    fprintf (fp, "\n");
  }else if (format == qh_PRINTsummary) {
    qh_countfacets (facetlist, facets, printall, &numfacets, &numsimplicial, 
      &totneighbors, &numridges, &numcoplanars);
    vertices= qh_facetvertices (facetlist, facets, printall); 
    fprintf (fp, "7 %d %d %d %d %d %d %d\n2 ", qh hull_dim, 
                qh num_points + qh_setsize (qh other_points),
                qh num_vertices, qh num_facets - qh num_visible,
                qh_setsize (vertices), numfacets, numcoplanars);
    fprintf (fp, qh_REAL_2n, qh max_outside + qh DISTround, 
                  qh min_vertex - qh DISTround);  /* agrees with qh_check_points */
    qh_settempfree (&vertices);
  }else if (format == qh_PRINTvneighbors)
    qh_printvneighbors (fp, facetlist, facets, printall);
  else if (qh VORONOI && format == qh_PRINToff)
    qh_printvoronoi (fp, format, facetlist, facets, printall);
  else if (qh VORONOI && format == qh_PRINTgeom) {
    qh_printbegin (fp, format, facetlist, facets, printall);
    qh_printvoronoi (fp, format, facetlist, facets, printall);
    qh_printend (fp, format, facetlist, facets, printall);
  }else {
    qh_printbegin (fp, format, facetlist, facets, printall);
    FORALLfacet_(facetlist)
      qh_printafacet (fp, format, facet, printall);
    FOREACHfacet_(facets) 
      qh_printafacet (fp, format, facet, printall);
    qh_printend (fp, format, facetlist, facets, printall);
  }
  qh RANDOMdist= qh old_randomdist;
} /* printfacets */


/*-----------------------------------------
-printhelp_degenerate- prints descriptive message
  no message if qh_QUICKhelp
*/
void qh_printhelp_degenerate(FILE *fp) {
  
  if (qh MERGEexact || qh PREmerge) 
    fprintf(fp, "\n\
A Qhull error has occurred.  It should have corrected the above\n\
precision error.  Please report the error to [e-mail]\n");
  else if (!qh_QUICKhelp) {
    fprintf(fp, "\n\
Precision problems were detected during construction of the convex hull.\n\
This occurs because convex hull algorithms assume that calculations are\n\
exact, but floating-point arithmetic has round-off errors.\n\
\n\
To correct for precision problems, use option 'C-0' (2-d to 4-d) or \n\
'Qx' (5-d and higher).  This will merge non-convex facets.  See \"Imprecision\n\
in Qhull\" (qh-impre.html).\n\
\n\
If you do not use options 'C-n', 'A-n' or 'Qx', the output may include\n\
coplanar ridges, concave ridges, and flipped facets.  In 4-d and higher,\n\
Qhull may produce a ridge with four neighbors or two facets with the same \n\
vertices.  Qhull reports these events when they occur.  It stops when a\n\
concave ridge, flipped facet, or duplicate facet occurs.\n");
#if REALfloat
    fprintf (fp, "\
\n\
Qhull is currently using single precision arithmetic.  The following\n\
will probably remove the precision problems:\n\
  - recompile qhull for double precision (#define REALfloat 0 in user.h).\n");
#endif
    fprintf(fp, "\
\n\
If you can not merge facets because you want simplicial facets,\n\
try one or more of the following options.  They can not guarantee an output.\n\
  - use 'Po' to produce output and prevent partitioning for flipped facets\n\
  - use 'V0' to set min. distance to visible facet as 0 instead of roundoff\n\
  - use 'En' to specify a maximum roundoff error less than %2.2g.\n\
  - options 'Qf', 'QR0', and 'QbB' may also help\n",
               qh DISTround);
    if (!qh ATinfinity && qh DELAUNAY)
      fprintf(fp, "\
  - remove 'Qu' to add a point \"at-infinity\"\n");
    fprintf(fp, "\
\n\
To guarantee simplicial output:\n\
  - use exact arithmetic (see \"Imprecision in Qhull\", qh-impre.html)\n\
\n\
To merge non-convex facets\n");
    if (qh hull_dim >= 5) 
      fprintf (fp, "\
  - use option 'Qx'\n");
    else
      fprintf (fp, "\
  - use option 'C-0'\n");
  }
} /* printhelp_degenerate */


/*-----------------------------------------
-printhelp_singular- prints descriptive message
  qh_QUICKhelp just prints the numbers.
*/
void qh_printhelp_singular(FILE *fp) {
  facetT *facet;
  vertexT *vertex, **vertexp;
  realT min, max, *coord, dist;
  int i,k;
  
  fprintf(fp, "\n\
The input to qhull appears to be less than %d dimensional, or a\n\
computation has overflowed.\n\n\
Qhull could not construct a clearly convex simplex from points:\n",
           qh hull_dim);
  qh_printvertexlist (fp, "", qh facet_list, NULL, qh_ALL);
  if (!qh_QUICKhelp)
    fprintf(fp, "\n\
The center point is coplanar with a facet, or a vertex is coplanar\n\
with a neighboring facet.  The maximum round off error for\n\
computing distances is %2.2g.  The center point, facets and distances\n\
to the center point are as follows:\n\n", qh DISTround);
  qh_printpoint (fp, "center point", qh interior_point);
  fprintf (fp, "\n");
  FORALLfacets {
    fprintf (fp, "facet");
    FOREACHvertex_(facet->vertices)
      fprintf (fp, " p%d", qh_pointid(vertex->point));
    zinc_(Zdistio);
    qh_distplane(qh interior_point, facet, &dist);
    fprintf (fp, " distance= %4.2g\n", dist);
  }
  if (!qh_QUICKhelp) {
    if (qh HALFspace) 
      fprintf (fp, "\n\
These points are the dual of the given halfspaces.  They indicate that\n\
the intersection is degenerate.\n");
    fprintf (fp,"\n\
These points either have a maximum or minimum x-coordinate, or\n\
they maximize the determinant for k coordinates.  Trial points\n\
are first selected from points that maximize a coordinate.\n");
    if (qh hull_dim >= qh_INITIALmax)
      fprintf (fp, "\n\
Because of the high dimension, the min x-coordinate and max-coordinate\n\
points are used if the determinant is non-zero.  The option 'Qs' will\n\
do a better, though much slower, job.  Instead of 'Qs', you can change\n\
the points by randomly rotating the input with 'QR0'.\n");
  }
  fprintf (fp, "\nThe min and max coordinates for each dimension are:\n");
  for (k=0; k<qh hull_dim; k++) {
    min= REALmax;
    max= -REALmin;
    for (i=qh num_points, coord= qh first_point+k; i--; coord += qh hull_dim) {
      maximize_(max, *coord);
      minimize_(min, *coord);
    }
    fprintf (fp, "  %d:  %8.4g  %8.4g  difference= %4.4g\n", k, min, max, max-min);
  }
  if (!qh_QUICKhelp) {
    fprintf (fp, "\n\
If the input should be full dimensional, you have several options that\n\
may determine an initial simplex:\n\
  - use 'QbB' to scale the points to the unit cube\n\
  - use 'QR0' to randomly rotate the input for different maximum points\n\
  - use the 'Qs' flag to search all points for the initial simplex\n\
  - use 'En' to specify a maximum roundoff error less than %2.2g.\n\
  - trace execution with 'T3' to see the determinant for each point.\n",
                     qh DISTround);
#if REALfloat
    fprintf (fp, "\
  - recompile qhull for double precision (#define REALfloat 0 in qhull.h).\n");
#endif
    fprintf (fp, "\n\
If the input is lower dimensional:\n\
   - use 'Qbk:0Bk:0' to delete coordinate k from the input.  You should\n\
     pick the coordinate with the least range.  The hull will have the\n\
     correct topology.\n\
   - determine the flat containing the points, rotate the points\n\
     into a coordinate plane, and delete the other coordinates.\n\
   - add one or more points to make the input full dimensional.\n");
  }
} /* printhelp_singular */

/*-----------------------------------------
-printhyperplaneintersection- print Geomview OFF or 4OFF for
   the intersection of two hyperplanes in 3-d or 4-d
*/
void qh_printhyperplaneintersection(FILE *fp, facetT *facet1, facetT *facet2,
		   setT *vertices, realT color[3]) {
  realT d1, d2, costheta, denominator, dist1, dist2, s, t, mindenom, p[4];
  vertexT *vertex, **vertexp;
  int i, k;
  boolT nearzero1, nearzero2;
  
  costheta= qh_getangle(facet1->normal, facet2->normal);
  denominator= 1 - costheta * costheta;
  d1= -facet1->offset;
  d2= -facet2->offset;
  i= qh_setsize(vertices);
  if (qh hull_dim == 3)
    fprintf(fp, "VECT 1 %d 1 %d 1 ", i, i);
  else if (qh hull_dim == 4 && qh DROPdim >= 0)
    fprintf(fp, "OFF 3 1 1 ");
  else
    qh printoutvar++;
  fprintf (fp, "# intersect f%d f%d\n", facet1->id, facet2->id);
  mindenom= 1 / (10.0 * qh maxmaxcoord);
  FOREACHvertex_(vertices) {
    zadd_(Zdistio, 2);
    qh_distplane(vertex->point, facet1, &dist1);
    qh_distplane(vertex->point, facet2, &dist2);
    s= qh_divzero (-dist1 + costheta * dist2, denominator,mindenom,&nearzero1);
    t= qh_divzero (-dist2 + costheta * dist1, denominator,mindenom,&nearzero2);
    if (nearzero1 || nearzero2)
      s= t= 0.0;
    for(k= qh hull_dim; k--; )
      p[k]= vertex->point[k] + facet1->normal[k] * s + facet2->normal[k] * t;
    if (qh PRINTdim <= 3) {
      qh_projectdim3 (p, p);
      fprintf(fp, "%8.4g %8.4g %8.4g # ", p[0], p[1], p[2]);
    }else 
      fprintf(fp, "%8.4g %8.4g %8.4g %8.4g # ", p[0], p[1], p[2], p[3]);
    if (nearzero1+nearzero2)
      fprintf (fp, "p%d (coplanar facets)\n", qh_pointid (vertex->point));
    else
      fprintf (fp, "projected p%d\n", qh_pointid (vertex->point));
  }
  if (qh hull_dim == 3)
    fprintf(fp, "%8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]); 
  else if (qh hull_dim == 4 && qh DROPdim >= 0)  
    fprintf(fp, "3 0 1 2 %8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]);
} /* printhyperplaneintersection */

/*----------------------------------------
-printline3geom- prints a line as a VECT
  0's for DROPdim
  if pointA == pointB, it's a 1 point VECT
*/
void qh_printline3geom (FILE *fp, pointT *pointA, pointT *pointB, realT color[3]) {
  int k;
  realT pA[4], pB[4];

  qh_projectdim3(pointA, pA);
  qh_projectdim3(pointB, pB);
  if ((fabs(pA[0] - pB[0]) > 1e-3) || 
      (fabs(pA[1] - pB[1]) > 1e-3) || 
      (fabs(pA[2] - pB[2]) > 1e-3)) {
    fprintf (fp, "VECT 1 2 1 2 1\n");
    for (k= 0; k < 3; k++)
       fprintf (fp, "%8.4g ", pB[k]);
    fprintf (fp, " # p%d\n", qh_pointid (pointB));
  }else
    fprintf (fp, "VECT 1 1 1 1 1\n");
  for (k=0; k< 3; k++)
    fprintf (fp, "%8.4g ", pA[k]);
  fprintf (fp, " # p%d\n", qh_pointid (pointA));
  fprintf (fp, "%8.4g %8.4g %8.4g 1\n", color[0], color[1], color[2]);
}

/*----------------------------------------
-printneighborhood- print neighborhood of one or two facets
  calls findgood_all if active and !printall
  bumps qh visit_id
*/
void qh_printneighborhood (FILE *fp, int format, facetT *facetA, facetT *facetB, boolT printall) {
  facetT *neighbor, **neighborp, *facet;
  setT *facets;

  if (format == qh_PRINTnone)
    return;
  qh_findgood_all (qh facet_list);
  if (facetA == facetB)
    facetB= NULL;
  facets= qh_settemp (2*(qh_setsize (facetA->neighbors)+1));
  qh visit_id++;
  for (facet= facetA; facet; facet= ((facet == facetA) ? facetB : NULL)) {
    if (facet->visitid != qh visit_id) {
      facet->visitid= qh visit_id;
      qh_setappend (&facets, facet);
    }
    FOREACHneighbor_(facet) {
      if (neighbor->visitid == qh visit_id)
        continue;
      neighbor->visitid= qh visit_id;
      if (printall || !qh_skipfacet (neighbor))
        qh_setappend (&facets, neighbor);
    }
  }
  qh_printfacets (fp, format, NULL, facets, printall);
  qh_settempfree (&facets);
} /* printneighborhood */

/*----------------------------------------
-printpoint- prints the coordinates of a point
  prints point id if string and valid point
  nop if point is NULL
*/
void qh_printpoint(FILE *fp, char *string, pointT *point) {
  int k,id;
  realT r; /*bug fix*/
  
  if (!point)
    return;
  if (string) {
    fputs (string, fp);
    if ((id= qh_pointid(point)) != -1)
      fprintf(fp, " p%d: ", id);
  }
  for(k= qh hull_dim; k--; ) {
    if (string)
      fprintf(fp, " %8.4g", r=*point++);
    else
      fprintf(fp, qh_REAL_1, r=*point++);
  }
  fprintf(fp, "\n");
} /* printpoint */

/*------------------------------------
-printpoint3- prints 2-d , 3-d, or 4-d point as Geomview 3-d coordinates
*/
void qh_printpoint3 (FILE *fp, pointT *point) {
  int k;
  realT p[4];
  
  qh_projectdim3 (point, p);
  for (k=0; k< 3; k++)
    fprintf (fp, "%8.4g ", p[k]);
  fprintf (fp, " # p%d\n", qh_pointid (point));
} /* printpoint3 */

/*----------------------------------------
-printpoints- print pointids for a set of points starting at index 
   see geom.c
*/

/*----------------------------------------
-printpoints_out- prints vertices for facets by their point coordinates
  same format as qhull input
  allows CDDoutput
*/
void qh_printpoints_out (FILE *fp, facetT *facetlist, setT *facets, int printall) {
  int allpoints= qh nu