text: /*<html><pre>  -<a                             href="qh-io_r.htm"
    >-------------------------------</a><a name="TOP">-</a>
  
     io_r.c
     Input/Output routines of qhull application
  
     see qh-io_r.htm and io_r.h
  
     see user_r.c for qh_errprint and qh_printfacetlist
  
     unix_r.c calls qh_readpoints and qh_produce_output
  
     unix_r.c and user_r.c are the only callers of io_r.c functions
     This allows the user to avoid loading io_r.o from qhull.a
  
     Copyright (c) 1993-2020 The Geometry Center.
     $Id: //main/2019/qhull/src/libqhull_r/io_r.c#12 $$Change: 2965 $
     $DateTime: 2020/06/04 15:37:41 $$Author: bbarber $
  */
  
  #include "qhull_ra.h"
  
  /*========= -functions in alphabetical order after qh_produce_output(qh)  =====*/
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="produce_output">-</a>
  
    qh_produce_output(qh )
    qh_produce_output2(qh )
      prints out the result of qhull in desired format
      qh_produce_output2 does not call qh_prepare_output
        qh_checkpolygon is valid for qh_prepare_output
      if qh.GETarea
        computes and prints area and volume
      qh.PRINTout[] is an array of output formats
  
    notes:
      prints output in qh.PRINTout order
  */
  void qh_produce_output(qhT *qh) {
      int tempsize= qh_setsize(qh, qh->qhmem.tempstack);
  
      qh_prepare_output(qh);
      qh_produce_output2(qh);
      if (qh_setsize(qh, qh->qhmem.tempstack) != tempsize) {
          qh_fprintf(qh, qh->ferr, 6206, "qhull internal error (qh_produce_output): temporary sets not empty(%d)\n",
              qh_setsize(qh, qh->qhmem.tempstack));
          qh_errexit(qh, qh_ERRqhull, NULL, NULL);
      }
  } /* produce_output */
  
  
  void qh_produce_output2(qhT *qh) {
    int i, tempsize= qh_setsize(qh, qh->qhmem.tempstack), d_1;
  
    fflush(NULL);
    if (qh->PRINTsummary)
      qh_printsummary(qh, qh->ferr);
    else if (qh->PRINTout[0] == qh_PRINTnone)
      qh_printsummary(qh, qh->fout);
    for (i=0; i < qh_PRINTEND; i++)
      qh_printfacets(qh, qh->fout, qh->PRINTout[i], qh->facet_list, NULL, !qh_ALL);
    fflush(NULL);
  
    qh_allstatistics(qh);
    if (qh->PRINTprecision && !qh->MERGING && (qh->JOGGLEmax > REALmax/2 || qh->RERUN))
      qh_printstats(qh, qh->ferr, qh->qhstat.precision, NULL);
    if (qh->VERIFYoutput && (zzval_(Zridge) > 0 || zzval_(Zridgemid) > 0))
      qh_printstats(qh, qh->ferr, qh->qhstat.vridges, NULL);
    if (qh->PRINTstatistics) {
      qh_printstatistics(qh, qh->ferr, "");
      qh_memstatistics(qh, qh->ferr);
      d_1= (int)sizeof(setT) + (qh->hull_dim - 1) * SETelemsize;
      qh_fprintf(qh, qh->ferr, 8040, "\
      size in bytes: merge %d ridge %d vertex %d facet %d\n\
           normal %d ridge vertices %d facet vertices or neighbors %d\n",
              (int)sizeof(mergeT), (int)sizeof(ridgeT),
              (int)sizeof(vertexT), (int)sizeof(facetT),
              qh->normal_size, d_1, d_1 + SETelemsize);
    }
    if (qh_setsize(qh, qh->qhmem.tempstack) != tempsize) {
      qh_fprintf(qh, qh->ferr, 6065, "qhull internal error (qh_produce_output2): temporary sets not empty(%d)\n",
               qh_setsize(qh, qh->qhmem.tempstack));
      qh_errexit(qh, qh_ERRqhull, NULL, NULL);
    }
  } /* produce_output2 */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="dfacet">-</a>
  
    qh_dfacet(qh, id )
      print facet by id, for debugging
  
  */
  void qh_dfacet(qhT *qh, unsigned int id) {
    facetT *facet;
  
    FORALLfacets {
      if (facet->id == id) {
        qh_printfacet(qh, qh->fout, facet);
        break;
      }
    }
  } /* dfacet */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="dvertex">-</a>
  
    qh_dvertex(qh, id )
      print vertex by id, for debugging
  */
  void qh_dvertex(qhT *qh, unsigned int id) {
    vertexT *vertex;
  
    FORALLvertices {
      if (vertex->id == id) {
        qh_printvertex(qh, qh->fout, vertex);
        break;
      }
    }
  } /* dvertex */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="compare_facetarea">-</a>
  
    qh_compare_facetarea( p1, p2 )
      used by qsort() to order facets by area
  */
  int qh_compare_facetarea(const void *p1, const void *p2) {
    const facetT *a= *((facetT *const*)p1), *b= *((facetT *const*)p2);
  
    if (!a->isarea)
      return -1;
    if (!b->isarea)
      return 1;
    if (a->f.area > b->f.area)
      return 1;
    else if (a->f.area == b->f.area)
      return 0;
    return -1;
  } /* compare_facetarea */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="compare_facetvisit">-</a>
  
    qh_compare_facetvisit( p1, p2 )
      used by qsort() to order facets by visit id or id
  */
  int qh_compare_facetvisit(const void *p1, const void *p2) {
    const facetT *a= *((facetT *const*)p1), *b= *((facetT *const*)p2);
    int i,j;
  
    if (!(i= (int)a->visitid))
      i= (int)(0 - a->id); /* sign distinguishes id from visitid */
    if (!(j= (int)b->visitid))
      j= (int)(0 - b->id);
    return(i - j);
  } /* compare_facetvisit */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="compare_nummerge">-</a>
  
    qh_compare_nummerge( p1, p2 )
      used by qsort() to order facets by number of merges
  
  notes:
      called by qh_markkeep ('PMerge-keep')
  */
  int qh_compare_nummerge(const void *p1, const void *p2) {
    const facetT *a= *((facetT *const*)p1), *b= *((facetT *const*)p2);
  
    return(a->nummerge - b->nummerge);
  } /* compare_nummerge */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="copyfilename">-</a>
  
    qh_copyfilename(qh, dest, size, source, length )
      copy filename identified by qh_skipfilename()
  
    notes:
      see qh_skipfilename() for syntax
  */
  void qh_copyfilename(qhT *qh, char *filename, int size, const char* source, int length) {
    char c= *source;
  
    if (length > size + 1) {
        qh_fprintf(qh, qh->ferr, 6040, "qhull error: filename is more than %d characters, %s\n",  size-1, source);
        qh_errexit(qh, qh_ERRinput, NULL, NULL);
    }
    strncpy(filename, source, (size_t)length);
    filename[length]= '\0';
    if (c == '\'' || c == '"') {
      char *s= filename + 1;
      char *t= filename;
      while (*s) {
        if (*s == c) {
            if (s[-1] == '\\')
                t[-1]= c;
        }else
            *t++= *s;
        s++;
      }
      *t= '\0';
    }
  } /* copyfilename */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="countfacets">-</a>
  
    qh_countfacets(qh, facetlist, facets, printall,
            numfacets, numsimplicial, totneighbors, numridges, numcoplanar, numtricoplanars  )
      count good facets for printing and set visitid
      if allfacets, ignores qh_skipfacet()
  
    notes:
      qh_printsummary and qh_countfacets must match counts
  
    returns:
      numfacets, numsimplicial, total neighbors, numridges, coplanars
      each facet with ->visitid indicating 1-relative position
        ->visitid==0 indicates not good
  
    notes
      numfacets >= numsimplicial
      if qh.NEWfacets,
        does not count visible facets (matches qh_printafacet)
  
    design:
      for all facets on facetlist and in facets set
        unless facet is skipped or visible (i.e., will be deleted)
          mark facet->visitid
          update counts
  */
  void qh_countfacets(qhT *qh, facetT *facetlist, setT *facets, boolT printall,
      int *numfacetsp, int *numsimplicialp, int *totneighborsp, int *numridgesp, int *numcoplanarsp, int *numtricoplanarsp) {
    facetT *facet, **facetp;
    int numfacets= 0, numsimplicial= 0, numridges= 0, totneighbors= 0, numcoplanars= 0, numtricoplanars= 0;
  
    FORALLfacet_(facetlist) {
      if ((facet->visible && qh->NEWfacets)
      || (!printall && qh_skipfacet(qh, facet)))
        facet->visitid= 0;
      else {
        facet->visitid= (unsigned int)(++numfacets);
        totneighbors += qh_setsize(qh, facet->neighbors);
        if (facet->simplicial) {
          numsimplicial++;
          if (facet->keepcentrum && facet->tricoplanar)
            numtricoplanars++;
        }else
          numridges += qh_setsize(qh, facet->ridges);
        if (facet->coplanarset)
          numcoplanars += qh_setsize(qh, facet->coplanarset);
      }
    }
  
    FOREACHfacet_(facets) {
      if ((facet->visible && qh->NEWfacets)
      || (!printall && qh_skipfacet(qh, facet)))
        facet->visitid= 0;
      else {
        facet->visitid= (unsigned int)(++numfacets);
        totneighbors += qh_setsize(qh, facet->neighbors);
        if (facet->simplicial){
          numsimplicial++;
          if (facet->keepcentrum && facet->tricoplanar)
            numtricoplanars++;
        }else
          numridges += qh_setsize(qh, facet->ridges);
        if (facet->coplanarset)
          numcoplanars += qh_setsize(qh, facet->coplanarset);
      }
    }
    qh->visit_id += (unsigned int)numfacets + 1;
    *numfacetsp= numfacets;
    *numsimplicialp= numsimplicial;
    *totneighborsp= totneighbors;
    *numridgesp= numridges;
    *numcoplanarsp= numcoplanars;
    *numtricoplanarsp= numtricoplanars;
  } /* countfacets */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="detvnorm">-</a>
  
    qh_detvnorm(qh, vertex, vertexA, centers, offset )
      compute separating plane of the Voronoi diagram for a pair of input sites
      centers= set of facets (i.e., Voronoi vertices)
        facet->visitid= 0 iff vertex-at-infinity (i.e., unbounded)
  
    assumes:
      qh_ASvoronoi and qh_vertexneighbors() already set
  
    returns:
      norm
        a pointer into qh.gm_matrix to qh.hull_dim-1 reals
        copy the data before reusing qh.gm_matrix
      offset
        if 'QVn'
          sign adjusted so that qh.GOODvertexp is inside
        else
          sign adjusted so that vertex is inside
  
      qh.gm_matrix= simplex of points from centers relative to first center
  
    notes:
      in io_r.c so that code for 'v Tv' can be removed by removing io_r.c
      returns pointer into qh.gm_matrix to avoid tracking of temporary memory
  
    design:
      determine midpoint of input sites
      build points as the set of Voronoi vertices
      select a simplex from points (if necessary)
        include midpoint if the Voronoi region is unbounded
      relocate the first vertex of the simplex to the origin
      compute the normalized hyperplane through the simplex
      orient the hyperplane toward 'QVn' or 'vertex'
      if 'Tv' or 'Ts'
        if bounded
          test that hyperplane is the perpendicular bisector of the input sites
        test that Voronoi vertices not in the simplex are still on the hyperplane
      free up temporary memory
  */
  pointT *qh_detvnorm(qhT *qh, vertexT *vertex, vertexT *vertexA, setT *centers, realT *offsetp) {
    facetT *facet, **facetp;
    int  i, k, pointid, pointidA, point_i, point_n;
    setT *simplex= NULL;
    pointT *point, **pointp, *point0, *midpoint, *normal, *inpoint;
    coordT *coord, *gmcoord, *normalp;
    setT *points= qh_settemp(qh, qh->TEMPsize);
    boolT nearzero= False;
    boolT unbounded= False;
    int numcenters= 0;
    int dim= qh->hull_dim - 1;
    realT dist, offset, angle, zero= 0.0;
  
    midpoint= qh->gm_matrix + qh->hull_dim * qh->hull_dim;  /* last row */
    for (k=0; k < dim; k++)
      midpoint[k]= (vertex->point[k] + vertexA->point[k])/2;
    FOREACHfacet_(centers) {
      numcenters++;
      if (!facet->visitid)
        unbounded= True;
      else {
        if (!facet->center)
          facet->center= qh_facetcenter(qh, facet->vertices);
        qh_setappend(qh, &points, facet->center);
      }
    }
    if (numcenters > dim) {
      simplex= qh_settemp(qh, qh->TEMPsize);
      qh_setappend(qh, &simplex, vertex->point);
      if (unbounded)
        qh_setappend(qh, &simplex, midpoint);
      qh_maxsimplex(qh, dim, points, NULL, 0, &simplex);
      qh_setdelnth(qh, simplex, 0);
    }else if (numcenters == dim) {
      if (unbounded)
        qh_setappend(qh, &points, midpoint);
      simplex= points;
    }else {
      qh_fprintf(qh, qh->ferr, 6216, "qhull internal error (qh_detvnorm): too few points(%d) to compute separating plane\n", numcenters);
      qh_errexit(qh, qh_ERRqhull, NULL, NULL);
    }
    i= 0;
    gmcoord= qh->gm_matrix;
    point0= SETfirstt_(simplex, pointT);
    FOREACHpoint_(simplex) {
      if (qh->IStracing >= 4)
        qh_printmatrix(qh, qh->ferr, "qh_detvnorm: Voronoi vertex or midpoint",
                                &point, 1, dim);
      if (point != point0) {
        qh->gm_row[i++]= gmcoord;
        coord= point0;
        for (k=dim; k--; )
          *(gmcoord++)= *point++ - *coord++;
      }
    }
    qh->gm_row[i]= gmcoord;  /* does not overlap midpoint, may be used later for qh_areasimplex */
    normal= gmcoord;
    qh_sethyperplane_gauss(qh, dim, qh->gm_row, point0, True,
                  normal, &offset, &nearzero);
    /* nearzero is true for axis-parallel hyperplanes (e.g., a bounding box).  Should detect degenerate hyperplanes.  See 'Tv' check following */
    if (qh->GOODvertexp == vertexA->point)
      inpoint= vertexA->point;
    else
      inpoint= vertex->point;
    zinc_(Zdistio);
    dist= qh_distnorm(dim, inpoint, normal, &offset);
    if (dist > 0) {
      offset= -offset;
      normalp= normal;
      for (k=dim; k--; ) {
        *normalp= -(*normalp);
        normalp++;
      }
    }
    if (qh->VERIFYoutput || qh->PRINTstatistics) {
      pointid= qh_pointid(qh, vertex->point);
      pointidA= qh_pointid(qh, vertexA->point);
      if (!unbounded) {
        zinc_(Zdiststat);
        dist= qh_distnorm(dim, midpoint, normal, &offset);
        if (dist < 0)
          dist= -dist;
        zzinc_(Zridgemid);
        wwmax_(Wridgemidmax, dist);
        wwadd_(Wridgemid, dist);
        trace4((qh, qh->ferr, 4014, "qh_detvnorm: points %d %d midpoint dist %2.2g\n",
                   pointid, pointidA, dist));
        for (k=0; k < dim; k++)
          midpoint[k]= vertexA->point[k] - vertex->point[k];  /* overwrites midpoint! */
        qh_normalize(qh, midpoint, dim, False);
        angle= qh_distnorm(dim, midpoint, normal, &zero); /* qh_detangle uses dim+1 */
        if (angle < 0.0)
          angle= angle + 1.0;
        else
          angle= angle - 1.0;
        if (angle < 0.0)
          angle= -angle;
        trace4((qh, qh->ferr, 4015, "qh_detvnorm: points %d %d angle %2.2g nearzero %d\n",
                   pointid, pointidA, angle, nearzero));
        if (nearzero) {
          zzinc_(Zridge0);
          wwmax_(Wridge0max, angle);
          wwadd_(Wridge0, angle);
        }else {
          zzinc_(Zridgeok)
          wwmax_(Wridgeokmax, angle);
          wwadd_(Wridgeok, angle);
        }
      }
      if (simplex != points) {
        FOREACHpoint_i_(qh, points) {
          if (!qh_setin(simplex, point)) {
            facet= SETelemt_(centers, point_i, facetT);
            zinc_(Zdiststat);
            dist= qh_distnorm(dim, point, normal, &offset);
            if (dist < 0)
              dist= -dist;
            zzinc_(Zridge);
            wwmax_(Wridgemax, dist);
            wwadd_(Wridge, dist);
            trace4((qh, qh->ferr, 4016, "qh_detvnorm: points %d %d Voronoi vertex %d dist %2.2g\n",
                               pointid, pointidA, facet->visitid, dist));
          }
        }
      }
    }
    *offsetp= offset;
    if (simplex != points)
      qh_settempfree(qh, &simplex);
    qh_settempfree(qh, &points);
    return normal;
  } /* detvnorm */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="detvridge">-</a>
  
    qh_detvridge(qh, vertexA )
      determine Voronoi ridge from 'seen' neighbors of vertexA
      include one vertex-at-infinite if an !neighbor->visitid
  
    returns:
      temporary set of centers (facets, i.e., Voronoi vertices)
      sorted by center id
  */
  setT *qh_detvridge(qhT *qh, vertexT *vertex) {
    setT *centers= qh_settemp(qh, qh->TEMPsize);
    setT *tricenters= qh_settemp(qh, qh->TEMPsize);
    facetT *neighbor, **neighborp;
    boolT firstinf= True;
  
    FOREACHneighbor_(vertex) {
      if (neighbor->seen) {
        if (neighbor->visitid) {
          if (!neighbor->tricoplanar || qh_setunique(qh, &tricenters, neighbor->center))
            qh_setappend(qh, &centers, neighbor);
        }else if (firstinf) {
          firstinf= False;
          qh_setappend(qh, &centers, neighbor);
        }
      }
    }
    qsort(SETaddr_(centers, facetT), (size_t)qh_setsize(qh, centers),
               sizeof(facetT *), qh_compare_facetvisit);
    qh_settempfree(qh, &tricenters);
    return centers;
  } /* detvridge */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="detvridge3">-</a>
  
    qh_detvridge3(qh, atvertex, vertex )
      determine 3-d Voronoi ridge from 'seen' neighbors of atvertex and vertex
      include one vertex-at-infinite for !neighbor->visitid
      assumes all facet->seen2= True
  
    returns:
      temporary set of centers (facets, i.e., Voronoi vertices)
      listed in adjacency order (!oriented)
      all facet->seen2= True
  
    design:
      mark all neighbors of atvertex
      for each adjacent neighbor of both atvertex and vertex
        if neighbor selected
          add neighbor to set of Voronoi vertices
  */
  setT *qh_detvridge3(qhT *qh, vertexT *atvertex, vertexT *vertex) {
    setT *centers= qh_settemp(qh, qh->TEMPsize);
    setT *tricenters= qh_settemp(qh, qh->TEMPsize);
    facetT *neighbor, **neighborp, *facet= NULL;
    boolT firstinf= True;
  
    FOREACHneighbor_(atvertex)
      neighbor->seen2= False;
    FOREACHneighbor_(vertex) {
      if (!neighbor->seen2) {
        facet= neighbor;
        break;
      }
    }
    while (facet) {
      facet->seen2= True;
      if (neighbor->seen) {
        if (facet->visitid) {
          if (!facet->tricoplanar || qh_setunique(qh, &tricenters, facet->center))
            qh_setappend(qh, &centers, facet);
        }else if (firstinf) {
          firstinf= False;
          qh_setappend(qh, &centers, facet);
        }
      }
      FOREACHneighbor_(facet) {
        if (!neighbor->seen2) {
          if (qh_setin(vertex->neighbors, neighbor))
            break;
          else
            neighbor->seen2= True;
        }
      }
      facet= neighbor;
    }
    if (qh->CHECKfrequently) {
      FOREACHneighbor_(vertex) {
        if (!neighbor->seen2) {
            qh_fprintf(qh, qh->ferr, 6217, "qhull internal error (qh_detvridge3): neighbors of vertex p%d are not connected at facet %d\n",
                   qh_pointid(qh, vertex->point), neighbor->id);
          qh_errexit(qh, qh_ERRqhull, neighbor, NULL);
        }
      }
    }
    FOREACHneighbor_(atvertex)
      neighbor->seen2= True;
    qh_settempfree(qh, &tricenters);
    return centers;
  } /* detvridge3 */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="eachvoronoi">-</a>
  
    qh_eachvoronoi(qh, fp, printvridge, vertex, visitall, innerouter, inorder )
      if visitall,
        visit all Voronoi ridges for vertex (i.e., an input site)
      else
        visit all unvisited Voronoi ridges for vertex
        all vertex->seen= False if unvisited
      assumes
        all facet->seen= False
        all facet->seen2= True (for qh_detvridge3)
        all facet->visitid == 0 if vertex_at_infinity
                           == index of Voronoi vertex
                           >= qh.num_facets if ignored
      innerouter:
        qh_RIDGEall--  both inner (bounded) and outer(unbounded) ridges
        qh_RIDGEinner- only inner
        qh_RIDGEouter- only outer
  
      if inorder
        orders vertices for 3-d Voronoi diagrams
  
    returns:
      number of visited ridges (does not include previously visited ridges)
  
      if printvridge,
        calls printvridge( fp, vertex, vertexA, centers)
          fp== any pointer (assumes FILE*)
               fp may be NULL for QhullQh::qh_fprintf which calls appendQhullMessage
          vertex,vertexA= pair of input sites that define a Voronoi ridge
          centers= set of facets (i.e., Voronoi vertices)
                   ->visitid == index or 0 if vertex_at_infinity
                   ordered for 3-d Voronoi diagram
    notes:
      uses qh.vertex_visit
  
    see:
      qh_eachvoronoi_all()
  
    design:
      mark selected neighbors of atvertex
      for each selected neighbor (either Voronoi vertex or vertex-at-infinity)
        for each unvisited vertex
          if atvertex and vertex share more than d-1 neighbors
            bump totalcount
            if printvridge defined
              build the set of shared neighbors (i.e., Voronoi vertices)
              call printvridge
  */
  int qh_eachvoronoi(qhT *qh, FILE *fp, printvridgeT printvridge, vertexT *atvertex, boolT visitall, qh_RIDGE innerouter, boolT inorder) {
    boolT unbounded;
    int count;
    facetT *neighbor, **neighborp, *neighborA, **neighborAp;
    setT *centers;
    setT *tricenters= qh_settemp(qh, qh->TEMPsize);
  
    vertexT *vertex, **vertexp;
    boolT firstinf;
    unsigned int numfacets= (unsigned int)qh->num_facets;
    int totridges= 0;
  
    qh->vertex_visit++;
    atvertex->seen= True;
    if (visitall) {
      FORALLvertices
        vertex->seen= False;
    }
    FOREACHneighbor_(atvertex) {
      if (neighbor->visitid < numfacets)
        neighbor->seen= True;
    }
    FOREACHneighbor_(atvertex) {
      if (neighbor->seen) {
        FOREACHvertex_(neighbor->vertices) {
          if (vertex->visitid != qh->vertex_visit && !vertex->seen) {
            vertex->visitid= qh->vertex_visit;
            count= 0;
            firstinf= True;
            qh_settruncate(qh, tricenters, 0);
            FOREACHneighborA_(vertex) {
              if (neighborA->seen) {
                if (neighborA->visitid) {
                  if (!neighborA->tricoplanar || qh_setunique(qh, &tricenters, neighborA->center))
                    count++;
                }else if (firstinf) {
                  count++;
                  firstinf= False;
                }
              }
            }
            if (count >= qh->hull_dim - 1) {  /* e.g., 3 for 3-d Voronoi */
              if (firstinf) {
                if (innerouter == qh_RIDGEouter)
                  continue;
                unbounded= False;
              }else {
                if (innerouter == qh_RIDGEinner)
                  continue;
                unbounded= True;
              }
              totridges++;
              trace4((qh, qh->ferr, 4017, "qh_eachvoronoi: Voronoi ridge of %d vertices between sites %d and %d\n",
                    count, qh_pointid(qh, atvertex->point), qh_pointid(qh, vertex->point)));
              if (printvridge) {
                if (inorder && qh->hull_dim == 3+1) /* 3-d Voronoi diagram */
                  centers= qh_detvridge3(qh, atvertex, vertex);
                else
                  centers= qh_detvridge(qh, vertex);
                (*printvridge)(qh, fp, atvertex, vertex, centers, unbounded);
                qh_settempfree(qh, &centers);
              }
            }
          }
        }
      }
    }
    FOREACHneighbor_(atvertex)
      neighbor->seen= False;
    qh_settempfree(qh, &tricenters);
    return totridges;
  } /* eachvoronoi */
  
  
  /*-<a                             href="qh-poly_r.htm#TOC"
    >-------------------------------</a><a name="eachvoronoi_all">-</a>
  
    qh_eachvoronoi_all(qh, fp, printvridge, isUpper, innerouter, inorder )
      visit all Voronoi ridges
  
      innerouter:
        see qh_eachvoronoi()
  
      if inorder
        orders vertices for 3-d Voronoi diagrams
  
    returns
      total number of ridges
  
      if isUpper == facet->upperdelaunay  (i.e., a Vornoi vertex)
        facet->visitid= Voronoi vertex index(same as 'o' format)
      else
        facet->visitid= 0
  
      if printvridge,
        calls printvridge( fp, vertex, vertexA, centers)
        [see qh_eachvoronoi]
  
    notes:
      Not used for qhull.exe
      same effect as qh_printvdiagram but ridges not sorted by point id
  */
  int qh_eachvoronoi_all(qhT *qh, FILE *fp, printvridgeT printvridge, boolT isUpper, qh_RIDGE innerouter, boolT inorder) {
    facetT *facet;
    vertexT *vertex;
    int numcenters= 1;  /* vertex 0 is vertex-at-infinity */
    int totridges= 0;
  
    qh_clearcenters(qh, qh_ASvoronoi);
    qh_vertexneighbors(qh);
    maximize_(qh->visit_id, (unsigned int)qh->num_facets);
    FORALLfacets {
      facet->visitid= 0;
      facet->seen= False;
      facet->seen2= True;
    }
    FORALLfacets {
      if (facet->upperdelaunay == isUpper)
        facet->visitid= (unsigned int)(numcenters++);
    }
    FORALLvertices
      vertex->seen= False;
    FORALLvertices {
      if (qh->GOODvertex > 0 && qh_pointid(qh, vertex->point)+1 != qh->GOODvertex)
        continue;
      totridges += qh_eachvoronoi(qh, fp, printvridge, vertex,
                     !qh_ALL, innerouter, inorder);
    }
    return totridges;
  } /* eachvoronoi_all */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="facet2point">-</a>
  
    qh_facet2point(qh, facet, point0, point1, mindist )
      return two projected temporary vertices for a 2-d facet
      may be non-simplicial
  
    returns:
      point0 and point1 oriented and projected to the facet
      returns mindist (maximum distance below plane)
  */
  void qh_facet2point(qhT *qh, facetT *facet, pointT **point0, pointT **point1, realT *mindist) {
    vertexT *vertex0, *vertex1;
    realT dist;
  
    if (facet->toporient ^ qh_ORIENTclock) {
      vertex0= SETfirstt_(facet->vertices, vertexT);
      vertex1= SETsecondt_(facet->vertices, vertexT);
    }else {
      vertex1= SETfirstt_(facet->vertices, vertexT);
      vertex0= SETsecondt_(facet->vertices, vertexT);
    }
    zadd_(Zdistio, 2);
    qh_distplane(qh, vertex0->point, facet, &dist);
    *mindist= dist;
    *point0= qh_projectpoint(qh, vertex0->point, facet, dist);
    qh_distplane(qh, vertex1->point, facet, &dist);
    minimize_(*mindist, dist);
    *point1= qh_projectpoint(qh, vertex1->point, facet, dist);
  } /* facet2point */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="facetvertices">-</a>
  
    qh_facetvertices(qh, facetlist, facets, allfacets )
      returns temporary set of vertices in a set and/or list of facets
      if allfacets, ignores qh_skipfacet()
  
    returns:
      vertices with qh.vertex_visit
  
    notes:
      optimized for allfacets of facet_list
  
    design:
      if allfacets of facet_list
        create vertex set from vertex_list
      else
        for each selected facet in facets or facetlist
          append unvisited vertices to vertex set
  */
  setT *qh_facetvertices(qhT *qh, facetT *facetlist, setT *facets, boolT allfacets) {
    setT *vertices;
    facetT *facet, **facetp;
    vertexT *vertex, **vertexp;
  
    qh->vertex_visit++;
    if (facetlist == qh->facet_list && allfacets && !facets) {
      vertices= qh_settemp(qh, qh->num_vertices);
      FORALLvertices {
        vertex->visitid= qh->vertex_visit;
        qh_setappend(qh, &vertices, vertex);
      }
    }else {
      vertices= qh_settemp(qh, qh->TEMPsize);
      FORALLfacet_(facetlist) {
        if (!allfacets && qh_skipfacet(qh, facet))
          continue;
        FOREACHvertex_(facet->vertices) {
          if (vertex->visitid != qh->vertex_visit) {
            vertex->visitid= qh->vertex_visit;
            qh_setappend(qh, &vertices, vertex);
          }
        }
      }
    }
    FOREACHfacet_(facets) {
      if (!allfacets && qh_skipfacet(qh, facet))
        continue;
      FOREACHvertex_(facet->vertices) {
        if (vertex->visitid != qh->vertex_visit) {
          vertex->visitid= qh->vertex_visit;
          qh_setappend(qh, &vertices, vertex);
        }
      }
    }
    return vertices;
  } /* facetvertices */
  
  /*-<a                             href="qh-geom_r.htm#TOC"
    >-------------------------------</a><a name="geomplanes">-</a>
  
    qh_geomplanes(qh, facet, outerplane, innerplane )
      return outer and inner planes for Geomview
      qh.PRINTradius is size of vertices and points (includes qh.JOGGLEmax)
  
    notes:
      assume precise calculations in io_r.c with roundoff covered by qh_GEOMepsilon
  */
  void qh_geomplanes(qhT *qh, facetT *facet, realT *outerplane, realT *innerplane) {
    realT radius;
  
    if (qh->MERGING || qh->JOGGLEmax < REALmax/2) {
      qh_outerinner(qh, facet, outerplane, innerplane);
      radius= qh->PRINTradius;
      if (qh->JOGGLEmax < REALmax/2)
        radius -= qh->JOGGLEmax * sqrt((realT)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-io_r.htm#TOC"
    >-------------------------------</a><a name="markkeep">-</a>
  
    qh_markkeep(qh, facetlist )
      restrict good facets for qh.KEEParea, qh.KEEPmerge, and qh.KEEPminArea
      ignores visible facets (!part of convex hull)
  
    returns:
      may clear facet->good
      recomputes qh.num_good
  
    notes:
      only called by qh_prepare_output after qh_findgood_all
      does not throw errors except memory/corruption of qset_r.c
  
    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(qhT *qh, facetT *facetlist) {
    facetT *facet, **facetp;
    setT *facets= qh_settemp(qh, qh->num_facets);
    int size, count;
  
    trace2((qh, qh->ferr, 2006, "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(qh, &facets, facet);
    }
    size= qh_setsize(qh, facets);
    if (qh->KEEParea) {
      qsort(SETaddr_(facets, facetT), (size_t)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_t)size,
               sizeof(facetT *), qh_compare_nummerge);
      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(qh, &facets);
    count= 0;
    FORALLfacet_(facetlist) {
      if (facet->good)
        count++;
    }
    qh->num_good= count;
  } /* markkeep */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="markvoronoi">-</a>
  
    qh_markvoronoi(qh, 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(qhT *qh, 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, qh_ASvoronoi);  /* in case, qh_printvdiagram2 called by user */
    qh_vertexneighbors(qh);
    vertices= qh_pointvertex(qh);
    if (qh->ATinfinity)
      SETelem_(vertices, qh->num_points-1)= NULL;
    qh->visit_id++;
    maximize_(qh->visit_id, (unsigned int)qh->num_facets);
    FORALLfacet_(facetlist) {
      if (printall || !qh_skipfacet(qh, facet)) {
        if (!facet->upperdelaunay) {
          isLower= True;
          break;
        }
      }
    }
    FOREACHfacet_(facets) {
      if (printall || !qh_skipfacet(qh, 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(qh, facet))
        facet->visitid= (unsigned int)(numcenters++);
    }
    FOREACHfacet_(facets) {
      if (printall || !qh_skipfacet(qh, facet))
        facet->visitid= (unsigned int)(numcenters++);
    }
    *isLowerp= isLower;
    *numcentersp= numcenters;
    trace2((qh, qh->ferr, 2007, "qh_markvoronoi: isLower %d numcenters %d\n", isLower, numcenters));
    return vertices;
  } /* markvoronoi */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="order_vertexneighbors">-</a>
  
    qh_order_vertexneighbors(qh, vertex )
      order facet neighbors of vertex by 2-d (orientation), 3-d (adjacency), or n-d (f.visitid,id)
  
    notes:
      error if qh_vertexneighbors not called beforehand
      only 2-d orients the neighbors
      for 4-d and higher
        set or clear f.visitid for qh_compare_facetvisit
        for example, use qh_markvoronoi (e.g., qh_printvornoi) or qh_countfacets (e.g., qh_printvneighbors)
  
    design (2-d):
      see qh_printextremes_2d
    design (3-d):
      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
    design (n-d):
      qsort by f.visitid, or f.facetid (qh_compare_facetvisit)
      facet_id is negated (sorted before visit_id facets)
  */
  void qh_order_vertexneighbors(qhT *qh, vertexT *vertex) {
    setT *newset;
    facetT *facet, *facetA, *facetB, *neighbor, **neighborp;
    vertexT *vertexA;
    int numneighbors;
  
    trace4((qh, qh->ferr, 4018, "qh_order_vertexneighbors: order facet neighbors of v%d by 2-d (orientation), 3-d (adjacency), or n-d (f.visitid,id)\n", vertex->id));
    if (!qh->VERTEXneighbors) {
      qh_fprintf(qh, qh->ferr, 6428, "qhull internal error (qh_order_vertexneighbors): call qh_vertexneighbors before calling qh_order_vertexneighbors\n");
      qh_errexit(qh, qh_ERRqhull, NULL, NULL);
    }
    if (qh->hull_dim == 2) {
      facetA= SETfirstt_(vertex->neighbors, facetT);
      if (facetA->toporient ^ qh_ORIENTclock)
        vertexA= SETfirstt_(facetA->vertices, vertexT);
      else
        vertexA= SETsecondt_(facetA->vertices, vertexT);
      if (vertexA!=vertex) {
        facetB= SETsecondt_(vertex->neighbors, facetT);
        SETfirst_(vertex->neighbors)= facetB;
        SETsecond_(vertex->neighbors)= facetA;
      }
    }else if (qh->hull_dim == 3) {
      newset= qh_settemp(qh, qh_setsize(qh, vertex->neighbors));
      facet= (facetT *)qh_setdellast(vertex->neighbors);
      qh_setappend(qh, &newset, facet);
      while (qh_setsize(qh, vertex->neighbors)) {
        FOREACHneighbor_(vertex) {
          if (qh_setin(facet->neighbors, neighbor)) {
            qh_setdel(vertex->neighbors, neighbor);
            qh_setappend(qh, &newset, neighbor);
            facet= neighbor;
            break;
          }
        }
        if (!neighbor) {
          qh_fprintf(qh, qh->ferr, 6066, "qhull internal error (qh_order_vertexneighbors): no neighbor of v%d for f%d\n",
            vertex->id, facet->id);
          qh_errexit(qh, qh_ERRqhull, facet, NULL);
        }
      }
      qh_setfree(qh, &vertex->neighbors);
      qh_settemppop(qh);
      vertex->neighbors= newset;
    }else { /* qh.hull_dim >= 4 */
      numneighbors= qh_setsize(qh, vertex->neighbors);
      qsort(SETaddr_(vertex->neighbors, facetT), (size_t)numneighbors,
          sizeof(facetT *), qh_compare_facetvisit);
    }
  } /* order_vertexneighbors */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="prepare_output">-</a>
  
    qh_prepare_output(qh )
      prepare for qh_produce_output2(qh) according to
        qh.KEEPminArea, KEEParea, KEEPmerge, GOODvertex, GOODthreshold, GOODpoint, ONLYgood, SPLITthresholds
      does not reset facet->good
  
    notes
      called by qh_produce_output, qh_new_qhull, Qhull.outputQhull
      except for PRINTstatistics, no-op if previously called with same options
  */
  void qh_prepare_output(qhT *qh) {
    if (qh->VORONOI) {
      qh_clearcenters(qh, qh_ASvoronoi);  /* must be before qh_triangulate */
      qh_vertexneighbors(qh);
    }
    if (qh->TRIangulate && !qh->hasTriangulation) {
      qh_triangulate(qh);
      if (qh->VERIFYoutput && !qh->CHECKfrequently)
        qh_checkpolygon(qh, qh->facet_list);
    }
    qh_findgood_all(qh, qh->facet_list);
    if (qh->GETarea)
      qh_getarea(qh, qh->facet_list);
    if (qh->KEEParea || qh->KEEPmerge || qh->KEEPminArea < REALmax/2)
      qh_markkeep(qh, qh->facet_list);
    if (qh->PRINTstatistics)
      qh_collectstatistics(qh);
  }
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printafacet">-</a>
  
    qh_printafacet(qh, 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(qhT *qh, FILE *fp, qh_PRINT 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(qh, facet))
      return;
    if (facet->visible && qh->NEWfacets && format != qh_PRINTfacets)
      return;
    qh->printoutnum++;
    switch (format) {
    case qh_PRINTarea:
      if (facet->isarea) {
        qh_fprintf(qh, fp, 9009, qh_REAL_1, facet->f.area);
        qh_fprintf(qh, fp, 9010, "\n");
      }else
        qh_fprintf(qh, fp, 9011, "0\n");
      break;
    case qh_PRINTcoplanars:
      qh_fprintf(qh, fp, 9012, "%d", qh_setsize(qh, facet->coplanarset));
      FOREACHpoint_(facet->coplanarset)
        qh_fprintf(qh, fp, 9013, " %d", qh_pointid(qh, point));
      qh_fprintf(qh, fp, 9014, "\n");
      break;
    case qh_PRINTcentrums:
      qh_printcenter(qh, fp, format, NULL, facet);
      break;
    case qh_PRINTfacets:
      qh_printfacet(qh, fp, facet);
      break;
    case qh_PRINTfacets_xridge:
      qh_printfacetheader(qh, 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(qh, color, color);
      if (qh->PRINTdim != qh->hull_dim)
        qh_normalize2(qh, color, 3, True, NULL, NULL);
      if (qh->hull_dim <= 2)
        qh_printfacet2geom(qh, fp, facet, color);
      else if (qh->hull_dim == 3) {
        if (facet->simplicial)
          qh_printfacet3geom_simplicial(qh, fp, facet, color);
        else
          qh_printfacet3geom_nonsimplicial(qh, fp, facet, color);
      }else {
        if (facet->simplicial)
          qh_printfacet4geom_simplicial(qh, fp, facet, color);
        else
          qh_printfacet4geom_nonsimplicial(qh, fp, facet, color);
      }
      break;
    case qh_PRINTids:
      qh_fprintf(qh, fp, 9015, "%d\n", facet->id);
      break;
    case qh_PRINTincidences:
    case qh_PRINToff:
    case qh_PRINTtriangles:
      if (qh->hull_dim == 3 && format != qh_PRINTtriangles)
        qh_printfacet3vertex(qh, fp, facet, format);
      else if (facet->simplicial || qh->hull_dim == 2 || format == qh_PRINToff)
        qh_printfacetNvertex_simplicial(qh, fp, facet, format);
      else
        qh_printfacetNvertex_nonsimplicial(qh, fp, facet, qh->printoutvar++, format);
      break;
    case qh_PRINTinner:
      qh_outerinner(qh, facet, NULL, &innerplane);
      offset= facet->offset - innerplane;
      goto LABELprintnorm;
      break; /* prevent warning */
    case qh_PRINTmerges:
      qh_fprintf(qh, fp, 9016, "%d\n", facet->nummerge);
      break;
    case qh_PRINTnormals:
      offset= facet->offset;
      goto LABELprintnorm;
      break; /* prevent warning */
    case qh_PRINTouter:
      qh_outerinner(qh, facet, &outerplane, NULL);
      offset= facet->offset - outerplane;
    LABELprintnorm:
      if (!facet->normal) {
        qh_fprintf(qh, fp, 9017, "no normal for facet f%d\n", facet->id);
        break;
      }
      if (qh->CDDoutput) {
        qh_fprintf(qh, fp, 9018, qh_REAL_1, -offset);
        for (k=0; k < qh->hull_dim; k++)
          qh_fprintf(qh, fp, 9019, qh_REAL_1, -facet->normal[k]);
      }else {
        for (k=0; k < qh->hull_dim; k++)
          qh_fprintf(qh, fp, 9020, qh_REAL_1, facet->normal[k]);
        qh_fprintf(qh, fp, 9021, qh_REAL_1, offset);
      }
      qh_fprintf(qh, fp, 9022, "\n");
      break;
    case qh_PRINTmathematica:  /* either 2 or 3-d by qh_printbegin */
    case qh_PRINTmaple:
      if (qh->hull_dim == 2)
        qh_printfacet2math(qh, fp, facet, format, qh->printoutvar++);
      else
        qh_printfacet3math(qh, fp, facet, format, qh->printoutvar++);
      break;
    case qh_PRINTneighbors:
      qh_fprintf(qh, fp, 9023, "%d", qh_setsize(qh, facet->neighbors));
      FOREACHneighbor_(facet)
        qh_fprintf(qh, fp, 9024, " %d",
                 neighbor->visitid ? neighbor->visitid - 1: 0 - neighbor->id);
      qh_fprintf(qh, fp, 9025, "\n");
      break;
    case qh_PRINTpointintersect:
      if (!qh->feasible_point) {
        qh_fprintf(qh, qh->ferr, 6067, "qhull input error (qh_printafacet): option 'Fp' needs qh->feasible_point\n");
        qh_errexit(qh, qh_ERRinput, NULL, NULL);
      }
      if (facet->offset > 0)
        goto LABELprintinfinite;
      point= coordp= (coordT *)qh_memalloc(qh, 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(qh, point, qh->normal_size);
            goto LABELprintinfinite;
          }
        }
      }
      qh_printpoint(qh, fp, NULL, point);
      qh_memfree(qh, point, qh->normal_size);
      break;
    LABELprintinfinite:
      for (k=qh->hull_dim; k--; )
        qh_fprintf(qh, fp, 9026, qh_REAL_1, qh_INFINITE);
      qh_fprintf(qh, fp, 9027, "\n");
      break;
    case qh_PRINTpointnearest:
      FOREACHpoint_(facet->coplanarset) {
        int id, id2;
        vertex= qh_nearvertex(qh, facet, point, &dist);
        id= qh_pointid(qh, vertex->point);
        id2= qh_pointid(qh, point);
        qh_fprintf(qh, fp, 9028, "%d %d %d " qh_REAL_1 "\n", id, id2, facet->id, dist);
      }
      break;
    case qh_PRINTpoints:  /* VORONOI only by qh_printbegin */
      if (qh->CDDoutput)
        qh_fprintf(qh, fp, 9029, "1 ");
      qh_printcenter(qh, fp, format, NULL, facet);
      break;
    case qh_PRINTvertices:
      qh_fprintf(qh, fp, 9030, "%d", qh_setsize(qh, facet->vertices));
      FOREACHvertex_(facet->vertices)
        qh_fprintf(qh, fp, 9031, " %d", qh_pointid(qh, vertex->point));
      qh_fprintf(qh, fp, 9032, "\n");
      break;
    default:
      break;
    }
  } /* printafacet */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printbegin">-</a>
  
    qh_printbegin(qh )
      prints header for all output formats
  
    returns:
      checks for valid format
  
    notes:
      uses qh.visit_id for 3/4off
      changes qh.interior_point if printing centrums
      qh_countfacets clears facet->visitid for non-good facets
  
    see
      qh_printend() and qh_printafacet()
  
    design:
      count facets and related statistics
      print header for format
  */
  void qh_printbegin(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall) {
    int numfacets, numsimplicial, numridges, totneighbors, numcoplanars, numtricoplanars;
    int i, num;
    facetT *facet, **facetp;
    vertexT *vertex, **vertexp;
    setT *vertices;
    pointT *point, **pointp, *pointtemp;
  
    qh->printoutnum= 0;
    qh_countfacets(qh, facetlist, facets, printall, &numfacets, &numsimplicial,
        &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
    switch (format) {
    case qh_PRINTnone:
      break;
    case qh_PRINTarea:
      qh_fprintf(qh, fp, 9033, "%d\n", numfacets);
      break;
    case qh_PRINTcoplanars:
      qh_fprintf(qh, fp, 9034, "%d\n", numfacets);
      break;
    case qh_PRINTcentrums:
      if (qh->CENTERtype == qh_ASnone)
        qh_clearcenters(qh, qh_AScentrum);
      qh_fprintf(qh, fp, 9035, "%d\n%d\n", qh->hull_dim, numfacets);
      break;
    case qh_PRINTfacets:
    case qh_PRINTfacets_xridge:
      if (facetlist)
        qh_printvertexlist(qh, 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))
        qh_fprintf(qh, qh->ferr, 7049, "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)))
        qh_fprintf(qh, qh->ferr, 7050, "qhull warning: output for outer/inner planes and centrums not implemented in 4-d\n");
      if (qh->PRINTdim == 4 && (qh->PRINTspheres))
        qh_fprintf(qh, qh->ferr, 7051, "qhull warning: output for vertices not implemented in 4-d\n");
      if (qh->PRINTdim == 4 && qh->DOintersections && qh->PRINTnoplanes)
        qh_fprintf(qh, qh->ferr, 7052, "qhull warning: 'Gnh' generates no output in 4-d\n");
      if (qh->PRINTdim == 2) {
        qh_fprintf(qh, fp, 9036, "{appearance {linewidth 3} LIST # %s | %s\n",
                qh->rbox_command, qh->qhull_command);
      }else if (qh->PRINTdim == 3) {
        qh_fprintf(qh, fp, 9037, "{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(qh, NULL, facet, &num, printall);
        FOREACHfacet_(facets)
          qh_printend4geom(qh, NULL, facet, &num, printall);
        qh->ridgeoutnum= num;
        qh->printoutvar= 0;  /* counts number of ridges in output */
        qh_fprintf(qh, fp, 9038, "LIST # %s | %s\n", qh->rbox_command, qh->qhull_command);
      }
  
      if (qh->PRINTdots) {
        qh->printoutnum++;
        num= qh->num_points + qh_setsize(qh, qh->other_points);
        if (qh->DELAUNAY && qh->ATinfinity)
          num--;
        if (qh->PRINTdim == 4)
          qh_fprintf(qh, fp, 9039, "4VECT %d %d 1\n", num, num);
        else
          qh_fprintf(qh, fp, 9040, "VECT %d %d 1\n", num, num);
  
        for (i=num; i--; ) {
          if (i % 20 == 0)
            qh_fprintf(qh, fp, 9041, "\n");
          qh_fprintf(qh, fp, 9042, "1 ");
        }
        qh_fprintf(qh, fp, 9043, "# 1 point per line\n1 ");
        for (i=num-1; i--; ) { /* num at least 3 for D2 */
          if (i % 20 == 0)
            qh_fprintf(qh, fp, 9044, "\n");
          qh_fprintf(qh, fp, 9045, "0 ");
        }
        qh_fprintf(qh, fp, 9046, "# 1 color for all\n");
        FORALLpoints {
          if (!qh->DELAUNAY || !qh->ATinfinity || qh_pointid(qh, point) != qh->num_points-1) {
            if (qh->PRINTdim == 4)
              qh_printpoint(qh, fp, NULL, point);
              else
                qh_printpoint3(qh, fp, point);
          }
        }
        FOREACHpoint_(qh->other_points) {
          if (qh->PRINTdim == 4)
            qh_printpoint(qh, fp, NULL, point);
          else
            qh_printpoint3(qh, fp, point);
        }
        qh_fprintf(qh, fp, 9047, "0 1 1 1  # color of points\n");
      }
  
      if (qh->PRINTdim == 4  && !qh->PRINTnoplanes)
        /* 4dview loads up multiple 4OFF objects slowly */
        qh_fprintf(qh, fp, 9048, "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->MAXabs_coord * qh_MINradius);
      if (qh->premerge_cos < REALmax/2) {
        maximize_(qh->PRINTradius, (1- qh->premerge_cos) * qh->MAXabs_coord);
      }else if (!qh->PREmerge && qh->POSTmerge && qh->postmerge_cos < REALmax/2) {
        maximize_(qh->PRINTradius, (1- qh->postmerge_cos) * qh->MAXabs_coord);
      }
      maximize_(qh->PRINTradius, qh->MINvisible);
      if (qh->JOGGLEmax < REALmax/2)
        qh->PRINTradius += qh->JOGGLEmax * sqrt((realT)qh->hull_dim);
      if (qh->PRINTdim != 4 &&
          (qh->PRINTcoplanar || qh->PRINTspheres || qh->PRINTcentrums)) {
        vertices= qh_facetvertices(qh, facetlist, facets, printall);
        if (qh->PRINTspheres && qh->PRINTdim <= 3)
          qh_printspheres(qh, fp, vertices, qh->PRINTradius);
        if (qh->PRINTcoplanar || qh->PRINTcentrums) {
          qh->firstcentrum= True;
          if (qh->PRINTcoplanar&& !qh->PRINTspheres) {
            FOREACHvertex_(vertices)
              qh_printpointvect2(qh, fp, vertex->point, NULL, qh->interior_point, qh->PRINTradius);
          }
          FORALLfacet_(facetlist) {
            if (!printall && qh_skipfacet(qh, facet))
              continue;
            if (!facet->normal)
              continue;
            if (qh->PRINTcentrums && qh->PRINTdim <= 3)
              qh_printcentrum(qh, fp, facet, qh->PRINTcradius);
            if (!qh->PRINTcoplanar)
              continue;
            FOREACHpoint_(facet->coplanarset)
              qh_printpointvect2(qh, fp, point, facet->normal, NULL, qh->PRINTradius);
            FOREACHpoint_(facet->outsideset)
              qh_printpointvect2(qh, fp, point, facet->normal, NULL, qh->PRINTradius);
          }
          FOREACHfacet_(facets) {
            if (!printall && qh_skipfacet(qh, facet))
              continue;
            if (!facet->normal)
              continue;
            if (qh->PRINTcentrums && qh->PRINTdim <= 3)
              qh_printcentrum(qh, fp, facet, qh->PRINTcradius);
            if (!qh->PRINTcoplanar)
              continue;
            FOREACHpoint_(facet->coplanarset)
              qh_printpointvect2(qh, fp, point, facet->normal, NULL, qh->PRINTradius);
            FOREACHpoint_(facet->outsideset)
              qh_printpointvect2(qh, fp, point, facet->normal, NULL, qh->PRINTradius);
          }
        }
        qh_settempfree(qh, &vertices);
      }
      qh->visit_id++; /* for printing hyperplane intersections */
      break;
    case qh_PRINTids:
      qh_fprintf(qh, fp, 9049, "%d\n", numfacets);
      break;
    case qh_PRINTincidences:
      if (qh->VORONOI && qh->PRINTprecision)
        qh_fprintf(qh, qh->ferr, 7053, "qhull warning: input sites of Delaunay regions (option 'i').  Use option 'p' or 'o' for Voronoi centers.  Disable warning with option 'Pp'\n");
      qh->printoutvar= (int)qh->vertex_id;  /* centrum id for 4-d+, non-simplicial facets */
      if (qh->hull_dim <= 3)
        qh_fprintf(qh, fp, 9050, "%d\n", numfacets);
      else
        qh_fprintf(qh, fp, 9051, "%d\n", numsimplicial+numridges);
      break;
    case qh_PRINTinner:
    case qh_PRINTnormals:
    case qh_PRINTouter:
      if (qh->CDDoutput)
        qh_fprintf(qh, fp, 9052, "%s | %s\nbegin\n    %d %d real\n", qh->rbox_command,
              qh->qhull_command, numfacets, qh->hull_dim+1);
      else
        qh_fprintf(qh, fp, 9053, "%d\n%d\n", qh->hull_dim+1, numfacets);
      break;
    case qh_PRINTmathematica:
    case qh_PRINTmaple:
      if (qh->hull_dim > 3)  /* qh_initbuffers also checks */
        goto LABELnoformat;
      if (qh->VORONOI)
        qh_fprintf(qh, qh->ferr, 7054, "qhull warning: output is the Delaunay triangulation\n");
      if (format == qh_PRINTmaple) {
        if (qh->hull_dim == 2)
          qh_fprintf(qh, fp, 9054, "PLOT(CURVES(\n");
        else
          qh_fprintf(qh, fp, 9055, "PLOT3D(POLYGONS(\n");
      }else
        qh_fprintf(qh, fp, 9056, "{\n");
      qh->printoutvar= 0;   /* counts number of facets for notfirst */
      break;
    case qh_PRINTmerges:
      qh_fprintf(qh, fp, 9057, "%d\n", numfacets);
      break;
    case qh_PRINTpointintersect:
      qh_fprintf(qh, fp, 9058, "%d\n%d\n", qh->hull_dim, numfacets);
      break;
    case qh_PRINTneighbors:
      qh_fprintf(qh, fp, 9059, "%d\n", numfacets);
      break;
    case qh_PRINToff:
    case qh_PRINTtriangles:
      if (qh->VORONOI)
        goto LABELnoformat;
      num= qh->hull_dim;
      if (format == qh_PRINToff || qh->hull_dim == 2)
        qh_fprintf(qh, fp, 9060, "%d\n%d %d %d\n", num,
          qh->num_points+qh_setsize(qh, qh->other_points), numfacets, totneighbors/2);
      else { /* qh_PRINTtriangles */
        qh->printoutvar= qh->num_points+qh_setsize(qh, qh->other_points); /* first centrum */
        if (qh->DELAUNAY)
          num--;  /* drop last dimension */
        qh_fprintf(qh, fp, 9061, "%d\n%d %d %d\n", num, qh->printoutvar
          + numfacets - numsimplicial, numsimplicial + numridges, totneighbors/2);
      }
      FORALLpoints
        qh_printpointid(qh, qh->fout, NULL, num, point, qh_IDunknown);
      FOREACHpoint_(qh->other_points)
        qh_printpointid(qh, qh->fout, NULL, num, point, qh_IDunknown);
      if (format == qh_PRINTtriangles && qh->hull_dim > 2) {
        FORALLfacets {
          if (!facet->simplicial && facet->visitid)
            qh_printcenter(qh, qh->fout, format, NULL, facet);
        }
      }
      break;
    case qh_PRINTpointnearest:
      qh_fprintf(qh, fp, 9062, "%d\n", numcoplanars);
      break;
    case qh_PRINTpoints:
      if (!qh->VORONOI)
        goto LABELnoformat;
      if (qh->CDDoutput)
        qh_fprintf(qh, fp, 9063, "%s | %s\nbegin\n%d %d real\n", qh->rbox_command,
             qh->qhull_command, numfacets, qh->hull_dim);
      else
        qh_fprintf(qh, fp, 9064, "%d\n%d\n", qh->hull_dim-1, numfacets);
      break;
    case qh_PRINTvertices:
      qh_fprintf(qh, fp, 9065, "%d\n", numfacets);
      break;
    case qh_PRINTsummary:
    default:
    LABELnoformat:
      qh_fprintf(qh, qh->ferr, 6068, "qhull internal error (qh_printbegin): can not use this format for dimension %d\n",
           qh->hull_dim);
      qh_errexit(qh, qh_ERRqhull, NULL, NULL);
    }
  } /* printbegin */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printcenter">-</a>
  
    qh_printcenter(qh, fp, string, facet )
      print facet->center as centrum or Voronoi center
      string may be NULL.  Don't include '%' codes.
      nop if qh->CENTERtype neither CENTERvoronoi nor CENTERcentrum
      if upper envelope of Delaunay triangulation and point at-infinity
        prints qh_INFINITE instead;
  
    notes:
      defines facet->center if needed
      if format=PRINTgeom, adds a 0 if would otherwise be 2-d
      Same as QhullFacet::printCenter
  */
  void qh_printcenter(qhT *qh, FILE *fp, qh_PRINT format, const char *string, facetT *facet) {
    int k, num;
  
    if (qh->CENTERtype != qh_ASvoronoi && qh->CENTERtype != qh_AScentrum)
      return;
    if (string)
      qh_fprintf(qh, fp, 9066, string);
    if (qh->CENTERtype == qh_ASvoronoi) {
      num= qh->hull_dim-1;
      if (!facet->normal || !facet->upperdelaunay || !qh->ATinfinity) {
        if (!facet->center)
          facet->center= qh_facetcenter(qh, facet->vertices);
        for (k=0; k < num; k++)
          qh_fprintf(qh, fp, 9067, qh_REAL_1, facet->center[k]);
      }else {
        for (k=0; k < num; k++)
          qh_fprintf(qh, fp, 9068, qh_REAL_1, qh_INFINITE);
      }
    }else /* qh.CENTERtype == qh_AScentrum */ {
      num= qh->hull_dim;
      if (format == qh_PRINTtriangles && qh->DELAUNAY)
        num--;
      if (!facet->center)
        facet->center= qh_getcentrum(qh, facet);
      for (k=0; k < num; k++)
        qh_fprintf(qh, fp, 9069, qh_REAL_1, facet->center[k]);
    }
    if (format == qh_PRINTgeom && num == 2)
      qh_fprintf(qh, fp, 9070, " 0\n");
    else
      qh_fprintf(qh, fp, 9071, "\n");
  } /* printcenter */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printcentrum">-</a>
  
    qh_printcentrum(qh, fp, facet, radius )
      print centrum for a facet in OOGL format
      radius defines size of centrum
      2-d or 3-d only
  
    returns:
      defines facet->center if needed
  */
  void qh_printcentrum(qhT *qh, 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(qh, facet);
      centrum= facet->center;
    }else {
      centrum= qh_getcentrum(qh, facet);
      tempcentrum= True;
    }
    qh_fprintf(qh, fp, 9072, "{appearance {-normal -edge normscale 0} ");
    if (qh->firstcentrum) {
      qh->firstcentrum= False;
      qh_fprintf(qh, fp, 9073, "{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
      qh_fprintf(qh, fp, 9074, "{INST geom { : centrum } transform { # f%d\n", facet->id);
    apex= SETfirstt_(facet->vertices, vertexT);
    qh_distplane(qh, apex->point, facet, &dist);
    projpt= qh_projectpoint(qh, 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(qh, xaxis, xaxis);
      qh_projectdim3(qh, normal, normal);
      qh_normalize2(qh, normal, qh->PRINTdim, True, NULL, NULL);
    }
    qh_crossproduct(3, xaxis, normal, yaxis);
    qh_fprintf(qh, fp, 9075, "%8.4g %8.4g %8.4g 0\n", xaxis[0], xaxis[1], xaxis[2]);
    qh_fprintf(qh, fp, 9076, "%8.4g %8.4g %8.4g 0\n", yaxis[0], yaxis[1], yaxis[2]);
    qh_fprintf(qh, fp, 9077, "%8.4g %8.4g %8.4g 0\n", normal[0], normal[1], normal[2]);
    qh_printpoint3(qh, fp, centrum);
    qh_fprintf(qh, fp, 9078, "1 }}}\n");
    qh_memfree(qh, projpt, qh->normal_size);
    qh_printpointvect(qh, fp, centrum, facet->normal, NULL, radius, green);
    if (tempcentrum)
      qh_memfree(qh, centrum, qh->normal_size);
  } /* printcentrum */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printend">-</a>
  
    qh_printend(qh, fp, format )
      prints trailer for all output formats
  
    see:
      qh_printbegin() and qh_printafacet()
  
  */
  void qh_printend(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall) {
    int num;
    facetT *facet, **facetp;
  
    if (!qh->printoutnum)
      qh_fprintf(qh, qh->ferr, 7055, "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(qh, fp, facet,&num, printall);
        FOREACHfacet_(facets)
          qh_printend4geom(qh, fp, facet, &num, printall);
        if (num != qh->ridgeoutnum || qh->printoutvar != qh->ridgeoutnum) {
          qh_fprintf(qh, qh->ferr, 6069, "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, qh_ERRqhull, NULL, NULL);
        }
      }else
        qh_fprintf(qh, fp, 9079, "}\n");
      break;
    case qh_PRINTinner:
    case qh_PRINTnormals:
    case qh_PRINTouter:
      if (qh->CDDoutput)
        qh_fprintf(qh, fp, 9080, "end\n");
      break;
    case qh_PRINTmaple:
      qh_fprintf(qh, fp, 9081, "));\n");
      break;
    case qh_PRINTmathematica:
      qh_fprintf(qh, fp, 9082, "}\n");
      break;
    case qh_PRINTpoints:
      if (qh->CDDoutput)
        qh_fprintf(qh, fp, 9083, "end\n");
      break;
    default:
      break;
    }
  } /* printend */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printend4geom">-</a>
  
    qh_printend4geom(qh, fp, facet, numridges, printall )
      helper function for qh_printbegin/printend
  
    returns:
      number of printed ridges
  
    notes:
      just counts printed ridges if fp=NULL
      uses facet->visitid
      must agree with qh_printfacet4geom...
  
    design:
      computes color for facet from its normal
      prints each ridge of facet
  */
  void qh_printend4geom(qhT *qh, 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(qh, 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)
            qh_fprintf(qh, fp, 9084, "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)
            qh_fprintf(qh, fp, 9085, "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 */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printextremes">-</a>
  
    qh_printextremes(qh, fp, facetlist, facets, printall )
      print extreme points for convex hulls or halfspace intersections
  
    notes:
      #points, followed by ids, one per line
  
      sorted by id
      same order as qh_printpoints_out if no coplanar/interior points
  */
  void qh_printextremes(qhT *qh, FILE *fp, facetT *facetlist, setT *facets, boolT printall) {
    setT *vertices, *points;
    pointT *point;
    vertexT *vertex, **vertexp;
    int id;
    int numpoints=0, point_i, point_n;
    int allpoints= qh->num_points + qh_setsize(qh, qh->other_points);
  
    points= qh_settemp(qh, allpoints);
    qh_setzero(qh, points, 0, allpoints);
    vertices= qh_facetvertices(qh, facetlist, facets, printall);
    FOREACHvertex_(vertices) {
      id= qh_pointid(qh, vertex->point);
      if (id >= 0) {
        SETelem_(points, id)= vertex->point;
        numpoints++;
      }
    }
    qh_settempfree(qh, &vertices);
    qh_fprintf(qh, fp, 9086, "%d\n", numpoints);
    FOREACHpoint_i_(qh, points) {
      if (point)
        qh_fprintf(qh, fp, 9087, "%d\n", point_i);
    }
    qh_settempfree(qh, &points);
  } /* printextremes */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printextremes_2d">-</a>
  
    qh_printextremes_2d(qh, fp, facetlist, facets, printall )
      prints point ids for facets in qh_ORIENTclock order
  
    notes:
      #points, followed by ids, one per line
      if facetlist/facets are disjoint than the output includes skips
      errors if facets form a loop
      does not print coplanar points
  */
  void qh_printextremes_2d(qhT *qh, FILE *fp, facetT *facetlist, setT *facets, boolT printall) {
    int numfacets, numridges, totneighbors, numcoplanars, numsimplicial, numtricoplanars;
    setT *vertices;
    facetT *facet, *startfacet, *nextfacet;
    vertexT *vertexA, *vertexB;
  
    qh_countfacets(qh, facetlist, facets, printall, &numfacets, &numsimplicial,
        &totneighbors, &numridges, &numcoplanars, &numtricoplanars); /* marks qh->visit_id */
    vertices= qh_facetvertices(qh, facetlist, facets, printall);
    qh_fprintf(qh, fp, 9088, "%d\n", qh_setsize(qh, vertices));
    qh_settempfree(qh, &vertices);
    if (!numfacets)
      return;
    facet= startfacet= facetlist ? facetlist : SETfirstt_(facets, facetT);
    qh->vertex_visit++;
    qh->visit_id++;
    do {
      if (facet->toporient ^ qh_ORIENTclock) {
        vertexA= SETfirstt_(facet->vertices, vertexT);
        vertexB= SETsecondt_(facet->vertices, vertexT);
        nextfacet= SETfirstt_(facet->neighbors, facetT);
      }else {
        vertexA= SETsecondt_(facet->vertices, vertexT);
        vertexB= SETfirstt_(facet->vertices, vertexT);
        nextfacet= SETsecondt_(facet->neighbors, facetT);
      }
      if (facet->visitid == qh->visit_id) {
        qh_fprintf(qh, qh->ferr, 6218, "qhull internal error (qh_printextremes_2d): loop in facet list.  facet %d nextfacet %d\n",
                   facet->id, nextfacet->id);
        qh_errexit2(qh, qh_ERRqhull, facet, nextfacet);
      }
      if (facet->visitid) {
        if (vertexA->visitid != qh->vertex_visit) {
          vertexA->visitid= qh->vertex_visit;
          qh_fprintf(qh, fp, 9089, "%d\n", qh_pointid(qh, vertexA->point));
        }
        if (vertexB->visitid != qh->vertex_visit) {
          vertexB->visitid= qh->vertex_visit;
          qh_fprintf(qh, fp, 9090, "%d\n", qh_pointid(qh, vertexB->point));
        }
      }
      facet->visitid= qh->visit_id;
      facet= nextfacet;
    }while (facet && facet != startfacet);
  } /* printextremes_2d */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printextremes_d">-</a>
  
    qh_printextremes_d(qh, fp, facetlist, facets, printall )
      print extreme points of input sites for Delaunay triangulations
  
    notes:
      #points, followed by ids, one per line
  
      unordered
  */
  void qh_printextremes_d(qhT *qh, FILE *fp, facetT *facetlist, setT *facets, boolT printall) {
    setT *vertices;
    vertexT *vertex, **vertexp;
    boolT upperseen, lowerseen;
    facetT *neighbor, **neighborp;
    int numpoints=0;
  
    vertices= qh_facetvertices(qh, facetlist, facets, printall);
    qh_vertexneighbors(qh);
    FOREACHvertex_(vertices) {
      upperseen= lowerseen= False;
      FOREACHneighbor_(vertex) {
        if (neighbor->upperdelaunay)
          upperseen= True;
        else
          lowerseen= True;
      }
      if (upperseen && lowerseen) {
        vertex->seen= True;
        numpoints++;
      }else
        vertex->seen= False;
    }
    qh_fprintf(qh, fp, 9091, "%d\n", numpoints);
    FOREACHvertex_(vertices) {
      if (vertex->seen)
        qh_fprintf(qh, fp, 9092, "%d\n", qh_pointid(qh, vertex->point));
    }
    qh_settempfree(qh, &vertices);
  } /* printextremes_d */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacet">-</a>
  
    qh_printfacet(qh, fp, facet )
      prints all fields of a facet to fp
  
    notes:
      ridges printed in neighbor order
  */
  void qh_printfacet(qhT *qh, FILE *fp, facetT *facet) {
  
    qh_printfacetheader(qh, fp, facet);
    if (facet->ridges)
      qh_printfacetridges(qh, fp, facet);
  } /* printfacet */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacet2geom">-</a>
  
    qh_printfacet2geom(qh, fp, facet, color )
      print facet as part of a 2-d VECT for Geomview
  
      notes:
        assume precise calculations in io_r.c with roundoff covered by qh_GEOMepsilon
        mindist is calculated within io_r.c.  maxoutside is calculated elsewhere
        so a DISTround error may have occurred.
  */
  void qh_printfacet2geom(qhT *qh, FILE *fp, facetT *facet, realT color[3]) {
    pointT *point0, *point1;
    realT mindist, innerplane, outerplane;
    int k;
  
    qh_facet2point(qh, facet, &point0, &point1, &mindist);
    qh_geomplanes(qh, facet, &outerplane, &innerplane);
    if (qh->PRINTouter || (!qh->PRINTnoplanes && !qh->PRINTinner))
      qh_printfacet2geom_points(qh, fp, point0, point1, facet, outerplane, color);
    if (qh->PRINTinner || (!qh->PRINTnoplanes && !qh->PRINTouter &&
                  outerplane - innerplane > 2 * qh->MAXabs_coord * qh_GEOMepsilon)) {
      for (k=3; k--; )
        color[k]= 1.0 - color[k];
      qh_printfacet2geom_points(qh, fp, point0, point1, facet, innerplane, color);
    }
    qh_memfree(qh, point1, qh->normal_size);
    qh_memfree(qh, point0, qh->normal_size);
  } /* printfacet2geom */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacet2geom_points">-</a>
  
    qh_printfacet2geom_points(qh, fp, point1, point2, facet, offset, color )
      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(qhT *qh, FILE *fp, pointT *point1, pointT *point2,
                                 facetT *facet, realT offset, realT color[3]) {
    pointT *p1= point1, *p2= point2;
  
    qh_fprintf(qh, fp, 9093, "VECT 1 2 1 2 1 # f%d\n", facet->id);
    if (offset != 0.0) {
      p1= qh_projectpoint(qh, p1, facet, -offset);
      p2= qh_projectpoint(qh, p2, facet, -offset);
    }
    qh_fprintf(qh, fp, 9094, "%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(qh, p1, qh->normal_size);
      qh_memfree(qh, p2, qh->normal_size);
    }
    qh_fprintf(qh, fp, 9095, "%8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]);
  } /* printfacet2geom_points */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacet2math">-</a>
  
    qh_printfacet2math(qh, fp, facet, format, notfirst )
      print 2-d Maple or Mathematica output for a facet
      may be non-simplicial
  
    notes:
      use %16.8f since Mathematica 2.2 does not handle exponential format
      see qh_printfacet3math
  */
  void qh_printfacet2math(qhT *qh, FILE *fp, facetT *facet, qh_PRINT format, int notfirst) {
    pointT *point0, *point1;
    realT mindist;
    const char *pointfmt;
  
    qh_facet2point(qh, facet, &point0, &point1, &mindist);
    if (notfirst)
      qh_fprintf(qh, fp, 9096, ",");
    if (format == qh_PRINTmaple)
      pointfmt= "[[%16.8f, %16.8f], [%16.8f, %16.8f]]\n";
    else
      pointfmt= "Line[{{%16.8f, %16.8f}, {%16.8f, %16.8f}}]\n";
    qh_fprintf(qh, fp, 9097, pointfmt, point0[0], point0[1], point1[0], point1[1]);
    qh_memfree(qh, point1, qh->normal_size);
    qh_memfree(qh, point0, qh->normal_size);
  } /* printfacet2math */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacet3geom_nonsimplicial">-</a>
  
    qh_printfacet3geom_nonsimplicial(qh, fp, facet, color )
      print Geomview OFF for a 3-d nonsimplicial facet.
      if DOintersections, prints ridges to unvisited neighbors(qh->visit_id)
  
    notes
      uses facet->visitid for intersections and ridges
  */
  void qh_printfacet3geom_nonsimplicial(qhT *qh, 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, outerplane, innerplane;
    int cntvertices, k;
    realT black[3]={0, 0, 0}, green[3]={0, 1, 0};
  
    qh_geomplanes(qh, facet, &outerplane, &innerplane);
    vertices= qh_facet3vertex(qh, facet); /* oriented */
    cntvertices= qh_setsize(qh, vertices);
    projectedpoints= qh_settemp(qh, cntvertices);
    FOREACHvertex_(vertices) {
      zinc_(Zdistio);
      qh_distplane(qh, vertex->point, facet, &dist);
      projpt= qh_projectpoint(qh, vertex->point, facet, dist);
      qh_setappend(qh, &projectedpoints, projpt);
    }
    if (qh->PRINTouter || (!qh->PRINTnoplanes && !qh->PRINTinner))
      qh_printfacet3geom_points(qh, fp, projectedpoints, facet, outerplane, color);
    if (qh->PRINTinner || (!qh->PRINTnoplanes && !qh->PRINTouter &&
                  outerplane - innerplane > 2 * qh->MAXabs_coord * qh_GEOMepsilon)) {
      for (k=3; k--; )
        color[k]= 1.0 - color[k];
      qh_printfacet3geom_points(qh, fp, projectedpoints, facet, innerplane, color);
    }
    FOREACHpoint_(projectedpoints)
      qh_memfree(qh, point, qh->normal_size);
    qh_settempfree(qh, &projectedpoints);
    qh_settempfree(qh, &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(qh, fp, facet, neighbor, ridge->vertices, black);
          if (qh->PRINTridges) {
            vertexA= SETfirstt_(ridge->vertices, vertexT);
            vertexB= SETsecondt_(ridge->vertices, vertexT);
            qh_printline3geom(qh, fp, vertexA->point, vertexB->point, green);
          }
        }
      }
    }
  } /* printfacet3geom_nonsimplicial */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacet3geom_points">-</a>
  
    qh_printfacet3geom_points(qh, fp, points, facet, offset )
      prints a 3-d facet as OFF Geomview object.
      offset is relative to the facet's hyperplane
      Facet is determined as a list of points
  */
  void qh_printfacet3geom_points(qhT *qh, FILE *fp, setT *points, facetT *facet, realT offset, realT color[3]) {
    int k, n= qh_setsize(qh, points), i;
    pointT *point, **pointp;
    setT *printpoints;
  
    qh_fprintf(qh, fp, 9098, "{ OFF %d 1 1 # f%d\n", n, facet->id);
    if (offset != 0.0) {
      printpoints= qh_settemp(qh, n);
      FOREACHpoint_(points)
        qh_setappend(qh, &printpoints, qh_projectpoint(qh, point, facet, -offset));
    }else
      printpoints= points;
    FOREACHpoint_(printpoints) {
      for (k=0; k < qh->hull_dim; k++) {
        if (k == qh->DROPdim)
          qh_fprintf(qh, fp, 9099, "0 ");
        else
          qh_fprintf(qh, fp, 9100, "%8.4g ", point[k]);
      }
      if (printpoints != points)
        qh_memfree(qh, point, qh->normal_size);
      qh_fprintf(qh, fp, 9101, "\n");
    }
    if (printpoints != points)
      qh_settempfree(qh, &printpoints);
    qh_fprintf(qh, fp, 9102, "%d ", n);
    for (i=0; i < n; i++)
      qh_fprintf(qh, fp, 9103, "%d ", i);
    qh_fprintf(qh, fp, 9104, "%8.4g %8.4g %8.4g 1.0 }\n", color[0], color[1], color[2]);
  } /* printfacet3geom_points */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacet3geom_simplicial">-</a>
  
    qh_printfacet3geom_simplicial(qh )
      print Geomview OFF for a 3-d simplicial facet.
  
    notes:
      may flip color
      uses facet->visitid for intersections and ridges
  
      assume precise calculations in io_r.c with roundoff covered by qh_GEOMepsilon
      innerplane may be off by qh->DISTround.  Maxoutside is calculated elsewhere
      so a DISTround error may have occurred.
  */
  void qh_printfacet3geom_simplicial(qhT *qh, FILE *fp, facetT *facet, realT color[3]) {
    setT *points, *vertices;
    vertexT *vertex, **vertexp, *vertexA, *vertexB;
    facetT *neighbor, **neighborp;
    realT outerplane, innerplane;
    realT black[3]={0, 0, 0}, green[3]={0, 1, 0};
    int k;
  
    qh_geomplanes(qh, facet, &outerplane, &innerplane);
    vertices= qh_facet3vertex(qh, facet);
    points= qh_settemp(qh, qh->TEMPsize);
    FOREACHvertex_(vertices)
      qh_setappend(qh, &points, vertex->point);
    if (qh->PRINTouter || (!qh->PRINTnoplanes && !qh->PRINTinner))
      qh_printfacet3geom_points(qh, fp, points, facet, outerplane, color);
    if (qh->PRINTinner || (!qh->PRINTnoplanes && !qh->PRINTouter &&
                outerplane - innerplane > 2 * qh->MAXabs_coord * qh_GEOMepsilon)) {
      for (k=3; k--; )
        color[k]= 1.0 - color[k];
      qh_printfacet3geom_points(qh, fp, points, facet, innerplane, color);
    }
    qh_settempfree(qh, &points);
    qh_settempfree(qh, &vertices);
    if ((qh->DOintersections || qh->PRINTridges)
    && (!facet->visible || !qh->NEWfacets)) {
      facet->visitid= qh->visit_id;
      FOREACHneighbor_(facet) {
        if (neighbor->visitid != qh->visit_id) {
          vertices= qh_setnew_delnthsorted(qh, facet->vertices, qh->hull_dim,
                            SETindex_(facet->neighbors, neighbor), 0);
          if (qh->DOintersections)
             qh_printhyperplaneintersection(qh, fp, facet, neighbor, vertices, black);
          if (qh->PRINTridges) {
            vertexA= SETfirstt_(vertices, vertexT);
            vertexB= SETsecondt_(vertices, vertexT);
            qh_printline3geom(qh, fp, vertexA->point, vertexB->point, green);
          }
          qh_setfree(qh, &vertices);
        }
      }
    }
  } /* printfacet3geom_simplicial */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacet3math">-</a>
  
    qh_printfacet3math(qh, fp, facet, notfirst )
      print 3-d Maple or Mathematica output for a facet
  
    notes:
      may be non-simplicial
      use %16.8f since Mathematica 2.2 does not handle exponential format
      see qh_printfacet2math
  */
  void qh_printfacet3math(qhT *qh, FILE *fp, facetT *facet, qh_PRINT format, int notfirst) {
    vertexT *vertex, **vertexp;
    setT *points, *vertices;
    pointT *point, **pointp;
    boolT firstpoint= True;
    realT dist;
    const char *pointfmt, *endfmt;
  
    if (notfirst)
      qh_fprintf(qh, fp, 9105, ",\n");
    vertices= qh_facet3vertex(qh, facet);
    points= qh_settemp(qh, qh_setsize(qh, vertices));
    FOREACHvertex_(vertices) {
      zinc_(Zdistio);
      qh_distplane(qh, vertex->point, facet, &dist);
      point= qh_projectpoint(qh, vertex->point, facet, dist);
      qh_setappend(qh, &points, point);
    }
    if (format == qh_PRINTmaple) {
      qh_fprintf(qh, fp, 9106, "[");
      pointfmt= "[%16.8f, %16.8f, %16.8f]";
      endfmt= "]";
    }else {
      qh_fprintf(qh, fp, 9107, "Polygon[{");
      pointfmt= "{%16.8f, %16.8f, %16.8f}";
      endfmt= "}]";
    }
    FOREACHpoint_(points) {
      if (firstpoint)
        firstpoint= False;
      else
        qh_fprintf(qh, fp, 9108, ",\n");
      qh_fprintf(qh, fp, 9109, pointfmt, point[0], point[1], point[2]);
    }
    FOREACHpoint_(points)
      qh_memfree(qh, point, qh->normal_size);
    qh_settempfree(qh, &points);
    qh_settempfree(qh, &vertices);
    qh_fprintf(qh, fp, 9110, "%s", endfmt);
  } /* printfacet3math */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacet3vertex">-</a>
  
    qh_printfacet3vertex(qh, fp, facet, format )
      print vertices in a 3-d facet as point ids
  
    notes:
      prints number of vertices first if format == qh_PRINToff
      the facet may be non-simplicial
  */
  void qh_printfacet3vertex(qhT *qh, FILE *fp, facetT *facet, qh_PRINT format) {
    vertexT *vertex, **vertexp;
    setT *vertices;
  
    vertices= qh_facet3vertex(qh, facet);
    if (format == qh_PRINToff)
      qh_fprintf(qh, fp, 9111, "%d ", qh_setsize(qh, vertices));
    FOREACHvertex_(vertices)
      qh_fprintf(qh, fp, 9112, "%d ", qh_pointid(qh, vertex->point));
    qh_fprintf(qh, fp, 9113, "\n");
    qh_settempfree(qh, &vertices);
  } /* printfacet3vertex */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacet4geom_nonsimplicial">-</a>
  
    qh_printfacet4geom_nonsimplicial(qh )
      print Geomview 4OFF file for a 4d nonsimplicial facet
      prints all ridges to unvisited neighbors (qh.visit_id)
      if qh.DROPdim
        prints in OFF format
  
    notes:
      must agree with printend4geom()
  */
  void qh_printfacet4geom_nonsimplicial(qhT *qh, FILE *fp, facetT *facet, realT color[3]) {
    facetT *neighbor;
    ridgeT *ridge, **ridgep;
    vertexT *vertex, **vertexp;
    pointT *point;
    int k;
    realT dist;
  
    facet->visitid= qh->visit_id;
    if (qh->PRINTnoplanes || (facet->visible && qh->NEWfacets))
      return;
    FOREACHridge_(facet->ridges) {
      neighbor= otherfacet_(ridge, facet);
      if (neighbor->visitid == qh->visit_id)
        continue;
      if (qh->PRINTtransparent && !neighbor->good)
        continue;
      if (qh->DOintersections)
        qh_printhyperplaneintersection(qh, fp, facet, neighbor, ridge->vertices, color);
      else {
        if (qh->DROPdim >= 0)
          qh_fprintf(qh, fp, 9114, "OFF 3 1 1 # f%d\n", facet->id);
        else {
          qh->printoutvar++;
          qh_fprintf(qh, fp, 9115, "# r%d between f%d f%d\n", ridge->id, facet->id, neighbor->id);
        }
        FOREACHvertex_(ridge->vertices) {
          zinc_(Zdistio);
          qh_distplane(qh, vertex->point,facet, &dist);
          point=qh_projectpoint(qh, vertex->point,facet, dist);
          for (k=0; k < qh->hull_dim; k++) {
            if (k != qh->DROPdim)
              qh_fprintf(qh, fp, 9116, "%8.4g ", point[k]);
          }
          qh_fprintf(qh, fp, 9117, "\n");
          qh_memfree(qh, point, qh->normal_size);
        }
        if (qh->DROPdim >= 0)
          qh_fprintf(qh, fp, 9118, "3 0 1 2 %8.4g %8.4g %8.4g\n", color[0], color[1], color[2]);
      }
    }
  } /* printfacet4geom_nonsimplicial */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacet4geom_simplicial">-</a>
  
    qh_printfacet4geom_simplicial(qh, fp, facet, color )
      print Geomview 4OFF file for a 4d simplicial facet
      prints triangles for unvisited neighbors (qh.visit_id)
  
    notes:
      must agree with printend4geom()
  */
  void qh_printfacet4geom_simplicial(qhT *qh, FILE *fp, facetT *facet, realT color[3]) {
    setT *vertices;
    facetT *neighbor, **neighborp;
    vertexT *vertex, **vertexp;
    int k;
  
    facet->visitid= qh->visit_id;
    if (qh->PRINTnoplanes || (facet->visible && qh->NEWfacets))
      return;
    FOREACHneighbor_(facet) {
      if (neighbor->visitid == qh->visit_id)
        continue;
      if (qh->PRINTtransparent && !neighbor->good)
        continue;
      vertices= qh_setnew_delnthsorted(qh, facet->vertices, qh->hull_dim,
                            SETindex_(facet->neighbors, neighbor), 0);
      if (qh->DOintersections)
        qh_printhyperplaneintersection(qh, fp, facet, neighbor, vertices, color);
      else {
        if (qh->DROPdim >= 0)
          qh_fprintf(qh, fp, 9119, "OFF 3 1 1 # ridge between f%d f%d\n",
                  facet->id, neighbor->id);
        else {
          qh->printoutvar++;
          qh_fprintf(qh, fp, 9120, "# ridge between f%d f%d\n", facet->id, neighbor->id);
        }
        FOREACHvertex_(vertices) {
          for (k=0; k < qh->hull_dim; k++) {
            if (k != qh->DROPdim)
              qh_fprintf(qh, fp, 9121, "%8.4g ", vertex->point[k]);
          }
          qh_fprintf(qh, fp, 9122, "\n");
        }
        if (qh->DROPdim >= 0)
          qh_fprintf(qh, fp, 9123, "3 0 1 2 %8.4g %8.4g %8.4g\n", color[0], color[1], color[2]);
      }
      qh_setfree(qh, &vertices);
    }
  } /* printfacet4geom_simplicial */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacetNvertex_nonsimplicial">-</a>
  
    qh_printfacetNvertex_nonsimplicial(qh, fp, facet, id, format )
      print vertices for an N-d non-simplicial facet
      triangulates each ridge to the id
  */
  void qh_printfacetNvertex_nonsimplicial(qhT *qh, FILE *fp, facetT *facet, int id, qh_PRINT format) {
    vertexT *vertex, **vertexp;
    ridgeT *ridge, **ridgep;
  
    if (facet->visible && qh->NEWfacets)
      return;
    FOREACHridge_(facet->ridges) {
      if (format == qh_PRINTtriangles)
        qh_fprintf(qh, fp, 9124, "%d ", qh->hull_dim);
      qh_fprintf(qh, fp, 9125, "%d ", id);
      if ((ridge->top == facet) ^ qh_ORIENTclock) {
        FOREACHvertex_(ridge->vertices)
          qh_fprintf(qh, fp, 9126, "%d ", qh_pointid(qh, vertex->point));
      }else {
        FOREACHvertexreverse12_(ridge->vertices)
          qh_fprintf(qh, fp, 9127, "%d ", qh_pointid(qh, vertex->point));
      }
      qh_fprintf(qh, fp, 9128, "\n");
    }
  } /* printfacetNvertex_nonsimplicial */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacetNvertex_simplicial">-</a>
  
    qh_printfacetNvertex_simplicial(qh, fp, facet, format )
      print vertices for an N-d simplicial facet
      prints vertices for non-simplicial facets
        2-d facets (orientation preserved by qh_mergefacet2d)
        PRINToff ('o') for 4-d and higher
  */
  void qh_printfacetNvertex_simplicial(qhT *qh, FILE *fp, facetT *facet, qh_PRINT format) {
    vertexT *vertex, **vertexp;
  
    if (format == qh_PRINToff || format == qh_PRINTtriangles)
      qh_fprintf(qh, fp, 9129, "%d ", qh_setsize(qh, facet->vertices));
    if ((facet->toporient ^ qh_ORIENTclock)
    || (qh->hull_dim > 2 && !facet->simplicial)) {
      FOREACHvertex_(facet->vertices)
        qh_fprintf(qh, fp, 9130, "%d ", qh_pointid(qh, vertex->point));
    }else {
      FOREACHvertexreverse12_(facet->vertices)
        qh_fprintf(qh, fp, 9131, "%d ", qh_pointid(qh, vertex->point));
    }
    qh_fprintf(qh, fp, 9132, "\n");
  } /* printfacetNvertex_simplicial */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacetheader">-</a>
  
    qh_printfacetheader(qh, fp, facet )
      prints header fields of a facet to fp
  
    notes:
      for 'f' output and debugging
      Same as QhullFacet::printHeader()
  */
  void qh_printfacetheader(qhT *qh, FILE *fp, facetT *facet) {
    pointT *point, **pointp, *furthest;
    facetT *neighbor, **neighborp;
    realT dist;
  
    if (facet == qh_MERGEridge) {
      qh_fprintf(qh, fp, 9133, " MERGEridge\n");
      return;
    }else if (facet == qh_DUPLICATEridge) {
      qh_fprintf(qh, fp, 9134, " DUPLICATEridge\n");
      return;
    }else if (!facet) {
      qh_fprintf(qh, fp, 9135, " NULLfacet\n");
      return;
    }
    qh->old_randomdist= qh->RANDOMdist;
    qh->RANDOMdist= False;
    qh_fprintf(qh, fp, 9136, "- f%d\n", facet->id);
    qh_fprintf(qh, fp, 9137, "    - flags:");
    if (facet->toporient)
      qh_fprintf(qh, fp, 9138, " top");
    else
      qh_fprintf(qh, fp, 9139, " bottom");
    if (facet->simplicial)
      qh_fprintf(qh, fp, 9140, " simplicial");
    if (facet->tricoplanar)
      qh_fprintf(qh, fp, 9141, " tricoplanar");
    if (facet->upperdelaunay)
      qh_fprintf(qh, fp, 9142, " upperDelaunay");
    if (facet->visible)
      qh_fprintf(qh, fp, 9143, " visible");
    if (facet->newfacet)
      qh_fprintf(qh, fp, 9144, " newfacet");
    if (facet->tested)
      qh_fprintf(qh, fp, 9145, " tested");
    if (!facet->good)
      qh_fprintf(qh, fp, 9146, " notG");
    if (facet->seen && qh->IStracing)
      qh_fprintf(qh, fp, 9147, " seen");
    if (facet->seen2 && qh->IStracing)
      qh_fprintf(qh, fp, 9418, " seen2");
    if (facet->isarea)
      qh_fprintf(qh, fp, 9419, " isarea");
    if (facet->coplanarhorizon)
      qh_fprintf(qh, fp, 9148, " coplanarhorizon");
    if (facet->mergehorizon)
      qh_fprintf(qh, fp, 9149, " mergehorizon");
    if (facet->cycledone)
      qh_fprintf(qh, fp, 9420, " cycledone");
    if (facet->keepcentrum)
      qh_fprintf(qh, fp, 9150, " keepcentrum");
    if (facet->dupridge)
      qh_fprintf(qh, fp, 9151, " dupridge");
    if (facet->mergeridge && !facet->mergeridge2)
      qh_fprintf(qh, fp, 9152, " mergeridge1");
    if (facet->mergeridge2)
      qh_fprintf(qh, fp, 9153, " mergeridge2");
    if (facet->newmerge)
      qh_fprintf(qh, fp, 9154, " newmerge");
    if (facet->flipped)
      qh_fprintf(qh, fp, 9155, " flipped");
    if (facet->notfurthest)
      qh_fprintf(qh, fp, 9156, " notfurthest");
    if (facet->degenerate)
      qh_fprintf(qh, fp, 9157, " degenerate");
    if (facet->redundant)
      qh_fprintf(qh, fp, 9158, " redundant");
    qh_fprintf(qh, fp, 9159, "\n");
    if (facet->isarea)
      qh_fprintf(qh, fp, 9160, "    - area: %2.2g\n", facet->f.area);
    else if (qh->NEWfacets && facet->visible && facet->f.replace)
      qh_fprintf(qh, fp, 9161, "    - replacement: f%d\n", facet->f.replace->id);
    else if (facet->newfacet) {
      if (facet->f.samecycle && facet->f.samecycle != facet)
        qh_fprintf(qh, fp, 9162, "    - shares same visible/horizon as f%d\n", facet->f.samecycle->id);
    }else if (facet->tricoplanar /* !isarea */) {
      if (facet->f.triowner)
        qh_fprintf(qh, fp, 9163, "    - owner of normal & centrum is facet f%d\n", facet->f.triowner->id);
    }else if (facet->f.newcycle)
      qh_fprintf(qh, fp, 9164, "    - was horizon to f%d\n", facet->f.newcycle->id);
    if (facet->nummerge == qh_MAXnummerge)
      qh_fprintf(qh, fp, 9427, "    - merges: %dmax\n", qh_MAXnummerge);
    else if (facet->nummerge)
      qh_fprintf(qh, fp, 9165, "    - merges: %d\n", facet->nummerge);
    qh_printpointid(qh, fp, "    - normal: ", qh->hull_dim, facet->normal, qh_IDunknown);
    qh_fprintf(qh, fp, 9166, "    - offset: %10.7g\n", facet->offset);
    if (qh->CENTERtype == qh_ASvoronoi || facet->center)
      qh_printcenter(qh, fp, qh_PRINTfacets, "    - center: ", facet);
  #if qh_MAXoutside
    if (facet->maxoutside > qh->DISTround) /* initial value */
      qh_fprintf(qh, fp, 9167, "    - maxoutside: %10.7g\n", facet->maxoutside);
  #endif
    if (!SETempty_(facet->outsideset)) {
      furthest= (pointT *)qh_setlast(facet->outsideset);
      if (qh_setsize(qh, facet->outsideset) < 6) {
        qh_fprintf(qh, fp, 9168, "    - outside set(furthest p%d):\n", qh_pointid(qh, furthest));
        FOREACHpoint_(facet->outsideset)
          qh_printpoint(qh, fp, "     ", point);
      }else if (qh_setsize(qh, facet->outsideset) < 21) {
        qh_printpoints(qh, fp, "    - outside set:", facet->outsideset);
      }else {
        qh_fprintf(qh, fp, 9169, "    - outside set:  %d points.", qh_setsize(qh, facet->outsideset));
        qh_printpoint(qh, fp, "  Furthest", furthest);
      }
  #if !qh_COMPUTEfurthest
      qh_fprintf(qh, fp, 9170, "    - furthest distance= %2.2g\n", facet->furthestdist);
  #endif
    }
    if (!SETempty_(facet->coplanarset)) {
      furthest= (pointT *)qh_setlast(facet->coplanarset);
      if (qh_setsize(qh, facet->coplanarset) < 6) {
        qh_fprintf(qh, fp, 9171, "    - coplanar set(furthest p%d):\n", qh_pointid(qh, furthest));
        FOREACHpoint_(facet->coplanarset)
          qh_printpoint(qh, fp, "     ", point);
      }else if (qh_setsize(qh, facet->coplanarset) < 21) {
        qh_printpoints(qh, fp, "    - coplanar set:", facet->coplanarset);
      }else {
        qh_fprintf(qh, fp, 9172, "    - coplanar set:  %d points.", qh_setsize(qh, facet->coplanarset));
        qh_printpoint(qh, fp, "  Furthest", furthest);
      }
      zinc_(Zdistio);
      qh_distplane(qh, furthest, facet, &dist);
      qh_fprintf(qh, fp, 9173, "      furthest distance= %2.2g\n", dist);
    }
    qh_printvertices(qh, fp, "    - vertices:", facet->vertices);
    qh_fprintf(qh, fp, 9174, "    - neighboring facets:");
    FOREACHneighbor_(facet) {
      if (neighbor == qh_MERGEridge)
        qh_fprintf(qh, fp, 9175, " MERGEridge");
      else if (neighbor == qh_DUPLICATEridge)
        qh_fprintf(qh, fp, 9176, " DUPLICATEridge");
      else
        qh_fprintf(qh, fp, 9177, " f%d", neighbor->id);
    }
    qh_fprintf(qh, fp, 9178, "\n");
    qh->RANDOMdist= qh->old_randomdist;
  } /* printfacetheader */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacetridges">-</a>
  
    qh_printfacetridges(qh, fp, facet )
      prints ridges of a facet to fp
  
    notes:
      ridges printed in neighbor order
      assumes the ridges exist
      for 'f' output
      same as QhullFacet::printRidges
  */
  void qh_printfacetridges(qhT *qh, FILE *fp, facetT *facet) {
    facetT *neighbor, **neighborp;
    ridgeT *ridge, **ridgep;
    int numridges= 0;
    int n;
  
    if (facet->visible && qh->NEWfacets) {
      qh_fprintf(qh, fp, 9179, "    - ridges (tentative ids):");
      FOREACHridge_(facet->ridges)
        qh_fprintf(qh, fp, 9180, " r%d", ridge->id);
      qh_fprintf(qh, fp, 9181, "\n");
    }else {
      qh_fprintf(qh, fp, 9182, "    - ridges:\n");
      FOREACHridge_(facet->ridges)
        ridge->seen= False;
      if (qh->hull_dim == 3) {
        ridge= SETfirstt_(facet->ridges, ridgeT);
        while (ridge && !ridge->seen) {
          ridge->seen= True;
          qh_printridge(qh, fp, ridge);
          numridges++;
          ridge= qh_nextridge3d(ridge, facet, NULL);
          }
      }else {
        FOREACHneighbor_(facet) {
          FOREACHridge_(facet->ridges) {
            if (otherfacet_(ridge, facet) == neighbor && !ridge->seen) {
              ridge->seen= True;
              qh_printridge(qh, fp, ridge);
              numridges++;
            }
          }
        }
      }
      n= qh_setsize(qh, facet->ridges);
      if (n == 1 && facet->newfacet && qh->NEWtentative) {
        qh_fprintf(qh, fp, 9411, "     - horizon ridge to visible facet\n");
      }
      if (numridges != n) {
        qh_fprintf(qh, fp, 9183, "     - all ridges:");
        FOREACHridge_(facet->ridges)
          qh_fprintf(qh, fp, 9184, " r%d", ridge->id);
        qh_fprintf(qh, fp, 9185, "\n");
      }
      /* non-3d ridges w/o non-simplicial neighbors */
      FOREACHridge_(facet->ridges) {
        if (!ridge->seen)
          qh_printridge(qh, fp, ridge);
      }
    }
  } /* printfacetridges */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printfacets">-</a>
  
    qh_printfacets(qh, fp, format, facetlist, facets, printall )
      prints facetlist and/or facet set in output format
  
    notes:
      also used for specialized formats ('FO' and summary)
      turns off 'Rn' option since want actual numbers
  */
  void qh_printfacets(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall) {
    int numfacets, numsimplicial, numridges, totneighbors, numcoplanars, numtricoplanars;
    facetT *facet, **facetp;
    setT *vertices;
    coordT *center;
    realT outerplane, innerplane;
  
    qh->old_randomdist= qh->RANDOMdist;
    qh->RANDOMdist= False;
    if (qh->CDDoutput && (format == qh_PRINTcentrums || format == qh_PRINTpointintersect || format == qh_PRINToff))
      qh_fprintf(qh, qh->ferr, 7056, "qhull warning: CDD format is not available for centrums, halfspace\nintersections, and OFF file format.\n");
    if (format == qh_PRINTnone)
      ; /* print nothing */
    else if (format == qh_PRINTaverage) {
      vertices= qh_facetvertices(qh, facetlist, facets, printall);
      center= qh_getcenter(qh, vertices);
      qh_fprintf(qh, fp, 9186, "%d 1\n", qh->hull_dim);
      qh_printpointid(qh, fp, NULL, qh->hull_dim, center, qh_IDunknown);
      qh_memfree(qh, center, qh->normal_size);
      qh_settempfree(qh, &vertices);
    }else if (format == qh_PRINTextremes) {
      if (qh->DELAUNAY)
        qh_printextremes_d(qh, fp, facetlist, facets, printall);
      else if (qh->hull_dim == 2)
        qh_printextremes_2d(qh, fp, facetlist, facets, printall);
      else
        qh_printextremes(qh, fp, facetlist, facets, printall);
    }else if (format == qh_PRINToptions)
      qh_fprintf(qh, fp, 9187, "Options selected for Qhull %s:\n%s\n", qh_version, qh->qhull_options);
    else if (format == qh_PRINTpoints && !qh->VORONOI)
      qh_printpoints_out(qh, fp, facetlist, facets, printall);
    else if (format == qh_PRINTqhull)
      qh_fprintf(qh, fp, 9188, "%s | %s\n", qh->rbox_command, qh->qhull_command);
    else if (format == qh_PRINTsize) {
      qh_fprintf(qh, fp, 9189, "0\n2 ");
      qh_fprintf(qh, fp, 9190, qh_REAL_1, qh->totarea);
      qh_fprintf(qh, fp, 9191, qh_REAL_1, qh->totvol);
      qh_fprintf(qh, fp, 9192, "\n");
    }else if (format == qh_PRINTsummary) {
      qh_countfacets(qh, facetlist, facets, printall, &numfacets, &numsimplicial,
        &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
      vertices= qh_facetvertices(qh, facetlist, facets, printall);
      qh_fprintf(qh, fp, 9193, "10 %d %d %d %d %d %d %d %d %d %d\n2 ", qh->hull_dim,
                  qh->num_points + qh_setsize(qh, qh->other_points),
                  qh->num_vertices, qh->num_facets - qh->num_visible,
                  qh_setsize(qh, vertices), numfacets, numcoplanars,
                  numfacets - numsimplicial, zzval_(Zdelvertextot),
                  numtricoplanars);
      qh_settempfree(qh, &vertices);
      qh_outerinner(qh, NULL, &outerplane, &innerplane);
      qh_fprintf(qh, fp, 9194, qh_REAL_2n, outerplane, innerplane);
    }else if (format == qh_PRINTvneighbors)
      qh_printvneighbors(qh, fp, facetlist, facets, printall);
    else if (qh->VORONOI && format == qh_PRINToff)
      qh_printvoronoi(qh, fp, format, facetlist, facets, printall);
    else if (qh->VORONOI && format == qh_PRINTgeom) {
      qh_printbegin(qh, fp, format, facetlist, facets, printall);
      qh_printvoronoi(qh, fp, format, facetlist, facets, printall);
      qh_printend(qh, fp, format, facetlist, facets, printall);
    }else if (qh->VORONOI
    && (format == qh_PRINTvertices || format == qh_PRINTinner || format == qh_PRINTouter))
      qh_printvdiagram(qh, fp, format, facetlist, facets, printall);
    else {
      qh_printbegin(qh, fp, format, facetlist, facets, printall);
      FORALLfacet_(facetlist)
        qh_printafacet(qh, fp, format, facet, printall);
      FOREACHfacet_(facets)
        qh_printafacet(qh, fp, format, facet, printall);
      qh_printend(qh, fp, format, facetlist, facets, printall);
    }
    qh->RANDOMdist= qh->old_randomdist;
  } /* printfacets */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printhyperplaneintersection">-</a>
  
    qh_printhyperplaneintersection(qh, fp, facet1, facet2, vertices, color )
      print Geomview OFF or 4OFF for the intersection of two hyperplanes in 3-d or 4-d
  */
  void qh_printhyperplaneintersection(qhT *qh, FILE *fp, facetT *facet1, facetT *facet2,
                     setT *vertices, realT color[3]) {
    realT costheta, denominator, dist1, dist2, s, t, mindenom, p[4];
    vertexT *vertex, **vertexp;
    int i, k;
    boolT nearzero1, nearzero2;
  
    costheta= qh_getangle(qh, facet1->normal, facet2->normal);
    denominator= 1 - costheta * costheta;
    i= qh_setsize(qh, vertices);
    if (qh->hull_dim == 3)
      qh_fprintf(qh, fp, 9195, "VECT 1 %d 1 %d 1 ", i, i);
    else if (qh->hull_dim == 4 && qh->DROPdim >= 0)
      qh_fprintf(qh, fp, 9196, "OFF 3 1 1 ");
    else
      qh->printoutvar++;
    qh_fprintf(qh, fp, 9197, "# intersect f%d f%d\n", facet1->id, facet2->id);
    mindenom= 1 / (10.0 * qh->MAXabs_coord);
    FOREACHvertex_(vertices) {
      zadd_(Zdistio, 2);
      qh_distplane(qh, vertex->point, facet1, &dist1);
      qh_distplane(qh, 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(qh, p, p);
        qh_fprintf(qh, fp, 9198, "%8.4g %8.4g %8.4g # ", p[0], p[1], p[2]);
      }else
        qh_fprintf(qh, fp, 9199, "%8.4g %8.4g %8.4g %8.4g # ", p[0], p[1], p[2], p[3]);
      if (nearzero1+nearzero2)
        qh_fprintf(qh, fp, 9200, "p%d(coplanar facets)\n", qh_pointid(qh, vertex->point));
      else
        qh_fprintf(qh, fp, 9201, "projected p%d\n", qh_pointid(qh, vertex->point));
    }
    if (qh->hull_dim == 3)
      qh_fprintf(qh, fp, 9202, "%8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]);
    else if (qh->hull_dim == 4 && qh->DROPdim >= 0)
      qh_fprintf(qh, fp, 9203, "3 0 1 2 %8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]);
  } /* printhyperplaneintersection */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printline3geom">-</a>
  
    qh_printline3geom(qh, fp, pointA, pointB, color )
      prints a line as a VECT
      prints 0's for qh.DROPdim
  
    notes:
      if pointA == pointB,
        it's a 1 point VECT
  */
  void qh_printline3geom(qhT *qh, FILE *fp, pointT *pointA, pointT *pointB, realT color[3]) {
    int k;
    realT pA[4], pB[4];
  
    qh_projectdim3(qh, pointA, pA);
    qh_projectdim3(qh, pointB, pB);
    if ((fabs(pA[0] - pB[0]) > 1e-3) ||
        (fabs(pA[1] - pB[1]) > 1e-3) ||
        (fabs(pA[2] - pB[2]) > 1e-3)) {
      qh_fprintf(qh, fp, 9204, "VECT 1 2 1 2 1\n");
      for (k=0; k < 3; k++)
         qh_fprintf(qh, fp, 9205, "%8.4g ", pB[k]);
      qh_fprintf(qh, fp, 9206, " # p%d\n", qh_pointid(qh, pointB));
    }else
      qh_fprintf(qh, fp, 9207, "VECT 1 1 1 1 1\n");
    for (k=0; k < 3; k++)
      qh_fprintf(qh, fp, 9208, "%8.4g ", pA[k]);
    qh_fprintf(qh, fp, 9209, " # p%d\n", qh_pointid(qh, pointA));
    qh_fprintf(qh, fp, 9210, "%8.4g %8.4g %8.4g 1\n", color[0], color[1], color[2]);
  }
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printneighborhood">-</a>
  
    qh_printneighborhood(qh, fp, format, facetA, facetB, printall )
      print neighborhood of one or two facets
  
    notes:
      calls qh_findgood_all()
      bumps qh.visit_id
  */
  void qh_printneighborhood(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetA, facetT *facetB, boolT printall) {
    facetT *neighbor, **neighborp, *facet;
    setT *facets;
  
    if (format == qh_PRINTnone)
      return;
    qh_findgood_all(qh, qh->facet_list);
    if (facetA == facetB)
      facetB= NULL;
    facets= qh_settemp(qh, 2*(qh_setsize(qh, 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(qh, &facets, facet);
      }
      FOREACHneighbor_(facet) {
        if (neighbor->visitid == qh->visit_id)
          continue;
        neighbor->visitid= qh->visit_id;
        if (printall || !qh_skipfacet(qh, neighbor))
          qh_setappend(qh, &facets, neighbor);
      }
    }
    qh_printfacets(qh, fp, format, NULL, facets, printall);
    qh_settempfree(qh, &facets);
  } /* printneighborhood */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printpoint">-</a>
  
    qh_printpoint(qh, fp, string, point )
    qh_printpointid(qh, fp, string, dim, point, id )
      prints the coordinates of a point
  
    returns:
      if string is defined
        prints 'string p%d'.  Skips p%d if id=qh_IDunknown(-1) or qh_IDnone(-3)
  
    notes:
      nop if point is NULL
      Same as QhullPoint's printPoint
  */
  void qh_printpoint(qhT *qh, FILE *fp, const char *string, pointT *point) {
    int id= qh_pointid(qh, point);
  
    qh_printpointid(qh, fp, string, qh->hull_dim, point, id);
  } /* printpoint */
  
  void qh_printpointid(qhT *qh, FILE *fp, const char *string, int dim, pointT *point, int id) {
    int k;
    realT r; /*bug fix*/
  
    if (!point)
      return;
    if (string) {
      qh_fprintf(qh, fp, 9211, "%s", string);
      if (id != qh_IDunknown && id != qh_IDnone)
        qh_fprintf(qh, fp, 9212, " p%d: ", id);
    }
    for (k=dim; k--; ) {
      r= *point++;
      if (string)
        qh_fprintf(qh, fp, 9213, " %8.4g", r);
      else
        qh_fprintf(qh, fp, 9214, qh_REAL_1, r);
    }
    qh_fprintf(qh, fp, 9215, "\n");
  } /* printpointid */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printpoint3">-</a>
  
    qh_printpoint3(qh, fp, point )
      prints 2-d, 3-d, or 4-d point as Geomview 3-d coordinates
  */
  void qh_printpoint3(qhT *qh, FILE *fp, pointT *point) {
    int k;
    realT p[4];
  
    qh_projectdim3(qh, point, p);
    for (k=0; k < 3; k++)
      qh_fprintf(qh, fp, 9216, "%8.4g ", p[k]);
    qh_fprintf(qh, fp, 9217, " # p%d\n", qh_pointid(qh, point));
  } /* printpoint3 */
  
  /*----------------------------------------
  -printpoints- print pointids for a set of points starting at index
     see geom_r.c
  */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printpoints_out">-</a>
  
    qh_printpoints_out(qh, fp, facetlist, facets, printall )
      prints vertices, coplanar/inside points, for facets by their point coordinates
      allows qh.CDDoutput
  
    notes:
      same format as qhull input
      if no coplanar/interior points,
        same order as qh_printextremes
  */
  void qh_printpoints_out(qhT *qh, FILE *fp, facetT *facetlist, setT *facets, boolT printall) {
    int allpoints= qh->num_points + qh_setsize(qh, qh->other_points);
    int numpoints=0, point_i, point_n;
    setT *vertices, *points;
    facetT *facet, **facetp;
    pointT *point, **pointp;
    vertexT *vertex, **vertexp;
    int id;
  
    points= qh_settemp(qh, allpoints);
    qh_setzero(qh, points, 0, allpoints);
    vertices= qh_facetvertices(qh, facetlist, facets, printall);
    FOREACHvertex_(vertices) {
      id= qh_pointid(qh, vertex->point);
      if (id >= 0)
        SETelem_(points, id)= vertex->point;
    }
    if (qh->KEEPinside || qh->KEEPcoplanar || qh->KEEPnearinside) {
      FORALLfacet_(facetlist) {
        if (!printall && qh_skipfacet(qh, facet))
          continue;
        FOREACHpoint_(facet->coplanarset) {
          id= qh_pointid(qh, point);
          if (id >= 0)
            SETelem_(points, id)= point;
        }
      }
      FOREACHfacet_(facets) {
        if (!printall && qh_skipfacet(qh, facet))
          continue;
        FOREACHpoint_(facet->coplanarset) {
          id= qh_pointid(qh, point);
          if (id >= 0)
            SETelem_(points, id)= point;
        }
      }
    }
    qh_settempfree(qh, &vertices);
    FOREACHpoint_i_(qh, points) {
      if (point)
        numpoints++;
    }
    if (qh->CDDoutput)
      qh_fprintf(qh, fp, 9218, "%s | %s\nbegin\n%d %d real\n", qh->rbox_command,
               qh->qhull_command, numpoints, qh->hull_dim + 1);
    else
      qh_fprintf(qh, fp, 9219, "%d\n%d\n", qh->hull_dim, numpoints);
    FOREACHpoint_i_(qh, points) {
      if (point) {
        if (qh->CDDoutput)
          qh_fprintf(qh, fp, 9220, "1 ");
        qh_printpoint(qh, fp, NULL, point);
      }
    }
    if (qh->CDDoutput)
      qh_fprintf(qh, fp, 9221, "end\n");
    qh_settempfree(qh, &points);
  } /* printpoints_out */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printpointvect">-</a>
  
    qh_printpointvect(qh, fp, point, normal, center, radius, color )
      prints a 2-d, 3-d, or 4-d point as 3-d VECT's relative to normal or to center point
  */
  void qh_printpointvect(qhT *qh, FILE *fp, pointT *point, coordT *normal, pointT *center, realT radius, realT color[3]) {
    realT diff[4], pointA[4];
    int k;
  
    for (k=qh->hull_dim; k--; ) {
      if (center)
        diff[k]= point[k]-center[k];
      else if (normal)
        diff[k]= normal[k];
      else
        diff[k]= 0;
    }
    if (center)
      qh_normalize2(qh, diff, qh->hull_dim, True, NULL, NULL);
    for (k=qh->hull_dim; k--; )
      pointA[k]= point[k]+diff[k] * radius;
    qh_printline3geom(qh, fp, point, pointA, color);
  } /* printpointvect */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printpointvect2">-</a>
  
    qh_printpointvect2(qh, fp, point, normal, center, radius )
      prints a 2-d, 3-d, or 4-d point as 2 3-d VECT's for an imprecise point
  */
  void qh_printpointvect2(qhT *qh, FILE *fp, pointT *point, coordT *normal, pointT *center, realT radius) {
    realT red[3]={1, 0, 0}, yellow[3]={1, 1, 0};
  
    qh_printpointvect(qh, fp, point, normal, center, radius, red);
    qh_printpointvect(qh, fp, point, normal, center, -radius, yellow);
  } /* printpointvect2 */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printridge">-</a>
  
    qh_printridge(qh, fp, ridge )
      prints the information in a ridge
  
    notes:
      for qh_printfacetridges()
      same as operator<< [QhullRidge.cpp]
  */
  void qh_printridge(qhT *qh, FILE *fp, ridgeT *ridge) {
  
    qh_fprintf(qh, fp, 9222, "     - r%d", ridge->id);
    if (ridge->tested)
      qh_fprintf(qh, fp, 9223, " tested");
    if (ridge->nonconvex)
      qh_fprintf(qh, fp, 9224, " nonconvex");
    if (ridge->mergevertex)
      qh_fprintf(qh, fp, 9421, " mergevertex");
    if (ridge->mergevertex2)
      qh_fprintf(qh, fp, 9422, " mergevertex2");
    if (ridge->simplicialtop)
      qh_fprintf(qh, fp, 9425, " simplicialtop");
    if (ridge->simplicialbot)
      qh_fprintf(qh, fp, 9423, " simplicialbot");
    qh_fprintf(qh, fp, 9225, "\n");
    qh_printvertices(qh, fp, "           vertices:", ridge->vertices);
    if (ridge->top && ridge->bottom)
      qh_fprintf(qh, fp, 9226, "           between f%d and f%d\n",
              ridge->top->id, ridge->bottom->id);
  } /* printridge */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printspheres">-</a>
  
    qh_printspheres(qh, fp, vertices, radius )
      prints 3-d vertices as OFF spheres
  
    notes:
      inflated octahedron from Stuart Levy earth/mksphere2
  */
  void qh_printspheres(qhT *qh, FILE *fp, setT *vertices, realT radius) {
    vertexT *vertex, **vertexp;
  
    qh->printoutnum++;
    qh_fprintf(qh, fp, 9227, "{appearance {-edge -normal normscale 0} {\n\
  INST geom {define vsphere OFF\n\
  18 32 48\n\
  \n\
  0 0 1\n\
  1 0 0\n\
  0 1 0\n\
  -1 0 0\n\
  0 -1 0\n\
  0 0 -1\n\
  0.707107 0 0.707107\n\
  0 -0.707107 0.707107\n\
  0.707107 -0.707107 0\n\
  -0.707107 0 0.707107\n\
  -0.707107 -0.707107 0\n\
  0 0.707107 0.707107\n\
  -0.707107 0.707107 0\n\
  0.707107 0.707107 0\n\
  0.707107 0 -0.707107\n\
  0 0.707107 -0.707107\n\
  -0.707107 0 -0.707107\n\
  0 -0.707107 -0.707107\n\
  \n\
  3 0 6 11\n\
  3 0 7 6 \n\
  3 0 9 7 \n\
  3 0 11 9\n\
  3 1 6 8 \n\
  3 1 8 14\n\
  3 1 13 6\n\
  3 1 14 13\n\
  3 2 11 13\n\
  3 2 12 11\n\
  3 2 13 15\n\
  3 2 15 12\n\
  3 3 9 12\n\
  3 3 10 9\n\
  3 3 12 16\n\
  3 3 16 10\n\
  3 4 7 10\n\
  3 4 8 7\n\
  3 4 10 17\n\
  3 4 17 8\n\
  3 5 14 17\n\
  3 5 15 14\n\
  3 5 16 15\n\
  3 5 17 16\n\
  3 6 13 11\n\
  3 7 8 6\n\
  3 9 10 7\n\
  3 11 12 9\n\
  3 14 8 17\n\
  3 15 13 14\n\
  3 16 12 15\n\
  3 17 10 16\n} transforms { TLIST\n");
    FOREACHvertex_(vertices) {
      qh_fprintf(qh, fp, 9228, "%8.4g 0 0 0 # v%d\n 0 %8.4g 0 0\n0 0 %8.4g 0\n",
        radius, vertex->id, radius, radius);
      qh_printpoint3(qh, fp, vertex->point);
      qh_fprintf(qh, fp, 9229, "1\n");
    }
    qh_fprintf(qh, fp, 9230, "}}}\n");
  } /* printspheres */
  
  
  /*----------------------------------------------
  -printsummary-
                  see libqhull_r.c
  */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printvdiagram">-</a>
  
    qh_printvdiagram(qh, fp, format, facetlist, facets, printall )
      print voronoi diagram
        # of pairs of input sites
        #indices site1 site2 vertex1 ...
  
      sites indexed by input point id
        point 0 is the first input point
      vertices indexed by 'o' and 'p' order
        vertex 0 is the 'vertex-at-infinity'
        vertex 1 is the first Voronoi vertex
  
    see:
      qh_printvoronoi()
      qh_eachvoronoi_all()
  
    notes:
      if all facets are upperdelaunay,
        prints upper hull (furthest-site Voronoi diagram)
  */
  void qh_printvdiagram(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall) {
    setT *vertices;
    int totcount, numcenters;
    boolT isLower;
    qh_RIDGE innerouter= qh_RIDGEall;
    printvridgeT printvridge= NULL;
  
    if (format == qh_PRINTvertices) {
      innerouter= qh_RIDGEall;
      printvridge= qh_printvridge;
    }else if (format == qh_PRINTinner) {
      innerouter= qh_RIDGEinner;
      printvridge= qh_printvnorm;
    }else if (format == qh_PRINTouter) {
      innerouter= qh_RIDGEouter;
      printvridge= qh_printvnorm;
    }else {
      qh_fprintf(qh, qh->ferr, 6219, "qhull internal error (qh_printvdiagram): unknown print format %d.\n", format);
      qh_errexit(qh, qh_ERRqhull, NULL, NULL);
    }
    vertices= qh_markvoronoi(qh, facetlist, facets, printall, &isLower, &numcenters);
    totcount= qh_printvdiagram2(qh, NULL, NULL, vertices, innerouter, False);
    qh_fprintf(qh, fp, 9231, "%d\n", totcount);
    totcount= qh_printvdiagram2(qh, fp, printvridge, vertices, innerouter, True /* inorder*/);
    qh_settempfree(qh, &vertices);
  #if 0  /* for testing qh_eachvoronoi_all */
    qh_fprintf(qh, fp, 9232, "\n");
    totcount= qh_eachvoronoi_all(qh, fp, printvridge, qh->UPPERdelaunay, innerouter, True /* inorder*/);
    qh_fprintf(qh, fp, 9233, "%d\n", totcount);
  #endif
  } /* printvdiagram */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printvdiagram2">-</a>
  
    qh_printvdiagram2(qh, fp, printvridge, vertices, innerouter, inorder )
      visit all pairs of input sites (vertices) for selected Voronoi vertices
      vertices may include NULLs
  
    innerouter:
      qh_RIDGEall   print inner ridges(bounded) and outer ridges(unbounded)
      qh_RIDGEinner print only inner ridges
      qh_RIDGEouter print only outer ridges
  
    inorder:
      print 3-d Voronoi vertices in order
  
    assumes:
      qh_markvoronoi marked facet->visitid for Voronoi vertices
      all facet->seen= False
      all facet->seen2= True
  
    returns:
      total number of Voronoi ridges
      if printvridge,
        calls printvridge( fp, vertex, vertexA, centers) for each ridge
        [see qh_eachvoronoi()]
  
    see:
      qh_eachvoronoi_all()
  */
  int qh_printvdiagram2(qhT *qh, FILE *fp, printvridgeT printvridge, setT *vertices, qh_RIDGE innerouter, boolT inorder) {
    int totcount= 0;
    int vertex_i, vertex_n;
    vertexT *vertex;
  
    FORALLvertices
      vertex->seen= False;
    FOREACHvertex_i_(qh, vertices) {
      if (vertex) {
        if (qh->GOODvertex > 0 && qh_pointid(qh, vertex->point)+1 != qh->GOODvertex)
          continue;
        totcount += qh_eachvoronoi(qh, fp, printvridge, vertex, !qh_ALL, innerouter, inorder);
      }
    }
    return totcount;
  } /* printvdiagram2 */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printvertex">-</a>
  
    qh_printvertex(qh, fp, vertex )
      prints the information in a vertex
      Duplicated as operator<< [QhullVertex.cpp]
  */
  void qh_printvertex(qhT *qh, FILE *fp, vertexT *vertex) {
    pointT *point;
    int k, count= 0;
    facetT *neighbor, **neighborp;
    realT r; /*bug fix*/
  
    if (!vertex) {
      qh_fprintf(qh, fp, 9234, "  NULLvertex\n");
      return;
    }
    qh_fprintf(qh, fp, 9235, "- p%d(v%d):", qh_pointid(qh, vertex->point), vertex->id);
    point= vertex->point;
    if (point) {
      for (k=qh->hull_dim; k--; ) {
        r= *point++;
        qh_fprintf(qh, fp, 9236, " %5.2g", r);
      }
    }
    if (vertex->deleted)
      qh_fprintf(qh, fp, 9237, " deleted");
    if (vertex->delridge)
      qh_fprintf(qh, fp, 9238, " delridge");
    if (vertex->newfacet)
      qh_fprintf(qh, fp, 9415, " newfacet");
    if (vertex->seen && qh->IStracing)
      qh_fprintf(qh, fp, 9416, " seen");
    if (vertex->seen2 && qh->IStracing)
      qh_fprintf(qh, fp, 9417, " seen2");
    qh_fprintf(qh, fp, 9239, "\n");
    if (vertex->neighbors) {
      qh_fprintf(qh, fp, 9240, "  neighbors:");
      FOREACHneighbor_(vertex) {
        if (++count % 100 == 0)
          qh_fprintf(qh, fp, 9241, "\n     ");
        qh_fprintf(qh, fp, 9242, " f%d", neighbor->id);
      }
      qh_fprintf(qh, fp, 9243, "\n");
    }
  } /* printvertex */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printvertexlist">-</a>
  
    qh_printvertexlist(qh, fp, string, facetlist, facets, printall )
      prints vertices used by a facetlist or facet set
      tests qh_skipfacet() if !printall
  */
  void qh_printvertexlist(qhT *qh, FILE *fp, const char* string, facetT *facetlist,
                           setT *facets, boolT printall) {
    vertexT *vertex, **vertexp;
    setT *vertices;
  
    vertices= qh_facetvertices(qh, facetlist, facets, printall);
    qh_fprintf(qh, fp, 9244, "%s", string);
    FOREACHvertex_(vertices)
      qh_printvertex(qh, fp, vertex);
    qh_settempfree(qh, &vertices);
  } /* printvertexlist */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printvertices">-</a>
  
    qh_printvertices(qh, fp, string, vertices )
      prints vertices in a set
      duplicated as printVertexSet [QhullVertex.cpp]
  */
  void qh_printvertices(qhT *qh, FILE *fp, const char* string, setT *vertices) {
    vertexT *vertex, **vertexp;
  
    qh_fprintf(qh, fp, 9245, "%s", string);
    FOREACHvertex_(vertices)
      qh_fprintf(qh, fp, 9246, " p%d(v%d)", qh_pointid(qh, vertex->point), vertex->id);
    qh_fprintf(qh, fp, 9247, "\n");
  } /* printvertices */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printvneighbors">-</a>
  
    qh_printvneighbors(qh, fp, facetlist, facets, printall )
      print vertex neighbors of vertices in facetlist and facets ('FN')
  
    notes:
      qh_countfacets clears facet->visitid for non-printed facets
  
    design:
      collect facet count and related statistics
      if necessary, build neighbor sets for each vertex
      collect vertices in facetlist and facets
      build a point array for point->vertex and point->coplanar facet
      for each point
        list vertex neighbors or coplanar facet
  */
  void qh_printvneighbors(qhT *qh, FILE *fp, facetT* facetlist, setT *facets, boolT printall) {
    int numfacets, numsimplicial, numridges, totneighbors, numneighbors, numcoplanars, numtricoplanars;
    setT *vertices, *vertex_points, *coplanar_points;
    int numpoints= qh->num_points + qh_setsize(qh, qh->other_points);
    vertexT *vertex, **vertexp;
    int vertex_i, vertex_n;
    facetT *facet, **facetp, *neighbor, **neighborp;
    pointT *point, **pointp;
  
    qh_countfacets(qh, facetlist, facets, printall, &numfacets, &numsimplicial,
        &totneighbors, &numridges, &numcoplanars, &numtricoplanars);  /* sets facet->visitid */
    qh_fprintf(qh, fp, 9248, "%d\n", numpoints);
    qh_vertexneighbors(qh);
    vertices= qh_facetvertices(qh, facetlist, facets, printall);
    vertex_points= qh_settemp(qh, numpoints);
    coplanar_points= qh_settemp(qh, numpoints);
    qh_setzero(qh, vertex_points, 0, numpoints);
    qh_setzero(qh, coplanar_points, 0, numpoints);
    FOREACHvertex_(vertices)
      qh_point_add(qh, vertex_points, vertex->point, vertex);
    FORALLfacet_(facetlist) {
      FOREACHpoint_(facet->coplanarset)
        qh_point_add(qh, coplanar_points, point, facet);
    }
    FOREACHfacet_(facets) {
      FOREACHpoint_(facet->coplanarset)
        qh_point_add(qh, coplanar_points, point, facet);
    }
    FOREACHvertex_i_(qh, vertex_points) {
      if (vertex) {
        numneighbors= qh_setsize(qh, vertex->neighbors);
        qh_fprintf(qh, fp, 9249, "%d", numneighbors);
        qh_order_vertexneighbors(qh, vertex);
        FOREACHneighbor_(vertex)
          qh_fprintf(qh, fp, 9250, " %d",
                   neighbor->visitid ? neighbor->visitid - 1 : 0 - neighbor->id);
        qh_fprintf(qh, fp, 9251, "\n");
      }else if ((facet= SETelemt_(coplanar_points, vertex_i, facetT)))
        qh_fprintf(qh, fp, 9252, "1 %d\n",
                    facet->visitid ? facet->visitid - 1 : 0 - facet->id);
      else
        qh_fprintf(qh, fp, 9253, "0\n");
    }
    qh_settempfree(qh, &coplanar_points);
    qh_settempfree(qh, &vertex_points);
    qh_settempfree(qh, &vertices);
  } /* printvneighbors */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printvoronoi">-</a>
  
    qh_printvoronoi(qh, fp, format, facetlist, facets, printall )
      print voronoi diagram in 'o' or 'G' format
      for 'o' format
        prints voronoi centers for each facet and for infinity
        for each vertex, lists ids of printed facets or infinity
        assumes facetlist and facets are disjoint
      for 'G' format
        prints an OFF object
        adds a 0 coordinate to center
        prints infinity but does not list in vertices
  
    see:
      qh_printvdiagram()
  
    notes:
      if 'o',
        prints a line for each point except "at-infinity"
      if all facets are upperdelaunay,
        reverses lower and upper hull
  */
  void qh_printvoronoi(qhT *qh, FILE *fp, qh_PRINT format, facetT *facetlist, setT *facets, boolT printall) {
    int k, numcenters, numvertices= 0, numneighbors, numinf, vid=1, vertex_i, vertex_n;
    facetT *facet, **facetp, *neighbor, **neighborp;
    setT *vertices;
    vertexT *vertex;
    boolT isLower;
    unsigned int numfacets= (unsigned int)qh->num_facets;
  
    vertices= qh_markvoronoi(qh, facetlist, facets, printall, &isLower, &numcenters);
    FOREACHvertex_i_(qh, vertices) {
      if (vertex) {
        numvertices++;
        numneighbors= numinf= 0;
        FOREACHneighbor_(vertex) {
          if (neighbor->visitid == 0)
            numinf= 1;
          else if (neighbor->visitid < numfacets)
            numneighbors++;
        }
        if (numinf && !numneighbors) {
          SETelem_(vertices, vertex_i)= NULL;
          numvertices--;
        }
      }
    }
    if (format == qh_PRINTgeom)
      qh_fprintf(qh, fp, 9254, "{appearance {+edge -face} OFF %d %d 1 # Voronoi centers and cells\n",
                  numcenters, numvertices);
    else
      qh_fprintf(qh, fp, 9255, "%d\n%d %d 1\n", qh->hull_dim-1, numcenters, qh_setsize(qh, vertices));
    if (format == qh_PRINTgeom) {
      for (k=qh->hull_dim-1; k--; )
        qh_fprintf(qh, fp, 9256, qh_REAL_1, 0.0);
      qh_fprintf(qh, fp, 9257, " 0 # infinity not used\n");
    }else {
      for (k=qh->hull_dim-1; k--; )
        qh_fprintf(qh, fp, 9258, qh_REAL_1, qh_INFINITE);
      qh_fprintf(qh, fp, 9259, "\n");
    }
    FORALLfacet_(facetlist) {
      if (facet->visitid && facet->visitid < numfacets) {
        if (format == qh_PRINTgeom)
          qh_fprintf(qh, fp, 9260, "# %d f%d\n", vid++, facet->id);
        qh_printcenter(qh, fp, format, NULL, facet);
      }
    }
    FOREACHfacet_(facets) {
      if (facet->visitid && facet->visitid < numfacets) {
        if (format == qh_PRINTgeom)
          qh_fprintf(qh, fp, 9261, "# %d f%d\n", vid++, facet->id);
        qh_printcenter(qh, fp, format, NULL, facet);
      }
    }
    FOREACHvertex_i_(qh, vertices) {
      numneighbors= 0;
      numinf=0;
      if (vertex) {
        qh_order_vertexneighbors(qh, vertex);
        FOREACHneighbor_(vertex) {
          if (neighbor->visitid == 0)
            numinf= 1;
          else if (neighbor->visitid < numfacets)
            numneighbors++;
        }
      }
      if (format == qh_PRINTgeom) {
        if (vertex) {
          qh_fprintf(qh, fp, 9262, "%d", numneighbors);
          FOREACHneighbor_(vertex) {
            if (neighbor->visitid && neighbor->visitid < numfacets)
              qh_fprintf(qh, fp, 9263, " %d", neighbor->visitid);
          }
          qh_fprintf(qh, fp, 9264, " # p%d(v%d)\n", vertex_i, vertex->id);
        }else
          qh_fprintf(qh, fp, 9265, " # p%d is coplanar or isolated\n", vertex_i);
      }else {
        if (numinf)
          numneighbors++;
        qh_fprintf(qh, fp, 9266, "%d", numneighbors);
        if (vertex) {
          FOREACHneighbor_(vertex) {
            if (neighbor->visitid == 0) {
              if (numinf) {
                numinf= 0;
                qh_fprintf(qh, fp, 9267, " %d", neighbor->visitid);
              }
            }else if (neighbor->visitid < numfacets)
              qh_fprintf(qh, fp, 9268, " %d", neighbor->visitid);
          }
        }
        qh_fprintf(qh, fp, 9269, "\n");
      }
    }
    if (format == qh_PRINTgeom)
      qh_fprintf(qh, fp, 9270, "}\n");
    qh_settempfree(qh, &vertices);
  } /* printvoronoi */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printvnorm">-</a>
  
    qh_printvnorm(qh, fp, vertex, vertexA, centers, unbounded )
      print one separating plane of the Voronoi diagram for a pair of input sites
      unbounded==True if centers includes vertex-at-infinity
  
    assumes:
      qh_ASvoronoi and qh_vertexneighbors() already set
  
    note:
      parameter unbounded is UNUSED by this callback
  
    see:
      qh_printvdiagram()
      qh_eachvoronoi()
  */
  void qh_printvnorm(qhT *qh, FILE *fp, vertexT *vertex, vertexT *vertexA, setT *centers, boolT unbounded) {
    pointT *normal;
    realT offset;
    int k;
    QHULL_UNUSED(unbounded);
  
    normal= qh_detvnorm(qh, vertex, vertexA, centers, &offset);
    qh_fprintf(qh, fp, 9271, "%d %d %d ",
        2+qh->hull_dim, qh_pointid(qh, vertex->point), qh_pointid(qh, vertexA->point));
    for (k=0; k< qh->hull_dim-1; k++)
      qh_fprintf(qh, fp, 9272, qh_REAL_1, normal[k]);
    qh_fprintf(qh, fp, 9273, qh_REAL_1, offset);
    qh_fprintf(qh, fp, 9274, "\n");
  } /* printvnorm */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="printvridge">-</a>
  
    qh_printvridge(qh, fp, vertex, vertexA, centers, unbounded )
      print one ridge of the Voronoi diagram for a pair of input sites
      unbounded==True if centers includes vertex-at-infinity
  
    see:
      qh_printvdiagram()
  
    notes:
      the user may use a different function
      parameter unbounded is UNUSED
  */
  void qh_printvridge(qhT *qh, FILE *fp, vertexT *vertex, vertexT *vertexA, setT *centers, boolT unbounded) {
    facetT *facet, **facetp;
    QHULL_UNUSED(unbounded);
  
    qh_fprintf(qh, fp, 9275, "%d %d %d", qh_setsize(qh, centers)+2,
         qh_pointid(qh, vertex->point), qh_pointid(qh, vertexA->point));
    FOREACHfacet_(centers)
      qh_fprintf(qh, fp, 9276, " %d", facet->visitid);
    qh_fprintf(qh, fp, 9277, "\n");
  } /* printvridge */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="projectdim3">-</a>
  
    qh_projectdim3(qh, source, destination )
      project 2-d 3-d or 4-d point to a 3-d point
      uses qh.DROPdim and qh.hull_dim
      source and destination may be the same
  
    notes:
      allocate 4 elements to destination just in case
  */
  void qh_projectdim3(qhT *qh, pointT *source, pointT *destination) {
    int i,k;
  
    for (k=0, i=0; k < qh->hull_dim; k++) {
      if (qh->hull_dim == 4) {
        if (k != qh->DROPdim)
          destination[i++]= source[k];
      }else if (k == qh->DROPdim)
        destination[i++]= 0;
      else
        destination[i++]= source[k];
    }
    while (i < 3)
      destination[i++]= 0.0;
  } /* projectdim3 */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="readfeasible">-</a>
  
    qh_readfeasible(qh, dim, curline )
      read feasible point from current line and qh.fin
  
    returns:
      number of lines read from qh.fin
      sets qh.feasible_point with malloc'd coordinates
  
    notes:
      checks for qh.HALFspace
      assumes dim > 1
  
    see:
      qh_setfeasible
  */
  int qh_readfeasible(qhT *qh, int dim, const char *curline) {
    boolT isfirst= True;
    int linecount= 0, tokcount= 0;
    const char *s;
    char *t, firstline[qh_MAXfirst+1];
    coordT *coords, value;
  
    if (!qh->HALFspace) {
      qh_fprintf(qh, qh->ferr, 6070, "qhull input error: feasible point(dim 1 coords) is only valid for halfspace intersection\n");
      qh_errexit(qh, qh_ERRinput, NULL, NULL);
    }
    if (qh->feasible_string)
      qh_fprintf(qh, qh->ferr, 7057, "qhull input warning: feasible point(dim 1 coords) overrides 'Hn,n,n' feasible point for halfspace intersection\n");
    if (!(qh->feasible_point= (coordT *)qh_malloc((size_t)dim * sizeof(coordT)))) {
      qh_fprintf(qh, qh->ferr, 6071, "qhull error: insufficient memory for feasible point\n");
      qh_errexit(qh, qh_ERRmem, NULL, NULL);
    }
    coords= qh->feasible_point;
    while ((s= (isfirst ?  curline : fgets(firstline, qh_MAXfirst, qh->fin)))) {
      if (isfirst)
        isfirst= False;
      else
        linecount++;
      while (*s) {
        while (isspace(*s))
          s++;
        value= qh_strtod(s, &t);
        if (s == t)
          break;
        s= t;
        *(coords++)= value;
        if (++tokcount == dim) {
          while (isspace(*s))
            s++;
          qh_strtod(s, &t);
          if (s != t) {
            qh_fprintf(qh, qh->ferr, 6072, "qhull input error: coordinates for feasible point do not finish out the line: %s\n",
                 s);
            qh_errexit(qh, qh_ERRinput, NULL, NULL);
          }
          return linecount;
        }
      }
    }
    qh_fprintf(qh, qh->ferr, 6073, "qhull input error: only %d coordinates.  Could not read %d-d feasible point.\n",
             tokcount, dim);
    qh_errexit(qh, qh_ERRinput, NULL, NULL);
    return 0;
  } /* readfeasible */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="readpoints">-</a>
  
    qh_readpoints(qh, numpoints, dimension, ismalloc )
      read points from qh.fin into qh.first_point, qh.num_points
      qh.fin is lines of coordinates, one per vertex, first line number of points
      if 'rbox D4',
        gives message
      if qh.ATinfinity,
        adds point-at-infinity for Delaunay triangulations
  
    returns:
      number of points, array of point coordinates, dimension, ismalloc True
      if qh.DELAUNAY & !qh.PROJECTinput, projects points to paraboloid
          and clears qh.PROJECTdelaunay
      if qh.HALFspace, reads optional feasible point, reads halfspaces,
          converts to dual.
  
    for feasible point in "cdd format" in 3-d:
      3 1
      coordinates
      comments
      begin
      n 4 real/integer
      ...
      end
  
    notes:
      dimension will change in qh_initqhull_globals if qh.PROJECTinput
      uses malloc() since qh_mem not initialized
      QH11012 FIX: qh_readpoints needs rewriting, too long
  */
  coordT *qh_readpoints(qhT *qh, int *numpoints, int *dimension, boolT *ismalloc) {
    coordT *points, *coords, *infinity= NULL;
    realT paraboloid, maxboloid= -REALmax, value;
    realT *coordp= NULL, *offsetp= NULL, *normalp= NULL;
    char *s= 0, *t, firstline[qh_MAXfirst+1];
    int diminput=0, numinput=0, dimfeasible= 0, newnum, k, tempi;
    int firsttext=0, firstshort=0, firstlong=0, firstpoint=0;
    int tokcount= 0, linecount=0, maxcount, coordcount=0;
    boolT islong, isfirst= True, wasbegin= False;
    boolT isdelaunay= qh->DELAUNAY && !qh->PROJECTinput;
  
    if (qh->CDDinput) {
      while ((s= fgets(firstline, qh_MAXfirst, qh->fin))) {
        linecount++;
        if (qh->HALFspace && linecount == 1 && isdigit(*s)) {
          dimfeasible= qh_strtol(s, &s);
          while (isspace(*s))
            s++;
          if (qh_strtol(s, &s) == 1)
            linecount += qh_readfeasible(qh, dimfeasible, s);
          else
            dimfeasible= 0;
        }else if (!memcmp(firstline, "begin", (size_t)5) || !memcmp(firstline, "BEGIN", (size_t)5))
          break;
        else if (!*qh->rbox_command)
          strncat(qh->rbox_command, s, sizeof(qh->rbox_command)-1);
      }
      if (!s) {
        qh_fprintf(qh, qh->ferr, 6074, "qhull input error: missing \"begin\" for cdd-formated input\n");
        qh_errexit(qh, qh_ERRinput, NULL, NULL);
      }
    }
    while (!numinput && (s= fgets(firstline, qh_MAXfirst, qh->fin))) {
      linecount++;
      if (!memcmp(s, "begin", (size_t)5) || !memcmp(s, "BEGIN", (size_t)5))
        wasbegin= True;
      while (*s) {
        while (isspace(*s))
          s++;
        if (!*s)
          break;
        if (!isdigit(*s)) {
          if (!*qh->rbox_command) {
            strncat(qh->rbox_command, s, sizeof(qh->rbox_command)-1);
            firsttext= linecount;
          }
          break;
        }
        if (!diminput)
          diminput= qh_strtol(s, &s);
        else {
          numinput= qh_strtol(s, &s);
          if (numinput == 1 && diminput >= 2 && qh->HALFspace && !qh->CDDinput) {
            linecount += qh_readfeasible(qh, diminput, s); /* checks if ok */
            dimfeasible= diminput;
            diminput= numinput= 0;
          }else
            break;
        }
      }
    }
    if (!s) {
      qh_fprintf(qh, qh->ferr, 6075, "qhull input error: short input file.  Did not find dimension and number of points\n");
      qh_errexit(qh, qh_ERRinput, NULL, NULL);
    }
    if (diminput > numinput) {
      tempi= diminput;    /* exchange dim and n, e.g., for cdd input format */
      diminput= numinput;
      numinput= tempi;
    }
    if (diminput < 2) {
      qh_fprintf(qh, qh->ferr, 6220, "qhull input error: dimension %d (first or smaller number) should be at least 2\n",
              diminput);
      qh_errexit(qh, qh_ERRinput, NULL, NULL);
    }
    if (numinput < 1 || numinput > qh_POINTSmax) {
      qh_fprintf(qh, qh->ferr, 6411, "qhull input error: expecting between 1 and %d points.  Got %d %d-d points\n",
        qh_POINTSmax, numinput, diminput);
      qh_errexit(qh, qh_ERRinput, NULL, NULL);
      /* same error message in qh_initqhull_globals */
    }
  
    if (isdelaunay && qh->HALFspace) {
      qh_fprintf(qh, qh->ferr, 6037, "qhull option error (qh_readpoints): can not use Delaunay('d') or Voronoi('v') with halfspace intersection('H')\n");
      qh_errexit(qh, qh_ERRinput, NULL, NULL);
      /* otherwise corrupted memory allocations, same error message as in qh_initqhull_globals */
    }else if (isdelaunay) {
      qh->PROJECTdelaunay= False;
      if (qh->CDDinput)
        *dimension= diminput;
      else
        *dimension= diminput+1;
      *numpoints= numinput;
      if (qh->ATinfinity)
        (*numpoints)++;
    }else if (qh->HALFspace) {
      *dimension= diminput - 1;
      *numpoints= numinput;
      if (diminput < 3) {
        qh_fprintf(qh, qh->ferr, 6221, "qhull input error: dimension %d (first number, includes offset) should be at least 3 for halfspaces\n",
              diminput);
        qh_errexit(qh, qh_ERRinput, NULL, NULL);
      }
      if (dimfeasible) {
        if (dimfeasible != *dimension) {
          qh_fprintf(qh, qh->ferr, 6222, "qhull input error: dimension %d of feasible point is not one less than dimension %d for halfspaces\n",
            dimfeasible, diminput);
          qh_errexit(qh, qh_ERRinput, NULL, NULL);
        }
      }else
        qh_setfeasible(qh, *dimension);
    }else {
      if (qh->CDDinput)
        *dimension= diminput-1;
      else
        *dimension= diminput;
      *numpoints= numinput;
    }
    qh->normal_size= *dimension * (int)sizeof(coordT); /* for tracing with qh_printpoint */
    if (qh->HALFspace) {
      qh->half_space= coordp= (coordT *)qh_malloc((size_t)qh->normal_size + sizeof(coordT));
      if (qh->CDDinput) {
        offsetp= qh->half_space;
        normalp= offsetp + 1;
      }else {
        normalp= qh->half_space;
        offsetp= normalp + *dimension;
      }
    }
    qh->maxline= diminput * (qh_REALdigits + 5);
    maximize_(qh->maxline, 500);
    qh->line= (char *)qh_malloc((size_t)(qh->maxline+1) * sizeof(char));
    *ismalloc= True;  /* use malloc since memory not setup */
    coords= points= qh->temp_malloc=  /* numinput and diminput >=2 by QH6220 */
          (coordT *)qh_malloc((size_t)((*numpoints)*(*dimension))*sizeof(coordT));
    if (!coords || !qh->line || (qh->HALFspace && !qh->half_space)) {
      qh_fprintf(qh, qh->ferr, 6076, "qhull error: insufficient memory to read %d points\n",
              numinput);
      qh_errexit(qh, qh_ERRmem, NULL, NULL);
    }
    if (isdelaunay && qh->ATinfinity) {
      infinity= points + numinput * (*dimension);
      for (k= (*dimension) - 1; k--; )
        infinity[k]= 0.0;
    }
    maxcount= numinput * diminput;
    paraboloid= 0.0;
    while ((s= (isfirst ?  s : fgets(qh->line, qh->maxline, qh->fin)))) {
      if (!isfirst) {
        linecount++;
        if (*s == 'e' || *s == 'E') {
          if (!memcmp(s, "end", (size_t)3) || !memcmp(s, "END", (size_t)3)) {
            if (qh->CDDinput )
              break;
            else if (wasbegin)
              qh_fprintf(qh, qh->ferr, 7058, "qhull input warning: the input appears to be in cdd format.  If so, use 'Fd'\n");
          }
        }
      }
      islong= False;
      while (*s) {
        while (isspace(*s))
          s++;
        value= qh_strtod(s, &t);
        if (s == t) {
          if (!*qh->rbox_command)
           strncat(qh->rbox_command, s, sizeof(qh->rbox_command)-1);
          if (*s && !firsttext)
            firsttext= linecount;
          if (!islong && !firstshort && coordcount)
            firstshort= linecount;
          break;
        }
        if (!firstpoint)
          firstpoint= linecount;
        s= t;
        if (++tokcount > maxcount)
          continue;
        if (qh->HALFspace) {
          if (qh->CDDinput)
            *(coordp++)= -value; /* both coefficients and offset */
          else
            *(coordp++)= value;
        }else {
          *(coords++)= value;
          if (qh->CDDinput && !coordcount) {
            if (value != 1.0) {
              qh_fprintf(qh, qh->ferr, 6077, "qhull input error: for cdd format, point at line %d does not start with '1'\n",
                     linecount);
              qh_errexit(qh, qh_ERRinput, NULL, NULL);
            }
            coords--;
          }else if (isdelaunay) {
            paraboloid += value * value;
            if (qh->ATinfinity) {
              if (qh->CDDinput)
                infinity[coordcount-1] += value;
              else
                infinity[coordcount] += value;
            }
          }
        }
        if (++coordcount == diminput) {
          coordcount= 0;
          if (isdelaunay) {
            *(coords++)= paraboloid;
            maximize_(maxboloid, paraboloid);
            paraboloid= 0.0;
          }else if (qh->HALFspace) {
            if (!qh_sethalfspace(qh, *dimension, coords, &coords, normalp, offsetp, qh->feasible_point)) {
              qh_fprintf(qh, qh->ferr, 8048, "The halfspace was on line %d\n", linecount);
              if (wasbegin)
                qh_fprintf(qh, qh->ferr, 8049, "The input appears to be in cdd format.  If so, you should use option 'Fd'\n");
              qh_errexit(qh, qh_ERRinput, NULL, NULL);
            }
            coordp= qh->half_space;
          }
          while (isspace(*s))
            s++;
          if (*s) {
            islong= True;
            if (!firstlong)
              firstlong= linecount;
          }
        }
      }
      if (!islong && !firstshort && coordcount)
        firstshort= linecount;
      if (!isfirst && s - qh->line >= qh->maxline) {
        qh_fprintf(qh, qh->ferr, 6078, "qhull input error: line %d contained more than %d characters\n",
                linecount, (int) (s - qh->line));   /* WARN64 */
        qh_errexit(qh, qh_ERRinput, NULL, NULL);
      }
      isfirst= False;
    }
    if (qh->rbox_command[0])
      qh->rbox_command[strlen(qh->rbox_command)-1]= '\0'; /* remove \n, previous qh_errexit's display command as two lines */
    if (tokcount != maxcount) {
      newnum= fmin_(numinput, tokcount/diminput);
      if (qh->ALLOWshort)
        qh_fprintf(qh, qh->ferr, 7073, "qhull warning: instead of %d points in %d-d, input contains %d points and %d extra coordinates.\n",
            numinput, diminput, tokcount/diminput, tokcount % diminput);
      else
        qh_fprintf(qh, qh->ferr, 6410, "qhull error: instead of %d points in %d-d, input contains %d points and %d extra coordinates.\n",
            numinput, diminput, tokcount/diminput, tokcount % diminput);
      if (firsttext)
        qh_fprintf(qh, qh->ferr, 8051, "    Line %d is the first comment.\n", firsttext);
      qh_fprintf(qh, qh->ferr, 8033,   "    Line %d is the first point.\n", firstpoint);
      if (firstshort)
        qh_fprintf(qh, qh->ferr, 8052, "    Line %d is the first short line.\n", firstshort);
      if (firstlong)
        qh_fprintf(qh, qh->ferr, 8053, "    Line %d is the first long line.\n", firstlong);
      if (qh->ALLOWshort)
        qh_fprintf(qh, qh->ferr, 8054, "    Continuing with %d points.\n", newnum);
      else {
        qh_fprintf(qh, qh->ferr, 8077, "    Override with option 'Qa' (allow-short)\n");
        qh_errexit(qh, qh_ERRinput, NULL, NULL);
      }
      numinput= newnum;
      if (isdelaunay && qh->ATinfinity) {
        for (k= tokcount % diminput; k--; )
          infinity[k] -= *(--coords);
        *numpoints= newnum+1;
      }else {
        coords -= tokcount % diminput;
        *numpoints= newnum;
      }
    }
    if (isdelaunay && qh->ATinfinity) {
      for (k= (*dimension) - 1; k--; )
        infinity[k] /= numinput;
      if (coords == infinity)
        coords += (*dimension) -1;
      else {
        for (k=0; k < (*dimension) - 1; k++)
          *(coords++)= infinity[k];
      }
      *(coords++)= maxboloid * 1.1;
    }
    if (!strcmp(qh->rbox_command, "./rbox D4"))
      qh_fprintf(qh, qh->ferr, 8055, "\n\
  This is the qhull test case.  If any errors or core dumps occur,\n\
  recompile qhull with 'make new'.  If errors still occur, there is\n\
  an incompatibility.  You should try a different compiler.  You can also\n\
  change the choices in user_r.h.  If you discover the source of the problem,\n\
  please send mail to qhull_bug@qhull.org.\n\
  \n\
  Type 'qhull' for a short list of options.\n");
    qh_free(qh->line);
    qh->line= NULL;
    if (qh->half_space) {
      qh_free(qh->half_space);
      qh->half_space= NULL;
    }
    qh->temp_malloc= NULL;
    trace1((qh, qh->ferr, 1008,"qh_readpoints: read in %d %d-dimensional points\n",
            numinput, diminput));
    return(points);
  } /* readpoints */
  
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="setfeasible">-</a>
  
    qh_setfeasible(qh, dim )
      set qh.feasible_point from qh.feasible_string in "n,n,n" or "n n n" format
  
    notes:
      "n,n,n" already checked by qh_initflags()
      see qh_readfeasible()
      called only once from qh_new_qhull, otherwise leaks memory
  */
  void qh_setfeasible(qhT *qh, int dim) {
    int tokcount= 0;
    char *s;
    coordT *coords, value;
  
    if (!(s= qh->feasible_string)) {
      qh_fprintf(qh, qh->ferr, 6223, "qhull input error: halfspace intersection needs a feasible point.  Either prepend the input with 1 point or use 'Hn,n,n'.  See manual.\n");
      qh_errexit(qh, qh_ERRinput, NULL, NULL);
    }
    if (!(qh->feasible_point= (pointT *)qh_malloc((size_t)dim * sizeof(coordT)))) {
      qh_fprintf(qh, qh->ferr, 6079, "qhull error: insufficient memory for 'Hn,n,n'\n");
      qh_errexit(qh, qh_ERRmem, NULL, NULL);
    }
    coords= qh->feasible_point;
    while (*s) {
      value= qh_strtod(s, &s);
      if (++tokcount > dim) {
        qh_fprintf(qh, qh->ferr, 7059, "qhull input warning: more coordinates for 'H%s' than dimension %d\n",
            qh->feasible_string, dim);
        break;
      }
      *(coords++)= value;
      if (*s)
        s++;
    }
    while (++tokcount <= dim)
      *(coords++)= 0.0;
  } /* setfeasible */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="skipfacet">-</a>
  
    qh_skipfacet(qh, facet )
      returns 'True' if this facet is not to be printed
  
    notes:
      based on the user provided slice thresholds and 'good' specifications
  */
  boolT qh_skipfacet(qhT *qh, facetT *facet) {
    facetT *neighbor, **neighborp;
  
    if (qh->PRINTneighbors) {
      if (facet->good)
        return !qh->PRINTgood;
      FOREACHneighbor_(facet) {
        if (neighbor->good)
          return False;
      }
      return True;
    }else if (qh->PRINTgood)
      return !facet->good;
    else if (!facet->normal)
      return True;
    return(!qh_inthresholds(qh, facet->normal, NULL));
  } /* skipfacet */
  
  /*-<a                             href="qh-io_r.htm#TOC"
    >-------------------------------</a><a name="skipfilename">-</a>
  
    qh_skipfilename(qh, string )
      returns pointer to character after filename
  
    notes:
      skips leading spaces
      ends with spacing or eol
      if starts with ' or " ends with the same, skipping \' or \"
      For qhull, qh_argv_to_command() only uses double quotes
  */
  char *qh_skipfilename(qhT *qh, char *filename) {
    char *s= filename;  /* non-const due to return */
    char c;
  
    while (*s && isspace(*s))
      s++;
    c= *s++;
    if (c == '\0') {
      qh_fprintf(qh, qh->ferr, 6204, "qhull input error: filename expected, none found.\n");
      qh_errexit(qh, qh_ERRinput, NULL, NULL);
    }
    if (c == '\'' || c == '"') {
      while (*s !=c || s[-1] == '\\') {
        if (!*s) {
          qh_fprintf(qh, qh->ferr, 6203, "qhull input error: missing quote after filename -- %s\n", filename);
          qh_errexit(qh, qh_ERRinput, NULL, NULL);
        }
        s++;
      }
      s++;
    }
    else while (*s && !isspace(*s))
        s++;
    return s;
  } /* skipfilename */
  
  