merge.c

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

/* merge.c - merges non-convex facets

   see README and merge.h
   
   other modules call qh_premerge() and qh_postmerge()
   
   the user may call qh_postmerge() to perform additional merges.

     To remove deleted facets and vertices (qhull() in qhull.c):
      qh_partitionvisible (!qh_ALL, &numoutside);  // visible_list, newfacet_list
      qh_deletevisible ();         // qh visible_list
      qh_resetlists (False);       // qh visible_list newvertex_list newfacet_list 

     and optionally execute:
      qh_check_maxout();

   to avoid loading merge functions, define qh_NOmerge in user.h 

   assumes qh CENTERtype= centrum
   
   merges occur in qh_mergefacet and in qh_mergecycle
   vertex->neighbors not set until the first merge occurs
      
   copyright (c) 1993-1995 The Geometry Center        
*/

#include "qhull_a.h"

#ifndef qh_NOmerge

static int qh_compareangle(const void *p1, const void *p2);
static int qh_comparemerge(const void *p1, const void *p2);
static int qh_comparevisit (const void *p1, const void *p2);

																														
/* ===== functions (alphabetical after premerge and postmerge) ======*/

/*-------------------------------------------------
-premerge- pre-merge nonconvex facets in newfacet_list for apex
  maxcentrum and maxangle define coplanar and concave (qh_test_appendmerge)
  uses globals, MERGEexact, PREmerge
returns:
  deleted facets added to visible_list with facet->visible set
*/
void qh_premerge (vertexT *apex, realT maxcentrum, realT maxangle) {
  boolT othermerge= False;
  facetT *newfacet;
  
  if (qh ZEROcentrum && qh_checkzero(!qh_ALL))
    return;    
  trace2((qh ferr, "qh_premerge: premerge centrum %2.2g angle %2.2g for apex v%d facetlist f%d\n",
	    maxcentrum, maxangle, apex->id, getid_(qh newfacet_list)));
  if (qh IStracing >= 4 && qh num_facets < 50)
    qh_printlists();
  qh centrum_radius= maxcentrum;
  qh cos_max= maxangle;
  qh degen_mergeset= qh_settemp (qh TEMPsize);
  qh facet_mergeset= qh_settemp (qh TEMPsize);
  if (qh hull_dim >=3) { 
    qh_mark_dupridges (qh newfacet_list); /* facet_mergeset */
    qh_mergecycle_all (qh newfacet_list, &othermerge);
    qh_forcedmerges (&othermerge /* qh facet_mergeset */); 
    FORALLnew_facets {  /* test samecycle merges */
      if (!newfacet->simplicial && !newfacet->mergeridge)
	qh_degen_redundant_neighbors (newfacet, NULL);
    }
    if (qh_merge_degenredundant())
      othermerge= True;
  }else /* qh hull_dim == 2 */
    qh_mergecycle_all (qh newfacet_list, &othermerge);
  qh_flippedmerges (qh newfacet_list, &othermerge);
  if (!qh MERGEexact || zzval_(Ztotmerge)) {
    zinc_(Zpremergetot);
    qh POSTmerging= False;
    qh_getmergeset_initial (qh newfacet_list);
    qh_all_merges (othermerge, False);
  }
  qh_settempfree(&qh facet_mergeset);
  qh_settempfree(&qh degen_mergeset);
} /* premerge */
  
/*-------------------------------------------------
-postmerge- post-merge nonconvex facets as defined by maxcentrum/maxangle
  if firsttime (visible list not set), 
    builds facet_newlist, newvertex_list
  calls 
  'reason' is for reporting progress
  if firstmerge, calls qh_reducevertices before  getmergeset
  if vneighbors, calls qh_test_vneighbors at end of qh_all_merge 
returns:
  deleted facets added to visible_list with facet->visible
  visible_list == facet_list
*/
void qh_postmerge (char *reason, realT maxcentrum, realT maxangle, 
                      boolT vneighbors) {
  facetT *newfacet;
  boolT othermerges= False;
  vertexT *vertex;

  if (qh REPORTfreq || qh IStracing) {
    qh_buildtracing (NULL, NULL);
    qh_printsummary (qh ferr);
    if (qh PRINTstatistics) 
      qh_printallstatistics (qh ferr, "reason");
    fprintf (qh ferr, "\n%s with 'C%.2g' and 'A%.2g'\n", 
        reason, maxcentrum, maxangle);
  }
  trace2((qh ferr, "qh_postmerge: postmerge.  test vneighbors? %d\n",
	    vneighbors));
  qh centrum_radius= maxcentrum;
  qh cos_max= maxangle;
  qh POSTmerging= True;
  qh degen_mergeset= qh_settemp (qh TEMPsize);
  qh facet_mergeset= qh_settemp (qh TEMPsize);
  if (qh visible_list != qh facet_list) {  /* first call */
    qh NEWfacets= True;
    qh visible_list= qh newfacet_list= qh facet_list;
    FORALLnew_facets {
      newfacet->newfacet= True;
       if (!newfacet->simplicial)
        newfacet->newmerge= True;
     zinc_(Zpostfacets);
    }
    qh newvertex_list= qh vertex_list;
    FORALLvertices
      vertex->newlist= True;
    if (qh VERTEXneighbors) { /* a merge has occurred */
      FORALLvertices
	vertex->delridge= True; /* test for redundant, needed? */
      if (qh MERGEexact) {
	if (qh hull_dim <= qh_DIMreduceBuild)
	  qh_reducevertices(); /* was skipped during pre-merging */
      }
    }
    if (!qh PREmerge && !qh MERGEexact) 
      qh_flippedmerges (qh newfacet_list, &othermerges);
  }
  qh_getmergeset_initial (qh newfacet_list);
  qh_all_merges (False, vneighbors);
  qh_settempfree(&qh facet_mergeset);
  qh_settempfree(&qh degen_mergeset);
} /* post_merge */

/*-------------------------------------------------
-all_merges- merge all non-convex facets
  othermerge set if already merged (for qh_reducevertices)
  if vneighbors, tests vertex neighbors for convexity at end
  qh facet_mergeset lists the non-convex ridges in qh_newfacet_list
  qh degen_mergeset is defined
  if MERGEexact && !POSTmerging, does not merge coplanar facets
returns:
  deleted facets added to visible_list with facet->visible
  deleted vertices added delvertex_list with vertex->delvertex
notes:
  unless !MERGEindependent, merges facets in independent sets
  uses qh newfacet_list as argument since merges call removefacet()
*/
void qh_all_merges (boolT othermerge, boolT vneighbors) {
  facetT *facet1, *facet2;
  mergeT *merge;
  boolT wasmerge= True, isreduce;
  void **freelistp;
  vertexT *vertex;
  mergeType mergetype;
  int numcoplanar=0, numconcave=0, numdegenredun= 0, numnewmerges= 0;
  
  trace2((qh ferr, "qh_merge_all: starting to merge facets beginning from f%d\n",
	    getid_(qh newfacet_list)));
  while (True) {
    wasmerge= False;
    while (qh_setsize (qh facet_mergeset)) {
      while ((merge= (mergeT*)qh_setdellast(qh facet_mergeset))) {
	facet1= merge->facet1;
	facet2= merge->facet2;
	mergetype= merge->type;
	qh_memfree_(merge, sizeof(mergeT), freelistp);
	if (facet1->visible || facet2->visible) /*deleted facet*/
	  continue;  
	if ((facet1->newfacet && !facet1->tested)
	        || (facet2->newfacet && !facet2->tested)) {
	  if (qh MERGEindependent && mergetype <= MRGanglecoplanar)
	    continue;      /* perform independent sets of merges */
	}
	qh_merge_nonconvex (facet1, facet2, mergetype);
        numdegenredun += qh_merge_degenredundant();
        numnewmerges++;
        wasmerge= True;
	if (mergetype == MRGconcave)
	  numconcave++;
	else /* MRGcoplanar or MRGanglecoplanar */
	  numcoplanar++;
      } /* while setdellast */
      if (qh POSTmerging && qh hull_dim <= qh_DIMreduceBuild 
      && numnewmerges > qh_MAXnewmerges) {
	numnewmerges= 0;
	qh_reducevertices();  /* otherwise large post merges too slow */
      }
      qh_getmergeset (qh newfacet_list); /* facet_mergeset */
    } /* while mergeset */
    if (qh VERTEXneighbors) {
      isreduce= False;
      if (qh hull_dim >=4 && qh POSTmerging) {
	FORALLvertices  
	  vertex->delridge= True;
	isreduce= True;
      }
      if ((wasmerge || othermerge) && (!qh MERGEexact || qh POSTmerging) 
	  && qh hull_dim <= qh_DIMreduceBuild) {
	othermerge= False;
	isreduce= True;
      }
      if (isreduce) {
	if (qh_reducevertices()) {
	  qh_getmergeset (qh newfacet_list); /* facet_mergeset */
	  continue;
	}
      }
    }
    if (vneighbors && qh_test_vneighbors(/* qh newfacet_list */)) 
      continue;
    break;
  } /* while (True) */
  if (qh CHECKfrequently && !qh MERGEexact)
    qh_checkconvex (qh newfacet_list, qh_ALGORITHMfault);
  /* qh_checkconnect (); [this is slow and it changes the facet order] */
  trace1((qh ferr, "qh_merge_all: merged %d coplanar facets %d concave facets and %d degen or redundant facets.\n",
    numcoplanar, numconcave, numdegenredun));
  if (qh IStracing >= 4 && qh num_facets < 50)
    qh_printlists ();
} /* merge_all */


/*-------------------------------------------------
-appendmergeset- appends an entry to facet_mergeset
  angle ignored if NULL or !qh ANGLEmerge
returns:
  merge appended to facet_mergeset or degen_mergeset
    sets ->degenerate or ->redundant if degen_mergeset
notes:
  see test_appendmerge()
*/
void qh_appendmergeset(facetT *facet, facetT *neighbor, mergeType mergetype, realT *angle) {
  mergeT *merge, *lastmerge;
  void **freelistp;

  if (facet->redundant)
    return;
  if (facet->degenerate && mergetype == MRGdegen)
    return;
  qh_memalloc_(sizeof(mergeT), freelistp, merge, mergeT);
  merge->facet1= facet;
  merge->facet2= neighbor;
  merge->type= mergetype;
  if (angle && qh ANGLEmerge)
    merge->angle= *angle;
  if (mergetype < MRGdegen)
    qh_setappend (&(qh facet_mergeset), merge);
  else if (mergetype == MRGdegen) {
    facet->degenerate= True;
    if (!(lastmerge= (mergeT*)qh_setlast (qh degen_mergeset)) 
    || lastmerge->type == MRGdegen)
      qh_setappend (&(qh degen_mergeset), merge);
    else
      qh_setaddnth (&(qh degen_mergeset), 0, merge);
  }else { /* mergetype == MRGredundant */
    facet->redundant= True;
    qh_setappend (&(qh degen_mergeset), merge);
  }
} /* appendmergeset */


/*-----------------------------------------------
-basevertices- return temporary set of base vertices for samecycle
  assumes apex is SETfirst
returns:
  vertices (settemp)
  all ->seen are cleared
notes:
  uses qh_vertex_visit;
*/
setT *qh_basevertices (facetT *samecycle) {
  facetT *same;
  vertexT *apex, *vertex, **vertexp;
  setT *vertices= qh_settemp (qh TEMPsize);
  
  apex= SETfirst_(samecycle->vertices);
  apex->visitid= ++qh vertex_visit;
  FORALLsame_cycle_(samecycle) {
    if (same->mergeridge)
      continue;
    FOREACHvertex_(same->vertices) {
      if (vertex->visitid != qh vertex_visit) {
        qh_setappend (&vertices, vertex);
        vertex->visitid= qh vertex_visit;
        vertex->seen= False;
      }
    }
  }
  trace4((qh ferr, "qh_basevertices: found %d vertices\n", 
         qh_setsize (vertices)));
  return vertices;
} /* basevertices */

/*-------------------------------------------------
-checkconnect- check that new facets are connected
notes:
  this is slow and it changes the order of the facets
  uses qh visit_id
*/
void qh_checkconnect (void /* qh newfacet_list */) {
  facetT *facet, *newfacet, *errfacet= NULL, *neighbor, **neighborp;

  facet= qh newfacet_list;
  qh_removefacet (facet);
  qh_appendfacet (facet);
  facet->visitid= ++qh visit_id;
  FORALLfacet_(facet) {
    FOREACHneighbor_(facet) {
      if (neighbor->visitid != qh visit_id) {
        qh_removefacet (neighbor);
        qh_appendfacet (neighbor);
        neighbor->visitid= qh visit_id;
      }
    }
  }
  FORALLnew_facets {
    if (newfacet->visitid == qh visit_id)
      break;
    fprintf(qh ferr, "qhull error: f%d is not attached to the new facets\n",
         newfacet->id);
    errfacet= newfacet;
  }
} /* checkconnect */

/*-------------------------------------------------
-checkzero- check that facets are clearly convex for DISTround
  if testall, test all facets for qh MERGEexact post-merging
  else test qh newfacet_list
    if qh MERGEexact, allows coplanar ridges
        skip convexity test while qh ZEROall_ok
returns:
  returns True if all facets !flipped, !dupridge, normal
      if all horizon facets are simplicial
      if all vertices are clearly below neighbor
      if all opposite vertices of horizon are below 
  clears qh ZEROall_ok if any problems or coplanar facets
notes:
  uses qh vertex_visit, horizon facets may define multiple new facets
*/
boolT qh_checkzero (boolT testall) {
  facetT *facet, *neighbor, **neighborp;
  facetT *horizon, *facetlist;
  int neighbor_i;
  vertexT *vertex, **vertexp;
  realT dist;

  if (testall) 
    facetlist= qh facet_list;
  else {
    facetlist= qh newfacet_list;
    FORALLfacet_(facetlist) {
      horizon= SETfirst_(facet->neighbors);
      if (!horizon->simplicial)
        goto LABELproblem;
      if (facet->flipped || facet->dupridge || !facet->normal)
        goto LABELproblem;
    }
    if (qh MERGEexact && qh ZEROall_ok) {
      trace2((qh ferr, "qh_checkzero: skip convexity check until first pre-merge\n"));
      return True;
    }
  }
  FORALLfacet_(facetlist) {
    qh vertex_visit++;
    neighbor_i= 0;
    horizon= NULL;
    FOREACHneighbor_(facet) {
      if (!neighbor_i && !testall) {
        horizon= neighbor;
	neighbor_i++;
        continue; /* horizon facet tested in qh_findhorizon */