io.c

关于网格剖分的

  see:
    qh_printbegin() and qh_printafacet()
      
*/
void qh_printend (FILE *fp, int format, facetT *facetlist, setT *facets, boolT printall) {
  int num;
  facetT *facet, **facetp;

  if (!qh printoutnum)
    fprintf (qh ferr, "qhull warning: no facets printed\n");
  switch (format) {
  case qh_PRINTgeom:
    if (qh hull_dim == 4 && qh DROPdim < 0  && !qh PRINTnoplanes) {
      qh visit_id++;
      num= 0;
      FORALLfacet_(facetlist)
        qh_printend4geom (fp, facet,&num, printall);
      FOREACHfacet_(facets) 
        qh_printend4geom (fp, facet, &num, printall);
      if (num != qh ridgeoutnum || qh printoutvar != qh ridgeoutnum) {
	fprintf (qh ferr, "qhull internal error (qh_printend): number of ridges %d != number printed %d and at end %d\n", qh ridgeoutnum, qh printoutvar, num);
	qh_errexit (qh_ERRqhull, NULL, NULL);
      }
    }else
      fprintf(fp, "}\n");
    break;
  case qh_PRINTinner:
  case qh_PRINTnormals:
  case qh_PRINTouter:
    if (qh CDDoutput) 
      fprintf (fp, "end\n");
    break;
  case qh_PRINTmathematica:
    fprintf(fp, "}\n");
    break;
  case qh_PRINTpoints:
    if (qh CDDoutput)
      fprintf (fp, "end\n");
    break;
  }
} /* printend */

/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="printend4geom">-</a>
  
  qh_printend4geom( 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 (FILE *fp, facetT *facet, int *nump, boolT printall) {
  realT color[3];
  int i, num= *nump;
  facetT *neighbor, **neighborp;
  ridgeT *ridge, **ridgep;
  
  if (!printall && qh_skipfacet(facet))
    return;
  if (qh PRINTnoplanes || (facet->visible && qh NEWfacets))
    return;
  if (!facet->normal)
    return;
  if (fp) {
    for (i=0; i < 3; i++) {
      color[i]= (facet->normal[i]+1.0)/2.0;
      maximize_(color[i], -1.0);
      minimize_(color[i], +1.0);
    }
  }
  facet->visitid= qh visit_id;
  if (facet->simplicial) {
    FOREACHneighbor_(facet) {
      if (neighbor->visitid != qh visit_id) {
	if (fp)
          fprintf (fp, "3 %d %d %d %8.4g %8.4g %8.4g 1 # f%d f%d\n",
		 3*num, 3*num+1, 3*num+2, color[0], color[1], color[2],
		 facet->id, neighbor->id);
	num++;
      }
    }
  }else {
    FOREACHridge_(facet->ridges) {
      neighbor= otherfacet_(ridge, facet);
      if (neighbor->visitid != qh visit_id) {
	if (fp)
          fprintf (fp, "3 %d %d %d %8.4g %8.4g %8.4g 1 #r%d f%d f%d\n",
		 3*num, 3*num+1, 3*num+2, color[0], color[1], color[2],
		 ridge->id, facet->id, neighbor->id);
	num++;
      }
    }
  }
  *nump= num;
} /* printend4geom */

/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="printextremes">-</a>
  
  qh_printextremes( 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 (FILE *fp, facetT *facetlist, setT *facets, int 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 other_points);

  points= qh_settemp (allpoints);
  qh_setzero (points, 0, allpoints);
  vertices= qh_facetvertices (facetlist, facets, printall);
  FOREACHvertex_(vertices) {
    id= qh_pointid (vertex->point);
    if (id >= 0) {
      SETelem_(points, id)= vertex->point;
      numpoints++;
    }
  }
  qh_settempfree (&vertices);
  fprintf (fp, "%d\n", numpoints);
  FOREACHpoint_i_(points) {
    if (point) 
      fprintf (fp, "%d\n", point_i);
  }
  qh_settempfree (&points);
} /* printextremes */

/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="printextremes_2d">-</a>
  
  qh_printextremes_2d( 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 (FILE *fp, facetT *facetlist, setT *facets, int printall) {
  int numfacets, numridges, totneighbors, numcoplanars, numsimplicial;
  setT *vertices;
  facetT *facet, *startfacet, *nextfacet;
  vertexT *vertexA, *vertexB;

  qh_countfacets (facetlist, facets, printall, &numfacets, &numsimplicial, 
      &totneighbors, &numridges, &numcoplanars); /* marks qh visit_id */
  vertices= qh_facetvertices (facetlist, facets, printall);
  fprintf(fp, "%d\n", qh_setsize (vertices));
  qh_settempfree (&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) {
      fprintf(qh ferr, "qh_printextremes_2d: loop in facet list.  facet %d nextfacet %d\n",
                 facet->id, nextfacet->id);
      qh_errexit2 (qh_ERRqhull, facet, nextfacet);
    }
    if (facet->visitid) {
      if (vertexA->visitid != qh vertex_visit) {
	vertexA->visitid= qh vertex_visit;
	fprintf(fp, "%d\n", qh_pointid (vertexA->point));
      }
      if (vertexB->visitid != qh vertex_visit) {
	vertexB->visitid= qh vertex_visit;
	fprintf(fp, "%d\n", qh_pointid (vertexB->point));
      }
    }
    facet->visitid= qh visit_id;
    facet= nextfacet;
  }while (facet && facet != startfacet);
} /* printextremes_2d */

/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="printextremes_d">-</a>
  
  qh_printextremes_d( 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 (FILE *fp, facetT *facetlist, setT *facets, int printall) {
  setT *vertices;
  vertexT *vertex, **vertexp;
  boolT upperseen, lowerseen;
  facetT *neighbor, **neighborp;
  int numpoints=0;

  vertices= qh_facetvertices (facetlist, facets, printall);
  qh_vertexneighbors();
  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;
  }
  fprintf (fp, "%d\n", numpoints);
  FOREACHvertex_(vertices) {
    if (vertex->seen)
      fprintf (fp, "%d\n", qh_pointid (vertex->point));
  }
  qh_settempfree (&vertices);
} /* printextremes_d */

/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="printfacet">-</a>
  
  qh_printfacet( fp, facet )
    prints all fields of a facet to fp

  notes:
    ridges printed in neighbor order
*/
void qh_printfacet(FILE *fp, facetT *facet) {

  qh_printfacetheader (fp, facet);
  if (facet->ridges)
    qh_printfacetridges (fp, facet);
} /* printfacet */


/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="printfacet2geom">-</a>
  
  qh_printfacet2geom( fp, facet, color )
    print facet as part of a 2-d VECT for Geomview
  
    notes:
      assume precise calculations in io.c with roundoff covered by qh_GEOMepsilon
      mindist is calculated within io.c.  maxoutside is calculated elsewhere
      so a DISTround error may have occured.
*/
void qh_printfacet2geom(FILE *fp, facetT *facet, realT color[3]) {
  pointT *point0, *point1;
  realT mindist, innerplane, outerplane;
  int k;

  qh_facet2point (facet, &point0, &point1, &mindist);
  qh_geomplanes (facet, &outerplane, &innerplane);
  if (qh PRINTouter || (!qh PRINTnoplanes && !qh PRINTinner))
    qh_printfacet2geom_points(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(fp, point0, point1, facet, innerplane, color);
  }
  qh_memfree (point1, qh normal_size);
  qh_memfree (point0, qh normal_size); 
} /* printfacet2geom */

/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="printfacet2geom_points">-</a>
  
  qh_printfacet2geom_points( 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(FILE *fp, pointT *point1, pointT *point2,
			       facetT *facet, realT offset, realT color[3]) {
  pointT *p1= point1, *p2= point2;

  fprintf(fp, "VECT 1 2 1 2 1 # f%d\n", facet->id);
  if (offset != 0.0) {
    p1= qh_projectpoint (p1, facet, -offset);
    p2= qh_projectpoint (p2, facet, -offset);
  }
  fprintf(fp, "%8.4g %8.4g %8.4g\n%8.4g %8.4g %8.4g\n",
           p1[0], p1[1], 0.0, p2[0], p2[1], 0.0);
  if (offset != 0.0) {
    qh_memfree (p1, qh normal_size);
    qh_memfree (p2, qh normal_size);
  }
  fprintf(fp, "%8.4g %8.4g %8.4g 1.0\n", color[0], color[1], color[2]);
} /* printfacet2geom_points */


/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="printfacet2math">-</a>
  
  qh_printfacet2math( fp, facet, notfirst )
    print 2-d Mathematica output for a facet
    may be non-simplicial

  notes:
    use %16.8f since Mathematica 2.2 does not handle exponential format
*/
void qh_printfacet2math(FILE *fp, facetT *facet, int notfirst) {
  pointT *point0, *point1;
  realT mindist;
  
  qh_facet2point (facet, &point0, &point1, &mindist);
  if (notfirst)
    fprintf(fp, ",");
  fprintf(fp, "Line[{{%16.8f, %16.8f}, {%16.8f, %16.8f}}]\n",
	  point0[0], point0[1], point1[0], point1[1]);
  qh_memfree (point1, qh normal_size);
  qh_memfree (point0, qh normal_size);
} /* printfacet2math */


/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="printfacet3geom_nonsimplicial">-</a>
  
  qh_printfacet3geom_nonsimplicial( 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(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 (facet, &outerplane, &innerplane); 
  vertices= qh_facet3vertex (facet); /* oriented */
  cntvertices= qh_setsize(vertices);
  projectedpoints= qh_settemp(cntvertices);
  FOREACHvertex_(vertices) {
    zinc_(Zdistio);
    qh_distplane(vertex->point, facet, &dist);
    projpt= qh_projectpoint(vertex->point, facet, dist);
    qh_setappend (&projectedpoints, projpt);
  }
  if (qh PRINTouter || (!qh PRINTnoplanes && !qh PRINTinner))
    qh_printfacet3geom_points(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(fp, projectedpoints, facet, innerplane, color);
  }
  FOREACHpoint_(projectedpoints)
    qh_memfree (point, qh normal_size);
  qh_settempfree(&projectedpoints);
  qh_settempfree(&vertices);
  if ((qh DOintersections || qh PRINTridges)
  && (!facet->visible || !qh NEWfacets)) {
    facet->visitid= qh visit_id;
    FOREACHridge_(facet->ridges) {
      neighbor= otherfacet_(ridge, facet);
      if (neighbor->visitid != qh visit_id) {
        if (qh DOintersections)
          qh_printhyperplaneintersection(fp, facet, neighbor, ridge->vertices, black);
        if (qh PRINTridges) {
          vertexA= SETfirstt_(ridge->vertices, vertexT);
          verte