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);