/* io.c - Input/Output routines of qhull application   see README and io.h      see user.c for qh_errprint and qh_printfacetlist   unix.c calls qh_readpoints and qh_produce_output      unix.c and user.c are the only callers of io.c functions   This allows the user to avoid loading io.o from qhull.a   copyright (c) 1993-1995 The Geometry Center        */#include "qhull_a.h"static int qh_compare_facetarea(const void *p1, const void *p2);static int qh_compare_facetmerge(const void *p1, const void *p2);static int qh_compare_facetvisit(const void *p1, const void *p2);int qh_compare_vertexpoint(const void *p1, const void *p2); /* not used *//*========= -functions in alphabetical order after produce_output  ========= *//*--------------------------------------------produce_output- prints out the result of qhull in desired format*/void qh_produce_output(void) {  int i, tempsize= qh_setsize (qhmem.tempstack), d_1;  if (qh VORONOI) {    qh_clearcenters (qh_ASvoronoi);    qh_vertexneighbors();  }  if (qh GETarea)    qh_getarea(qh facet_list);  qh_findgood_all (qh facet_list);   if (qh KEEParea || qh KEEPmerge || qh KEEPminArea < REALmax/2)    qh_markkeep (qh facet_list);  if (qh PRINTsummary)    qh_printsummary(qh ferr);  else if (qh PRINTout[0] == qh_PRINTnone)    qh_printsummary(qh fout);  for (i= 0; i< qh_PRINTEND; i++)    qh_printfacets (qh fout, qh PRINTout[i], qh facet_list, NULL, !qh_ALL);  qh_allstatistics();  if (qh PRINTprecision && !qh MERGING)    qh_printstats (qh ferr, qhstat precision, NULL);  if (qh PRINTstatistics) {    qh_collectstatistics();    qh_printstatistics(qh ferr, "");    qh_memstatistics (qh ferr);    d_1= sizeof(setT) + (qh hull_dim - 1) * SETelemsize;    fprintf(qh ferr, "\    size in bytes: hashentry %d merge %d ridge %d vertex %d facet %d\n\         normal %d ridge vertices %d facet vertices or neighbors %d\n",	    sizeof(hashentryT), sizeof(mergeT), sizeof(ridgeT),	    sizeof(vertexT), sizeof(facetT),	    qh normal_size, d_1, d_1 + SETelemsize);  }  if (qh_setsize (qhmem.tempstack) != tempsize) {    fprintf (qh ferr, "qhull internal error (qh_produce_output): temporary sets not empty (%d)\n",	     qh_setsize (qhmem.tempstack));    qh_errexit (qh_ERRqhull, NULL, NULL);  }} /* produce_output *//*--------------------------------------------------compare_vertexpoint- used by qsort() to order vertices by point id  them*/int qh_compare_vertexpoint(const void *p1, const void *p2) {  vertexT *a= *((vertexT **)p1), *b= *((vertexT **)p2);   return ((qh_pointid(a->point) > qh_pointid(b->point)?1:-1));} /* compare_vertexpoint *//*--------------------------------------------------compare_facetarea- used by qsort() to order facets by area*/static int qh_compare_facetarea(const void *p1, const void *p2) {  facetT *a= *((facetT **)p1), *b= *((facetT **)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 *//*--------------------------------------------------compare_facetmerge- used by qsort() to order facets by number of merges*/static int qh_compare_facetmerge(const void *p1, const void *p2) {  facetT *a= *((facetT **)p1), *b= *((facetT **)p2);   return (a->nummerge - b->nummerge);} /* compare_facetvisit *//*--------------------------------------------------compare_facetvisit- used by qsort() to order facets by visit id or id*/static int qh_compare_facetvisit(const void *p1, const void *p2) {  facetT *a= *((facetT **)p1), *b= *((facetT **)p2);  int i,j;  if (!(i= a->visitid))    i= - a->id;  if (!(j= b->visitid))    j= - b->id;  return (i - j);} /* compare_facetvisit *//*--------------------------------------------------countfacets- count good facets for printing and set visitidreturns:  each facet with ->visitid indicating 1-relative position      ->visitid==0 indicates not good  numfacets, numsimplicial, total neighbors, numridges, coplanars  if NEWfacets, does not count visible facets (matches qh_printafacet)*/void qh_countfacets (facetT *facetlist, setT *facets, boolT printall,    int *numfacetsp, int *numsimplicialp, int *totneighborsp, int *numridgesp, int *numcoplanarsp) {  facetT *facet, **facetp;  int numfacets= 0, numsimplicial= 0, numridges= 0, totneighbors= 0, numcoplanars= 0;  FORALLfacet_(facetlist) {    if ((facet->visible && qh NEWfacets)    || (!printall && qh_skipfacet(facet)))      facet->visitid= 0;    else {      facet->visitid= ++numfacets;      totneighbors += qh_setsize (facet->neighbors);      if (facet->simplicial)         numsimplicial++;      else        numridges += qh_setsize (facet->ridges);      if (facet->coplanarset)        numcoplanars += qh_setsize (facet->coplanarset);    }  }  FOREACHfacet_(facets) {    if ((facet->visible && qh NEWfacets)    || (!printall && qh_skipfacet(facet)))      facet->visitid= 0;    else {      facet->visitid= ++numfacets;      totneighbors += qh_setsize (facet->neighbors);      if (facet->simplicial)        numsimplicial++;      else        numridges += qh_setsize (facet->ridges);      if (facet->coplanarset)        numcoplanars += qh_setsize (facet->coplanarset);    }  }  qh visit_id += numfacets+1;  *numfacetsp= numfacets;  *numsimplicialp= numsimplicial;  *totneighborsp= totneighbors;  *numridgesp= numridges;  *numcoplanarsp= numcoplanars;} /* countfacets *//*----------------------------------------dfacet- print facet by id, for debugging*/void dfacet (int id) {  facetT *facet;  FORALLfacets {    if (facet->id == id) {      qh_printfacet (qh fout, facet);      break;    }  }} /* dfacet *//*----------------------------------------dvertex- print vertex by id, for debugging*/void dvertex (int id) {  vertexT *vertex;  FORALLvertices {    if (vertex->id == id) {      qh_printvertex (qh fout, vertex);      break;    }  }} /* dvertex *//*-----------------------------------------facet2point- return two temporary projected points for a 2-d facet  may be non-simplicial, returns mindist*/void qh_facet2point(facetT *facet, pointT **point0, pointT **point1, realT *mindist) {  vertexT *vertex0, *vertex1;  realT dist;    if (facet->toporient ^ qh_ORIENTclock) {    vertex0= SETfirst_(facet->vertices);    vertex1= SETsecond_(facet->vertices);  }else {    vertex1= SETfirst_(facet->vertices);    vertex0= SETsecond_(facet->vertices);  }  zadd_(Zdistio, 2);  qh_distplane(vertex0->point, facet, &dist);  *mindist= dist;  *point0= qh_projectpoint(vertex0->point, facet, dist);  qh_distplane(vertex1->point, facet, &dist);  minimize_(*mindist, dist);		  *point1= qh_projectpoint(vertex1->point, facet, dist);} /* facet2point *//*--------------------------------------------facetvertices- returns temporary set of vertices in a set or list of facets  optimized for allfacets of facet_list*/setT *qh_facetvertices (facetT *facetlist, setT *facets, boolT allfacets) {  setT *vertices;  facetT *facet, **facetp;  vertexT *vertex, **vertexp;  if (facetlist == qh facet_list && allfacets && !facets) {    vertices= qh_settemp (qh num_vertices);    FORALLvertices       qh_setappend (&vertices, vertex);  }else {    vertices= qh_settemp (qh TEMPsize);    qh vertex_visit++;    FORALLfacet_(facetlist) {      if (!allfacets && qh_skipfacet (facet))        continue;      FOREACHvertex_(facet->vertices) {        if (vertex->visitid != qh vertex_visit) {          vertex->visitid= qh vertex_visit;          qh_setappend (&vertices, vertex);        }      }    }  }  FOREACHfacet_(facets) {    if (!allfacets && qh_skipfacet (facet))      continue;    FOREACHvertex_(facet->vertices) {      if (vertex->visitid != qh vertex_visit) {        vertex->visitid= qh vertex_visit;        qh_setappend (&vertices, vertex);      }    }  }  return vertices;} /* facetvertices *//*------------------------------------------markkeep- mark facets that meet qh KEEParea, qh KEEPmerge, and qh KEEPminArea  clears ->good  recomputes qh num_good*/void qh_markkeep (facetT *facetlist) {  facetT *facet, **facetp;  setT *facets= qh_settemp (qh num_facets);  int size, count;  trace2((qh ferr, "qh_markkeep: only keep %d largest and/or %d most merged facets and/or min area %.2g\n",          qh KEEParea, qh KEEPmerge, qh KEEPminArea));  FORALLfacet_(facetlist) {    if (!facet->visible && facet->good)      qh_setappend (&facets, facet);  }  size= qh_setsize (facets);  if (qh KEEParea) {    qsort (SETaddr_(facets, facetT), size,             sizeof (facetT *), qh_compare_facetarea);    if ((count= size - qh KEEParea) > 0) {      FOREACHfacet_(facets) {        facet->good= False;        if (--count == 0)          break;      }    }  }  if (qh KEEPmerge) {    qsort (SETaddr_(facets, facetT), size,             sizeof (facetT *), qh_compare_facetmerge);    if ((count= size - qh KEEPmerge) > 0) {      FOREACHfacet_(facets) {        facet->good= False;        if (--count == 0)          break;      }    }  }  if (qh KEEPminArea < REALmax/2) {    FOREACHfacet_(facets) {      if (!facet->isarea || facet->f.area < qh KEEPminArea)	facet->good= False;    }  }  qh_settempfree (&facets);  count= 0;  FORALLfacet_(facetlist) {    if (facet->good)      count++;  }  qh num_good= count;} /* markkeep *//*------------------------------------------order_vertexneighbors- order neighbors for a 2-d or 3-d vertex by adjacency  does not orient the neighbors*/void qh_order_vertexneighbors(vertexT *vertex) {  setT *newset;  facetT *facet, *neighbor, **neighborp;  trace4((qh ferr, "qh_order_vertexneighbors: order neighbors of v%d for 3-d\n", vertex->id));  newset= qh_settemp (qh_setsize (vertex->neighbors));  facet= (facetT*)qh_setdellast (vertex->neighbors);  qh_setappend (&newset, facet);  while (qh_setsize (vertex->neighbors)) {    FOREACHneighbor_(vertex) {      if (qh_setin (facet->neighbors, neighbor)) {        qh_setdel(vertex->neighbors, neighbor);        qh_setappend (&newset, neighbor);        facet= neighbor;        break;      }    }    if (!neighbor) {      fprintf (qh ferr, "qhull internal error (qh_order_vertexneighbors): no neighbor of v%d for f%d\n",        vertex->id, facet->id);      qh_errexit (qh_ERRqhull, facet, NULL);    }  }  qh_setfree (&vertex->neighbors);  qh_settemppop ();  vertex->neighbors= newset;} /* order_vertexneighbors *//*------------------------------------------printafacet- print facet in given output format (see qh PRINTout)  nop if skipfacet() unless printall    nop if visible facet and NEWfacets and format != PRINTfacets    must match qh_countfacets  preserves qh visit_id  facet->normal may be null if PREmerge/MERGEexact and STOPcone before merge*/void qh_printafacet(FILE *fp, int format, facetT *facet, boolT printall) {  realT color[4], mindist, offset, dist;  boolT zerodiv;  coordT *point, *normp, *coordp, **pointp, *feasiblep;  int k;  vertexT *vertex, **vertexp;  facetT *neighbor, **neighborp;  if (!printall && qh_skipfacet (facet))    return;  if (facet->visible && qh NEWfacets && format != qh_PRINTfacets)    return;  qh printoutnum++;  switch (format) {  case qh_PRINTarea:    if (facet->isarea) {      fprintf (fp, qh_REAL_1, facet->f.area);      fprintf (fp, "\n");    }    break;  case qh_PRINTcoplanars:    fprintf (fp, "%d", qh_setsize (facet->coplanarset));    FOREACHpoint_(facet->coplanarset)      fprintf (fp, " %d", qh_pointid (point));    fprintf (fp, "\n");    break;  case qh_PRINTcentrums:    qh_printcenter (fp, format, NULL, facet);    break;  case qh_PRINTfacets:    qh_printfacet (fp, facet);    break;  case qh_PRINTfacets_xridge: