merge.c

关于网格剖分的

#ifndef qh_NOtrace
  if (qh IStracing >= 4) {
    fprintf (qh ferr, "qh_find_newvertex: find new vertex for v%d from ",
	     oldvertex->id);
    FOREACHvertex_(vertices) 
      fprintf (qh ferr, "v%d ", vertex->id);
    FOREACHridge_(ridges)
      fprintf (qh ferr, "r%d ", ridge->id);
    fprintf (qh ferr, "\n");
  }
#endif
  FOREACHvertex_(vertices) 
    vertex->visitid= 0;
  FOREACHridge_(ridges) {
    FOREACHvertex_(ridge->vertices) 
      vertex->visitid++;
  }
  FOREACHvertex_(vertices) {
    if (!vertex->visitid) {
      qh_setdelnth (vertices, SETindex_(vertices,vertex));
      vertexp--; /* repeat since deleted this vertex */
    }
  }
  qh vertex_visit += qh_setsize (ridges);
  if (!qh_setsize (vertices)) {
    trace4((qh ferr, "qh_find_newvertex: vertices not in ridges for v%d\n",
	    oldvertex->id));
    return NULL;
  }
  qsort (SETaddr_(vertices, vertexT), qh_setsize (vertices),
	        sizeof (vertexT *), qh_comparevisit);
  /* can now use qh vertex_visit */
  if (qh PRINTstatistics) {
    size= qh_setsize (vertices);
    zinc_(Zintersect);
    zadd_(Zintersecttot, size);
    zmax_(Zintersectmax, size);
  }
  hashsize= qh_newhashtable (qh_setsize (ridges));
  FOREACHridge_(ridges)
    qh_hashridge (qh hash_table, hashsize, ridge, oldvertex);
  FOREACHvertex_(vertices) {
    newridges= qh_vertexridges (vertex);
    FOREACHridge_(newridges) {
      if (qh_hashridge_find (qh hash_table, hashsize, ridge, vertex, oldvertex, &hash)) {
	zinc_(Zdupridge);
	break;
      }
    }
    qh_settempfree (&newridges);
    if (!ridge)
      break;  /* found a rename */
  }
  if (vertex) {
    /* counted in qh_renamevertex */
    trace2((qh ferr, "qh_find_newvertex: found v%d for old v%d from %d vertices and %d ridges.\n",
      vertex->id, oldvertex->id, qh_setsize (vertices), qh_setsize (ridges)));
  }else {
    zinc_(Zfindfail);
    trace0((qh ferr, "qh_find_newvertex: no vertex for renaming v%d (all duplicated ridges) during p%d\n",
      oldvertex->id, qh furthest_id));
  }
  qh_setfree (&qh hash_table);
  return vertex;
} /* find_newvertex */

/*-<a                             href="qh-c.htm#merge"
  >-------------------------------</a><a name="findbest_test">-</a>
  
  qh_findbest_test( testcentrum, facet, neighbor, bestfacet, dist, mindist, maxdist )
    test neighbor of facet for qh_findbestneighbor()
    if testcentrum,
      tests centrum (assumes it is defined)
    else 
      tests vertices

  returns:
    if a better facet (i.e., vertices/centrum of facet closer to neighbor)
      updates bestfacet, dist, mindist, and maxdist
*/
void qh_findbest_test (boolT testcentrum, facetT *facet, facetT *neighbor,
      facetT **bestfacet, realT *distp, realT *mindistp, realT *maxdistp) {
  realT dist, mindist, maxdist;

  if (testcentrum) {
    zzinc_(Zbestdist);
    qh_distplane(facet->center, neighbor, &dist);
    dist *= qh hull_dim; /* estimate furthest vertex */
    if (dist < 0) {
      maxdist= 0;
      mindist= dist;
      dist= -dist;
    }else
      maxdist= dist;
  }else
    dist= qh_getdistance (facet, neighbor, &mindist, &maxdist);
  if (dist < *distp) {
    *bestfacet= neighbor;
    *mindistp= mindist;
    *maxdistp= maxdist;
    *distp= dist;
  }
} /* findbest_test */

/*-<a                             href="qh-c.htm#merge"
  >-------------------------------</a><a name="findbestneighbor">-</a>
  
  qh_findbestneighbor( facet, dist, mindist, maxdist )
    finds best neighbor (least dist) of a facet for merging

  returns:
    returns min and max distances and their max absolute value
  
  notes:
    avoids merging old into new
    assumes ridge->nonconvex only set on one ridge between a pair of facets
    could use an early out predicate but not worth it

  design:
    if a large facet
      will test centrum
    else
      will test vertices
    if a large facet
      test nonconvex neighbors for best merge
    else
      test all neighbors for the best merge
    if testing centrum
      get distance information
*/
facetT *qh_findbestneighbor(facetT *facet, realT *distp, realT *mindistp, realT *maxdistp) {
  facetT *neighbor, **neighborp, *bestfacet= NULL;
  ridgeT *ridge, **ridgep;
  boolT nonconvex= True, testcentrum= False;
  int size= qh_setsize (facet->vertices);

  *distp= REALmax;
  if (size > qh_BESTcentrum2 * qh hull_dim + qh_BESTcentrum) {
    testcentrum= True;
    zinc_(Zbestcentrum);
    if (!facet->center)
       facet->center= qh_getcentrum (facet);
  }
  if (size > qh hull_dim + qh_BESTnonconvex) {
    FOREACHridge_(facet->ridges) {
      if (ridge->nonconvex) {
        neighbor= otherfacet_(ridge, facet);
	qh_findbest_test (testcentrum, facet, neighbor,
			  &bestfacet, distp, mindistp, maxdistp);
      }
    }
  }
  if (!bestfacet) {     
    nonconvex= False;
    FOREACHneighbor_(facet)
      qh_findbest_test (testcentrum, facet, neighbor,
			&bestfacet, distp, mindistp, maxdistp);
  }
  if (!bestfacet) {
    fprintf (qh ferr, "qhull internal error (qh_findbestneighbor): no neighbors for f%d\n", facet->id);
    
    qh_errexit (qh_ERRqhull, facet, NULL);
  }
  if (testcentrum) 
    qh_getdistance (facet, bestfacet, mindistp, maxdistp);
  trace3((qh ferr, "qh_findbestneighbor: f%d is best neighbor for f%d testcentrum? %d nonconvex? %d dist %2.2g min %2.2g max %2.2g\n",
     bestfacet->id, facet->id, testcentrum, nonconvex, *distp, *mindistp, *maxdistp));
  return(bestfacet);
} /* findbestneighbor */


/*-<a                             href="qh-c.htm#merge"
  >-------------------------------</a><a name="flippedmerges">-</a>
  
  qh_flippedmerges( facetlist, wasmerge )
    merge flipped facets into best neighbor
    assumes qh.facet_mergeset at top of temporary stack

  returns:
    no flipped facets on facetlist
    sets wasmerge if merge occurred
    degen/redundant merges passed through

  notes:
    othermerges not needed since qh.facet_mergeset is empty before & after
      keep it in case of change

  design:
    append flipped facets to qh.facetmergeset
    for each flipped merge
      find best neighbor
      merge facet into neighbor
      merge degenerate and redundant facets
    remove flipped merges from qh.facet_mergeset
*/
void qh_flippedmerges(facetT *facetlist, boolT *wasmerge) {
  facetT *facet, *neighbor, *facet1;
  realT dist, mindist, maxdist;
  mergeT *merge, **mergep;
  setT *othermerges;
  int nummerge=0;

  trace4((qh ferr, "qh_flippedmerges: begin\n"));
  FORALLfacet_(facetlist) {
    if (facet->flipped && !facet->visible) 
      qh_appendmergeset (facet, facet, MRGflip, NULL);
  }
  othermerges= qh_settemppop(); /* was facet_mergeset */
  qh facet_mergeset= qh_settemp (qh TEMPsize);
  qh_settemppush (othermerges);
  FOREACHmerge_(othermerges) {
    facet1= merge->facet1;
    if (merge->type != MRGflip || facet1->visible) 
      continue;
    if (qh TRACEmerge-1 == zzval_(Ztotmerge))
      qhmem.IStracing= qh IStracing= qh TRACElevel;
    neighbor= qh_findbestneighbor (facet1, &dist, &mindist, &maxdist);
    trace0((qh ferr, "qh_flippedmerges: merge flipped f%d into f%d dist %2.2g during p%d\n",
      facet1->id, neighbor->id, dist, qh furthest_id));
    qh_mergefacet (facet1, neighbor, &mindist, &maxdist, !qh_MERGEapex);
    nummerge++;
    if (qh PRINTstatistics) {
      zinc_(Zflipped);
      wadd_(Wflippedtot, dist);
      wmax_(Wflippedmax, dist);
    }
    qh_merge_degenredundant();
  }
  FOREACHmerge_(othermerges) {
    if (merge->facet1->visible || merge->facet2->visible)
      qh_memfree (merge, sizeof(mergeT));
    else
      qh_setappend (&qh facet_mergeset, merge);
  }
  qh_settempfree (&othermerges);
  if (nummerge)
    *wasmerge= True;
  trace1((qh ferr, "qh_flippedmerges: merged %d flipped facets into a good neighbor\n", nummerge));
} /* flippedmerges */


/*-<a                             href="qh-c.htm#merge"
  >-------------------------------</a><a name="forcedmerges">-</a>
  
  qh_forcedmerges( wasmerge )
    merge duplicated ridges

  returns:
    removes all duplicate ridges on facet_mergeset
    wasmerge set if merge
    qh.facet_mergeset may include non-forced merges (none for now)
    qh.degen_mergeset includes degen/redun merges

  notes: 
    duplicate ridges occur when the horizon is pinched,
        i.e. a subridge occurs in more than two horizon ridges.
     could rename vertices that pinch the horizon
    assumes qh_merge_degenredundant() has not be called
    othermerges isn't needed since facet_mergeset is empty afterwards
      keep it in case of change

  design:
    for each duplicate ridge
      find current facets by chasing f.replace links
      determine best direction for facet
      merge one facet into the other
      remove duplicate ridges from qh.facet_mergeset
*/
void qh_forcedmerges(boolT *wasmerge) {
  facetT *facet1, *facet2;
  mergeT *merge, **mergep;
  realT dist1, dist2, mindist1, mindist2, maxdist1, maxdist2;
  setT *othermerges;
  int nummerge=0, numflip=0;

  if (qh TRACEmerge-1 == zzval_(Ztotmerge))
    qhmem.IStracing= qh IStracing= qh TRACElevel;
  trace4((qh ferr, "qh_forcedmerges: begin\n"));  
  othermerges= qh_settemppop(); /* was facet_mergeset */
  qh facet_mergeset= qh_settemp (qh TEMPsize);
  qh_settemppush (othermerges);
  FOREACHmerge_(othermerges) {
    if (merge->type != MRGridge) 
    	continue;
    facet1= merge->facet1;
    facet2= merge->facet2;
    while (facet1->visible)    	 /* must exist, no qh_merge_degenredunant */
      facet1= facet1->f.replace; /* previously merged facet */
    while (facet2->visible)
      facet2= facet2->f.replace; /* previously merged facet */
    if (facet1 == facet2)
      continue;
    if (!qh_setin (facet2->neighbors, facet1)) {
      fprintf (qh ferr, "qhull internal error (qh_forcedmerges): f%d and f%d had a duplicate ridge but as f%d and f%d they are no longer neighbors\n",
	       merge->facet1->id, merge->facet2->id, facet1->id, facet2->id);
      qh_errexit2 (qh_ERRqhull, facet1, facet2);
    }
    if (qh TRACEmerge-1 == zzval_(Ztotmerge))
      qhmem.IStracing= qh IStracing= qh TRACElevel;
    dist1= qh_getdistance (facet1, facet2, &mindist1, &maxdist1);
    dist2= qh_getdistance (facet2, facet1, &mindist2, &maxdist2);
    trace0((qh ferr, "qh_forcedmerges: duplicate ridge between f%d and f%d, dist %2.2g and reverse dist %2.2g during p%d\n",
	    facet1->id, facet2->id, dist1, dist2, qh furthest_id));
    if (dist1 < dist2) 
      qh_mergefacet (facet1, facet2, &mindist1, &maxdist1, !qh_MERGEapex);
    else {
      qh_mergefacet (facet2, facet1, &mindist2, &maxdist2, !qh_MERGEapex);
      dist1= dist2;
      facet1= facet2;
    }
    if (facet1->flipped) {
      zinc_(Zmergeflipdup);
      numflip++;
    }else
      nummerge++;
    if (qh PRINTstatistics) {
      zinc_(Zduplicate);
      wadd_(Wduplicatetot, dist1);
      wmax_(Wduplicatemax, dist1);
    }
  }
  FOREACHmerge_(othermerges) {
    if (merge->type == MRGridge)
      qh_memfree (merge, sizeof(mergeT));
    else
      qh_setappend (&qh facet_mergeset, merge);
  }
  qh_settempfree (&othermerges);
  if (nummerge)
    *wasmerge= True;
  trace1((qh ferr, "qh_forcedmerges: merged %d facets and %d flipped facets across duplicated ridges\n", 
                nummerge, numflip));
} /* forcedmerges */


/*-<a                             href="qh-c.htm#merge"
  >-------------------------------</a><a name="getmergeset">-</a>
  
  qh_getmergeset( facetlist )
    determines nonconvex facets on facetlist
    tests !tested ridges and nonconvex ridges of !tested facets

  returns:
    returns sorted qh.facet_mergeset of facet-neighbor pairs to be merged
    all ridges tested
  
  notes:
    assumes no nonconvex ridges with both facets tested
    uses facet->tested/ridge->tested to prevent duplicate tests
    can not limit tests to modified ridges since the centrum changed
    uses qh.visit_id
  
  see:
    qh_getmergeset_initial()

  design:
    for each facet on facetlist
      for each ridge of facet
        if untested ridge
          test ridge for convexity
          if non-convex
            append ridge to qh.facet_mergeset
    sort qh.facet_mergeset by angle  
*/
void qh_getmergeset(facetT *facetlist) {
  facetT *facet, *neighbor, **neighborp;
  ridgeT *ridge, **ridgep;
  int nummerges;
  
  nummerges= qh_setsize (qh facet_mergeset);
  trace4((qh ferr, "qh_getmergeset: started.\n"));
  qh visit_id++;
  FORALLfacet_(facetlist) {
    if (facet->tested)
      continue;
    facet->visitid= qh visit_id;
    facet->tested= True;  /* must be non-simplicial due to merge */
    FOREACHneighbor_(facet)
      neighbor->seen= False;
    FOREACHridge_(facet->ridges) {
      if (ridge->tested && !ridge->nonconvex)
	continue;
      /* if tested & nonconvex, need to append merge */
      neighbor= otherfacet_(ridge, facet);
      if (neighbor->seen) {