io.c

关于网格剖分的

      }
    }
    if (simplex != points) {
      FOREACHpoint_i_(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 ferr, "qh_detvnorm: points %d %d Voronoi vertex %d dist %2.2g\n",
                             pointid, pointidA, facet->visitid, dist));
        }
      }
    }
  }
  *offsetp= offset;
  if (simplex != points)
    qh_settempfree (&simplex);
  qh_settempfree (&points);
  return normal;
} /* detvnorm */

/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="detvridge">-</a>

  qh_detvridge( 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 (vertexT *vertex) {
  setT *centers= qh_settemp (qh TEMPsize);
  facetT *neighbor, **neighborp;
  boolT firstinf= True;
  
  FOREACHneighbor_(vertex) {
    if (neighbor->seen) {
      if (neighbor->visitid)
        qh_setappend (&centers, neighbor);
      else if (firstinf) {
        firstinf= False;
        qh_setappend (&centers, neighbor);
      }
    }
  }
  qsort (SETaddr_(centers, facetT), qh_setsize (centers),
             sizeof (facetT *), qh_compare_facetvisit);
  return centers;
} /* detvridge */      

/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="detvridge3">-</a>

  qh_detvridge3( 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 (not 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 (vertexT *atvertex, vertexT *vertex) {
  setT *centers= qh_settemp (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)
        qh_setappend (&centers, facet);
      else if (firstinf) {
        firstinf= False;
        qh_setappend (&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) {
	fprintf (stderr, "qh_detvridge3: neigbors of vertex p%d are not connected at facet %d\n",
	         qh_pointid (vertex->point), neighbor->id);
	qh_errexit (qh_ERRqhull, neighbor, NULL);
      }
    }
  }
  FOREACHneighbor_(atvertex) 
    neighbor->seen2= True;
  return centers;
} /* detvridge3 */      

/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="eachvoronoi">-</a>
  
  qh_eachvoronoi( 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*)
        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 (FILE *fp, printvridgeT printvridge, vertexT *atvertex, boolT visitall, qh_RIDGE innerouter, boolT inorder) {
  boolT unbounded;
  int count;
  facetT *neighbor, **neighborp, *neighborA, **neighborAp;
  setT *centers;
  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;
          FOREACHneighborA_(vertex) {
            if (neighborA->seen) {
	      if (neighborA->visitid)
                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 ferr, "qh_eachvoronoi: Voronoi ridge of %d vertices between sites %d and %d\n",
                  count, qh_pointid (atvertex->point), qh_pointid (vertex->point)));
            if (printvridge) { 
	      if (inorder && qh hull_dim == 3+1) /* 3-d Voronoi diagram */
                centers= qh_detvridge3 (atvertex, vertex);
	      else
                centers= qh_detvridge (vertex);
              (*printvridge) (fp, atvertex, vertex, centers, unbounded);
              qh_settempfree (&centers);
            }
          }
        }
      }
    }
  }
  FOREACHneighbor_(atvertex) 
    neighbor->seen= False;
  return totridges;
} /* eachvoronoi */
  

/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="eachvoronoi_all">-</a>
  
  qh_eachvoronoi_all( 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:
    same effect as qh_printvdiagram but ridges not sorted by point id
*/
int qh_eachvoronoi_all (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_ASvoronoi);
  qh_vertexneighbors();
  maximize_(qh visit_id, (unsigned) qh num_facets);
  FORALLfacets {
    facet->visitid= 0;
    facet->seen= False;
    facet->seen2= True;
  }
  FORALLfacets {
    if (facet->upperdelaunay == isupper)
      facet->visitid= numcenters++;
  }
  FORALLvertices 
    vertex->seen= False;
  FORALLvertices 
    totridges += qh_eachvoronoi (fp, printvridge, vertex, 
                    !qh_ALL, innerouter, inorder);
  return totridges;
} /* eachvoronoi_all */
      
/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="facet2point">-</a>
  
  qh_facet2point( 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(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(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 */


/*-<a                             href="qh-c.htm#io"
  >-------------------------------</a><a name="facetvertices">-</a>
  
  qh_facetvertices( 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 (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 num_vertices);
    FORALLvertices {
      vertex->visitid= qh vertex_visit; 
      qh_setappend (&vertices, vertex);
    }
  }else {
    vertices= qh_settemp (qh TEMPsize);
    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 */

/*-<a                             href="qh-c.htm#geom"
  >-------------------------------</a><a name="geomplanes">-</a>
  
  qh_geomplanes( 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.c with roundoff covered by qh_GEOMepsilon
*/
void qh_geomplanes (facetT *facet, realT *outerplane, realT *innerplane) {