poly.c

关于网格剖分的

  notes:
    assumes at least firstindex+1 elements
    assumes skipelem is NULL, in set, or part of hash
    
    hashes memory addresses which may change over different runs of the same data
    using sum for hash does badly in high d
*/
unsigned qh_gethash (int hashsize, setT *set, int size, int firstindex, void *skipelem) {
  void **elemp= SETelemaddr_(set, firstindex, void);
  ptr_intT hash = 0, elem;
  int i;

  switch (size-firstindex) {
  case 1:
    hash= (ptr_intT)(*elemp) - (ptr_intT) skipelem;
    break;
  case 2:
    hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] - (ptr_intT) skipelem;
    break;
  case 3:
    hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] + (ptr_intT)elemp[2]
      - (ptr_intT) skipelem;
    break;
  case 4:
    hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] + (ptr_intT)elemp[2]
      + (ptr_intT)elemp[3] - (ptr_intT) skipelem;
    break;
  case 5:
    hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] + (ptr_intT)elemp[2]
      + (ptr_intT)elemp[3] + (ptr_intT)elemp[4] - (ptr_intT) skipelem;
    break;
  case 6:
    hash= (ptr_intT)(*elemp) + (ptr_intT)elemp[1] + (ptr_intT)elemp[2]
      + (ptr_intT)elemp[3] + (ptr_intT)elemp[4]+ (ptr_intT)elemp[5]
      - (ptr_intT) skipelem;
    break;
  default:
    hash= 0;
    i= 3;
    do {     /* this is about 10% in 10-d */
      if ((elem= (ptr_intT)*elemp++) != (ptr_intT)skipelem) {
        hash ^= (elem << i) + (elem >> (32-i));
	i += 3;
	if (i >= 32)
	  i -= 32;
      }
    }while(*elemp);
    break;
  }
  hash %= (ptr_intT) hashsize;
  /* hash= 0; for debugging purposes */
  return hash;
} /* gethash */

/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="makenewfacet">-</a>
  
  qh_makenewfacet( vertices, toporient, horizon )
    creates a toporient? facet from vertices

  returns:
    returns newfacet
      adds newfacet to qh.facet_list 
      newfacet->neighbor= horizon, but not vice versa
      newfacet->vertices= vertices
    newvertex_list updated with vertices
*/
facetT *qh_makenewfacet(setT *vertices, boolT toporient,facetT *horizon) {
  facetT *newfacet;
  vertexT *vertex, **vertexp;

  FOREACHvertex_(vertices) {
    if (!vertex->newlist) {
      qh_removevertex (vertex);
      qh_appendvertex (vertex);
    }
  }
  newfacet= qh_newfacet();
  newfacet->vertices= vertices;
  newfacet->toporient= toporient;
  qh_setappend(&(newfacet->neighbors), horizon);
  qh_appendfacet(newfacet);
  return(newfacet);
} /* makenewfacet */


/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="makenewplanes">-</a>
  
  qh_makenewplanes()
    make new hyperplanes for facets on qh.newfacet_list

  returns:
    all facets have hyperplanes or are marked for   merging
    doesn't create hyperplane if horizon is coplanar (will merge)
    updates qh.min_vertex if qh.JOGGLEmax

  notes:
    facet->f.samecycle is defined for facet->mergehorizon facets
*/
void qh_makenewplanes (void /* newfacet_list */) {
  facetT *newfacet;

  FORALLnew_facets {
    if (!newfacet->mergehorizon)
      qh_setfacetplane (newfacet);  
  }
  if (qh JOGGLEmax < REALmax/2)  
    minimize_(qh min_vertex, -wwval_(Wnewvertexmax));
} /* makenewplanes */

/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="makenew_nonsimplicial">-</a>
  
  qh_makenew_nonsimplicial( visible, apex, numnew )
    make new facets for ridges of a visible facet
    
  returns:
    first newfacet, bumps numnew as needed
    attaches new facets if !qh.ONLYgood
    marks ridge neighbors for simplicial visible
    if (qh.ONLYgood)
      ridges on newfacet, horizon, and visible
    else
      ridge and neighbors between newfacet and   horizon
      visible facet's ridges are deleted    

  notes:
    qh.visit_id if visible has already been processed
    sets neighbor->seen for building f.samecycle
      assumes all 'seen' flags initially false
    
  design:
    for each ridge of visible facet
      get neighbor of visible facet
      if neighbor was already processed
        delete the ridge (will delete all visible facets later)
      if neighbor is a horizon facet
        create a new facet
        if neighbor coplanar
          adds newfacet to f.samecycle for later merging
        else 
          updates neighbor's neighbor set
          (checks for non-simplicial facet with multiple ridges to visible facet)
        updates neighbor's ridge set
        (checks for simplicial neighbor to non-simplicial visible facet)
          
*/
#ifndef qh_NOmerge
facetT *qh_makenew_nonsimplicial (facetT *visible, vertexT *apex, int *numnew) {
  void **freelistp; /* used !qh_NOmem */
  ridgeT *ridge, **ridgep;
  facetT *neighbor, *newfacet= NULL, *samecycle;
  setT *vertices;
  boolT toporient;
  int ridgeid;

  FOREACHridge_(visible->ridges) {
    ridgeid= ridge->id;
    neighbor= otherfacet_(ridge, visible);
    if (neighbor->visible) {
      if (!qh ONLYgood) {
        if (neighbor->visitid == qh visit_id) {
          qh_setfree (&(ridge->vertices));  /* delete on 2nd visit */
	  qh_memfree_(ridge, sizeof(ridgeT), freelistp);
	}
      }
    }else {  /* neighbor is an horizon facet */
      toporient= (ridge->top == visible);
      vertices= qh_setnew (qh hull_dim); /* makes sure this is quick */
      qh_setappend (&vertices, apex);
      qh_setappend_set (&vertices, ridge->vertices);
      newfacet= qh_makenewfacet(vertices, toporient, neighbor);
      (*numnew)++;
      if (neighbor->coplanar) {
	newfacet->mergehorizon= True;
        if (!neighbor->seen) {
          newfacet->f.samecycle= newfacet;
          neighbor->f.newcycle= newfacet;
        }else {
          samecycle= neighbor->f.newcycle;
          newfacet->f.samecycle= samecycle->f.samecycle;
          samecycle->f.samecycle= newfacet;
	}
      }
      if (qh ONLYgood) {
        if (!neighbor->simplicial)
 	  qh_setappend(&(newfacet->ridges), ridge);
      }else {  /* qh_attachnewfacets */
        if (neighbor->seen) {
	  if (neighbor->simplicial) {
	    fprintf (qh ferr, "qhull internal error (qh_makenew_nonsimplicial): simplicial f%d sharing two ridges with f%d\n", 
	           neighbor->id, visible->id);
	    qh_errexit2 (qh_ERRqhull, neighbor, visible);
	  }
	  qh_setappend (&(neighbor->neighbors), newfacet);
	}else
          qh_setreplace (neighbor->neighbors, visible, newfacet);
        if (neighbor->simplicial) {
          qh_setdel (neighbor->ridges, ridge);
          qh_setfree (&(ridge->vertices)); 
	  qh_memfree (ridge, sizeof(ridgeT));
	}else {
 	  qh_setappend(&(newfacet->ridges), ridge);
 	  if (toporient)
 	    ridge->top= newfacet;
 	  else
 	    ridge->bottom= newfacet;
 	}
      trace4((qh ferr, "qh_makenew_nonsimplicial: created facet f%d from v%d and r%d of horizon f%d\n",
	    newfacet->id, apex->id, ridgeid, neighbor->id));
      }
    }
    neighbor->seen= True;        
  } /* for each ridge */
  if (!qh ONLYgood)
    SETfirst_(visible->ridges)= NULL;
  return newfacet;
} /* makenew_nonsimplicial */
#else /* qh_NOmerge */
facetT *qh_makenew_nonsimplicial (facetT *visible, vertexT *apex, int *numnew) {
  return NULL;
}
#endif /* qh_NOmerge */

/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="makenew_simplicial">-</a>
  
  qh_makenew_simplicial( visible, apex, numnew )
    make new facets for simplicial visible facet and apex

  returns:
    attaches new facets if (!qh.ONLYgood)
      neighbors between newfacet and horizon

  notes:
    nop if neighbor->seen or neighbor->visible (see qh_makenew_nonsimplicial)

  design:
    locate neighboring horizon facet for visible facet
    determine vertices and orientation
    create new facet
    if coplanar,
      add new facet to f.samecycle
    update horizon facet's neighbor list        
*/
facetT *qh_makenew_simplicial (facetT *visible, vertexT *apex, int *numnew) {
  facetT *neighbor, **neighborp, *newfacet= NULL;
  setT *vertices;
  boolT flip, toporient;
  int horizonskip, visibleskip;

  FOREACHneighbor_(visible) {
    if (!neighbor->seen && !neighbor->visible) {
      vertices= qh_facetintersect(neighbor,visible, &horizonskip, &visibleskip, 1);
      SETfirst_(vertices)= apex;
      flip= ((horizonskip & 0x1) ^ (visibleskip & 0x1));
      if (neighbor->toporient)         
	toporient= horizonskip & 0x1;
      else
	toporient= (horizonskip & 0x1) ^ 0x1;
      newfacet= qh_makenewfacet(vertices, toporient, neighbor);
      (*numnew)++;
      if (neighbor->coplanar && (qh PREmerge || qh MERGEexact)) {
#ifndef qh_NOmerge
	newfacet->f.samecycle= newfacet;
	newfacet->mergehorizon= True;
#endif
      }
      if (!qh ONLYgood)
        SETelem_(neighbor->neighbors, horizonskip)= newfacet;
      trace4((qh ferr, "qh_makenew_simplicial: create facet f%d top %d from v%d and horizon f%d skip %d top %d and visible f%d skip %d, flip? %d\n",
	    newfacet->id, toporient, apex->id, neighbor->id, horizonskip,
	      neighbor->toporient, visible->id, visibleskip, flip));
    }
  }
  return newfacet;
} /* makenew_simplicial */

/*-<a                             href="qh-c.htm#poly"
  >-------------------------------</a><a name="matchneighbor">-</a>
  
  qh_matchneighbor( newfacet, newskip, hashsize, hashcount )
    either match subridge of newfacet with neighbor or add to hash_table

  returns:
    duplicate ridges are unmatched and marked by qh_DUPLICATEridge

  notes:
    ridge is newfacet->vertices w/o newskip vertex
    do not allocate memory (need to free hash_table cleanly)
    uses linear hash chains
  
  see also:
    qh_matchduplicates

  design:
    for each possible matching facet in qh.hash_table
      if vertices match
        set ismatch, if facets have opposite orientation
        if ismatch and matching facet doesn't have a match
          match the facets by updating their neighbor sets
        else
          indicate a duplicate ridge
          set facet hyperplane for later testing
          add facet to hashtable
          unless the other facet was already a duplicate ridge
            mark both facets with a duplicate ridge
            add other facet (if defined) to hash table
*/
void qh_matchneighbor (facetT *newfacet, int newskip, int hashsize, int *hashcount) {
  boolT newfound= False;   /* True, if new facet is already in hash chain */
  boolT same, ismatch;
  int hash, scan;
  facetT *facet, *matchfacet;
  int skip, matchskip;

  hash= (int)qh_gethash (hashsize, newfacet->vertices, qh hull_dim, 1, 
                     SETelem_(newfacet->vertices, newskip));
  trace4((qh ferr, "qh_matchneighbor: newfacet f%d skip %d hash %d hashcount %d\n",
	  newfacet->id, newskip, hash, *hashcount));
  zinc_(Zhashlookup);
  for (scan= hash; (facet= SETelemt_(qh hash_table, scan, facetT)); 
       scan= (++scan >= hashsize ? 0 : scan)) {
    if (facet == newfacet) {
      newfound= True;
      continue;
    }
    zinc_(Zhashtests);
    if (qh_matchvertices (1, newfacet->vertices, newskip, facet->vertices, &skip, &same)) {
      if (SETelem_(newfacet->vertices, newskip) == 
          SETelem_(facet->vertices, skip)) {
        qh_precision ("two facets with the same vertices");
        fprintf (qh ferr, "qhull precision error: Vertex sets are the same for f%d and f%d.  Can not force output.\n",
          facet->id, newfacet->id);
        qh_errexit2 (qh_ERRprec, facet, newfacet);
      }
      ismatch= (same == (newfacet->toporient ^ facet->toporient));
      matchfacet= SETelemt_(facet->neighbors, skip, facetT);
      if (ismatch && !matchfacet) {
        SETelem_(facet->neighbors, skip)= newfacet;
        SETelem_(newfacet->neighbors, newskip)= facet;
        (*hashcount)--;
        trace4((qh ferr, "qh_matchneighbor: f%d skip %d matched with new f%d skip %d\n",
           facet->id, skip, newfacet->id, newskip));
        return;
      }
      if (!qh PREmerge && !qh MERGEexact) {
        qh_precision ("a ridge with more than two neighbors");
	fprintf (qh ferr, "qhull precision error: facets f%d, f%d and f%d meet at a ridge with more than 2 neighbors.  Can not continue.\n",
		 facet->id, newfacet->id, getid_(matchfacet));
	qh_errexit2 (qh_ERRprec, facet, newfacet);
      }
      SETelem_(newfacet->neighbors, newskip)= qh_DUPLICATEridge;
      newfacet->dupridge= True;
      if (!newfacet->normal)
	qh_setfacetplane (newfacet);
      qh_addhash (newfacet, qh hash_table, hashsize, hash);
      (*hashcount)++;
      if (!facet->normal)
	qh_setfacetplane (facet);
      if (matchfacet != qh_DUPLICATEridge) {
	SETelem_(facet->neighbors, skip)= qh_DUPLICATEridge;
	facet->dupridge= True;
	if (!facet->normal)
	  qh_setfacetplane (facet);
	if (matchfacet) {
	  matchskip= qh_setindex (matchfacet->neighbors, facet);
	  SETelem_(matchfacet->neighbors, matchskip)= qh_DUPLICATEridge;
	  matchfacet->dupridge= True;
	  if (!matchfacet->normal)
	    qh_setfacetplane (matchfacet);
	  qh_addhash (matchfacet, qh hash_table, hashsize, hash);
	  *hashcount += 2;
	}
      }
      trace4((qh ferr, "qh_matchneighbor: new f%d skip %d duplicates ridge for f%d skip %d matching f%d ismatch %d at hash %d\n",
	   newfacet->id, newskip, facet->id, skip, 
	   (matchfacet == qh_DUPLICATEridge ? -2 : getid_(matchfacet)), 
	   ismatch, hash));
      return; /* end of duplicate ridge */
    }
  }
  if (!newfound) 
    SETelem_(qh hash_table, scan)= newfacet;  /* same as qh_addhash */
  (*hashcount)++;
  trace4((qh ferr, "qh_matchneighbor: no match for f%d skip %d at hash %d\n",
           newfacet->id, newskip, hash));
} /* matchneighbor */