poly.c

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

    get vertices with nth neighbor by deleting nth vertex
  if PREmerge/MERGEexact or FORCEoutput 
    all facets check for flipped (also prevents point partitioning)
  if duplicate ridges and PREmerge/MERGEexact
    facet->dupridge set
    missing neighbor links identifies extra ridges to be merging
notes:
  do not allocate memory after hash_table (need to free it cleanly)
*/
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= SETelem_(facet->neighbors, k);
	  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 */

    
/*----------------------------------------
-matchvertices- tests whether vertices match with a single skip
  starts match at firstindex since all new facets have a common vertex
  assumes skipA is in A and both sets are the same size
returns:
  skip index
  sets same iff vertices have the same 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 */

/*----------------------------------------
-newfacet- creates and allocates space for a facet
returns:
    all fields initialized or cleared (NULL)
    preallocates neighbors
*/
facetT *qh_newfacet(void) {
  facetT *facet;
  void **freelistp;
  
  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 */


/*----------------------------------------
-newridge- creates and allocates space for a ridge
*/
ridgeT *qh_newridge(void) {
  ridgeT *ridge;
  void **freelistp;

  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 */


/*-------------------------------------------------
-pointid- return id for a point, -3 if null, -2 if interior, or -1 if not known
notes:
  alternative code:
  unsigned long id;
  id= ((unsigned long)point - (unsigned long)qh first_point)/qh normal_size;
*/
int qh_pointid (pointT *point) {
  long offset, id;

  if (!point)
    return -3;
  offset= point - qh first_point;
  id= offset / qh hull_dim;
  if (id >= qh num_points || id < 0) {
    if (point == qh interior_point)
      id= -2;
    else if ((id= qh_setindex (qh other_points, point)) != -1)
      id += qh num_points;
  }
  return (int) id;
} /* pointid */
  
/*-------------------------------------------------
-removefacet- unlinks facet from qh facet_list,
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 */


/*-------------------------------------------------
-removevertex- unlinks vertex from qh vertex_list,
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 */


/*-----------------------------------------------
-updatevertices - update vertex neighbors and delete interior vertices
  if qh VERTEXneighbors, update neighbors for each vertex
  interior vertices added to qh del_vertices for later partitioning
*/
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 */