realT radius;

  if (qh MERGING || qh JOGGLEmax < REALmax/2) {
    qh_outerinner (facet, outerplane, innerplane);
    radius= qh PRINTradius;
    if (qh JOGGLEmax < REALmax/2)
      radius -= qh JOGGLEmax * sqrt (qh hull_dim);  /* already accounted for in qh_outerinner() */
    *outerplane += radius;
    *innerplane -= radius;
    if (qh PRINTcoplanar || qh PRINTspheres) {
      *outerplane += qh MAXabs_coord * qh_GEOMepsilon;
      *innerplane -= qh MAXabs_coord * qh_GEOMepsilon;
    }
  }else 
    *innerplane= *outerplane= 0;
} /* geomplanes */


/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="markkeep">-</a>
  
  qh_markkeep( facetlist )
    mark good facets that meet qh.KEEParea, qh.KEEPmerge, and qh.KEEPminArea
    ignores visible facets (not part of convex hull)

  returns:
    may clear facet->good
    recomputes qh.num_good

  design:
    get set of good facets
    if qh.KEEParea
      sort facets by area
      clear facet->good for all but n largest facets
    if qh.KEEPmerge
      sort facets by merge count
      clear facet->good for all but n most merged facets
    if qh.KEEPminarea
      clear facet->good if area too small
    update qh.num_good    
*/
void qh_markkeep (facetT *facetlist) {
  facetT *facet, **facetp;
  setT *facets= qh_settemp (qh num_facets);
  int size, count;

  trace2((qh ferr, "qh_markkeep: only keep %d largest and/or %d most merged facets and/or min area %.2g\n",
          qh KEEParea, qh KEEPmerge, qh KEEPminArea));
  FORALLfacet_(facetlist) {
    if (!facet->visible && facet->good)
      qh_setappend (&facets, facet);
  }
  size= qh_setsize (facets);
  if (qh KEEParea) {
    qsort (SETaddr_(facets, facetT), size,
             sizeof (facetT *), qh_compare_facetarea);
    if ((count= size - qh KEEParea) > 0) {
      FOREACHfacet_(facets) {
        facet->good= False;
        if (--count == 0)
          break;
      }
    }
  }
  if (qh KEEPmerge) {
    qsort (SETaddr_(facets, facetT), size,
             sizeof (facetT *), qh_compare_facetmerge);
    if ((count= size - qh KEEPmerge) > 0) {
      FOREACHfacet_(facets) {
        facet->good= False;
        if (--count == 0)
          break;
      }
    }
  }
  if (qh KEEPminArea < REALmax/2) {
    FOREACHfacet_(facets) {
      if (!facet->isarea || facet->f.area < qh KEEPminArea)
	facet->good= False;
    }
  }
  qh_settempfree (&facets);
  count= 0;
  FORALLfacet_(facetlist) {
    if (facet->good)
      count++;
  }
  qh num_good= count;
} /* markkeep */


/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="markvoronoi">-</a>
  
  qh_markvoronoi( facetlist, facets, printall, islower, numcenters )
    mark voronoi vertices for printing by site pairs
  
  returns:
    temporary set of vertices indexed by pointid
    islower set if printing lower hull (i.e., at least one facet is lower hull)
    numcenters= total number of Voronoi vertices
    bumps qh.printoutnum for vertex-at-infinity
    clears all facet->seen and sets facet->seen2
    
    if selected
      facet->visitid= Voronoi vertex id
    else if upper hull (or 'Qu' and lower hull)
      facet->visitid= 0
    else
      facet->visitid >= qh num_facets
  
  notes:
    ignores qh ATinfinity, if defined
*/
setT *qh_markvoronoi (facetT *facetlist, setT *facets, boolT printall, boolT *islowerp, int *numcentersp) {
  int numcenters=0;
  facetT *facet, **facetp;
  setT *vertices;
  boolT islower= False;

  qh printoutnum++;
  qh_clearcenters (qh_ASvoronoi);  /* in case, qh_printvdiagram2 called by user */
  qh_vertexneighbors();
  vertices= qh_pointvertex();
  if (qh ATinfinity) 
    SETelem_(vertices, qh num_points-1)= NULL;
  qh visit_id++;
  maximize_(qh visit_id, (unsigned) qh num_facets);
  FORALLfacet_(facetlist) {  /* FIXUP: could merge with below */
    if (printall || !qh_skipfacet (facet)) {
      if (!facet->upperdelaunay) {
        islower= True;
	break;
      }
    }
  }
  FOREACHfacet_(facets) {
    if (printall || !qh_skipfacet (facet)) {
      if (!facet->upperdelaunay) {
        islower= True;
	break;
      }
    }
  }
  FORALLfacets {
    if (facet->normal && (facet->upperdelaunay == islower))
      facet->visitid= 0;  /* facetlist or facets may overwrite */
    else
      facet->visitid= qh visit_id;
    facet->seen= False;
    facet->seen2= True;
  }
  numcenters++;  /* qh_INFINITE */
  FORALLfacet_(facetlist) {
    if (printall || !qh_skipfacet (facet))
      facet->visitid= numcenters++;
  }
  FOREACHfacet_(facets) {
    if (printall || !qh_skipfacet (facet))
      facet->visitid= numcenters++;  
  }
  *islowerp= islower;
  *numcentersp= numcenters;
  trace2((qh ferr, "qh_markvoronoi: islower %d numcenters %d\n", islower, numcenters));
  return vertices;
} /* markvoronoi */

/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="order_vertexneighbors">-</a>
  
  qh_order_vertexneighbors( vertex )
    order facet neighbors of a 2-d or 3-d vertex by adjacency

  notes:
    does not orient the neighbors

  design:
    initialize a new neighbor set with the first facet in vertex->neighbors
    while vertex->neighbors non-empty
      select next neighbor in the previous facet's neighbor set
    set vertex->neighbors to the new neighbor set
*/
void qh_order_vertexneighbors(vertexT *vertex) {
  setT *newset;
  facetT *facet, *neighbor, **neighborp;

  trace4((qh ferr, "qh_order_vertexneighbors: order neighbors of v%d for 3-d\n", vertex->id));
  newset= qh_settemp (qh_setsize (vertex->neighbors));
  facet= (facetT*)qh_setdellast (vertex->neighbors);
  qh_setappend (&newset, facet);
  while (qh_setsize (vertex->neighbors)) {
    FOREACHneighbor_(vertex) {
      if (qh_setin (facet->neighbors, neighbor)) {
        qh_setdel(vertex->neighbors, neighbor);
        qh_setappend (&newset, neighbor);
        facet= neighbor;
        break;
      }
    }
    if (!neighbor) {
      fprintf (qh ferr, "qhull internal error (qh_order_vertexneighbors): no neighbor of v%d for f%d\n",
        vertex->id, facet->id);
      qh_errexit (qh_ERRqhull, facet, NULL);
    }
  }
  qh_setfree (&vertex->neighbors);
  qh_settemppop ();
  vertex->neighbors= newset;
} /* order_vertexneighbors */

/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="printafacet">-</a>
  
  qh_printafacet( fp, format, facet, printall )
    print facet to fp in given output format (see qh.PRINTout)

  returns:
    nop if !printall and qh_skipfacet()
    nop if visible facet and NEWfacets and format != PRINTfacets
    must match qh_countfacets

  notes
    preserves qh.visit_id
    facet->normal may be null if PREmerge/MERGEexact and STOPcone before merge

  see
    qh_printbegin() and qh_printend()

  design:
    test for printing facet
    call appropriate routine for format
    or output results directly
*/
void qh_printafacet(FILE *fp, int format, facetT *facet, boolT printall) {
  realT color[4], offset, dist, outerplane, innerplane;
  boolT zerodiv;
  coordT *point, *normp, *coordp, **pointp, *feasiblep;
  int k;
  vertexT *vertex, **vertexp;
  facetT *neighbor, **neighborp;

  if (!printall && qh_skipfacet (facet))
    return;
  if (facet->visible && qh NEWfacets && format != qh_PRINTfacets)
    return;
  qh printoutnum++;
  switch (format) {
  case qh_PRINTarea:
    if (facet->isarea) {
      fprintf (fp, qh_REAL_1, facet->f.area);
      fprintf (fp, "\n");
    }else
      fprintf (fp, "0\n");
    break;
  case qh_PRINTcoplanars:
    fprintf (fp, "%d", qh_setsize (facet->coplanarset));
    FOREACHpoint_(facet->coplanarset)
      fprintf (fp, " %d", qh_pointid (point));
    fprintf (fp, "\n");
    break;
  case qh_PRINTcentrums:
    qh_printcenter (fp, format, NULL, facet);
    break;
  case qh_PRINTfacets:
    qh_printfacet (fp, facet);
    break;
  case qh_PRINTfacets_xridge:
    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_PRINTincidences:
  case qh_PRINToff:
  case qh_PRINTtriangles:
    if (qh hull_dim == 3 && format != qh_PRINTtriangles) 
      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++, format);
    break;
  case qh_PRINTinner:
    qh_outerinner (facet, NULL, &innerplane);
    offset= facet->offset - innerplane;
    goto LABELprintnorm;
    break; /* prevent warning */
  case qh_PRINTmerges:
    fprintf (fp, "%d\n", facet->nummerge);
    break;
  case qh_PRINTnormals:
    offset= facet->offset;
    goto LABELprintnorm;
    break; /* prevent warning */
  case qh_PRINTouter:
    qh_outerinner (facet, &outerplane, NULL);
    offset= facet->offset - outerplane;
  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));