poly.c

关于网格剖分的

/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="matchnewfacets">-</a>
  
  qh_matchnewfacets()
    match newfacets in qh.newfacet_list to their newfacet neighbors

  returns:
    qh.newfacet_list with full neighbor sets
      get vertices with nth neighbor by deleting nth vertex
    if qh.PREmerge/MERGEexact or qh.FORCEoutput 
      all facets check for flipped (also prevents point partitioning)
    if duplicate ridges and qh.PREmerge/MERGEexact
      facet->dupridge set
      missing neighbor links identifies extra ridges to be merging

  notes:
    newfacets already have neighbor[0] (horizon facet)
    assumes qh.hash_table is NULL
    vertex->neighbors has not been updated yet
    do not allocate memory after qh.hash_table (need to free it cleanly)

  design:
    delete neighbor sets for all new facets
    initialize a hash table
    for all new facets
      match facet with neighbors
    if unmatched facets (due to duplicate ridges)
      for each new facet with a duplicate ridge
        match it with a facet
    check for flipped facets
*/
void qh_matchnewfacets (void) {
  int numnew=0, hashcount=0, newskip;
  facetT *newfacet, *neighbor;
  int dim= qh hull_dim, hashsize, neighbor_i, neighbor_n;
  setT *neighbors;
#ifndef qh_NOtrace
  int facet_i, facet_n, numfree= 0;
  facetT *facet;
#endif
  
  trace1((qh ferr, "qh_matchnewfacets: match neighbors for new facets.\n"));
  FORALLnew_facets {
    numnew++;
    {  /* inline qh_setzero (newfacet->neighbors, 1, qh hull_dim); */
      neighbors= newfacet->neighbors;
      neighbors->e[neighbors->maxsize].i= dim+1; /*may be overwritten*/
      memset ((char *)SETelemaddr_(neighbors, 1, void), 0, dim * SETelemsize);
    }    
  }
  qh_newhashtable (numnew*(qh hull_dim-1)); /* twice what is normally needed,
                                     but every ridge could be DUPLICATEridge */
  hashsize= qh_setsize (qh hash_table);
  FORALLnew_facets {
    for (newskip=1; newskip<qh hull_dim; newskip++) /* furthest/horizon already matched */
      qh_matchneighbor (newfacet, newskip, hashsize, &hashcount);
#if 0   /* use the following to trap hashcount errors */
    {
      int count= 0, k;
      facetT *facet, *neighbor;

      count= 0;
      FORALLfacet_(qh newfacet_list) {  /* newfacet already in use */
	for (k=1; k < qh hull_dim; k++) {
	  neighbor= SETelemt_(facet->neighbors, k, facetT);
	  if (!neighbor || neighbor == qh_DUPLICATEridge)
	    count++;
	}
	if (facet == newfacet)
	  break;
      }
      if (count != hashcount) {
	fprintf (qh ferr, "qh_matchnewfacets: after adding facet %d, hashcount %d != count %d\n",
		 newfacet->id, hashcount, count);
	qh_errexit (qh_ERRqhull, newfacet, NULL);
      }
    }
#endif  /* end of trap code */
  }
  if (hashcount) {
    FORALLnew_facets {
      if (newfacet->dupridge) {
        FOREACHneighbor_i_(newfacet) {
          if (neighbor == qh_DUPLICATEridge) {
            qh_matchduplicates (newfacet, neighbor_i, hashsize, &hashcount);
         	    /* this may report MERGEfacet */
	  }
        }
      }
    }
  }
  if (hashcount) {
    fprintf (qh ferr, "qhull internal error (qh_matchnewfacets): %d neighbors did not match up\n",
        hashcount);
    qh_printhashtable (qh ferr);
    qh_errexit (qh_ERRqhull, NULL, NULL);
  }
#ifndef qh_NOtrace
  if (qh IStracing >= 2) {
    FOREACHfacet_i_(qh hash_table) {
      if (!facet)
        numfree++;
    }
    fprintf (qh ferr, "qh_matchnewfacets: %d new facets, %d unused hash entries .  hashsize %d\n",
	     numnew, numfree, qh_setsize (qh hash_table));
  }
#endif /* !qh_NOtrace */
  qh_setfree (&qh hash_table);
  if (qh PREmerge || qh MERGEexact) {
    if (qh IStracing >= 4)
      qh_printfacetlist (qh newfacet_list, NULL, qh_ALL);
    FORALLnew_facets {
      if (newfacet->normal)
	qh_checkflipped (newfacet, NULL, qh_ALL);
    }
  }else if (qh FORCEoutput)
    qh_checkflipped_all (qh newfacet_list);  /* prints warnings for flipped */
} /* matchnewfacets */

    
/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="matchvertices">-</a>
  
  qh_matchvertices( firstindex, verticesA, skipA, verticesB, skipB, same )
    tests whether vertices match with a single skip
    starts match at firstindex since all new facets have a common vertex

  returns:
    true if matched vertices
    skip index for each set
    sets same iff vertices have the same orientation

  notes:
    assumes skipA is in A and both sets are the same size

  design:
    set up pointers
    scan both sets checking for a match
    test orientation
*/
boolT qh_matchvertices (int firstindex, setT *verticesA, int skipA, 
       setT *verticesB, int *skipB, boolT *same) {
  vertexT **elemAp, **elemBp, **skipBp=NULL, **skipAp;

  elemAp= SETelemaddr_(verticesA, firstindex, vertexT);
  elemBp= SETelemaddr_(verticesB, firstindex, vertexT);
  skipAp= SETelemaddr_(verticesA, skipA, vertexT);
  do if (elemAp != skipAp) {
    while (*elemAp != *elemBp++) {
      if (skipBp)
        return False;
      skipBp= elemBp;  /* one extra like FOREACH */
    }
  }while(*(++elemAp));
  if (!skipBp)
    skipBp= ++elemBp;
  *skipB= SETindex_(verticesB, skipB);
  *same= !(((ptr_intT)skipA & 0x1) ^ ((ptr_intT)*skipB & 0x1));
  trace4((qh ferr, "qh_matchvertices: matched by skip %d (v%d) and skip %d (v%d) same? %d\n",
	  skipA, (*skipAp)->id, *skipB, (*(skipBp-1))->id, *same));
  return (True);
} /* matchvertices */

/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="newfacet">-</a>
  
  qh_newfacet()
    return a new facet 

  returns:
    all fields initialized or cleared   (NULL)
    preallocates neighbors set
*/
facetT *qh_newfacet(void) {
  facetT *facet;
  void **freelistp; /* used !qh_NOmem */
  
  qh_memalloc_(sizeof(facetT), freelistp, facet, facetT);
  memset ((char *)facet, 0, sizeof(facetT));
  if (qh facet_id == qh tracefacet_id)
    qh tracefacet= facet;
  facet->id= qh facet_id++;
  facet->neighbors= qh_setnew(qh hull_dim);
#if !qh_COMPUTEfurthest
  facet->furthestdist= 0.0;
#endif
#if qh_MAXoutside
  if (qh FORCEoutput && qh APPROXhull)
    facet->maxoutside= qh MINoutside;
  else
    facet->maxoutside= qh DISTround;
#endif
  facet->simplicial= True;
  facet->good= True;
  facet->newfacet= True;
  trace4((qh ferr, "qh_newfacet: created facet f%d\n", facet->id));
  return (facet);
} /* newfacet */


/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="newridge">-</a>
  
  qh_newridge()
    return a new ridge
*/
ridgeT *qh_newridge(void) {
  ridgeT *ridge;
  void **freelistp;   /* used !qh_NOmem */

  qh_memalloc_(sizeof(ridgeT), freelistp, ridge, ridgeT);
  memset ((char *)ridge, 0, sizeof(ridgeT));
  zinc_(Ztotridges);
  if (qh ridge_id == 0xFFFFFF) {
    fprintf(qh ferr, "\
qhull warning: more than %d ridges.  Id field overflows and two ridges\n\
may have the same identifier.  Otherwise output ok.\n", 0xFFFFFF);
  }
  ridge->id= qh ridge_id++;     
  trace4((qh ferr, "qh_newridge: created ridge r%d\n", ridge->id));
  return (ridge);
} /* newridge */


/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="pointid">-</a>
  
  qh_pointid(  )
    return id for a point, 
    returns -3 if null, -2 if interior, or -1 if not known

  alternative code:
    unsigned long id;
    id= ((unsigned long)point - (unsigned long)qh.first_point)/qh.normal_size;

  notes:
    if point not in point array
      the code does a comparison of unrelated pointers.
*/
int qh_pointid (pointT *point) {
  long offset, id;

  if (!point)
    id= -3;
  else if (point == qh interior_point)
    id= -2;
  else if (point >= qh first_point
  && point < qh first_point + qh num_points * qh hull_dim) {
    offset= point - qh first_point;
    id= offset / qh hull_dim;
  }else if ((id= qh_setindex (qh other_points, point)) != -1)
    id += qh num_points;
  else
    id= -1;
  return (int) id;
} /* pointid */
  
/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="removefacet">-</a>
  
  qh_removefacet( facet )
    unlinks facet from qh.facet_list,

  returns:
    updates qh.facet_list .newfacet_list .facet_next visible_list
    decrements qh.num_facets
*/
void qh_removefacet(facetT *facet) {
  facetT *next= facet->next, *previous= facet->previous;
  
  if (facet == qh newfacet_list)
    qh newfacet_list= next;
  if (facet == qh facet_next)
    qh facet_next= next;
  if (facet == qh visible_list)
    qh visible_list= next; 
  if (previous) {
    previous->next= next;
    next->previous= previous;
  }else {  /* 1st facet in qh facet_list */
    qh facet_list= next;
    qh facet_list->previous= NULL;
  }
  qh num_facets--;
  trace4((qh ferr, "qh_removefacet: remove f%d from facet_list\n", facet->id));
} /* removefacet */


/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="removevertex">-</a>
  
  qh_removevertex( vertex )
    unlinks vertex from qh.vertex_list,

  returns:
    updates qh.vertex_list .newvertex_list 
    decrements qh.num_vertices
*/
void qh_removevertex(vertexT *vertex) {
  vertexT *next= vertex->next, *previous= vertex->previous;
  
  if (vertex == qh newvertex_list)
    qh newvertex_list= next;
  if (previous) {
    previous->next= next;
    next->previous= previous;
  }else {  /* 1st vertex in qh vertex_list */
    qh vertex_list= vertex->next;
    qh vertex_list->previous= NULL;
  }
  qh num_vertices--;
  trace4((qh ferr, "qh_removevertex: remove v%d from vertex_list\n", vertex->id));
} /* removevertex */


/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="updatevertices">-</a>
  
  qh_updatevertices()
    update vertex neighbors and delete interior vertices

  returns:
    if qh.VERTEXneighbors
      updates neighbors for each vertex
    interior vertices added to qh.del_vertices for later partitioning

  design:
    if qh.VERTEXneighbors
      deletes references to visible facets from vertex neighbors
      appends new facets to the neighbor list for each vertex
      checks all vertices of visible facets
        removes visible facets from neighbor lists
        marks unused vertices for deletion
*/
void qh_updatevertices (void) {
  facetT *newfacet= NULL, *neighbor, **neighborp, *visible;
  vertexT *vertex, **vertexp;

  trace3((qh ferr, "qh_updatevertices: delete interior vertices and update vertex->neighbors\n"));
  if (qh VERTEXneighbors) {
    FORALLvertex_(qh newvertex_list) {
      FOREACHneighbor_(vertex) {
	if (neighbor->visible) 
	  SETref_(neighbor)= NULL;
      }
      qh_setcompact (vertex->neighbors);
    }
    FORALLnew_facets {
      FOREACHvertex_(newfacet->vertices)
        qh_setappend (&vertex->neighbors, newfacet);
    }
    FORALLvisible_facets {
      FOREACHvertex_(visible->vertices) {
        if (!vertex->newlist && !vertex->deleted) {
  	  FOREACHneighbor_(vertex) { /* this can happen under merging */
	    if (!neighbor->visible)
	      break;
	  }
	  if (neighbor)
	    qh_setdel (vertex->neighbors, visible);
	  else {
	    vertex->deleted= True;
	    qh_setappend (&qh del_vertices, vertex);
	    trace2((qh ferr, "qh_updatevertices: delete vertex p%d (v%d) in f%d\n",
		  qh_pointid(vertex->point), vertex->id, visible->id));
  	  }
        }
      }
    }
  }else {  /* !VERTEXneighbors */
    FORALLvisible_facets {
      FOREACHvertex_(visible->vertices) {
        if (!vertex->newlist && !vertex->deleted) {
          vertex->deleted= True;
	  qh_setappend (&qh del_vertices, vertex);
	  trace2((qh ferr, "qh_updatevertices: delete vertex p%d (v%d) in f%d\n",
		  qh_pointid(vertex->point), vertex->id, visible->id));
  	}
      }
    }
  }
} /* updatevertices */