poly.c

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

/* poly.c -- implements polygons and simplices

   see README, poly.h and qhull.h
   
   copyright (c) 1993-1995, The Geometry Center

   infrequent code is in poly2.c 
        (all but top 50 and their callers 12/3/95)
*/

#include "qhull_a.h"

/*======== functions in alphabetical order ==========*/

/*-------------------------------------------------
-appendfacet- appends facet to end of qh facet_list,
  updates qh facet_list, facet_tail, newfacet_list, facet_next
  increments qh numfacets
  assumes qh facet_list/facet_tail is defined (createsimplex)
*/
void qh_appendfacet(facetT *facet) {
  facetT *tail= qh facet_tail;

  if (tail == qh newfacet_list)
    qh newfacet_list= facet;
  if (tail == qh facet_next)
    qh facet_next= facet;
  facet->previous= tail->previous;
  facet->next= tail;
  if (tail->previous)
    tail->previous->next= facet;
  else
    qh facet_list= facet;
  tail->previous= facet;
  qh num_facets++;
  trace4((qh ferr, "qh_appendfacet: append f%d to facet_list\n", facet->id));
} /* appendfacet */


/*-------------------------------------------------
-appendvertex- appends vertex to end of qh vertex_list,
  updates qh vertex_list, vertex_tail, newvertex_list
  increments qh num_vertices
  sets vertex->newlist
  assumes qh vertex_list/vertex_tail is defined (createsimplex)
*/
void qh_appendvertex (vertexT *vertex) {
  vertexT *tail= qh vertex_tail;

  if (tail == qh newvertex_list)
    qh newvertex_list= vertex;
  vertex->newlist= True;
  vertex->previous= tail->previous;
  vertex->next= tail;
  if (tail->previous)
    tail->previous->next= vertex;
  else
    qh vertex_list= vertex;
  tail->previous= vertex;
  qh num_vertices++;
  trace4((qh ferr, "qh_appendvertex: append v%d to vertex_list\n", vertex->id));
} /* appendvertex */


/*-------------------------------------------------
-attachnewfacets- attach horizon facets to new facets in qh newfacet_list
  only needed for qh ONLYgood
  newfacets have neighbor and ridge links to horizon but not vice versa
  will set NEWfacets
returns:
  horizon facets linked to new facets 
     ridges changed from visible facets to new facets
     simplicial ridges deleted
  qh visible_list, no ridges valid
     ->f.replace is a newfacet (if any)
*/
void qh_attachnewfacets (void ) {
  facetT *newfacet= NULL, *neighbor, **neighborp, *horizon, *visible;
  ridgeT *ridge, **ridgep;

  qh NEWfacets= True;
  trace3((qh ferr, "qh_attachnewfacets: delete interior ridges\n"));
  qh visit_id++;
  FORALLvisible_facets {
    visible->visitid= qh visit_id;
    if (visible->ridges) {
      FOREACHridge_(visible->ridges) {
	neighbor= otherfacet_(ridge, visible);
	if (neighbor->visitid == qh visit_id
	    || (!neighbor->visible && neighbor->simplicial)) {
	  if (!neighbor->visible)  /* delete ridge for simplicial horizon */
	    qh_setdel (neighbor->ridges, ridge);
	  qh_setfree (&(ridge->vertices)); /* delete on 2nd visit */
	  qh_memfree (ridge, sizeof(ridgeT));
	}
      }
      SETfirst_(visible->ridges)= NULL;
    }
    SETfirst_(visible->neighbors)= NULL;
  }
  trace1((qh ferr, "qh_attachnewfacets: attach horizon facets to new facets\n"));
  FORALLnew_facets {
    horizon= SETfirst_(newfacet->neighbors);
    if (horizon->simplicial) {
      visible= NULL;
      FOREACHneighbor_(horizon) {   /* may have more than one horizon ridge */
	if (neighbor->visible) {
	  if (visible) {
	    if (qh_setequal_skip (newfacet->vertices, 0, horizon->vertices,
				  SETindex_(horizon->neighbors, neighbor))) {
	      visible= neighbor;
	      break;
	    }
	  }else
	    visible= neighbor;
	}
      }
      if (visible) {
	visible->f.replace= newfacet;
	qh_setreplace (horizon->neighbors, visible, newfacet);
      }else {
	fprintf (qh ferr, "qhull internal error (qh_attachnewfacets): couldn't find visible facet for horizon f%d of newfacet f%d\n",
		 horizon->id, newfacet->id);
	qh_errexit2 (qh_ERRqhull, horizon, newfacet);
      }
    }else { /* non-simplicial, with a ridge for newfacet */
      FOREACHneighbor_(horizon) {    /* may hold for many new facets */
	if (neighbor->visible) {
	  neighbor->f.replace= newfacet;
	  qh_setdelnth (horizon->neighbors,
			SETindex_(horizon->neighbors, neighbor));
	  neighborp--; /* repeat */
	}
      }
      qh_setappend (&horizon->neighbors, newfacet);
      ridge= SETfirst_(newfacet->ridges);
      if (ridge->top == horizon)
	ridge->bottom= newfacet;
      else
	ridge->top= newfacet;
      }
  } /* newfacets */
  if (qh PRINTstatistics) {
    FORALLvisible_facets {
      if (!visible->f.replace) 
	zinc_(Zinsidevisible);
    }
  }
} /* attachnewfacets */

/*-------------------------------------------------
-checkflipped- checks facet orientation to interior point
  tests against 0 if !allerror since tested against DISTround before
returns:
  False if flipped orientation (sets facet->flipped)
  distance if non-NULL
*/
boolT qh_checkflipped (facetT *facet, realT *distp, boolT allerror) {
  realT dist;

  if (facet->flipped && !distp)
    return False;
  zzinc_(Zdistcheck);
  qh_distplane(qh interior_point, facet, &dist);
  if (distp)
    *distp= dist;
  if ((allerror && dist > -qh DISTround)|| (!allerror && dist >= 0.0)) {
    facet->flipped= True;
    zzinc_(Zflippedfacets);
    trace0((qh ferr, "qh_checkflipped: facet f%d is flipped, distance= %6.12g during p%d\n",
              facet->id, dist, qh furthest_id));
    return False;
  }
  return True;
} /* checkflipped */

/*-------------------------------------------------
-delfacet- removes facet from facet_list and frees up its memory
   assumes vertices and ridges already freed
*/
void qh_delfacet(facetT *facet) {
  void **freelistp;

  trace5((qh ferr, "qh_delfacet: delete f%d\n", facet->id));
  if (facet == qh tracefacet)
    qh tracefacet= NULL;
  qh_removefacet(facet);
  qh_memfree_(facet->normal, qh normal_size, freelistp);
  if (qh CENTERtype == qh_ASvoronoi) {   /* uses macro calls */
    qh_memfree_(facet->center, qh center_size, freelistp);
  }else /* AScentrum */ {
    qh_memfree_(facet->center, qh normal_size, freelistp);
  }
  qh_setfree(&(facet->neighbors));
  if (facet->ridges)
    qh_setfree(&(facet->ridges));
  qh_setfree(&(facet->vertices));
  if (facet->outsideset)
    qh_setfree(&(facet->outsideset));
  if (facet->coplanarset)
    qh_setfree(&(facet->coplanarset));
  qh_memfree_(facet, sizeof(facetT), freelistp);
} /* delfacet */


/*-------------------------------------------------
-deletevisible- delete visible facets and vertices
    ridges already deleted
    horizon facets do not reference facets on qh visible_list
    new facets in qh newfacet_list
returns:
    deletes each facet and removes from facetlist
    uses qh visit_id;
    qh visible_list empty (== qh newfacet_list)
*/
void qh_deletevisible (/*qh visible_list*/) {
  facetT *visible, *nextfacet;
  vertexT *vertex, **vertexp;
  int numvisible= 0, numdel= qh_setsize(qh del_vertices);

  trace1((qh ferr, "qh_deletevisible: delete %d visible facets and %d vertices\n",
         qh num_visible, numdel));
  for (visible= qh visible_list; visible && visible->visible; 
                visible= nextfacet) { /* deleting current */
    nextfacet= visible->next;        
    numvisible++;
    qh_delfacet(visible);
  }
  if (numvisible != qh num_visible) {
    fprintf (qh ferr, "qhull internal error (qh_deletevisible): qh num_visible %d is not number of visible facets %d\n",
             qh num_visible, numvisible);
    qh_errexit (qh_ERRqhull, NULL, NULL);
  }
  qh num_visible= 0;
  zadd_(Zvisfacettot, numvisible);
  zmax_(Zvisfacetmax, numvisible);
  zadd_(Zdelvertextot, numdel);
  zmax_(Zdelvertexmax, numdel);
  FOREACHvertex_(qh del_vertices) 
    qh_delvertex (vertex);
  qh_settruncate (qh del_vertices, 0);
} /* deletevisible */

/*-------------------------------------------------
-facetintersect- return vertices for intersection of two simplicial facets
  may include 1 prepended entry (if more, need to settemppush)
returns:
  returns set of hull_dim-1 + optional extra
  returns skipped index for each test and checks for exactly one
notes:
  does not need settemp since set in quick memory
  see also qh_vertexintersect and qh_vertexintersect_new
  use qh_setnew_delnthsorted to get nth ridge (no skip information)
*/
setT *qh_facetintersect (facetT *facetA, facetT *facetB,
			 int *skipA,int *skipB, int prepend) {
  setT *intersect;
  int dim= qh hull_dim, i, j;
  facetT **neighborsA, **neighborsB;

  neighborsA= SETaddr_(facetA->neighbors, facetT);
  neighborsB= SETaddr_(facetB->neighbors, facetT);
  i= j= 0;
  if (facetB == *neighborsA++)
    *skipA= 0;
  else if (facetB == *neighborsA++)
    *skipA= 1;
  else if (facetB == *neighborsA++)
    *skipA= 2;
  else {
    for (i= 3; i < dim; i++) {
      if (facetB == *neighborsA++) {
        *skipA= i;
        break;
      }
    }
  }
  if (facetA == *neighborsB++)
    *skipB= 0;
  else if (facetA == *neighborsB++)
    *skipB= 1;
  else if (facetA == *neighborsB++)
    *skipB= 2;
  else {
    for (j= 3; j < dim; j++) {
      if (facetA == *neighborsB++) {
        *skipB= j;
        break;
      }
    }
  }
  if (i >= dim || j >= dim) {
    fprintf (qh ferr, "qhull internal error (qh_facetintersect): f%d or f%d not in others neighbors\n",
            facetA->id, facetB->id);
    qh_errexit2 (qh_ERRqhull, facetA, facetB);
  }
  intersect= qh_setnew_delnthsorted (facetA->vertices, qh hull_dim, *skipA, prepend);
  trace4((qh ferr, "qh_facetintersect: f%d skip %d matches f%d skip %d\n",
	  facetA->id, *skipA, facetB->id, *skipB));
  return(intersect);
} /* facetintersect */

/*----------------------------------------
-gethash- return hashvalue for a set with firstindex and skipelem
  assumes at least firstindex+1 elements
  sum of elements does badly in high d
  assumes skipelem is NULL, in set, or part of hash
*/
unsigned qh_gethash (int hashsize, setT *set, int size, int firstindex, void *skipelem) {
  void **elemp= SETelemaddr_(set, firstindex, void);
  ptr_intT hash, elem;
  int i;