poly2.c

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

	totvneighbors += qh_setsize (vertex->neighbors);
      }
      FORALLfacet_(facetlist)
	totvertices += qh_setsize (facet->vertices);
      if (totvneighbors != totvertices) {
	fprintf(qh ferr, "qhull internal error (qh_checkpolygon): vertex neighbors inconsistent.  Totvneighbors %d, totvertices %d\n",
		totvneighbors, totvertices);
	waserror= True;
      }
    }
    if (numvertices != qh num_vertices - qh_setsize(qh del_vertices)) {
      fprintf(qh ferr, "qhull internal error (qh_checkpolygon): actual number of vertices is %d, cumulative vertex count is %d\n",
	      numvertices, qh num_vertices - qh_setsize(qh del_vertices));
      waserror= True;
    }
    if (qh hull_dim == 2 && numvertices != numfacets) {
      fprintf (qh ferr, "qhull internal error (qh_checkpolygon): #vertices %d != #facets %d\n",
        numvertices, numfacets);
      waserror= True;
    }
    if (qh hull_dim == 3 && numvertices + numfacets - numridges/2 != 2) {
      fprintf (qh ferr, "qhull internal error (qh_checkpolygon): #vertices %d + #facets %d - #edges %d != 2\n",
        numvertices, numfacets, numridges/2);
      waserror= True;
    }
  }
  if (waserror) 
    qh_errexit(qh_ERRqhull, NULL, NULL);
} /* checkpolygon */


/*-------------------------------------------------
-checkvertex- check vertex for consistency
notes:
  neighbors checked efficiently in checkpolygon
*/
void qh_checkvertex (vertexT *vertex) {
  boolT waserror= False;
  facetT *neighbor, **neighborp, *errfacet=NULL;
  int size;

  if (qh_pointid (vertex->point) == -1) {
    fprintf (qh ferr, "qhull internal error (qh_checkvertex): unknown point id %p\n", vertex->point);
    waserror= True;
  }
  if (vertex->id >= qh vertex_id) {
    fprintf (qh ferr, "qhull internal error (qh_checkvertex): unknown vertex id %d\n", vertex->id);
    waserror= True;
  }
  if (!waserror && !vertex->deleted) {
    if ((size= qh_setsize (vertex->neighbors))) {
      FOREACHneighbor_(vertex) {
        if (!qh_setin (neighbor->vertices, vertex)) {
          fprintf (qh ferr, "qhull internal error (qh_checkvertex): neighbor f%d does not contain v%d\n", neighbor->id, vertex->id);
	  errfacet= neighbor;
	  waserror= True;
	}
      }
    }
  }
  if (waserror) {
    qh_errprint ("ERRONEOUS", NULL, NULL, NULL, vertex);
    qh_errexit (qh_ERRqhull, errfacet, NULL);
  }
} /* checkvertex */
  
/*-------------------------------------------------
-clearcenters- clear old data from facet->center
  sets new centertype
  nop if CENTERtype is the same
*/
void qh_clearcenters (qh_CENTER type) {
  facetT *facet;
  
  if (qh CENTERtype != type) {
    FORALLfacets {
      if (qh CENTERtype == qh_ASvoronoi){
        if (facet->center) {
          qh_memfree (facet->center, qh center_size);
          facet->center= NULL;
        }
      }else /* qh CENTERtype == qh_AScentrum */ {
        if (facet->center) {
          qh_memfree (facet->center, qh normal_size);
	  facet->center= NULL;
        }
      }
    }
    qh CENTERtype= type;
  }
  trace2((qh ferr, "clearcenters: switched to center type %d\n", type));
} /* clearcenters */

/*----------------------------------------
-createsimplex- creates a simplex from a set of vertices
returns:
    initializes qh facet_list to the simplex
*/
void qh_createsimplex(setT *vertices) {
  facetT *facet= NULL, *newfacet;
  boolT toporient= True;
  int vertex_i, vertex_n, nth;
  setT *newfacets= qh_settemp (qh hull_dim+1);
  vertexT *vertex;
  
  qh facet_list= qh newfacet_list= qh facet_tail= qh_newfacet();
  qh vertex_list= qh newvertex_list= qh vertex_tail= qh_newvertex(NULL);
  FOREACHvertex_i_(vertices) {
    newfacet= qh_newfacet();
    newfacet->vertices= qh_setnew_delnthsorted (vertices, vertex_n,
						vertex_i, 0);
    newfacet->toporient= toporient;
    qh_appendfacet(newfacet);
    newfacet->newfacet= True;
    qh_appendvertex (vertex);
    qh_setappend (&newfacets, newfacet);
    toporient ^= True;
  }
  FORALLnew_facets {
    nth= 0;
    FORALLfacet_(qh newfacet_list) {
      if (facet != newfacet) 
        SETelem_(newfacet->neighbors, nth++)= facet;
    }
    qh_settruncate (newfacet->neighbors, qh hull_dim);
  }
  qh_settempfree (&newfacets);
  trace1((qh ferr, "qh_createsimplex: created simplex\n"));
} /* createsimplex */

/*-------------------------------------------------
-delridge- deletes ridge from data structures it belongs to and frees up the
memory occupied by it
notes:
  in merge.c, caller sets vertex->delridge for each vertex
  also freed in qh_freeqhull
*/
void qh_delridge(ridgeT *ridge) {
  void **freelistp;
  
  qh_setdel(ridge->top->ridges, ridge);
  qh_setdel(ridge->bottom->ridges, ridge);
  qh_setfree(&(ridge->vertices));
  qh_memfree_(ridge, sizeof(ridgeT), freelistp);
} /* delridge */


/*-------------------------------------------------
-delvertex- deletes a vertex and frees its memory
  assumes vertex->adjacencies have been updated if needed
  unlinks for vertex_list
*/
void qh_delvertex (vertexT *vertex) {

  if (vertex == qh tracevertex)
    qh tracevertex= NULL;
  qh_removevertex (vertex);
  qh_setfree (&vertex->neighbors);
  qh_memfree(vertex, sizeof(vertexT));
} /* delvertex */


/*----------------------------------------
-facet3vertex- return temporary set of 3-d vertices
  in qh_ORIENTclock order
*/
setT *qh_facet3vertex (facetT *facet) {
  ridgeT *ridge, *firstridge;
  vertexT *vertex;
  int cntvertices, cntprojected=0;
  setT *vertices;

  cntvertices= qh_setsize(facet->vertices);
  vertices= qh_settemp (cntvertices);
  if (facet->simplicial) {
    if (cntvertices != 3) {
      fprintf (qh ferr, "qhull internal error (qh_facet3vertex): only %d vertices for simplicial facet f%d\n", 
                  cntvertices, facet->id);
      qh_errexit(qh_ERRqhull, facet, NULL);
    }
    qh_setappend (&vertices, SETfirst_(facet->vertices));
    if (facet->toporient ^ qh_ORIENTclock)
      qh_setappend (&vertices, SETsecond_(facet->vertices));
    else
      qh_setaddnth (&vertices, 0, SETsecond_(facet->vertices));
    qh_setappend (&vertices, SETelem_(facet->vertices, 2));
  }else {
    ridge= firstridge= SETfirst_(facet->ridges);   /* no infinite */
    while ((ridge= qh_nextridge3d (ridge, facet, &vertex))) {
      qh_setappend (&vertices, vertex);
      if (++cntprojected > cntvertices || ridge == firstridge)
        break;
    }
    if (!ridge || cntprojected != cntvertices) {
      fprintf (qh ferr, "qhull internal error (qh_facet3vertex): ridges for facet %d don't match up.  got at least %d\n", 
                  facet->id, cntprojected);
      qh_errexit(qh_ERRqhull, facet, ridge);
    }
  }
  return vertices;
} /* facet3vertex */


/*-------------------------------------------------
-findfacet- find facet that is furthest below a point 
  if 'facet' defined, starts search at 'facet'
  if NULL 'facet', starts at first non-flipped non-upperDelaunay facet
  searches neighbors of coplanar and most flipped facets
  does not search upper envelope of Delaunay triangulations
returns:
  best facet
  distance to facet
  isoutside if point is outside of the hull
  numpart counts the number of distance tests
notes:
  this works for Delaunay triangulations.
  this works if the initial facet is not above the point.
  warning: a lens-shaped hull may mistakenly return a facet
      that is above the point
  uses qh visit_id, qh searchset
*/
facetT *qh_findfacet (pointT *point, facetT *facet, 
           realT *dist, boolT *isoutside, int *numpart) {
  facetT *bestfacet;
  
  if (!facet)
    facet= qh_nonupper( qh facet_list);
  bestfacet= qh_findbest (point, facet, qh_ALL, False, dist, isoutside, numpart);
  return bestfacet;
} /* findfacet */ 
 
/*-------------------------------------------------
-findgood- identify good facets for qh ONLYgood
  GOODvertex>0 - facet includes point as vertex
    if !match, returns goodhorizon
    inactive if qh MERGING
  GOODpoint - facet is visible or coplanar (>0) or not visible (<0) 
  GOODthreshold - facet->normal matches threshold
    if !goodhorizon and !match, selects facet with closest angle
       and sets GOODclosest
returns:
  number of new, good facets found
  determins facet->good
  may update GOODclosest
notes:
  findgood_all further reduces the good region
*/
int qh_findgood (facetT *facetlist, int goodhorizon) {
  facetT *facet, *bestfacet;
  realT angle, bestangle, dist;
  int  numgood=0;

  if (qh GOODclosest) {
    bestfacet= qh GOODclosest;
    qh_inthresholds (bestfacet->normal, &bestangle);
  }else {
    bestfacet= NULL;
    bestangle= -REALmax;
  }
  FORALLfacet_(facetlist) {
    if (facet->good)
      numgood++;
  }
  if (qh GOODvertex>0 && !qh MERGING) {
    FORALLfacet_(facetlist) {
      if (!qh_isvertex (qh GOODvertexp, facet->vertices)) {
        facet->good= False;
        numgood--;
      }
    }
  }
  if (qh GOODpoint && numgood) {
    FORALLfacet_(facetlist) {
      if (facet->good && facet->normal) {
        zinc_(Zdistgood);
        qh_distplane (qh GOODpointp, facet, &dist);
        if ((qh GOODpoint > 0) ^ (dist > 0.0)) {
          facet->good= False;
          numgood--;
        }
      }
    }
  }
  if (qh GOODthreshold && (numgood || goodhorizon)) {
    FORALLfacet_(facetlist) {
      if (facet->good && facet->normal) {
        if (!qh_inthresholds (facet->normal, &angle)) {
          facet->good= False;
          numgood--;
          angle= fabs_(angle);
          if (angle > bestangle) {
            bestangle= angle;
            bestfacet= facet;
          }
        }
      }
    }
    if (!numgood && bestfacet && bestfacet != qh GOODclosest) {
      if (qh GOODclosest)
	qh GOODclosest->good= False;
      qh GOODclosest= bestfacet;
      bestfacet->good= True;
      numgood++;
      trace2((qh ferr, "qh_findgood: f%d is closest (%2.2g) to thresholds\n", 
           bestfacet->id, bestangle));
      return numgood;
    }else if (numgood && qh GOODclosest)
      qh GOODclosest->good= False;
  }
  zadd_(Zgoodfacet, numgood);
  trace2((qh ferr, "qh_findgood: found %d good facets\n", numgood));
  if (!numgood && qh GOODvertex>0 && !qh MERGING) 
    return goodhorizon;
  return numgood;
} /* findgood */

/*-------------------------------------------------
-findgood_all- apply other constraints for good facets (used by qh PRINTgood)
  GOODvertex - facet includes (>0) or doesn't include (<0) point as vertex
    if last good facet, prints warning and continues
  SPLITthresholds- facet->normal matches threshold, or if none, the closest one
  calls findgood if !ONLYgood
  nop if good not used
returns:
  clears facet->good if not good
  sets qh num_good
notes:
  this is like findgood but more restrictive
*/
void qh_findgood_all (facetT *facetlist) {
  facetT *facet, *bestfacet=NULL;
  realT angle, bestangle= REALmax;
  int  numgood=0, startgood;

  if (!qh GOODvertex && !qh GOODthreshold && !qh GOODpoint 
  && !qh SPLITthresholds)
    return;
  if (!qh ONLYgood)
    qh_findgood (qh facet_list, 0);
  FORALLfacet_(facetlist) {
    if (facet->good)
      numgood++;
  }
  if (qh GOODvertex <0 || (qh GOODvertex > 0 && qh MERGING)) {
    FORALLfacet_(facetlist) {
      if (facet->good && ((qh GOODvertex > 0) ^ !!qh_isvertex (qh GOODvertexp, facet->vertices))) {
        if (!--numgood) {
          fprintf (qh ferr, "qhull warning: good vertex p%d does not match last good facet f%d.  Ignored.\n",
             qh_pointid(qh GOODvertexp), facet->id);
          return;
        }
        facet->good= False;
      }
    }
  }
  startgood= numgood;
  if (qh SPLITthresholds) {
    FORALLfacet_(facetlist) {
      if (facet->good) {
        if (!qh_inthresholds (facet->normal, &angle)) {
          facet->good= False;
          numgood--;
          angle= fabs_(angle);
          if (angle < bestangle) {
            bestangle= angle;
            bestfacet= facet;
          }
        }
      }
    }
    if (!numgood) {
      bestfacet->good= True;
      numgood++;
      trace0((qh ferr, "qh_findgood_all: f%d is closest (%2.2g) to thresholds\n", 
           bestfacet->id, bestangle));