io.c

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

/* 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 visitid
returns:
  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: